00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <math.h>
00018 #include <fstream>
00019 #include <iomanip>
00020 #include <iostream>
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024
00025 #include "../def.h"
00026 #include "../model.h"
00027 #include "../nr.h"
00028
00029 using namespace std;
00030
00031
00032 inline void _dcopy (double * dest, double * src, long n)
00033 {
00034 memcpy (dest + 1, src + 1, n * sizeof (double));
00035 }
00036
00037 inline void _dfill0 (double * x, long n)
00038 {
00039 memset (x + 1, 0, n * sizeof (double));
00040 }
00041
00042 void _initYt (double *Yt, double *state)
00043 {
00044
00045 if (state)
00046 _dcopy (Yt, state, NEQNS);
00047 _dfill0 (Yt + NEQNS, (NPARAMS + NEQNS) * NEQNS);
00048 for (long i = 1; i <= NEQNS; i++)
00049 Yt[(NPARAMS + i) * NEQNS + i] = 1;
00050 }
00051
00052
00053 extern "C" int odessa_ (int NEQ[],
00054 double Y[],
00055 double PAR[],
00056 double *T,
00057 double *TOUT,
00058 int *ITOL,
00059 double *RTOL,
00060 double *ATOL,
00061 int *ITASK,
00062 int *ISTATE,
00063 int IOPT[],
00064 double RWORK[],
00065 int *LRW, int IWORK[], int *LIW, int *MF);
00066
00067 extern "C" void dop_komponente_ (int NEQ[], int &J_OPT, int &J_GES);
00068
00069 void call_odessa(Glob *globs,GlobExp *ex,int N, int M_PAR_GES, char *doP,
00070 int D_FLAG, int MAX_NR_STEP, int STIFF, int inhomogen, int nmesh,
00071 double *tmesh, double t0, double t1, double *zustand, double **y,
00072 double ***dmds, double ***dmdp, double **dyds, double **dydp,
00073 double *parameter, double rtol, double atol,
00074 double min_stepsize, double max_stepsize, double initial_stepsize)
00075 {
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
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
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
00133
00134
00135 int IPAR_offset = 10;
00136
00137
00138
00139 int *NEQ;
00140 NEQ = ivector (1, IPAR_offset + nPExp);
00141
00142
00143
00144 for (i = 1; i <= nPExp; i++)
00145 NEQ[IPAR_offset + i] = (doP[i] != '0' ? 1 : 0);
00146
00147
00148
00149
00150
00151 NEQ[1] = nVar;
00152 NEQ[2] = M_GES;
00153 NEQ[3] = M_PAR_SENS;
00154 NEQ[4] = D_FLAG;
00155 NEQ[5] = nPExp;
00156
00157
00158
00159
00160
00161 double **Yt = dmatrix (0, M_GES, 1, nVar);
00162
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
00170
00171 if (j <= M_PAR_SENS)
00172 for (i = 1; i <= nVar; i++)
00173 Yt[j][i] = 0.0;
00174
00175
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 }
00184
00185
00186
00187
00188
00189
00190 double *PARAMETER;
00191 PARAMETER = dvector (1, nPExp + nVar * nVar);
00192
00193
00194
00195 int zaehler = 1;
00196 for (i = 1; i <= nPExp; i++)
00197 PARAMETER[i] = parameter[i];
00198
00199
00200
00201
00202
00203 int SqualCheck=FALSE;
00204 int ITOL = SqualCheck ? 1 : 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;
00214 ATOL[i] = atol;
00215 }
00216 for (; i <= nTOL; i++)
00217 {
00218 RTOL[i] = rtol * 1e8;
00219 ATOL[i] = atol * 1e8;
00220 }
00221 }
00222 else
00223 {
00224 RTOL[1] = rtol;
00225 ATOL[1] = atol;
00226 }
00227
00228
00229
00230
00231
00232 int ITASK = 1;
00233
00234 int ISTATE = 1;
00235
00236 int *IOPT;
00237 IOPT = ivector (1, 3);
00238 IOPT[1] = 1;
00239 IOPT[2] = 1;
00240 IOPT[3] = inhomogen;
00241
00242 int MITER = 1;
00243 int METH;
00244 int MAXORD;
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
00258
00259
00260
00261
00262
00263 int LRW;
00264 int LIW;
00265
00266 if (IOPT[1] == 0)
00267 {
00268 LRW = 20 + nVar * (MAXORD + 1) + 3 * nVar + nVar * nVar + 2 + 1;
00269 LIW = 20 + nVar;
00270 }
00271 if (IOPT[1] == 1)
00272 {
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
00286
00287
00288
00289
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;
00308
00309
00310
00311
00312
00313
00314 double T = t0;
00315 double TOUT;
00316
00317
00318 if (dmds)
00319 {
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;
00331 else
00332 TOUT = tmesh[k];
00333
00334
00335
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;
00347
00348
00349 int jlong;
00350 if (k > nmesh)
00351 {
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
00359 for (i = 1; i <= nVar; i++)
00360 dydp[i][jlong] = Yt[j][i];
00361
00362 }
00363
00364 for (j = 1; j <= nVar; j++)
00365
00366 for (i = 1; i <= nVar; i++)
00367 {
00368 dyds[i][j] = Yt[M_PAR_SENS + j][i];
00369 }
00370 }
00371 else
00372 {
00373 _dcopy (y[k], Yt[0], nVar);
00374
00375
00376
00377
00378
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 }
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 }
00411 }
00412 }
00413
00414 T = TOUT;
00415 }
00416
00417 if (dgdy)
00418 free_dmatrix (dgdy, 1, NOBS, 1, nVar);
00419 if (dgdp)
00420 free_dmatrix (dgdp, 1, NOBS, 1, nPExp);
00421
00422 }
00423 else
00424 {
00425
00426 for (k = 1; k <= nmesh + 1; k++)
00427 {
00428 if (k > nmesh)
00429 TOUT = t1;
00430 else
00431 TOUT = tmesh[k];
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;
00442
00443
00444 if (k > nmesh)
00445 _dcopy (zustand, Yt[0], nVar);
00446 else
00447 {
00448 _dcopy (y[k], Yt[0], nVar);
00449
00450
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;
00456 }
00457
00458 }
00459
00460
00462
00463
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
00474
00475 free_ivector (NEQ, 1, IPAR_offset + nPExp);
00476 free_ivector (IWORK, 1, LIW);
00477 free_ivector (IOPT, 1, 3);
00478
00479
00480
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 }
00500
00501
00502
00503
00504
00505
00506
00507
00508 extern "C" void ode_ (int &n, int &m, double &x, double *y, double *dydx, double *p)
00509 {
00510 ode (x, y - 1, dydx - 1, p - 1);
00511 }
00512
00513 extern "C" void jac_ (int *NEQ, double &t, double *y, int &n, double *p, int &M, double *pd)
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 }
00521
00522 extern "C" void od_ode_ (int *NEQ, double &t, double *y, int &N, double *p,
00523 int &M, double *dfdp, int &j)
00524 {
00525 if (!j)
00526 return;
00527 inhomo1 (t, y - 1, dfdp - 1, p - 1, j);
00528 }
00529
00530 extern "C" void ad_ode_ (int &, int &, int &, double &, double *, double *, double *,
00531 int &, double *, double *, int &)
00532 {
00533 cerr << "unexpected call to AD_ODE" <<endl;
00534 }