/* Zeiraphera two site Parasitoid & Quality model */ /**************************************************************************/ /* This template file is for setting up a difference equation model for */ /* fitting to data: */ /* The user supplies routines map(), initst(), initfit(), and can call */ /* uc(), and uf() for unknown constants and functions. */ /**************************************************************************/ /* Essential include files */ #include #include #include #include "ddefit95.h" #include "rangen.h" #include "spline.h" #include "map.h" // This file is for discrete time models only /***************************************************************************/ /* Put your global CONSTANTS next. e.g. double IamAglobalVariable=2.0; */ /***************************************************************************/ /***************************************************************************/ /* Problem specific routines */ /***************************************************************************/ void map(s,c,t,swno) double *s,*c,t;int swno; /* Routine which takes the state vector forward 1 time unit Nt - hosts St - survival percentages Qt - Quality Lt - Needle length b+c*Qt Pt - Parasitoid density - this is not a state variable - just an intermediate */ { double Nte,Qte,Ste,Nts,Qts,Sts, // states Nte1,Qte1,Ste1,Nts1,Qts1,Sts1, // updated states Lte,Lts,Lte1,Lts1, // intermediate `states' a,delta,r,Ke,Ks,alpha,gamma,b,C,h,w,SUR,kick; // parameters int i; for (i=0;i<6;i++) if (s[i]<0.0) s[i]=0.0; Nte=s[0];Lte=s[1];Ste=s[2];Nts=s[3];Lts=s[4];Sts=s[5]; a=uc(0);delta=uc(1);r=uc(2);Ke=uc(3);Ks=uc(4); alpha=uc(5);gamma=uc(6);b=uc(7);C=uc(8);h=uc(9);w=uc(10);SUR=uc(17); kick=uc(18); if (Nte<0.0) Nte=0.0; if (Nts<0.0) Nts=0.0; // host untransform Nte=pow(Nte,10.0);Nts=pow(Nts,10.0); // host untransform Ste/=SUR;Sts/=SUR; // deal with under-reporting of Parasitism rate // Turn proportion parasitised into proportion surviving parasitism .... Ste=1.0-Ste;Sts=1.0-Sts; if (Ste>1) Ste=1.0; if (Sts>1.0) Sts=1.0; // Now get Qt from Lt Qte=(Lte-b)/C;Qts=(Lts-b)/C; // Engadine update if (delta>0.0) Lte=r*Qte/(delta+Qte); else Lte=r; // food quality effect on host Nte1=Nte*Lte*exp( - Nte/Ke)*Ste; // Host Ste1=exp(-a*Nte*(1-Ste)/(1+a*w*Nte*(1-Ste)+a*h*Nte1)); // Survival of parasitism if (Nte>0.0) Nte=Nte/(gamma+Nte); else Nte=0.0; // population effect on quality Qte1=(1-alpha)*(1-Nte)+alpha*Qte; // quality // Sils update if (delta>0.0) Lts=r*Qts/(delta+Qts); else Lts=r; // food quality effect on host Nts1=Nts*Lts*exp( - Nts/Ks)*Sts; // Host if (t>37.5&&t<38.5) Nts1+=kick; Sts1=exp(-a*Nts*(1-Sts)/(1+a*w*Nts*(1-Sts)+a*h*Nts1)); // Survival of parasitism if (Nts>0.0) Nts=Nts/(gamma+Nts); else Nts=0.0; // population effect on quality Qts1=(1-alpha)*(1-Nts)+alpha*Qts; // quality // It is now necessary to extensively untransform things..... Ste1=1.0-Ste1;Sts1=1.0-Sts1; Ste1*=SUR;Sts1*=SUR; if (Ste1>1.0) Ste1=1.0; if (Sts1>1.0) Sts1=1.0; if (Ste1<0.0) Ste1=0.0; if (Sts1<0.0) Sts1=0.0; if (Nte1<0.0) Nte1=0.0;if(Nts1<0.0) Nts1=0.0; Nte1=pow(Nte1,0.1);Nts1=pow(Nts1,0.1); // host transform ET /* transform quality back to needle length */ Lte1=b + C*Qte1;Lts1=b+C*Qts1; s[0]=Nte1;s[1]=Lte1;s[2]=Ste1;s[3]=Nts1;s[4]=Lts1;s[5]=Sts1; } void initst(s,c,t) double *s,*c,t; /* initialise state variables and any global constants here, using c or uc*/ { double SUR; s[0]=uc(11); // Nte s[1]=uc(12); // Lte s[2]=uc(13); // Ste s[3]=uc(14); // Nts s[4]=uc(15); // Lts s[5]=uc(16); // Sts SUR=uc(17); if (s[0]<0.0) s[0]=0.0; if (s[3]<0.0) s[3]=0.0; s[0]=pow(s[0],0.1);s[3]=pow(s[3],0.1); // transform s[2]*=SUR;s[5]*=SUR; if (s[2]>1.0) s[2]=1.0; if (s[5]>1.0) s[5]=1.0; if (s[2]<0.0) s[2]=0.0; if (s[5]<0.0) s[5]=0.0; } /**************************************************************************/ /* Fit specification routines */ /**************************************************************************/ void initfit(d) ddefit_control *d; /* This is where you have to supply details about the model fit, and the graphical display of the model fit. * This routine must not use: uc(),uf(),guf(),iuf() or call any functions which do. If you call these functions in initfit() you will be comitting memory errors. * It is a good idea not to edit the comments out of this routine. * NOTE if you are supplying a floating point number (a double) it must contain a decimal point or e. If you are supplying an integer it must not! C is very unforgiving on this point. */ { /***************** The following values must be supplied ****************/ d->discrete=1; // MUST be set to 1 to indicate a discrete time model /************************************************************************/ /* The following define numbers of things ....... */ /************************************************************************/ d->no_uc=19; // Number of unknown parameters (integer<5000) d->no_uf=0; // Number of unknown functions d->no_c=0; // Number of (known) model coefficients (integer<5000) d->no_s=6; // Number of state variables (integer<5000) d->no_fit=4; // Number of variables for which there is data to fit d->t0=0.0; // iteration start time d->t1=400.0; // Default end time, when simulating d->dout=1.0; // don't change this in a discrete model! /**************************************************************************/ /* Now choose the error model and fitting method. */ /* There are 3 fitting methods (all constrained): */ /* 0 - Gauss-Newton with steepest descent backup. This is good on small */ /* residual problems. It tries both GN and steepest descent direction */ /* at each step. Can be slow on large residual problems. */ /* 1 - Quasi-Newton. Usually faster than Gauss-Newton on moderate to */ /* large residual problems. */ /* 2 - Iterative least squares/quadratic programming. Very fast on well */ /* behaved problems - can diverge otherwise. */ /* There are 3 error models: */ /* 0 - Error independent of mean (normal) - single shot fitting - fastest */ /* and most robust. ONLY this one works when estimating smoothing */ /* parameters. */ /* 1 - Variance proportional to mean (Poisson). Iteratively re-weighted */ /* least sqaures objective used. Convergence can be slow (and is not */ /* guaranteed - these are not GLMs). */ /* 2 - Standard deviation proportional to mean (constant CV). Same */ /* method and comments as model 1. */ /* Finally specify the number of bootstrap restarts to use, this is very */ /* useful with difficult objectives, but is ignored whenthere are */ /* smoothing parameters to be estimated. */ /**************************************************************************/ d->fitmethod=0; d->errors=0; d->bsr_reps=0; /**************************************************************************/ /* Now supply default values and names for the unknown parameters... */ /**************************************************************************/ d->uco[0]=0.684442; d->cname[0]="a"; d->uco[1]=44.5989; d->cname[1]="delta"; d->uco[2]=410.744; d->cname[2]="r"; d->uco[3]=83.7881; d->cname[3]="Ke"; d->uco[4]=90.4732; d->cname[4]="Ks"; d->uco[5]=0.240882; d->cname[5]= "alpha"; d->uco[6]=260.551; d->cname[6]= "gamma"; d->uco[7]=1e-4; d->cname[7]= "b"; d->uco[8]=1.503; d->cname[8]= "C"; d->uco[9]=1e-5; d->cname[9]= "h"; d->uco[10]=0.139018; d->cname[10]= "w"; d->uco[11]=5.2578; d->cname[11]= "Nte"; d->uco[12]=1.20055; d->cname[12]= "Lte"; d->uco[13]=0.000506; d->cname[13]= "Ste"; d->uco[14]=1.51478; d->cname[14]= "Nts"; d->uco[15]=1.49966; d->cname[15]= "Lts"; d->uco[16]=0.533476; d->cname[16]= "Sts"; d->uco[17]=0.706272; d->cname[17]= "SUR"; d->uco[18]=0.0; d->cname[18]="ignore"; /**************************************************************************/ /* Supply the type of constraints and any upper and lower bounds on the */ /* uc's. An unknown coefficient can be of the following types: */ /* * UNBOUND - no constraints */ /* * B_BELOW - bounded below: lower bound must be supplied in d->uclb[i] */ /* * B_ABOVE - bounded above: upper bound must be supplied in d->uclb[i] */ /* * B_ABOVE|B_BELOW - bounded above and below. */ /* * FIXED - will not be changed in fitting. */ /**************************************************************************/ d->uctype[0]=B_BELOW|B_ABOVE;d->uclb[0]= 0.0;d->ucub[0]=100.0; // a d->uctype[1]=B_BELOW|B_ABOVE;d->uclb[1]=0.0;d->ucub[1]=1000.0; // delta d->uctype[2]=B_BELOW|B_ABOVE;d->uclb[2]=0.0;d->ucub[2]=10000.0; // r d->uctype[3]=B_BELOW|B_ABOVE;d->uclb[3]=0.001;d->ucub[3]=10000.0; // Ke d->uctype[4]=B_BELOW|B_ABOVE;d->uclb[4]=0.001;d->ucub[4]=10000.0; // Ks d->uctype[5]=B_BELOW|B_ABOVE; d->uclb[5]= 0.00001; d->ucub[5]=.9999; // alpha d->uctype[6]=B_BELOW|B_ABOVE; d->uclb[6]=0.00001; d->ucub[6]=1000.0; // gamma d->uctype[7]=B_BELOW|B_ABOVE; d->uclb[7]= 0.0; d->ucub[7]=100.0; // b d->uctype[8]=B_BELOW|B_ABOVE; d->uclb[8]=0.5; d->ucub[8]=100.0; // C d->uctype[9]=B_BELOW|B_ABOVE; d->uclb[9]=0.0; d->ucub[9]=1000.0; // h d->uctype[10]=B_BELOW|B_ABOVE;d->uclb[10]=0.0; d->ucub[10]=1000.0; // w d->uctype[11]=B_BELOW|B_ABOVE;d->uclb[11]=0.001; d->ucub[11]=1000.0; // Nte d->uctype[12]=B_BELOW|B_ABOVE;d->uclb[12]=1.2; d->ucub[12]=1.5; // Lte d->uctype[13]=B_BELOW|B_ABOVE;d->uclb[13]=0.0; d->ucub[13]=1.0; // Ste d->uctype[14]=B_BELOW|B_ABOVE;d->uclb[14]=0.001; d->ucub[14]=1000.0; // Nts d->uctype[15]=B_BELOW|B_ABOVE;d->uclb[15]=1.2; d->ucub[15]=1.5; // Lts d->uctype[16]=B_BELOW|B_ABOVE;d->uclb[16]=0.0; d->ucub[16]=1.0; // Sts d->uctype[17]=B_BELOW|B_ABOVE;d->uclb[17]=0.7; d->ucub[17]=1.0; // SUR d->uctype[18]=FIXED; /**************************************************************************/ /* Definitions for unknown functions. */ /* d->ufdf[i] is maximum degrees of freedom for ith u.f. */ /* d->uft0[i] to d->uft0[i] defines the range of the argument u.f. */ /* d->uftype[i] can be UNBOUND, INCREASING, DECREASING, B_BELOW, B_ABOVE */ /* RE1, ORIGIN or FIXED or combinations of these using |. */ /* d->ufub[i] and d->uflb[i] give the optional upper an lower bounds. */ /* d->ufsp[i] is the smoothing parameter. If any of these are -ve then */ /* automatic selection will be used. */ /**************************************************************************/ /**************************************************************************/ /* This part deals with the equivalence of state variables and columns */ /* of data in the input file. */ /* index[i] is the state variable corresponding to input file column i */ /* where it is assumed that column 0 of the input file is time. */ /**************************************************************************/ d->index[1]=0; d->index[2]=2;d->index[3]=3;d->index[4]=4; d->dfile="zwrong.dat"; // Data file name d->wfile="zwrong.w"; // Weight file name /**************************************************************************/ /* It is possible to specify extra statistics of the input series that */ /* should be appended to the data vector as "extra" data to be fitted. */ /* The weight to be given to these statistics can also be specified. */ /* Available options are: */ /* NO_STATS - nothing (option to use if formal statistical stuff is to */ /* be done) */ /* STDEV - standard deviation of series. */ /* ABOLUTEDEV - Mean absulte deviation of data from series mean. */ /* MEANGRAD - Mean gradient between data points. */ /* MEANFREQ - Mean frequency of series. */ /* ACF - Autocorrelation function to first 10 lags. */ /* These can be combined, e.g. d->statistics=STDEV|ACF|MEANFREQ. */ /* In fact only ACF is really useful, in my experience. It is useful for */ /* finding parameters in the right ballpark for long cyclic series. */ /* Weights to give to statistics are specified in d->wstats. e.g. for the */ /* above example: */ /* d->wstats[STDEV]=100.0;d->wstats[MEANFREQ]=1.0;d->wstats[ACF]=30.0; */ /* would specify weights 100, 1 and 30 for the selected statistics. */ /**************************************************************************/ d->statistics=NO_STATS; d->wstats[ACF]=100.0; /**************************************************************************/ /* Now you get to specify how the model and data are to be displayed */ /**************************************************************************/ d->no_windows=2; // Number of output windows for model and data series /**************************************************************************/ /* Number of model series to be displayed in each window. If a series */ /* is to be fitted to data then that data is automatically displayed */ /* with it, in the same colour. */ /**************************************************************************/ d->lines[0]=3; d->lines[1]=3; /**************************************************************************/ /* Now specify exactly where each displayed model state variable series */ /* will be displayed, by filling in an element of the d->windex[] array */ /* for each variable to be displayed. For example if state variable 2 is */ /* to be displayed as the 0th curve in the 0th window and state */ /* variable 0 is the 1st curve in that window: */ /* d->windex[2].win=0;d->windex[2].cur=0; */ /* d->windex[0].win=0;d->windex[0].cur=1; */ /* If a state variable is not to be displayed then don't fill out its */ /* element of d->windex[]. */ /**************************************************************************/ d->windex[0].win=0;d->windex[0].cur=0; d->windex[1].win=0;d->windex[1].cur=1; d->windex[2].win=0;d->windex[2].cur=2; d->windex[3].win=1;d->windex[3].cur=0; d->windex[4].win=1;d->windex[4].cur=1; d->windex[5].win=1;d->windex[5].cur=2; /**************************************************************************/ /* You must specify initial y-axis ranges for each display window. */ /**************************************************************************/ d->range[0].y0=0.0;d->range[0].y1=2.0; d->range[1].y0=0.0;d->range[1].y1=2.0; /**************************************************************************/ /* Labels for state variables. You can label any model state variable, */ /* and these labels will be used to label the displayed model series. */ /**************************************************************************/ d->label[0]="E. Hosts";d->label[1]="E. Quality";d->label[2]="E. Parasites"; d->label[3]="S. Hosts";d->label[4]="S. Quality";d->label[5]="S. Parasites"; /**************************************************************************/ /* Labels for windows. You can give names to each output display window. */ /**************************************************************************/ d->wname[0]="Engadine"; d->wname[1]="Sils"; } double inuf(int i,double t) /* This is where unknown function intial values must be set. * i is the index of the unknown function and t is its argument. * Given this information you must return the value of the ith unknown function at t from this function. */ { double mu=30.0,sig=20.0; switch (i) { case 0: return(0.0); case 1: return(0.0); } return(0.0); } void inboot(bsctype *bsc) /* This is the routine for setting up bootstrapping, (not bootstrap restarting) * You can leave it blank, or set bsc->reps=0, if you don't want to bootstrap. * bsc->parametric: set to 0 for non-parametric, 1 for parametric or 2 for semi-parametric. * bsc->lumped: set to 0 to treat each fitted series separately for bootstrapping purposes, or to 1 to lump the stages in bootstrapping. * bsc->iocontrol: 1 to write replicate data to a file, 2 to read replicate data from a file 0 for no i/o of data. * bsc->fdname: name of the file to which replicate data (actually weights for no-para) will be written. * bsc->fpname: name of file to which replicate parameter estimates will be written. * bsc->SSname: (weighted) sum of squares is written to this file, for each rep. * bsc->reps: number of bootstrap replicates to generate and fit. * bsc->nextra: number of extra data series (see documentation - usually set to zero or leave out) * bsc->ignore: user initialised array of fit variables to ignore (see documentation) * bsc->restarts[0/1/2]: Number of bootstrap restarts per fit on first second and subsequent reps. * bsc->carry_p: set to 1 to carry the best fit parameter vector from the last rep on as the starting parameter vector for the next fit. Set to 0 to use best fit parameters to original data as starting parameters for each fit. Set to -1 to always use the original starting parameters. * bsc->n_start_files: you can read start parameters from a file: this specifies the number of start files to use. * bsc->start_file[]: array for start file names - you must initialise this, e.g. for a list of 3 files: bsc->start_file=(char **)calloc(3,sizeof(char*)); A * in front of the the first filename indicates that this is a parameter file generated by a previous set of bootstrap runs, with a different set of parameters for each b.s. rep. e.g. bsc->start_file[0]="*zd0.bsp"; tells the program that "zd0.bsp" is a text file with one set of starting parameters on each line. For more information please see the written documentation. */ { bsc->parametric=2; // semi-para bsc->lumped=0; bsc->iocontrol=2; bsc->fdname="zd0.bsy"; bsc->fpname="zfd0.bsp"; bsc->SSname="zfd0.SS"; bsc->reps=231; bsc->nextra=0; bsc->ignore=(int *)calloc(4,sizeof(int)); bsc->restarts[0]=40; bsc->restarts[1]=40; bsc->restarts[2]=40; bsc->carry_p=0; bsc->n_start_files=2; bsc->start_file=(char **)calloc(3,sizeof(char*)); bsc->start_file[0]="*zd0.bsp"; bsc->start_file[1]="estet.p"; bsc->start_file[2]="zfull.p"; }