/* Routines for fitting DDE models to observed data. This file contains 
the routines used by the threads that do fitting and numerical solution 
of models. All i/o is via the routines in the GUI thread (se w95dde.c).
 
Copyright (C) 1991-2000 Simon N. Wood  snw@st-and.ac.uk

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License   
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
(www.gnu.org/copyleft/gpl.html)

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
USA. */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h> // used to avoid too frequent data output
#include <sys\timeb.h> // ditto - not completely portable

#define INCL_WIN
#define INCL_DOSPROCESS
#define INCL_WINDIALOGS
#define INCL_DOS
#define STRICT
#include <process.h>
#include <windows.h>
#include <ctype.h> /* white space testing */

#define round(a) ( (a)-floor(a)>0.5 ? (1+(int)floor(a)) : (int) floor(a))

#include "rangen.h"
#include "matrix.h"

#include "ddeq.h"
#include "qp.h"
#include "spline.h"
#include "ddefit95.h"
#include "w95dde.h"
#include "gcv.h"
#include "stochmin.h"

extern qpoutdatatype qpoutdata;
extern short stop_the_fit; // used for signalling that fitting should be halted

double (*splp)(); /* This declares a pointer to a function which is used to
		     assign either the spl() routine or user defined routines
		     to the above macros at run time */

void (*ddep)(); /* this declares a pointer to a function which will be used, either for the 
                   dde() function from ddeq.c, or for the ddem() function from map.c,
                   depending on the type of model to be fitted */



void ErrorMessage(char *msg,int fatal);

void output(double *s,double t,void *dd,int init)

/* This is the routine that updates the output structures. It is called by
   dde() when required (with dd pointing to nothing). If init!=0 then the 
   object in dd is examined for information about the data structure into 
   which output should be directed. This is used by first calling output()
   with init==1, and dd pointing to a display_data_op_type structure 
   (already initialized). Then dde() is called, and will call output() with
   init=0. 
*/



{ static display_data_op_type *d;
  int i;
  if (init) 
  { d=(display_data_op_type *)dd;
    d->n_t=0;
    return;
  }
  for (i=0;i<d->n_disp;i++)  // output data
  { d->display[i][d->n_t]=s[d->distos[i]];}
  d->t[d->n_t]=t;
  d->n_t++; // update output counter
  if (d->n_t>=d->max_n) // expand update buffer if needed
  { for (i=0;i<d->n_disp;i++)
    d->display[i]=(double *)realloc(d->display[i],(size_t)(d->max_n+100)*sizeof(double));
    d->t=(double *)realloc(d->t,(size_t)(d->max_n+100)*sizeof(double));
    d->max_n+=100;
  } 
  
}






void getdifferrors(matrix f,matrix J,matrix a,fit_control_type *fc,double h,double eA,int i,double *df,double *db,
                   double *dc,double *Cf,double *Cb,double *Cp,double *p)


/* finds tedious set of error bounds required to check finite difference
   intervals for acceptability */

{ matrix ff,fb,temp;
  int j;
  ff=initmat(f.r,f.c);fb=initmat(f.r,f.c);temp=initmat(f.r,f.c);
  a.V[i]+=h;F(ff,J,J,a,a,(void *)fc,0,1.0);a.V[i]-=2*h;
  F(fb,J,J,a,a,(void *)fc,0,1.0);a.V[i]+=h;
  for (j=0;j<f.r;j++) temp.V[j]=(ff.V[j]-f.V[j])/h; /*forward diffs */
  *df=enorm(temp)/f.r;*Cf=2*eA/(h*(*df));
  for (j=0;j<f.r;j++) temp.V[j]=(f.V[j]-fb.V[j])/h; /*backward diffs */
  *db=enorm(temp)/f.r;*Cb=2*eA/(h*(*df));
  for (j=0;j<f.r;j++) temp.V[j]=0.5*(ff.V[j]-fb.V[j])/h; /*backward diffs */
  *dc=enorm(temp)/f.r;
  for (j=0;j<f.r;j++)                  /* Curvature */
  temp.V[j]=(ff.V[j]-2*f.V[j]+fb.V[j])/(h*h);
  *p=enorm(temp)/f.r;*Cp=4*eA/(h*h*(*p));
  freemat(ff);freemat(fb);freemat(temp);
}





void Jacobian(int (*F)(matrix,matrix,matrix,matrix,matrix,void*,int,double),
              matrix D,matrix J,matrix a,matrix dav,fit_control_type *fc,matrix f,
                                                                  double tol_factor)

/* This routine finds the partial derivatives of the state variables w.r.t
   the parameters at a series of times the result is a matrix:
	_                                                                 _
       |  ds[0](t[0])/da[0], ds[0](t[0])/da[1], .... ,ds[0](t[0])/da[i]..  |
   J = |  ds[1](t[0])/da[0], ds[1](t[0])/da[1], .... ,ds[1](t[0])/da[i]..  |
       |        .                    .                        .            |
       |_       .                    .                        .           _|
   the routine calls F() which in turn uses the problem specific routines
   defined by the user. The routine also returns a matrix D the first a.r
   columns of which are the fitted values given +ve perturbations, the 2nd a.r
   cols of which are the fitted values for negative pert., and the vector da of
   perturbations -this allows direct FD estimation of the gradient of the
   objective function.
*/

{ matrix ap,fp,fb;
  static matrix dam;
  double da,mean_ce,mean_te;
  long i,j;
  int stop;
  static int first=1;
  if (first)
  { dam=initmat(a.r,1L); /* fractional f.d. intervals */
    for (i=0;i<a.r;i++) dam.V[i]=1e-8;
    first=0;
  }
  ap=initmat(a.r,1L);
  fp=initmat(f.r,1L); /* holds results of forward differences */
  fb=initmat(f.r,1L); /* holds results of backward differences */
  for (i=0;i<a.r;i++) ap.V[i]=a.V[i];
  for (i=0;i<a.r;i++)  /* loop through all parameters */
  { da=dam.V[i]*fabs(a.V[i]); /* the fabs() matters here! */
    if (da<1e-8*dam.V[i]) da=1e-8*dam.V[i]; /* should probably be made proportional to int tol */
    dav.V[i]=da;  /* store the intervals actually used */
    ap.V[i]+=da;  /* forward difference */
    stop=F(fp,J,D,ap,dav,(void *)fc,-1,tol_factor);
    if (stop==1) break; // abort calculation as user specified stop
    ap.V[i]-=da;  /* restore parameter to previous value */
    for (j=0;j<fp.r;j++)
    { D.M[j][i]=fp.V[j];
    }
    { ap.V[i]-=da;  /* backward difference */
      stop=F(fb,J,D,ap,dav,(void *)fc,-1,tol_factor);
      if (stop==1) break;
      ap.V[i]+=da;  /* restore to previous value */
      mean_ce=0.0;mean_te=0.0;
      for (j=0;j<fp.r;j++)
      { mean_te+=0.5*(fp.V[j]-2.0*f.V[j]+fb.V[j])/da; /* estimate of truncation error 0.5*h*f'' */
        mean_ce+=2.0*max(f.V[j],fp.V[j])*1.0e-15/da;  /* estimate of condition error */
	    J.M[j][i]=(fp.V[j]-fb.V[j])/(2.0*da);
	    D.M[j][i+a.r]=fb.V[j];
      }
      if ((dam.V[i]>=1e-10)&&(mean_te>10.0*mean_ce))
      dam.V[i]/=10.0; /* the interval is too large */
      /* limited because too small an interval can cause no change */
      if ((dam.V[i]<=0.001)&&(mean_ce>10.0*mean_te))
      dam.V[i]*=10.0; /* the interval is too small */
      /* but is limited to 1% at most (intervals that are too large can crash the RK scheme */
    } 
  }
  freemat(ap);freemat(fp);freemat(fb);
}




int F(matrix f,matrix J,matrix D,matrix a,matrix da,void *problem_data,int getJ,double tol_factor)


/* this is the routine required by NonLinLS() and optNLLS() prototype in
   qp.h
   getJ should be set to 2 to force output without calculating Jacobian
			 1 if the Jacobian is required
			-1 if being called from Jacobian()
			 0 otherwise
   these codes ensure that finite differencing does not lose accuracy as
   a result of adaptive stepping.

   tol_factor is a number used to multiply the default integration tolerance - usually set to 1
   , it may be set to a smaller (or larger) value as part of the estimation of integration 
   accuracy.

   problem_data is a pointer used for passing detailed problem specific information to this 
   routine. It should immediately be cast to type fit_control_type.

   This routine returns 0 on normal termination and 1 when the user has requested a stop.
*/


{ fit_control_type *fc;
  double dout,*stats,*s,*c,*cp,t0,t1,eps,dt,*pp;
  int j,i,nlag,nsw,hbsize,ns,z,step,nhv,reset,op=0,nt;
  long k,l;
  matrix y,min,max;
  static double ot0=0.0,ot1,tup=0.0,tnow;
  struct timeb tb;
  display_data_op_type *outp;
  status_op_type *outs;
  /* new code */
  fc=(fit_control_type *)problem_data; // fc now contains the detailed information needed by F()
  if (stop_the_fit) return(1); // return without doing anything, because user has requested a stop
  ns=fc->n_s;
  if (fc->trial||getJ>=1) op=1; // then output is required
  s=(double *)calloc((size_t)ns,sizeof(double));
  cp=c=(double *)calloc((size_t)(fc->n_uc+fc->n_c),sizeof(double));
  eps=fc->tol*tol_factor;
  nhv=fc->n_hv;
  nlag=fc->n_lag;
  nsw=fc->n_sw;
  hbsize=(long)fc->hbsize;
  //ioff=20+round(afix.V[7]);
  if (getJ==1)
  step=1;
  else if (getJ==-1) step=2;else step=0;
  /* set up constants */
  c+=fc->n_uc;
  //coff=20+round(afix.V[7])*(1+round(afix.V[5]));
  for (j=0;j<fc->n_c;j++) c[j]=fc->c[j];
  for (j=0;j<fc->n_uc;j++) c[-1-j]=a.V[j];
  /* Setting up splines */
  k=fc->n_uc;
  /* z is offset for uf df */
  //z=(long)((fc->n_fit+1)*fc->n_st+fc->n_c)+20L;
  for (j=0;j<fc->n_uf;j++)   // setting up splines for current run
  { l=fc->ufdf[j];
    if (l<0) l=-l;
    y=initmat((long)l,1L);
    for (i=0;i<l;i++)
    { y.V[i]=a.V[k];k++;}
    spl(j,y,y,fc->t0,0.0,1,0,1);
    freemat(y);
  }
  // Send message to fitting routine that something is happening iff it's more than 1 
  // second since last such message........
  ftime(&tb); tnow=(double)tb.time+tb.millitm*1e-3;
  if (!fc->trial&&tnow-tup>1.0) // tell parent thread is still alive and kicking
  { PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(0,0),(LPARAM)0);
    tup=tnow;
  }
  if (op) // then check that enough time has elapsed since last o/p 
  { if (!fc->trial&&getJ!=2) // it's being called from a fitting routine
    { if (tnow-ot0>0.5) // only output if it's long enough since last time
      { ot0=tnow;
      } else op=0; // cancel output as it's too soon after last output
    } 
  }
  if (op) // then set up for output
  { outp=(display_data_op_type *)calloc(1,sizeof(display_data_op_type));
    splp=audituf; // assigning pointer to audit uf so that uf() usable and function domains recorded
    min=initmat((long)fc->n_uf,1L);max=initmat((long)fc->n_uf,1L); // storage for use by audituf
    audituf(fc->n_uf,min,max,0.0,0.0,1,1,0); // clearing the range (and indicating number of ufs)
    dout=fc->dout; // approx. o/p timestep
    nt=(int)floor((fc->ts[fc->n_st-1]-fc->t0)/dout)+1; // number of output times
    // set up output data structures
    init_display_data_op_type(outp,nt,fc->n_dis,fc->n_uf,fc->distos);
    output(s,0.0,(void *)outp,1);   // initialise function output() for output from dde
  } else 
  dout=0.0; // don't output graphics
  initst(s,c,fc->t0);
  if (fc->t0>fc->ts[0])
  { ErrorMessage("ERROR: Start time after first sample time!",1);}
  t1=fc->t0;dt=fc->dt/*(afix.V[21]-afix.V[20])*0.0005*/;reset=1;
  k=0L; /* f counter */
  for (i=0;i<fc->n_st;i++) /* working through sample times */
  { t0=t1;t1=fc->ts[i];
    if (t1>t0)
    { ddep(s,c,t0,t1,&dt,eps,dout,ns,nsw,nhv,hbsize,nlag,reset,step);reset=0;}
    /* Time to load up vector f (carefully) */
    for (j=0;j<fc->n_fit;j++) /* work through all fit variables */
    { z=i*fc->n_fit+j;
      if (fc->yts[z]>-0.5) { f.V[k]=s[fc->yts[z]];k++;}
    }
  }
  if (op) // output required
  { splp=spl; // re-assign splp to spl instead of audituf
    
    //  post output to the o/p parent window
    audituf(0,min,max,0.0,0.0,1,2,0); // read information out of audituf
    for (i=0;i<fc->n_uf;i++)  // NOTE: could change this so that matrices are not allocated using initmat()
    { outp->ufxmax[i]=max.V[i];
      outp->ufxmin[i]=min.V[i]; 
    }
    PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(3,0),(LPARAM)outp);
    if (!fc->trial)
    { pp=(double *)calloc((size_t)fc->n_p,sizeof(double));
      for (i=0;i<fc->n_p;i++) pp[i]=a.V[i];
      PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(1,0),(LPARAM)pp);  
    } 
    outs=(status_op_type *)calloc((size_t)1,sizeof(status_op_type));
    outs->F=qpoutdata.obj;
    outs->dF=qpoutdata.obj_change;
    outs->cons=qpoutdata.constraints; 
    PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(2,0),(LPARAM)outs);
 
    freemat(max);freemat(min);
  }
  /* now add  values for random effects */
  l=fc->n_uc;
  for (j=0;j<fc->n_uf;j++)   // work through uf's looking for r.e.s
  { //b=afix.V[z+j];   // df's of uf -ve for random effect
    z=fc->ufdf[j];
    if (z<0) // it's a random effect (being negative)
    { z=-z;
      for (i=0;i<z;i++)
      { f.V[k]=a.V[l++];k++;} // load random effect nuisance params straight into f
    }
  }
 


  /* Now add the statistics */
  j=statistics(fc->stats,f,fc,&stats);
  for (i=0;i<j;i++) f.V[f.r-j+i]=stats[i];
  if (j) free(stats);
  free(cp);
  free(s);
  if (getJ==1) Jacobian(F,D,J,a,da,fc,f,tol_factor);
  if (stop_the_fit) return(1);
  return(0); // regular return
}


