\section{/home/peifer/diffit/diffit.cc File Reference}
\label{diffit_8cc}\index{/home/peifer/diffit/diffit.cc@{/home/peifer/diffit/diffit.cc}}
Root program for {\bf diffit} (contains function main). 

{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$stdlib.h$>$}\par
{\tt \#include \char`\"{}def.h\char`\"{}}\par
{\tt \#include \char`\"{}model.h\char`\"{}}\par
{\tt \#include \char`\"{}nr.h\char`\"{}}\par
\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 read\-Data} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs, long exp\-Nr)
\item 
void {\bf set\-Mesh} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs, long exp\-Nr)
\item 
void {\bf out\-Fit} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs)
\item 
void {\bf free\-Mem} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs, int simit)
\item 
void {\bf initialise} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs, int simit)
\item 
void {\bf sim\-Init} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs)
\item 
void {\bf fit\-It} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs)
\item 
void {\bf glob\-Opt} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs)
\item 
{\bf main} (int argc, char $\ast$argv[$\,$])
\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{Detailed Description}
Root program for {\bf diffit} (contains function main). 

Diffit is has a modular structure. This file can be considered as root module for the entire C stand-alone version. A matlab interface is also available under: {\bf diffit\_\-mex.cc}{\rm (p.\,\pageref{diffit__mex_8cc})}. The structure of the matlab interface is similar than this file, apart from parsing the data and the program parameters. Note that all global variables are capsuled into two C structures named {\bf Glob}{\rm (p.\,\pageref{structGlob})} and {\bf Glob\-Exp}{\rm (p.\,\pageref{structGlobExp})}. The first structure {\bf Glob}{\rm (p.\,\pageref{structGlob})} stores all experiment unspecific parameters such as the used integrator, or the integration accuracy, etc. Experiment specific program variables are held in {\bf Glob\-Exp}{\rm (p.\,\pageref{structGlobExp})}, these are, e.g., the number of multiple shooting intervals. Modules called by function {\bf main} is in order of occurrence:

\begin{enumerate}
\item {\bf parseopts()}{\rm (p.\,\pageref{diffit_8cc_3422340b235f8185bdc6533ff0df9921})}: Parsing all program parameters either from the commant line or a file.\item {\bf read\-Data()}{\rm (p.\,\pageref{diffit_8cc_2bf00ef122e7627b0af66c27e8adf570})}: Reading the data.\item {\bf set\-Mesh()}{\rm (p.\,\pageref{diffit_8cc_343605cc3bf29ae8d91094465c787e95})}: This module sets up the multiple shooting intervals for each observation and experiment.\item {\bf initialise()}{\rm (p.\,\pageref{diffit_8cc_c66882342aff1604543a41022eb8a3f7})}: Allocates memory and initializes vectors and matrices.\item {\bf sim\-Init()}{\rm (p.\,\pageref{diffit_8cc_1c8ff8f8413349876df41981e468ea78})}: If some state variables are not directly obseved (non trivial observation function), this function integrates a trajectory using the initial guess to set up the initial guess for the state variables at multiple shooting intervals. This is activated by the command line argument {\bf \char`\"{}-siminit\char`\"{}}. If not used, the state variables are initialized with 1. Since the integrated trajectory for the initial guess is continuous, rather than discontinuous the {\bf \char`\"{}-pert $<$value$>$\char`\"{}} perturbs this guess using mean zero Gaussian noise. The standard deviation of this noise can be chosen by the value of the -pert argument.\item {\bf glob\-Opt()}{\rm (p.\,\pageref{diffit_8cc_d83d1b2a175d158bafd1cb79d3aeb172})}: Experimental !! Global optimizer instead of multiple shooting. {\bf Currently} {\bf not} {\bf working} {\bf properly}.\item {\bf fit\-It()}{\rm (p.\,\pageref{diffit_8cc_f78b0cce8098e0edba6bd140ffdef395})}: {\bf General} {\bf fitting} {\bf function}, all the numerics is done in here.\item {\bf out\-Fit()}{\rm (p.\,\pageref{diffit_8cc_89e77e765b35e43fc6de7e0d52e83ae7})}: This routine writes the output such as parameter estimates and its standard error.\item {\bf free\-Mem()}{\rm (p.\,\pageref{diffit_8cc_c39e0b1a98465f5380027ff10e613058})}: Deallocation of vectors and matrices. \end{enumerate}


Definition in file {\bf diffit.cc}.

\subsection{Function Documentation}
\index{diffit.cc@{diffit.cc}!fitIt@{fitIt}}
\index{fitIt@{fitIt}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void fit\-It ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs})}\label{diffit_8cc_f78b0cce8098e0edba6bd140ffdef395}




