/home/peifer/diffit/simit.cc File Reference

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include "def.h"
#include "model.h"
#include "nr.h"

Go to the source code of this file.

Functions

GlobExpparseopts (int argc, char *argv[], Glob *globs, char *outstr)
void initialise (GlobExp ex[], Glob *globs, int simit)
void tabulateValues (Glob *globs, GlobExp *ex, double t0, double t1, double *t, long n, double *state, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp)
void outSimit (GlobExp ex[], Glob *globs, double *t, long n, double **y)
void freeMem (GlobExp *ex, Glob *globs, int simit)
 main (int argc, char *argv[])

Variables

ofstream * dbg
 debug stream (C++)


Function Documentation

void freeMem GlobExp ex,
Glob globs,
int  simit
 

Definition at line 11 of file freeMem.cc.

00012 {
00013   long nExp,i;
00014   long nrExp=globs->nrExp;
00015 
00016 
00017   delete globs->gnuFp;
00018   
00019   if(simit==FALSE)
00020     {
00021       for (nExp=1;nExp<=nrExp;++nExp) 
00022         { 
00023           long nPoints=ex[nExp].nPoints;
00024           long nP=globs->npar;
00025           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
00026           long nMeas=ex[nExp].nMeasure;
00027           long  me=ex[nExp].me, mg=ex[nExp].mg;
00028           free_dvector(ex[nExp].errP,1,nP);
00029           free_dvector(ex[nExp].errY0,1,nvar);
00030           free_dmatrix(ex[nExp].yTry,1,nPoints,1,nvar);
00031           free_dmatrix(ex[nExp].yTrySave,1,nPoints,1,nvar);
00032           free_dmatrix(ex[nExp].yComp,1,nPoints,1,nvar);
00033           free_dmatrix(ex[nExp].yPred,1,nMeas,1,nobs);
00034           free_dmatrix(ex[nExp].h,1,nPoints,1,nvar);   
00035           free_dmatrix(ex[nExp].residues,1,nMeas,1,nobs);
00036           free_d3tensor(ex[nExp].dyds,1,nPoints,1,nvar,1,nvar);
00037           free_d3tensor(ex[nExp].dydp,1,nPoints,1,nvar,1,nP);
00038           free_d3tensor(ex[nExp].dmds,1,nMeas,1,nobs,1,nvar);
00039           free_d3tensor(ex[nExp].dmdp,1,nMeas,1,nobs,1,nP);
00040           if (me>0) 
00041             {
00042               free_dvector(ex[nExp].r2,1,me);
00043               free_d3tensor(ex[nExp].dR2ds,1,me,1,nPoints,1,nvar);
00044               free_dmatrix(ex[nExp].dR2dp,1,me,1,nP);
00045             }
00046           if (mg>0) 
00047             {
00048               free_dvector(ex[nExp].r3,1,mg);
00049               free_d3tensor(ex[nExp].dR3ds,1,mg,1,nPoints,1,nvar);
00050               free_dmatrix(ex[nExp].dR3dp,1,mg,1,nP);
00051             }
00052           
00053           free_dvector(ex[nExp].ua,1,nMeas*nobs);
00054           free_dmatrix(ex[nExp].Ea,1,nMeas*nobs,1,nvar);
00055           free_dmatrix(ex[nExp].Pa,1,nMeas*nobs,1,nP);
00056           free_dvector(ex[nExp].ue,1,me);
00057           free_dmatrix(ex[nExp].Ee,1,me,1,nvar);
00058           free_dmatrix(ex[nExp].Pe,1,me,1,nP);
00059           free_dvector(ex[nExp].ug,1,mg);
00060           free_dmatrix(ex[nExp].Eg,1,mg,1,nvar);
00061           free_dmatrix(ex[nExp].Pg,1,mg,1,nP);
00062           free_dmatrix(ex[nExp].dS,1,nPoints,1,nvar);
00063           free_dvector(ex[nExp].dP,1,nP);
00064           
00065         }
00066     }
00067   for (nExp=1;nExp<=nrExp;++nExp) 
00068     {
00069       free(ex[nExp].splineNodes);
00070       free(ex[nExp].splineY);
00071       free(ex[nExp].splineGam);
00072       free(ex[nExp].nNodes);
00073       for(i=1;i<=NSPLINES;i++)
00074         {
00075           free(ex[nExp].splineNodes[i]);
00076           free(ex[nExp].splineY[i]);
00077           free(ex[nExp].splineGam[i]);
00078         }
00079     }
00080   if(simit==FALSE && globs->covar!=NULL)
00081     free_dmatrix(globs->covar,1,globs->fitdim,1,globs->fitdim);
00082 }

void initialise GlobExp  ex[],
Glob globs,
int  simit
 

Definition at line 98 of file initialise.cc.

