- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/modules/intODE.cc File Reference

`#include <iostream>`
`#include <fstream>`
`#include <math.h>`
`#include <string.h>`
`#include <stdio.h>`
`#include "../def.h"`
`#include "../model.h"`
`#include "../nr.h"`

Go to the source code of this file.

## Functions
void | integrateRK (Glob *globs, double ystart[], double p[], long nvar, double x1, double x2, double eps, double h1, double hmin, long *nok, long *nbad, void(*derivType)(double, double[], double[], double[]))
void | call_odessa (Glob *globs, GlobExp *ex, int N, int M_PAR_GES, char *doP, int D_FLAG, int MAX_NR_STEP, int STIFF, int inhomogen, int nmesh, double *tmesh, double t0, double t1, double *zustand, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp, double *parameter, double rtol, double atol, double min_stepsize, double max_stepsize, double initial_stepsize)
long | maxTableLen (double *mesh, long nPoints, double *xMeasure, long nMeasure)
void | dcopy (double *dest, double *src, long n)
void | dfill0 (double *x, long n)
void | initYt (double *Yt, double *state)
double | max (double a, double b)
void | sensDerivs (double t, double *Yt, double *Ft, double *p)
void | tabulateValues (Glob *globs, GlobExp *ex, double t0, double t1, double *t, long n, double *state, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp)
void | intODE (GlobExp *ex, Glob *globs, int doDerivs, int doPlot, long expNr)

## Function Documentation

void call_odessa | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | int | N,
|  | int | M_PAR_GES,
|  | char * | doP,
|  | int | D_FLAG,
|  | int | MAX_NR_STEP,
|  | int | STIFF,
|  | int | inhomogen,
|  | int | nmesh,
|  | double * | tmesh,
|  | double | t0,
|  | double | t1,
|  | double * | zustand,
|  | double ** | y,
|  | double *** | dmds,
|  | double *** | dmdp,
|  | double ** | dyds,
|  | double ** | dydp,
|  | double * | parameter,
|  | double | rtol,
|  | double | atol,
|  | double | min_stepsize,
|  | double | max_stepsize,
|  | double | initial_stepsize
| ) |

