/* This template file is for setting up a delay differential equation model
   that you would like to fit to data. The template file is based on the
   Solver program written by Bill Gurney, Pete Maas, Roger Nisbet and
   Steve Blythe, but there are important differences in the details:
   0) indeces start at 0 in the usual "C" style.
   1) Its all in C, not pascal.
   2) you have to store gradients of lagged variables, in order to maintain
      numerical probity in the adaptive stepsize control that is used.
      (to help, a pastgradient() function is provided)
*/

/* Essential include files */
#include <math.h>
#include <stdio.h>
#include "ddeq.h"
#include "ddefit95.h"
#include "spline.h"

/***************************************************************************/
/* Put your global CONSTANTS next. e.g. double IamAglobalVariable=2.0;     */
/* (Don't be tempted to use global variables - they will almost certainly  */
/* not work in the way you expect, - in particular grad is often called    */
/* with a value of t earlier than some previously used value)              */
/***************************************************************************/


/***************************************************************************/
/*             Problem specific routines                                   */
/***************************************************************************/


void switchfunctions(sw,s,c,t)
double *sw,*s,*c,t;

/* This routine sets the values of the switch functions. When the switch
   functions pass through zero from positive to negative the state variables
   may be reset in function map(). The switch functions should pass smoothly
   through 0 and should never have both value and first derivative zero. The
   same switch must not pass through zero from positive to negative more than
   once in an integration timestep. An example of a switch function is:
			   sw[0]=sin(pi*time/30.0)
   which passes through zero every 60 time units. Switches may include state
   variables provided the above conditions are met. Note that to get 'Solver'
   style switches define twice as many switches and let e.g. sw[1]=-sw[0] */

{ sw[0]=50-t;
}

void map(s,c,t,swno)
double *s,*c,t;int swno;

/* This routine is called whenever one of the switch functions passes through
   zero. 'swno' is the number of the switch function. The state variables
   can be changed discontinuously within this routine. eg:
     if (swno==1)
     { state[0]=coeff[1]*(state[0]);}
   time and the coefficients should not be changed.
*/

{ s[0]=s[0]*0.5;
}




void grad(g,s,c,t)
double *g,*s,*c,t;

/* This routine must provide the gradients g for the state variables s.
   So ds[i]/dt=g[i].
   c is the coefficient vector.
   s is the stat variable vector.
   t is time.
   You can access lagged variable that you have stored in storehistory()
   with the functions pastvalue(i,t-T,j) and pastgradient(i,t-T,j).
   The final argument in these calls is the "lagpointer index". You
   can always set this to zero, but if your model includes references to
   the same variable at two different lags, then it is computationally
   more efficient to define one lagpointer for each. So, for example,
   you might call pastvalue(0,t-T1,0) and pastvalue(0,t-T2,1). Use
   of this facility requires that you have specified the number of
   such pointers in initfit().
   Unknown coefficients (i.e. coeffs being fitted) are referenced uc(i)
   Unknown functions are called by:
   uf(i,t) - ith unknown function at any t
   guf(i,t) - gradient of ith unknown function at any t
   iuf(i,t0,t1) - integral from t0 to t1 of ith unknown function
   uc,uf,guf & iuf return doubles (and are implemented as macros)
*/

{ double T,b,d;
  T=uc(0);b=uc(1);d=uc(2);
  if (t>T)
  g[0]=pastvalue(0,t-T,0)*b-d*s[0];
  else
  g[0]= -d*s[0];
}


void storehistory(his,ghis,g,s,c,t)
double *his,*ghis,*g,*s,*c,t;

/* This is the routine in which the values of the history variables at time
   t are calculated and put in the array his, along with gradients in ghis,
   using state variables s gradients of s, g, and coefficients c:
   unknown functions and coeffs can be used in the manner described in
   grad()*/

{ his[0]=s[0];
  ghis[0]=g[0];
}

void statescale(double *scale)

/* In this routine you can set scale factors for error control. For each state
   variable the maximum permisable error will be bounded below by the tolerance
   multiplied by scale[i]. If you don't supply values then zero will be used.
   Non-zero scale values are useful for variables that start at zero and leave
   zero without 3rd order continuity. */

{ scale[0]=1e-5;
}

void initst(s,c,t) double *s,*c,t;

/* initialise state variables and any global constants here, using c or uc*/


{ s[0]=uc(3);
}

/**************************************************************************/
/* Fit specification routines                                             */
/**************************************************************************/

void initfit(d)
ddefit_control *d;

/* This is where you have to supply details about 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.
   You have to fill a control structure for the fit. NOTE if its a double
   it must contain a decimal point or e. If its an integer it must not!
*/