00099 {
00100   long i,j,k,nExp;
00101   long nrExp=globs->nrExp;
00102   
00103   ifstream inSplines;
00104 
00105   if(simit==FALSE)
00106     {
00107       for (nExp=1;nExp<=nrExp;++nExp) 
00108         { 
00109           
00110           //non-local parameters are equal for each experiment
00111           for(i=1;i<=globs->npar;i++)
00112             {
00113               if(globs->doP[i]==TRUE)
00114                 ex[nExp].par[i]=ex[1].par[i];
00115             }
00116           
00117           // set ex->me and ex->mg given in model.cc
00118           setNrConstraints(&ex[nExp],globs->npar,ex[nExp].par);
00119           
00120           long nPoints=ex[nExp].nPoints;
00121           long nP=globs->npar;
00122           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
00123           long nMeas=ex[nExp].nMeasure;
00124           long  me=ex[nExp].me, mg=ex[nExp].mg;
00125           
00126 #ifdef PRINTINITVALUES
00127           *dbg << "initialize: exp. #" << nExp << ", nPoints=" << nPoints
00128                << ", nvar=" << nvar << ", nobs=" << nobs << "\nnMeas="
00129                << nMeas << ", nVal=" << nVal << ", me=" << me
00130                << ", mg=" << mg << ", nP=" << nP << "\n";
00131           dbg->flush();
00132 #endif
00133           //std. error -> output
00134           ex[nExp].errP=dvector(1,nP);
00135           ex[nExp].errY0=dvector(1,nvar);         
00136           // initial guesses at mesh points
00137           ex[nExp].yTry=dmatrix(1,nPoints,1,nvar); 
00138           ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar); 
00139           // computed values at mesh points   
00140           ex[nExp].yComp=dmatrix(1,nPoints,1,nvar);
00141           // computed values at measuring points
00142           ex[nExp].yPred=dmatrix(1,nMeas,1,nobs);
00143           
00144           ex[nExp].h=dmatrix(1,nPoints,1,nvar);      // discrepancies
00145           ex[nExp].residues=dmatrix(1,nMeas,1,nobs); // residues
00146           
00147           // derivatives
00148           // d(yComp[i][j]) / d(yTry[i-1][k])
00149           ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar);
00150           // d(yComp[i][j]) / d(p[k])
00151           ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP);
00152           // d(yPred[i][j]) / d(yTry[**][k])
00153           ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar);
00154           // d(yPred[i][j]) / d(p[k])
00155           ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP);
00156           if (me>0) 
00157             {
00158               ex[nExp].r2=dvector(1,me);                 // equal. constr.
00159               // d(yR2[i]) / d(yTry[j][k])
00160               ex[nExp].dR2ds=d3tensor(1,me,1,nPoints,1,nvar);
00161               for(i=1;i<=me;i++)
00162                 for(j=1;j<=nPoints;j++)
00163                   for(k=1;k<=nvar;k++)
00164                     ex[nExp].dR2ds[i][j][k]=0.;
00165               
00166               // d(yR2[i]) / d(p[j])
00167               ex[nExp].dR2dp=dmatrix(1,me,1,nP);
00168               for(i=1;i<=me;i++)
00169                 for(j=1;j<=nP;j++)
00170                   ex[nExp].dR2dp[i][j]=0.;
00171             } 
00172           else 
00173             {
00174               ex[nExp].r2=NULL; ex[nExp].dR2ds=NULL; ex[nExp].dR2dp=NULL;
00175             }
00176           if (mg>0) 
00177             {
00178               ex[nExp].r3=dvector(1,mg);                 // ineq.  constr.
00179               // d(yR3[i]) / d(yTry[j][k])
00180               ex[nExp].dR3ds=d3tensor(1,mg,1,nPoints,1,nvar);
00181               for(i=1;i<=mg;i++)
00182                 for(j=1;j<=nPoints;j++)
00183                   for(k=1;k<=nvar;k++)
00184                     ex[nExp].dR3ds[i][j][k]=0.;
00185               // d(yR3[i]) / d(p[j])
00186               ex[nExp].dR3dp=dmatrix(1,mg,1,nP);
00187               for(i=1;i<=mg;i++)
00188                 for(j=1;j<=nP;j++)
00189                   ex[nExp].dR3dp[i][j]=0.;
00190             } 
00191           else 
00192             {
00193               ex[nExp].r3=NULL; ex[nExp].dR3ds=NULL; ex[nExp].dR3dp=NULL;
00194             }
00195           
00196           // compute yTry from ex (all ex read from file still present)
00197           setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0);
00198           // check if initial values are compatible with constraints
00199           // CURRENTLY DISABLED (1. Dez. 2004)
00200 //        if (!initValuesOK(globs,&ex[nExp])) 
00201 //          {
00202 //            cerr << "Experiment #" << nExp 
00203 //                 << ": initial values are not compatible with constraints\n";
00204 //            cerr << "R2:";
00205 //            for (i=1;i<=ex[nExp].me;++i) 
00206 //              cerr << "\t" << ex[nExp].r2[i];
00207 //            cerr << "\nR3:";
00208 //            for (i=1;i<=ex[nExp].mg;++i) 
00209 //              cerr << "\t" << ex[nExp].r3[i];
00210 //            cerr << "\n";
00211 //            exit(1);
00212 //          }
00213           
00214 #ifdef PRINTEX
00215           *dbg << "Ex used for fitting: \n\n";
00216           for (i=1; i<=ex[nExp].nMeasure; ++i) 
00217             {
00218               *dbg << ex[nExp].xMeasure[i];
00219               for (j=1; j<=ex[nExp].nobs; ++j)
00220                 *dbg << "\t" << ex[nExp].yMeasure[i][j] << "\t" << ex[nExp].sigma[i][j];
00221               *dbg << '\n';
00222             }
00223           *dbg << '\n';
00224           dbg->flush();
00225 #endif
00226           //Allocate memory for condensation
00227           //least squares
00228           ex[nExp].ua=dvector(1,nMeas*nobs);
00229           ex[nExp].Ea=dmatrix(1,nMeas*nobs,1,nvar);
00230           ex[nExp].Pa=dmatrix(1,nMeas*nobs,1,nP);
00231           //equality constraints
00232           ex[nExp].ue=dvector(1,me);
00233           ex[nExp].Ee=dmatrix(1,me,1,nvar);
00234           ex[nExp].Pe=dmatrix(1,me,1,nP);
00235           //inequality constraints
00236           ex[nExp].ug=dvector(1,mg);
00237           ex[nExp].Eg=dmatrix(1,mg,1,nvar);
00238           ex[nExp].Pg=dmatrix(1,mg,1,nP);
00239           
00240           //update steps
00241           ex[nExp].dS=dmatrix(1,nPoints,1,nvar);
00242           ex[nExp].dP=dvector(1,nP);
00243         }// matches: for(nExp=1;nExp<=nrExp;++nExp) loop over experiments
00244       //covar in solvLin allocated
00245       globs->covar=NULL;
00246     } // matches with if(simit==FALSE)...
00247 
00248   //Initialise GNUplot
00249   if(globs->noGnu==FALSE)
00250     setupGnuFp(globs,ex);
00251   
00252   globs->cond=0;
00253   globs->Lambda=1;
00254   
00255   char line[1000];
00256   //Initialise Splines
00257   if(globs->initSpline==TRUE)
00258     {
00259       for (nExp=1;nExp<=nrExp;++nExp) 
00260         {
00261           ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00262           ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00263           ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00264           ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
00265           
00266           for(i=1;i<=NSPLINES;i++)
00267             {
00268               j=0;
00269               while(ex[nExp].splineFile[i][j]!='\0')
00270                 {
00271                   line[j]=ex[nExp].splineFile[i][j];
00272                   j++;
00273                 }
00274               line[j]='\0';
00275               inSplines.open(line);
00276               if(inSplines==NULL)
00277                 {
00278                   cerr << "Cannot open " << ex[nExp].splineFile[i] << ".\n";
00279                   exit(1);
00280                 }
00281               inSplines >>  ex[nExp].nNodes[i];
00282               ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00283               ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00284               ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00285               for(j=1;j<= ex[nExp].nNodes[i];j++)
00286                 {
00287                   inSplines >>  ex[nExp].splineNodes[i][j];
00288                   inSplines >>  ex[nExp].splineY[i][j];
00289                   inSplines >>  ex[nExp].splineGam[i][j];             
00290                 }
00291               inSplines.close();
00292             }
00293         }
00294     }
00295 }