Definition at line 33 of file fit\-It.cc.

\footnotesize\begin{verbatim}34 {
35   long nExp,i,j,k,l,n;
36   long tryOnceMore,ok;
37   double objF, sum, scale, lambda;
38   long maxPoints=0,maxVar=0;
39   double minImprovement=globs->minimp;
40 
41   double uS=1.;
42 
43   for (i=1;i<=globs->nrExp;++i) 
44     {
45       if (ex[i].nPoints>maxPoints) 
46         maxPoints=ex[i].nPoints;
47       if (ex[i].nvar>maxVar)       
48         maxVar=ex[i].nvar;
49     }
50   
51   double **dPSave=dmatrix(1,globs->nrExp,1,globs->npar);      
52   double ***dSSave=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
53   double **P0=dmatrix(1,globs->nrExp,1,globs->npar);
54   double ***S0=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
55   double **dP=dmatrix(1,globs->nrExp,1,globs->npar);      
56   double ***dS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
57   
58   *dbg << "Minimum norm for Gauss-Newton update: minImprovement = " << minImprovement << "\n\n";
59   
60   // GGN iteration
61   // do the actual fitting
62   
63   globs->nIter=1; 
64 
65   tryOnceMore=TRUE;
66   while (globs->nIter < globs->maxit && tryOnceMore) 
67     {
68       globs->gnuindex=1;
69       if(globs->silent!=TRUE)
70         {
71           cout << "Integration # " << globs->nIter << '\n';
72           for(j=1;j<=globs->nrExp;j++)
73             {
74               cout << "Parameter value(s) of experiment " << j << ":\n";
75               for (i=1; i<=globs->npar; ++i) 
76                 {
77                   cout << ex[j].par[i] << " ";
78                 }
79               cout << endl;
80               cout << "Initial value(s) of experiment " << j << ":\n";
81               for (i=1; i<=ex[j].nvar; ++i) 
82                 {
83                   cout << ex[j].yTry[1][i] << " ";
84                 }
85               cout << endl;
86             }
87           cout << "\n";
88           cout.flush();
89         }
90       *dbg << "Integration # " << globs->nIter << '\n';
91       *dbg << "Parameter value(s):";
92       for (i=1; i<=globs->npar; ++i) 
93         {
94           if(globs->doP[i]!='L')
95             *dbg << '\t' << ex[1].par[i];
96         }
97       *dbg << "\n\n";
98       
99       // all experiments
100       sum=0;
101       for (nExp=1;nExp<=globs->nrExp;++nExp) 
102         {
103 #ifdef PRINTINITVALUES
104           *dbg << "Experiment #" << nExp << "\n";
105           *dbg << "Initial values at mesh points:\n";
106           for (i=1; i<ex[nExp].nPoints; ++i) 
107             {
108               *dbg << ex[nExp].mesh[i] << ':';
109               for (j=1; j<=ex[nExp].nvar; ++j) {
110                 *dbg << '\t' << ex[nExp].yTry[i][j];
111               }
112               *dbg << '\n';
113             }
114           dbg->flush();
115 #endif    
116           intODE(&ex[nExp],globs,TRUE,TRUE,nExp);
117           
118           // set up internal data (residues, etc.)
119           objF=computeRight(globs,&ex[nExp]);
120           sum += objF;
121           
122           *dbg << "Experiment #" << nExp << ": objective f ="
123                << objF << "\n";
124           *dbg << "Discrepancies: ";
125           for (i=1; i<ex[nExp].nPoints; ++i) 
126             {
127               for (j=1; j<=ex[nExp].nvar; ++j)
128                 *dbg << ' ' << ex[nExp].h[i][j];
129               *dbg << ",";
130             }
131           *dbg << "\n\n";
132           dbg->flush();
133         }
134 
135       if(globs->silent!=TRUE)
136         cout << "Objective function f = " << sum << "\n";
137       *dbg << "Global objective function f = " << sum << "\n";
138       
139           for (nExp=1;nExp<=globs->nrExp;++nExp) 
140             {
141               for (i=1; i<ex[nExp].nPoints; ++i)
142                 {
143                   for (j=1; j<=ex[nExp].nvar; ++j) 
144                     {
145                       ex[nExp].yTrySave[i][j] = ex[nExp].yTry[i][j];
146                     }
147                 }
148             }
149       // solve linearized minimisation problem
150       solvLin(globs,ex,FALSE);
151 
152       //save old values and new updates
153       for (nExp=1;nExp<=globs->nrExp;nExp++) 
154         {
155           for (i=1; i<ex[nExp].nPoints; i++)
156             {
157               for (j=1; j<=ex[nExp].nvar; j++)
158                 {
159                   S0[nExp][i][j] = ex[nExp].yTry[i][j];
160                   dSSave[nExp][i][j]=ex[nExp].dS[i][j];
161                   dS[nExp][i][j]=ex[nExp].dS[i][j];
162                 }
163             }
164           for (i=1; i<=globs->npar; i++)
165             {
166               P0[nExp][i] = ex[nExp].par[i];
167               dPSave[nExp][i]=ex[nExp].dP[i];
168               dP[nExp][i]=ex[nExp].dP[i];
169             }
170         }
171       
172       uS=1.;
173       //damping
174       if(!globs->nodamp)
175         lambda=dampIt(globs,ex,P0,S0,dP,dS,&uS);
176       else
177         lambda=1.;
178 
179       globs->Lambda=lambda;
180 
181       globs->chisq=sum;
182       *dbg << "Chisq = " << globs->chisq << endl;
183       
184       // next guess
185       *dbg << "using lambda=" << lambda << "\n";
186       
187       sum=0;
188       for (nExp=1;nExp<=globs->nrExp;++nExp) 
189         {
190           for (i=1; i<ex[nExp].nPoints; ++i)
191             {
192               for (j=1; j<=ex[nExp].nvar; ++j) 
193                 {
194                   sum += (lambda*uS*dSSave[nExp][i][j])*(lambda*uS*dSSave[nExp][i][j]);
195                 }
196             }
197           for (i=1; i<=globs->npar; ++i) 
198             {
199               if(globs->doP[i]=='L')
200                 sum  += (lambda*dPSave[nExp][i])*(lambda*dPSave[nExp][i]);
201             }
202         }
203       for (i=1; i<=globs->npar; ++i) 
204         {
205           if(globs->doP[i]!='L')
206             sum  += (lambda*dPSave[1][i])*(lambda*dPSave[1][i]);
207         }
208       
209       tryOnceMore = (sqrt(sum)/lambda > minImprovement);
210       
211       if(uS==0)
212         tryOnceMore=TRUE;
213 
214       if (tryOnceMore) 
215         {
216           // set up new initial values for next iteration
217           for (nExp=1;nExp<=globs->nrExp;++nExp) 
218             {
219               for (i=1; i<ex[nExp].nPoints; ++i)
220                 {
221                   for (j=1; j<=ex[nExp].nvar; ++j) 
222                     {
223                       
224                       ex[nExp].yTry[i][j] = S0[nExp][i][j]+lambda*uS*dSSave[nExp][i][j];
225                       // positiveness of initial cond --> test
226                       if(ex[nExp].yTry[i][j]<0.0)
227                         ex[nExp].yTry[i][j]=0;
228                     }
229                 }
230               for (i=1; i<=globs->npar; ++i) 
231                 {
232                   ex[nExp].par[i]=P0[nExp][i]+lambda*dPSave[nExp][i];
233                 }
234             }
235         } 
236 
237       // else:
238       // keep old values; parCovar is covariance matrix at solution point
239       // (update is too small to be worth iterating once more)
240       
241         ++globs->nIter;
242         if(globs->silent!=TRUE)
243           {
244             cout << "Norm of update vector: " << sqrt(sum) <<"\n\n";
245             cout.flush();
246           }
247         *dbg << "Norm of update vector: " << sqrt(sum) <<"\n\n";
248         
249         dbg->flush();
250         if(globs->wait)
251           {
252             cout << "<Press any key to continue>\n";
253             getchar();
254           }
255     }
256   
257   if (globs->nIter>=globs->maxit) 
258     {
259       if(globs->silent!=TRUE)
260         {
261           cout << "fitit: no convergence after " << globs->nIter << " iterations\n";
262           cout << "PRELIMINARY RESULTS! \n";
263         }
264       globs->fitConverged=FALSE;
265       
266     }
267   else
268     globs->fitConverged=TRUE;
269   //to obtain the covariance matrix
270   *dbg << "last call of solvLin to obtain the covariance matrix\n";
271   if(globs->reg=TRUE)
272     {
273       globs->minimiser=1;
274       //globs->reg=FALSE;
275       solvLin(globs,ex,TRUE);
276     }
277   else
278     {
279       globs->minimiser=1;
280       //globs->reg=FALSE;
281       solvLin(globs,ex,TRUE);
282     }
283 
284   if(!globs->wait && !globs->nowait)
285     {
286       cout << "<Press any key to continue>\n";
287       getchar();
288     }
289   
290   
291   if(globs->noGnu==FALSE)
292     {
293       for(j=1;j<=globs->ngnu;j++)
294         pclose(globs->gnuFp[j]);
295     }
296 
297   free_dmatrix(P0,1,globs->nrExp,1,globs->npar);
298   free_d3tensor(S0,1,globs->nrExp,1,maxPoints,1,maxVar);
299   free_dmatrix(dPSave,1,globs->nrExp,1,globs->npar);
300   free_d3tensor(dSSave,1,globs->nrExp,1,maxPoints,1,maxVar);       
301   free_dmatrix(dP,1,globs->nrExp,1,globs->npar);
302   free_d3tensor(dS,1,globs->nrExp,1,maxPoints,1,maxVar);
303 }
\end{verbatim}\normalsize 


