\section{/home/peifer/diffit/simit\_\-mex.cc File Reference}
\label{simit__mex_8cc}\index{/home/peifer/diffit/simit_mex.cc@{/home/peifer/diffit/simit\_\-mex.cc}}
{\tt \#include \char`\"{}mex.h\char`\"{}}\par
{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$stdlib.h$>$}\par
{\tt \#include $<$math.h$>$}\par
{\tt \#include \char`\"{}def.h\char`\"{}}\par
{\tt \#include \char`\"{}model.h\char`\"{}}\par
{\tt \#include \char`\"{}nr.h\char`\"{}}\par
\subsection*{Defines}
\begin{CompactItemize}
\item 
\#define {\bf IN}~prhs[0]
\item 
\#define {\bf N}~plhs[0]
\item 
\#define {\bf YOUT}~plhs[1]
\item 
\#define {\bf DMDS}~plhs[2]
\item 
\#define {\bf DMDP}~plhs[3]
\end{CompactItemize}
\subsection*{Functions}
\begin{CompactItemize}
\item 
{\bf Glob\-Exp} $\ast$ {\bf parseopts} (int argc, char $\ast$argv[$\,$], {\bf Glob} $\ast$globs, char $\ast$outstr)
\item 
void {\bf initialise} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs, int simit)
\item 
void {\bf tabulate\-Values} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, double t0, double t1, double $\ast$t, long n, double $\ast$state, double $\ast$$\ast$y, double $\ast$$\ast$$\ast$dmds, double $\ast$$\ast$$\ast$dmdp, double $\ast$$\ast$dyds, double $\ast$$\ast$dydp)
\item 
void {\bf out\-Simit} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs, double $\ast$t, long n, double $\ast$$\ast$y)
\item 
void {\bf free\-Mem} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs, int simit)
\item 
void {\bf mex\-Function} (int nlhs, mx\-Array $\ast$plhs[$\,$], int nrhs, const mx\-Array $\ast$prhs[$\,$])
\end{CompactItemize}
\subsection*{Variables}
\begin{CompactItemize}
\item 
ofstream $\ast$ {\bf dbg}
\begin{CompactList}\small\item\em debug stream (C++) \item\end{CompactList}\end{CompactItemize}


\subsection{Define Documentation}
\index{simit_mex.cc@{simit\_\-mex.cc}!DMDP@{DMDP}}
\index{DMDP@{DMDP}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define DMDP~plhs[3]}\label{simit__mex_8cc_ed9e5c2d3c49b69907af6932fc537c1f}




Definition at line 38 of file simit\_\-mex.cc.\index{simit_mex.cc@{simit\_\-mex.cc}!DMDS@{DMDS}}
\index{DMDS@{DMDS}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define DMDS~plhs[2]}\label{simit__mex_8cc_d177e0e99147f38a34ce540c29d2e91a}




Definition at line 37 of file simit\_\-mex.cc.\index{simit_mex.cc@{simit\_\-mex.cc}!IN@{IN}}
\index{IN@{IN}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define IN~prhs[0]}\label{simit__mex_8cc_c2bbd6d630a06a980d9a92ddb9a49928}




Definition at line 31 of file simit\_\-mex.cc.\index{simit_mex.cc@{simit\_\-mex.cc}!N@{N}}
\index{N@{N}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define N~plhs[0]}\label{simit__mex_8cc_0240ac851181b84ac374872dc5434ee4}




Definition at line 35 of file simit\_\-mex.cc.

Referenced by mex\-Function().\index{simit_mex.cc@{simit\_\-mex.cc}!YOUT@{YOUT}}
\index{YOUT@{YOUT}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define YOUT~plhs[1]}\label{simit__mex_8cc_9fd5fc158d6d9066d0ec47a1d646a8f9}




Definition at line 36 of file simit\_\-mex.cc.

Referenced by mex\-Function().

\subsection{Function Documentation}
\index{simit_mex.cc@{simit\_\-mex.cc}!freeMem@{freeMem}}
\index{freeMem@{freeMem}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void free\-Mem ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs}, int {\em simit})}\label{simit__mex_8cc_c39e0b1a98465f5380027ff10e613058}




Definition at line 11 of file free\-Mem.cc.

References FALSE, free\_\-d3tensor(), free\_\-dmatrix(), free\_\-dvector(), Glob::gnu\-Fp, Glob\-Exp::me, Glob\-Exp::mg, Glob\-Exp::n\-Measure, Glob\-Exp::nobs, Glob::npar, Glob::nr\-Exp, and Glob\-Exp::nvar.

\footnotesize\begin{verbatim}12 {
13   long nExp,i;
14   long nrExp=globs->nrExp;
15 
16 
17   delete globs->gnuFp;
18   
19   if(simit==FALSE)
20     {
21       for (nExp=1;nExp<=nrExp;++nExp) 
22         { 
23           long nPoints=ex[nExp].nPoints;
24           long nP=globs->npar;
25           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
26           long nMeas=ex[nExp].nMeasure;
27           long  me=ex[nExp].me, mg=ex[nExp].mg;
28           free_dvector(ex[nExp].errP,1,nP);
29           free_dvector(ex[nExp].errY0,1,nvar);
30           free_dmatrix(ex[nExp].yTry,1,nPoints,1,nvar);
31           free_dmatrix(ex[nExp].yTrySave,1,nPoints,1,nvar);
32           free_dmatrix(ex[nExp].yComp,1,nPoints,1,nvar);
33           free_dmatrix(ex[nExp].yPred,1,nMeas,1,nobs);
34           free_dmatrix(ex[nExp].h,1,nPoints,1,nvar);   
35           free_dmatrix(ex[nExp].residues,1,nMeas,1,nobs);
36           free_d3tensor(ex[nExp].dyds,1,nPoints,1,nvar,1,nvar);
37           free_d3tensor(ex[nExp].dydp,1,nPoints,1,nvar,1,nP);
38           free_d3tensor(ex[nExp].dmds,1,nMeas,1,nobs,1,nvar);
39           free_d3tensor(ex[nExp].dmdp,1,nMeas,1,nobs,1,nP);
40           if (me>0) 
41             {
42               free_dvector(ex[nExp].r2,1,me);
43               free_d3tensor(ex[nExp].dR2ds,1,me,1,nPoints,1,nvar);
44               free_dmatrix(ex[nExp].dR2dp,1,me,1,nP);
45             }
46           if (mg>0) 
47             {
48               free_dvector(ex[nExp].r3,1,mg);
49               free_d3tensor(ex[nExp].dR3ds,1,mg,1,nPoints,1,nvar);
50               free_dmatrix(ex[nExp].dR3dp,1,mg,1,nP);
51             }
52           
53           free_dvector(ex[nExp].ua,1,nMeas*nobs);
54           free_dmatrix(ex[nExp].Ea,1,nMeas*nobs,1,nvar);
55           free_dmatrix(ex[nExp].Pa,1,nMeas*nobs,1,nP);
56           free_dvector(ex[nExp].ue,1,me);
57           free_dmatrix(ex[nExp].Ee,1,me,1,nvar);
58           free_dmatrix(ex[nExp].Pe,1,me,1,nP);
59           free_dvector(ex[nExp].ug,1,mg);
60           free_dmatrix(ex[nExp].Eg,1,mg,1,nvar);
61           free_dmatrix(ex[nExp].Pg,1,mg,1,nP);
62           free_dmatrix(ex[nExp].dS,1,nPoints,1,nvar);
63           free_dvector(ex[nExp].dP,1,nP);
64           
65         }
66     }
67   for (nExp=1;nExp<=nrExp;++nExp) 
68     {
69       free(ex[nExp].splineNodes);
70       free(ex[nExp].splineY);
71       free(ex[nExp].splineGam);
72       free(ex[nExp].nNodes);
73       for(i=1;i<=NSPLINES;i++)
74         {
75           free(ex[nExp].splineNodes[i]);
76           free(ex[nExp].splineY[i]);
77           free(ex[nExp].splineGam[i]);
78         }
79     }
80   if(simit==FALSE && globs->covar!=NULL)
81     free_dmatrix(globs->covar,1,globs->fitdim,1,globs->fitdim);
82 }
\end{verbatim}\normalsize 


