\section{/home/peifer/diffit/modules/int\-ODE.cc File Reference}
\label{intODE_8cc}\index{/home/peifer/diffit/modules/intODE.cc@{/home/peifer/diffit/modules/intODE.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*{Functions}
\begin{CompactItemize}
\item 
void {\bf integrate\-RK} ({\bf Glob} $\ast$globs, double ystart[$\,$], double p[$\,$], long nvar, double x1, double x2, double eps, double h1, double hmin, long $\ast$nok, long $\ast$nbad, void($\ast$deriv\-Type)(double, double[$\,$], double[$\,$], double[$\,$]))
\item 
void {\bf call\_\-odessa} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, int N, int M\_\-PAR\_\-GES, char $\ast$do\-P, int D\_\-FLAG, int MAX\_\-NR\_\-STEP, int STIFF, int inhomogen, int nmesh, double $\ast$tmesh, double t0, double t1, double $\ast$zustand, double $\ast$$\ast$y, double $\ast$$\ast$$\ast$dmds, double $\ast$$\ast$$\ast$dmdp, double $\ast$$\ast$dyds, double $\ast$$\ast$dydp, double $\ast$parameter, double rtol, double atol, double min\_\-stepsize, double max\_\-stepsize, double initial\_\-stepsize)
\item 
long {\bf max\-Table\-Len} (double $\ast$mesh, long n\-Points, double $\ast$x\-Measure, long n\-Measure)
\item 
void {\bf dcopy} (double $\ast$dest, double $\ast$src, long n)
\item 
void {\bf dfill0} (double $\ast$x, long n)
\item 
void {\bf init\-Yt} (double $\ast$Yt, double $\ast$state)
\item 
double {\bf max} (double a, double b)
\item 
void {\bf sens\-Derivs} (double t, double $\ast$Yt, double $\ast$Ft, double $\ast$p)
\item 
void {\bf tabulate\-Values} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, double t0, double t1, double $\ast$t, long n, double $\ast$state, double $\ast$$\ast$y, double $\ast$$\ast$$\ast$dmds, double $\ast$$\ast$$\ast$dmdp, double $\ast$$\ast$dyds, double $\ast$$\ast$dydp)
\item 
void {\bf int\-ODE} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs, int do\-Derivs, int do\-Plot, long exp\-Nr)
\end{CompactItemize}


\subsection{Function Documentation}
\index{intODE.cc@{int\-ODE.cc}!call_odessa@{call\_\-odessa}}
\index{call_odessa@{call\_\-odessa}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void call\_\-odessa ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex}, int {\em N}, int {\em M\_\-PAR\_\-GES}, char $\ast$ {\em do\-P}, int {\em D\_\-FLAG}, int {\em MAX\_\-NR\_\-STEP}, int {\em STIFF}, int {\em inhomogen}, int {\em nmesh}, double $\ast$ {\em tmesh}, double {\em t0}, double {\em t1}, double $\ast$ {\em zustand}, double $\ast$$\ast$ {\em y}, double $\ast$$\ast$$\ast$ {\em dmds}, double $\ast$$\ast$$\ast$ {\em dmdp}, double $\ast$$\ast$ {\em dyds}, double $\ast$$\ast$ {\em dydp}, double $\ast$ {\em parameter}, double {\em rtol}, double {\em atol}, double {\em min\_\-stepsize}, double {\em max\_\-stepsize}, double {\em initial\_\-stepsize})}\label{intODE_8cc_77d593ebea400902d8585a26dff9e89e}




Definition at line 69 of file call\_\-odessa.cc.

References dmatrix(), dvector(), ivector(), Glob::npar, Glob\-Exp::nvar, and TRUE.

Referenced by tabulate\-Values().