| Definition at line 69 of file call_odessa.cc. References dmatrix(), dvector(), ivector(), Glob::npar, GlobExp::nvar, and TRUE. Referenced by tabulateValues().
    00075 {
    00076
    00077 /* from DOKUMENTATION/call_odessa.doc
    00078 obsolet: N      nVar=anzahl der gleichungen
    00079 obsolet: M_PAR_GES, nPExp=gesamt_anzahl der parameter
    00080 char *doP       [1..nPExp] parameter fixen ("0") oder optimieren (sonst)
    00081 hieraus berechnet:
    00082     M_PAR_SENS  anzahl zu variierender parameter
    00083     M_GES =     nVar+M_PAR_SENS; anzahl fit-variablen
    00084 int D_FLAG      verwende OD_ODE (1) oder AD_ODE (0)
    00085 int MAX_NR_STEP maximale anzahl der integrationsschritte
    00086 int STIFF       steifes / nicht-steifes problem (1 / 0)
    00087 int inhomogen   liegt Inhomogenitätsmatrix vor (1)? bestimme die ableitungen
    00088                 ansonsten über finite differenzen (0)
    00089 int nmesh       anzahl der punkte, auf denen der zustand berechnet werden soll
    00090 double *tmesh           zeitpunkte, an denen der zustand berechnet werden soll
    00091 obsolet: sens_mesh;     Sensitivitäten werden entweder gar nicht berechnet
    00092                         (dmds==NULL) oder am Endpunkt UND an den Datenpunkten
    00093                         (dmds,dmdp,dyds,dydp!=NULL)
    00094 double t0               Startpunkt, kann < tmesh[1] sein
    00095 double t1               Endpunkt, kann > tmesh[nmesh] sein
    00096 double **y              [k=1..nmesh][1..nVar] output: y(tmesh[k])
    00097 double *zustand         [1..nVar] input: y(t0), output: y(t1)
    00098 wie in tabVal:
    00099 double ***dmds          [k=1..nmesh][i=1..nVar][j=1..nVar] dy[i]/dy0[j](jacobi)
    00100 double ***dmdp          [k=1..nmesh][i=1..nVar][j=1..nPExp] (hiess sensitiv)
    00101                         output: Sensitivitäten dy[i]/dp[j] bei t=tmesh[k]
    00102 double **dyds           [i=1..nVar][j=1..nVar]         dy[i]/dy0[j] bei t=t1
    00103 double **dydp           [i=1..nVar][j=1..nPExp] dy[i]/dp[j]  bei t=t1
    00104 double *parameter       [1..nPExp] parameter-vektor
    00105 double rtol             error test passes if yerr<rtol*|y|+atol
    00106 double atol
    00107 double min_stepsize     minimale schrittweite
    00108 double max_stepsize     maximale schrittweite
    00109 double initial_stepsize
    00110  */
    00111
    00112   /*==================*/
    00113   /* lokale variablen */
    00114   /*==================*/
    00115
    00116   int i = 0, j = 0, k = 0, M_GES, M_PAR_SENS = 0;
    00117   long nPExp=globs->npar;
    00118   long nVar=ex->nvar;
    00119   double * gy=dvector(1,NOBS);
    00120   int generic=TRUE;
    00121
    00122   /*===============================*/
    00123   /* berechne M_GES und M_PAR_SENS */
    00124   /*===============================*/
    00125
    00126   for (i = 1; i <= nPExp; i++)
    00127     if (doP[i] != '0')
    00128       M_PAR_SENS++;
    00129   M_GES = nVar + M_PAR_SENS;
    00130
    00131   /*=======================*/
    00132   /* stelle NEQ() zusammen */
    00133   /*=======================*/
    00134
    00135   int IPAR_offset = 10;
    00136   // wenn diese größe verändert wird, muß berechne_doP_komponente.f
    00137   // ebenfalls angepaßt werden !!
    00138
    00139   int *NEQ;
    00140   NEQ = ivector (1, IPAR_offset + nPExp);
    00141
    00142   // stecke den gesamten doP-vektor (in int umgewandelt) in NEQ,
    00143   // sodass dies global zugänglich ist!
    00144   for (i = 1; i <= nPExp; i++)
    00145     NEQ[IPAR_offset + i] = (doP[i] != '0' ? 1 : 0);
    00146
    00147   // alle NEQ-einträge die startwerte betreffend werden
    00148   // in berechne_komponente auf null gesetzt, sodass keine
    00149   // inhomogenitäten für diese "Parameter" berechnet werden
    00150
    00151   NEQ[1] = nVar;
    00152   NEQ[2] = M_GES;               // = nVar + M_PAR_SENS
    00153   NEQ[3] = M_PAR_SENS;          // Anzahl Fitparameter, <=nPExp
    00154   NEQ[4] = D_FLAG;
    00155   NEQ[5] = nPExp;
    00156
    00157   /*===================*/
    00158   /* stelle Y zusammen */
    00159   /*===================*/
    00160
    00161   double **Yt = dmatrix (0, M_GES, 1, nVar);
    00162   //Y = dvector(1,nVar + nVar*M_GES);
    00163   for (i = 1; i <= nVar; i++)
    00164     Yt[0][i] = zustand[i];
    00165
    00166   for (j = 1; j <= M_PAR_SENS + nVar; j++)
    00167     {
    00168
    00169       // initialisiere: Sij = 0 für dy/dp
    00170
    00171       if (j <= M_PAR_SENS)
    00172         for (i = 1; i <= nVar; i++)
    00173           Yt[j][i] = 0.0;
    00174
    00175       // initialisiere: Sii = 1 für dyi/dyi°
    00176
    00177       else
    00178         for (i = 1; i <= nVar; i++)
    00179           if ((j - M_PAR_SENS) == i)
    00180             Yt[j][i] = 1.0;
    00181           else
    00182             Yt[j][i] = 0.0;
    00183     }                           //nächster zustand i, der abgeleitet werden soll
    00184   //Yt[2][2]=0;//!!
    00185
    00186   /*=====================*/
    00187   /* stelle PAR zusammen */
    00188   /*=====================*/
    00189
    00190   double *PARAMETER;
    00191   PARAMETER = dvector (1, nPExp + nVar * nVar);
    00192   //(1,nPExp) wird von ODESSA erwartet, der rest ist zur vereinfachung
    00193
    00194   //übergebe an PARAMETER alle Parameter!
    00195   int zaehler = 1;
    00196   for (i = 1; i <= nPExp; i++)
    00197     PARAMETER[i] = parameter[i];
    00198
    00199   /*=================================*/
    00200   /* stelle fehlerparameter zusammen */
    00201   /*=================================*/
    00202
    00203   int SqualCheck=FALSE;
    00204   int ITOL = SqualCheck ? 1 : 4;        // erlaubt: 1 or 4
    00205   int nTOL = (ITOL == 4) ? nVar * (M_GES + 1) : 1;
    00206
    00207   double *RTOL = dvector (1, nTOL);
    00208   double *ATOL = dvector (1, nTOL);
    00209   if (ITOL == 4)
    00210     {
    00211       for (i = 1; i <= nVar; i++)
    00212         {
    00213           RTOL[i] = rtol;       // 1..nVar: rtol
    00214           ATOL[i] = atol;
    00215         }
    00216       for (; i <= nTOL; i++)
    00217         {
    00218           RTOL[i] = rtol * 1e8; // nVar+1..: fast ignoriert
    00219           ATOL[i] = atol * 1e8;
    00220         }
    00221     }
    00222   else
    00223     {
    00224       RTOL[1] = rtol;           // alle: rtol
    00225       ATOL[1] = atol;
    00226     }
    00227
    00228   /*===================================*/
    00229   /* stelle kontrollparameter zusammen */
    00230   /*===================================*/
    00231
    00232   int ITASK = 1;                //siehe ODESSA-beschreibung -> normale integration von T bis TOUT
    00233   // variable der status-abfrage
    00234   int ISTATE = 1;               //-> erster aufruf von ODESSA!
    00235
    00236   int *IOPT;
    00237   IOPT = ivector (1, 3);
    00238   IOPT[1] = 1;                  //zusätzliche angaben über RWORK-feld
    00239   IOPT[2] = 1;                  //sensitivitäten werden berechnet !
    00240   IOPT[3] = inhomogen;          //verwende inhomogenitäten, sonst finite differenzen
    00241
    00242   int MITER = 1;                // -> volle Jacobi-Matrix vorhanden
    00243   int METH;                     // normaler / steifer Integrator
    00244   int MAXORD;                   // maximale ordnung der verwendeten approximationen
    00245
    00246   if (STIFF == 0)
    00247     {
    00248       METH = 1;
    00249       MAXORD = 12;
    00250     }
    00251   else
    00252     {
    00253       METH = 2;
    00254       MAXORD = 5;
    00255     }
    00256   int MF = 10 * METH + MITER;
    00257   //cout<<"mf = "<<MF<<LF;
    00258
    00259   /*=================================*/
    00260   /* stelle arbeitsspeicher zusammen */
    00261   /*=================================*/
    00262
    00263   int LRW;
    00264   int LIW;
    00265
    00266   if (IOPT[1] == 0)
    00267     {                           //ohne sensitivitäten
    00268       LRW = 20 + nVar * (MAXORD + 1) + 3 * nVar + nVar * nVar + 2 + 1;
    00269       LIW = 20 + nVar;
    00270     }
    00271   if (IOPT[1] == 1)
    00272     {                           //mit sensitivitäten
    00273       LRW = 20 + nVar * (M_GES + 1) * (MAXORD + 1) + 2 * nVar * (M_GES + 1) +
    00274         nVar * nVar + 2 + nVar + 1;
    00275       LIW = 21 + nVar + M_GES;
    00276     }
    00277
    00278   double *RWORK;
    00279   RWORK = dvector (1, LRW);
    00280
    00281   int *IWORK;
    00282   IWORK = ivector (1, LIW);
    00283
    00284   /*==================================*/
    00285   /* stelle optionalen input zusammen */
    00286   /*==================================*/
    00287
    00288   // odeint detects stepsize underflow independently of hmin
    00289   // odessa relies on min_stepsize being large enough to avoid prevent it
    00290   if (t0 + min_stepsize == t0)
    00291     {
    00292       cerr<< t0 << "\t" << min_stepsize << endl;
    00293       throw 1;
    00294     }
    00295   if (t1 + min_stepsize == t0)
    00296     {
    00297       cerr<< t0 << "\t" << min_stepsize << endl;
    00298       throw 1;
    00299     }
    00300
    00301   RWORK[5] = initial_stepsize;
    00302   RWORK[6] = max_stepsize;
    00303   RWORK[7] = min_stepsize;
    00304
    00305   IWORK[5] = MAXORD;
    00306   IWORK[6] = MAX_NR_STEP;
    00307   IWORK[7] = 3;                 // not more than 3 warnings of type T=T+H
    00308
    00309   /*==================================*/
    00310   /* bestimme T und TOUT mittels tmesh */
    00311   /*==================================*/
    00312
    00313
    00314   double T = t0;
    00315   double TOUT;                  //wird mittels gitter vorgegeben
    00316
    00317
    00318   if (dmds)
    00319     {                           // Sensitivitäten speichern
    00320       double **dgdy = NULL, **dgdp = NULL;
    00321       dgdy = dmatrix (1, NOBS, 1, nVar);
    00322       dgdp = dmatrix (1, NOBS, 1, nPExp);
    00323       _dfill0 (dgdy[1], NOBS * nVar);
    00324       _dfill0 (dgdp[1], NOBS * nPExp);
    00325
    00326
    00327       for (k = 1; k <= nmesh + 1; k++)
    00328         {
    00329           if (k > nmesh)
    00330             TOUT = t1;          // Endpunkt
    00331           else
    00332             TOUT = tmesh[k];    // Datenpunkt
    00333
    00334           /*=================*/
    00335           /* rufe ODESSA auf */
    00336           /*=================*/
    00337
    00338           odessa_ (&NEQ[1],
    00339                    &Yt[0][1],
    00340                    &PARAMETER[1],
    00341                    &T, &TOUT, &ITOL, &RTOL[1], &ATOL[1],
    00342                    &ITASK, &ISTATE, &IOPT[1], &RWORK[1], &LRW, &IWORK[1],
    00343                    &LIW, &MF);
    00344
    00345           if (ISTATE < 0)
    00346             break;              // skip remaining k
    00347
    00348           // Zustand speichern
    00349           int jlong;            // 1..nPExp
    00350           if (k > nmesh)
    00351             {                   // Endpunkt
    00352               _dcopy (zustand, Yt[0], nVar);
    00353
    00354               for (j = 1; j <= M_PAR_SENS; j++)
    00355                 {
    00356                   dop_komponente_ (&NEQ[1], j, jlong);
    00357
    00358                   //übernehme die ableitung nach den parametern
    00359                   for (i = 1; i <= nVar; i++)
    00360                     dydp[i][jlong] = Yt[j][i];
    00361                   // jlong and j were exchanged until 3/02
    00362                 }
    00363
    00364               for (j = 1; j <= nVar; j++)
    00365                 //übernehme die ableitung nach den startwerten
    00366                 for (i = 1; i <= nVar; i++)
    00367                   {
    00368                     dyds[i][j] = Yt[M_PAR_SENS + j][i];
    00369                   }
    00370             }
    00371           else
    00372             {                   // k<=nmesh: Datenpunkt
    00373               _dcopy (y[k], Yt[0], nVar);
    00374
    00375               // write sensitivities to dm/dp and dm/ds
    00376               // support for extra observables as in odeExpment::tabVal
    00377               //if (nVarObs > nVar)
    00378               //  observation (globs,ex,TOUT, y[k], globs->par, dgdy, dgdp);
    00379               generic=observation(globs,ex,TOUT, y[k],gy, ex->par, dgdy, dgdp);
    00380               for(i=1;i<=NOBS;i++)
    00381                 y[k][i]=gy[i];
    00382
    00383               for (i = 1; i <= NOBS; i++)
    00384                 {
    00385                   int ivar = i;
    00386                   for (j = 1; j <= M_PAR_SENS; j++)
    00387                     {
    00388                       dop_komponente_ (&NEQ[1], j, jlong);
    00389                       double &dest = dmdp[k][i][jlong];
    00390                       if (!generic)
    00391                         {
    00392                           dest = dgdp[i][j];
    00393                           for (long l = 1; l <= nVar; l++)
    00394                             dest += dgdy[i][l] * Yt[j][l];
    00395                         }
    00396                       else
    00397                         dest = Yt[j][i];
    00398                     }           // for j
    00399                   for (j = 1; j <= nVar; j++)
    00400                     {
    00401                       double &dest = dmds[k][i][j];
    00402                       if (!generic)
    00403                         {
    00404                           dest = 0;
    00405                           for (long l = 1; l <= nVar; l++)
    00406                             dest +=  dgdy[ivar][l] * Yt[M_PAR_SENS + j][l];
    00407                         }
    00408                       else
    00409                         dest = Yt[M_PAR_SENS + j][ivar];
    00410                     }           // for j
    00411                 }               // for i
    00412             }                   // if k>nmesh else
    00413
    00414           T = TOUT;             // für nächstes k
    00415         }                       //nächster zeitpunkt k
    00416       //BEGIN eingefuegt von Felix
    00417       if (dgdy)
    00418         free_dmatrix (dgdy, 1, NOBS, 1, nVar);
    00419       if (dgdp)
    00420         free_dmatrix (dgdp, 1, NOBS, 1, nPExp);
    00421       //END
    00422     }
    00423   else
    00424     {                           // !dmds
    00425
    00426       for (k = 1; k <= nmesh + 1; k++)
    00427         {
    00428           if (k > nmesh)
    00429             TOUT = t1;          // Endpunkt
    00430           else
    00431             TOUT = tmesh[k];    // Datenpunkt
    00432
    00433           odessa_ (&NEQ[1],
    00434                    &Yt[0][1],
    00435                    &PARAMETER[1],
    00436                    &T, &TOUT, &ITOL, &RTOL[1], &ATOL[1],
    00437                    &ITASK, &ISTATE, &IOPT[1], &RWORK[1], &LRW, &IWORK[1],
    00438                    &LIW, &MF);
    00439
    00440           if (ISTATE < 0)
    00441             break;              // skip remaining k
    00442
    00443           // Zustand speichern
    00444           if (k > nmesh)        // Endpunkt
    00445             _dcopy (zustand, Yt[0], nVar);
    00446           else
    00447             {                   // k<=nmesh: Datenpunkt
    00448               _dcopy (y[k], Yt[0], nVar);
    00449
    00450               // compute extra observables
    00451               generic=observation (globs,ex,TOUT, y[k],gy,ex->par, NULL, NULL);
    00452               for(i=1;i<=NOBS;i++)
    00453                 y[k][i]=gy[i];
    00454             }
    00455           T = TOUT;             // für nächstes k
    00456         }                       //nächster zeitpunkt k
    00457
    00458     }                           // if dmds else
    00459
    00460
    00462
    00463   // dvector
    00464
    00465   free_dmatrix (Yt, 0, M_GES, 1, nVar);
    00466   free_dvector (PARAMETER, 1, nPExp + nVar * nVar);
    00467
    00468   free_dvector (RTOL, 1, nTOL);
    00469   free_dvector (ATOL, 1, nTOL);
    00470   free_dvector (RWORK, 1, LRW);
    00471   free_dvector (gy,1,NOBS);
    00472
    00473   // ivector
    00474
    00475   free_ivector (NEQ, 1, IPAR_offset + nPExp);
    00476   free_ivector (IWORK, 1, LIW);
    00477   free_ivector (IOPT, 1, 3);
    00478
    00479   /*===================================================*/
    00480   /* abfang eventueller fehler / auswertung von ISTATE */
    00481   /*===================================================*/
    00482
    00483   char *odessaErrStr[6] = {
    00484     "excess work done on this call (perhaps wrong MF)",
    00485     "excess accuracy requested (tolerances too small)",
    00486     "illegal input detected (see printed message)",
    00487     "repeated error test failures (check all inputs)",
    00488     "repeated convergence failures\n"
    00489       "(perhaps bad jacobian supplied or wrong choice of MF or tolerances)",
    00490     "error weight became zero during problem\n"
    00491       "(solution component I,J vanished, and ATOL or ATOL(I,J) = 0.0)"
    00492   };
    00493
    00494   if (ISTATE < 0)
    00495     {
    00496       cerr << odessaErrStr[-ISTATE - 1] <<endl;
    00497       throw 1;
    00498     }
    00499 }                              // odeExpment::call_odessa