\index{diffit.cc@{diffit.cc}!freeMem@{freeMem}}
\index{freeMem@{freeMem}!diffit.cc@{diffit.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{diffit_8cc_c39e0b1a98465f5380027ff10e613058}




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

\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{diffit.cc@{diffit.cc}!globOpt@{globOpt}}
\index{globOpt@{globOpt}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void glob\-Opt ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs})}\label{diffit_8cc_d83d1b2a175d158bafd1cb79d3aeb172}




Definition at line 89 of file glob\-Opt.cc.

References ex\-Ptr, globs\-Ptr, lvector(), Glob::n\-Iter, Glob\-Exp::n\-Points, and Glob::nr\-Exp.

\footnotesize\begin{verbatim}90 {
91   globsPtr=globs;
92   exPtr=ex;
93 
94 
95   globsPtr->nIter=0;
96   long i,j;
97   //variables for libSRES
98   ESParameter *param;
99   ESPopulation *population;
100   ESStatistics *stats;
101   ESfcnTrsfm *trsfm;
102   unsigned int seed;
103   int es;
104   int dim=0;
105   double *ub, *lb;
106   int miu, lambda, gen;
107   double gamma, alpha, varphi;
108   int retry;
109   double pf; 
110 
111   //keeping old values
112   long unsigned *nms=lvector(1,globsPtr->nrExp);
113   for(i=1;i<=globsPtr->nrExp;i++)
114     {
115       nms[i]=exPtr[i].nPoints;
116       exPtr[i].nPoints=2;
117     }
118 
119   seed = shareDefSeed;
120   gamma = esDefGamma;
121   alpha = esDefAlpha;
122   varphi = esDefVarphi;
123   retry = esDefRetry;
124   pf = essrDefPf;
125   es = esDefESSlash;
126   miu = 30;
127   lambda = 200;
128   gen = 1750;
129   trsfm = NULL;
130 
131   //determine # of decision variables
132   for(i=1;i<=globsPtr->nrExp;i++)
133     {
134       for(j=1;j<=NEQNS;j++)
135         {
136           if(globsPtr->y0fix[j]==TRUE)
137             dim++;
138         }
139       for(j=1;j<=globsPtr->npar;j++)
140         {
141           if(globsPtr->doP[j]=='L')
142             dim++;
143         }
144     }      
145   for(j=1;j<=globsPtr->npar;j++)
146     {
147       if(globsPtr->doP[j]==TRUE)
148         dim++;
149     }
150 
151 
152   //bounds ...for test only
153 
154   lb=dvector(0,dim);
155   ub=dvector(0,dim);
156 
157   for(i=0;i<dim;i++)
158     {
159       lb[i]=0;
160       ub[i]=2;
161     }
162 
163   ESInitial(seed, &param, trsfm, fitness, es,0, dim, ub, lb, \
164             miu, lambda, gen, gamma, alpha, varphi,  \
165             retry, &population, &stats);
166   
167   for(int i=1;i<=globsPtr->maxit;i++)
168     { 
169       globsPtr->nIter++;
170       ESStep(population, param, stats, pf);
171       if(stats->curgen >= param->gen)
172         break;
173     }
174   
175   //copying old stuff back
176   for(i=1;i<=globsPtr->nrExp;i++)
177     {
178       exPtr[i].nPoints=nms[i];
179     }
180 
181   cout << endl;
182 
183   free_dvector(lb,0,dim);
184   free_dvector(ub,0,dim);
185   free_lvector(nms,1,globsPtr->nrExp);
186 }
\end{verbatim}\normalsize 


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




