- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/modules/simInit.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 | 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 | simInit (GlobExp *ex, Glob *globs)

## Function Documentation

void simInit | ( | GlobExp * | ex,
|  | Glob * | globs
| ) |

| Definition at line 21 of file simInit.cc. References dmatrix(), dvector(), GlobExp::nPoints, GlobExp::nvar, GlobExp::y0, and GlobExp::yTry.
    00022 {
    00023   long i,j;
    00024   long nvar=ex->nvar;
    00025   long nPoints=ex->nPoints;
    00026   long idum=-42;
    00027   double *t=dvector(1,2);
    00028   double **y=dmatrix(1,2,1,nvar);
    00029
    00030    if (ex->y0)
    00031     {
    00032       for (i=1; i<=nvar; ++i)
    00033         ex->yTry[1][i]=ex->y0[i];
    00034     }
    00035    else
    00036      {
    00037         for (i=1; i<=nvar; ++i)
    00038           {
    00039             ex->yTry[1][i]=DefYValues[i-1];
    00040           }
    00041      }
    00042
    00043    //integrating
    00044
    00045    for(i=2;i<=nPoints;i++)
    00046      {
    00047        for(j=1;j<=nvar;j++)
    00048          ex->yTry[i][j]=ex->yTry[i-1][j];
    00049        t[1]=ex->mesh[i-1];
    00050        t[2]=ex->mesh[i];
    00051        tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
    00052       for(j=1;j<=nvar;j++)
    00053         ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
    00054      }
    00055
    00056
    00057
    00058    free_dvector(t,1,2);
    00059    free_dmatrix(y,1,2,1,nvar);
    00060 }

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 }