\index{simit_mex.cc@{simit\_\-mex.cc}!initialise@{initialise}}
\index{initialise@{initialise}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void initialise ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs}, int {\em simit})}\label{simit__mex_8cc_c66882342aff1604543a41022eb8a3f7}




Definition at line 98 of file initialise.cc.

References d3tensor(), dbg, dmatrix(), Glob\-Exp::dmdp, Glob\-Exp::dmds, Glob::do\-P, Glob\-Exp::d\-R2dp, Glob\-Exp::d\-R2ds, Glob\-Exp::d\-R3dp, Glob\-Exp::d\-R3ds, dvector(), Glob\-Exp::dydp, Glob\-Exp::dyds, Glob\-Exp::err\-P, Glob\-Exp::err\-Y0, FALSE, Glob\-Exp::h, Glob\-Exp::me, Glob\-Exp::mg, Glob\-Exp::n\-Measure, Glob\-Exp::nobs, Glob::npar, Glob\-Exp::n\-Points, Glob::nr\-Exp, Glob\-Exp::nvar, Glob\-Exp::par, Glob\-Exp::r2, Glob\-Exp::residues, TRUE, Glob\-Exp::y\-Comp, Glob\-Exp::y\-Pred, Glob\-Exp::y\-Try, and Glob\-Exp::y\-Try\-Save.

\footnotesize\begin{verbatim}99 {
100   long i,j,k,nExp;
101   long nrExp=globs->nrExp;
102   
103   ifstream inSplines;
104 
105   if(simit==FALSE)
106     {
107       for (nExp=1;nExp<=nrExp;++nExp) 
108         { 
109           
110           //non-local parameters are equal for each experiment
111           for(i=1;i<=globs->npar;i++)
112             {
113               if(globs->doP[i]==TRUE)
114                 ex[nExp].par[i]=ex[1].par[i];
115             }
116           
117           // set ex->me and ex->mg given in model.cc
118           setNrConstraints(&ex[nExp],globs->npar,ex[nExp].par);
119           
120           long nPoints=ex[nExp].nPoints;
121           long nP=globs->npar;
122           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
123           long nMeas=ex[nExp].nMeasure;
124           long  me=ex[nExp].me, mg=ex[nExp].mg;
125           
126 #ifdef PRINTINITVALUES
127           *dbg << "initialize: exp. #" << nExp << ", nPoints=" << nPoints
128                << ", nvar=" << nvar << ", nobs=" << nobs << "\nnMeas="
129                << nMeas << ", nVal=" << nVal << ", me=" << me
130                << ", mg=" << mg << ", nP=" << nP << "\n";
131           dbg->flush();
132 #endif
133           //std. error -> output
134           ex[nExp].errP=dvector(1,nP);
135           ex[nExp].errY0=dvector(1,nvar);         
136           // initial guesses at mesh points
137           ex[nExp].yTry=dmatrix(1,nPoints,1,nvar); 
138           ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar); 
139           // computed values at mesh points   
140           ex[nExp].yComp=dmatrix(1,nPoints,1,nvar);
141           // computed values at measuring points
142           ex[nExp].yPred=dmatrix(1,nMeas,1,nobs);
143           
144           ex[nExp].h=dmatrix(1,nPoints,1,nvar);      // discrepancies
145           ex[nExp].residues=dmatrix(1,nMeas,1,nobs); // residues
146           
147           // derivatives
148           // d(yComp[i][j]) / d(yTry[i-1][k])
149           ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar);
150           // d(yComp[i][j]) / d(p[k])
151           ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP);
152           // d(yPred[i][j]) / d(yTry[**][k])
153           ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar);
154           // d(yPred[i][j]) / d(p[k])
155           ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP);
156           if (me>0) 
157             {
158               ex[nExp].r2=dvector(1,me);                 // equal. constr.
159               // d(yR2[i]) / d(yTry[j][k])
160               ex[nExp].dR2ds=d3tensor(1,me,1,nPoints,1,nvar);
161               for(i=1;i<=me;i++)
162                 for(j=1;j<=nPoints;j++)
163                   for(k=1;k<=nvar;k++)
164                     ex[nExp].dR2ds[i][j][k]=0.;
165               
166               // d(yR2[i]) / d(p[j])
167               ex[nExp].dR2dp=dmatrix(1,me,1,nP);
168               for(i=1;i<=me;i++)
169                 for(j=1;j<=nP;j++)
170                   ex[nExp].dR2dp[i][j]=0.;
171             } 
172           else 
173             {
174               ex[nExp].r2=NULL; ex[nExp].dR2ds=NULL; ex[nExp].dR2dp=NULL;
175             }
176           if (mg>0) 
177             {
178               ex[nExp].r3=dvector(1,mg);                 // ineq.  constr.
179               // d(yR3[i]) / d(yTry[j][k])
180               ex[nExp].dR3ds=d3tensor(1,mg,1,nPoints,1,nvar);
181               for(i=1;i<=mg;i++)
182                 for(j=1;j<=nPoints;j++)
183                   for(k=1;k<=nvar;k++)
184                     ex[nExp].dR3ds[i][j][k]=0.;
185               // d(yR3[i]) / d(p[j])
186               ex[nExp].dR3dp=dmatrix(1,mg,1,nP);
187               for(i=1;i<=mg;i++)
188                 for(j=1;j<=nP;j++)
189                   ex[nExp].dR3dp[i][j]=0.;
190             } 
191           else 
192             {
193               ex[nExp].r3=NULL; ex[nExp].dR3ds=NULL; ex[nExp].dR3dp=NULL;
194             }
195           
196           // compute yTry from ex (all ex read from file still present)
197           setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0);
198           // check if initial values are compatible with constraints
199           // CURRENTLY DISABLED (1. Dez. 2004)
200 //        if (!initValuesOK(globs,&ex[nExp])) 
201 //          {
202 //            cerr << "Experiment #" << nExp 
203 //                 << ": initial values are not compatible with constraints\n";
204 //            cerr << "R2:";
205 //            for (i=1;i<=ex[nExp].me;++i) 
206 //              cerr << "\t" << ex[nExp].r2[i];
207 //            cerr << "\nR3:";
208 //            for (i=1;i<=ex[nExp].mg;++i) 
209 //              cerr << "\t" << ex[nExp].r3[i];
210 //            cerr << "\n";
211 //            exit(1);
212 //          }
213           
214 #ifdef PRINTEX
215           *dbg << "Ex used for fitting: \n\n";
216           for (i=1; i<=ex[nExp].nMeasure; ++i) 
217             {
218               *dbg << ex[nExp].xMeasure[i];
219               for (j=1; j<=ex[nExp].nobs; ++j)
220                 *dbg << "\t" << ex[nExp].yMeasure[i][j] << "\t" << ex[nExp].sigma[i][j];
221               *dbg << '\n';
222             }
223           *dbg << '\n';
224           dbg->flush();
225 #endif
226           //Allocate memory for condensation
227           //least squares
228           ex[nExp].ua=dvector(1,nMeas*nobs);
229           ex[nExp].Ea=dmatrix(1,nMeas*nobs,1,nvar);
230           ex[nExp].Pa=dmatrix(1,nMeas*nobs,1,nP);
231           //equality constraints
232           ex[nExp].ue=dvector(1,me);
233           ex[nExp].Ee=dmatrix(1,me,1,nvar);
234           ex[nExp].Pe=dmatrix(1,me,1,nP);
235           //inequality constraints
236           ex[nExp].ug=dvector(1,mg);
237           ex[nExp].Eg=dmatrix(1,mg,1,nvar);
238           ex[nExp].Pg=dmatrix(1,mg,1,nP);
239           
240           //update steps
241           ex[nExp].dS=dmatrix(1,nPoints,1,nvar);
242           ex[nExp].dP=dvector(1,nP);
243         }// matches: for(nExp=1;nExp<=nrExp;++nExp) loop over experiments
244       //covar in solvLin allocated
245       globs->covar=NULL;
246     } // matches with if(simit==FALSE)...
247 
248   //Initialise GNUplot
249   if(globs->noGnu==FALSE)
250     setupGnuFp(globs,ex);
251   
252   globs->cond=0;
253   globs->Lambda=1;
254   
255   char line[1000];
256   //Initialise Splines
257   if(globs->initSpline==TRUE)
258     {
259       for (nExp=1;nExp<=nrExp;++nExp) 
260         {
261           ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
262           ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
263           ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
264           ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
265           
266           for(i=1;i<=NSPLINES;i++)
267             {
268               j=0;
269               while(ex[nExp].splineFile[i][j]!='\0')
270                 {
271                   line[j]=ex[nExp].splineFile[i][j];
272                   j++;
273                 }
274               line[j]='\0';
275               inSplines.open(line);
276               if(inSplines==NULL)
277                 {
278                   cerr << "Cannot open " << ex[nExp].splineFile[i] << ".\n";
279                   exit(1);
280                 }
281               inSplines >>  ex[nExp].nNodes[i];
282               ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
283               ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
284               ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
285               for(j=1;j<= ex[nExp].nNodes[i];j++)
286                 {
287                   inSplines >>  ex[nExp].splineNodes[i][j];
288                   inSplines >>  ex[nExp].splineY[i][j];
289                   inSplines >>  ex[nExp].splineGam[i][j];             
290                 }
291               inSplines.close();
292             }
293         }
294     }
295 }
\end{verbatim}\normalsize 