Definition at line 98 of file initialise.cc.

Referenced by mex\-Function().

\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{diffit.cc@{diffit.cc}!main@{main}}
\index{main@{main}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}main (int {\em argc}, char $\ast$ {\em argv}[$\,$])}\label{diffit_8cc_e0665038b72011f5c680c660fcb59459}


function main 

Definition at line 81 of file diffit.cc.

References dbg, Glob::nr\-Exp, and parseopts().

\footnotesize\begin{verbatim}82 {
83   long k,i,j;
84   char outstr[100];
85   Glob globs;
86   GlobExp *ex;
87   
88   //  parseopt
89   ex=parseopts(argc,argv,&globs,outstr);
90 
91   //open DEBUG stream
92   dbg=new ofstream("diffit.dbg");
93   if (!dbg) 
94     {
95       cerr << "Error opening DEBUG file.\n";
96       exit(1);
97     }
98   dbg->precision(4);     // set output format for floating point data
99   dbg->unsetf(ios::floatfield);
100 
101   //print debugging information
102   *dbg << DefModelDescription << "\n";
103   if (globs.nrExp==1) 
104     *dbg << "1 experiment\n";
105   else
106     *dbg << globs.nrExp << " experiments\n";
107    *dbg << "\n";
108    for(i=1;i<=globs.nrExp;i++)
109      {
110        *dbg << "Experiment: " << i << "\n";
111        *dbg << NPARAMS << " parameter(s):";
112        for (k=1; k<=NPARAMS; ++k) 
113          *dbg << " " << ex[i].par[k];
114        *dbg << "\n";
115      }
116    
117   dbg->flush();
118  
119   for(i=1;i<=globs.nrExp;i++)
120     {
121       // readData()
122       readData(ex,&globs,i);
123       //set mesh
124       setMesh(ex,&globs,i);
125 
126 #ifdef PRINTDATA
127       *dbg << "\nData: \n";
128       for (j=1; j<=ex[i].nMeasure; j++) {
129         *dbg << ex[i].xMeasure[j];
130         for (k=1; k<=ex[i].nobs; k++)
131           *dbg << "\t" << ex[i].yMeasure[j][k] << "\t" << ex[i].sigma[j][k];
132         *dbg << "\n";
133       }
134       *dbg << "\n";
135 #endif
136       *dbg << "Mesh:";
137       for (k=1;k<=ex[i].nPoints;++k) 
138         *dbg << " " << ex[i].mesh[k];
139       *dbg << "\n";
140     }
141 
142   //initial values
143   for(i=1;i<=globs.nrExp;i++)
144     {
145       *dbg << "\nExperiment " << i << ":\n";
146       if (ex[i].y0) 
147         {
148           *dbg << NEQNS << " starting value(s):";
149           for (k=1; k<=NEQNS; ++k) 
150             *dbg << ' ' << ex[i].y0[k];
151           *dbg << "\n";
152         } 
153       else 
154         {
155           *dbg << "No starting values specified\n";
156         }
157       *dbg << "\n";
158       dbg->flush();
159     }
160   
161    initialise(ex,&globs,FALSE);
162    
163    //simulate initial state
164    if(globs.simInit==TRUE)
165      {
166        for (i=1;i<=globs.nrExp;++i) 
167          simInit(&ex[i],&globs);
168      }
169    
170    // starting of the numerics
171    //*************************
172       
173    if(globs.strategy==2)
174      globOpt(ex,&globs);
175    else
176      {
177        try 
178          {            
179             fitIt(ex,&globs);
180          }
181        catch(int i)  //Exeption handling
182          {
183            if(globs.noGnu==FALSE)
184              {
185                for(j=1;j<=globs.ngnu;j++)
186                  pclose(globs.gnuFp[j]);
187              }       
188            exit(1);
189          }
190      }
191 
192    
193    //*************************
194 
195   //Output after convergence
196   outFit(ex,&globs);
197   if(!globs.noGnu)
198     system("rm -f gnuout.dat");
199 
200   freeMem(ex,&globs,FALSE);
201 }
\end{verbatim}\normalsize 


