\section{/home/peifer/diffit/diffit\_\-mex.cc File Reference}
\label{diffit__mex_8cc}\index{/home/peifer/diffit/diffit_mex.cc@{/home/peifer/diffit/diffit\_\-mex.cc}}
{\tt \#include \char`\"{}mex.h\char`\"{}}\par
{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$stdlib.h$>$}\par
{\tt \#include $<$string$>$}\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 INEX}~prhs[1]
\item 
\#define {\bf OUT}~plhs[0]
\item 
\#define {\bf OUTEX}~plhs[1]
\end{CompactItemize}
\subsection*{Functions}
\begin{CompactItemize}
\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 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{diffit_mex.cc@{diffit\_\-mex.cc}!IN@{IN}}
\index{IN@{IN}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define IN~prhs[0]}\label{diffit__mex_8cc_c2bbd6d630a06a980d9a92ddb9a49928}




Definition at line 34 of file diffit\_\-mex.cc.

Referenced by mex\-Function().\index{diffit_mex.cc@{diffit\_\-mex.cc}!INEX@{INEX}}
\index{INEX@{INEX}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define INEX~prhs[1]}\label{diffit__mex_8cc_fbe7ff759541b50d26b4dd0738caf2a0}




Definition at line 35 of file diffit\_\-mex.cc.\index{diffit_mex.cc@{diffit\_\-mex.cc}!OUT@{OUT}}
\index{OUT@{OUT}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define OUT~plhs[0]}\label{diffit__mex_8cc_ec78e7a9e90a406a56f859ee456e8eae}




Definition at line 37 of file diffit\_\-mex.cc.\index{diffit_mex.cc@{diffit\_\-mex.cc}!OUTEX@{OUTEX}}
\index{OUTEX@{OUTEX}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define OUTEX~plhs[1]}\label{diffit__mex_8cc_55dabb70c26bc03a3a88ff33b013f876}




Definition at line 38 of file diffit\_\-mex.cc.

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




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

References Glob::minimp, and Glob::nr\-Exp.

\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_mex.cc@{diffit\_\-mex.cc}!freeMem@{freeMem}}
\index{freeMem@{freeMem}!diffit_mex.cc@{diffit\_\-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{diffit__mex_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_mex.cc@{diffit\_\-mex.cc}!initialise@{initialise}}
\index{initialise@{initialise}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void initialise ({\bf Glob\-Exp} {\em ex}[$\,$], {\bf Glob} $\ast$ {\em globs}, int {\em simit})}\label{diffit__mex_8cc_c66882342aff1604543a41022eb8a3f7}




Definition at line 98 of file initialise.cc.

\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_mex.cc@{diffit\_\-mex.cc}!mexFunction@{mexFunction}}
\index{mexFunction@{mexFunction}!diffit_mex.cc@{diffit\_\-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{diffit__mex_8cc_6a215cbfde54f82a3ce599228fc3fce5}




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

References Glob::do\-P, Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob::gnu\-Fp, IN, Glob::init\-Spline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, Glob::maxstp, Glob::minimiser, Glob::minimp, Glob::nodamp, Glob::no\-Gnu, Glob::no\-Measurements, Glob::nowait, Glob::npar, Glob::nr\-Exp, Glob::pert, Glob::reg, Glob::silent, Glob::sim\-Init, Glob::stiff, TRUE, Glob::wait, and Glob::wquer.

\footnotesize\begin{verbatim}45 {
46   long k,i,j,nExp;
47   long index,nglob=0;
48   unsigned long *ndat;
49   Glob globs;
50   GlobExp *ex;
51   mxArray *tmp,*sig;
52   double *tmp_,*sig_,*spline_;
53   char tmpstr[1000];
54   mxArray *splinedat;
55   int except=0;
56   int init=0;
57 
58   //for output
59   const char *fieldNamesOUT[]={"wquer","converged","chisq","except","Lambda","cond"};
60   const char *fieldNamesOUTEX[]={"p","y0","mesh","errP","errY0"};
61 
62   int dims[2];
63 
64   //GET DATA FOR GLOBS
65 
66   globs.npar=NPARAMS;
67   globs.wait=FALSE;
68   globs.nowait=TRUE;
69   globs.elastic=1.;
70   globs.initSpline=FALSE;
71   globs.maxstp=5000;
72   globs.gnuFp=NULL;
73   globs.initSpline=TRUE; 
74   globs.wquer=100.;
75   globs.silent=FALSE;
76   globs.minimiser=2;
77 
78   //eps
79   tmp=mxGetField(IN,0,"eps");
80   tmp_=mxGetPr(tmp);
81   globs.eps=fabs(tmp_[0]);
82   //eps
83   tmp=mxGetField(IN,0,"nognu");
84   tmp_=mxGetPr(tmp);
85   globs.noGnu=(int)tmp_[0];  
86   //nomeasure
87   tmp=mxGetField(IN,0,"nomeasure");
88   tmp_=mxGetPr(tmp);
89   globs.noMeasurements=(int)tmp_[0];   
90   //maxit
91   tmp=mxGetField(IN,0,"maxit");
92   tmp_=mxGetPr(tmp);
93   globs.maxit=abs((int)tmp_[0]);  
94   //odessa
95   tmp=mxGetField(IN,0,"int");
96   tmp_=mxGetPr(tmp);
97   globs.integrator=abs((int)tmp_[0]);  
98   globs.stiff=TRUE;
99   //minimp
100   tmp=mxGetField(IN,0,"minimp");
101   tmp_=mxGetPr(tmp);
102   globs.minimp=fabs(tmp_[0]);
103   //reg
104   tmp=mxGetField(IN,0,"reg");
105   tmp_=mxGetPr(tmp);
106   globs.reg=(int)tmp_[0];
107   //epsilon
108   tmp=mxGetField(IN,0,"epsilon");
109   tmp_=mxGetPr(tmp);
110   globs.epsilon=fabs(tmp_[0]); 
111   //lambda
112   tmp=mxGetField(IN,0,"lambda");
113   tmp_=mxGetPr(tmp);
114   globs.lambda=fabs(tmp_[0]); 
115   //siminit
116   tmp=mxGetField(IN,0,"siminit");
117   tmp_=mxGetPr(tmp);
118   globs.simInit=(int)tmp_[0];
119   //pert
120   tmp=mxGetField(IN,0,"pert");
121   tmp_=mxGetPr(tmp);
122   globs.pert=fabs(tmp_[0]);
123   //nodamp
124   tmp=mxGetField(IN,0,"nodamp");
125   tmp_=mxGetPr(tmp);
126   globs.nodamp=(int)tmp_[0];
127   //nrExp
128   tmp=mxGetField(IN,0,"nexp");
129   tmp_=mxGetPr(tmp);
130   globs.nrExp=(long)tmp_[0];
131   //doP
132   globs.doP=ivector(1,globs.npar);
133   for(k=1;k<=globs.npar;k++)
134     globs.doP[k]=TRUE;
135   tmp=mxGetField(IN,0,"doP");
136   if(mxGetString(tmp,tmpstr,1000))
137     {
138       cerr << "Parsing doP failed.\n";
139       return;
140     }
141   if(tmpstr[0]!='N')
142     {
143       if(strlen(tmpstr)!=globs.npar)
144         {
145           cerr << "doP must be of length " << globs.npar << ".\n";
146           return;
147         }
148       else
149         {
150           for(k=0;k < globs.npar;k++)
151             {
152               if(tmpstr[k]=='0')
153                 globs.doP[k+1]=FALSE;
154               else if (tmpstr[k]=='L')
155                 globs.doP[k+1]='L'; 
156               else
157                 globs.doP[k+1]=TRUE;
158             }
159         }
160     }
161   //y0fix
162   globs.y0fix=ivector(1,NEQNS);
163   for(k=1;k<=NEQNS;k++)
164     globs.y0fix[k]=TRUE;
165   tmp=mxGetField(IN,0,"y0fix");
166   if(mxGetString(tmp,tmpstr,1000))
167     {
168       cerr << "Parsing y0fix failed.\n";
169       return;
170     }
171   if(tmpstr[0]!='N')
172     {
173       if(strlen(tmpstr)!=NEQNS)
174         {
175           cerr << "y0fix must be of length " << NEQNS << ".\n";
176           return;
177         }
178       else
179         {
180           for(k=0;k < NEQNS ;k++)
181             {
182               if(tmpstr[k]=='0')
183                 globs.y0fix[k+1]=FALSE;
184               else
185                 globs.y0fix[k+1]=TRUE;
186             }
187         }
188     }
189   //wquer
190   tmp=mxGetField(IN,0,"wquer");
191   tmp_=mxGetPr(tmp);
192   globs.wquer=tmp_[0];
193   //init
194   tmp=mxGetField(IN,0,"init");
195   tmp_=mxGetPr(tmp);
196   init=(int)tmp_[0];
197    //silent
198   tmp=mxGetField(IN,0,"silent");
199   tmp_=mxGetPr(tmp);
200   globs.silent=(int)tmp_[0]; 
201    //opt
202   tmp=mxGetField(IN,0,"opt");
203   tmp_=mxGetPr(tmp);
204   globs.minimiser=(int)tmp_[0]; 
205   //Lconst
206   tmp=mxGetField(IN,0,"Lconst");      
207   tmp_=mxGetPr(tmp);
208   if(mxIsEmpty(tmp))
209   {
210      globs.faktorLexist=FALSE;
211   }
212   else
213   {
214      globs.faktorL=dvector(1,NPARAMS);
215      for(k=1;k<=NPARAMS;k++)
216         globs.faktorL[k]=-1;
217      for(k=1;k<=mxGetNumberOfElements(tmp) && k<=NPARAMS;k++)
218      {
219         globs.faktorL[k]=tmp_[k-1];
220      }
221      globs.faktorLexist=TRUE;
222   }     
223 
224 
225 
226   // GET DATA FOR EX
227   ndat=lvector(1,globs.nrExp);
228   ex=new GlobExp[globs.nrExp+1];
229 
230   for(nExp=1;nExp<=globs.nrExp;nExp++)
231     {
232       //general stuff
233       ex[nExp].nvar=NEQNS;
234       ex[nExp].splinesDefined=FALSE;
235       //t0
236       tmp=mxGetField(INEX,nExp-1,"t0");
237       tmp_=mxGetPr(tmp);
238       ex[nExp].fitstart=tmp_[0];
239     
240       //t1
241       tmp=mxGetField(INEX,nExp-1,"t1");
242       tmp_=mxGetPr(tmp);
243       ex[nExp].fitend=tmp_[0];
244       if(ex[nExp].fitstart >= ex[nExp].fitend)
245         {
246           cerr << "Fitstart >= fitend in exp. " << nExp << ".\n";
247           return;
248         }
249       //nobs
250       tmp=mxGetField(INEX,nExp-1,"nobs");      
251       tmp_=mxGetPr(tmp);
252       ex[nExp].nobs=(long)tmp_[0];
253       if(ex[nExp].nobs > NOBS)
254         {
255           cerr << "More observables than required in exp. " << nExp << ".\n";
256           return;
257         }
258       //nms
259       tmp=mxGetField(INEX,nExp-1,"nms");      
260       tmp_=mxGetPr(tmp);
261       ex[nExp].nPoints=(long)tmp_[0]+1;
262       //p
263       ex[nExp].par=dvector(1,NPARAMS);
264       tmp=mxGetField(INEX,nExp-1,"p");  
265       tmp_=mxGetPr(tmp);
266       if(mxIsEmpty(tmp))
267         {
268           for(k=1;k<=globs.npar;k++)
269             ex[nExp].par[k]=DefParameters[k-1];
270         }
271       else
272         {
273           if(mxGetNumberOfElements(tmp)!=globs.npar)
274             {
275               cerr << globs.npar << " parameter(s) required in exp. " << nExp <<".\n";
276               return;
277             }
278           else
279             {
280               for(k=1;k<=globs.npar;k++)
281                 ex[nExp].par[k]= tmp_[k-1];
282             }
283         }
284       //y0
285       tmp=mxGetField(INEX,nExp-1,"y0");      
286       tmp_=mxGetPr(tmp);
287       if(mxIsEmpty(tmp))
288         {
289           ex[nExp].y0=NULL;
290         }
291       else
292         {
293           if(mxGetNumberOfElements(tmp)!=ex[nExp].nvar)
294             {
295               cerr << ex[nExp].nvar << " initial value(s) required in exp. " << nExp <<".\n";
296               return;
297             }
298           else
299             {
300               ex[nExp].y0=dvector(1,ex[nExp].nvar);
301               for(k=1;k<=ex[nExp].nvar;k++)
302                 ex[nExp].y0[k]= tmp_[k-1];
303             }
304         }     
305       //n
306       tmp=mxGetField(INEX,nExp-1,"n");      
307       tmp_=mxGetPr(tmp);
308       ndat[nExp]=(long)tmp_[0];
309       
310     } //end for(nExp=1;nExp<=globs.nrExp;nExp++)
311 
312   //open DEBUG stream
313   dbg=new ofstream("diffit.dbg");
314   if (!dbg) 
315     {
316       cerr << "Error opening DEBUG file.\n";
317       return;
318     }
319   dbg->precision(4);     // set output format for floating point data
320   dbg->unsetf(ios::floatfield);
321   
322   //print debugging information
323   *dbg << DefModelDescription << "\n";
324   if (globs.nrExp==1) 
325     *dbg << "1 experiment\n";
326   else
327     *dbg << globs.nrExp << " experiments\n";
328   *dbg << "\n";
329   for(i=1;i<=globs.nrExp;i++)
330     {
331       *dbg << "Experiment: " << i << "\n";
332       *dbg << NPARAMS << " parameter(s):";
333       for (k=1; k<=NPARAMS; ++k) 
334         *dbg << " " << ex[i].par[k];
335       *dbg << "\n";
336     }
337    
338   dbg->flush();
339   
340  
341   for(i=1;i<=globs.nrExp;i++)
342     {
343       //read data and check for for consistency
344 
345       //data 
346       tmp=mxGetField(INEX,i-1,"data");      
347       tmp_=mxGetPr(tmp);
348       sig=mxGetField(INEX,i-1,"sig");      
349       sig_=mxGetPr(sig);
350       //deternine nMeasure
351       index=1;
352       for(j=1;j<=ndat[i];j++)
353         {
354           if(tmp_[j-1]>= ex[i].fitstart && tmp_[j-1]<= ex[i].fitend)
355             index++;
356         }
357       index--; 
358       if(index==0)
359         {
360           cerr << "Nothing to fit in exp " << i << ".\n";
361           return;
362         }
363       ex[i].nMeasure=index;
364       //memory allocation
365       ex[i].xMeasure=dvector(1,ex[i].nMeasure);
366       ex[i].yMeasure=dmatrix(1,ex[i].nMeasure,1,ex[i].nobs);
367       ex[i].sigma=dmatrix(1,ex[i].nMeasure,1,ex[i].nobs);
368       //copy data
369       index=1;
370       for(j=1;j<=ndat[i];j++)
371         {
372           if(tmp_[j-1]>= ex[i].fitstart && tmp_[j-1]<= ex[i].fitend)
373             {
374               ex[i].xMeasure[index]=tmp_[j-1];
375               for(k=1;k<=ex[i].nobs;k++)
376                 {
377                   ex[i].yMeasure[index][k]=tmp_[ndat[i]*k+j-1];
378                   ex[i].sigma[index][k]=sig_[ndat[i]*(k-1)+j-1];
379                 }
380               index++;
381             }
382         }
383       //check if time is in ascending order
384       for(j=2;j<=ex[i].nMeasure;j++)
385         {
386           if(ex[i].xMeasure[j-1] >= ex[i].xMeasure[j])
387             {
388               cerr << "Data not in temporal ascending order in exp " << i << ",\n";
389               return;
390             }
391         }
392       //setting some variables
393       ex[i].firstMeasure=1;
394       ex[i].lastMeasure=ex[i].nMeasure;
395       
396       //set mesh
397       tmp=mxGetField(INEX,i-1,"mesh");      
398       tmp_=mxGetPr(tmp);
399       if(mxIsEmpty(tmp))
400         setMesh(ex,&globs,i);
401       else
402         {
403           ex[i].nPoints=mxGetM(tmp)+2;
404           if (ex[i].nMeasure < ex[i].nPoints) 
405             {
406               cerr << "Too few measurements to construct mesh\n";
407               return;
408             }
409 
410           if(mxGetN(tmp) <= ex[i].nvar+1)
411             {
412               ex[i].mesh=dvector(1,ex[i].nPoints);
413               ex[i].mesh[1]= ex[i].fitstart;
414               ex[i].mesh[ex[i].nPoints]= ex[i].fitend;
415               //reading mesh data
416               for(j=1;j<=ex[i].nPoints-2;j++)
417                 {
418                   ex[i].mesh[j+1]=tmp_[j-1];
419                 }
420               //check if mesh is in ascending order
421               for(j=2;j<=ex[i].nPoints;j++)
422                 {
423                   if(ex[i].mesh[j-1] >= ex[i].mesh[j])
424                     {
425                       cerr << "Mesh not in ascending order in exp " << i << ",\n";
426                       return;
427                     }
428                 }
429             }
430           else
431             {
432               cerr << "Incorrect format of mesh in exp " << i << ",\n"; 
433               return;
434             }
435         }
436       
437 #ifdef PRINTDATA
438       *dbg << "\nData: \n";
439       for (j=1; j<=ex[i].nMeasure; j++) {
440         *dbg << ex[i].xMeasure[j];
441         for (k=1; k<=ex[i].nobs; k++)
442           *dbg << "\t" << ex[i].yMeasure[j][k] << "\t" << ex[i].sigma[j][k];
443         *dbg << "\n";
444       }
445       *dbg << "\n";
446 #endif
447       *dbg << "Mesh:";
448       for (k=1;k<=ex[i].nPoints;++k) 
449         *dbg << " " << ex[i].mesh[k];
450       *dbg << "\n";
451     }
452   
453   //initial values
454   for(i=1;i<=globs.nrExp;i++)
455     {
456       *dbg << "\nExperiment " << i << ":\n";
457       if (ex[i].y0) 
458         {
459           *dbg << NEQNS << " starting value(s):";
460           for (k=1; k<=NEQNS; ++k) 
461             *dbg << ' ' << ex[i].y0[k];
462           *dbg << "\n";
463         } 
464       else 
465         {
466           *dbg << "No starting values specified\n";
467         }
468       *dbg << "\n";
469       dbg->flush();
470     } // end for(i=1;i<=globs.nrExp;i++)
471   // starting of the numerics
472 
473   
474    // spline section and initialse
475   
476   if(NSPLINES==0)
477     {
478       globs.initSpline=TRUE;
479       initialise(ex,&globs,FALSE);
480     }
481   else
482     {
483       for(nExp=1;nExp<=globs.nrExp;nExp++)
484         {
485           globs.initSpline=FALSE;
486           initialise(ex,&globs,FALSE);  
487           ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
488           ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
489           ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
490           ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
491           tmp=mxGetField(INEX,nExp-1,"spline"); 
492           for(i=1;i<=NSPLINES;i++)
493             {
494               if(mxGetNumberOfElements(tmp) < NSPLINES)
495                 {
496                   cerr << "Please specify " << NSPLINES << " spline(s).\n";
497                   return;
498                 }
499               
500               splinedat=mxGetCell(tmp,i-1);
501               spline_=mxGetPr(splinedat);
502               ex[nExp].nNodes[i]=mxGetM(splinedat);
503               if(mxGetN(splinedat)!=3)
504                 {
505                   cerr << "Invalid data for spline " << i << ".\n";
506                   return;
507                 }
508               ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
509               ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
510               ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));         
511               
512               //copy spline data
513               for(j=1;j<=ex[nExp].nNodes[i];j++)
514                 {
515                   ex[nExp].splineNodes[i][j]=spline_[j-1];
516                   ex[nExp].splineY[i][j]=spline_[ex[nExp].nNodes[i]+j-1];
517                   ex[nExp].splineGam[i][j]=spline_[2*ex[nExp].nNodes[i]+j-1];       
518                 }
519             }
520         }
521     }
522  
523   //simulate initial state
524   try
525     {
526       if(globs.simInit==TRUE)
527         {
528           for (i=1;i<=globs.nrExp;++i) 
529             simInit(&ex[i],&globs);
530         }
531     }
532   catch (int i)
533     {
534       cerr << "Cannot integrate trajectory, no siminit!\n";
535     }
536 
537   //setting initial values from mesh data
538   for(nExp=1;nExp<=globs.nrExp;nExp++)
539     {
540       tmp=mxGetField(INEX,nExp-1,"mesh");      
541       tmp_=mxGetPr(tmp);
542       if(!mxIsEmpty(tmp))
543         {
544           index=mxGetN(tmp);
545           for(i=1;i<=ex[nExp].nPoints-2;i++)
546             {
547               for(j=1;j<=index-1;j++)
548                 {
549                   ex[nExp].yTry[i+1][j]=tmp_[j*(ex[nExp].nPoints-2)+i-1];
550                 }
551             }
552         }
553     }
554   
555   //*************************
556    try 
557       {
558 
559         if(init!=1)
560           fitIt(ex,&globs);
561         else
562           globs.fitConverged=0;
563       }
564    catch(int i)  //Exception handling
565      {
566        dims[0]=1;
567        dims[1]=1;
568        OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT);
569          //error
570        tmp = mxCreateDoubleMatrix(1,1,mxREAL);
571        tmp_= mxGetPr(tmp);
572        tmp_[0]= 1;
573        mxSetField(OUT,0,"except",tmp);    
574 
575        dims[0]=1;
576        dims[1]=globs.nrExp;
577        OUTEX = mxCreateStructArray(2,dims,sizeof(fieldNamesOUTEX)/sizeof(*fieldNamesOUTEX),fieldNamesOUTEX);
578        if(globs.noGnu==FALSE)
579          {
580            for(j=1;j<=globs.ngnu;j++)
581              pclose(globs.gnuFp[j]);
582          }
583        
584        freeMem(ex,&globs,FALSE);
585        delete ex;
586        return;
587      }
588   
589   //*************************
590       
591   //Output after convergence
592   if(init!=1)
593       {    
594         outFit(ex,&globs);
595         if(!globs.noGnu)
596           system("rm gnuout.dat");
597       }
598   
599   //Mex Output
600   
601   dims[0]=1;
602   dims[1]=1;
603   
604   
605   OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT);
606   
607   //wquer
608   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
609   tmp_= mxGetPr(tmp);
610   tmp_[0]= globs.wquer;
611   mxSetField(OUT,0,"wquer",tmp); 
612   //converged
613   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
614   tmp_= mxGetPr(tmp);
615   tmp_[0]= globs.fitConverged;
616   mxSetField(OUT,0,"converged",tmp);   
617   //chisq
618   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
619   tmp_= mxGetPr(tmp);
620   tmp_[0]= globs.chisq;
621   mxSetField(OUT,0,"chisq",tmp); 
622   //error
623   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
624   tmp_= mxGetPr(tmp);
625   tmp_[0]= 0;
626   mxSetField(OUT,0,"except",tmp); 
627   //Lambda
628   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
629   tmp_= mxGetPr(tmp);
630   tmp_[0]= globs.Lambda;
631   mxSetField(OUT,0,"Lambda",tmp);  
632   //cond
633   tmp = mxCreateDoubleMatrix(1,1,mxREAL);
634   tmp_= mxGetPr(tmp);
635   tmp_[0]= globs.cond;
636   mxSetField(OUT,0,"cond",tmp);  
637 
638   //Experiment Specific output -> OUTEX
639 
640   dims[0]=1;
641   dims[1]=globs.nrExp;
642 
643   OUTEX = mxCreateStructArray(2,dims,sizeof(fieldNamesOUTEX)/sizeof(*fieldNamesOUTEX),fieldNamesOUTEX);
644 
645 
646   for(nExp=1;nExp<=globs.nrExp;nExp++)
647     {
648       //p
649       tmp = mxCreateDoubleMatrix(1,globs.npar,mxREAL);
650       tmp_= mxGetPr(tmp);
651       for(i=1;i<=globs.npar;i++)
652         tmp_[i-1]= ex[nExp].par[i];
653       mxSetField(OUTEX,nExp-1,"p",tmp); 
654       //errP
655       tmp = mxCreateDoubleMatrix(1,globs.npar,mxREAL);
656       tmp_= mxGetPr(tmp);
657       for(i=1;i<=globs.npar;i++)
658         tmp_[i-1]= ex[nExp].errP[i];
659       mxSetField(OUTEX,nExp-1,"errP",tmp); 
660       //y0
661       tmp = mxCreateDoubleMatrix(1,ex[nExp].nvar,mxREAL);
662       tmp_= mxGetPr(tmp);
663       for(i=1;i<=ex[nExp].nvar;i++)
664         tmp_[i-1]= ex[nExp].yTry[1][i];
665       mxSetField(OUTEX,nExp-1,"y0",tmp); 
666       //errY0
667       tmp = mxCreateDoubleMatrix(1,ex[nExp].nvar,mxREAL);
668       tmp_= mxGetPr(tmp);
669       for(i=1;i<=ex[nExp].nvar;i++)
670         tmp_[i-1]= ex[nExp].errY0[i];
671       mxSetField(OUTEX,nExp-1,"errY0",tmp); 
672       //mesh
673       tmp = mxCreateDoubleMatrix(ex[nExp].nPoints-2,ex[nExp].nvar+1,mxREAL);
674       tmp_= mxGetPr(tmp);
675       for(i=1;i<=ex[nExp].nPoints-2;i++)
676         {
677           for(j=1;j<=ex[nExp].nvar+1;j++)
678             {
679               if(j==1)
680                 tmp_[i-1]=ex[nExp].mesh[i+1];
681               else
682                 tmp_[(j-1)*(ex[nExp].nPoints-2)+i-1]=ex[nExp].yTry[i+1][j-1];
683             }
684         }
685       mxSetField(OUTEX,nExp-1,"mesh",tmp); 
686     }
687 
688   freeMem(ex,&globs,FALSE);
689   free_lvector(ndat,1,globs.nrExp);
690   delete ex;
691 
692 }
\end{verbatim}\normalsize 


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




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