\index{simit_mex.cc@{simit\_\-mex.cc}!mexFunction@{mexFunction}}
\index{mexFunction@{mexFunction}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void mex\-Function (int {\em nlhs}, mx\-Array $\ast$ {\em plhs}[$\,$], int {\em nrhs}, const mx\-Array $\ast$ {\em prhs}[$\,$])}\label{simit__mex_8cc_6a215cbfde54f82a3ce599228fc3fce5}




Definition at line 41 of file simit\_\-mex.cc.

References dmatrix(), Glob::do\-P, Glob::dt, dvector(), Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob\-Exp::fitend, Glob\-Exp::fitstart, Glob::gnu\-Fp, IN, initialise(), Glob::init\-Spline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, Glob::maxstp, Glob::minimp, N, Glob\-Exp::n\-Nodes, Glob\-Exp::nobs, Glob::no\-Gnu, Glob::no\-Measurements, Glob::nowait, Glob::npar, Glob::nr\-Exp, Glob\-Exp::nvar, Glob\-Exp::par, Glob::reg, Glob::sig, SPLINE, Glob\-Exp::spline\-Gam, Glob\-Exp::spline\-Nodes, Glob\-Exp::spline\-Y, Glob::stiff, TRUE, Glob::usesig, Glob::wait, Glob\-Exp::y0, and YOUT.