\footnotesize\begin{verbatim}75 {
76 
77 /* from DOKUMENTATION/call_odessa.doc
78 obsolet: N      nVar=anzahl der gleichungen
79 obsolet: M_PAR_GES, nPExp=gesamt_anzahl der parameter
80 char *doP       [1..nPExp] parameter fixen ("0") oder optimieren (sonst)
81 hieraus berechnet:
82     M_PAR_SENS  anzahl zu variierender parameter
83     M_GES =     nVar+M_PAR_SENS; anzahl fit-variablen
84 int D_FLAG      verwende OD_ODE (1) oder AD_ODE (0)
85 int MAX_NR_STEP maximale anzahl der integrationsschritte
86 int STIFF       steifes / nicht-steifes problem (1 / 0)
87 int inhomogen   liegt Inhomogenitätsmatrix vor (1)? bestimme die ableitungen
88                 ansonsten über finite differenzen (0)
89 int nmesh       anzahl der punkte, auf denen der zustand berechnet werden soll
90 double *tmesh           zeitpunkte, an denen der zustand berechnet werden soll
91 obsolet: sens_mesh;     Sensitivitäten werden entweder gar nicht berechnet
92                         (dmds==NULL) oder am Endpunkt UND an den Datenpunkten
93                         (dmds,dmdp,dyds,dydp!=NULL)
94 double t0               Startpunkt, kann < tmesh[1] sein
95 double t1               Endpunkt, kann > tmesh[nmesh] sein
96 double **y              [k=1..nmesh][1..nVar] output: y(tmesh[k])
97 double *zustand         [1..nVar] input: y(t0), output: y(t1)
98 wie in tabVal:
99 double ***dmds          [k=1..nmesh][i=1..nVar][j=1..nVar] dy[i]/dy0[j](jacobi)
100 double ***dmdp          [k=1..nmesh][i=1..nVar][j=1..nPExp] (hiess sensitiv)
101                         output: Sensitivitäten dy[i]/dp[j] bei t=tmesh[k]
102 double **dyds           [i=1..nVar][j=1..nVar]         dy[i]/dy0[j] bei t=t1
103 double **dydp           [i=1..nVar][j=1..nPExp] dy[i]/dp[j]  bei t=t1
104 double *parameter       [1..nPExp] parameter-vektor
105 double rtol             error test passes if yerr<rtol*|y|+atol
106 double atol
107 double min_stepsize     minimale schrittweite
108 double max_stepsize     maximale schrittweite
109 double initial_stepsize
110  */
111 
112   /*==================*/
113   /* lokale variablen */
114   /*==================*/
115 
116   int i = 0, j = 0, k = 0, M_GES, M_PAR_SENS = 0;
117   long nPExp=globs->npar;
118   long nVar=ex->nvar;
119   double * gy=dvector(1,NOBS);
120   int generic=TRUE;
121   
122   /*===============================*/
123   /* berechne M_GES und M_PAR_SENS */
124   /*===============================*/
125 
126   for (i = 1; i <= nPExp; i++)
127     if (doP[i] != '0')
128       M_PAR_SENS++;
129   M_GES = nVar + M_PAR_SENS;
130 
131   /*=======================*/
132   /* stelle NEQ() zusammen */
133   /*=======================*/
134 
135   int IPAR_offset = 10;
136   // wenn diese größe verändert wird, muß berechne_doP_komponente.f
137   // ebenfalls angepaßt werden !!
138 
139   int *NEQ;
140   NEQ = ivector (1, IPAR_offset + nPExp);
141 
142   // stecke den gesamten doP-vektor (in int umgewandelt) in NEQ,
143   // sodass dies global zugänglich ist!
144   for (i = 1; i <= nPExp; i++)
145     NEQ[IPAR_offset + i] = (doP[i] != '0' ? 1 : 0);
146 
147   // alle NEQ-einträge die startwerte betreffend werden
148   // in berechne_komponente auf null gesetzt, sodass keine
149   // inhomogenitäten für diese "Parameter" berechnet werden
150 
151   NEQ[1] = nVar;
152   NEQ[2] = M_GES;               // = nVar + M_PAR_SENS
153   NEQ[3] = M_PAR_SENS;          // Anzahl Fitparameter, <=nPExp
154   NEQ[4] = D_FLAG;
155   NEQ[5] = nPExp;
156 
157   /*===================*/
158   /* stelle Y zusammen */
159   /*===================*/
160 
161   double **Yt = dmatrix (0, M_GES, 1, nVar);
162   //Y = dvector(1,nVar + nVar*M_GES);
163   for (i = 1; i <= nVar; i++)
164     Yt[0][i] = zustand[i];
165 
166   for (j = 1; j <= M_PAR_SENS + nVar; j++)
167     {
168 
169       // initialisiere: Sij = 0 für dy/dp
170 
171       if (j <= M_PAR_SENS)
172         for (i = 1; i <= nVar; i++)
173           Yt[j][i] = 0.0;
174 
175       // initialisiere: Sii = 1 für dyi/dyi°
176 
177       else
178         for (i = 1; i <= nVar; i++)
179           if ((j - M_PAR_SENS) == i)
180             Yt[j][i] = 1.0;
181           else
182             Yt[j][i] = 0.0;
183     }                           //nächster zustand i, der abgeleitet werden soll
184   //Yt[2][2]=0;//!!
185 
186   /*=====================*/
187   /* stelle PAR zusammen */
188   /*=====================*/
189 
190   double *PARAMETER;
191   PARAMETER = dvector (1, nPExp + nVar * nVar);
192   //(1,nPExp) wird von ODESSA erwartet, der rest ist zur vereinfachung
193 
194   //übergebe an PARAMETER alle Parameter!
195   int zaehler = 1;
196   for (i = 1; i <= nPExp; i++)
197     PARAMETER[i] = parameter[i];
198 
199   /*=================================*/
200   /* stelle fehlerparameter zusammen */
201   /*=================================*/
202 
203   int SqualCheck=FALSE;
204   int ITOL = SqualCheck ? 1 : 4;        // erlaubt: 1 or 4
205   int nTOL = (ITOL == 4) ? nVar * (M_GES + 1) : 1;
206 
207   double *RTOL = dvector (1, nTOL);
208   double *ATOL = dvector (1, nTOL);
209   if (ITOL == 4)
210     {
211       for (i = 1; i <= nVar; i++)
212         {
213           RTOL[i] = rtol;       // 1..nVar: rtol
214           ATOL[i] = atol;
215         }
216       for (; i <= nTOL; i++)
217         {
218           RTOL[i] = rtol * 1e8; // nVar+1..: fast ignoriert
219           ATOL[i] = atol * 1e8;
220         }
221     }
222   else
223     {
224       RTOL[1] = rtol;           // alle: rtol
225       ATOL[1] = atol;
226     }
227 
228   /*===================================*/
229   /* stelle kontrollparameter zusammen */
230   /*===================================*/
231 
232   int ITASK = 1;                //siehe ODESSA-beschreibung -> normale integration von T bis TOUT
233   // variable der status-abfrage
234   int ISTATE = 1;               //-> erster aufruf von ODESSA!
235 
236   int *IOPT;
237   IOPT = ivector (1, 3);
238   IOPT[1] = 1;                  //zusätzliche angaben über RWORK-feld
239   IOPT[2] = 1;                  //sensitivitäten werden berechnet !
240   IOPT[3] = inhomogen;          //verwende inhomogenitäten, sonst finite differenzen
241 
242   int MITER = 1;                // -> volle Jacobi-Matrix vorhanden
243   int METH;                     // normaler / steifer Integrator
244   int MAXORD;                   // maximale ordnung der verwendeten approximationen
245 
246   if (STIFF == 0)
247     {
248       METH = 1;
249       MAXORD = 12;
250     }
251   else
252     {
253       METH = 2;
254       MAXORD = 5;
255     }
256   int MF = 10 * METH + MITER;
257   //cout<<"mf = "<<MF<<LF;
258 
259   /*=================================*/
260   /* stelle arbeitsspeicher zusammen */
261   /*=================================*/
262 
263   int LRW;
264   int LIW;
265 
266   if (IOPT[1] == 0)
267     {                           //ohne sensitivitäten
268       LRW = 20 + nVar * (MAXORD + 1) + 3 * nVar + nVar * nVar + 2 + 1;
269       LIW = 20 + nVar;
270     }
271   if (IOPT[1] == 1)
272     {                           //mit sensitivitäten
273       LRW = 20 + nVar * (M_GES + 1) * (MAXORD + 1) + 2 * nVar * (M_GES + 1) +
274         nVar * nVar + 2 + nVar + 1;
275       LIW = 21 + nVar + M_GES;
276     }
277 
278   double *RWORK;
279   RWORK = dvector (1, LRW);
280 
281   int *IWORK;
282   IWORK = ivector (1, LIW);
283 
284   /*==================================*/
285   /* stelle optionalen input zusammen */
286   /*==================================*/
287 
288   // odeint detects stepsize underflow independently of hmin
289   // odessa relies on min_stepsize being large enough to avoid prevent it
290   if (t0 + min_stepsize == t0)
291     {
292       cerr<< t0 << "\t" << min_stepsize << endl;
293       throw 1;
294     }
295   if (t1 + min_stepsize == t0)
296     {
297       cerr<< t0 << "\t" << min_stepsize << endl;
298       throw 1;
299     }
300 
301   RWORK[5] = initial_stepsize;
302   RWORK[6] = max_stepsize;
303   RWORK[7] = min_stepsize;
304 
305   IWORK[5] = MAXORD;
306   IWORK[6] = MAX_NR_STEP;
307   IWORK[7] = 3;                 // not more than 3 warnings of type T=T+H
308 
309   /*==================================*/
310   /* bestimme T und TOUT mittels tmesh */
311   /*==================================*/
312 
313 
314   double T = t0;
315   double TOUT;                  //wird mittels gitter vorgegeben
316 
317 
318   if (dmds)
319     {                           // Sensitivitäten speichern
320       double **dgdy = NULL, **dgdp = NULL;
321       dgdy = dmatrix (1, NOBS, 1, nVar);
322       dgdp = dmatrix (1, NOBS, 1, nPExp);
323       _dfill0 (dgdy[1], NOBS * nVar);
324       _dfill0 (dgdp[1], NOBS * nPExp);
325 
326 
327       for (k = 1; k <= nmesh + 1; k++)
328         {
329           if (k > nmesh)
330             TOUT = t1;          // Endpunkt
331           else
332             TOUT = tmesh[k];    // Datenpunkt
333 
334           /*=================*/
335           /* rufe ODESSA auf */
336           /*=================*/
337 
338           odessa_ (&NEQ[1],
339                    &Yt[0][1],
340                    &PARAMETER[1],
341                    &T, &TOUT, &ITOL, &RTOL[1], &ATOL[1],
342                    &ITASK, &ISTATE, &IOPT[1], &RWORK[1], &LRW, &IWORK[1],
343                    &LIW, &MF);
344 
345           if (ISTATE < 0)
346             break;              // skip remaining k
347 
348           // Zustand speichern
349           int jlong;            // 1..nPExp
350           if (k > nmesh)
351             {                   // Endpunkt
352               _dcopy (zustand, Yt[0], nVar);
353 
354               for (j = 1; j <= M_PAR_SENS; j++)
355                 {
356                   dop_komponente_ (&NEQ[1], j, jlong);
357 
358                   //übernehme die ableitung nach den parametern
359                   for (i = 1; i <= nVar; i++)
360                     dydp[i][jlong] = Yt[j][i];
361                   // jlong and j were exchanged until 3/02
362                 }
363 
364               for (j = 1; j <= nVar; j++)
365                 //übernehme die ableitung nach den startwerten
366                 for (i = 1; i <= nVar; i++)
367                   {
368                     dyds[i][j] = Yt[M_PAR_SENS + j][i];
369                   }
370             }
371           else
372             {                   // k<=nmesh: Datenpunkt
373               _dcopy (y[k], Yt[0], nVar);
374 
375               // write sensitivities to dm/dp and dm/ds
376               // support for extra observables as in odeExpment::tabVal
377               //if (nVarObs > nVar)
378               //  observation (globs,ex,TOUT, y[k], globs->par, dgdy, dgdp);
379               generic=observation(globs,ex,TOUT, y[k],gy, ex->par, dgdy, dgdp);
380               for(i=1;i<=NOBS;i++)
381                 y[k][i]=gy[i];
382               
383               for (i = 1; i <= NOBS; i++)
384                 {
385                   int ivar = i;
386                   for (j = 1; j <= M_PAR_SENS; j++)
387                     {
388                       dop_komponente_ (&NEQ[1], j, jlong);
389                       double &dest = dmdp[k][i][jlong];
390                       if (!generic)
391                         {
392                           dest = dgdp[i][j];
393                           for (long l = 1; l <= nVar; l++)
394                             dest += dgdy[i][l] * Yt[j][l];
395                         }
396                       else
397                         dest = Yt[j][i];
398                     }           // for j
399                   for (j = 1; j <= nVar; j++)
400                     {
401                       double &dest = dmds[k][i][j];
402                       if (!generic)
403                         {
404                           dest = 0;
405                           for (long l = 1; l <= nVar; l++)
406                             dest +=  dgdy[ivar][l] * Yt[M_PAR_SENS + j][l];
407                         }
408                       else
409                         dest = Yt[M_PAR_SENS + j][ivar];
410                     }           // for j
411                 }               // for i
412             }                   // if k>nmesh else
413 
414           T = TOUT;             // für nächstes k
415         }                       //nächster zeitpunkt k
416       //BEGIN eingefuegt von Felix
417       if (dgdy)
418         free_dmatrix (dgdy, 1, NOBS, 1, nVar);
419       if (dgdp)
420         free_dmatrix (dgdp, 1, NOBS, 1, nPExp);
421       //END
422     }
423   else
424     {                           // !dmds
425 
426       for (k = 1; k <= nmesh + 1; k++)
427         {
428           if (k > nmesh)
429             TOUT = t1;          // Endpunkt
430           else
431             TOUT = tmesh[k];    // Datenpunkt
432 
433           odessa_ (&NEQ[1],
434                    &Yt[0][1],
435                    &PARAMETER[1],
436                    &T, &TOUT, &ITOL, &RTOL[1], &ATOL[1],
437                    &ITASK, &ISTATE, &IOPT[1], &RWORK[1], &LRW, &IWORK[1],
438                    &LIW, &MF);
439 
440           if (ISTATE < 0)
441             break;              // skip remaining k
442 
443           // Zustand speichern
444           if (k > nmesh)        // Endpunkt
445             _dcopy (zustand, Yt[0], nVar);
446           else
447             {                   // k<=nmesh: Datenpunkt
448               _dcopy (y[k], Yt[0], nVar);
449               
450               // compute extra observables
451               generic=observation (globs,ex,TOUT, y[k],gy,ex->par, NULL, NULL);
452               for(i=1;i<=NOBS;i++)
453                 y[k][i]=gy[i];
454             }
455           T = TOUT;             // für nächstes k
456         }                       //nächster zeitpunkt k
457 
458     }                           // if dmds else
459 
460 
462 
463   // dvector
464 
465   free_dmatrix (Yt, 0, M_GES, 1, nVar);
466   free_dvector (PARAMETER, 1, nPExp + nVar * nVar);
467 
468   free_dvector (RTOL, 1, nTOL);
469   free_dvector (ATOL, 1, nTOL);
470   free_dvector (RWORK, 1, LRW);
471   free_dvector (gy,1,NOBS);
472 
473   // ivector
474 
475   free_ivector (NEQ, 1, IPAR_offset + nPExp);
476   free_ivector (IWORK, 1, LIW);
477   free_ivector (IOPT, 1, 3);
478 
479   /*===================================================*/
480   /* abfang eventueller fehler / auswertung von ISTATE */
481   /*===================================================*/
482 
483   char *odessaErrStr[6] = {
484     "excess work done on this call (perhaps wrong MF)",
485     "excess accuracy requested (tolerances too small)",
486     "illegal input detected (see printed message)",
487     "repeated error test failures (check all inputs)",
488     "repeated convergence failures\n"
489       "(perhaps bad jacobian supplied or wrong choice of MF or tolerances)",
490     "error weight became zero during problem\n"
491       "(solution component I,J vanished, and ATOL or ATOL(I,J) = 0.0)"
492   };
493 
494   if (ISTATE < 0)
495     {
496       cerr << odessaErrStr[-ISTATE - 1] <<endl;
497       throw 1;
498     }
499 }                              // odeExpment::call_odessa
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!dcopy@{dcopy}}
\index{dcopy@{dcopy}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void dcopy (double $\ast$ {\em dest}, double $\ast$ {\em src}, long {\em n})\hspace{0.3cm}{\tt  [inline]}}\label{intODE_8cc_bcca6016d0bf32216f925814497616b4}




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

