00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<string.h>
00005 #include<stdio.h>
00006
00007 #include "../def.h"
00008 #include "../model.h"
00009 #include "../nr.h"
00010
00011 using namespace std;
00012
00013
00014
00015
00016
00017
00018 void integrateRK(Glob *globs,double ystart[], double p[], long nvar, double x1, double x2,double eps,
00019 double h1, double hmin, long *nok, long *nbad,
00020 void (*derivType) (double, double[], double[], double[]));
00021
00022
00023
00024 void call_odessa(Glob *globs,GlobExp *ex,int N, int M_PAR_GES, char *doP,
00025 int D_FLAG, int MAX_NR_STEP, int STIFF, int inhomogen, int nmesh,
00026 double *tmesh, double t0, double t1, double *zustand, double **y,
00027 double ***dmds, double ***dmdp, double **dyds, double **dydp,
00028 double *parameter, double rtol, double atol,
00029 double min_stepsize, double max_stepsize, double initial_stepsize);
00030
00031
00032
00033
00034
00035 long maxTableLen(double *mesh, long nPoints, double *xMeasure, long nMeasure)
00036 {
00037 int i,j, max=0;
00038
00039 return(nPoints > nMeasure ? nPoints : nMeasure);
00040 }
00041
00042 inline void dcopy (double * dest, double * src, long n)
00043 {
00044 memcpy (dest + 1, src + 1, n * sizeof (double));
00045 }
00046
00047 inline void dfill0 (double * x, long n)
00048 {
00049 memset (x + 1, 0, n * sizeof (double));
00050 }
00051
00052 void initYt (double *Yt, double *state)
00053 {
00054
00055 if (state)
00056 dcopy (Yt, state, NEQNS);
00057 dfill0 (Yt + NEQNS, (NPARAMS + NEQNS) * NEQNS);
00058 for (long i = 1; i <= NEQNS; i++)
00059 Yt[(NPARAMS + i) * NEQNS + i] = 1;
00060 }
00061
00062 inline double max(double a,double b)
00063 {
00064 if(a>=b)
00065 return(a);
00066 else
00067 return(b);
00068 }
00069
00070
00071
00072 void sensDerivs (double t, double *Yt, double *Ft, double *p)
00073 {
00074
00075
00076
00077
00078
00079
00080
00081 long i, j, k;
00082 dfill0 (Ft, (1 + NPARAMS + NEQNS) * NEQNS);
00083
00084
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);
00089 inhomo (t, Yt, dfdpt, p);
00090 delete[]++ dfdpt;
00091
00092
00093
00094 double **M = dmatrix (1, NEQNS, 1, NEQNS);
00095 dfill0 (M[1], NEQNS * NEQNS);
00096 jacobi (t, Yt, M, p);
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 }
00104 free_dmatrix (M, 1, NEQNS, 1, NEQNS);
00105 }
00106
00107
00108 void tabulateValues (Glob *globs,GlobExp *ex,double t0, double t1, double *t, long n, double *state,
00109 double **y, double ***dmds, double ***dmdp, double **dyds,double **dydp)
00110 {
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
00122
00123
00124
00125
00126
00127
00128
00129
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 {
00149 #ifdef ODESSA
00150 double &rtol = eps;
00151
00152
00153 double atol = eps;
00154
00155
00156
00157
00158
00159
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 }
00180 double *gy=dvector(1,nobs);
00181 long i, j, k;
00182 double xx1 = t0, xx2;
00183
00184 if (dmds)
00185 {
00186
00187 double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
00188 initYt (Yt[0], state);
00189
00190
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 :
00202 t[k];
00203
00204 globs->rkqs_ign = (nPExp + nVar) * nVar;
00205 if (xx2 > xx1 && nVar > 0)
00206 {
00207 if(globs->integrator==1)
00208 {
00209
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
00219 }
00220 }
00221
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 }
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 }
00255 }
00256 }
00257 xx1 = xx2;
00258 }
00259
00260
00261 dcopy (state, Yt[0], nVar);
00262
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 {
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
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
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 }
00307 }
00308
00309 free_dvector(gy,1,nobs);
00310 }
00311
00312
00313 void intODE(GlobExp *ex,Glob *globs,int doDerivs,int doPlot,long expNr)
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
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
00333
00334
00335
00336
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
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
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 {
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
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
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];
00420
00421 for (j=ksav; j<k; ++j)
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
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
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
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 }
00518