\footnotesize\begin{verbatim}46 {
47   long k,i,j,n,ldum;
48   char outstr[100];
49   Glob globs;
50   GlobExp *ex;
51   double time,**y,*state;
52   double t0; //init. time
53   double *t,*y0,*p;
54   double **yout,*yout_,**ydum;
55   double ***dmds,***dmdp,*dmds_,*dmdp_;
56   double **dyds,**dydp;
57   double *n_,*spline_;
58   double *tdum=dvector(1,1);
59   mxArray *splinedat;
60   mxArray *tmp;
61   mxArray *SPLINE;
62   double *tmp_;
63   int hidden;
64 
65   // input -> initial time
66   tmp=mxGetField(IN,0,"t");
67   tmp_=mxGetPr(tmp);
68   n=mxGetNumberOfElements(tmp);
69   t=dvector(1,n);
70   for(i=1;i<=n;i++)
71     t[i]=tmp_[i-1];
72  
73   //allocation of init. state and parameters
74   y0=dvector(1,NEQNS);
75   p=dvector(1,NPARAMS);
76 
77   // input -> initial time
78   tmp=mxGetField(IN,0,"t0");
79   tmp_=mxGetPr(tmp);
80   t0=tmp_[0];
81 
82   // input -> y0
83   tmp=mxGetField(IN,0,"y0");
84   tmp_=mxGetPr(tmp);
85   if(mxGetNumberOfElements(tmp)!=NEQNS)
86     {
87       if(mxGetNumberOfElements(tmp)!=0)
88         mexWarnMsgTxt("Insufficient number of initial values.\nUse predefined.\n");
89       for(i=1;i<=NEQNS;i++)
90         y0[i]=DefYValues[i-1];
91     }
92   else
93     {
94       for(i=1;i<=NEQNS;i++)
95         y0[i]=tmp_[i-1];
96     }
97   
98   // input -> p
99   tmp=mxGetField(IN,0,"p");
100   tmp_=mxGetPr(tmp);
101   if(mxGetNumberOfElements(tmp)!=NPARAMS)
102     {
103       if(mxGetNumberOfElements(tmp)!=0)
104         mexWarnMsgTxt("Insufficient number of parameters.\nUse predefined.\n");
105       for(i=1;i<=NPARAMS;i++)
106         p[i]=DefParameters[i-1];
107     }  
108   else
109     {
110       for(i=1;i<=NPARAMS;i++)
111         p[i]=tmp_[i-1];
112     }
113 
114   // input -> eps
115   tmp=mxGetField(IN,0,"eps");
116   tmp_=mxGetPr(tmp);
117   globs.eps=fabs(tmp_[0]);
118   
119   // input -> integrator
120   tmp=mxGetField(IN,0,"int");
121   tmp_=mxGetPr(tmp);
122   globs.integrator=abs((int)tmp_[0]);
123   globs.stiff=TRUE;
124   // input -> spline
125   SPLINE=mxGetField(IN,0,"spline");
126   
127   // input -> hidden
128   tmp=mxGetField(IN,0,"hidden");
129   tmp_=mxGetPr(tmp);
130   hidden=(int)tmp_[0];
131 
132   //allocation of trajectory vect. for output
133   if(hidden==TRUE)
134     YOUT=mxCreateDoubleMatrix(n,NOBS+NEQNS,mxREAL);
135   else
136     YOUT=mxCreateDoubleMatrix(n,NOBS,mxREAL);
137   yout_=mxGetPr(YOUT);  
138   yout=dmatrix(1,n,1,NOBS+NEQNS);
139   ydum=dmatrix(1,2,1,NOBS+NEQNS);
140 
141   //initialise some global parameters in globs
142   globs.noGnu=TRUE;
143   globs.npar=NPARAMS;
144   globs.noMeasurements=FALSE;
145   globs.doP=ivector(1,globs.npar);
146   for(k=1;k<=globs.npar;k++)
147     globs.doP[k]=TRUE;
148   globs.maxit=1000;
149   globs.gnuFp=NULL;
150   globs.wait=FALSE;
151   globs.usesig=FALSE;
152   globs.stiff=TRUE;
153   globs.maxstp=5000;
154   globs.minimp=0.05;
155   globs.nowait=FALSE;
156   globs.elastic=1.;
157   globs.reg=FALSE;
158   globs.epsilon=1e-10;
159   globs.lambda=1e6;
160   globs.dt=0.1;
161   globs.sig=0.;
162   globs.nrExp=1;  
163 
164   ex=new GlobExp[globs.nrExp+1];   
165   //initialise some global parameters
166 
167   ex[1].nobs=NOBS;
168   ex[1].nvar=NEQNS;
169   ex[1].par=dvector(1,NPARAMS);
170   ex[1].y0=dvector(1,NEQNS);
171 
172 
173   for(k=1;k<=NPARAMS;k++)
174     ex[1].par[k]=p[k];
175 
176  
177   ex[1].fitstart=t[1];
178   ex[1].fitend=t[n];
179   //initial values
180 
181   state=dvector(1,ex[1].nvar);
182   for(k=1;k<=NEQNS;k++)
183     state[k]=y0[k];
184 
185 
186   N=mxCreateDoubleMatrix(3,1,mxREAL);
187   n_=mxGetPr(N);
188   
189   n_[0]=NEQNS;
190   n_[1]=NOBS;
191   n_[2]=NPARAMS;
192   
193   // spline section
194   if(NSPLINES==0)
195     {
196       globs.initSpline=TRUE;
197       initialise(ex,&globs,TRUE);
198     }
199   else
200     {
201       globs.initSpline=FALSE;      
202       initialise(ex,&globs,TRUE);  
203       ex[1].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
204       ex[1].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
205       ex[1].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
206       ex[1].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
207 
208       for(i=1;i<=NSPLINES;i++)
209         {
210           if(mxGetNumberOfElements(SPLINE) < NSPLINES)
211             {
212               cerr << "Please specify " << NSPLINES << " spline(s).\n";
213               return;
214             }
215 
216           splinedat=mxGetCell(SPLINE,i-1);
217           spline_=mxGetPr(splinedat);
218           ex[1].nNodes[i]=mxGetM(splinedat);
219           if(mxGetN(splinedat)!=3)
220             {
221               cerr << "Invalid data for spline " << i << ".\n";
222               return;
223             }
224           ex[1].splineNodes[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
225           ex[1].splineY[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
226           ex[1].splineGam[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));           
227 
228           //copy spline data
229           for(j=1;j<=ex[1].nNodes[i];j++)
230              {
231                ex[1].splineNodes[i][j]=spline_[j-1];
232                ex[1].splineY[i][j]=spline_[ex[1].nNodes[i]+j-1];
233                ex[1].splineGam[i][j]=spline_[2*ex[1].nNodes[i]+j-1];        
234              }
235         }
236     }
237 
238   //... and integrate
239   try
240     {
241       if(nlhs!=2)
242         {
243           dmds=d3tensor(1,n,1,NEQNS,1,NEQNS);
244           dmdp=d3tensor(1,n,1,NEQNS,1,NPARAMS);
245           dyds=dmatrix(1,NEQNS,1,NEQNS);
246           dydp=dmatrix(1,NEQNS,1,NPARAMS);
247           
248           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,dmds,dmdp,dyds,dydp);
249           
250           DMDS=mxCreateDoubleMatrix(NOBS*NEQNS,n,mxREAL);
251           dmds_=mxGetPr(DMDS);
252           DMDP=mxCreateDoubleMatrix(NOBS*NPARAMS,n,mxREAL);
253           dmdp_=mxGetPr(DMDP);
254           
255 
256           //copying 
257           
258           for(k=1;k<=NEQNS;k++) 
259             for(i=1;i<=NOBS;i++)
260               for(j=1;j<=n;j++)
261                 {
262                   dmds_[(j-1)*NOBS*NEQNS+(k-1)*NOBS+i-1]=dmds[j][i][k];
263                 }
264           
265           
266           for(k=1;k<=NPARAMS;k++)
267             for(i=1;i<=NOBS;i++)
268               for(j=1;j<=n;j++)
269                 {
270                   dmdp_[(j-1)*NOBS*NPARAMS+(k-1)*NOBS+i-1]=dmdp[j][i][k];
271                 }
272           
273           free_d3tensor(dmds,1,n,1,NEQNS,1,NEQNS);
274           free_d3tensor(dmdp,1,n,1,NEQNS,1,NPARAMS);
275           free_dmatrix(dyds,1,NEQNS,1,NEQNS);
276           free_dmatrix(dydp,1,NEQNS,1,NPARAMS);
277         }
278       else
279         {
280           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,NULL,NULL,NULL,NULL);
281           
282         }
283 
284       if(hidden==TRUE)
285         {
286           ldum=NOBS+NEQNS;
287           for(i=1;i<=NEQNS;i++)
288             state[i]=y0[i];
289           tdum[1]=t[1];
290           if(t0!=t[1])
291             tabulateValues(&globs,&ex[1],t0,t[1],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
292           for(j=1;j<=NEQNS;j++)
293             yout[1][NOBS+j]=state[j];
294 
295           for(i=2;i<=n;i++)
296             {
297               tdum[1]=t[i]; 
298               tabulateValues(&globs,&ex[1],t[i-1],t[i],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
299               for(j=1;j<=NEQNS;j++)
300                 yout[i][NOBS+j]=state[j];
301              
302             }
303         }
304       else
305         ldum=NOBS;
306 
307       for(i=1;i<=ldum;i++)
308         for(j=1;j<=n;j++)
309           yout_[(i-1)*n+j-1]=yout[j][i];
310     }
311   catch (int i)
312     {
313       return;
314     }
315   //free mem
316   freeMem(ex,&globs,TRUE);
317   free_dvector(state,1,ex[1].nvar);
318   free_dmatrix(yout,1,n,1,NEQNS+NOBS);
319   free_dvector(t,1,n);
320   free_dvector(y0,1,NEQNS);
321   free_dvector(p,1,NPARAMS);
322   free_dvector(ex[1].par,1,NPARAMS);
323   free_dvector(ex[1].y0,1,NEQNS);
324   free_dvector(tdum,1,1);
325   free_dmatrix(ydum,1,2,1,NOBS+NEQNS);
326 
327   delete ex;
328 }
\end{verbatim}\normalsize 


\index{simit_mex.cc@{simit\_\-mex.cc}!outSimit@{outSimit}}
\index{outSimit@{outSimit}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void out\-Simit ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs}, double $\ast$ {\em t}, long {\em n}, double $\ast$$\ast$ {\em y})}\label{simit__mex_8cc_47814cfd9dc34050f085fe35b880db5f}