{ int i;

  d->discrete=0;
/***************** The following values must be supplied ******************/
/* The following define numbers of things .......                         */
/**************************************************************************/
  d->no_uc=4;    /* Number of unknown parameters (integer<5000)          */
  d->no_uf=0;    /* Number of unknown functions (integer<=125)           */
  d->no_c=0;     /* Number of (known) model coefficients (integer<5000)  */
  d->no_hv=1;    /* Number of lagged (history) variables (integer)       */
  d->no_s=1;     /* Number of state variables (integer<5000)             */
  d->no_fit=1;   /* Number of variables for which there is data to fit  */
  d->no_sw=1;    /* Number of switches in model (integer)                */
  d->hbsize=1000L;   /* size of history buffer                         */
  d->nlag=1;     /* number of delays per lag variable.(can be left at 1) */
  d->tol=1e-5;   /* relative tolerance used in stepsize control          */
  d->t0=0.0;     /* integration start time                               */
  d->t1=100.0;   /* Default end time, when simulating                    */
  d->dt=1.0;   /* initial time step and 1% of MAXIMUM time step*/
  d->dout=1.0;   /* Output timestep   */
  

/**************************************************************************/
/* 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.                                                   */
/* 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.                                    */
/**************************************************************************/

  d->fitmethod=1;
  d->errors=0;
  d->bsr_reps=0;

/**************************************************************************/
/*  Supply initial values and names for the unknown constants. Any names  */
/* specified beyond the number required will be taken as names of unknown */
/* functions.                                                             */
/**************************************************************************/

 d->uco[0]=4.0;d->cname[0]="Lag - T";
 d->uco[1]=0.3;d->cname[1]="pcbr";
 d->uco[2]=0.28;d->cname[2]="pcdr";
 d->uco[3]=2.0;d->cname[3]="N(0)";



/**************************************************************************/
/* Specify nature of constraints on unknown constants. A constant can be  */
/* UNBOUND, B_BELOW, B_ABOVE or B_BELOW|B_ABOVE. If bounded, you must     */
/* supply bounds here (bounds ignored if type UNBOUND). Default UNBOUND.  */
/**************************************************************************/

 d->uctype[0]=B_BELOW|B_ABOVE;d->uclb[0]=0.001;d->ucub[0]=6.0;
 d->uctype[1]=B_BELOW|B_ABOVE;d->uclb[1]=0.001;d->ucub[1]=1.0;
 d->uctype[2]=B_BELOW|B_ABOVE;d->uclb[2]=0.001;d->ucub[2]=1.0;
 d->uctype[3]=B_BELOW|B_ABOVE;d->uclb[3]=0.001;d->ucub[3]=20.0;

/**************************************************************************/
/* Supply values of known coefficients.                                   */
/**************************************************************************/

/**************************************************************************/
/* 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 INCREASING, DECREASING, B_BELOW, B_ABOVE or        */
/*		combinations of these using |.                            */
/* d->ufstart[i] gives the initial starting value for the u.f.            */
/* 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.                           */
/**************************************************************************/

 /* d->ufdf[0]=7;d->uft0[0]=d->t0;d->uft1[0]=90.0;
  d->uftype[0]=B_BELOW;
  d->uflb[0]=0.0;d->ufub[0]=30.0;
  d->ufsp[0]=-1.0;*/

/**************************************************************************/
/* Specify data and weight file names. The 0th column of the data file    */
/* must be sample time. Then fill out the array specifying which input    */
/* column is which state variable. d->index[i] should give the state      */
/* variable for the ith input column (col 0 is time).                     */
/* Specifify number of rows to read in d0->rows- reads all if <4          */
/* d->cfname is the optional name of a file from which to read initial    */
/* estimates and to which to write final estimates: defaults are used if  */
/* file is empty.                                                         */
/**************************************************************************/

  d->index[1]=0;
  d->dfile="dexp.dat"; /* Data file name */
  d->wfile="dexp.w";   /* Weight file name */
  

/**************************************************************************/
/* You can specify additional statistics of the time series that should   */
/* be fitted. These are STDEV, ABSOLUTEDEV, MEANGRAD, MEANFREQ, NO_STATS  */
/* or ACF. The latter is probably the most useful!                        */
/* They can all be combined using the | bit operator. You get to specify  */
/* what weight they should get. e.g. to give a weight of 1000.0 to the    */
/* square deviation of model from data standard deviation set             */
/* d->wstats[STDEV]=1000.0.                                               */
/**************************************************************************/
  d->statistics=NO_STATS;

/**************************************************************************/
/* Now tell the program how you would like the output to look on the      */
/* screen. Start with the number of windows and the lines in each window. */
/* Next fill in the initial y axis range for each window (d->range[i].yj).*/
/* Then fill in the structure saying where each plotted state variable is */
/* to go. e.g. to make state variable i the curve 2 in  window 1 enter    */
/* d->windex[i].win=1;d->windex[i].cur=1; As usual all indeces start at 0.*/
/**************************************************************************/

  d->no_windows=1;     /* Number of output windows */
  d->lines[0]=1;       /* Number of lines in each window */

  d->windex[0].win=0;d->windex[0].cur=0;

  d->range[0].y0=0.0;d->range[0].y1=150.0;

/**************************************************************************/
/* Now provide labels for the state variables (don't bother with ones that*/
/* will not be plotted) and names for the plotting windows.               */
/**************************************************************************/

  d->label[0]="Population";

  d->wname[0]="Lagged exponential growth model";

}

/*************************************************************************/
/* Routine to initialise unknown functions, for fitting and simulating.  */
/*************************************************************************/



double inuf(int i,double t)

/* Initialises ufs. Called with t as the argument of the ith u.f. should return
   intial value for uf(i,t). */

{
}



void inboot(bsctype *bsc)

/* this is the routine for setting up bootstrapping.
   Leave it blank, or set bsc->reps=0, if you don't want any.
*/

{ 
}