main int  argc,
char *  argv[]
 

Definition at line 31 of file simit.cc.

References dbg, and parseopts().

00032 {
00033   long k,i,j;
00034   char outstr[100];
00035   Glob globs;
00036   GlobExp *ex;
00037   double *t,time,**y,*state;
00038   long n;
00039 
00040   ex=parseopts(argc,argv,&globs,outstr);
00041 
00042   //open DEBUG stream
00043   dbg=new ofstream("simit.dbg");
00044   if (!dbg) 
00045     {
00046       cerr << "Error opening DEBUG file.\n";
00047       exit(1);
00048     }
00049   dbg->precision(4);     // set output format for floating point data
00050   dbg->unsetf(ios::floatfield);
00051 
00052   //print debugging information
00053   *dbg << DefModelDescription << "\n";
00054   *dbg << NPARAMS << " parameter(s):";
00055   for (k=1; k<=NPARAMS; ++k) 
00056     {
00057       *dbg << " " << ex[1].par[k];
00058     }
00059   *dbg << "\n";
00060   dbg->flush();
00061  
00062    
00063   //initial values
00064   *dbg << "\n\n";
00065   if (ex[1].y0) 
00066     {
00067       *dbg << NEQNS << " starting value(s):";
00068       for (k=1; k<=NEQNS; ++k) 
00069         *dbg << ' ' << ex[1].y0[k];
00070       *dbg << "\n";
00071     } 
00072   else 
00073     {
00074       *dbg << "No starting values specified\n";
00075       ex[1].y0=dvector(1,ex[1].nvar);
00076       for (k=1; k<=NEQNS; ++k) 
00077         ex[1].y0[k]=DefYValues[k-1];
00078     }
00079   *dbg << "\n";
00080   dbg->flush();
00081 
00082   //setting up time points
00083   if(globs.dt <= 0.)
00084     {
00085       cerr << "illegal integration step: " << globs.dt <<  " \n";
00086       exit(1);
00087     }
00088   *dbg << "Integration step: " << globs.dt << endl;
00089 
00090   time=ex[1].fitstart;
00091   n=0;
00092   
00093   while(time < ex[1].fitend-globs.eps)
00094     {
00095       time+=globs.dt;
00096       n++;
00097     }
00098   n++;
00099   *dbg << "Number of time points : " << n << endl;
00100   
00101   t=dvector(1,n);
00102   t[1]=ex[1].fitstart;
00103   t[n]=ex[1].fitend;
00104 
00105   for(k=2;k<=n-1;k++)
00106     t[k]=t[k-1]+globs.dt;
00107 #ifdef PLOTTIMEPOINTS
00108   *dbg << "Time points : ";
00109   for(k=1;k<=n;k++)
00110     {
00111       *dbg << t[k];
00112       if(k!=n)
00113         *dbg << ";";
00114     }
00115   *dbg << endl;
00116 #endif
00117   dbg->flush();
00118 
00119   y=dmatrix(1,n,1,ex[1].nvar);
00120   state=dvector(1,ex[1].nvar);
00121   for(k=1;k<=ex[1].nvar;k++)
00122     state[k]=ex[1].y0[k];
00123 
00124   initialise(ex,&globs,TRUE);
00125 
00126   try
00127     {
00128       tabulateValues(&globs,&ex[1],t[1],t[n],t,n,state,y,NULL,NULL,NULL,NULL);
00129     }
00130   catch (int i)
00131     {
00132       exit(1);
00133     }
00134 
00135   outSimit(ex,&globs,t,n,y);
00136   cout << "creating " << ex[1].fileName << endl;
00137 
00138   freeMem(ex,&globs,TRUE);
00139 }

void outSimit GlobExp  ex[],
Glob globs,
double *  t,
long  n,
double **  y
 

Definition at line 12 of file outSimit.cc.