Referenced by init\-Yt().

\footnotesize\begin{verbatim}43 {
44   memcpy (dest + 1, src + 1, n * sizeof (double));
45 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!dfill0@{dfill0}}
\index{dfill0@{dfill0}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void dfill0 (double $\ast$ {\em x}, long {\em n})\hspace{0.3cm}{\tt  [inline]}}\label{intODE_8cc_04d1f32a396be53490aaf01dae30c3d3}




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

Referenced by init\-Yt(), sens\-Derivs(), and tabulate\-Values().

\footnotesize\begin{verbatim}48 {
49   memset (x + 1, 0, n * sizeof (double));
50 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!initYt@{initYt}}
\index{initYt@{initYt}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void init\-Yt (double $\ast$ {\em Yt}, double $\ast$ {\em state})}\label{intODE_8cc_da9617ac05c6a795eb7aaf2a714f2bd5}




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

References dcopy(), and dfill0().

Referenced by tabulate\-Values().

\footnotesize\begin{verbatim}53 {
54   // initialise Y
55   if (state)
56     dcopy (Yt, state, NEQNS);    // original state vector
57   dfill0 (Yt + NEQNS, (NPARAMS + NEQNS) * NEQNS);    // 0 in middle and rear part
58   for (long i = 1; i <= NEQNS; i++)       // identity matrix in rear part dy/dy0
59     Yt[(NPARAMS + i) * NEQNS + i] = 1;
60 }                       
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!integrateRK@{integrateRK}}
\index{integrateRK@{integrateRK}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void integrate\-RK ({\bf Glob} $\ast$ {\em globs}, double {\em ystart}[$\,$], double {\em p}[$\,$], long {\em nvar}, double {\em x1}, double {\em x2}, double {\em eps}, double {\em h1}, double {\em hmin}, long $\ast$ {\em nok}, long $\ast$ {\em nbad}, void($\ast$)(double, double[$\,$], double[$\,$], double[$\,$]) {\em deriv\-Type})}\label{intODE_8cc_1794f960789c9073b1fbe822e40728e1}




