/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
 


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