00013 {
00014   long k,l;
00015   long idum=-time(NULL);
00016   long nobs=ex[1].nobs;
00017   ofstream out;
00018   
00019   out.open(ex[1].fileName);
00020 
00021   if(globs->sig==0)
00022     {
00023       for(k=1;k<=n;k++)
00024         {
00025           out << t[k];
00026           for(l=1;l<=nobs;l++)
00027             out << " " << y[k][l] << " 1";
00028           out << endl;
00029         }
00030 
00031     }
00032   else
00033     {
00034       double *max=dvector(1,nobs),*min=dvector(1,nobs);
00035       for(l=1;l<=nobs;l++)
00036         {
00037           max[l]=y[1][l];
00038           min[l]=y[1][l];
00039         }
00040       for(k=2;k<=n;k++)
00041         {
00042           for(l=1;l<=nobs;l++)
00043             {
00044               if(max[l]<y[k][l])
00045                 max[l]=y[k][l];
00046               if(min[l]>y[k][l])
00047                 min[l]=y[k][l];
00048             }
00049         }
00050       for(k=1;k<=n;k++)
00051         {
00052           out << t[k] ;
00053           for(l=1;l<=nobs;l++)
00054             out << " " << y[k][l]+globs->sig*fabs(max[l]-min[l])*gasdev(&idum) << " " << globs->sig*fabs(max[l]-min[l]);
00055           out << endl;
00056         }
00057     }
00058 
00059 }

GlobExp* parseopts int  argc,
char *  argv[],
Glob globs,
char *  outstr
 

Main module routine doing all the pasing.

Parameters:
argc The number of command line arguments as passed by main.
*argv[] String of the arguments.
*globs The structure Glob.
*outstr Depreciated, not used.

Definition at line 120 of file parse.cc.