Definition at line 19 of file integrate\-RK.cc.

References dvector(), Glob::maxstp, nrerror(), and SIGN.

Referenced by tabulate\-Values().

\footnotesize\begin{verbatim}22 {
23   long ode_fkmax = 0, ode_fkount = 0, ode_kmax = 0, ode_kount = 0;
24   double ode_dxsav = 0, *ode_fxp = 0, *ode_xp = 0, **ode_yp = 0;
25   long rkqs_ign = 0;
26   long nstp,i;
27   double xsav,x,hnext,hdid,h;
28   double *yscal,*y,*dydx;
29   long odeint_MAXSTP=globs->maxstp;
30 
31   yscal=dvector(1,nvar);
32   y=dvector(1,nvar);
33   dydx=dvector(1,nvar);
34   x=x1;
35   h=SIGN(h1,x2-x1);
36   *nok = (*nbad) = ode_kount = 0;
37   for (i=1;i<=nvar;i++) y[i]=ystart[i];
38   if (ode_dxsav<0) nrerror("odeint:dxsav<0");
39   if (ode_kmax>0 && ode_kount<ode_kmax-1) {
40      ode_xp[++ode_kount]=x;
41      for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
42      xsav=x;
43   } /* save first point, w.h. 19.03.99 */
44   for (nstp=1;nstp<=odeint_MAXSTP;nstp++) {
45     (*derivs)(x,y,dydx,p);
46     for (i=1;i<=nvar;i++)
47       yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
48     if (ode_fkmax) {
49       if (ode_fkount>=ode_fkmax) nrerror("ode_fkount exceeds ode_fkmax");
50       ode_fxp[++ode_fkount]=x;
51     }
52     if (ode_kmax>0 && ode_kount<ode_kmax-1 && fabs(x-xsav)>ode_dxsav) {
53       ode_xp[++ode_kount]=x;
54       for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
55       xsav=x;
56     }
57     if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
58     rkqs(y,dydx,p,nvar,&x,h,eps,yscal,&hdid,&hnext,globs->rkqs_ign,derivs);
59     if (hdid == h) ++(*nok); 
60     else ++(*nbad);
61     if ((x-x2)*(x2-x1) >= 0.0) {
62       for (i=1;i<=nvar;i++) ystart[i]=y[i];
63       if (ode_kmax>0 && ode_kount<ode_kmax-1) {
64         ode_xp[++ode_kount]=x;
65         for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
66       }
67       free_dvector(dydx,1,nvar);
68       free_dvector(y,1,nvar);
69       free_dvector(yscal,1,nvar);
70       return;
71     }
72     if (fabs(hnext) <= hmin)
73       {
74         cerr << "Step size too small in odeint\n";
75         throw 1;
76       }
77     h=hnext;
78   }
79   cerr << "Too many steps in routine odeint; try -maxstp <n> greater n\n";
80   throw 1;
81 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!intODE@{intODE}}
\index{intODE@{intODE}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void int\-ODE ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs}, int {\em do\-Derivs}, int {\em do\-Plot}, long {\em exp\-Nr})}\label{intODE_8cc_f476965d47230b44b3b9d785a278fd27}




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