\index{diffit.cc@{diffit.cc}!outFit@{outFit}}
\index{outFit@{outFit}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void out\-Fit ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs})}\label{diffit_8cc_89e77e765b35e43fc6de7e0d52e83ae7}




Definition at line 11 of file out\-Fit.cc.

\footnotesize\begin{verbatim}12 {
13   long i,j;
14   long nExp;
15   long nvar,npar,nglob=0;
16   long ind=1;
17 
18   npar=globs->npar;     
19   nvar=ex[1].nvar;
20 
21   double **errorS=dmatrix(1,globs->nrExp,1,nvar);
22   double **errorP=dmatrix(1,globs->nrExp,1,npar);
23   
24 
25   if(globs->strategy!=2)
26     {
27       //prepare output
28       for(nExp=1;nExp<=globs->nrExp;++nExp) 
29         {
30           for(i=1;i<=nvar;i++)
31             errorS[nExp][i]=-1;
32           for(i=1;i<=npar;i++)
33             errorP[nExp][i]=-1;
34           
35           for(i=1;i<=nvar;i++)
36             {
37               if(globs->y0fix[i]!=FALSE)
38                 {
39                   errorS[nExp][i]=sqrt(globs->covar[ind][ind]);
40                   ind++;
41                 }
42             }
43           for(i=1;i<=npar;i++)
44             {
45               if(globs->doP[i]=='L')
46                 {
47                   errorP[nExp][i]=sqrt(globs->covar[ind][ind]);
48                   ind++;
49                   nglob++;
50                 }
51             }
52         }
53       
54       for(i=1;i<=npar;i++)
55         {
56           if(globs->doP[i]==TRUE)
57             {
58               for(nExp=1;nExp<=globs->nrExp;++nExp) 
59                 {
60                   errorP[nExp][i]=sqrt(globs->covar[ind][ind]);
61                 }
62               ind++;
63             }
64         } 
65     }  //end if(globs->strategy..)
66 
67   /* print results */
68   if(globs->silent!=TRUE)
69     {
70       cout << "Number of iterations: " << globs->nIter << ", chi^2 = " << globs->chisq << "\n\n";
71       cout << "\nBest fit parameters +/- standard errors:\n";
72       cout << "----------------------------------------\n\n"; 
73       
74       cout << "Global Parameters:\n";
75       for(j=1;j<=npar;j++)
76         { 
77           if(globs->doP[j]!='L')
78             {
79               if(globs->strategy!=2)
80                 {
81                   for(nExp=1;nExp<=globs->nrExp;++nExp) 
82                     ex[nExp].errP[j]=errorP[1][j];
83                   cout << ParameterNames[j-1];
84                   if(errorP[1][j]!=-1)
85                     cout << " = " << ex[1].par[j] <<  " +/- " << errorP[1][j] << endl;
86                   else
87                     cout << " (fixed)\n";
88                 }
89               else
90                 cout << ParameterNames[j-1] << " = " << ex[1].par[j] << endl;
91                 
92                 
93             }
94         }
95       cout << endl;
96       for(i=1;i<=globs->nrExp;++i) 
97         {
98           cout << "Experiment " << i << ":";
99           if(nglob!=0)
100             cout << "\n\nLocal Parameters:\n";
101           else
102             cout << endl;
103           
104           for(j=1;j<=npar;j++)
105             {
106               if(globs->doP[j]=='L')
107                 {
108                   if(globs->strategy!=2)
109                     {
110                       ex[i].errP[j]=errorP[i][j];
111                       cout << ParameterNames[j-1] << " = " << ex[i].par[j];
112                       if(errorP[i][j]!=-1)
113                         cout <<  " +/- " << errorP[i][j] << endl;
114                       else
115                         cout << " (fixed)\n";
116                     }
117                   else
118                     cout << ParameterNames[j-1] << " = " << ex[i].par[j] << endl;
119                 }
120             }
121           cout << "Initial Values" << endl;
122           for(j=1;j<=nvar;j++)
123             {
124               if(globs->strategy!=2)
125                 {
126                   ex[i].errY0[j]=errorS[i][j];
127                   cout << VariableNames[j-1] << " = " << ex[i].yTry[1][j];
128                   if(errorS[i][j]!=-1)
129                     cout <<  " +/- " << errorS[i][j] << endl;
130                   else
131                     cout << " (fixed)\n";
132                 }
133               else
134                 cout << VariableNames[j-1] << " = " << ex[i].yTry[1][j] << endl;
135             }
136           cout << endl;
137         }
138     } //end of --> if(silent!=TRUE)
139 
140   free_dmatrix(errorS,1,globs->nrExp,1,nvar);
141   free_dmatrix(errorP,1,globs->nrExp,1,npar);
142 }
\end{verbatim}\normalsize 