void G(matrix f,matrix J,matrix D,matrix a,matrix da,void *problem_data,int getJ,
       double tol_factor)


/* this is the routine required by NonLinLS() and optNLLS() prototype in
   qp.h.
   getJ should be set to 2 to force output without calculating Jacobian
			 1 if the Jacobian is required
			-1 if being called from Jacobian()
			 0 otherwise
   these codes are not critical for this version.
   This version is for use when trying to match gradients at sample points,
   predicted by the model on the basis of smoothed raw data, to gradients
   at the sample points obtained from the smooth of the raw data.
   This routine would benefit from some address optimization......
*/


{ fit_control_type *fc;
  double *s,*c,*g,*cp,t1;
  int j,i,ns,z,first_time;
  long k,l;
  matrix y;
  /* new code */
  fc=(fit_control_type *)problem_data;
  ns=fc->n_s;
  s=(double *)calloc(ns,sizeof(double)); /* state variable array */
  g=(double *)calloc(ns,sizeof(double)); /* gradient variable array */
  cp=c=(double *)calloc((size_t)(fc->n_uc+fc->n_c),sizeof(double));
//  ioff=20+round(afix.V[7]);
  /* set up constants */
  c+=fc->n_uc; /* offsetting pointer to c to allow ucs in -ve elements */
  //coff=20+round(afix.V[7])*(1+round(afix.V[5]));
  //soff=coff+round(afix.V[10]+afix.V[9]);
  for (j=0;j<fc->n_c;j++) c[j]=fc->c[j]; /* fixed constants */
  for (j=0;j<fc->n_uc;j++) c[-1-j]=a.V[j];  /* ucs */
  /* Setting up splines representing ufs. */
  k=fc->n_uc;
  //z=(long)(round(afix.V[5]+1)*round(afix.V[7])+round(afix.V[8]))+19L;
  for (j=0;j<fc->n_uf;j++)
  { l=fc->ufdf[j];
    y=initmat((long)l,1L);
    for (i=0;i<l;i++)
    { y.V[i]=a.V[k];k++;}
    spl(j,y,y,fc->t0,0.0,1,0,1);
    freemat(y);
  }
  //if (getJ>=1){ newoutput();}
  k=0L; /* f counter */
  first_time=(int) fc->dt; // NOTE: don't understand this line (while translating to new control structure)
  for (i=first_time;i<fc->n_st;i++) /* working through sample times */
  { t1=fc->ts[i]; /* time at which grad required */
    /* load up the state vector for a call to grad */
    for (j=0;j<ns;j++)
    s[j]= fc->smooth[j][i];//afix.V[soff+i*ns+j];

    /* call grad */

    grad(g,s,c,t1);
// Following commented out for compatibility of re-written code
//    if (getJ>=1) poutput(g,t1); /* put output in window buffers */

    /* unload g into fit vector f */

    for (j=0;j<ns;j++) /* work through all state variables */
    { z=i*ns+j;
      if (fc->yts[z]>-0.5) { f.V[k]=g[fc->yts[z]];k++;}
    }
  }
//  if (getJ>=1) forceoutput(a);
  /* Now add the statistics */
  free(cp);
  free(s);
  free(g);
  if (getJ==1) Jacobian(G,D,J,a,da,fc,f,1.0); 
}


double redirect(no,x,y,x0,x1,mode,resetx,resety)
int no,mode,resetx,resety;matrix x,y;double x0,x1;
/* this routine redirects calls to the macros:
   uf,guf,iuf
   to user defined functions for the purposes of simulating data.
   This is done by setting global pointer (to function) splp to
   redirect, instead of spl.
*/

{ switch (mode)
  { case 0 :; //return(simifunc(no,x0,x1)); /* integral of function */
    case 1 :; //return(simfunc(no,x0));     /* function itself */
    case 2 :; //return(simgfunc(no,x0));    /* gradient of function */
  }
  return(0.0);
}





void profiler(char *name,matrix y,matrix w,matrix f,matrix a,fit_control_type *fc,
              int gradfit,int n,matrix amin,matrix amax)

/* Routine added 20/2/99, to output cross-sections through the un-penalized
   objective to a file. For each parameter  in a the routine steps from amin
   to amax. Parameters with amax<=amin are not output. n steps are taken from
   amin to amax. Output is in pairs of columns - first of pair is
   always the parameter - the second is always the objective. */

{ FILE *fi;
  int i,j,k;
  double aold,z,r;
  matrix da,J,D;
  fi=fopen(name,"wt");
  for (i=0;i<n;i++)
  { for (j=0;j<a.r;j++)
    if (amax.V[j]>amin.V[j])
    { aold=a.V[j];     // store in order to restore after run
      a.V[j]=amin.V[j]+(amax.V[j]-amin.V[j])*i/((double)n-1.0);
      //if (gradfit)
      //G(f,J,D,a,da,afix,2); // NOTE: must be modified when G modified
      //else
      F(f,J,D,a,da,fc,2,1.0);
      r=0.0; for (k=0;k<y.r;k++) { z=y.V[k]-f.V[k];r+=w.V[k]*z*z;}
      r/=2;        // least squares objective
      fprintf(fi,"%8.5g  %8.5g  ",a.V[j],r);
      a.V[j]=aold; // restore original value
    }
    fprintf(fi,"\n");
  }
  fclose(fi);
  //if (gradfit)
  //G(f,J,D,a,da,afix,2); // NOTE: as above
  //else
  F(f,J,D,a,da,fc,2,1.0);
}





double audituf(int no,matrix x,matrix y,double x0,double x1,
               int mode,int resetx,int resety)

/* This function is provided as a means of assessing the range of input values to the 
   model unknown functions. 
   * If it's called with resetx==1 then a new recording is started, and no must contain the 
     number of uf's
   * If it's called with resetx==2 then the max and min x values since last call with
     resetx==1 are returned in matrices min and max.
   * Otherwise max and min values are updated and spl called.
   
   Example of use: 

   splp=audituf;
   min=initmat((long)data.no_uf,1L);max=initmat((long)data.no_uf,1L);
   audituf(no_uf,min,max,0.0,0.0,1,1,0); // clearing the range (and indicating number of ufs)
   
   // code that calls calls the uf's  

   audituf(0,min,max,0.0,0.0,1,2,0); // finding out what the range was

*/

{ static double *low,*hi;
  int i;
  static int initialised=0,size;
  if (resetx==1) /* then it is time to reset the records of uf range limits */
  { if (initialised) {free(low);free(hi);initialised=0;}
    initialised=1;
    low=(double *)calloc((size_t)no,sizeof(double));
    hi=(double *)calloc((size_t)no,sizeof(double));
    for (i=0;i<no;i++) {hi[i]=0.0;low[i]=1.0;}
    size=no;
    return(0.0);
  } else
  if (resetx==2) /* readout requested - this is done into vectors x and y */
  { if ((x.r!=size)||(y.r!=size))
    ErrorMessage("Error in audituf",0);
    for (i=0;i<size;i++)
    { x.V[i]=low[i];
      y.V[i]=hi[i];
    }
  } else /* regular call */
  { if (hi[no]<low[no]) /* first recording call to this uf */
    { if (mode==0) { low[no]=x0;hi[no]=x1;}
      else { low[no]=hi[no]=x0;}
    } else   /* regular call */
    { if (mode==0)      /* integrating */
      { if (x1>hi[no]) hi[no]=x1;}
      else if (x0>hi[no]) hi[no]=x0;
      if (x0<low[no]) low[no]=x0;
    }
    return(spl(no,x,y,x0,x1,mode,0,0));
  }
  return(0.0);
}

void finalrun(matrix y,matrix f,matrix a,fit_control_type *fc,int gradfit)

/* This is a function that does a final model run once the model is fitted.
   Its purpose is to allow monitoring of the final fit - for example by checking
   the ranges of the u.f.s */

{ FILE *fi;
  int i;
  matrix J,D,da,min,max;
  double r2,m,em,e2;
  /* first assign range audit function to splp() */
  splp=audituf;
  min=initmat((long)fc->n_uf,1L);max=initmat((long)fc->n_uf,1L);
  audituf(0,J,J,0.0,0.0,1,1,0);
  if (gradfit) G(f,J,D,a,da,(void *)fc,2,1.0);
  else F(f,J,D,a,da,(void *)fc,2,1.0);
  audituf(0,min,max,0.0,0.0,1,2,0);
  fi=fopen("ddefit.log","at");
  if (fc->n_uf)
  { fprintf(fi,"u.f. input range at best fit\n");
    fprintf(fi,"----------------------------\n");
    fprintf(fi,"Index        Min.          Max.\n");
    for (i=0;i<fc->n_uf;i++)
    { fprintf(fi,"%5d      %8.3g     %8.3g\n",i,min.V[i],max.V[i]);}
  }
  m=0.0;em=0.0;
  for (i=0;i<y.r;i++) { m+=y.V[i];em+=(y.V[i]-f.V[i]);} m/=y.r;em/=y.r;
  r2=0.0;e2=0.0;
  for (i=0;i<y.r;i++)
  { r2+=(y.V[i]-m)*(y.V[i]-m);
    e2+=(y.V[i]-f.V[i]-em)*(y.V[i]-f.V[i]-em);
  }
  r2=1.0-e2/r2;
  fprintf(fi,"\n r^2 = %g\n",r2);
  fclose(fi);
  freemat(min);freemat(max);
}






double *setupcoeffs(ddefit_control *d,int create)
/* setup for coeffs */

{ static double *c;
  double *cc;
  int i;
  if (create)
  { c=(double *)calloc((size_t)d->no_c+d->no_uc+1,sizeof(double));
    cc=c+d->no_uc;
    for (i=0;i<d->no_uc;i++) cc[-i-1]=d->uco[i];
    for (i=0;i<d->no_c;i++) cc[i]=d->c[i];
    return(cc);
  }
  free(c);
  return (double *)NULL;
}
#undef isspace

int IsCharNumber(char c)

{ if (c=='-') return(1);
  if (c=='+') return(1);
  if (c=='e') return(1);
  if (c=='E') return(1);
  if (c=='0') return(1);
  if (c=='1') return(1);
  if (c=='2') return(1);
  if (c=='3') return(1);
  if (c=='4') return(1);
  if (c=='5') return(1);
  if (c=='6') return(1);
  if (c=='7') return(1);
  if (c=='8') return(1);
  if (c=='9') return(1);
  if (c=='.') return(1);
  return(0); 
}

int getno(FILE *f,double *x)

/* returns 0 for number, 
           1 missing value, 
           2 for newline & number, 
           3 for newline & missing value
           5 for eof or newline with no number
 
   - Modified 24/5/2000 - previous version did something very stupid with 
     static variables that I would rather not admit to, but meant that you 
     could get the last character of the previous file returned as the first 
     number from the current file (although this did happily cause a flagged 
     error in the routines that used this).

*/

{ int i,started,ok,newline=0,missing;
  char c,st[200];
  started=0;ok=1;i=0;missing=0;
  while (ok) 
  { if (feof(f)) { c=' ';ok=0;} else c=(char)fgetc(f);
    if (started)
    { if (IsCharNumber(c)) st[i++]=c; // add character to number string
      else // other character terminate number string
      { ok=0;st[i]=(char)NULL;
        if (c=='\n'||c=='\r') newline=1;
      } 
    } else // number not started yet
    { if (IsCharNumber(c)) { started=1;st[i++]=c;} // add character to number string
      else 
      if (c=='\n'||c=='\r') { return(5);} // end of line and no new number
      else 
      if (c=='*')  // missing value
      { missing=1;ok=0;i++;}
    } 
  }
  if (!i) return(5); // eof terminated 
  else  // read forward to next number or newline or missing
  { ok=1;
    while (ok)
    { if (feof(f)) {ok=0;newline=1;} else c=(char)fgetc(f);
      if (IsCharNumber(c)||c=='*') { ok=0;ungetc(c,f);}
      else if (c=='\n'||c=='\r') {newline=1;ok=0;}
    } 
  } 
  if (missing) if (newline) return(3); else return(1);
  // At this stage there is a supposed number in st
  sscanf(st,"%lf",x);
  if (newline) return(2); // number and newline
  else return(0); // number, no newline
  
}

#define ACFLAG 10

int countstats(int statype)

/* returns number of additional statistics being used */

{ int i=0;
  if (statype&STDEV) i++; if (statype&ABSOLUTEDEV) i++;
  if (statype&MEANGRAD) i++; if (statype&MEANFREQ) i++;
  if (statype&ACF) i+=ACFLAG;
  return(i);
}

double rmsgrad(matrix s,matrix t)

{ int i;
  double m=0.0,z;
  for (i=0;i<t.r-1;i++)
  z=fabs(s.V[i+1]-s.V[i])/(t.V[i+1]-t.V[i]);
  z=z*z;
  m+=z;
  m/=t.r-1;
  return(sqrt(m));
}


double SFTPower(t,h,p) matrix t,h,*p;
/* Does slow fourier transform on data in h and t, puts power in p
and returns frequency interval */

{ static double twopi, Hr,Hi;
  static int first=1;
  double *hdt,f,df;
  long i,k,n;
  if (first)
  { first=0;twopi=asin(1.0)*4.0;}
  k=t.r/2+1;
  (*p)=initmat(k,1L);
  hdt=(double *)malloc(t.r*sizeof(double));
  hdt[0]=h.V[0]*(t.V[1]-t.V[0]);
  for (i=1;i<t.r-1;i++) hdt[i]=h.V[i]*(t.V[i+1]-t.V[i-1]);
  hdt[t.r-1]=h.V[t.r-1]*(t.V[t.r-1]-t.V[t.r-2]);
  df=(t.r-1.0)/(2.0*(k-1)*(t.V[t.r-1]-t.V[0]));
  for (n=0;n<k;n++)
  { f=n*df*twopi;
    Hr=hdt[0]*cos(f*t.V[0]);
    Hi=hdt[0]*sin(f*t.V[0]);
    for (i=1;i<t.r-1;i++)
    { Hr+=hdt[i]*cos(f*t.V[i]);
      Hi+=hdt[i]*sin(f*t.V[i]);
    }
    Hr+=hdt[t.r-1]*cos(f*t.V[t.r-1]);
    Hi+=hdt[t.r-1]*sin(f*t.V[t.r-1]);
    Hr*=0.5;Hi*=0.5;
    p->V[n]=Hr*Hr+Hi*Hi;
  }
  free(hdt);
  return(df);
}