Definition at line 12 of file out\-Simit.cc.

References Glob\-Exp::nobs, and Glob::sig.

\footnotesize\begin{verbatim}13 {
14   long k,l;
15   long idum=-time(NULL);
16   long nobs=ex[1].nobs;
17   ofstream out;
18   
19   out.open(ex[1].fileName);
20 
21   if(globs->sig==0)
22     {
23       for(k=1;k<=n;k++)
24         {
25           out << t[k];
26           for(l=1;l<=nobs;l++)
27             out << " " << y[k][l] << " 1";
28           out << endl;
29         }
30 
31     }
32   else
33     {
34       double *max=dvector(1,nobs),*min=dvector(1,nobs);
35       for(l=1;l<=nobs;l++)
36         {
37           max[l]=y[1][l];
38           min[l]=y[1][l];
39         }
40       for(k=2;k<=n;k++)
41         {
42           for(l=1;l<=nobs;l++)
43             {
44               if(max[l]<y[k][l])
45                 max[l]=y[k][l];
46               if(min[l]>y[k][l])
47                 min[l]=y[k][l];
48             }
49         }
50       for(k=1;k<=n;k++)
51         {
52           out << t[k] ;
53           for(l=1;l<=nobs;l++)
54             out << " " << y[k][l]+globs->sig*fabs(max[l]-min[l])*gasdev(&idum) << " " << globs->sig*fabs(max[l]-min[l]);
55           out << endl;
56         }
57     }
58 
59 }
\end{verbatim}\normalsize 


\index{simit_mex.cc@{simit\_\-mex.cc}!parseopts@{parseopts}}
\index{parseopts@{parseopts}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}{\bf Glob\-Exp}$\ast$ parseopts (int {\em argc}, char $\ast$ {\em argv}[$\,$], {\bf Glob} $\ast$ {\em globs}, char $\ast$ {\em outstr})}\label{simit__mex_8cc_f7cc992592a2b95c9c3dec59ccf66584}


Main module routine doing all the pasing. \begin{Desc}
\item[Parameters:]
\begin{description}
\item[{\em argc}]The number of command line arguments as passed by {\bf main}. \item[{\em $\ast$argv\mbox{[}$\,$\mbox{]}}]String of the arguments. \item[{\em $\ast$globs}]The structure {\bf Glob}{\rm (p.\,\pageref{structGlob})}. \item[{\em $\ast$outstr}]Depreciated, not used. \end{description}
\end{Desc}


Definition at line 120 of file parse.cc.

References Glob::do\-P, Glob::dt, Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob::gnu\-Fp, Glob::init\-Spline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, MAXLENGTH, Glob::maxstp, Glob::minimp, Glob::no\-Gnu, Glob::no\-Measurements, Glob::nowait, Glob::npar, Glob::nr\-Exp, Glob::reg, Glob::sig, Glob::stiff, TRUE, Glob::usesig, and Glob::wait.