\index{diffit.cc@{diffit.cc}!parseopts@{parseopts}}
\index{parseopts@{parseopts}!diffit.cc@{diffit.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{diffit_8cc_3422340b235f8185bdc6533ff0df9921}


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.

Referenced by main().

\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{diffit.cc@{diffit.cc}!readData@{readData}}
\index{readData@{readData}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void read\-Data ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs}, long {\em exp\-Nr})}\label{diffit_8cc_2bf00ef122e7627b0af66c27e8adf570}




Definition at line 17 of file read\-Data.cc.

References Glob\-Exp::file\-Name, and Glob\-Exp::nobs.

\footnotesize\begin{verbatim}18 {
19   int lineNr=0;
20   long k,l,index;
21   long nLines=0,maxLeng=0;
22   long nobsmax=ex[expNr].nobs;
23   char c,*line,*arg;
24   double *x;
25   double **y;
26   double **sig;
27 
28   ifstream in;
29 
30   FILE *fp;
31 
32   fp=fopen(ex[expNr].fileName,"r");
33   if (fp==NULL) {
34     cerr << "Error opening data file " << ex[expNr].fileName << "\n";
35     exit(1);
36   }
37   
38   //read total number of lines
39   while((c=fgetc(fp))!=EOF)
40         {
41           for(k=1;c!='\n';k++)
42             {
43               if(k>maxLeng)
44                 maxLeng=k;
45               c=fgetc(fp);
46             }
47           nLines++;
48         }
49 #ifdef DEBUGREADDATA
50   *dbg << "Total number of lines in " <<  ex[expNr].fileName << " : " << nLines << endl;
51   *dbg << "Maximal length of line in " << ex[expNr].fileName << " : " << maxLeng << endl;
52 #endif
53   fclose(fp);
54 
55   in.open(ex[expNr].fileName);
56 
57   //allocate memory
58   line=new char[maxLeng+2];
59   arg=new char[maxLeng+2];
60   x=dvector(1,nLines);
61   y=dmatrix(1,nLines,1,ex[expNr].nobs);
62   sig=dmatrix(1,nLines,1,ex[expNr].nobs);
63 
64   //extracting the data
65   while(!in.eof())
66         {
67           in.getline(line,maxLeng+2,'\n');
68           if(isspace(line[0])==FALSE && line[0]!= '#' && line[0] != '\0' )
69             {
70               lineNr++;
71 #ifdef DEBUGREADDATA
72               *dbg << "Extracting data line no. : " << lineNr << endl;
73 #endif
74               l=0;
75               index=0;
76               for(k=0;k<=strlen(line);k++)
77                 {
78                   if(isspace(line[k])!=FALSE)
79                     {
80                       arg[l]='\0';
81                       l=0;
82                       if(index==0)
83                         x[lineNr]=atof(arg);
84                       else if(index <=2*ex[expNr].nobs && (index % 2)==1)
85                         y[lineNr][(index+1)/2]=atof(arg);
86                       else if(index <=2*ex[expNr].nobs && (index % 2)==0)
87                         sig[lineNr][index/2]=atof(arg);
88                       index++;
89                     }
90                   else
91                     {
92                       arg[l]=line[k];
93                       l++;
94                     }
95                 }
96               arg[l]='\0';      
97               if(index==0)
98                 x[lineNr]=atof(arg);
99               else if(index <=2*ex[expNr].nobs && (index % 2)==1)
100                 y[lineNr][(index+1)/2]=atof(arg);
101               else if(index <=2*ex[expNr].nobs && (index % 2)==0)
102                 sig[lineNr][index/2]=atof(arg);
103               
104               //determ. max. complete obs.
105               if((index % 2)==1)
106                 index--;
107               if(index==0)
108                 {
109                   cerr << "Data line " << lineNr << " contains no data." << endl;
110                   exit(1);
111                 }
112               if(index/2 <=  nobsmax)
113                 nobsmax=index/2;
114             }
115                     
116         }
117   in.close();
118 
119   ex[expNr].nobs=nobsmax;
120   ex[expNr].nMeasure=lineNr;
121   if(ex[expNr].nMeasure==0)
122     {
123       cerr << "File " << ex[expNr].fileName << " contains no data." << endl;
124       exit(1);
125     }
126 
127   //check if the time is in ascending order
128     for(k=2;k<=ex[expNr].nMeasure;k++)
129       {
130         if(x[k-1] >= x[k])
131           {
132             cerr << "Data not in temporal ascending order at dataline " << k << "\n";
133             exit(1);
134           }
135       }
136     //fit range stuff
137     if(ex[expNr].fitstart == -1e99)
138       ex[expNr].fitstart=x[1];
139     if(ex[expNr].fitend ==  1e99)
140       ex[expNr].fitend=x[ex[expNr].nMeasure];
141     
142     //for memory allocation
143     index=1;
144     for(k=1;k<=ex[expNr].nMeasure;k++)
145       {
146         if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
147           { 
148             index++;
149           }
150       }
151     index--; 
152   
153     //the actual amount of data
154     ex[expNr].nMeasure=index;
155 
156     //memory allocation
157     ex[expNr].xMeasure=dvector(1,ex[expNr].nMeasure);
158     ex[expNr].yMeasure=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
159     ex[expNr].sigma=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
160 
161     // copy data to data structure
162     index=1;
163     for(k=1;k<=lineNr;k++)
164       {
165         if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
166           {
167             ex[expNr].xMeasure[index]=x[k];
168             for(l=1;l<=ex[expNr].nobs;l++)
169               {
170                 ex[expNr].yMeasure[index][l]=y[k][l];
171                 ex[expNr].sigma[index][l]=sig[k][l];
172               }
173             index++;
174           }
175       }
176     index--;
177 
178     //setting some variables
179     ex[expNr].firstMeasure=1;
180     ex[expNr].lastMeasure=ex[expNr].nMeasure;
181     
182     //free memory
183     free_dvector(x,1,nLines);
184     free_dmatrix(y,1,nLines,1,ex[expNr].nobs);
185     free_dmatrix(sig,1,nLines,1,ex[expNr].nobs);
186 }
\end{verbatim}\normalsize 