double meanf(y,t)matrix y,t;

/* returns the average frequency of data y taken at points t. This is done by
   a slow fourier transform - which makes dealing with unevenly spaced data
   with non power of 2 elements easy, but expensive. */

{ int i;
  matrix p;
  double df,mean=0.0, mp=0.0,w;
  df=SFTPower(t,y,&p);
  for (i=0;i<p.r;i++)
  { w=p.V[i];mean+=i*w;mp+=w;
  }
  freemat(p);
  return(df*mean/mp);
}


void testSFT()

{ matrix t,y;
  int i;
  double pi;
  pi=4*asin(1.0);
  t=initmat(100L,1L);y=initmat(t.r,1L);
  for (i=0;i<t.r;i++)
  { t.V[i]=i*0.1;
    y.V[i]=sin(pi*t.V[i]*0.2);
  }
  pi=meanf(y,t);
  freemat(t);freemat(y);
}



int statistics(int statype,matrix y,fit_control_type *fc,double **res)

/* this routine calculates statistics for the vector y, which is disaggregated
   using afix into data relating to separate state variables. The stats for each
   state variable are summed - although this could be altered
   NOTE that res is allocated memory in this routine, that must be freed
   elsewhere */

{ long k,*kk;
  int i,j,z,q;
  matrix t,*s;
  if (statype==NO_STATS) return(0);
  t.V=fc->ts;t.c=1L;t.r=(long)fc->n_st;
  s=(matrix *)calloc((size_t)fc->n_fit,sizeof(matrix));
  kk=(long *)calloc((size_t)fc->n_fit,sizeof(long));
  for (i=0;i<fc->n_fit;i++) s[i]=initmat(t.r,1L);
  k=0L; /* Next apportion the data betwixt the s vectors */
  for (i=0;i<fc->n_st;i++) /* working through sample times */
  for (j=0;j<fc->n_fit;j++) /* work through all fit variables */
  { z=i*fc->n_fit+j;
    if (fc->yts[z]>-0.5)
    {s[j].V[kk[j]]=y.V[k];kk[j]++;k++;} else s[j].r--;
  }
  free(kk);
  /* count the statistics */
  z=countstats(statype);
  i=fc->n_fit;
  (*res)=(double *)calloc((size_t)(z*i),sizeof(double)); /*location for results*/

  k=0L;
  /* work through statistics */
  if (statype&STDEV)
  { for (j=0;j<i;j++) {(*res)[k]+=sqrt(variance(s[j]));k++;}}
  if (statype&ABSOLUTEDEV)
  { for (j=0;j<i;j++) {(*res)[k]+=absdev(s[j]);k++;}}
  if (statype&MEANGRAD)
  { for (j=0;j<i;j++) {(*res)[k]+=rmsgrad(s[j],t);k++;}}
  if (statype&MEANFREQ)
  { for (j=0;j<i;j++) {(*res)[k]+=meanf(s[j],t);k++;}}
  if (statype&ACF)
  { for (j=0;j<i;j++) // stage loop
    for (q=0;q<ACFLAG;q++) {(*res)[k]=acf(s[j],q+1);k++;} // lag loop
  }
  /*clear up*/
  for (j=0;j<fc->n_fit;j++) freemat(s[j]);
  free(s);
  return(z*i);  /* returns the number of statistics */
}

void statweight(fit_control_type *fc,double **ws)
/* routine to sort out the weight to be given to different statistics
   that are to be fitted, assumes the order of working through stats,
   used in statistics() */
{ int i=0,j,z;
  (*ws)=(double *)calloc((size_t)(countstats(fc->stats)*fc->n_fit),sizeof(double));
  if (fc->stats&STDEV)
  for (j=0;j<fc->n_fit;j++) { (*ws)[i]=fc->wstats[STDEV]; i++;}
  if (fc->stats&ABSOLUTEDEV)
  for (j=0;j<fc->n_fit;j++) { (*ws)[i]=fc->wstats[ABSOLUTEDEV];i++;}
  if (fc->stats&MEANGRAD)
  for (j=0;j<fc->n_fit;j++) { (*ws)[i]=fc->wstats[MEANGRAD];i++;}
  if (fc->stats&MEANFREQ)
  for (j=0;j<fc->n_fit;j++) { (*ws)[i]=fc->wstats[MEANFREQ];i++;}
  if (fc->stats&ACF)
  { for (j=0;j<fc->n_fit;j++) // stage loop
    for (z=0;z<ACFLAG;z++) {(*ws)[i]=fc->wstats[ACF];i++;} // lag loop
  }
}


double datastore(int n,int stages,double *t,int *index,double *y,int state,double time)

/* if n!=0 stores data state by state for later access from the model
   specification routines */

{ static double *ts,*ys;
  static int i,j,k,ns,nt;
  if (n)
  { ns=stages;nt=n;
    ys=(double *)calloc((size_t)n*ns,sizeof(double));
    ts=(double *)calloc((size_t)n,sizeof(double));
    k=0;
    for (i=0;i<n;i++)
    { ts[i]=t[i];
      for (j=0;j<ns;j++)
      { if (index[i*ns+j]!= -1) ys[i*ns+j]=y[k++];
        else ys[i*ns+j]=-1.0;
      }
    }
  } else
  { if (time<ts[0]) return(-1.0);
    k=0; while ((time!=ts[k])&&(k<nt)) k++;
    if (k==nt)
    ErrorMessage("Over-run time in datastore!",0);
    i=k*ns+state;
    return(ys[i]);
  }
}

double obs(int i,double s,double t)

/* gives map model routines access to observed state variables: The user gets
   to know which fit variables are which state variables, s must contain the
   corresponding state variable, in case the observed is missing */

{ double r;
  int dum;
  r=datastore(0,0,&s,&dum,&s,i,t);
  if (r==-1.0) return(s); else return(r);
}


void readin(matrix *y,matrix *w,fit_control_type *fc)

/* This routine reads data into y and weights into w.
   * Data are read from fc->dname and weights from fc->wname.
   * fc->index[i] contains the state variable index for data in column i of
     the data file (where col 0 is time).
   * fc->n_st contains the number of rows of the data file to read in
     - if it's 0 then all will be read in. fc->n_st gets reset to number of rows 
     actually read on exit.
   * fc->raw and fc->rawcount are also set up to contain x,y data for each fitted
     series and number of data in each fitted series respectively. 
   * NOTE: not sure about start time initialization yet 
   
   * Missing values in the input data file are ok: but not missing sample times 
     (first column of data file), or missing weights. The missing value character 
     is '*' (without the quotes!).
   
   * Here are the elements of fc to be initialized on input:
     n_st, dfile, wfile, n_fit, n_uf, nfdf, uftype, stats, index, re_vmax, wstats
   
   * Here are the elements that are initialized in here:
     raw, raw_n, yts, n_st (reset), ts, n
 
*/

{ double x,yy,*stats,*ws;
  int i,k,j=0,rows,cols=0,l;
  FILE *f,*fw;
  f=fopen(fc->dfile,"rt");
  if (!f) 
  ErrorMessage("Input data file not found",1); // NOTE: eventually prompt for this if not found
  i=0;rows=0;
  while (!feof(f)) // counting the 
  { k=getno(f,&x);
    if (k%2==0) i++;   /* A number rather than a missing value */
    if (k!=5) // then it's not a lone newline
    { if (!rows) cols++; /* counting columns */
      else j++;
    } 
    if ((k/2>0)||(feof(f)))         /* there was a newline char after last no. */
    { if (rows&&(j!=cols)&&(j!=0))
      { ErrorMessage("\nInput file has unequal row lengths\n",1);}
      if (!rows||j>0) rows++;j=0;
    }
  }
  if (fc->n_st>3 && fc->n_st<rows) rows=fc->n_st; // how many rows to read??
  fc->n_st=rows;  // rows actually read
  if (fc->n_fit!=cols-1) ErrorMessage("Number of cols doesn't match number of fit variables",1);
  fc->raw=(xytype **)calloc((size_t)fc->n_fit,sizeof(xytype *));
  for (j=0;j<fc->n_fit;j++)
  fc->raw[j]=(xytype *)calloc((size_t)rows,sizeof(xytype));/* temp storage for raw data */
  fc->raw_n=(int *)calloc((size_t)fc->n_fit,sizeof(int));  /* counters for raw data */
  fseek(f,0,0);
  /* setting up y - three elements to it:
     1. the actual data (always present)
     2. the expected values of the random effects (optional)
     3. the statistics of the data (optional)
  */
  k=0;
  for (l=0;l<fc->n_uf;l++) /* count up random effect terms */
  { j=fc->uftype[l];if (j&RE1) k+=abs(fc->ufdf[l]);}
  k += countstats(fc->stats)*fc->n_fit; // extra statistics to be fitted
  (*y)=initmat((long)(i-rows+k),1L);
  fc->ts=(double *)calloc((size_t)fc->n_st,sizeof(double)); // space for sample times
  //(*afix)=initmat((long)(rows*cols+20+dfc->no_uf+dfc->no_c),1L);
  (*w)=initmat(y->r,1L);
  fc->n=y->r; // number of data to be fitted
  fw=fopen(fc->wfile,"rt");
  if (!fw) ErrorMessage("Weight file not found, weights set to unity",0);
  fc->yts=(int *)calloc((size_t)rows*(cols-1),sizeof(int));
  l=0;
  for (i=0;i<rows;i++) /* read data into Y */
  { k=getno(f,&x);fc->ts[i]=x; /* sample times */
    for (j=0;j<cols-1;j++)
    { k=getno(f,&yy);
      if (k%2==0) // then data is not missing
      { y->V[l++]=yy; // loading the data vector
        fc->raw[j][fc->raw_n[j]].x=(float)x; // load data for display
	    fc->raw[j][fc->raw_n[j]].y=(float)yy;
        fc->raw_n[j]++; // counting data
      } 
      // now fill in array used for unpacking state variables into fitted value vector
      // in routine F() - see above.
      if (k%2) fc->yts[i*(cols-1)+j]=-1;else  /* -1 for no data */
      fc->yts[i*(cols-1)+j]=fc->index[j+1];     /* state variable no. */
    }
  }
  /* Now append the expected values of any random effects to the data vector */
  //matrixintegritycheck();
  for (i=0;i<fc->n_uf;i++)
  { k=fc->uftype[i];
    if (k&RE1)
    for (j=0;j<fc->ufdf[i];j++) y->V[l++]=0.0;
  }
  l=0;
  for (i=0;i<rows;i++)
  { for (j=0;j<cols-1;j++)
    { if (fw) getno(fw,&x);else x=1.0;
      if (fc->yts[i*(cols-1)+j]>=0) { w->V[l]=x;l++;}
    }
  }
  if (fw) fclose(fw);if (f) fclose(f);
  /* weights for any random effects */

  for (i=0;i<fc->n_uf;i++)
  { k=fc->uftype[i];
    if (k&RE1)
    for (j=0;j<fc->ufdf[i];j++) w->V[l++]=1.0/fc->re_vmax[i];
  }

  datastore(rows,cols-1,fc->ts,fc->yts,y->V,0,0.0);

  /* setting up rest of elements in afix 13 through 19 currently free
   NOTE: this is left here for reference about what was in afix!!!
   only remove when all references to afix have been purged (hooray).
   */
 // afix->V[0]=dfc->tol;
 // afix->V[1]=(double) dfc->no_hv;
 // afix->V[2]=(double) dfc->no_s;
 // afix->V[3]=(double) dfc->nlag;
 // afix->V[4]=(double) dfc->no_sw;
 // afix->V[5]=(double) dfc->no_fit; /* cols-1 */
 // afix->V[6]=dfc->t0;
 // afix->V[7]=(double) rows;
 // afix->V[8]=(double) dfc->no_uc;
 // afix->V[9]=(double) dfc->no_uf;
 // afix->V[10]=(double)dfc->no_c;
 // afix->V[11]=(double)dfc->hbsize;
 // afix->V[12]=(double)dfc->statistics;
 // afix->V[13]=dfc->dt;
 // dfc->int0=dfc->t0;      /* start of model integration */
 // dfc->t0=afix->V[20];    /* start of data */
 // dfc->t1=afix->V[19+round(afix->V[7])];
 // for (i=0;i<dfc->no_c;i++)
/*  afix->V[afix->r-dfc->no_uf-dfc->no_c+i]=dfc->c[i];*/
//  afix->V[20+(dfc->no_fit+1)*rows+i]=dfc->c[i];
  for (i=0;i<fc->n_uf;i++)
  /* afix->V[afix->r-dfc->no_uf+i]=dfc->ufdf[i];*/
  { k=fc->uftype[i];
    if (k&RE1) fc->ufdf[i]= -abs(fc->ufdf[i]); // signals random effect - not convinced this should be here
  }
  i=statistics(fc->stats,*y,fc,&stats);
  statweight(fc,&ws);
  for (j=0;j<i;j++)
  { y->V[y->r-i+j]=stats[j];
    w->V[w->r-i+j]=ws[j];
  }
  if (i) { free(stats);free(ws);}
}

void taper(afix,w,w1,c,leader) matrix afix,w,w1;double c;long leader;

/* Takes an initial vector of weights w and multiplies those that relate
   to observations of time dependent state variables, by a factor
   exp(-c*t) where t is time of the observation after sample 'leader' and c
   is a decay constant. weights for the first 'leader' sample times are
   unaltered. The weights given to statistics
   of the time series are unaffected. The new weights are returned in w1
   which the user must initialise to be the same length as w. */

{ long rows,cols,i,j,l;
  rows=(long)round(afix.V[7]);
  cols=(long)round(afix.V[5]);
  l=0L;
  for (i=0;i<rows;i++)
  { for (j=0;j<cols;j++)
    { if (round(afix.V[20+rows+i*cols+j])>=0)
      { if (i<leader) w1.V[l]=w.V[l]; else
	     w1.V[l]=w.V[l]*exp(-c*(afix.V[20+i]-afix.V[20+leader]));
	     l++;
      }
    }
  }
  for (i=l;i<w.r;i++) w1.V[i]=w.V[i]; /* statistics unchanged */
}

void setupunknowns(fit_control_type *fc,matrix *S,int *off,matrix *A,matrix *Af,
                   matrix *b,matrix *a)

/* Sets up overall constraint and penalty matrices. Note that individual
   smoothness matrices are generated here, and are not y scaled at present,
   before being added (although they are x scaled).
   
   S and off are arrays - memory must be initialized before call, but matrices
   initialised in this routine.

   A, Af, b and a are initialized in this routine, and must be freed later
   with freemat();

   Don't forget that unknown parameter vector starts with unknown coeffs and
   then moves on to uf's.
*/