void dcopy | ( | double * | dest,
|  | double * | src,
|  | long | n
| ) | `[inline]`

| Definition at line 42 of file intODE.cc. Referenced by initYt().
    00043 {
    00044   memcpy (dest + 1, src + 1, n * sizeof (double));
    00045 }

void dfill0 | ( | double * | x,
|  | long | n
| ) | `[inline]`

| Definition at line 47 of file intODE.cc. Referenced by initYt(), sensDerivs(), and tabulateValues().
    00048 {
    00049   memset (x + 1, 0, n * sizeof (double));
    00050 }

void initYt | ( | double * | Yt,
|  | double * | state
| ) |

| Definition at line 52 of file intODE.cc. References dcopy(), and dfill0(). Referenced by tabulateValues().
    00053 {
    00054   // initialise Y
    00055   if (state)
    00056     dcopy (Yt, state, NEQNS);    // original state vector
    00057   dfill0 (Yt + NEQNS, (NPARAMS + NEQNS) * NEQNS);    // 0 in middle and rear part
    00058   for (long i = 1; i <= NEQNS; i++)       // identity matrix in rear part dy/dy0
    00059     Yt[(NPARAMS + i) * NEQNS + i] = 1;
    00060 }

void integrateRK | ( | Glob * | globs,
|  | double | ystart[],
|  | double | p[],
|  | long | nvar,
|  | double | x1,
|  | double | x2,
|  | double | eps,
|  | double | h1,
|  | double | hmin,
|  | long * | nok,
|  | long * | nbad,
|  | void(*)(double, double[], double[], double[]) | derivType
| ) |

