- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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

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

Go to the source code of this file.

## Functions
void | _dcopy (double *dest, double *src, long n)
void | _dfill0 (double *x, long n)
void | _initYt (double *Yt, double *state)
int | odessa_ (int NEQ[], double Y[], double PAR[], double *T, double *TOUT, int *ITOL, double *RTOL, double *ATOL, int *ITASK, int *ISTATE, int IOPT[], double RWORK[], int *LRW, int IWORK[], int *LIW, int *MF)
void | dop_komponente_ (int NEQ[], int &J_OPT, int &J_GES)
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)
void | ode_ (int &n, int &m, double &x, double *y, double *dydx, double *p)
void | jac_ (int *NEQ, double &t, double *y, int &n, double *p, int &M, double *pd)
void | od_ode_ (int *NEQ, double &t, double *y, int &N, double *p, int &M, double *dfdp, int &j)
void | ad_ode_ (int &, int &, int &, double &, double *, double *, double *, int &, double *, double *, int &)

## Function Documentation

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

| Definition at line 32 of file call_odessa.cc. Referenced by _initYt().
    00033 {
    00034   memcpy (dest + 1, src + 1, n * sizeof (double));
    00035 }

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

| Definition at line 37 of file call_odessa.cc. Referenced by _initYt().
    00038 {
    00039   memset (x + 1, 0, n * sizeof (double));
    00040 }

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

| Definition at line 42 of file call_odessa.cc. References _dcopy(), and _dfill0().
    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 }

void ad_ode_ | ( | int & | ,
|  | int & | ,
|  | int & | ,
|  | double & | ,
|  | double * | ,
|  | double * | ,
|  | double * | ,
|  | int & | ,
|  | double * | ,
|  | double * | ,
|  | int & |
| ) |

| Definition at line 530 of file call_odessa.cc.
    00532 {
    00533   cerr << "unexpected call to AD_ODE" <<endl;
    00534 }

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 dop_komponente_ | ( | int | NEQ[],
|  | int & | J_OPT,
|  | int & | J_GES
| ) |

void jac_ | ( | int * | NEQ,
|  | double & | t,
|  | double * | y,
|  | int & | n,
|  | double * | p,
|  | int & | M,
|  | double * | pd
| ) |

| Definition at line 513 of file call_odessa.cc.
    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 }

void od_ode_ | ( | int * | NEQ,
|  | double & | t,
|  | double * | y,
|  | int & | N,
|  | double * | p,
|  | int & | M,
|  | double * | dfdp,
|  | int & | j
| ) |

| Definition at line 522 of file call_odessa.cc.
    00524 {
    00525   if (!j)
    00526     return;
    00527   inhomo1 (t, y - 1, dfdp - 1, p - 1, j);
    00528 }                               // od_ode_

void ode_ | ( | int & | n,
|  | int & | m,
|  | double & | x,
|  | double * | y,
|  | double * | dydx,
|  | double * | p
| ) |

| Definition at line 508 of file call_odessa.cc.
    00509 {
    00510   ode (x, y - 1, dydx - 1, p - 1);
    00511 }                               // ode_

int odessa_ | ( | int | NEQ[],
|  | double | Y[],
|  | double | PAR[],
|  | double * | T,
|  | double * | TOUT,
|  | int * | ITOL,
|  | double * | RTOL,
|  | double * | ATOL,
|  | int * | ITASK,
|  | int * | ISTATE,
|  | int | IOPT[],
|  | double | RWORK[],
|  | int * | LRW,
|  | int | IWORK[],
|  | int * | LIW,
|  | int * | MF
| ) |