{ int i,j,k,r,c,offi,up;
  long coff;
  matrix D,x,*Ai,*bi,Ac,bc,temp,temp1,bt,y,h;
  k=0;for (i=0;i<fc->n_uf;i++) k+=fc->ufdf[i];
  k+=fc->n_uc; /* total no. params. to be fitted */
  (*a)=initmat((long)k,1L);
  Ai=(matrix *)calloc((long)fc->n_uf,sizeof(matrix));
  bi=(matrix *)calloc((long)fc->n_uf,sizeof(matrix));
  *Af=initmat(k,k); /* the fixed constraint matrix */
  A->r=0L;offi=fc->n_uc;Af->r=0L;
  for (i=0;i<fc->n_uc;i++) a->V[i]=fc->p[i];
 
  for (i=0;i<fc->n_uf;i++)
  { off[i]=offi;
    //dt=(dfc->uft1[i]-dfc->uft0[i])/(dfc->ufdf[i]-1);
    x=initmat((long)fc->ufdf[i],1L);
    y=initmat(x.r,1L);
    for (j=0;j<x.r;j++)
    { x.V[j]=fc->ufx[i][j];  
      a->V[offi+j]=y.V[j]=fc->p[offi+j];
    }
    spl(i,x,y,x.V[0],0.0,1,1,1);
    k=fc->uftype[i];
    if (k&RE1)  /* then this u.f. is a random effect and has identity matrix penalty */
    { S[i]=initmat(x.r,x.r);
      for (j=0;j<x.r;j++) S[i].M[j][j]=1.0;
      //WarningQuery("Penalties for r.e. in setupunknowns()");
    }
    else getSmooth(S+i,x,1);  /* regular spline based uf */
//    D=initmat(x.r,x.r); -appeared to be here for no good reason

    if (k&ORIGIN) /* the first knot of the function is to be fixed */
    { Af->M[Af->r][offi]=1.0;
      Af->r++;
    }
    offi+=S[i].r;
    if ((k&INCREASING)||(k&DECREASING))   /* monotonic */
    { if (k&INCREASING) up=1;else up=0;
      h=initmat(x.r-1,1L);
      for (j=0;j<h.r;j++) h.V[j]=x.V[j+1]-x.V[j];
      D=getD(h,0);
      createA(D,&temp,h,up,0);
      freemat(D);freemat(h);
      j=temp.r;
      if (k&B_BELOW) j++;if (k&B_ABOVE) j++;
      Ai[i]=initmat(j,temp.c);
      bi[i]=initmat(j,1L);
      for (r=0;r<temp.r;r++) for (c=0;c<temp.c;c++)
      Ai[i].M[r][c]=temp.M[r][c];
      if (k&B_BELOW)
      if (k&INCREASING)
      { Ai[i].M[temp.r][0]=1.0;bi[i].V[temp.r]=fc->uflb[i];}
      else { Ai[i].M[temp.r][temp.c-1]=1.0;bi[i].V[temp.r]=fc->uflb[i];}
      if (k&B_ABOVE)
	   if (k&INCREASING)
      { Ai[i].M[Ai[i].r-1][temp.c-1]= -1.0;bi[i].V[Ai[i].r-1]= -fc->ufub[i];}
      else { Ai[i].M[Ai[i].r-1][0]= -1.0;bi[i].V[Ai[i].r-1]= -fc->ufub[i];}
      freemat(temp);
    } else
    { temp.r=temp1.r=0L;
      if (k&B_BELOW) { poscon(&temp,&bt,x,5L,1);freemat(bt);}
      if (k&B_ABOVE) { poscon(&temp1,&bt,x,5L,1);freemat(bt);}
      Ai[i]=initmat(temp.r+temp1.r,x.r);
      bi[i]=initmat(temp.r+temp1.r,1L);
      if (k&B_BELOW)
      for (r=0;r<temp.r;r++)
      { bi[i].V[r]=fc->uflb[i];
        for (c=0;c<x.r;c++) Ai[i].M[r][c]=temp.M[r][c];
      }
      if (k&B_ABOVE)
      for (r=0;r<temp1.r;r++)
      { bi[i].V[r+temp.r]=-fc->ufub[i];
        for (c=0;c<x.r;c++) Ai[i].M[r+temp.r][c]=-temp1.M[r][c];
      }
      if (k&B_BELOW) freemat(temp);
      if (k&B_ABOVE) freemat(temp1);
    }
    freemat(x);freemat(y);
    A->r+=Ai[i].r;
  }
  /* now deal with constraints on coeffs */
  Ac=initmat(2*fc->n_uc,a->r);
  bc=initmat(2*fc->n_uc,1L);
  k=0;
  for (i=0;i<fc->n_uc;i++)
  { if (fc->uctype[i]&FIXED) // coefficient is fixed at its initial value
    { Af->M[Af->r][i]=1.0; Af->r++; // updating fixed constraint matrix
    } else                    // it might be bounded ....
    { if (fc->uctype[i]&B_BELOW)
      { Ac.M[k][i]=1.0;bc.V[k]=fc->uclb[i];k++;}
      if (fc->uctype[i]&B_ABOVE)
      { Ac.M[k][i]= -1.0;bc.V[k]= -fc->ucub[i];k++;}
    }
  }
  /* Now collect together the matrices (and free the local ones)*/
  (*A)=initmat(A->r+k,a->r);
  (*b)=initmat(A->r,1L);
  offi=0;coff=(long)fc->n_uc;
  for (i=0;i<fc->n_uf;i++)
  { for (r=0;r<Ai[i].r;r++)
    { b->V[offi+r]=bi[i].V[r];
      for (c=0;c<Ai[i].c;c++)
      A->M[offi+r][coff+c]=Ai[i].M[r][c];
    }
    offi+=Ai[i].r;coff+=Ai[i].c;
    freemat(Ai[i]);freemat(bi[i]);
  }
  for (r=0;r<k;r++)
  { b->V[offi+r]=bc.V[r];
    for (c=0;c<Ac.c;c++)
    A->M[offi+r][c]=Ac.M[r][c];
  }
  freemat(Ac);freemat(bc);
}

void ufsetup(double *ufx,double *ufy,double x0,double x1,int nk,int i)

/* Routine for setting up knot values ufx[] and user defined initial values at the knots 
   ufy[] of the ith spline, which is assumed to have nk knots. x0 is lowest knot position
   while x1 is highest. */

{ int j;
  double dx;
  dx=(x1-x0)/(nk-1);
  for (j=0;j<nk;j++)
  { ufx[j]=x0+j*dx;ufy[j]=inuf(i,ufx[j]);
  }
}

double pvspliner(no_his,no_uf,x,y,t,mode,reset)
int no_his,no_uf,mode,reset;matrix x,y;double t;

/* this is the interface routine for calling pastvalues of spline
   smoothed data via pastvalue() ( the reason for setting this up is so that
   models do not need to be re-specified in order to use fitting to smoothed
   gradients) The routine must be called in set up mode for each lagged
   variable in turn. no_uf is the number of unknown functions used in the
   modelling (these get stored first). mode tells the routine whether to return
   gradient (2) or value (1)*/

{ static int total_uf;
  if (reset) /* a new smooth is being stored */
  { total_uf=no_uf;
    return(spl(no_his+no_uf,x,y,t,0.0,mode,1,1));
  } else
  return(spl(no_his+total_uf,x,y,t,0.0,mode,0,0));
}

double pvscaller(no,t,markno) int no,markno;double t;

/* this routine is a dummy for assigning to ppastvalue(), it calls
   the spline routines substituted for pastvalue routines when fitting
   to smoothed gradients */

{ matrix dum;
  return(pvspliner(no,0,dum,dum,t,1,0));
}

double pgcaller(no,t,markno) int no,markno;double t;

{ matrix dum;
  return(pvspliner(no,0,dum,dum,t,2,0));
}

void gethindex(dfc,hindex) ddefit_control *dfc;int *hindex;

/* gets the array whose ith element is the data variable correspoding to
   the ith history variable - this is achieved using dfc->index and by
   interrogation of storehistory() */

{ double *s,*h,*gh,*g,*c,t=0.0;
  int *hsmap,i,j; /* hsmap[i] is state variable for ith history variable */
  s=(double *)calloc((size_t)dfc->no_s,sizeof(double));
  h=(double *)calloc((size_t)dfc->no_hv,sizeof(double));
  c=(double *)calloc((size_t)dfc->no_c,sizeof(double));
  g=(double *)calloc((size_t)dfc->no_s,sizeof(double));
  gh=(double *)calloc((size_t)dfc->no_hv,sizeof(double));
  hsmap=(int *)calloc((size_t)dfc->no_hv,sizeof(int));
  for (j=0;j<dfc->no_hv;j++) hsmap[j]=-1;
  for (i=0;i<dfc->no_s;i++)
  { s[i]=1000.0;      /* probe storehistory with a known vector */
    storehistory(h,gh,g,s,c,t);
    for (j=0;j<dfc->no_hv;j++)
    if (h[j]==1000.0) hsmap[j]=i;
    s[i]=0.0;
  }
  for (j=0;j<dfc->no_hv;j++)
  if (h[j]==-1)
  ErrorMessage("A history variable is missing or is not a simple state vector, can't gradient fit.",1);
  for (j=0;j<dfc->no_hv;j++)
  { i=1;while((hsmap[j]!=dfc->index[i])&&(i<dfc->no_fit+1)) i++;
    if (i==dfc->no_fit+1)
    ErrorMessage("Can not use gradient fit with unobserved state variables.",1);
    hindex[j]=i-1; /* the -1 is 'cos of column notation */
  }
  free(s);free(h);free(c);free(g);free(gh);free(hsmap);
}



void setupgradfit(dfc,Y,W,Co,afix)
ddefit_control *dfc;matrix *Y,*W,*afix,*Co;

/* horrible routine for setting up gradient fit in user transparent manner
   (user transparancy == coding opacity)
   Returns the gradients to be fitted in Y and the corresponding Covariance
   matrix in Co.
*/

{ int *hindex; /* gives data variable of ith history variable */
  long k,*kk,**Z;
  int i,j,z,ioff,soff,hi;
  matrix D,h,t,*s,*b,*x,*w,a,c,d,tempM,A,*Cog;
  double lam,dd,sig2,df;
  hindex=(int *)calloc((size_t)dfc->no_hv,sizeof(int));
  gethindex(dfc,hindex);
  /* set pastvalue functions to point to splines, not lags */
  ppastvalue=pvscaller;
  ppastgradient=pgcaller;

  /* now find first sample time */
  i=0;while (afix->V[20+i]-afix->V[20]<dfc->biglag) i++;
  afix->V[13]=(double) i;

  /* sort out data for smoothing, into vectors s[i] */

  t.V=afix->V+20;t.c=1L;t.r=(long)round(afix->V[7]);
  s=(matrix *)calloc((size_t)round(afix->V[5]),sizeof(matrix));
  w=(matrix *)calloc((size_t)round(afix->V[5]),sizeof(matrix));
  x=(matrix *)calloc((size_t)round(afix->V[5]),sizeof(matrix));
  b=(matrix *)calloc((size_t)round(afix->V[5]),sizeof(matrix));
  Cog=(matrix *)calloc((size_t)round(afix->V[5]),sizeof(matrix));
  kk=(long *)calloc((size_t)round(afix->V[5]),sizeof(long));
  for (i=0;i<round(afix->V[5]);i++)
  { s[i]=initmat(t.r,1L);
    w[i]=initmat(t.r,1L);
    x[i]=initmat(t.r,1L);
    Cog[i]=initmat(t.r,t.r);
  }
  ioff=20+(int)round(afix->V[7]);
  k=0L; /* Next apportion the data betwixt the s vectors */
  for (i=0;i<round(afix->V[7]);i++) /* working through sample times */
  for (j=0;j<round(afix->V[5]);j++) /* work through all state variables */
  { z=ioff+i*(int) round(afix->V[5])+j;
    if (afix->V[z]>-0.5)
    { s[j].V[kk[j]]=Y->V[k];
      w[j].V[kk[j]]=W->V[k];
      x[j].V[kk[j]]=t.V[i];
      kk[j]++;k++;
    } else { s[j].r--;w[j].r--;x[j].r--;}
  }
  for (i=0;i<round(afix->V[5]);i++)
  b[i]=initmat(x[i].r,1L);
  free(kk);

  /* now do state by state smoothing */

  tempM=(*afix);
  (*afix)=initmat(tempM.r+round(afix->V[7]*afix->V[5]),1L);
  for (i=0;i<tempM.r;i++) afix->V[i]=tempM.V[i];
  freemat(tempM);t.V=afix->V+20;
  for (i=0;i<dfc->no_s;i++) /* work through all data variables */
  { a=initmat(x[i].r,1L);c=initmat(a.r,1L);d=initmat(a.r,1L);
    lam=ss_fit(x[i],s[i],w[i],a,b[i],c,d,2,0); /* 2,0 => gcv all coeffs returned */
    /* debug code to check GetInf matches ss_fit */
    A=initmat(x[i].r,x[i].r);
    GetInf(x[i],w[i],A,1.0/lam); /* Influence matrix seems ok - but there is
				     small mismatch ss_fit to GetInf + big coeff
				     problem - needs fixing with ref to pascal*/
    /* AW^{-1} sig2 is covariance of fit */
    sig2=0.0;df=0.0;
    for (j=0;j<s[i].r;j++)
    { sig2+=w[i].V[j]*(s[i].V[j]-a.V[j])*(s[i].V[j]-a.V[j]);
      df+=(1.0-A.M[j][j]);
    }
    sig2/=df;
    for (j=0;j<A.r;j++)
    for (k=0;k<A.r;k++) A.M[j][k]*= sig2/w[i].V[j];
    h=initmat(x[i].r-1,1L);
    for (j=0;j<x[i].r-1;j++) h.V[j]=x[i].V[j+1]-x[i].V[j];
    D=getD(h,0);
    freemat(h);
    Cog[i]=initmat(A.r,A.r);
    symproduct(D,A,Cog[i],0,0); /* gets cov matrix for gradients */
    freemat(D);freemat(A);
    /* Now it gets nasty - the covariance matrix for a stage has
       to be incorporated in the overall covariance matrix for all
       gradients (after the discarded data segment */

    /* covariance matrix ought to be returned here */
    /* work through all sample times getting smoothed state variable */
    dd=20.0+afix->V[7]*(afix->V[5]+1)+afix->V[10]+afix->V[9];
    soff=round(dd);soff+=dfc->index[i+1]; /* ensuring goes in correct state
    vector location */
    k=0L; /* counter for location in x vector */
    for (j=0;j<round(afix->V[7]);j++)
    { if (t.V[j]==x[i].V[k])
      { dd=a.V[k];k++;}
      else
      { if (k) k--;
	dd=t.V[j]-x[i].V[k];
	dd=a.V[k]+b[i].V[k]*dd+c.V[k]*dd*dd+d.V[k]*dd*dd*dd;
	if (k) k++;
      }
      afix->V[soff+j*dfc->no_s]=dd;
    }

    /* history variable splines initialised here */

    hi=-1;
    for (j=0;j<dfc->no_hv;j++)
    if (hindex[j]==i) hi=j;   /* is there a his variable for state i*/
    if (hi>=0)
    dd=pvspliner(hi,dfc->no_uf,x[i],a,t.V[0],1,1); /* if so: set it up */
    freemat(a);freemat(c);freemat(d);
  }
  /* gradients must now be packed into revised Y vector, using b[i],x[i]
     and t and starting at the first sample time (afiv->V[13]) */
  k=0L;
  kk=(long *)calloc((size_t)round(afix->V[5]),sizeof(long));
  for (j=0;j<round(afix->V[5]);j++) /* move counters to first x>=t.V[afix->V[13]] */
  { kk[j]=0L;while (x[j].V[kk[j]]<t.V[round(afix->V[13])]) kk[j]++;
    dfc->rawcount[j]=0;
  }
  /* following loop fills Y and dfc structures that get transfered to o/p
     windows as raw data ....... */
  /* initialize the state data -> global vector map */
  Z=(long **)calloc((size_t)round(afix->V[5]),sizeof(long *));
  for (j=0;j<round(afix->V[5]);j++)
  { Z[j]=(long *) calloc((size_t)b[j].r,sizeof(long));
    for (i=0;i<b[j].r;i++) Z[j][i]=-1L;
  }
  for (i=round(afix->V[13]);i<t.r;i++) /* work through samples */
  for (j=0;j<round(afix->V[5]);j++) /* .... and states */
  { if (x[j].V[kk[j]]==t.V[i])
    { Y->V[k]=b[j].V[kk[j]];
      Z[j][kk[j]]=k; /* Z[j][i] is location in Y of ith element of b[j] */
      W->V[k]=w[j].V[kk[j]];
      dfc->raw[j][dfc->rawcount[j]].y=(float)Y->V[k];
      dfc->raw[j][dfc->rawcount[j]].x=(float)x[j].V[kk[j]];
      dfc->rawcount[j]++;
      k++;
      kk[j]++; /* move to next gradient for variable j */
    }
  }
  (*Co)=initmat(k,k);
  free(kk);  Y->r=k;W->r=k;
  /* now fill global Covariance matrix */
  for (j=0;j<round(afix->V[5]);j++) /*loop through states */
  { for (k=0;k<b[j].r;k++)
    if (Z[j][k]>=0L)
    for (i=0;i<b[j].r;i++)
    { if (Z[j][i]>=0L) Co->M[Z[j][i]][Z[j][k]]=Cog[j].M[i][k];
    } /* all done now wash your memory */
  }
}