| Definition at line 19 of file integrateRK.cc. References dvector(), Glob::maxstp, nrerror(), and SIGN. Referenced by tabulateValues().
    00022 {
    00023   long ode_fkmax = 0, ode_fkount = 0, ode_kmax = 0, ode_kount = 0;
    00024   double ode_dxsav = 0, *ode_fxp = 0, *ode_xp = 0, **ode_yp = 0;
    00025   long rkqs_ign = 0;
    00026   long nstp,i;
    00027   double xsav,x,hnext,hdid,h;
    00028   double *yscal,*y,*dydx;
    00029   long odeint_MAXSTP=globs->maxstp;
    00030
    00031   yscal=dvector(1,nvar);
    00032   y=dvector(1,nvar);
    00033   dydx=dvector(1,nvar);
    00034   x=x1;
    00035   h=SIGN(h1,x2-x1);
    00036   *nok = (*nbad) = ode_kount = 0;
    00037   for (i=1;i<=nvar;i++) y[i]=ystart[i];
    00038   if (ode_dxsav<0) nrerror("odeint:dxsav<0");
    00039   if (ode_kmax>0 && ode_kount<ode_kmax-1) {
    00040      ode_xp[++ode_kount]=x;
    00041      for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
    00042      xsav=x;
    00043   } /* save first point, w.h. 19.03.99 */
    00044   for (nstp=1;nstp<=odeint_MAXSTP;nstp++) {
    00045     (*derivs)(x,y,dydx,p);
    00046     for (i=1;i<=nvar;i++)
    00047       yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
    00048     if (ode_fkmax) {
    00049       if (ode_fkount>=ode_fkmax) nrerror("ode_fkount exceeds ode_fkmax");
    00050       ode_fxp[++ode_fkount]=x;
    00051     }
    00052     if (ode_kmax>0 && ode_kount<ode_kmax-1 && fabs(x-xsav)>ode_dxsav) {
    00053       ode_xp[++ode_kount]=x;
    00054       for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
    00055       xsav=x;
    00056     }
    00057     if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
    00058     rkqs(y,dydx,p,nvar,&x,h,eps,yscal,&hdid,&hnext,globs->rkqs_ign,derivs);
    00059     if (hdid == h) ++(*nok);
    00060     else ++(*nbad);
    00061     if ((x-x2)*(x2-x1) >= 0.0) {
    00062       for (i=1;i<=nvar;i++) ystart[i]=y[i];
    00063       if (ode_kmax>0 && ode_kount<ode_kmax-1) {
    00064         ode_xp[++ode_kount]=x;
    00065         for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
    00066       }
    00067       free_dvector(dydx,1,nvar);
    00068       free_dvector(y,1,nvar);
    00069       free_dvector(yscal,1,nvar);
    00070       return;
    00071     }
    00072     if (fabs(hnext) <= hmin)
    00073       {
    00074         cerr << "Step size too small in odeint\n";
    00075         throw 1;
    00076       }
    00077     h=hnext;
    00078   }
    00079   cerr << "Too many steps in routine odeint; try -maxstp <n> greater n\n";
    00080   throw 1;
    00081 }