\footnotesize\begin{verbatim}121 {
122   int longindex,opt;
123   int parlistSpecified=FALSE;
124   long k,l,nExp,sp;
125   char name[MAXLENGTH],line[MAXLENGTH];
126 
127   // option definition
128   static struct option longopts[]={
129     {"help", 0, 0,  'h'},
130     {"help", 0, 0,  '?'},
131     {"x0"  , 1, 0,   1 },
132     {"x1"  , 1, 0,   2 },
133     {"nms"  ,1, 0,   3 },
134     {"f"    ,1, 0,   4 },
135     {"p"    ,1, 0,   5 },
136     {"y0"   ,1, 0,   6 },
137     {"nobs" ,1, 0,   7 },
138     {"eps"  ,1, 0,   8 },
139     {"nognu",0, 0,   9 },
140     {"nomeasure",0,0,10},
141     {"doP"  ,1, 0,   11},
142     {"maxit",1, 0,   12},
143     {"wait" ,0, 0,   13},
144     {"usesig",0,0,   14},
145     {"int"   ,1,0,   15},
146     {"maxstp",1,0,   17},
147     {"minimp",1,0,   18},
148     {"nowait",0,0,   19},
149     {"elastic" ,1,0, 21},
150     {"reg"     ,0,0, 22},
151     {"epsilon" ,1,0, 23},
152     {"lambda"  ,1,0, 24},
153     {"spline"  ,1,0, 25},
154     {"siminit" ,0,0, 26},
155     {"pert"    ,1,0, 27},
156     {"y0fix"   ,1,0, 28},
157     {"nodamp"  ,0,0, 29},
158     {"strat"   ,1,0, 30},
159     {"opt"     ,1,0, 31},
160     {"Lconst"  ,1,0, 32},
161     {0, 0, 0, 0}
162   };
163   //initialise some global parameters in globs
164   globs->noGnu=FALSE;
165   globs->eps=1e-6;
166   globs->npar=NPARAMS;
167   globs->noMeasurements=FALSE;
168   globs->doP=ivector(1,globs->npar);
169   for(k=1;k<=globs->npar;k++)
170     globs->doP[k]=TRUE;
171   globs->y0fix=ivector(1,NEQNS);
172   for(k=1;k<=NEQNS;k++)
173     globs->y0fix[k]=TRUE;
174   globs->maxit=1000;
175   globs->gnuFp=NULL;
176   globs->wait=FALSE;
177   globs->usesig=FALSE;
178   globs->integrator=1;
179   globs->stiff=TRUE;
180   globs->maxstp=5000;
181   globs->minimp=1e-4;
182   globs->nowait=FALSE;
183   globs->elastic=1.;
184   globs->reg=FALSE;
185   globs->epsilon=1e-10;
186   globs->lambda=1e6;
187   globs->simInit=FALSE;
188   globs->pert=0.;
189   globs->nodamp=FALSE;
190   globs->initSpline=TRUE;
191   globs->wquer=-1;
192   globs->silent=FALSE;
193   globs->strategy=1;
194   globs->minimiser=2;
195   globs->faktorL=dvector(1,NPARAMS);
196   globs->faktorLexist=FALSE;
197   for(k=1;k<=NPARAMS;k++)
198   {
199         globs->faktorL[k]=-1;
200   }     
201   while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
202     {
203       switch(opt)
204         {
205         case 4:
206           parlistSpecified=TRUE;
207           strncpy(name,optarg,MAXLENGTH);
208           break;
209         }
210     }
211   optind=0;
212   
213   GlobExp *ex;
214   if(parlistSpecified==FALSE)
215     {
216       globs->nrExp=1;
217       ex=new GlobExp[globs->nrExp+1];   
218       //initialise some global parameters
219       ex[1].fitstart=-1e99;
220       ex[1].fitend=1e99;
221       ex[1].nobs=NOBS;
222       ex[1].nvar=NEQNS;
223       ex[1].nPoints=2;
224       ex[1].splinesDefined=FALSE;
225       ex[1].y0=NULL;
226       // initialisation of the parameters
227       ex[1].par=dvector(1,NPARAMS);
228       for(k=1;k<=globs->npar;k++)
229         ex[1].par[k]=DefParameters[k-1];
230       
231       if(NSPLINES!=0)
232         {
233           ex[1].nNodes=lvector(1,NSPLINES);
234         }
235     }
236   else //if parlist specified
237     {
238       ifstream in;      ofstream out;
239       char str[MAXLENGTH];
240 
241       in.open(name);
242 
243       if(!in)
244         {
245           cerr << "parameter list " << name << "  not found\n";
246           exit(1);
247         }
248       
249       //preprocessing of data list
250       globs->nrExp=0;
251       sprintf(str,"%s.%d",name,globs->nrExp);
252       out.open(str);
253       while(!in.eof())
254         {
255           in.getline(str,MAXLENGTH,'\n');
256           if(str[0]=='@')
257             {
258               out.flush();
259               out.close();
260               globs->nrExp++;
261               sprintf(str,"%s.%d",name,globs->nrExp);
262               out.open(str);
263             }
264           else if(str[0]!='#')
265             {
266               out << str << endl;
267             }
268         }
269       out.flush();
270       out.close();
271       in.close();
272 
273       if(globs->nrExp==0)
274         {
275           cerr << "insuffient number of experiments, @ missing\n";
276           exit(1);
277         }
278       
279       long _argc;
280       char *_argv[MAXLENGTH];
281 
282       for(k=0;k<=MAXLENGTH;k++)
283         _argv[k]=(char *)malloc((size_t) (150*sizeof(char)));
284 
285       {
286         ifstream inp;
287         
288         _argc=0;
289         sprintf(str,"%s.0",name);
290         inp.open(str);
291         
292         
293         while(!inp.eof())
294           {
295             inp >> _argv[_argc+1];
296             _argc++;
297           }
298         inp.close();
299       }
300    
301       //PARSING LIST
302 
303       while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
304         {
305           switch(opt) 
306             { 
307             case 'h':
308               cout << usage << endl;
309               exit(1);
310               break;
311             case '?':
312               cerr << usage << endl;
313               exit(1);  
314               break;
315             case 8:
316               globs->eps=fabs(atof(optarg));
317               break;
318             case 9:
319               globs->noGnu=TRUE;
320               break;
321             case 10:
322               globs->noMeasurements=TRUE;
323               break;
324             case 11:
325               if(strlen(optarg)!=globs->npar)
326                 {
327                   cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
328                   exit(1);
329                 }
330               for(k=0;k < globs->npar;k++)
331                 {
332                   if(optarg[k]=='0')
333                     globs->doP[k+1]=FALSE;
334                   else if (optarg[k]=='L')
335                     globs->doP[k+1]='L'; 
336                   else
337                     globs->doP[k+1]=TRUE;
338                 }
339               break;
340             case 12:
341               globs->maxit=abs(atol(optarg));
342               break;
343             case 13:
344               globs->wait=TRUE;
345               break;
346             case 14:
347               globs->usesig=TRUE;
348               break;
349             case 15:
350               globs->integrator=abs(atoi(optarg));
351               break;
352             case 17:
353               globs->maxstp=abs(atoi(optarg));
354               break;
355             case 18:
356               globs->minimp=fabs(atof(optarg));
357               break;
358             case 19:
359               globs->nowait=TRUE;
360               break;
361             case 21:
362               globs->elastic=fabs(atof(optarg));
363               if(globs->elastic > 1. || globs->elastic == 0.)
364                 {
365                   cerr << "Insufficient range of -elast <val> \n";
366                   exit(1);
367                 }
368               break;
369             case 22:
370               globs->reg=TRUE;
371               break;
372             case 23:
373               globs->epsilon=fabs(atof(optarg));
374               break;
375             case 24:
376               globs->lambda=fabs(atof(optarg));
377               break;
378             case 26:
379               globs->simInit=TRUE;
380               break;
381             case 27:
382               globs->pert=fabs(atof(optarg));
383               break;
384             case 28:
385               if(strlen(optarg)!=NEQNS)
386                 {
387                   cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
388                   exit(1);
389                 }
390               for(k=0;k < NEQNS;k++)
391                 {
392                   if(optarg[k]=='0')
393                     globs->y0fix[k+1]=FALSE;
394                   else
395                     globs->y0fix[k+1]=TRUE;
396                 }
397               break;
398             case 29:
399               globs->nodamp=TRUE;
400               break; 
401             case 30:
402               globs->strategy=abs(atoi(optarg));
403               break;
404             case 31:
405               globs->minimiser=abs(atoi(optarg));
406               break;
407             case 32:
408               get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
409               globs->faktorLexist=TRUE;
410               break;                  
411             default:
412               cerr << endl;
413               cerr << "Parsing parameter list produced errors.\n\n";
414               exit(1);
415             }
416         }
417 
418       optind=0;
419       ex=new GlobExp[globs->nrExp+1];   
420   
421       //PARSING experiment specific list
422       for(nExp=1;nExp<=globs->nrExp;nExp++)
423         {
424           //initialise some global parameters
425           ex[nExp].fitstart=-1e99;
426           ex[nExp].fitend=1e99;
427           ex[nExp].nobs=NOBS;
428           ex[nExp].nvar=NEQNS;
429           ex[nExp].nPoints=2;
430           ex[nExp].splinesDefined=FALSE;
431           ex[nExp].y0=NULL;
432           // initialisation of the parameters
433           ex[nExp].par=dvector(1,NPARAMS);
434           for(k=1;k<=globs->npar;k++)
435             ex[nExp].par[k]=DefParameters[k-1];
436           
437           if(NSPLINES!=0)
438             {
439               ex[nExp].nNodes=lvector(1,NSPLINES);
440             }
441 
442           _argc=0;
443 
444           ifstream inp;
445           sprintf(str,"%s.%d",name,nExp);
446           inp.open(str);
447 
448           while(!inp.eof())
449             {
450               inp >> _argv[_argc+1];
451               _argc++;
452             }
453 
454           inp.close();
455           while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
456              {
457                switch(opt)
458                  {
459                  case 1:
460                    ex[nExp].fitstart=atof(optarg);
461                    break;
462                  case 2:
463                    ex[nExp].fitend=atof(optarg);
464                    break;
465                  case 3:
466                    ex[nExp].nPoints=abs(atol(optarg));
467                    ex[nExp].nPoints+=1;
468                    break;
469                  case 5:
470                    get_list(globs->npar,ex[nExp].par,optarg," parameter ");
471                    break;
472                  case 6:
473                    ex[nExp].y0=dvector(1,ex[nExp].nvar);
474                    for(k=1;k<=ex[nExp].nvar;k++)
475                      ex[nExp].y0[k]=0.;
476                    k=get_list(ex[nExp].nvar,ex[nExp].y0,optarg," initial values ");
477                    if(k!=ex[nExp].nvar)
478                      {
479                        cerr << ex[nExp].nvar << " initial values required.\n";
480                        exit(1);
481                      }
482                    break;
483                  case 7:
484                    ex[nExp].nobs=atol(optarg);
485                    break;
486                  case 25:
487                    k=0;
488                    l=1;
489                    ex[nExp].splinesDefined=TRUE;
490                    while(optarg[k]!='\0')
491                      {
492                        if(optarg[k]==',')
493                          l++;
494                        k++;
495                      }
496                    if(l!=NSPLINES)
497                      {
498                        cerr << "Incompatible number of splines.\n";
499                        exit(1);
500                      }
501                    ex[nExp].splineFile=new string[l+1];
502                    l=1;
503                    sp=0;
504                    k=0;
505                    while(optarg[k]!='\0')
506                      {
507                        if(optarg[k]==',')
508                          {
509                            line[sp]='\0';
510                            ex[nExp].splineFile[l]=(string)line;
511                            sp=0;
512                            l++;
513                          }
514                        else
515                          {
516                            line[sp]=optarg[k];
517                            sp++;
518                          }
519                        k++;
520                      }
521                    line[sp]='\0';
522                    ex[nExp].splineFile[l]=(string)line;
523                    break;
524                  default:
525                    cerr << endl;
526                    cerr << "Parsing parameter list produced errors.\n\n";
527                    exit(1);
528                  }
529              }
530 
531            //Parsing file name
532            if((_argc-optind)!=1)
533              {
534                cerr << "No/Too many datafile(s) specified for experiment " << nExp << " .  \n\n";
535                exit(1);
536              }
537            else
538              {
539                ex[nExp].fileName=new char[strlen(_argv[optind])+1];
540                strcpy(ex[nExp].fileName,_argv[optind]);
541              }
542            //checking if x0 < x1
543            if(ex[nExp].fitstart >= ex[nExp].fitend)
544              {
545                cerr << "x0 must be smaller than x1.\n";
546                exit(1);
547              }
548            optind=0;
549         }
550       
551      
552       //for(k=0;k<=MAXLENGTH;k++)
553       //free(_argv[k]);
554 
555       sprintf(str,"rm -f %s.*",name);
556       system(str);
557     }
558 
559   optind=0;
560   
561   while((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
562     {
563       switch(opt)
564         {
565         case 'h':
566           cout << usage << endl;
567           exit(1);
568           break;
569         case '?':
570           cerr << usage << endl;
571           exit(1);      
572           break;
573         case 1:
574           ex[1].fitstart=atof(optarg);
575           break;
576         case 2:
577           ex[1].fitend=atof(optarg);
578           break;
579         case 3:
580           ex[1].nPoints=abs(atol(optarg));
581           ex[1].nPoints+=1;
582           break;
583         case 4:
584           break;
585         case 5:
586           get_list(globs->npar,ex[1].par,optarg," parameter ");
587           break;
588         case 6:
589           ex[1].y0=dvector(1,ex[1].nvar);
590           for(k=1;k<=ex[1].nvar;k++)
591             ex[1].y0[k]=0.;
592           k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
593           if(k!=ex[1].nvar)
594             {
595               cerr << ex[1].nvar << " initial values required.\n";
596               exit(1);
597             }
598           break;
599         case 7:
600           ex[1].nobs=atol(optarg);
601           break;
602         case 8:
603           globs->eps=fabs(atof(optarg));
604           break;
605         case 9:
606           globs->noGnu=TRUE;
607           break;
608         case 10:
609           globs->noMeasurements=TRUE;
610           break;
611         case 11:
612           if(strlen(optarg)!=globs->npar)
613             {
614               cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
615               exit(1);
616             }
617           for(k=0;k < globs->npar;k++)
618             {
619               if(optarg[k]=='0')
620                 globs->doP[k+1]=FALSE;
621               else if (optarg[k]=='L')
622                 globs->doP[k+1]='L'; 
623               else
624                 globs->doP[k+1]=TRUE;
625             }
626           break;
627         case 12:
628           globs->maxit=abs(atol(optarg));
629           break;
630         case 13:
631           globs->wait=TRUE;
632           break;
633         case 14:
634           globs->usesig=TRUE;
635           break;
636         case 15:
637           globs->integrator=abs(atoi(optarg));
638           break;
639         case 17:
640           globs->maxstp=abs(atoi(optarg));
641           break;
642         case 18:
643           globs->minimp=fabs(atof(optarg));
644           break;
645         case 19:
646           globs->nowait=TRUE;
647           break;
648         case 21:
649           globs->elastic=fabs(atof(optarg));
650           if(globs->elastic > 1. || globs->elastic == 0.)
651             {
652               cerr << "Insufficient range of -elast <val> \n";
653               exit(1);
654             }
655           break;
656         case 22:
657           globs->reg=TRUE;
658           break;
659         case 23:
660           globs->epsilon=fabs(atof(optarg));
661           break;
662         case 24:
663           globs->lambda=fabs(atof(optarg));
664           break;
665         case 25:
666           k=0;
667           l=1;
668           ex[1].splinesDefined=TRUE;
669           while(optarg[k]!='\0')
670             {
671               if(optarg[k]==',')
672                 l++;
673               k++;
674             }
675           if(l!=NSPLINES)
676             {
677               cerr << "Incompatible number of splines.\n";
678               exit(1);
679             }
680           ex[1].splineFile=new string[l+1];
681           l=1;
682           sp=0;
683           k=0;
684           while(optarg[k]!='\0')
685             {
686               if(optarg[k]==',')
687                 {
688                   line[sp]='\0';
689                   ex[1].splineFile[l]=(string)line;
690                   sp=0;
691                   l++;
692                 }
693               else
694                 {
695                   line[sp]=optarg[k];
696                   sp++;
697                 }
698               k++;
699             }
700           line[sp]='\0';
701           ex[1].splineFile[l]=(string)line;
702           break;           
703         case 26:
704           globs->simInit=TRUE;
705           break;
706         case 27:
707           globs->pert=fabs(atof(optarg));
708           break;
709         case 28:
710           if(strlen(optarg)!=NEQNS)
711             {
712               cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
713               exit(1);
714             }
715           for(k=0;k < NEQNS;k++)
716             {
717               if(optarg[k]=='0')
718                 globs->y0fix[k+1]=FALSE;
719               else
720                 globs->y0fix[k+1]=TRUE;
721             }
722           break;
723         case 29:
724           globs->nodamp=TRUE;
725           break;
726         case 30:
727           globs->strategy=abs(atoi(optarg));
728           break;
729         case 31:
730           globs->minimiser=abs(atoi(optarg));
731           break;
732         case 32:
733           get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
734           globs->faktorLexist=TRUE;
735           break;
736         default:
737           cerr << endl;
738           cerr << "Parsing command line options produced errors \n\n";
739           exit(1);
740         }
741     }
742   //checking if x0 < x1
743   if(ex[1].fitstart >= ex[1].fitend)
744     {
745       cerr << "x0 must be smaller than x1.\n";
746       exit(1);
747     }
748   
749   //Parsing file name for a single experiment
750   if((argc-optind)!=1 && parlistSpecified==FALSE)
751     {
752       cerr << "No/Too many datafile(s) specified.  \n\n";
753       exit(1);
754     }
755   else if(parlistSpecified==FALSE)
756     {
757       ex[1].fileName=new char[strlen(argv[optind])+1];
758       strcpy(ex[1].fileName,argv[optind]);
759     }
760   
761   if(NSPLINES!=0)
762     {
763       for(k=1;k<=globs->nrExp;k++)
764         {
765           if(ex[k].splinesDefined==FALSE)
766             {
767               cerr << "Please specify spline data in experiment " << k << ".\n";
768               exit(1);
769             }
770         }
771     }
772   
773   return(ex);
774 }
\end{verbatim}\normalsize 


\index{simit_mex.cc@{simit\_\-mex.cc}!tabulateValues@{tabulateValues}}
\index{tabulateValues@{tabulateValues}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void tabulate\-Values ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex}, double {\em t0}, double {\em t1}, double $\ast$ {\em t}, long {\em n}, double $\ast$ {\em state}, double $\ast$$\ast$ {\em y}, double $\ast$$\ast$$\ast$ {\em dmds}, double $\ast$$\ast$$\ast$ {\em dmdp}, double $\ast$$\ast$ {\em dyds}, double $\ast$$\ast$ {\em dydp})}\label{simit__mex_8cc_49291cd052f6e6ad7da733e5175315f5}