ddefit_control initddefit_control(int maxco,int maxfunc,int maxvar)

{ ddefit_control dfc;
  int i;
  dfc.uco=(double *)calloc((size_t)maxco,sizeof(double));
  dfc.c=(double *)calloc((size_t)maxco,sizeof(double));
  dfc.uft0=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.uft1=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.uflb=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.ufub=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.ufsp0=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.ufdsp=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.ufspmax=(double *)calloc((size_t)maxfunc,sizeof(double));
  for (i=0;i<maxfunc;i++) dfc.ufspmax[i] = -1.0;
  dfc.ufdf=(int *)calloc((size_t)maxfunc,sizeof(int));
  dfc.ufnsp=(int *)calloc((size_t)maxfunc,sizeof(int));
  dfc.index=(int *)calloc((size_t)maxco,sizeof(int));
  dfc.uftype=(int *)calloc((size_t)maxfunc,sizeof(int));
  dfc.revmax=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.uctype=(int *)calloc((size_t)maxfunc,sizeof(int));
  dfc.uclb=(double *)calloc((size_t)maxco,sizeof(double));
  dfc.ucub=(double *)calloc((size_t)maxco,sizeof(double));
  dfc.ufsp=(double *)calloc((size_t)maxfunc,sizeof(double));
  dfc.wstats=(double *)calloc((size_t)MAXSTAT+1,sizeof(double));
  dfc.lines=(int *)calloc((size_t)maxvar,sizeof(int));
  dfc.windex=(wctype *)calloc((size_t)maxvar,sizeof(wctype));
  dfc.label=(char **)calloc((size_t)maxvar,sizeof(char *));
  dfc.wname=(char **)calloc((size_t)maxvar,sizeof(char *));
  dfc.cname=(char **)calloc((size_t)maxco,sizeof(char *));
  dfc.range=(y0y1type *)calloc((size_t)maxvar,sizeof(y0y1type));
  dfc.gradfit=0;
  dfc.discrete=0;
  dfc.gcvmethod=0;
  dfc.fitmethod=2;
  dfc.errors=0;
  dfc.ppnoise=0; // user can set to unity to generate process errors estimates
  for (i=0;i<maxvar;i++)
  { dfc.wname[i]=dfc.label[i]=" ";
    dfc.windex[i].win=-1;dfc.windex[i].cur=-1;
  }
  for (i=0;i<maxco;i++)
  { dfc.cname[i]=" ";}
  for (i=0;i<maxfunc;i++) dfc.ufsp[i]=-1.0;
  dfc.bsr_reps=0; 
  return(dfc);
}

void freeddefit_control(dfc) ddefit_control dfc;

{ 
  free(dfc.uco);
  free(dfc.c);
  free(dfc.uft0);
  free(dfc.uft1);
  free(dfc.uflb);
  free(dfc.ufub);
  free(dfc.uclb);
  free(dfc.ucub);
  free(dfc.ufsp0);
  free(dfc.ufdsp);
  free(dfc.ufspmax);
  free(dfc.ufnsp);
  free(dfc.index);
  free(dfc.ufdf);
  free(dfc.ufsp);
  free(dfc.uftype);
  free(dfc.revmax);
  free(dfc.uctype);
  free(dfc.wstats);
  free(dfc.lines);
  free(dfc.label);
  free(dfc.wname);
  free(dfc.cname);
  free(dfc.range);
  free(dfc.windex); // added 21/4/98
}

void outputfit(fit_control_type *fc, char *name)


/* writes an output text file called whatever is in name:
   time ,    fit , data   ,   fit , data ......
   missing data is coded "*"
*/

{ int i,j,k,z;
  double t1;
  FILE *ff;
  ff=fopen(name,"wt");
  k=0L; /* f counter */
  for (i=0;i<fc->n_st;i++) /* working through sample times */
  { t1=fc->ts[i];
    fprintf(ff,"%8.4g",t1);
    for (j=0;j<fc->n_fit;j++) /* work through all fit variables */
    { z=i*fc->n_fit+j;
      if (fc->yts[z]>-0.5) // indicates missing or not
      { fprintf(ff,"      %8.4g    %8.4g",fc->f[k],fc->y[k]);k++;}
      else fprintf(ff,"             *         *");
    }
    fprintf(ff,"\n");
  }
  fclose(ff);
}

void output_r2(fit_control_type *fc,char *name, double dferr)

/* calculates r-squared values for fitted series and outputs the results to file 
   called name.
*/
  
{ int i,j,k,z,ki,fi;
  double sf,fm,sy,ym,sfy,r,sf1,r12,r22;
  FILE *ff;
  matrix f,y;
  ff=fopen(name,"wt");
  f.V=fc->f;y.V=fc->y; // dummy for cheap re-use of old code!! 
  fprintf(ff,"State  Correlation  Explained Variance  Adjusted Explained Variance");
  fprintf(ff,"\n-----  -----------  ------------------  ---------------------------\n");
  fm=0.0;ym=0.0;k=0;
  for (i=0;i<fc->n_st;i++)  // work through sample times
  for (j=0;j<fc->n_fit;j++)
  { z=i*fc->n_fit+j;
    if (fc->yts[z]>-0.5) { fm+=f.V[k];ym+=y.V[k];k++;}
  }
  fm/=k;ym/=k;
  sf1=sf=sy=sfy=0.0;k=0;
  for (i=0;i<fc->n_st;i++)
  for (j=0;j<fc->n_fit;j++)
  { z=i*fc->n_fit+j;
    if (fc->yts[z]>-0.5)
    { sfy+=(f.V[k]-fm)*(y.V[k]-ym);
      sf+=(f.V[k]-fm)*(f.V[k]-fm);
      sy+=(y.V[k]-ym)*(y.V[k]-ym);
      sf1+=(f.V[k]-y.V[k])*(f.V[k]-y.V[k]);
      k++;
    }
  }
  r12=1.0-sf1/sy; // variance explained
  sy/=k-1; sf/=k-1;sfy/=k-1; sf1/=dferr;r22=1.0-sf1/sy; // corrected variance explained
  r=sfy/sqrt(sy*sf); // correlation coefficient
  fprintf(ff,"  All    %9.3g           %9.3g                    %9.3g\n",r*r,r12,r22);

  for (fi=0;fi<fc->n_fit;fi++) // work through fit variables
  { fm=0.0;ym=0.0;ki=0;k=0;
    for (i=0;i<fc->n_st;i++)  // work through sample times
    for (j=0;j<fc->n_fit;j++)
    { z=i*fc->n_fit+j;
      //k=i*(int) round(afix.V[5])+j;
      if (fc->yts[z]>-0.5)
      { if (j==fi)
        { fm+=f.V[k];ym+=y.V[k];ki++;}
        k++;
      }
    }
    fm/=ki;ym/=ki;
    sf1=sf=sy=sfy=0.0;k=0;
    for (i=0;i<fc->n_st;i++)
    for (j=0;j<fc->n_fit;j++)
    { z=i*fc->n_fit+j;
      //k=i*(int) round(afix.V[5])+j;
      if (fc->yts[z]>-0.5)
      { if (j==fi)
        { sfy+=(f.V[k]-fm)*(y.V[k]-ym);
          sf+=(f.V[k]-fm)*(f.V[k]-fm);
          sy+=(y.V[k]-ym)*(y.V[k]-ym);
          sf1+=(f.V[k]-y.V[k])*(f.V[k]-y.V[k]);
        }
        k++;
      }
    }
    r12=1.0-sf1/sy; // variance explained
    sy/=ki-1;
    sf/=ki-1;sfy/=ki-1;
    sf1/=dferr;
    r22=1.0-sf1/sy; // corrected variance explained
    r=sfy/sqrt(sy*sf); // correlation coefficient
    fprintf(ff,"   %2d    %9.3g           %9.3g                   *\n",fi,r*r,r12);
  }
  fclose(ff);
}









double slowgcv(int setup,matrix lam,int method, int error,
               int F(matrix,matrix,matrix,matrix,matrix,void*,int,double),
               matrix J,matrix *Z,matrix f, matrix a,fit_control_type *fc,matrix
               w,matrix *S,matrix y,matrix A, matrix Af,matrix b, double ro,
               double *maxviol,int bsreps,double *trinf,long *off,matrix *Coa)

/* designed for bruteforce GCV minimisation. Call first with all arguments and
   setup set to one. Subsequent calls should have setup set to zero, in which
   case all arguments except lam are ignored, and only lam.V is used.
   *Coa should be initialised on input - it doesn't matter what dimensions.
   */

{ static int method1,error1,bsreps1;
  static long *off1;
  static double *trinf1,*maxviol1,ro1;
  static int (*F1)(matrix,matrix,matrix,matrix,matrix,void*,int,double);
  static matrix J1,*Z1,f1,a1,w1,*S1,y1,A1,Af1,b1,lam1,*Coa1;
  static fit_control_type *fc1;
  matrix St;
  double gcv=0.0;
  int i,j,k;
  if (setup) // store arguments for later use
  { method1=method;error1=error;bsreps1=bsreps;
    trinf1=trinf;maxviol1=maxviol;ro1=ro;F1=F;off1=off;Coa1=Coa;
    J1=J;Z1=Z;f1=f;a1=a;fc1=fc;w1=w;S1=S;y1=y;A1=A;Af1=Af;b1=b;lam1=lam;
  } else     // call fitter() and calculate GCV, using current lam
  { St=initmat(a1.r,a1.r);
    { for (k=0;k<lam1.r;k++) for (i=0;i<S1[k].r;i++) for (j=0;j<S1[k].c;j++)
      St.M[i+off1[k]][j+off1[k]]+=S1[k].M[i][j]*lam.V[k];
    }
    matrixintegritycheck();
   
    fitter(method1,error1,F1,J1,Z1,f1,a1,fc1,w1,St,y1,A1,Af1,b1,ro1,maxviol1,bsreps1,trinf1);
    for (i=0;i<y1.r;i++) gcv+=w1.V[i]*(y1.V[i]-f1.V[i])*(y1.V[i]-f1.V[i]);
    gcv/=(y1.r- *trinf1)*(y1.r- *trinf1);
    freemat(St);
    return(gcv);
  }
}

double slowgcvcall(double *p)

/* this routine is used to call slowgcv from a minimising algorithm that expects
   to call double f(double*).... slowgcv must be called with setup=1, before
   using this */

{ matrix lam;
  double a;long i;
  fit_control_type *dull;
  lam.V=p;
  a=slowgcv(0,lam,0,0,NULL,lam,&lam,lam,lam,dull,lam,&lam,lam,lam,lam,
          lam,0,&a,0,&a,&i,&lam);
  return(a);
}

void fitter(method,error,F,J,Z,f,a,fc,w,S,y,A,Af,b,ro,maxviol,bsreps,trinf)
int (*F)(matrix,matrix,matrix,matrix,matrix,void*,int,double);

matrix J,*Z,f,a,w,S,y,A,Af,b;double ro,*maxviol,*trinf;fit_control_type *fc;
int bsreps,method,error;/* This routine calls one of 3 non-linear fit routines according to method:
   0 - NonLinLS a hybrid Guass-Newton steepest descent method
   1 - CQN a quasi-Newton method
   2 - IRLS an iterative least squares method.
   There are 3 error structure options:
   0 - normal errors. fitting routine called once only.
   1 - variance proportional to expectation
   2 - constant CV model
   Options other than 0 require iterative re-weighting, with weights
   set to reciprocal of estimated variance.

   bsreps should contain the number of bootstrap restarts required. b.s.
   restarts perturb the fitting objective, by using non-parametric bootstrap
   resamples in place of the original data. The perturbed objective usually
   has different small scale local minima than the original, whilst maintaining
   roughly the same large scale topography (statistically the differences are
   only significant at level p with probability p!). set bsreps to zero for no
   restarts.

   * returns *trinf == -1.0 for user stop
*/

