#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "../def.h"
#include "../model.h"
#include "../nr.h"
Go to the source code of this file.
Functions | |
void | integrateRK (Glob *globs, double ystart[], double p[], long nvar, double x1, double x2, double eps, double h1, double hmin, long *nok, long *nbad, void(*derivType)(double, double[], double[], double[])) |
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) |
long | maxTableLen (double *mesh, long nPoints, double *xMeasure, long nMeasure) |
void | dcopy (double *dest, double *src, long n) |
void | dfill0 (double *x, long n) |
void | initYt (double *Yt, double *state) |
double | max (double a, double b) |
void | sensDerivs (double t, double *Yt, double *Ft, double *p) |
void | tabulateValues (Glob *globs, GlobExp *ex, double t0, double t1, double *t, long n, double *state, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp) |
void | intODE (GlobExp *ex, Glob *globs, int doDerivs, int doPlot, long expNr) |
|
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 42 of file intODE.cc. Referenced by initYt().
|
|
Definition at line 47 of file intODE.cc. Referenced by initYt(), sensDerivs(), and tabulateValues().
|
|
Definition at line 52 of file intODE.cc. References dcopy(), and dfill0(). Referenced by tabulateValues(). 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 }
|
|
Definition at line 19 of file integrateRK.cc. References dvector(), Glob::maxstp, nrerror(), and SIGN. Referenced by tabulateValues(). 00022 { 00023 long ode_fkmax = 0, ode_fkount = 0, ode_kmax = 0, ode_kount = 0; 00024 double ode_dxsav = 0, *ode_fxp = 0, *ode_xp = 0, **ode_yp = 0; 00025 long rkqs_ign = 0; 00026 long nstp,i; 00027 double xsav,x,hnext,hdid,h; 00028 double *yscal,*y,*dydx; 00029 long odeint_MAXSTP=globs->maxstp; 00030 00031 yscal=dvector(1,nvar); 00032 y=dvector(1,nvar); 00033 dydx=dvector(1,nvar); 00034 x=x1; 00035 h=SIGN(h1,x2-x1); 00036 *nok = (*nbad) = ode_kount = 0; 00037 for (i=1;i<=nvar;i++) y[i]=ystart[i]; 00038 if (ode_dxsav<0) nrerror("odeint:dxsav<0"); 00039 if (ode_kmax>0 && ode_kount<ode_kmax-1) { 00040 ode_xp[++ode_kount]=x; 00041 for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i]; 00042 xsav=x; 00043 } /* save first point, w.h. 19.03.99 */ 00044 for (nstp=1;nstp<=odeint_MAXSTP;nstp++) { 00045 (*derivs)(x,y,dydx,p); 00046 for (i=1;i<=nvar;i++) 00047 yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY; 00048 if (ode_fkmax) { 00049 if (ode_fkount>=ode_fkmax) nrerror("ode_fkount exceeds ode_fkmax"); 00050 ode_fxp[++ode_fkount]=x; 00051 } 00052 if (ode_kmax>0 && ode_kount<ode_kmax-1 && fabs(x-xsav)>ode_dxsav) { 00053 ode_xp[++ode_kount]=x; 00054 for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i]; 00055 xsav=x; 00056 } 00057 if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x; 00058 rkqs(y,dydx,p,nvar,&x,h,eps,yscal,&hdid,&hnext,globs->rkqs_ign,derivs); 00059 if (hdid == h) ++(*nok); 00060 else ++(*nbad); 00061 if ((x-x2)*(x2-x1) >= 0.0) { 00062 for (i=1;i<=nvar;i++) ystart[i]=y[i]; 00063 if (ode_kmax>0 && ode_kount<ode_kmax-1) { 00064 ode_xp[++ode_kount]=x; 00065 for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i]; 00066 } 00067 free_dvector(dydx,1,nvar); 00068 free_dvector(y,1,nvar); 00069 free_dvector(yscal,1,nvar); 00070 return; 00071 } 00072 if (fabs(hnext) <= hmin) 00073 { 00074 cerr << "Step size too small in odeint\n"; 00075 throw 1; 00076 } 00077 h=hnext; 00078 } 00079 cerr << "Too many steps in routine odeint; try -maxstp <n> greater n\n"; 00080 throw 1; 00081 }
|
|
Definition at line 313 of file intODE.cc. References d3tensor(), dmatrix(), GlobExp::dmdp, GlobExp::dmds, dvector(), GlobExp::dydp, GlobExp::dyds, maxTableLen(), GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, GlobExp::nvar, and GlobExp::xMeasure. 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 }
|
|
Definition at line 62 of file intODE.cc. Referenced by maxTableLen(), and tabulateValues().
|
|
Definition at line 35 of file intODE.cc. References max(). Referenced by intODE(). 00036 { 00037 int i,j, max=0; 00038 00039 return(nPoints > nMeasure ? nPoints : nMeasure); // to be improved 00040 }
|
|
Definition at line 72 of file intODE.cc. References dfill0(), and dmatrix(). Referenced by f(), and tabulateValues(). 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
|
|
Definition at line 108 of file intODE.cc. References call_odessa(), dfill0(), dmatrix(), Glob::doP, dvector(), Glob::eps, initYt(), integrateRK(), Glob::integrator, max(), Glob::maxstp, Glob::npar, GlobExp::nvar, GlobExp::par, Glob::rkqs_ign, sensDerivs(), Glob::stiff, and TRUE. 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 }
|