- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /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
GlobExp * | parseopts (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().