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

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

Go to the source code of this file.

Defines

#define IN   prhs[0]
#define N   plhs[0]
#define YOUT   plhs[1]
#define DMDS   plhs[2]
#define DMDP   plhs[3]

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)
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

Variables

ofstream * dbg
 debug stream (C++)


Define Documentation

#define DMDP   plhs[3]
 

Definition at line 38 of file simit_mex.cc.

#define DMDS   plhs[2]
 

Definition at line 37 of file simit_mex.cc.

#define IN   prhs[0]
 

Definition at line 31 of file simit_mex.cc.

#define N   plhs[0]
 

Definition at line 35 of file simit_mex.cc.

Referenced by mexFunction().

#define YOUT   plhs[1]
 

Definition at line 36 of file simit_mex.cc.

Referenced by mexFunction().


Function Documentation

void freeMem GlobExp ex,
Glob globs,
int  simit
 

Definition at line 11 of file freeMem.cc.

References FALSE, free_d3tensor(), free_dmatrix(), free_dvector(), Glob::gnuFp, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, Glob::nrExp, and GlobExp::nvar.

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.

References d3tensor(), dbg, dmatrix(), GlobExp::dmdp, GlobExp::dmds, Glob::doP, GlobExp::dR2dp, GlobExp::dR2ds, GlobExp::dR3dp, GlobExp::dR3ds, dvector(), GlobExp::dydp, GlobExp::dyds, GlobExp::errP, GlobExp::errY0, FALSE, GlobExp::h, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, GlobExp::par, GlobExp::r2, GlobExp::residues, TRUE, GlobExp::yComp, GlobExp::yPred, GlobExp::yTry, and GlobExp::yTrySave.

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 }

void mexFunction int  nlhs,
mxArray *  plhs[],
int  nrhs,
const mxArray *  prhs[]
 

Definition at line 41 of file simit_mex.cc.

References dmatrix(), Glob::doP, Glob::dt, dvector(), Glob::elastic, Glob::eps, Glob::epsilon, FALSE, GlobExp::fitend, GlobExp::fitstart, Glob::gnuFp, IN, initialise(), Glob::initSpline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, Glob::maxstp, Glob::minimp, N, GlobExp::nNodes, GlobExp::nobs, Glob::noGnu, Glob::noMeasurements, Glob::nowait, Glob::npar, Glob::nrExp, GlobExp::nvar, GlobExp::par, Glob::reg, Glob::sig, SPLINE, GlobExp::splineGam, GlobExp::splineNodes, GlobExp::splineY, Glob::stiff, TRUE, Glob::usesig, Glob::wait, GlobExp::y0, and YOUT.