References Glob::covar, dmatrix(), FALSE, Glob::npar, Glob::nr\-Exp, Glob\-Exp::nvar, Glob::strategy, and Glob::y0fix.

\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_mex.cc@{diffit\_\-mex.cc}!setMesh@{setMesh}}
\index{setMesh@{setMesh}!diffit_mex.cc@{diffit\_\-mex.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__mex_8cc_343605cc3bf29ae8d91094465c787e95}




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

References dvector(), FALSE, Glob\-Exp::fitend, Glob\-Exp::fitstart, Glob\-Exp::mesh, Glob::no\-Measurements, and Glob\-Exp::n\-Points.

\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_mex.cc@{diffit\_\-mex.cc}!simInit@{simInit}}
\index{simInit@{simInit}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void sim\-Init ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs})}\label{diffit__mex_8cc_1c8ff8f8413349876df41981e468ea78}




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

References dmatrix(), dvector(), Glob\-Exp::n\-Points, Glob\-Exp::nvar, Glob\-Exp::y0, and Glob\-Exp::y\-Try.

\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_mex.cc@{diffit\_\-mex.cc}!dbg@{dbg}}
\index{dbg@{dbg}!diffit_mex.cc@{diffit\_\-mex.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}ofstream$\ast$ {\bf dbg}}\label{diffit__mex_8cc_c6b594e76b642beb861b55fce626b3f1}


debug stream (C++) 



Definition at line 31 of file diffit\_\-mex.cc.