\index{diffit.cc@{diffit.cc}!setMesh@{setMesh}}
\index{setMesh@{setMesh}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void set\-Mesh ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs}, long {\em exp\-Nr})}\label{diffit_8cc_343605cc3bf29ae8d91094465c787e95}




Definition at line 11 of file set\-Mesh.cc.

\footnotesize\begin{verbatim}12 {
13   long i,nPoints=ex[expNr].nPoints;
14   
15   ex[expNr].mesh=dvector(1,nPoints);
16   ex[expNr].mesh[1]= ex[expNr].fitstart;
17   ex[expNr].mesh[nPoints]= ex[expNr].fitend;
18 
19   if (globs->noMeasurements==FALSE) 
20     {         
21       // mesh at data points
22       if (ex[expNr].nMeasure < nPoints) 
23         {
24           cerr << "Too few measurements to construct mesh\n";
25           exit(1);
26         }
27       for (i=2; i < nPoints; ++i)
28         ex[expNr].mesh[i]=ex[expNr].xMeasure[ex[expNr].firstMeasure+int(double((i-1)*(ex[expNr].nMeasure-1))/double(nPoints-1))];
29     } 
30   else 
31     {        
32       // evenly spaced mesh
33       for (i=2; i < nPoints; ++i)
34         ex[expNr].mesh[i]=ex[expNr].fitstart+(ex[expNr].fitend-ex[expNr].fitstart)*double(i-1)/double(nPoints-1);
35     }
36 }
\end{verbatim}\normalsize 


