#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 | 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 | simInit (GlobExp *ex, Glob *globs) |
|
Definition at line 21 of file simInit.cc. References dmatrix(), dvector(), GlobExp::nPoints, GlobExp::nvar, GlobExp::y0, and GlobExp::yTry. 00022 { 00023 long i,j; 00024 long nvar=ex->nvar; 00025 long nPoints=ex->nPoints; 00026 long idum=-42; 00027 double *t=dvector(1,2); 00028 double **y=dmatrix(1,2,1,nvar); 00029 00030 if (ex->y0) 00031 { 00032 for (i=1; i<=nvar; ++i) 00033 ex->yTry[1][i]=ex->y0[i]; 00034 } 00035 else 00036 { 00037 for (i=1; i<=nvar; ++i) 00038 { 00039 ex->yTry[1][i]=DefYValues[i-1]; 00040 } 00041 } 00042 00043 //integrating 00044 00045 for(i=2;i<=nPoints;i++) 00046 { 00047 for(j=1;j<=nvar;j++) 00048 ex->yTry[i][j]=ex->yTry[i-1][j]; 00049 t[1]=ex->mesh[i-1]; 00050 t[2]=ex->mesh[i]; 00051 tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL); 00052 for(j=1;j<=nvar;j++) 00053 ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum); 00054 } 00055 00056 00057 00058 free_dvector(t,1,2); 00059 free_dmatrix(y,1,2,1,nvar); 00060 }
|
|
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 }
|