void intODE | ( | GlobExp * | ex,
|  | Glob * | globs,
|  | int | doDerivs,
|  | int | doPlot,
|  | long | expNr
| ) |

| Definition at line 313 of file intODE.cc. References d3tensor(), dmatrix(), GlobExp::dmdp, GlobExp::dmds, dvector(), GlobExp::dydp, GlobExp::dyds, maxTableLen(), GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, GlobExp::nvar, and GlobExp::xMeasure.
    00314 {
    00315   long i,j,k,ksav,l,m,n,nok,nbad;
    00316   long nobs=ex->nobs;
    00317   long nPoints=ex->nPoints, nvar=ex->nvar, nMeasure=ex->nMeasure;
    00318   double *mesh=ex->mesh, *xMeasure=ex->xMeasure;
    00319   long tableLen=maxTableLen(mesh, nPoints, xMeasure, nMeasure)+1;
    00320   double temp, h,ymin,ymax;
    00321
    00322   // allocate memory
    00323   double *ys=dvector(1,nvar);
    00324   double *xTab=dvector(1,tableLen);
    00325   double **yTab=dmatrix(1,tableLen,1,nvar);
    00326   double ***dTabds, ***dTabdp;
    00327   char name[50];
    00328   ofstream gnuout;
    00329   ofstream outout;
    00330
    00331
    00332   //gnuout.open("gnuout.dat");
    00333
    00334   // compute derivatives?
    00335   // if not, don't write to graphics either
    00336   // (no genuine integration, but only estimation of omega)
    00337   if (doDerivs)
    00338     {
    00339       dTabds=d3tensor(1,tableLen,1,nvar,1,nvar);
    00340       dTabdp=d3tensor(1,tableLen,1,nvar,1,globs->npar);
    00341     }
    00342   else
    00343     {
    00344       dTabds=NULL;
    00345       dTabdp=NULL;
    00346     }
    00347
    00348
    00349   // first mesh point
    00350   if (doDerivs) {
    00351     for (i=1; i<=nvar; ++i)
    00352       {
    00353         for (j=1; j<=nvar; ++j)
    00354           {
    00355             ex->dyds[1][i][j] = (i==j) ? 1 : 0;
    00356             ex->dmds[1][i][j] = (i==j) ? 1 : 0;
    00357           }
    00358         for (j=1; j<=globs->npar; ++j)
    00359           {
    00360             ex->dydp[1][i][j] = 0;
    00361             ex->dmdp[1][i][j] = 0;
    00362           }
    00363       }
    00364   }
    00365   for (j=1; j<=nvar; ++j)
    00366     ex->yComp[1][j]=ex->yTry[1][j];
    00367
    00368
    00369   if (doDerivs && mesh[1]==xMeasure[1])
    00370     {
    00371       observation(globs,ex,mesh[1],ex->yTry[1],ex->yPred[1],ex->par,ex->dmds[1],ex->dmdp[1]);
    00372     }
    00373
    00374   // integrate
    00375   k=1;
    00376   ksav=1;
    00377   if (mesh[1]==xMeasure[1])
    00378     {
    00379       k=2;
    00380       ksav=2;
    00381     }
    00382   for (i=1; i< nPoints; ++i)
    00383     {   // given starting values
    00384       for (j=1; j<=nvar; ++j)
    00385         ys[j]=ex->yTry[i][j];
    00386       if (xMeasure[k]<mesh[i])
    00387         {
    00388           cerr << "fatal: intODE: xMeasure table corrupted (mesh: "
    00389                << mesh[i] << ", xMeasure: " << xMeasure[k] << ")\n";
    00390           throw 1;
    00391         }
    00392       // build table of measuring points
    00393       n=1;
    00394       while (k <= nMeasure && xMeasure[k]<=mesh[i+1])
    00395         xTab[n++]=xMeasure[k++];
    00396       xTab[n]=mesh[i+1];
    00397 #ifdef PRINTMESH
    00398       *dbg << "Mesh point " << mesh[i] << ": xTab = (";
    00399       for (j=1; j<=n; ++j)
    00400         *dbg << '\t' << xTab[j];
    00401       *dbg << ")\n";
    00402       dbg->flush();
    00403 #endif
    00404       //integration
    00405       tabulateValues(globs,ex,mesh[i],mesh[i+1],xTab,n,ys,yTab, dTabds, dTabdp,ex->dyds[i+1],ex->dydp[i+1]);
    00406
    00407 #ifdef PRINTINTEGRATE
    00408       *dbg << "Integration returned\t  (";
    00409       for (j=1; j<=n; ++j)
    00410         {
    00411           for (l=1; l<=nvar; ++l)
    00412             *dbg << " " << yTab[j][l];
    00413           *dbg << ",";
    00414         }
    00415       *dbg << ")\n";
    00416       dbg->flush();
    00417 #endif
    00418       for (j=1; j<=ex->nvar; ++j)
    00419         ex->yComp[i+1][j]=ys[j];   // store value at next mesh point
    00420
    00421       for (j=ksav; j<k; ++j)          // store value at measuring points
    00422         {
    00423           for (l=1; l<=ex->nobs; ++l)
    00424             {
    00425               ex->yPred[j][l]=yTab[j-ksav+1][l];
    00426               if (doDerivs)
    00427                 {
    00428                   for (m=1; m<=ex->nvar; ++m)
    00429                     ex->dmds[j][l][m]=dTabds[j-ksav+1][l][m];
    00430                   for (m=1; m<=globs->npar; ++m)
    00431                     ex->dmdp[j][l][m]=dTabdp[j-ksav+1][l][m];
    00432                 }
    00433             }
    00434         }
    00435       ksav=k;
    00436     }
    00437
    00438       //GnuPlotting
    00439   long mind=2;
    00440   if(doDerivs && (globs->noGnu==FALSE) && doPlot==TRUE)
    00441     {
    00442       for(l=1;l<=ex->nobs;l++)
    00443         {
    00444           if(globs->gnuindex <= globs->ngnu)
    00445             {
    00446               mind=2;
    00447               ofstream gnuout;
    00448               gnuout.open("gnuout.dat");
    00449               ymax=ex->yPred[1][l];
    00450               ymin=ex->yPred[1][l];
    00451               for(j=1;j<=ex->nMeasure;j++)
    00452                 {
    00453                   gnuout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
    00454                   if(mesh[mind] <= ex->xMeasure[j])
    00455                     {
    00456                       observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
    00457                       gnuout << mesh[mind] << "\t" << ys[l] << endl << endl;
    00458                       observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
    00459                       gnuout << mesh[mind] << "\t" << ys[l] << endl;
    00460                       mind++;
    00461                     }
    00462                   //determine y-range
    00463                   if(ymax <= ex->yPred[j][l])
    00464                     ymax=ex->yPred[j][l];
    00465                   if(ymin >= ex->yPred[j][l])
    00466                     ymin=ex->yPred[j][l];
    00467                   if(ymax <= ex->yMeasure[j][l])
    00468                     ymax=ex->yMeasure[j][l];
    00469                   if(ymin >= ex->yMeasure[j][l])
    00470                     ymin=ex->yMeasure[j][l];
    00471                 }
    00472               gnuout.flush();
    00473               gnuout.close();
    00474
    00475               fprintf(globs->gnuFp[globs->gnuindex],"reset\n");
    00476               fprintf(globs->gnuFp[globs->gnuindex],"set yrange[%f:%f]\n",ymin-0.01*ymin,ymax+0.01*ymax);
    00477               fprintf(globs->gnuFp[globs->gnuindex],"set xrange[%f:%f]\n",ex->xMeasure[1],ex->xMeasure[ex->nMeasure]);
    00478               fflush(globs->gnuFp[globs->gnuindex]);
    00479               fprintf(globs->gnuFp[globs->gnuindex],"plot \"gnuout.dat\" title \"\"w l 3,\"gnuout.dat\" u 1:3 title \"\" 1 \n");
    00480               fflush(globs->gnuFp[globs->gnuindex]);
    00481               //DELAY
    00482               for(long dd=1;dd<=1000000;dd++)
    00483                1+1;
    00484             }
    00485           sprintf(name,"Gnuout_Exp%d_Obs%d.dat",expNr,l);
    00486           outout.open(name);
    00487
    00488           globs->gnuindex++;
    00489           mind=2;
    00490           for(j=1;j<=ex->nMeasure;j++)
    00491             {
    00492
    00493               outout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
    00494               if(mesh[mind] <= ex->xMeasure[j])
    00495                 {
    00496                   observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
    00497                   outout << mesh[mind] << "\t" << ys[l] << endl << endl;
    00498                   observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
    00499                   outout << mesh[mind] << "\t" << ys[l] << endl;
    00500                   mind++;
    00501                 }
    00502             }
    00503           outout.flush();
    00504           outout.close();
    00505           for(long dd=1;dd<=1000000;dd++)
    00506             1+1;
    00507         }
    00508     }
    00509
    00510   if (doDerivs) {
    00511     free_d3tensor(dTabdp,1,tableLen,1,nvar,1,globs->npar);
    00512     free_d3tensor(dTabds,1,tableLen,1,nvar,1,nvar);
    00513   }
    00514   free_dmatrix(yTab,1,tableLen,1,nvar);
    00515   free_dvector(xTab,1,tableLen);
    00516   free_dvector(ys,1,nvar);
    00517 }