Definition at line 108 of file int\-ODE.cc.

\footnotesize\begin{verbatim}110 {
111   // integrate from t0 to t1, storing values at t[1..n]
112   // t[n] may be ==t1 or <t1
113   // in the latter case we will integrate from t[n] to t1 in the
114   // final loop k=n+1 without storing any data points
115   // If dmds!=NULL, compute derivatives as well
116   // and store them in dmds,dmdp,dyds,dydp
117   // with respect to observation, see odin.doc: Observation function
118   
119   //BEGIN von Felix veraendert
120   //   double hmin=max(fabs(t0),fabs(t1))*1e-8;
121   double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
122   // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest.
123   // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern
124   // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion.
125   // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht
126   // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich
127   // von 1e-18 noch einigermassen weit entfernt bleibt. 
128   // Man koennte das in Zukunft noch in einen externen Parameter stecken.
129   //END von Felix veraendert
130 
131   long nPExp=globs->npar;
132   long nVar=ex->nvar;
133   long nobs=NOBS;
134   int generic;
135   double *p=ex->par;
136   double eps=globs->eps;  
137   double h1 = (t1 - t0) * .01;
138   long nok,nbad;
139 
140   initInt(globs,ex);
141   if (n > 0 && (t0 > t[1] || t1 < t[n]))
142     {
143       cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl;
144       throw 1;
145     }
146 
147   if (globs->integrator==3)
148     {                           // said -odessa
149 #ifdef ODESSA
150       double &rtol = eps;
151 //BEGIN von Felix veraendert
152 //       double atol=eps*1e-6;  
153       double atol = eps;
154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!)
155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu
156 // siehe odessa.f: ATOL und RTOL
157 //END von Felix veraendert
158       // odeint uses TINY=1e-30 in a similar way as odessa's atol, 
159       // but atol<1e-16 or so results in step size underflow
160       
161       int odeint_MAXSTP=globs->maxstp;
162       int d_flag=1;
163       int stif=globs->stiff;
164       char doPStr[globs->npar+1];
165       
166       for(long i=1;i<=globs->npar;i++)
167         doPStr[i]=(char)globs->doP[i];
168       
169       call_odessa(globs,ex,nVar, nPExp, doPStr,
170          d_flag, odeint_MAXSTP, stif,
171          TRUE, n, t, t0, t1, state, y,
172          dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1);
173 
174       return;
175 #else
176       cerr << "-odessa used but not linked" << endl;
177       throw 1;
178 #endif // ODESSA
179     }                           // if odessa
180       double *gy=dvector(1,nobs);
181       long i, j, k;
182       double xx1 = t0, xx2;
183       
184       if (dmds)
185         {                           // compute sensitivities
186           // create extended state vector Y
187           double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
188           initYt (Yt[0], state);
189           
190           // chain rule for extra observables
191           double **dgdy, **dgdp;
192           
193           dgdy = dmatrix (1, nobs, 1, nVar);
194           dgdp = dmatrix (1, nobs, 1, nPExp);
195           dfill0 (dgdy[1], nobs * nVar);
196           dfill0 (dgdp[1], nobs * nPExp);
197           
198           
199           for (k = 1; k <= n + 1; k++)
200             {
201               xx2 = (k > n) ? t1 :  // end point
202                 t[k];               // data point
203               
204               globs->rkqs_ign = (nPExp + nVar) * nVar;
205               if (xx2 > xx1 && nVar > 0)
206                 {                 
207                   if(globs->integrator==1)
208                     {
209                       //Runge-Kutta integation
210                       integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,
211                                   eps, h1, hmin, &nok, &nbad, sensDerivs);
212                     }
213                   else if(globs->integrator==2)
214                     {
215                       cerr << "CVODES not optimized skip integrator" << endl;
216                       throw 1;
217             
218                       //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE);
219                     }
220                 }
221               // write state to y and sensitivities to dm/dp and dm/ds
222               if (k <= n)
223                 {
224                   dcopy (y[k], Yt[0], nVar);
225                   generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp);
226                   for(i=1;i<=nobs;i++)
227                     y[k][i]=gy[i];
228                   
229                   for (i = 1; i <= nobs; i++)
230                     {
231                       for (j = 1; j <= nPExp; j++)
232                         {
233                           double &dest = dmdp[k][i][j];
234                           if (!generic)
235                             {
236                               dest = dgdp[i][j];
237                               for (long l = 1; l <= nVar; l++)
238                                 dest += dgdy[i][l] * Yt[j][l];
239                             }
240                           else
241                             dest = Yt[j][i];
242                         }           // for j
243                       for (j = 1; j <= nVar; j++)
244                         {
245                           double &dest = dmds[k][i][j];
246                           if (!generic)
247                             {
248                               dest = 0;
249                               for (long l = 1; l <= nVar; l++)
250                                 dest += dgdy[i][l] * Yt[nPExp + j][l];
251                             }
252                           else
253                             dest = Yt[nPExp + j][i];
254                         }           // for j
255                     }               // for i
256                 }                   // if k<=n
257               xx1 = xx2;
258             }                       // for k
259           
260           // copy state and its derivatives to state, dydp and dyds
261           dcopy (state, Yt[0], nVar);       // original state vector
262           // write sensitivities to dy/dp and dy/ds
263           for (i = 1; i <= nVar; i++)
264             {
265               for (j = 1; j <= nPExp; j++)
266                 dydp[i][j] = Yt[j][i];
267               for (j = 1; j <= nVar; j++)
268                 dyds[i][j] = Yt[nPExp + j][i];
269             }
270           
271           free_dmatrix (dgdp,1 ,nobs , 1, nPExp);
272           free_dmatrix (dgdy,1 ,nobs, 1, nVar);
273           free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar);
274         }
275       else
276         {                           // !dmds
277           for (k = 1; k <= n + 1; k++)
278             {
279               xx2 = (k <= n) ? t[k] : t1;
280               
281               globs->rkqs_ign = 0;
282               if (xx2 > xx1 && nVar > 0)
283                 { 
284                   if(globs->integrator==1)
285                     {
286                       //Runge-Kutta integration
287                       integrateRK(globs,state, p, nVar, xx1, xx2,
288                                   eps, h1, hmin, &nok, &nbad, ode);
289                     }
290                   else if(globs->integrator==2)
291                     {
292                       cerr << "CVODES not optimized skip integrator" << endl;
293                       throw 1;
294 
295                       //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE);
296                     }
297                 }
298               if (k <= n)
299                 {
300                   dcopy (y[k], state, nVar);
301                   generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL);
302                   for(i=1;i<=nobs;i++)
303                     y[k][i]=gy[i];
304                 }
305               xx1 = xx2;
306             }                       // for k
307         }           
308       // if dmds else
309       free_dvector(gy,1,nobs); 
310 }
\end{verbatim}\normalsize 




\subsection{Variable Documentation}
\index{simit_mex.cc@{simit\_\-mex.cc}!dbg@{dbg}}
\index{dbg@{dbg}!simit_mex.cc@{simit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}ofstream$\ast$ {\bf dbg}}\label{simit__mex_8cc_c6b594e76b642beb861b55fce626b3f1}


debug stream (C++) 



Definition at line 29 of file simit\_\-mex.cc.

Referenced by initialise(), main(), and objective\-F().