\index{diffit.cc@{diffit.cc}!simInit@{simInit}}
\index{simInit@{simInit}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void sim\-Init ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs})}\label{diffit_8cc_1c8ff8f8413349876df41981e468ea78}




Definition at line 21 of file sim\-Init.cc.

\footnotesize\begin{verbatim}22 {
23   long i,j;
24   long nvar=ex->nvar;
25   long nPoints=ex->nPoints;
26   long idum=-42;
27   double *t=dvector(1,2);
28   double **y=dmatrix(1,2,1,nvar);
29   
30    if (ex->y0) 
31     {
32       for (i=1; i<=nvar; ++i) 
33         ex->yTry[1][i]=ex->y0[i];
34     } 
35    else
36      {
37         for (i=1; i<=nvar; ++i) 
38           {
39             ex->yTry[1][i]=DefYValues[i-1];
40           }
41      }
42 
43    //integrating 
44 
45    for(i=2;i<=nPoints;i++)
46      {
47        for(j=1;j<=nvar;j++)
48          ex->yTry[i][j]=ex->yTry[i-1][j];
49        t[1]=ex->mesh[i-1];
50        t[2]=ex->mesh[i];
51        tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
52       for(j=1;j<=nvar;j++)
53         ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
54      }
55    
56 
57 
58    free_dvector(t,1,2);
59    free_dmatrix(y,1,2,1,nvar);
60 }
\end{verbatim}\normalsize 




\subsection{Variable Documentation}
\index{diffit.cc@{diffit.cc}!dbg@{dbg}}
\index{dbg@{dbg}!diffit.cc@{diffit.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}ofstream$\ast$ {\bf dbg}}\label{diffit_8cc_c6b594e76b642beb861b55fce626b3f1}


debug stream (C++) 



Definition at line 78 of file diffit.cc.