double max | ( | double | a,
|  | double | b
| ) | `[inline]`

| Definition at line 62 of file intODE.cc. Referenced by maxTableLen(), and tabulateValues().
    00063 {
    00064   if(a>=b)
    00065     return(a);
    00066   else
    00067     return(b);
    00068 }

long maxTableLen | ( | double * | mesh,
|  | long | nPoints,
|  | double * | xMeasure,
|  | long | nMeasure
| ) |

| Definition at line 35 of file intODE.cc. References max(). Referenced by intODE().
    00036 {
    00037     int i,j, max=0;
    00038
    00039     return(nPoints > nMeasure ? nPoints : nMeasure);  // to be improved
    00040 }

void sensDerivs | ( | double | t,
|  | double * | Yt,
|  | double * | Ft,
|  | double * | p
| ) |

| Definition at line 72 of file intODE.cc. References dfill0(), and dmatrix(). Referenced by f(), and tabulateValues().
    00073 {
    00074   //  sensitivity equations
    00075   // P=(p,y0)
    00076   // S=dy/dP
    00077   // Y=(y,S)
    00078   // F=dY/dt=(f,df/dy*S+df/dP)
    00079   // Y and F are stored in column major format and rolled out as Yt and Ft
    00080
    00081   long i, j, k;
    00082   dfill0 (Ft, (1 + NPARAMS + NEQNS) * NEQNS);       // initialize to zero
    00083   // ode, inhomo and jacobi set only non-zero elements
    00084   // inhome expects a dmatrix; we create one who's data lies in Ft
    00085   double **dfdpt = new double *[NPARAMS] - 1;
    00086   for (i = 1; i <= NPARAMS; i++)
    00087     dfdpt[i] = Ft + i * NEQNS;
    00088   ode (t, Yt, Ft, p);           // first part of Y,F are y,f
    00089   inhomo (t, Yt, dfdpt, p);     // fill df/dp into F
    00090   delete[]++ dfdpt;
    00091
    00092   // now F=df/dP
    00093   // add M*S, M=df/dy
    00094   double **M = dmatrix (1, NEQNS, 1, NEQNS);
    00095   dfill0 (M[1], NEQNS * NEQNS);
    00096   jacobi (t, Yt, M, p);         // get M
    00097   for (i = 1; i <= NEQNS; i++)
    00098     for (j = 1; j <= NPARAMS + NEQNS; j++)
    00099       {
    00100         double &sum = Ft[j * NEQNS + i];
    00101         for (k = 1; k <= NEQNS; k++)
    00102           sum += M[k][i] * Yt[j * NEQNS + k];
    00103       }                         // for i,j
    00104   free_dmatrix (M, 1, NEQNS, 1, NEQNS);
    00105  }                               // sensDerivs

