- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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

Go to the documentation of this file.

    00001 #include<iostream>
    00002 #include<fstream>
    00003 #include<math.h>
    00004 #include<string.h>
    00005 #include<stdio.h>
    00006
    00007 #include "../def.h"
    00008 #include "../model.h"
    00009 #include "../nr.h"
    00010
    00011 using namespace std;
    00012
    00013
    00014 //Begin definition of sub-module prototypes
    00015
    00016 //Module modules/setInitialValues.cc
    00017
    00018 void integrateRK(Glob *globs,double ystart[], double p[], long nvar, double x1, double x2,double eps,
    00019                  double h1, double hmin, long *nok, long *nbad,
    00020                  void (*derivType) (double, double[], double[], double[]));
    00021
    00022 //void integrateCVODES(Glob *globs,double ystart[], double p[], long nvar, double x1,  double x2,double eps,int sensi);
    00023
    00024 void call_odessa(Glob *globs,GlobExp *ex,int N, int M_PAR_GES, char *doP,
    00025    int D_FLAG, int MAX_NR_STEP, int STIFF, int inhomogen, int nmesh,
    00026    double *tmesh, double t0, double t1, double *zustand, double **y,
    00027    double ***dmds, double ***dmdp, double **dyds, double **dydp,
    00028    double *parameter, double rtol, double atol,
    00029                  double min_stepsize, double max_stepsize, double initial_stepsize);
    00030 //End definition of sub-module prototypes
    00031
    00032
    00033 /* compute maximum number of measuring points between two mesh points */
    00034
    00035 long maxTableLen(double *mesh, long nPoints, double *xMeasure, long nMeasure)
    00036 {
    00037     int i,j, max=0;
    00038
    00039     return(nPoints > nMeasure ? nPoints : nMeasure);  // to be improved
    00040 }
    00041
    00042 inline void dcopy (double * dest, double * src, long n)
    00043 {
    00044   memcpy (dest + 1, src + 1, n * sizeof (double));
    00045 }
    00046
    00047 inline void dfill0 (double * x, long n)
    00048 {
    00049   memset (x + 1, 0, n * sizeof (double));
    00050 }
    00051
    00052 void initYt (double *Yt, double *state)
    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 }
    00061
    00062 inline double max(double a,double b)
    00063 {
    00064   if(a>=b)
    00065     return(a);
    00066   else
    00067     return(b);
    00068 }
    00069
    00070
    00071
    00072 void sensDerivs (double t, double *Yt, double *Ft, double *p)
    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
    00106
    00107
    00108 void tabulateValues (Glob *globs,GlobExp *ex,double t0, double t1, double *t, long n, double *state,
    00109                      double **y, double ***dmds, double ***dmdp, double **dyds,double **dydp)
    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 }
    00311
    00312
    00313 void intODE(GlobExp *ex,Glob *globs,int doDerivs,int doPlot,long expNr)
    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 }
    00518