\section{/home/peifer/diffit/modules/initialise.cc File Reference}
\label{initialise_8cc}\index{/home/peifer/diffit/modules/initialise.cc@{/home/peifer/diffit/modules/initialise.cc}}
{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$math.h$>$}\par
{\tt \#include $<$string.h$>$}\par
{\tt \#include $<$stdio.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 EPSTOL}~1e-8
\end{CompactItemize}
\subsection*{Functions}
\begin{CompactItemize}
\item 
void {\bf set\-Initial\-Values} ({\bf Glob\-Exp} $\ast$ex, int n\-P, double $\ast$parameters, double $\ast$y\-Values)
\item 
int {\bf init\-Values\-OK} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex)
\item 
void {\bf setup\-Gnu\-Fp} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} ex[$\,$])
\item 
void {\bf initialise} ({\bf Glob\-Exp} ex[$\,$], {\bf Glob} $\ast$globs, int simit)
\end{CompactItemize}


\subsection{Define Documentation}
\index{initialise.cc@{initialise.cc}!EPSTOL@{EPSTOL}}
\index{EPSTOL@{EPSTOL}!initialise.cc@{initialise.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define EPSTOL~1e-8}\label{initialise_8cc_dc7470da578c77fbfb7aa0151a363483}




Definition at line 20 of file initialise.cc.

Referenced by init\-Values\-OK().

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




Definition at line 98 of file initialise.cc.

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

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


\index{initialise.cc@{initialise.cc}!initValuesOK@{initValuesOK}}
\index{initValuesOK@{initValuesOK}!initialise.cc@{initialise.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}int init\-Values\-OK ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex})}\label{initialise_8cc_dfd0db7554654863d446e22b09f59a27}




Definition at line 23 of file initialise.cc.

References EPSTOL, FALSE, Glob\-Exp::me, Glob\-Exp::mg, Glob\-Exp::r2, Glob\-Exp::r3, and TRUE.

\footnotesize\begin{verbatim}24 {
25   int i,good=TRUE;
26   
27   if (ex->me > 0) 
28     {
29       R2(globs,ex, FALSE);
30       for (i=1; i<=ex->me; ++i)
31         if (good) 
32           good = fabs(ex->r2[i]) < EPSTOL;
33     }
34   
35   if (ex->mg>0) 
36     {
37       R3(globs,ex, FALSE);
38       for (i=1; i<=ex->mg; ++i) 
39         good = good && (ex->r3[i]>=0);
40     }
41   
42   return(good);
43 }
\end{verbatim}\normalsize 


\index{initialise.cc@{initialise.cc}!setInitialValues@{setInitialValues}}
\index{setInitialValues@{setInitialValues}!initialise.cc@{initialise.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void set\-Initial\-Values ({\bf Glob\-Exp} $\ast$ {\em ex}, int {\em n\-P}, double $\ast$ {\em parameters}, double $\ast$ {\em y\-Values})}\label{initialise_8cc_62b20da50932eea5fe90a1d268cf7ea2}




Definition at line 8 of file set\-Initial\-Values.cc.

References Glob\-Exp::mesh, Glob\-Exp::n\-Measure, Glob\-Exp::nobs, Glob\-Exp::n\-Points, Glob\-Exp::nvar, Glob\-Exp::x\-Measure, Glob\-Exp::y\-Measure, and Glob\-Exp::y\-Try.

\footnotesize\begin{verbatim}9 {
10   long i,j,k;
11   double xi, d, w1, w2;
12   long nobs=ex->nobs, nvar=ex->nvar, nMeasure=ex->nMeasure;
13   long nPoints=ex->nPoints;
14   
15   /* set initial values */
16   if (yValues) 
17     {
18       for (j=1; j<=nvar; ++j) 
19         ex->yTry[1][j]=yValues[j];
20       i=2;
21     } 
22   else
23     i=1;
24   k=1;
25   
26   while (i<=nPoints) 
27     {
28       xi=ex->mesh[i];
29       // make educated guess for yTry(xi)
30       
31       // could be improved (polynomial interpolation/extrapolation etc.)
32       
33       if (ex->xMeasure[k]>=xi) 
34         {
35           // before first measurement: take first measured value
36           for (j=1; j<=nvar; ++j)
37             if (j<=nobs)
38               ex->yTry[i][j]=ex->yMeasure[k][j];
39             else
40               ex->yTry[i][j]=1;
41         } else if (k==nMeasure) 
42           {
43             // after last measurement: take last measured value
44             for (j=1; j<=nvar; ++j)
45               if (j<=nobs)
46                 ex->yTry[i][j]=ex->yMeasure[k][j];
47               else
48                 ex->yTry[i][j]=1;
49           } 
50       else 
51         {
52           // interpolate linearly
53           while (k<nMeasure && ex->xMeasure[k]<xi) 
54             ++k;
55           d=ex->xMeasure[k]-ex->xMeasure[k-1];
56           w1=(ex->xMeasure[k]-xi)/d;
57           w2=(xi-ex->xMeasure[k-1])/d;
58           for (j=1; j<=nvar; ++j) 
59             if (j<=nobs)
60               ex->yTry[i][j]=w1*ex->yMeasure[k-1][j] + w2*ex->yMeasure[k][j];
61             else
62               ex->yTry[i][j]=1;
63         }
64       ++i;
65     }  
66 }
\end{verbatim}\normalsize 


\index{initialise.cc@{initialise.cc}!setupGnuFp@{setupGnuFp}}
\index{setupGnuFp@{setupGnuFp}!initialise.cc@{initialise.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void setup\-Gnu\-Fp ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} {\em ex}[$\,$])}\label{initialise_8cc_3ffbf02a86ded32b1308844a24a141a1}




Definition at line 46 of file initialise.cc.

References Glob::nr\-Exp.

\footnotesize\begin{verbatim}47 {
48   long gnu_n=20,ncomp,k;
49   int fullx=700; //max. geometry
50   int fully=500;
51   int nwi,dx=fullx/2,dy,xpos,ypos;
52   char outstr[100];
53 
54 
55   ncomp=0;
56   for(k=1;k<=globs->nrExp;k++)
57     {
58       ncomp+=ex[k].nobs;
59     }
60 
61   // determine nr of windows in y-direction
62   gnu_n = LMIN (gnu_n, ncomp);
63   globs->ngnu=gnu_n;
64 
65   // allocate streams
66   globs->gnuFp = new FILE *[gnu_n+1];
67   
68   if((gnu_n % 2)==1)
69     nwi=(gnu_n+1)/2;
70   else
71     nwi=gnu_n/2;
72   dy=fully/nwi;
73  
74   for(k=1;k<=gnu_n;k++)
75     {
76       if(k==1)
77         {
78           xpos=0;
79           ypos=0;
80         }
81       else if((k % 2)==0)
82         {
83           xpos=dx+10;
84         }
85       else
86         {
87           xpos=0;
88           ypos+=dy+25;
89         }
90       sprintf(outstr,"gnuplot -noraise -geometry %ix%i+%i+%i -title \"%i\"\0",dx,dy,xpos,ypos,k);
91       globs->gnuFp[k]=popen(outstr,"w");
92     }
93 }                               // setupGnuFp
\end{verbatim}\normalsize 