void tabulateValues | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | double | t0,
|  | double | t1,
|  | double * | t,
|  | long | n,
|  | double * | state,
|  | double ** | y,
|  | double *** | dmds,
|  | double *** | dmdp,
|  | double ** | dyds,
|  | double ** | dydp
| ) |

| Definition at line 108 of file intODE.cc. References call_odessa(), dfill0(), dmatrix(), Glob::doP, dvector(), Glob::eps, initYt(), integrateRK(), Glob::integrator, max(), Glob::maxstp, Glob::npar, GlobExp::nvar, GlobExp::par, Glob::rkqs_ign, sensDerivs(), Glob::stiff, and TRUE.
    00110 {
    00111   // integrate from t0 to t1, storing values at t[1..n]
    00112   // t[n] may be ==t1 or <t1
    00113   // in the latter case we will integrate from t[n] to t1 in the
    00114   // final loop k=n+1 without storing any data points
    00115   // If dmds!=NULL, compute derivatives as well
    00116   // and store them in dmds,dmdp,dyds,dydp
    00117   // with respect to observation, see odin.doc: Observation function
    00118
    00119   //BEGIN von Felix veraendert
    00120   //   double hmin=max(fabs(t0),fabs(t1))*1e-8;
    00121   double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
    00122   // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest.
    00123   // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern
    00124   // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion.
    00125   // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht
    00126   // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich
    00127   // von 1e-18 noch einigermassen weit entfernt bleibt.
    00128   // Man koennte das in Zukunft noch in einen externen Parameter stecken.
    00129   //END von Felix veraendert
    00130
    00131   long nPExp=globs->npar;
    00132   long nVar=ex->nvar;
    00133   long nobs=NOBS;
    00134   int generic;
    00135   double *p=ex->par;
    00136   double eps=globs->eps;
    00137   double h1 = (t1 - t0) * .01;
    00138   long nok,nbad;
    00139
    00140   initInt(globs,ex);
    00141   if (n > 0 && (t0 > t[1] || t1 < t[n]))
    00142     {
    00143       cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl;
    00144       throw 1;
    00145     }
    00146
    00147   if (globs->integrator==3)
    00148     {                           // said -odessa
    00149 #ifdef ODESSA
    00150       double &rtol = eps;
    00151 //BEGIN von Felix veraendert
    00152 //       double atol=eps*1e-6;
    00153       double atol = eps;
    00154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!)
    00155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu
    00156 // siehe odessa.f: ATOL und RTOL
    00157 //END von Felix veraendert
    00158       // odeint uses TINY=1e-30 in a similar way as odessa's atol,
    00159       // but atol<1e-16 or so results in step size underflow
    00160
    00161       int odeint_MAXSTP=globs->maxstp;
    00162       int d_flag=1;
    00163       int stif=globs->stiff;
    00164       char doPStr[globs->npar+1];
    00165
    00166       for(long i=1;i<=globs->npar;i++)
    00167         doPStr[i]=(char)globs->doP[i];
    00168
    00169       call_odessa(globs,ex,nVar, nPExp, doPStr,
    00170          d_flag, odeint_MAXSTP, stif,
    00171          TRUE, n, t, t0, t1, state, y,
    00172          dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1);
    00173
    00174       return;
    00175 #else
    00176       cerr << "-odessa used but not linked" << endl;
    00177       throw 1;
    00178 #endif // ODESSA
    00179     }                           // if odessa
    00180       double *gy=dvector(1,nobs);
    00181       long i, j, k;
    00182       double xx1 = t0, xx2;
    00183
    00184       if (dmds)
    00185         {                           // compute sensitivities
    00186           // create extended state vector Y
    00187           double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
    00188           initYt (Yt[0], state);
    00189
    00190           // chain rule for extra observables
    00191           double **dgdy, **dgdp;
    00192
    00193           dgdy = dmatrix (1, nobs, 1, nVar);
    00194           dgdp = dmatrix (1, nobs, 1, nPExp);
    00195           dfill0 (dgdy[1], nobs * nVar);
    00196           dfill0 (dgdp[1], nobs * nPExp);
    00197
    00198
    00199           for (k = 1; k <= n + 1; k++)
    00200             {
    00201               xx2 = (k > n) ? t1 :  // end point
    00202                 t[k];               // data point
    00203
    00204               globs->rkqs_ign = (nPExp + nVar) * nVar;
    00205               if (xx2 > xx1 && nVar > 0)
    00206                 {
    00207                   if(globs->integrator==1)
    00208                     {
    00209                       //Runge-Kutta integation
    00210                       integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,
    00211                                   eps, h1, hmin, &nok, &nbad, sensDerivs);
    00212                     }
    00213                   else if(globs->integrator==2)
    00214                     {
    00215                       cerr << "CVODES not optimized skip integrator" << endl;
    00216                       throw 1;
    00217
    00218                       //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE);
    00219                     }
    00220                 }
    00221               // write state to y and sensitivities to dm/dp and dm/ds
    00222               if (k <= n)
    00223                 {
    00224                   dcopy (y[k], Yt[0], nVar);
    00225                   generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp);
    00226                   for(i=1;i<=nobs;i++)
    00227                     y[k][i]=gy[i];
    00228
    00229                   for (i = 1; i <= nobs; i++)
    00230                     {
    00231                       for (j = 1; j <= nPExp; j++)
    00232                         {
    00233                           double &dest = dmdp[k][i][j];
    00234                           if (!generic)
    00235                             {
    00236                               dest = dgdp[i][j];
    00237                               for (long l = 1; l <= nVar; l++)
    00238                                 dest += dgdy[i][l] * Yt[j][l];
    00239                             }
    00240                           else
    00241                             dest = Yt[j][i];
    00242                         }           // for j
    00243                       for (j = 1; j <= nVar; j++)
    00244                         {
    00245                           double &dest = dmds[k][i][j];
    00246                           if (!generic)
    00247                             {
    00248                               dest = 0;
    00249                               for (long l = 1; l <= nVar; l++)
    00250                                 dest += dgdy[i][l] * Yt[nPExp + j][l];
    00251                             }
    00252                           else
    00253                             dest = Yt[nPExp + j][i];
    00254                         }           // for j
    00255                     }               // for i
    00256                 }                   // if k<=n
    00257               xx1 = xx2;
    00258             }                       // for k
    00259
    00260           // copy state and its derivatives to state, dydp and dyds
    00261           dcopy (state, Yt[0], nVar);       // original state vector
    00262           // write sensitivities to dy/dp and dy/ds
    00263           for (i = 1; i <= nVar; i++)
    00264             {
    00265               for (j = 1; j <= nPExp; j++)
    00266                 dydp[i][j] = Yt[j][i];
    00267               for (j = 1; j <= nVar; j++)
    00268                 dyds[i][j] = Yt[nPExp + j][i];
    00269             }
    00270
    00271           free_dmatrix (dgdp,1 ,nobs , 1, nPExp);
    00272           free_dmatrix (dgdy,1 ,nobs, 1, nVar);
    00273           free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar);
    00274         }
    00275       else
    00276         {                           // !dmds
    00277           for (k = 1; k <= n + 1; k++)
    00278             {
    00279               xx2 = (k <= n) ? t[k] : t1;
    00280
    00281               globs->rkqs_ign = 0;
    00282               if (xx2 > xx1 && nVar > 0)
    00283                 {
    00284                   if(globs->integrator==1)
    00285                     {
    00286                       //Runge-Kutta integration
    00287                       integrateRK(globs,state, p, nVar, xx1, xx2,
    00288                                   eps, h1, hmin, &nok, &nbad, ode);
    00289                     }
    00290                   else if(globs->integrator==2)
    00291                     {
    00292                       cerr << "CVODES not optimized skip integrator" << endl;
    00293                       throw 1;
    00294
    00295                       //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE);
    00296                     }
    00297                 }
    00298               if (k <= n)
    00299                 {
    00300                   dcopy (y[k], state, nVar);
    00301                   generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL);
    00302                   for(i=1;i<=nobs;i++)
    00303                     y[k][i]=gy[i];
    00304                 }
    00305               xx1 = xx2;
    00306             }                       // for k
    00307         }
    00308       // if dmds else
    00309       free_dvector(gy,1,nobs);
    00310 }