00121 {
00122   int longindex,opt;
00123   int parlistSpecified=FALSE;
00124   long k,l,nExp,sp;
00125   char name[MAXLENGTH],line[MAXLENGTH];
00126 
00127   // option definition
00128   static struct option longopts[]={
00129     {"help", 0, 0,  'h'},
00130     {"help", 0, 0,  '?'},
00131     {"x0"  , 1, 0,   1 },
00132     {"x1"  , 1, 0,   2 },
00133     {"nms"  ,1, 0,   3 },
00134     {"f"    ,1, 0,   4 },
00135     {"p"    ,1, 0,   5 },
00136     {"y0"   ,1, 0,   6 },
00137     {"nobs" ,1, 0,   7 },
00138     {"eps"  ,1, 0,   8 },
00139     {"nognu",0, 0,   9 },
00140     {"nomeasure",0,0,10},
00141     {"doP"  ,1, 0,   11},
00142     {"maxit",1, 0,   12},
00143     {"wait" ,0, 0,   13},
00144     {"usesig",0,0,   14},
00145     {"int"   ,1,0,   15},
00146     {"maxstp",1,0,   17},
00147     {"minimp",1,0,   18},
00148     {"nowait",0,0,   19},
00149     {"elastic" ,1,0, 21},
00150     {"reg"     ,0,0, 22},
00151     {"epsilon" ,1,0, 23},
00152     {"lambda"  ,1,0, 24},
00153     {"spline"  ,1,0, 25},
00154     {"siminit" ,0,0, 26},
00155     {"pert"    ,1,0, 27},
00156     {"y0fix"   ,1,0, 28},
00157     {"nodamp"  ,0,0, 29},
00158     {"strat"   ,1,0, 30},
00159     {"opt"     ,1,0, 31},
00160     {"Lconst"  ,1,0, 32},
00161     {0, 0, 0, 0}
00162   };
00163   //initialise some global parameters in globs
00164   globs->noGnu=FALSE;
00165   globs->eps=1e-6;
00166   globs->npar=NPARAMS;
00167   globs->noMeasurements=FALSE;
00168   globs->doP=ivector(1,globs->npar);
00169   for(k=1;k<=globs->npar;k++)
00170     globs->doP[k]=TRUE;
00171   globs->y0fix=ivector(1,NEQNS);
00172   for(k=1;k<=NEQNS;k++)
00173     globs->y0fix[k]=TRUE;
00174   globs->maxit=1000;
00175   globs->gnuFp=NULL;
00176   globs->wait=FALSE;
00177   globs->usesig=FALSE;
00178   globs->integrator=1;
00179   globs->stiff=TRUE;
00180   globs->maxstp=5000;
00181   globs->minimp=1e-4;
00182   globs->nowait=FALSE;
00183   globs->elastic=1.;
00184   globs->reg=FALSE;
00185   globs->epsilon=1e-10;
00186   globs->lambda=1e6;
00187   globs->simInit=FALSE;
00188   globs->pert=0.;
00189   globs->nodamp=FALSE;
00190   globs->initSpline=TRUE;
00191   globs->wquer=-1;
00192   globs->silent=FALSE;
00193   globs->strategy=1;
00194   globs->minimiser=2;
00195   globs->faktorL=dvector(1,NPARAMS);
00196   globs->faktorLexist=FALSE;
00197   for(k=1;k<=NPARAMS;k++)
00198   {
00199         globs->faktorL[k]=-1;
00200   }     
00201   while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00202     {
00203       switch(opt)
00204         {
00205         case 4:
00206           parlistSpecified=TRUE;
00207           strncpy(name,optarg,MAXLENGTH);
00208           break;
00209         }
00210     }
00211   optind=0;
00212   
00213   GlobExp *ex;
00214   if(parlistSpecified==FALSE)
00215     {
00216       globs->nrExp=1;
00217       ex=new GlobExp[globs->nrExp+1];   
00218       //initialise some global parameters
00219       ex[1].fitstart=-1e99;
00220       ex[1].fitend=1e99;
00221       ex[1].nobs=NOBS;
00222       ex[1].nvar=NEQNS;
00223       ex[1].nPoints=2;
00224       ex[1].splinesDefined=FALSE;
00225       ex[1].y0=NULL;
00226       // initialisation of the parameters
00227       ex[1].par=dvector(1,NPARAMS);
00228       for(k=1;k<=globs->npar;k++)
00229         ex[1].par[k]=DefParameters[k-1];
00230       
00231       if(NSPLINES!=0)
00232         {
00233           ex[1].nNodes=lvector(1,NSPLINES);
00234         }
00235     }
00236   else //if parlist specified
00237     {
00238       ifstream in;      ofstream out;
00239       char str[MAXLENGTH];
00240 
00241       in.open(name);
00242 
00243       if(!in)
00244         {
00245           cerr << "parameter list " << name << "  not found\n";
00246           exit(1);
00247         }
00248       
00249       //preprocessing of data list
00250       globs->nrExp=0;
00251       sprintf(str,"%s.%d",name,globs->nrExp);
00252       out.open(str);
00253       while(!in.eof())
00254         {
00255           in.getline(str,MAXLENGTH,'\n');
00256           if(str[0]=='@')
00257             {
00258               out.flush();
00259               out.close();
00260               globs->nrExp++;
00261               sprintf(str,"%s.%d",name,globs->nrExp);
00262               out.open(str);
00263             }
00264           else if(str[0]!='#')
00265             {
00266               out << str << endl;
00267             }
00268         }
00269       out.flush();
00270       out.close();
00271       in.close();
00272 
00273       if(globs->nrExp==0)
00274         {
00275           cerr << "insuffient number of experiments, @ missing\n";
00276           exit(1);
00277         }
00278       
00279       long _argc;
00280       char *_argv[MAXLENGTH];
00281 
00282       for(k=0;k<=MAXLENGTH;k++)
00283         _argv[k]=(char *)malloc((size_t) (150*sizeof(char)));
00284 
00285       {
00286         ifstream inp;
00287         
00288         _argc=0;
00289         sprintf(str,"%s.0",name);
00290         inp.open(str);
00291         
00292         
00293         while(!inp.eof())
00294           {
00295             inp >> _argv[_argc+1];
00296             _argc++;
00297           }
00298         inp.close();
00299       }
00300    
00301       //PARSING LIST
00302 
00303       while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00304         {
00305           switch(opt) 
00306             { 
00307             case 'h':
00308               cout << usage << endl;
00309               exit(1);
00310               break;
00311             case '?':
00312               cerr << usage << endl;
00313               exit(1);  
00314               break;
00315             case 8:
00316               globs->eps=fabs(atof(optarg));
00317               break;
00318             case 9:
00319               globs->noGnu=TRUE;
00320               break;
00321             case 10:
00322               globs->noMeasurements=TRUE;
00323               break;
00324             case 11:
00325               if(strlen(optarg)!=globs->npar)
00326                 {
00327                   cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
00328                   exit(1);
00329                 }
00330               for(k=0;k < globs->npar;k++)
00331                 {
00332                   if(optarg[k]=='0')
00333                     globs->doP[k+1]=FALSE;
00334                   else if (optarg[k]=='L')
00335                     globs->doP[k+1]='L'; 
00336                   else
00337                     globs->doP[k+1]=TRUE;
00338                 }
00339               break;
00340             case 12:
00341               globs->maxit=abs(atol(optarg));
00342               break;
00343             case 13:
00344               globs->wait=TRUE;
00345               break;
00346             case 14:
00347               globs->usesig=TRUE;
00348               break;
00349             case 15:
00350               globs->integrator=abs(atoi(optarg));
00351               break;
00352             case 17:
00353               globs->maxstp=abs(atoi(optarg));
00354               break;
00355             case 18:
00356               globs->minimp=fabs(atof(optarg));
00357               break;
00358             case 19:
00359               globs->nowait=TRUE;
00360               break;
00361             case 21:
00362               globs->elastic=fabs(atof(optarg));
00363               if(globs->elastic > 1. || globs->elastic == 0.)
00364                 {
00365                   cerr << "Insufficient range of -elast <val> \n";
00366                   exit(1);
00367                 }
00368               break;
00369             case 22:
00370               globs->reg=TRUE;
00371               break;
00372             case 23:
00373               globs->epsilon=fabs(atof(optarg));
00374               break;
00375             case 24:
00376               globs->lambda=fabs(atof(optarg));
00377               break;
00378             case 26:
00379               globs->simInit=TRUE;
00380               break;
00381             case 27:
00382               globs->pert=fabs(atof(optarg));
00383               break;
00384             case 28:
00385               if(strlen(optarg)!=NEQNS)
00386                 {
00387                   cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
00388                   exit(1);
00389                 }
00390               for(k=0;k < NEQNS;k++)
00391                 {
00392                   if(optarg[k]=='0')
00393                     globs->y0fix[k+1]=FALSE;
00394                   else
00395                     globs->y0fix[k+1]=TRUE;
00396                 }
00397               break;
00398             case 29:
00399               globs->nodamp=TRUE;
00400               break; 
00401             case 30:
00402               globs->strategy=abs(atoi(optarg));
00403               break;
00404             case 31:
00405               globs->minimiser=abs(atoi(optarg));
00406               break;
00407             case 32:
00408               get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00409               globs->faktorLexist=TRUE;
00410               break;                  
00411             default:
00412               cerr << endl;
00413               cerr << "Parsing parameter list produced errors.\n\n";
00414               exit(1);
00415             }
00416         }
00417 
00418       optind=0;
00419       ex=new GlobExp[globs->nrExp+1];   
00420   
00421       //PARSING experiment specific list
00422       for(nExp=1;nExp<=globs->nrExp;nExp++)
00423         {
00424           //initialise some global parameters
00425           ex[nExp].fitstart=-1e99;
00426           ex[nExp].fitend=1e99;
00427           ex[nExp].nobs=NOBS;
00428           ex[nExp].nvar=NEQNS;
00429           ex[nExp].nPoints=2;
00430           ex[nExp].splinesDefined=FALSE;
00431           ex[nExp].y0=NULL;
00432           // initialisation of the parameters
00433           ex[nExp].par=dvector(1,NPARAMS);
00434           for(k=1;k<=globs->npar;k++)
00435             ex[nExp].par[k]=DefParameters[k-1];
00436           
00437           if(NSPLINES!=0)
00438             {
00439               ex[nExp].nNodes=lvector(1,NSPLINES);
00440             }
00441 
00442           _argc=0;
00443 
00444           ifstream inp;
00445           sprintf(str,"%s.%d",name,nExp);
00446           inp.open(str);
00447 
00448           while(!inp.eof())
00449             {
00450               inp >> _argv[_argc+1];
00451               _argc++;
00452             }
00453 
00454           inp.close();
00455           while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00456              {
00457                switch(opt)
00458                  {
00459                  case 1:
00460                    ex[nExp].fitstart=atof(optarg);
00461                    break;
00462                  case 2:
00463                    ex[nExp].fitend=atof(optarg);
00464                    break;
00465                  case 3:
00466                    ex[nExp].nPoints=abs(atol(optarg));
00467                    ex[nExp].nPoints+=1;
00468                    break;
00469                  case 5:
00470                    get_list(globs->npar,ex[nExp].par,optarg," parameter ");
00471                    break;
00472                  case 6:
00473                    ex[nExp].y0=dvector(1,ex[nExp].nvar);
00474                    for(k=1;k<=ex[nExp].nvar;k++)
00475                      ex[nExp].y0[k]=0.;
00476                    k=get_list(ex[nExp].nvar,ex[nExp].y0,optarg," initial values ");
00477                    if(k!=ex[nExp].nvar)
00478                      {
00479                        cerr << ex[nExp].nvar << " initial values required.\n";
00480                        exit(1);
00481                      }
00482                    break;
00483                  case 7:
00484                    ex[nExp].nobs=atol(optarg);
00485                    break;
00486                  case 25:
00487                    k=0;
00488                    l=1;
00489                    ex[nExp].splinesDefined=TRUE;
00490                    while(optarg[k]!='\0')
00491                      {
00492                        if(optarg[k]==',')
00493                          l++;
00494                        k++;
00495                      }
00496                    if(l!=NSPLINES)
00497                      {
00498                        cerr << "Incompatible number of splines.\n";
00499                        exit(1);
00500                      }
00501                    ex[nExp].splineFile=new string[l+1];
00502                    l=1;
00503                    sp=0;
00504                    k=0;
00505                    while(optarg[k]!='\0')
00506                      {
00507                        if(optarg[k]==',')
00508                          {
00509                            line[sp]='\0';
00510                            ex[nExp].splineFile[l]=(string)line;
00511                            sp=0;
00512                            l++;
00513                          }
00514                        else
00515                          {
00516                            line[sp]=optarg[k];
00517                            sp++;
00518                          }
00519                        k++;
00520                      }
00521                    line[sp]='\0';
00522                    ex[nExp].splineFile[l]=(string)line;
00523                    break;
00524                  default:
00525                    cerr << endl;
00526                    cerr << "Parsing parameter list produced errors.\n\n";
00527                    exit(1);
00528                  }
00529              }
00530 
00531            //Parsing file name
00532            if((_argc-optind)!=1)
00533              {
00534                cerr << "No/Too many datafile(s) specified for experiment " << nExp << " .  \n\n";
00535                exit(1);
00536              }
00537            else
00538              {
00539                ex[nExp].fileName=new char[strlen(_argv[optind])+1];
00540                strcpy(ex[nExp].fileName,_argv[optind]);
00541              }
00542            //checking if x0 < x1
00543            if(ex[nExp].fitstart >= ex[nExp].fitend)
00544              {
00545                cerr << "x0 must be smaller than x1.\n";
00546                exit(1);
00547              }
00548            optind=0;
00549         }
00550       
00551      
00552       //for(k=0;k<=MAXLENGTH;k++)
00553       //free(_argv[k]);
00554 
00555       sprintf(str,"rm -f %s.*",name);
00556       system(str);
00557     }
00558 
00559   optind=0;
00560   
00561   while((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00562     {
00563       switch(opt)
00564         {
00565         case 'h':
00566           cout << usage << endl;
00567           exit(1);
00568           break;
00569         case '?':
00570           cerr << usage << endl;
00571           exit(1);      
00572           break;
00573         case 1:
00574           ex[1].fitstart=atof(optarg);
00575           break;
00576         case 2:
00577           ex[1].fitend=atof(optarg);
00578           break;
00579         case 3:
00580           ex[1].nPoints=abs(atol(optarg));
00581           ex[1].nPoints+=1;
00582           break;
00583         case 4:
00584           break;
00585         case 5:
00586           get_list(globs->npar,ex[1].par,optarg," parameter ");
00587           break;
00588         case 6:
00589           ex[1].y0=dvector(1,ex[1].nvar);
00590           for(k=1;k<=ex[1].nvar;k++)
00591             ex[1].y0[k]=0.;
00592           k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
00593           if(k!=ex[1].nvar)
00594             {
00595               cerr << ex[1].nvar << " initial values required.\n";
00596               exit(1);
00597             }
00598           break;
00599         case 7:
00600           ex[1].nobs=atol(optarg);
00601           break;
00602         case 8:
00603           globs->eps=fabs(atof(optarg));
00604           break;
00605         case 9:
00606           globs->noGnu=TRUE;
00607           break;
00608         case 10:
00609           globs->noMeasurements=TRUE;
00610           break;
00611         case 11:
00612           if(strlen(optarg)!=globs->npar)
00613             {
00614               cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
00615               exit(1);
00616             }
00617           for(k=0;k < globs->npar;k++)
00618             {
00619               if(optarg[k]=='0')
00620                 globs->doP[k+1]=FALSE;
00621               else if (optarg[k]=='L')
00622                 globs->doP[k+1]='L'; 
00623               else
00624                 globs->doP[k+1]=TRUE;
00625             }
00626           break;
00627         case 12:
00628           globs->maxit=abs(atol(optarg));
00629           break;
00630         case 13:
00631           globs->wait=TRUE;
00632           break;
00633         case 14:
00634           globs->usesig=TRUE;
00635           break;
00636         case 15:
00637           globs->integrator=abs(atoi(optarg));
00638           break;
00639         case 17:
00640           globs->maxstp=abs(atoi(optarg));
00641           break;
00642         case 18:
00643           globs->minimp=fabs(atof(optarg));
00644           break;
00645         case 19:
00646           globs->nowait=TRUE;
00647           break;
00648         case 21:
00649           globs->elastic=fabs(atof(optarg));
00650           if(globs->elastic > 1. || globs->elastic == 0.)
00651             {
00652               cerr << "Insufficient range of -elast <val> \n";
00653               exit(1);
00654             }
00655           break;
00656         case 22:
00657           globs->reg=TRUE;
00658           break;
00659         case 23:
00660           globs->epsilon=fabs(atof(optarg));
00661           break;
00662         case 24:
00663           globs->lambda=fabs(atof(optarg));
00664           break;
00665         case 25:
00666           k=0;
00667           l=1;
00668           ex[1].splinesDefined=TRUE;
00669           while(optarg[k]!='\0')
00670             {
00671               if(optarg[k]==',')
00672                 l++;
00673               k++;
00674             }
00675           if(l!=NSPLINES)
00676             {
00677               cerr << "Incompatible number of splines.\n";
00678               exit(1);
00679             }
00680           ex[1].splineFile=new string[l+1];
00681           l=1;
00682           sp=0;
00683           k=0;
00684           while(optarg[k]!='\0')
00685             {
00686               if(optarg[k]==',')
00687                 {
00688                   line[sp]='\0';
00689                   ex[1].splineFile[l]=(string)line;
00690                   sp=0;
00691                   l++;
00692                 }
00693               else
00694                 {
00695                   line[sp]=optarg[k];
00696                   sp++;
00697                 }
00698               k++;
00699             }
00700           line[sp]='\0';
00701           ex[1].splineFile[l]=(string)line;
00702           break;           
00703         case 26:
00704           globs->simInit=TRUE;
00705           break;
00706         case 27:
00707           globs->pert=fabs(atof(optarg));
00708           break;
00709         case 28:
00710           if(strlen(optarg)!=NEQNS)
00711             {
00712               cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
00713               exit(1);
00714             }
00715           for(k=0;k < NEQNS;k++)
00716             {
00717               if(optarg[k]=='0')
00718                 globs->y0fix[k+1]=FALSE;
00719               else
00720                 globs->y0fix[k+1]=TRUE;
00721             }
00722           break;
00723         case 29:
00724           globs->nodamp=TRUE;
00725           break;
00726         case 30:
00727           globs->strategy=abs(atoi(optarg));
00728           break;
00729         case 31:
00730           globs->minimiser=abs(atoi(optarg));
00731           break;
00732         case 32:
00733           get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00734           globs->faktorLexist=TRUE;
00735           break;
00736         default:
00737           cerr << endl;
00738           cerr << "Parsing command line options produced errors \n\n";
00739           exit(1);
00740         }
00741     }
00742   //checking if x0 < x1
00743   if(ex[1].fitstart >= ex[1].fitend)
00744     {
00745       cerr << "x0 must be smaller than x1.\n";
00746       exit(1);
00747     }
00748   
00749   //Parsing file name for a single experiment
00750   if((argc-optind)!=1 && parlistSpecified==FALSE)
00751     {
00752       cerr << "No/Too many datafile(s) specified.  \n\n";
00753       exit(1);
00754     }
00755   else if(parlistSpecified==FALSE)
00756     {
00757       ex[1].fileName=new char[strlen(argv[optind])+1];
00758       strcpy(ex[1].fileName,argv[optind]);
00759     }
00760   
00761   if(NSPLINES!=0)
00762     {
00763       for(k=1;k<=globs->nrExp;k++)
00764         {
00765           if(ex[k].splinesDefined==FALSE)
00766             {
00767               cerr << "Please specify spline data in experiment " << k << ".\n";
00768               exit(1);
00769             }
00770         }
00771     }
00772   
00773   return(ex);
00774 }

void tabulateValues Glob globs,
GlobExp ex,
double  t0,
double  t1,
double *  t,
long  n,
double *  state,
double **  y,
double ***  dmds,
double ***  dmdp,
double **  dyds,
double **  dydp
 

Definition at line 108 of file intODE.cc.

00110 {
00111   // integrate from t0 to t1, storing values at t[1..n]
00112   // t[n] may be ==t1 or <t1
00113   // in the latter case we will integrate from t[n] to t1 in the
00114   // final loop k=n+1 without storing any data points
00115   // If dmds!=NULL, compute derivatives as well
00116   // and store them in dmds,dmdp,dyds,dydp
00117   // with respect to observation, see odin.doc: Observation function
00118   
00119   //BEGIN von Felix veraendert
00120   //   double hmin=max(fabs(t0),fabs(t1))*1e-8;
00121   double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
00122   // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest.
00123   // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern
00124   // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion.
00125   // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht
00126   // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich
00127   // von 1e-18 noch einigermassen weit entfernt bleibt. 
00128   // Man koennte das in Zukunft noch in einen externen Parameter stecken.
00129   //END von Felix veraendert
00130 
00131   long nPExp=globs->npar;
00132   long nVar=ex->nvar;
00133   long nobs=NOBS;
00134   int generic;
00135   double *p=ex->par;
00136   double eps=globs->eps;  
00137   double h1 = (t1 - t0) * .01;
00138   long nok,nbad;
00139 
00140   initInt(globs,ex);
00141   if (n > 0 && (t0 > t[1] || t1 < t[n]))
00142     {
00143       cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl;
00144       throw 1;
00145     }
00146 
00147   if (globs->integrator==3)
00148     {                           // said -odessa
00149 #ifdef ODESSA
00150       double &rtol = eps;
00151 //BEGIN von Felix veraendert
00152 //       double atol=eps*1e-6;  
00153       double atol = eps;
00154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!)
00155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu
00156 // siehe odessa.f: ATOL und RTOL
00157 //END von Felix veraendert
00158       // odeint uses TINY=1e-30 in a similar way as odessa's atol, 
00159       // but atol<1e-16 or so results in step size underflow
00160       
00161       int odeint_MAXSTP=globs->maxstp;
00162       int d_flag=1;
00163       int stif=globs->stiff;
00164       char doPStr[globs->npar+1];
00165       
00166       for(long i=1;i<=globs->npar;i++)
00167         doPStr[i]=(char)globs->doP[i];
00168       
00169       call_odessa(globs,ex,nVar, nPExp, doPStr,
00170          d_flag, odeint_MAXSTP, stif,
00171          TRUE, n, t, t0, t1, state, y,
00172          dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1);
00173 
00174       return;
00175 #else
00176       cerr << "-odessa used but not linked" << endl;
00177       throw 1;
00178 #endif // ODESSA
00179     }                           // if odessa
00180       double *gy=dvector(1,nobs);
00181       long i, j, k;
00182       double xx1 = t0, xx2;
00183       
00184       if (dmds)
00185         {                           // compute sensitivities
00186           // create extended state vector Y
00187           double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
00188           initYt (Yt[0], state);
00189           
00190           // chain rule for extra observables
00191           double **dgdy, **dgdp;
00192           
00193           dgdy = dmatrix (1, nobs, 1, nVar);
00194           dgdp = dmatrix (1, nobs, 1, nPExp);
00195           dfill0 (dgdy[1], nobs * nVar);
00196           dfill0 (dgdp[1], nobs * nPExp);
00197           
00198           
00199           for (k = 1; k <= n + 1; k++)
00200             {
00201               xx2 = (k > n) ? t1 :  // end point
00202                 t[k];               // data point
00203               
00204               globs->rkqs_ign = (nPExp + nVar) * nVar;
00205               if (xx2 > xx1 && nVar > 0)
00206                 {                 
00207                   if(globs->integrator==1)
00208                     {
00209                       //Runge-Kutta integation
00210                       integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,
00211                                   eps, h1, hmin, &nok, &nbad, sensDerivs);
00212                     }
00213                   else if(globs->integrator==2)
00214                     {
00215                       cerr << "CVODES not optimized skip integrator" << endl;
00216                       throw 1;
00217             
00218                       //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE);
00219                     }
00220                 }
00221               // write state to y and sensitivities to dm/dp and dm/ds
00222               if (k <= n)
00223                 {
00224                   dcopy (y[k], Yt[0], nVar);
00225                   generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp);
00226                   for(i=1;i<=nobs;i++)
00227                     y[k][i]=gy[i];
00228                   
00229                   for (i = 1; i <= nobs; i++)
00230                     {
00231                       for (j = 1; j <= nPExp; j++)
00232                         {
00233                           double &dest = dmdp[k][i][j];
00234                           if (!generic)
00235                             {
00236                               dest = dgdp[i][j];
00237                               for (long l = 1; l <= nVar; l++)
00238                                 dest += dgdy[i][l] * Yt[j][l];
00239                             }
00240                           else
00241                             dest = Yt[j][i];
00242                         }           // for j
00243                       for (j = 1; j <= nVar; j++)
00244                         {
00245                           double &dest = dmds[k][i][j];
00246                           if (!generic)
00247                             {
00248                               dest = 0;
00249                               for (long l = 1; l <= nVar; l++)
00250                                 dest += dgdy[i][l] * Yt[nPExp + j][l];
00251                             }
00252                           else
00253                             dest = Yt[nPExp + j][i];
00254                         }           // for j
00255                     }               // for i
00256                 }                   // if k<=n
00257               xx1 = xx2;
00258             }                       // for k
00259           
00260           // copy state and its derivatives to state, dydp and dyds
00261           dcopy (state, Yt[0], nVar);       // original state vector
00262           // write sensitivities to dy/dp and dy/ds
00263           for (i = 1; i <= nVar; i++)
00264             {
00265               for (j = 1; j <= nPExp; j++)
00266                 dydp[i][j] = Yt[j][i];
00267               for (j = 1; j <= nVar; j++)
00268                 dyds[i][j] = Yt[nPExp + j][i];
00269             }
00270           
00271           free_dmatrix (dgdp,1 ,nobs , 1, nPExp);
00272           free_dmatrix (dgdy,1 ,nobs, 1, nVar);
00273           free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar);
00274         }
00275       else
00276         {                           // !dmds
00277           for (k = 1; k <= n + 1; k++)
00278             {
00279               xx2 = (k <= n) ? t[k] : t1;
00280               
00281               globs->rkqs_ign = 0;
00282               if (xx2 > xx1 && nVar > 0)
00283                 { 
00284                   if(globs->integrator==1)
00285                     {
00286                       //Runge-Kutta integration
00287                       integrateRK(globs,state, p, nVar, xx1, xx2,
00288                                   eps, h1, hmin, &nok, &nbad, ode);
00289                     }
00290                   else if(globs->integrator==2)
00291                     {
00292                       cerr << "CVODES not optimized skip integrator" << endl;
00293                       throw 1;
00294 
00295                       //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE);
00296                     }
00297                 }
00298               if (k <= n)
00299                 {
00300                   dcopy (y[k], state, nVar);
00301                   generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL);
00302                   for(i=1;i<=nobs;i++)
00303                     y[k][i]=gy[i];
00304                 }
00305               xx1 = xx2;
00306             }                       // for k
00307         }           
00308       // if dmds else
00309       free_dvector(gy,1,nobs); 
00310 }


Variable Documentation

ofstream* dbg
 

debug stream (C++)

Definition at line 29 of file simit.cc.


Generated on Mon Jan 29 17:09:15 2007 for Diffit by  doxygen 1.4.6