{ int (*fr)(int (*F)(matrix,matrix,matrix,matrix,matrix,void*,int,double),
     matrix INF,matrix* ,matrix f,matrix a, void *prob_dat,matrix W,matrix S,matrix y,
     matrix A,matrix Af,matrix b,double ro,double *maxviol,int reset);
  int iter=0,i,j,k,notcon,reset=1,r,stop;
  double maxw,mean,fmin,ftry;
  matrix nw,olda,fold,an,aSa,wr;
  if (method==2) fr=IRLS;
  else if (method==1) fr=CQN;
  else fr=NonLinLS;
  if (error) /* then iterative re-weighting is required */
  { notcon=1;
    fold=initmat(f.r,1L);
    olda=initmat(a.r,1L);
    for (i=0;i<a.r;i++) olda.V[i]=a.V[i];
    nw=initmat(w.r,1L);
    while (notcon)
    { if (iter)
      for (i=0;i<f.r;i++) fold.V[i]=f.V[i];
       /* move parameters away from last optimum */
      for (i=0;i<a.r;i++) a.V[i]=a.V[i]*0.98+0.02*olda.V[i];
      stop=fr(F,J,Z,f,a,(void *) fc,w,S,y,A,Af,b,ro,maxviol,reset);
      if (stop) // user stop
      { freemat(*Z);
        break;
      }
      /* set new weights */
      if (!iter) for (i=0;i<f.r;i++) fold.V[i]=f.V[i];
      iter++;
      maxw=0.0;
      for (i=0;i<f.r;i++)
      { if (f.V[i]>0.0)
	     { if (error==1) nw.V[i]=2.0/(f.V[i]+fold.V[i]);
	       else nw.V[i]=1.0/sqrt(0.5*(f.V[i]+fold.V[i]));
	     }
	     if (maxw<nw.V[i]) maxw=nw.V[i];
      }
      for (i=0;i<f.r;i++) if (f.V[i]<=0.0) nw.V[i]=maxw;
      /* test for weight convergence */
      mean=0.0;
      for (i=0;i<f.r;i++) mean+=fabs(w.V[i]-nw.V[i])/(w.V[i]+nw.V[i]);
      mean/=f.r;
      if (mean<0.03) notcon=0;
      else {freemat(*Z);}
      for (i=0;i<f.r;i++) w.V[i]=nw.V[i];
    }
    freemat(nw);freemat(olda);freemat(fold);
  } else
  { /* The following implements bootstrap restarts. The idea is that one
       perturbs the objective function by a statistically insignificant amount,
       by non-parametric bootstrap. This should allow descent through small
       scale obstacles that can not be significant at the scale that is
       significant for fitting..... */
    if (bsreps)
    { if (S.r) aSa=initmat(1L,1L);wr=initmat(w.r,1L);an=initmat(a.r,1L);
      stop=fr(F,J,Z,f,a,(void *) fc,w,S,y,A,Af,b,ro,maxviol,reset);
      freemat(*Z);
      fmin=0.0;for (i=0;i<w.r;i++) fmin+=(f.V[i]-y.V[i])*w.V[i]*(f.V[i]-y.V[i]);
      if (S.r) {multi(3,aSa,a,S,a,1,0,0);fmin+=ro*y.r*aSa.V[0];}
      if (stop==0) // then user has not tried to stop fit
      for (r=0;r<bsreps;r++)
      { for (k=0;k<wr.r;k++) wr.V[k]=0.0; // clearing resampled weight vector
        for (k=0;k<w.r;k++)
        { j=(int)floor(wr.r*ran1()); if (j>wr.r-1) j=wr.r-1;
          wr.V[j]+=w.V[j];
        }
        for (i=0;i<a.r;i++) an.V[i]=a.V[i];
        stop=fr(F,J,Z,f,an,(void *) fc,wr,S,y,A,Af,b,ro,maxviol,reset); // b.s. trial        
        freemat(*Z);
        if (stop) break;
        stop=fr(F,J,Z,f,an,(void *) fc,w,S,y,A,Af,b,ro,maxviol,reset);  // polish it !
        freemat(*Z);
        if (stop) break;
        ftry=0.0;for (i=0;i<w.r;i++) ftry+=(f.V[i]-y.V[i])*w.V[i]*(f.V[i]-y.V[i]);
        if (S.r) {multi(3,aSa,an,S,an,1,0,0);ftry+=ro*y.r*aSa.V[0];}
        if (ftry<fmin)
        { fmin=ftry;
          for (i=0;i<a.r;i++) a.V[i]=an.V[i];
        }
      }
      if (S.r) freemat(aSa);freemat(wr);freemat(an);
    }
    /*.... end of bootstrap restarts, do a final restart with original weights */
  }
  for (i=0;i<1;i++) 
  { stop=fr(F,J,Z,f,a,(void *) fc,w,S,y,A,Af,b,ro,maxviol,reset); // final fit, using final weights
    if (stop) break;
  }
  if (stop) *trinf=-1.0;
}


double SimAn(matrix a,fit_control_type *fc,matrix amin,matrix amax,matrix y,matrix w,
           matrix f,int gradfit,int N,double rho,double T0)

/* Simulated annealing routine following Brooks and Morgan 1995. */

{ int s=1,i,j,k;
  matrix J,D,da,at;
  double T,F0,F1,z,aold,Fmin;
  T=T0;
  /* get initial function ..... */
  if (gradfit) G(f,J,D,a,da,fc,2,1.0); else F(f,J,D,a,da,fc,2,1.0);
  F0=0.0;for (i=0;i<f.r;i++) { z=f.V[i]-y.V[i];F0+=w.V[i]*z*z;} F0/=2;
  Fmin=F0;
  at=initmat(a.r,1L);mcopy(&a,&at);
  while (s>0)
  { s=0;
    for (i=0;i<N;i++)
    { j=(int)floor(ran1()*a.r);if (j==a.r) j--;        // pick a random element
      aold=at.V[j];
      at.V[j]=amin.V[j]+(amax.V[j]-amin.V[j])*ran1();  // pick a random value
      if (gradfit) G(f,J,D,at,da,fc,0,1.0); else F(f,J,D,at,da,fc,0,1.0);
      F1=0.0;for (k=0;k<f.r;k++) { z=f.V[k]-y.V[k];F1+=w.V[k]*z*z;} F1/=2;
      if (F1<Fmin) { mcopy(&at,&a);Fmin=F1;}
      if ((F1<F0)||(ran1()<=exp(-(F1-F0)/T))) { F0=F1;s++;} // accept
      else at.V[j]=aold;                                    // reject
    }
    T*=rho;
    qpoutdata.obj=F0;
    if (gradfit) G(f,J,D,at,da,fc,2,1.0); else F(f,J,D,at,da,fc,2,1.0);
  }
  if (Fmin!=F0)
  { qpoutdata.obj=Fmin;
    if (gradfit) G(f,J,D,a,da,fc,2,1.0); else F(f,J,D,a,da,fc,2,1.0);
  }
  return(Fmin);
}


int BootStrap(bsctype bsc,int (*F)(matrix,matrix,matrix,matrix,matrix,void*,int,double),
               matrix J,
               matrix f,matrix a,fit_control_type *fc,matrix W,matrix S,matrix y,matrix A,
               matrix Af,matrix b,double ro,double *maxviol)

/* Bootstrapping routine to get intervals or replicate parameter sets. f must
   contain E(y). y contains corresponding observations. User fills out the
   bootstrapping control structure bsc:

   bsc.method     - the fitting method
   bsc.error      - the error model
   bsc.parametric - set to 1 for parametric bootstrap. 0 for non-parametric. 
                    2 for semi-parametric - E(y) para, resids non-para.
   bsc.lumped     - set to 1 to treat all stages the same. 0 for seperate stages.
   bsc.iocontrol  - 0 for no data io. 1 to output rep data to file. 2 to read
                    rep data from a file.
   bsc.fdname     - name for data rep file name. output or input depending on iocontrol.
   bsc.fpname     - parameter rep file name. Always output.
   bsc.SSname     - name of file to receive weighted error SS for each rep.
   bsc.reps       - number of bootstrap replicates to perform
   bsc.carry_p    - set to zero to use the same starting paraemter values for each rep.
   bsc.n_start_file - number of start parameter files to use as starting values
   bsc.start_file  - array of start file names. Any file name starting with a "*" is taken 
                     to refer to a text file which contains one line of start parameters
                     for each replicate. There must only be one such file.
                     Otherwise the same start parameters are used for each rep.
   
   There is a facility to introduce data on extra state variables into the
   stored bootstrap replicates. It is not designed to operate with extra
   statistics (a technical rather than a conceptual limitation). This facility
   allows model comparison when the more complicated model also has more data
   to fit. The control structure items are:

   bsc.nextra     - the number of supplementary states.
   bsc.extra[j][i]- the ith observation for the jth extra observed state.
                    extra[j] must be of length afix.V[7] (the number of sample
                    times). For non-parametric bootstrapping this contains the
                    weight for the ith observation of the jth stage.
   bsc.exloc[j]   - which fitted variable the jth extra variable will
                    belong to in the output. So if bsc.exloc[0]=2 and
                    bsc.exloc[1]=4, then the original fitted variable 2 will
                    become fit variable 3, and the original 3 will become
                    fit variable 5....
   bsc.missing    - the double that will code for missing data in bsc.extra[][]
   bsc.ignore[]   - bsc.ignore[i]==1 if fit variable i is to be ignored in ss
                    calculation, 0 if it is to be included. This is to allow
                    inter model comprability when models fit different numbers
                    of series.
   ... implemented 15/1/99

   returns 0 if everything went ok and 1 for a user stop.
*/



