/home/peifer/diffit/modules/call_odessa.cc

Go to the documentation of this file.
00001 /****************************************************************/
00002 /*                                                              */
00003 /* Programm zum Aufruf von ODESSA                               */
00004 /* Aufruf-Syntax ist diesselbe wie bei CALL_DOP853              */
00005 /* Ziel: Sensitivitätsanalyse der Differentialgleichung ODE     */
00006 /* bereitgestellte Routinen:                                    */
00007 /*                          ODE                                 */
00008 /*                          AD_ODE, OD_ODE                      */
00009 /*                          JAC                                 */
00010 /* Thorsten Mueller, Apr, 21, 1999                              */
00011 /* und September, 18, 2001 (Krakau, Polen)                      */
00012 /*                                                              */
00013 /****************************************************************/
00014 /* modified for inclusion in odin by w.h. 11/01
00015  */
00016 
00017 #include <math.h>
00018 #include <fstream>
00019 #include <iomanip>
00020 #include <iostream>
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024 
00025 #include "../def.h"
00026 #include "../model.h"
00027 #include "../nr.h"
00028 
00029 using namespace std;
00030 
00031 
00032 inline void _dcopy (double * dest, double * src, long n)
00033 {
00034   memcpy (dest + 1, src + 1, n * sizeof (double));
00035 }
00036 
00037 inline void _dfill0 (double * x, long n)
00038 {
00039   memset (x + 1, 0, n * sizeof (double));
00040 }
00041 
00042 void _initYt (double *Yt, double *state)
00043 {
00044   // initialise Y
00045   if (state)
00046     _dcopy (Yt, state, NEQNS);    // original state vector
00047   _dfill0 (Yt + NEQNS, (NPARAMS + NEQNS) * NEQNS);    // 0 in middle and rear part
00048   for (long i = 1; i <= NEQNS; i++)       // identity matrix in rear part dy/dy0
00049     Yt[(NPARAMS + i) * NEQNS + i] = 1;
00050 }                       
00051 
00052 
00053 extern "C" int odessa_ (int NEQ[],
00054                         double Y[],
00055                         double PAR[],
00056                         double *T,
00057                         double *TOUT,
00058                         int *ITOL,
00059                         double *RTOL,
00060                         double *ATOL,
00061                         int *ITASK,
00062                         int *ISTATE,
00063                         int IOPT[],
00064                         double RWORK[],
00065                         int *LRW, int IWORK[], int *LIW, int *MF);
00066 
00067 extern "C" void dop_komponente_ (int NEQ[], int &J_OPT, int &J_GES);
00068 
00069 void call_odessa(Glob *globs,GlobExp *ex,int N, int M_PAR_GES, char *doP,
00070    int D_FLAG, int MAX_NR_STEP, int STIFF, int inhomogen, int nmesh,
00071    double *tmesh, double t0, double t1, double *zustand, double **y,
00072    double ***dmds, double ***dmdp, double **dyds, double **dydp,
00073    double *parameter, double rtol, double atol,
00074    double min_stepsize, double max_stepsize, double initial_stepsize)
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
00500 
00501 
00502 /* wrapper routines between odessa and the dynamical equations generated
00503  * from the .mdf
00504  */
00505 
00506 
00507 
00508 extern "C" void ode_ (int &n, int &m, double &x, double *y, double *dydx, double *p)
00509 {
00510   ode (x, y - 1, dydx - 1, p - 1);
00511 }                               // ode_
00512 
00513 extern "C" void jac_ (int *NEQ, double &t, double *y, int &n, double *p, int &M, double *pd)
00514 {
00515   double **dfdy = new double *[NEQNS] - 1;
00516   for (int i = 1; i <= NEQNS; i++)
00517     dfdy[i] = pd - 1 + (i - 1) * NEQNS;
00518   jacobi (t, y - 1, dfdy, p - 1);
00519   delete[]++ dfdy;
00520 }
00521 
00522 extern "C" void od_ode_ (int *NEQ, double &t, double *y, int &N, double *p,
00523          int &M, double *dfdp, int &j)
00524 {
00525   if (!j)
00526     return;
00527   inhomo1 (t, y - 1, dfdp - 1, p - 1, j);
00528 }                               // od_ode_
00529 
00530 extern "C" void ad_ode_ (int &, int &, int &, double &, double *, double *, double *,
00531          int &, double *, double *, int &)
00532 {
00533   cerr << "unexpected call to AD_ODE" <<endl;
00534 }

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