#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 &) |
|
Definition at line 32 of file call_odessa.cc. Referenced by _initYt().
|
|
Definition at line 37 of file call_odessa.cc. Referenced by _initYt().
|
|
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 }
|
|
Definition at line 530 of file call_odessa.cc.
|
|
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
|
|
|
|
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 }
|
|
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_
|
|
Definition at line 508 of file call_odessa.cc.
|
|
|