00046 {
00047   long k,i,j,n,ldum;
00048   char outstr[100];
00049   Glob globs;
00050   GlobExp *ex;
00051   double time,**y,*state;
00052   double t0; //init. time
00053   double *t,*y0,*p;
00054   double **yout,*yout_,**ydum;
00055   double ***dmds,***dmdp,*dmds_,*dmdp_;
00056   double **dyds,**dydp;
00057   double *n_,*spline_;
00058   double *tdum=dvector(1,1);
00059   mxArray *splinedat;
00060   mxArray *tmp;
00061   mxArray *SPLINE;
00062   double *tmp_;
00063   int hidden;
00064 
00065   // input -> initial time
00066   tmp=mxGetField(IN,0,"t");
00067   tmp_=mxGetPr(tmp);
00068   n=mxGetNumberOfElements(tmp);
00069   t=dvector(1,n);
00070   for(i=1;i<=n;i++)
00071     t[i]=tmp_[i-1];
00072  
00073   //allocation of init. state and parameters
00074   y0=dvector(1,NEQNS);
00075   p=dvector(1,NPARAMS);
00076 
00077   // input -> initial time
00078   tmp=mxGetField(IN,0,"t0");
00079   tmp_=mxGetPr(tmp);
00080   t0=tmp_[0];
00081 
00082   // input -> y0
00083   tmp=mxGetField(IN,0,"y0");
00084   tmp_=mxGetPr(tmp);
00085   if(mxGetNumberOfElements(tmp)!=NEQNS)
00086     {
00087       if(mxGetNumberOfElements(tmp)!=0)
00088         mexWarnMsgTxt("Insufficient number of initial values.\nUse predefined.\n");
00089       for(i=1;i<=NEQNS;i++)
00090         y0[i]=DefYValues[i-1];
00091     }
00092   else
00093     {
00094       for(i=1;i<=NEQNS;i++)
00095         y0[i]=tmp_[i-1];
00096     }
00097   
00098   // input -> p
00099   tmp=mxGetField(IN,0,"p");
00100   tmp_=mxGetPr(tmp);
00101   if(mxGetNumberOfElements(tmp)!=NPARAMS)
00102     {
00103       if(mxGetNumberOfElements(tmp)!=0)
00104         mexWarnMsgTxt("Insufficient number of parameters.\nUse predefined.\n");
00105       for(i=1;i<=NPARAMS;i++)
00106         p[i]=DefParameters[i-1];
00107     }  
00108   else
00109     {
00110       for(i=1;i<=NPARAMS;i++)
00111         p[i]=tmp_[i-1];
00112     }
00113 
00114   // input -> eps
00115   tmp=mxGetField(IN,0,"eps");
00116   tmp_=mxGetPr(tmp);
00117   globs.eps=fabs(tmp_[0]);
00118   
00119   // input -> integrator
00120   tmp=mxGetField(IN,0,"int");
00121   tmp_=mxGetPr(tmp);
00122   globs.integrator=abs((int)tmp_[0]);
00123   globs.stiff=TRUE;
00124   // input -> spline
00125   SPLINE=mxGetField(IN,0,"spline");
00126   
00127   // input -> hidden
00128   tmp=mxGetField(IN,0,"hidden");
00129   tmp_=mxGetPr(tmp);
00130   hidden=(int)tmp_[0];
00131 
00132   //allocation of trajectory vect. for output
00133   if(hidden==TRUE)
00134     YOUT=mxCreateDoubleMatrix(n,NOBS+NEQNS,mxREAL);
00135   else
00136     YOUT=mxCreateDoubleMatrix(n,NOBS,mxREAL);
00137   yout_=mxGetPr(YOUT);  
00138   yout=dmatrix(1,n,1,NOBS+NEQNS);
00139   ydum=dmatrix(1,2,1,NOBS+NEQNS);
00140 
00141   //initialise some global parameters in globs
00142   globs.noGnu=TRUE;
00143   globs.npar=NPARAMS;
00144   globs.noMeasurements=FALSE;
00145   globs.doP=ivector(1,globs.npar);
00146   for(k=1;k<=globs.npar;k++)
00147     globs.doP[k]=TRUE;
00148   globs.maxit=1000;
00149   globs.gnuFp=NULL;
00150   globs.wait=FALSE;
00151   globs.usesig=FALSE;
00152   globs.stiff=TRUE;
00153   globs.maxstp=5000;
00154   globs.minimp=0.05;
00155   globs.nowait=FALSE;
00156   globs.elastic=1.;
00157   globs.reg=FALSE;
00158   globs.epsilon=1e-10;
00159   globs.lambda=1e6;
00160   globs.dt=0.1;
00161   globs.sig=0.;
00162   globs.nrExp=1;  
00163 
00164   ex=new GlobExp[globs.nrExp+1];   
00165   //initialise some global parameters
00166 
00167   ex[1].nobs=NOBS;
00168   ex[1].nvar=NEQNS;
00169   ex[1].par=dvector(1,NPARAMS);
00170   ex[1].y0=dvector(1,NEQNS);
00171 
00172 
00173   for(k=1;k<=NPARAMS;k++)
00174     ex[1].par[k]=p[k];
00175 
00176  
00177   ex[1].fitstart=t[1];
00178   ex[1].fitend=t[n];
00179   //initial values
00180 
00181   state=dvector(1,ex[1].nvar);
00182   for(k=1;k<=NEQNS;k++)
00183     state[k]=y0[k];
00184 
00185 
00186   N=mxCreateDoubleMatrix(3,1,mxREAL);
00187   n_=mxGetPr(N);
00188   
00189   n_[0]=NEQNS;
00190   n_[1]=NOBS;
00191   n_[2]=NPARAMS;
00192   
00193   // spline section
00194   if(NSPLINES==0)
00195     {
00196       globs.initSpline=TRUE;
00197       initialise(ex,&globs,TRUE);
00198     }
00199   else
00200     {
00201       globs.initSpline=FALSE;      
00202       initialise(ex,&globs,TRUE);  
00203       ex[1].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00204       ex[1].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00205       ex[1].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00206       ex[1].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
00207 
00208       for(i=1;i<=NSPLINES;i++)
00209         {
00210           if(mxGetNumberOfElements(SPLINE) < NSPLINES)
00211             {
00212               cerr << "Please specify " << NSPLINES << " spline(s).\n";
00213               return;
00214             }
00215 
00216           splinedat=mxGetCell(SPLINE,i-1);
00217           spline_=mxGetPr(splinedat);
00218           ex[1].nNodes[i]=mxGetM(splinedat);
00219           if(mxGetN(splinedat)!=3)
00220             {
00221               cerr << "Invalid data for spline " << i << ".\n";
00222               return;
00223             }
00224           ex[1].splineNodes[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
00225           ex[1].splineY[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
00226           ex[1].splineGam[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));           
00227 
00228           //copy spline data
00229           for(j=1;j<=ex[1].nNodes[i];j++)
00230              {
00231                ex[1].splineNodes[i][j]=spline_[j-1];
00232                ex[1].splineY[i][j]=spline_[ex[1].nNodes[i]+j-1];
00233                ex[1].splineGam[i][j]=spline_[2*ex[1].nNodes[i]+j-1];        
00234              }
00235         }
00236     }
00237 
00238   //... and integrate
00239   try
00240     {
00241       if(nlhs!=2)
00242         {
00243           dmds=d3tensor(1,n,1,NEQNS,1,NEQNS);
00244           dmdp=d3tensor(1,n,1,NEQNS,1,NPARAMS);
00245           dyds=dmatrix(1,NEQNS,1,NEQNS);
00246           dydp=dmatrix(1,NEQNS,1,NPARAMS);
00247           
00248           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,dmds,dmdp,dyds,dydp);
00249           
00250           DMDS=mxCreateDoubleMatrix(NOBS*NEQNS,n,mxREAL);
00251           dmds_=mxGetPr(DMDS);
00252           DMDP=mxCreateDoubleMatrix(NOBS*NPARAMS,n,mxREAL);
00253           dmdp_=mxGetPr(DMDP);
00254           
00255 
00256           //copying 
00257           
00258           for(k=1;k<=NEQNS;k++) 
00259             for(i=1;i<=NOBS;i++)
00260               for(j=1;j<=n;j++)
00261                 {
00262                   dmds_[(j-1)*NOBS*NEQNS+(k-1)*NOBS+i-1]=dmds[j][i][k];
00263                 }
00264           
00265           
00266           for(k=1;k<=NPARAMS;k++)
00267             for(i=1;i<=NOBS;i++)
00268               for(j=1;j<=n;j++)
00269                 {
00270                   dmdp_[(j-1)*NOBS*NPARAMS+(k-1)*NOBS+i-1]=dmdp[j][i][k];
00271                 }
00272           
00273           free_d3tensor(dmds,1,n,1,NEQNS,1,NEQNS);
00274           free_d3tensor(dmdp,1,n,1,NEQNS,1,NPARAMS);
00275           free_dmatrix(dyds,1,NEQNS,1,NEQNS);
00276           free_dmatrix(dydp,1,NEQNS,1,NPARAMS);
00277         }
00278       else
00279         {
00280           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,NULL,NULL,NULL,NULL);
00281           
00282         }
00283 
00284       if(hidden==TRUE)
00285         {
00286           ldum=NOBS+NEQNS;
00287           for(i=1;i<=NEQNS;i++)
00288             state[i]=y0[i];
00289           tdum[1]=t[1];
00290           if(t0!=t[1])
00291             tabulateValues(&globs,&ex[1],t0,t[1],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
00292           for(j=1;j<=NEQNS;j++)
00293             yout[1][NOBS+j]=state[j];
00294 
00295           for(i=2;i<=n;i++)
00296             {
00297               tdum[1]=t[i]; 
00298               tabulateValues(&globs,&ex[1],t[i-1],t[i],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
00299               for(j=1;j<=NEQNS;j++)
00300                 yout[i][NOBS+j]=state[j];
00301              
00302             }
00303         }
00304       else
00305         ldum=NOBS;
00306 
00307       for(i=1;i<=ldum;i++)
00308         for(j=1;j<=n;j++)
00309           yout_[(i-1)*n+j-1]=yout[j][i];
00310     }
00311   catch (int i)
00312     {
00313       return;
00314     }
00315   //free mem
00316   freeMem(ex,&globs,TRUE);
00317   free_dvector(state,1,ex[1].nvar);
00318   free_dmatrix(yout,1,n,1,NEQNS+NOBS);
00319   free_dvector(t,1,n);
00320   free_dvector(y0,1,NEQNS);
00321   free_dvector(p,1,NPARAMS);
00322   free_dvector(ex[1].par,1,NPARAMS);
00323   free_dvector(ex[1].y0,1,NEQNS);
00324   free_dvector(tdum,1,1);
00325   free_dmatrix(ydum,1,2,1,NOBS+NEQNS);
00326 
00327   delete ex;
00328 }

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

Definition at line 12 of file outSimit.cc.

References GlobExp::nobs, and Glob::sig.

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.

References Glob::doP, Glob::dt, Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob::gnuFp, Glob::initSpline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, MAXLENGTH, Glob::maxstp, Glob::minimp, Glob::noGnu, Glob::noMeasurements, Glob::nowait, Glob::npar, Glob::nrExp, Glob::reg, Glob::sig, Glob::stiff, TRUE, Glob::usesig, and Glob::wait.

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_mex.cc.

Referenced by initialise(), main(), and objectiveF().


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