{ double *e,**es,*mean,*stdev,*er,**ers,se,se2,ss,trinf,ssmin;
  int zo,ko,ke,z,k,i,j,*label,stages,samples,*nstage,*cstage,*index,bscreps;
  matrix f0,yr,ar,wr,fmin,amin,Z;
  FILE *fp,*fd,*fss,*fstart;
  char estr[100];
  f0=initmat(f.r,1L);yr=initmat(y.r,1L);ar=initmat(a.r,1L);wr=initmat(W.r,1L);
  fmin=initmat(f.r,1L);amin=initmat(a.r,1L);
  if (bsc.parametric&&bsc.error!=0)
  { ErrorMessage("Parametric bootstrapping only implemented for normal errors. Setting to normal",0);
    bsc.error=0;
  }
  fp=fopen(bsc.fpname,"wt");fclose(fp);
  ss=0.0;for (i=0;i<yr.r;i++) ss+=W.V[i]*(y.V[i]-f.V[i])*(y.V[i]-f.V[i]);
  fss=fopen(bsc.SSname,"wt");fprintf(fss,"%g (actual) \n",ss);fclose(fss);
  
  // process start file names.....
  fstart=(FILE *) NULL;
  for (i=0;i<bsc.n_start_files;i++)
  if (bsc.start_file[i][0]=='*') // then it's a multiple start value file
  { fstart=fopen(bsc.start_file[i]+1,"rt");
    if(!fstart)
    ErrorMessage("Multiple start parameter file not found",1);
    break;
  }

  // ioff=(int)round(afix.V[7])+20; redundant
  k=0;
  label=(int *)calloc((size_t) y.r, sizeof(int)); /* each element's stage */
  stages=fc->n_fit; /* that is stages to be fitted */
  samples=fc->n_st;
  nstage=(int *)calloc((size_t) stages,sizeof(int)); /* data per stage */
  /* the index in afix refers to state variables not fit stages - need to
       convert this index to a fitted stage index */
  index=(int *)calloc((size_t)fc->n_s,sizeof(int));
  for (i=0;i<samples;i++)
  for (j=0;j<stages;j++)
  { z=fc->yts[i*stages+j]; //   round(afix.V[ioff+i*stages+j]);
    if (z> -1) index[z]=j;
  } /* .... so index[z] is the data stage to which the zth state variable applies */
  for (i=0;i<samples;i++) for (j=0;j<stages;j++) /* filling label[] */
  { z=index[fc->yts[i*stages+j]];    // (int)round(afix.V[ioff+i*stages+j])];
    if (z > -1) /* then there is a datum for this i,j */
    { label[k]=z;k++;nstage[z]++;}
  }  /* label[k] is the data stage corresponding to the kth datum in y */
  if (bsc.iocontrol!=2) /* then routine is generating resamples */
  { e=(double *)calloc((size_t) y.r, sizeof(double)); /* residual array */
    er=(double *)calloc((size_t) y.r, sizeof(double)); /* bootstrap resid. array */
    cstage=(int *)calloc((size_t) stages,sizeof(int)); /* counters for stages */
    es=(double **)calloc((size_t) stages, sizeof(double *)); /* stage resid. array */
    ers=(double **)calloc((size_t) stages, sizeof(double *)); /* b.s. stage resids*/
    ers[0]=er;es[0]=e;
    for (i=0;i<stages-1;i++) { es[i+1]=es[i]+nstage[i];ers[i+1]=ers[i]+nstage[i];}
    /* now fill up the residual array(s) */
    for (i=0;i<y.r;i++) es[label[i]][cstage[label[i]]++]=y.V[i]-f.V[i];
    if (bsc.parametric) /* then parametric parameters must be generated */
    { mean=(double *)calloc((size_t)stages,sizeof(double));
      stdev=(double *)calloc((size_t)stages,sizeof(double));
      /* first get mean and variance */
      for (i=0;i<stages;i++)
      if (nstage[i]<2)
      { bsc.lumped=1;
        ErrorMessage("Insufficient data for unlumped bootstrap. Lumping...",0);
      }
      if (bsc.lumped) /* only one set of parameters */
      { se2=se=0.0;
        for (i=0;i<y.r;i++)
        { se+=e[i];se2+=e[i]*e[i];}
        se/=y.r;se2 = (se2 - y.r*se*se) / (y.r-1);
        se2=sqrt(se2);
        for (i=0;i<stages;i++) { mean[i]=se;stdev[i]=se2;}
      } else
      for (i=0;i<stages;i++)
      { se2=se=0.0;
        for (j=0;j<nstage[i];j++)
        { se+=es[i][j];se2+=es[i][j]*es[i][j];}
        se/=nstage[i];se2=(se2-nstage[i]*se*se)/(nstage[i]-1);se2=sqrt(se2);
        mean[i]=se;stdev[i]=se2;
      }
    }
    if (bsc.iocontrol==1) { fd=fopen(bsc.fdname,"wt");fclose(fd);}
  } else
  { fd=fopen(bsc.fdname,"rt");
    if (!fd)
    { sprintf(estr,"Input bootstrapping file %s could not be found.",bsc.fdname);
      ErrorMessage(estr,1);
    }
  }
  /* everything is now set up to actually bootstrap */
  for (z=0;z<bsc.reps;z++)
  { if (bsc.iocontrol==2) /* reps to be obtained from a file */
    { for (i=0;i<y.r;i++)
      { if (bsc.parametric)
        { fscanf(fd,"%lf",yr.V+i);wr.V[i]=W.V[i];}
        else
        { fscanf(fd,"%lf",wr.V+i);yr.V[i]=y.V[i];}
        if (feof(fd)&&i!=y.r-1&&z!=bsc.reps-1)
        { sprintf(estr,"Bootstrap data file %s does not contain enough data.",bsc.fdname);
          ErrorMessage(estr,1);
        }
      }
    } else /* reps to be generated */
    { if (bsc.parametric)
      { if (bsc.parametric==2) /* semi-parametric: E(y)- para; resids non-para */
        { if (bsc.lumped)
          for (i=0;i<y.r;i++)
          { k=(int)floor(ran1()*y.r); if (k>nstage[i]) k--;er[i]=e[k];}
          else
          for (i=0;i<stages;i++) for (j=0;j<nstage[i];j++)
          { k=(int)floor(ran1()*nstage[i]); if (k>nstage[i]) k--;
            ers[i][j]=es[i][k];
          }
        } else
        for (i=0;i<stages;i++) for (j=0;j<nstage[i];j++)
        ers[i][j]=ndev()*stdev[i]+mean[i];
        for (i=0;i<stages;i++) cstage[i]=0;
        for (i=0;i<y.r;i++) yr.V[i]=f.V[i]+ers[label[i]][cstage[label[i]]++];
        for (i=0;i<W.r;i++) wr.V[i]=W.V[i];
      } else /* non-parametric bootstrap */
      {  for (k=0;k<wr.r;k++) wr.V[k]=0.0; // clearing resampled weight vector
         for (k=0;k<W.r;k++)
         { j=(int)floor(wr.r*ran1()); if (j>wr.r-1) j=wr.r-1;
           wr.V[j]+=W.V[j];
         }
         for (i=0;i<y.r;i++) yr.V[i]=y.V[i];
      }
    }
    if (bsc.n_start_files) // then each rep to be repeated from several start values
    { for (k=0;k<bsc.n_start_files;k++)
      { // read start values into ar
        if (bsc.start_file[k][0]=='*') // then read from fstart
        { for (i=0;i<ar.r;i++) fscanf(fstart,"%lf",ar.V+i);}
        else
        { fp=fopen(bsc.start_file[k],"rt");
          if (!fp) ErrorMessage("Bootstrap start parameter file not found",1);
          for (i=0;i<ar.r;i++) fscanf(fp,"%lf",ar.V+i);
          fclose(fp);
        }
        // run fitter
        if (z>1) bscreps=bsc.restarts[2]; else bscreps=bsc.restarts[z];
        fitter(bsc.method,bsc.error,F,J,&Z,f0,ar,fc,wr,S,yr,A,Af,b,ro,maxviol,bscreps,&trinf);
        freemat(Z);
        if (trinf==-1.0) // then user has cancelled...... DO SOMETHING!!
        { break;
        }
        // obtain SS value
        ss=0.0;
        for (i=0;i<yr.r;i++)
        ss+= wr.V[i]*(yr.V[i]-f0.V[i])*(yr.V[i]-f0.V[i])*
         (1-bsc.ignore[index[label[i]]]);
        // store f0 and p if SS value is lowest yet 
        if (!k||ss<ssmin)
        { ssmin=ss;
          for (i=0;i<f0.r;i++) fmin.V[i]=f0.V[i];
          for (i=0;i<ar.r;i++) amin.V[i]=ar.V[i];
        } 
      }
      // check for user break
      if (trinf==-1.0) // then user has cancelled...... DO SOMETHING!!
      { break;
      }
      for (i=0;i<f0.r;i++) f0.V[i]=fmin.V[i];
      for (i=0;i<ar.r;i++) ar.V[i]=amin.V[i];      
    } else // just one start value per rep.
    { if (!bsc.carry_p||!z) for (i=0;i<a.r;i++) ar.V[i]=a.V[i];
      /* at this stage there is resampled vector in yr, call fitter with yr and f0 */
      if (z>1) bscreps=bsc.restarts[2]; else bscreps=bsc.restarts[z];
    
      fitter(bsc.method,bsc.error,F,J,&Z,f0,ar,fc,wr,S,yr,A,Af,b,ro,maxviol,bscreps,&trinf);
      freemat(Z);
      if (trinf==-1.0) // then user has cancelled...... DO SOMETHING!!
      { break;
      }
    } 
    /* now save the resampled data vector */
    if (bsc.iocontrol==1) /* output the current resample */
    { fd=fopen(bsc.fdname,"at");
      if (bsc.nextra)
      { zo=0;
        for (i=0;i<samples;i++)
        { ko=ke=0;
          for (k=0;k<bsc.nextra+stages;k++)
          { if (bsc.exloc[ke]!=k) /* write out original data */
            { if (fc->yts[i*stages+ko]> -1)  // afix.V[ioff+i*stages+ko]> -0.5)
              { fprintf(fd,"%g  ",(bsc.parametric==0)*wr.V[zo]+(bsc.parametric>0)*yr.V[zo]);
                zo++;
              }
              ko++;
            } else /* write out extra data */
            { if (bsc.extra[ke][i]!=bsc.missing)
              { fprintf(fd,"%g  ",bsc.extra[ke][i]);}
              ke++;if (ke==bsc.nextra) ke=0;
            }
          }
        }
      } else
      for (i=0;i<y.r;i++) /* outputs weights for non-parametric, data otherwise */
      fprintf(fd,"%g  ",(bsc.parametric==0)*wr.V[i]+(bsc.parametric>0)*yr.V[i]);
      fprintf(fd,"\n");fclose(fd);
    }
    /* now output current parameter vector and residual ss  */
    ss=0.0;
    for (i=0;i<yr.r;i++)
    ss+= wr.V[i]*(yr.V[i]-f0.V[i])*(yr.V[i]-f0.V[i])*
         (1-bsc.ignore[index[label[i]]]);
    fss=fopen(bsc.SSname,"at");fprintf(fss,"%g\n",ss);fclose(fss);
    fp=fopen(bsc.fpname,"at");
    for (i=0;i<a.r;i++) fprintf(fp,"%g  ",ar.V[i]);
    fprintf(fp,"\n");
    fclose(fp);

  }
  if (fstart) fclose(fstart);
  /* clear up */
  freemat(amin);freemat(fmin);
  if (bsc.iocontrol==2) fclose(fd);
  freemat(f0);freemat(yr);freemat(ar);
  if (bsc.parametric&&(bsc.iocontrol!=2)) {free(mean);free(stdev);}
  free(index);free(label);free(nstage);
  if (bsc.iocontrol!=2)
  { free(e);free(er);free(es);free(ers);free(cstage);
  }
  if (trinf==-1.0) return(1);
  else return(0);
}

typedef struct
{ int no,*knots,*type;
  double **x;
} splinestore;

splinestore StoreSplineKnots(ddefit_control *dfc)

/* stores knot positions of ufs (if there are any) */

{ splinestore sk;int i,j;
  double dx;
  { sk.no=(int)dfc->no_uf;
    if (sk.no)
    { sk.knots=(int *)calloc((size_t)sk.no,sizeof(int));
      sk.x=(double **)calloc((size_t)sk.no,sizeof(double *));
      sk.type=(int *)calloc((size_t)sk.no,sizeof(int));
    }
    for (i=0;i<sk.no;i++)
    { sk.knots[i]=dfc->ufdf[i];
      sk.x[i]=(double *)calloc((size_t)sk.knots[i],sizeof(double));
      dx=dfc->uft1[i]-dfc->uft0[i];dx/=sk.knots[i]-1;
      for (j=0;j<sk.knots[i];j++)
      sk.x[i][j]=dfc->uft0[i]+j*dx;
      sk.type[i]=dfc->uftype[i];
    }
  }
  return(sk);
}





/* the following are test function ....... */
void frignoise(double *e,double *y,double *s,double *c,double t0,double t1)

/* this is a test function designed to work only with single state data series
   it figures out the process noise that would have given exactly the observed
   data*/

{ FILE *f;
  int i=0;
  double t;
  initst(s,c,t0);
  map(s,c,t0,0);
  e[0]=y[0]-s[0];t=t0;
  f=fopen("dde.err","wt");
  fprintf(f,"%d  %g  %g\n",0,e[0],s[0]);
  while (t<t1)
  { t+=1.0;t=(double)round(t);
    s[0]=y[i];i++;
    map(s,c,t,0);
    e[i]=y[i]-s[0];
    fprintf(f,"%d  %g  %g\n",0,e[i],s[0]);
  }

  fclose(f);
}



void frig(matrix y,matrix a,matrix afix)

/* Routine for calling frignoise, once a fit has been obtained
*/


{ double *e,*s,*c,t0,b,t1,*cp;
  int i,j,ns,coff,ioff,z;
  long k;
  matrix p;
  /* new code */
  ns=round(afix.V[2]);
  s=(double *)calloc((size_t)ns,sizeof(double));
  cp=c=(double *)calloc((size_t)(round(afix.V[8])+round(afix.V[10])),sizeof(double));
  ioff=20+round(afix.V[7]);

  /* set up constants */
  c+=round(afix.V[8]);
  coff=20+round(afix.V[7])*(1+round(afix.V[5]));
  for (j=0;j<round(afix.V[10]);j++) c[j]=afix.V[coff+j];
  for (j=0;j<round(afix.V[8]);j++) c[-1-j]=a.V[j];
  /* Setting up splines */
  k=round(afix.V[8]);
  /* z is offset for uf df */
  z=(long)(round(afix.V[5]+1)*round(afix.V[7])+round(afix.V[/*8*/10]))+20L;
  for (j=0;j<round(afix.V[9]);j++)
  { b=afix.V[z+j];
    p=initmat((long)round(b),1L);
    for (i=0;i<round(b);i++)
    { p.V[i]=a.V[k];k++;}
    spl(j,p,p,afix.V[6],0.0,1,0,1);
    freemat(p);
  }
  t0=afix.V[20];t1=afix.V[20+(int)round(afix.V[7])-1];
  e=(double *)calloc((size_t)y.r,sizeof(double));
  frignoise(e,y.V,s,c,t0,t1);
  free(s);
  c-=round(afix.V[8]);free(c);
}

/* end */

void filecoeffs(char *name,matrix a,int out)

{ FILE *f;
  int i;
  double x;
  if (out)
  { f=fopen(name,"wt");
    for (i=0;i<a.r;i++)
    { fprintf(f,"%30.18g\n",a.V[i]);
    }
    fclose(f);
  } else
  { f=fopen(name,"rt");
    i=0;
    if (f)
    { while(!feof(f)) { fscanf(f,"%lf",&x);i++;} i--;
      if (i<a.r)
      ErrorMessage("Initial parameter file has too few rows, initializing with compiled defaults!",0);
      else if (i>a.r)
      ErrorMessage("Initial parameter file has too many rows, initializing with compiled defaults",0);
      else
      { fseek(f,0,0);
        for (i=0;i<a.r;i++)
        { fscanf(f,"%lf",&x);a.V[i]=x;
        }
      }
	  fclose(f);
    }
  }
}

void bscdefaults(int n_fit,bsctype *bsc)
/* does default initialisation of bsc control */

{ bsc->exloc=(int *)calloc((size_t)n_fit,sizeof(int));
  bsc->ignore= (int *)calloc((size_t)n_fit,sizeof(int));
  bsc->reps=0;
  bsc->parametric=0;
  bsc->lumped=1;
  bsc->iocontrol=0;
  bsc->fdname="dde.bsd";
  bsc->fpname="dde.bsp";
  bsc->SSname="dde.bss";
  bsc->restarts[0]=bsc->restarts[1]=bsc->restarts[2]=0;
  bsc->carry_p=0;
  bsc->nextra=0;
  bsc->missing=-1.2325335664454578e45;
}

void init_display_data_op_type(display_data_op_type *dd,int nt,int ns,int nuf,int *distos)

/* routine to initialise memory in structure of type used for exporting data to GUI thread */

{ int i;
  dd->n_disp=ns;
  dd->max_n=dd->n_t=nt;
  dd->ufxmax=(double *)calloc((size_t)nuf,sizeof(double));
  dd->ufxmin=(double *)calloc((size_t)nuf,sizeof(double));
  dd->distos=(int *) calloc((size_t)ns,sizeof(int));
  for (i=0;i<ns;i++) dd->distos[i]=distos[i];
  dd->t=(double *) calloc((size_t) nt,sizeof(double));
  dd->display=(double **) calloc((size_t) ns,sizeof(double *));
  for (i=0;i<ns;i++) dd->display[i]=(double *) calloc((size_t) nt,sizeof(double));
}

void free_display_data_op_type(display_data_op_type *dd)

/* routine to free memory in structure of type used for exporting data to GUI thread
   - will usually be called from GUI thread.
 */

{ int i;
  free(dd->distos);
  for (i=0;i<dd->n_disp;i++) free(dd->display[i]);
  free(dd->display);
  free(dd->t);
  free(dd->ufxmax);free(dd->ufxmin);
  free(dd);
}

void simulation(PVOID pvoid)

/* This thread is designed to execute the model once and post the results back to 
   the parent window that called it. It should be called with:
   _beginthread(simulation,0,&simcontrol);
   where simcontrol is of type simulation_control_type as defined in ddefit95.h
   Simulated data  is stored in a structure of type display_data_op_type defined
   in ddefit.h. This type is initialised here, but will be destroyed by the receiving 
   window after use. Same structure is used to post simulated, but noise free fit variable
   data to parent. March 2000.
*/