References d3tensor(), dmatrix(), Glob\-Exp::dmdp, Glob\-Exp::dmds, dvector(), Glob\-Exp::dydp, Glob\-Exp::dyds, max\-Table\-Len(), Glob\-Exp::mesh, Glob\-Exp::n\-Measure, Glob\-Exp::nobs, Glob::npar, Glob\-Exp::n\-Points, Glob\-Exp::nvar, and Glob\-Exp::x\-Measure.

\footnotesize\begin{verbatim}314 {
315   long i,j,k,ksav,l,m,n,nok,nbad;
316   long nobs=ex->nobs;
317   long nPoints=ex->nPoints, nvar=ex->nvar, nMeasure=ex->nMeasure;
318   double *mesh=ex->mesh, *xMeasure=ex->xMeasure;
319   long tableLen=maxTableLen(mesh, nPoints, xMeasure, nMeasure)+1; 
320   double temp, h,ymin,ymax;
321   
322   // allocate memory
323   double *ys=dvector(1,nvar);
324   double *xTab=dvector(1,tableLen);
325   double **yTab=dmatrix(1,tableLen,1,nvar);
326   double ***dTabds, ***dTabdp;
327   char name[50];
328   ofstream gnuout;
329   ofstream outout;
330 
331 
332   //gnuout.open("gnuout.dat");
333   
334   // compute derivatives?
335   // if not, don't write to graphics either
336   // (no genuine integration, but only estimation of omega)
337   if (doDerivs) 
338     {
339       dTabds=d3tensor(1,tableLen,1,nvar,1,nvar);
340       dTabdp=d3tensor(1,tableLen,1,nvar,1,globs->npar);
341     } 
342   else 
343     {
344       dTabds=NULL; 
345       dTabdp=NULL;
346     }
347   
348 
349   // first mesh point
350   if (doDerivs) {
351     for (i=1; i<=nvar; ++i) 
352       {
353         for (j=1; j<=nvar; ++j)
354           { 
355             ex->dyds[1][i][j] = (i==j) ? 1 : 0;
356             ex->dmds[1][i][j] = (i==j) ? 1 : 0;
357           }
358         for (j=1; j<=globs->npar; ++j)  
359           {
360             ex->dydp[1][i][j] = 0;
361             ex->dmdp[1][i][j] = 0;
362           }
363       }
364   }
365   for (j=1; j<=nvar; ++j) 
366     ex->yComp[1][j]=ex->yTry[1][j];
367 
368 
369   if (doDerivs && mesh[1]==xMeasure[1]) 
370     {
371       observation(globs,ex,mesh[1],ex->yTry[1],ex->yPred[1],ex->par,ex->dmds[1],ex->dmdp[1]);
372     }
373   
374   // integrate
375   k=1;
376   ksav=1;
377   if (mesh[1]==xMeasure[1]) 
378     {
379       k=2; 
380       ksav=2;
381     }
382   for (i=1; i< nPoints; ++i) 
383     {   // given starting values
384       for (j=1; j<=nvar; ++j) 
385         ys[j]=ex->yTry[i][j];  
386       if (xMeasure[k]<mesh[i]) 
387         {
388           cerr << "fatal: intODE: xMeasure table corrupted (mesh: "
389                << mesh[i] << ", xMeasure: " << xMeasure[k] << ")\n";
390           throw 1;
391         }
392       // build table of measuring points
393       n=1;
394       while (k <= nMeasure && xMeasure[k]<=mesh[i+1]) 
395         xTab[n++]=xMeasure[k++];
396       xTab[n]=mesh[i+1];
397 #ifdef PRINTMESH
398       *dbg << "Mesh point " << mesh[i] << ": xTab = (";
399       for (j=1; j<=n; ++j) 
400         *dbg << '\t' << xTab[j];
401       *dbg << ")\n";
402       dbg->flush();
403 #endif      
404       //integration 
405       tabulateValues(globs,ex,mesh[i],mesh[i+1],xTab,n,ys,yTab, dTabds, dTabdp,ex->dyds[i+1],ex->dydp[i+1]);
406 
407 #ifdef PRINTINTEGRATE
408       *dbg << "Integration returned\t  (";
409       for (j=1; j<=n; ++j) 
410         {
411           for (l=1; l<=nvar; ++l) 
412             *dbg << " " << yTab[j][l];
413           *dbg << ",";
414         }
415       *dbg << ")\n";
416       dbg->flush();
417 #endif      
418       for (j=1; j<=ex->nvar; ++j)
419         ex->yComp[i+1][j]=ys[j];   // store value at next mesh point
420 
421       for (j=ksav; j<k; ++j)          // store value at measuring points
422         {
423           for (l=1; l<=ex->nobs; ++l) 
424             {
425               ex->yPred[j][l]=yTab[j-ksav+1][l];
426               if (doDerivs) 
427                 {
428                   for (m=1; m<=ex->nvar; ++m)
429                     ex->dmds[j][l][m]=dTabds[j-ksav+1][l][m];
430                   for (m=1; m<=globs->npar; ++m)
431                     ex->dmdp[j][l][m]=dTabdp[j-ksav+1][l][m];
432                 }
433             }
434         }
435       ksav=k;
436     }
437 
438       //GnuPlotting
439   long mind=2;
440   if(doDerivs && (globs->noGnu==FALSE) && doPlot==TRUE)
441     { 
442       for(l=1;l<=ex->nobs;l++)
443         {
444           if(globs->gnuindex <= globs->ngnu)
445             {
446               mind=2;
447               ofstream gnuout;
448               gnuout.open("gnuout.dat");
449               ymax=ex->yPred[1][l];
450               ymin=ex->yPred[1][l];           
451               for(j=1;j<=ex->nMeasure;j++)
452                 {
453                   gnuout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
454                   if(mesh[mind] <= ex->xMeasure[j])
455                     {
456                       observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
457                       gnuout << mesh[mind] << "\t" << ys[l] << endl << endl;
458                       observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
459                       gnuout << mesh[mind] << "\t" << ys[l] << endl;
460                       mind++;
461                     }
462                   //determine y-range
463                   if(ymax <= ex->yPred[j][l])
464                     ymax=ex->yPred[j][l];
465                   if(ymin >= ex->yPred[j][l])
466                     ymin=ex->yPred[j][l];
467                   if(ymax <= ex->yMeasure[j][l])
468                     ymax=ex->yMeasure[j][l];
469                   if(ymin >= ex->yMeasure[j][l])
470                     ymin=ex->yMeasure[j][l];
471                 }
472               gnuout.flush();
473               gnuout.close();
474               
475               fprintf(globs->gnuFp[globs->gnuindex],"reset\n");
476               fprintf(globs->gnuFp[globs->gnuindex],"set yrange[%f:%f]\n",ymin-0.01*ymin,ymax+0.01*ymax);
477               fprintf(globs->gnuFp[globs->gnuindex],"set xrange[%f:%f]\n",ex->xMeasure[1],ex->xMeasure[ex->nMeasure]);        
478               fflush(globs->gnuFp[globs->gnuindex]);
479               fprintf(globs->gnuFp[globs->gnuindex],"plot \"gnuout.dat\" title \"\"w l 3,\"gnuout.dat\" u 1:3 title \"\" 1 \n");
480               fflush(globs->gnuFp[globs->gnuindex]);
481               //DELAY
482               for(long dd=1;dd<=1000000;dd++)
483                1+1;
484             }
485           sprintf(name,"Gnuout_Exp%d_Obs%d.dat",expNr,l);
486           outout.open(name);
487           
488           globs->gnuindex++;
489           mind=2;
490           for(j=1;j<=ex->nMeasure;j++)
491             {
492               
493               outout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
494               if(mesh[mind] <= ex->xMeasure[j])
495                 {
496                   observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
497                   outout << mesh[mind] << "\t" << ys[l] << endl << endl;
498                   observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
499                   outout << mesh[mind] << "\t" << ys[l] << endl;
500                   mind++;
501                 }
502             }
503           outout.flush();
504           outout.close();
505           for(long dd=1;dd<=1000000;dd++)
506             1+1; 
507         }
508     }
509   
510   if (doDerivs) {
511     free_d3tensor(dTabdp,1,tableLen,1,nvar,1,globs->npar);
512     free_d3tensor(dTabds,1,tableLen,1,nvar,1,nvar);
513   }
514   free_dmatrix(yTab,1,tableLen,1,nvar);
515   free_dvector(xTab,1,tableLen);
516   free_dvector(ys,1,nvar);
517 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!max@{max}}
\index{max@{max}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}double max (double {\em a}, double {\em b})\hspace{0.3cm}{\tt  [inline]}}\label{intODE_8cc_452d2724d2f82cc2c17f54b388232200}




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

