/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 

Generated on Mon Jan 29 17:09:12 2007 for Diffit by  doxygen 1.4.6