{ matrix min,max;
  simulation_control_type *sc;
  display_data_op_type *outp,*outf;
  double *s,*c,*cp,t0,dout,dt;
  int nt,i,j,k,reset;
  matrix x,y;
  sc=(simulation_control_type *) pvoid;
  outp=(display_data_op_type *)calloc(1,sizeof(display_data_op_type));
  outf=(display_data_op_type *)calloc(1,sizeof(display_data_op_type));
  s=(double *)calloc((size_t)sc->n_s,sizeof(double));
  cp=c=(double *)calloc((size_t)(sc->n_uc+sc->n_c),sizeof(double));
  
   /* set up constants */
  c+= sc->n_uc;
  for (j=0;j<sc->n_c;j++) c[j]=sc->c[j];
  for (j=0;j<sc->n_uc;j++) c[-1-j]=sc->p[j];

  k=sc->n_uc;
  for (j=0;j<sc->n_uf;j++)   // setting up splines for current run
  { y=initmat((long)sc->ufdf[j],1L);
    x=initmat((long)sc->ufdf[j],1L);
    for (i=0;i<sc->ufdf[j];i++)
    { y.V[i]=sc->p[k];x.V[i]=sc->ufx[j][i];k++;}
    spl(j,x,y,x.V[1],0.0,1,1,1);
    freemat(y);freemat(x);
  }
  //splp=spl; // assigning pointer to splines so that uf() becomes usable
  splp=audituf; // assigning pointer to audit uf so that uf() usable and function domains recorded
  min=initmat((long)sc->n_uf,1L);max=initmat((long)sc->n_uf,1L);
  audituf(sc->n_uf,min,max,0.0,0.0,1,1,0); // clearing the range (and indicating number of ufs)
   
  
  initst(s,c,sc->t0);
  reset=1;
  dout=sc->dt; /* setting o/p timestep for defualt output routine */
  k=0;
  nt=(int)floor((sc->t1-sc->t0)/dout)+1; // number of output times
  // set up output data structures
  init_display_data_op_type(outp,nt,sc->n_dis,sc->n_uf,sc->distos);
  init_display_data_op_type(outf,nt,sc->n_dis,sc->n_uf,sc->distos);
  t0=sc->t0;
  outf->t[0]=t0;
  for (j=0;j<sc->n_dis;j++) if (sc->fitdv[j]) outf->display[j][0]=s[sc->distos[j]];
  output(s,0.0,(void *)outp,1);   // initialise function output() for output from dde
  dt=sc->dt/100;
  for (i=0;i<nt-1;i++) /* working through sample times */
  { ddep(s,c,t0,t0+sc->dt,&dt,sc->tol,dout,sc->n_s,sc->n_sw,sc->n_hv,sc->hbsize,sc->n_lag,reset,0);
    reset=0;
    t0+=sc->dt;
    outf->t[i+1]=t0;
    for (j=0;j<sc->n_dis;j++) if (sc->fitdv[j]) 
    outf->display[j][i+1]=s[sc->distos[j]];
  }
  if (sc->noise>0.0) // then add noise to output data in outf
  { for (j=0;j<sc->n_dis;j++) if (sc->fitdv[j]) 
    for (i=0;i<outf->n_t;i++)
    outf->display[j][i]+=ndev()*sc->noise;
  } 
  audituf(0,min,max,0.0,0.0,1,2,0); // finding out what the x ranges of uf's  were
  for (i=0;i<sc->n_uf;i++) 
  { outp->ufxmax[i]=max.V[i];
    outp->ufxmin[i]=min.V[i]; 
  }
  freemat(min);freemat(max);
  // the following might be better done by using SendMessage() ....
  // (SendMessage calls window function and only returns when done)....
  PostMessage(sc->parent,WM_COMMAND,(WPARAM)MAKELONG(3,0),(LPARAM)outp);
  PostMessage(sc->parent,WM_COMMAND,(WPARAM)MAKELONG(4,0),(LPARAM)outf);
  free(s);free(cp);
  
  _endthread();
}

matrix covariance(matrix *J,matrix *Z,matrix *y, matrix *w, matrix *f,matrix *S,
                  double *lam,int *off,int m,int *edf)

/* Routine for finding estimated covariance matrix from fitted model.

   J is jacobian J_{ij}=df_i/dp_j
   Z is null space (actually stored as Householder rotations)
   y is data
   w is weight vector
   f is fitted values
   S[] is array of smoothness matrices
   lam[] is array of smoothing parameters
   off[] - start point offsets. S[k].M[0][0] should be in element off[k],off[k] of
           the overall penalty matrix.
   m is number of penalties
   edf - the estimated degrees of freedom for the model will be returned in this.

   Let A=JZ[Z'(J'W'J + \sum \lam_i S_i) Z]^{-1} Z'J'W

   Then define s^2=(y-f)'W(y-f)/tr(I-A)

   The estimated covariance matrix of the parameters is given by:

   s^2 Z[Z'(J'W'J + \sum \lambda_i S_i) Z]^{-1} Z'

*/
{ matrix Vp,A,B;
  double x,tr,rss,s2;
  int i,j,k;
  Vp=initmat(J->c,J->c);
  for (i=0;i<J->c;i++) for (j=0;j<=i;j++)
  { x=0.0;for (k=0;k<J->r;k++) x+=J->M[k][j]*w->V[k]*J->M[k][i];
    Vp.M[i][j]=Vp.M[j][i]=x;
  } // Vp contains J'WJ
  for (k=0;k<m;k++)
  { for (i=0;i<S[k].r;i++)
    for (j=0;j<S[k].c;j++) 
    Vp.M[i+off[k]][j+off[k]]+=lam[k]*S[k].M[i][j];
  } // Vp contains J'WJ + \sum \lambda_i S_i
   
  HQmult(Vp,*Z,1,1);
  HQmult(Vp,*Z,0,0);
  Vp.r += -Z->r;Vp.c += -Z->r;
  // Vp contains  Z'(J'WJ + \sum \lambda_i S_i)Z
  A=initmat(Vp.r,Z->c);A.c=Vp.r;
  if (chol(Vp,A,1,0)) // then Z'(J'WJ+ \sum \lam_i S_i) Z is +ve def - Choleski factor in A
  { A.c=Z->c;
    HQmult(A,*Z,0,1);  
    freemat(Vp);
    Vp=initmat(Z->c,Z->c);
    for (i=0;i<Vp.c;i++) for (j=0;j<=i;j++)
    { x=0.0;for (k=0;k<A.r;k++) x+=A.M[k][i]*A.M[k][j];
      Vp.M[i][j]=Vp.M[j][i]=x;
    } // Now Vp contains Z[Z'(J'WJ+ \sum \lam_i S_i) Z]^{-1} Z'
    if (m) // then there are penalties and edf is not so easy
    { B=initmat(J->r,A.r); 
      matmult(B,*J,A,0,1); // .. so B'B=JZ[Z'(J'WJ+ \sum \lam_i S_i) Z]^{-1} Z'J
      freemat(A);
      A=initmat(J->r,1L); // for storing the elements of the l.d. of BB'
      for (i=0;i<A.r;i++)
      { A.V[i]=0.0; for (j=0;j<B.c;j++) A.V[i]+=B.M[i][j]*B.M[i][j];}
      freemat(B);
      tr=0.0;
      for (i=0;i<A.r;i++) 
      tr+=A.V[i]*w->V[i];
    } else
    { tr=J->c-Z->r; // edf for model
    }
    rss=0.0;
    for (i=0;i<y->r;i++) rss+= (y->V[i]-f->V[i])*w->V[i]*(y->V[i]-f->V[i]);
    s2=rss/(y->r-tr);
    for (i=0;i<Vp.r;i++) for (j=0;j<Vp.c;j++) Vp.M[i][j]*=s2;  
  } else 
  { ErrorMessage("Singularity problem in covariance matrix estimation. Ignore variance estimates.",0);
    freemat(Vp);Vp=initmat(J->c,J->c);
    for (i=0;i<J->c;i++) Vp.M[i][i]=1.0;
  }
  freemat(A);
  return(Vp);
}

void fit_thread(PVOID pvoid)

/* This is the fitting thread function. Its purpose is to perform a single model fit 
   (including any bootstrap restarting). It is called with something like:
   _beginthread(fit_thread,0,&fc);
   where fitcontrol is a structure of fit_control_type. 
   Single trials of parameter vectors are also performed by this thread.
   The thread does not read in data, but assumes that it is already in fc->y and fc->w.
   The thread does the following:

   * assign splp 
   * set up splines
   * set up penalty matrices and constraint matrices
   * warn if parameters do not meet constraints
   * perform fit or trial 
   * clean up splines
   * clean up penalty and constraint matrices 
   
   Output is controlled from F() and happens whenever fc->trial is non-zero or after 
   the calculation of a new Jacobian.
*/
{ status_op_type *statop;
  fit_control_type *fc; 
  int i,j,k,*off;
  matrix x,*S,A,Af,b,p,y,w,f,J,Z,St,Vp,bt;
  double ro=0.0,maxviol,trinf,obj,*pp;
  fc=(fit_control_type *) pvoid; // cast pvoid to the correct type
  splp=spl;                      // assign splp to the correct function
  
  /* now set up Penalty S[i] and Constraint, A, Af and b, matrices along with the splines */
  if (fc->n_uf) 
  { S = (matrix *)calloc((size_t)fc->n_uf,sizeof(matrix));
    off = (int *)calloc((size_t)fc->n_uf,sizeof(int));
  }
  setupunknowns(fc,S,off,&A,&Af,&b,&p); // splines set up by this - p initialised from fc->p
  /* Now need to set up remaining matrices for fitting: w, y, p, f, J   */
  y.r=w.r=fc->n;y.c=w.c=1L;y.V=fc->y;w.V=fc->w; // data and weight vectors
 // p.r=fc->n_p;p.c=1L;p.V=fc->p;                 // parameter vector - incorrect - set up in setupunknowns()
  f=initmat(y.r,1L);                            // fitted value vector
  J=initmat(y.r,p.r);                           // Jacobian J_ij = df_i/dp_j
  // Now check whether smoothing parameters are to be estimated?
  for (i=0;i<fc->n_uf;i++) if (fc->ufsp[i]<0.0) fc->get_sp=1;
  if (fc->n_uf==0) fc->get_sp=0;  
  // Check that initial parameters satisfy the inequality constraints
  bt=initmat(b.r,1L);
  matmult(bt,A,p,0,0);
  for (i=0;i<b.r;i++) if (TOL*enorm(A)<b.V[i]-bt.V[i])
  { ErrorMessage("Trial parameters violate constraints",0);break;}
  freemat(bt);
  // Ready for trial/fitting .....
  if (fc->trial)  // then run model once before returning
  { k=0; // signals regular precision
    St.r=St.c=0L; // signalling to ignore penalties  
    obj=objective(F,f,J,J,p,(void *)fc,J,y,w,St,J,J,0.0,0.0,0L,2,&k);
    // ... repeated J's are for matrices not needed for a call that is not calculating J.
    // display should be automatic, but now post objective information to the o/p thread
    // NOTE CODE NEEDED HERE!!!!!! 
    if (k!=-2) // objective returned ok without user termination, so send updated objective
    { statop=(status_op_type *)calloc(1,sizeof(status_op_type));
      statop->F=obj;
      PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(2,0),(LPARAM)statop);
    }
  } else // we're fitting 
  if (fc->get_sp==0) // Don't estimate any smoothing parameters
  { // need to sum any penalty matrices to get one penalty St.....
    if (fc->n_uf) St=initmat(p.r,p.r);else St.r=St.c=0L;
    for (k=0;k<fc->n_uf;k++) for (i=0;i<S[k].r;i++) for (j=0;j<S[k].c;j++)
    St.M[i+off[k]][j+off[k]]+=S[k].M[i][j]*fc->ufsp[k];
    // Now fit the model......
    fitter(fc->fitmethod,0,F,J,&Z,f,p,fc,w,St,y,A,Af,b,ro,&maxviol,fc->bsr_reps,&trinf);
    
    if (fc->n_uf) freemat(St); // free total penalty matrix
  } else // estimating smoothing parameters
  { ro=1.0;
    if (fc->gcv_method==0)
    { // NOTE that fc->ufsp passed in unaltered
      optNLLS(F,J,&Z,f,p,(void *)fc,w,S,y,A,Af,b,off,&trinf,fc->ufsp,fc->ufspmax,
                  fc->n_uf,fc->fitmethod);
      // trinf == -1.0 now indicates user termination
    } else
    { // GCV by grid search
      /*Vp=initmat(1L,1L); // sacrificial
      slowgcv(1,lam,dfc.fitmethod,dfc.errors,F,&ZGZ,J,&Z,f,a,afix,W,S,Y,A,Af,b,ro,
              &maxviol,
              dfc.bsr_reps,
              &trinf,
              off,
              &Vp);
      if (dfc.gcvmethod==1)
      fastcrawl(slowgcvcall,lam.r,nk,p0,dp,lam.V,1,1); else
      slowcrawl(slowgcvcall,lam.r,nk,p0,dp,lam.V,1,1);
      free(nk);free(p0);free(dp);*/
    }
  }
  if (!fc->trial && trinf!=-1)
  { Vp=covariance(&J,&Z,&y,&w,&f,S,fc->ufsp,off,fc->n_uf,&(fc->edf)); // get the parameter covariance matrix
    freemat(Z); // initialized in fitter/optNlLS
    // Now unload Vp into fc->Vp
    for (i=0;i<Vp.r;i++) for (j=0;j<Vp.c;j++) fc->Vp[i][j]=Vp.M[i][j];
    freemat(Vp);
  }
  // NOTE: should unload fitted values into fc at this point, before cleaning up (unless trial)
  for (i=0;i<f.r;i++) fc->f[i]=f.V[i];
  if (!fc->trial) for (i=0;i<fc->n_p;i++) fc->p[i]=p.V[i];
  // perform bootstrapping if this has been requested...... 
  if (fc->bsc.reps)  BootStrap(fc->bsc,F,J,f,p,fc,w,St,y,A,Af,b,ro,&maxviol); 
  // need to force a final run with the best fit parameters and fc->trial==1 in order to 
  // force output....... CODE NEEDED ..... 
  
  /* clean up before returning */
  x.r=0;for (j=0;j<fc->n_uf;j++) spl(j,x,x,0.0,0.0,1,1,1); // clean up spline space
  // clean up matrices 
  for (i=0;i<fc->n_uf;i++) freemat(S[i]);if (fc->n_uf) { free(S);free(off);}
  freemat(f);freemat(J);
  // these 4 were initialised by setupunknowns().....
  freemat(A);
  freemat(Af);
  freemat(b);
  
  freemat(p);
  
  // now tell parent it's over as last action before ending thread
  if (fc->trial) PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(111,0),(LPARAM)0);
  else 
  { pp=(double *)calloc((size_t)fc->n_p,sizeof(double));
    for (i=0;i<fc->n_p;i++) pp[i]=fc->p[i];
    PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(1,0),(LPARAM)pp);
    PostMessage(fc->parent,WM_COMMAND,(WPARAM)MAKELONG(121,0),(LPARAM)0);
  }
  _endthread(); 
}