Referenced by max\-Table\-Len(), and tabulate\-Values().

\footnotesize\begin{verbatim}63 {
64   if(a>=b)
65     return(a);
66   else
67     return(b);
68 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!maxTableLen@{maxTableLen}}
\index{maxTableLen@{maxTableLen}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}long max\-Table\-Len (double $\ast$ {\em mesh}, long {\em n\-Points}, double $\ast$ {\em x\-Measure}, long {\em n\-Measure})}\label{intODE_8cc_1f76784558ce9eb9fbd2db9f50c68fe7}




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

References max().

Referenced by int\-ODE().

\footnotesize\begin{verbatim}36 {
37     int i,j, max=0;
38 
39     return(nPoints > nMeasure ? nPoints : nMeasure);  // to be improved
40 }
\end{verbatim}\normalsize 


\index{intODE.cc@{int\-ODE.cc}!sensDerivs@{sensDerivs}}
\index{sensDerivs@{sensDerivs}!intODE.cc@{int\-ODE.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void sens\-Derivs (double {\em t}, double $\ast$ {\em Yt}, double $\ast$ {\em Ft}, double $\ast$ {\em p})}\label{intODE_8cc_04f1a9061c83a6e0996603630ded0547}




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

References dfill0(), and dmatrix().

Referenced by f(), and tabulate\-Values().

\footnotesize\begin{verbatim}73 {
74   //  sensitivity equations
75   // P=(p,y0)
76   // S=dy/dP
77   // Y=(y,S)
78   // F=dY/dt=(f,df/dy*S+df/dP)
79   // Y and F are stored in column major format and rolled out as Yt and Ft
80 
81   long i, j, k;
82   dfill0 (Ft, (1 + NPARAMS + NEQNS) * NEQNS);       // initialize to zero
83   // ode, inhomo and jacobi set only non-zero elements
84   // inhome expects a dmatrix; we create one who's data lies in Ft
85   double **dfdpt = new double *[NPARAMS] - 1;
86   for (i = 1; i <= NPARAMS; i++)
87     dfdpt[i] = Ft + i * NEQNS;
88   ode (t, Yt, Ft, p);           // first part of Y,F are y,f
89   inhomo (t, Yt, dfdpt, p);     // fill df/dp into F
90   delete[]++ dfdpt;
91 
92   // now F=df/dP
93   // add M*S, M=df/dy
94   double **M = dmatrix (1, NEQNS, 1, NEQNS);
95   dfill0 (M[1], NEQNS * NEQNS);
96   jacobi (t, Yt, M, p);         // get M
97   for (i = 1; i <= NEQNS; i++)
98     for (j = 1; j <= NPARAMS + NEQNS; j++)
99       {
100         double &sum = Ft[j * NEQNS + i];
101         for (k = 1; k <= NEQNS; k++)
102           sum += M[k][i] * Yt[j * NEQNS + k];
103       }                         // for i,j
104   free_dmatrix (M, 1, NEQNS, 1, NEQNS);
105  }                               // sensDerivs
\end{verbatim}\normalsize 


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




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

References call\_\-odessa(), dfill0(), dmatrix(), Glob::do\-P, dvector(), Glob::eps, init\-Yt(), integrate\-RK(), Glob::integrator, max(), Glob::maxstp, Glob::npar, Glob\-Exp::nvar, Glob\-Exp::par, Glob::rkqs\_\-ign, sens\-Derivs(), Glob::stiff, and TRUE.

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


