00001
00011 #include<string.h>
00012 #include<stdlib.h>
00013 #include<stddef.h>
00014 #include<iostream>
00015 #include<fstream>
00016 #include<getopt.h>
00017 #include<string>
00018 #include<math.h>
00019
00020 #include "../nr.h"
00021 #include "../def.h"
00022 #include "../model.h"
00023
00025 #define MAXLENGTH 1000
00026
00027 using namespace std;
00028
00033 char usage[]="\n\t --- diffit : Fitting parameters in differential equations ---\n \
00034 diffit, v 2.3 12/Apr./2006 Univ. of Freiburg, \n \
00035 Authors: Martin Peifer, Kilian Bartholome\n\n \
00036 SYNOPSIS \n \t diffit [Options] <File>\n\n \
00037 DESCRIPTION\n \
00038 \t -h -? -help \t help\n \
00039 \t -x0 <val> \t fit starts at time <val>\n \
00040 \t -x1 <val> \t fit stops at time <val>\n\
00041 \t -nms <num> \t number of multiple shooting intervals \n\
00042 \t -f <File> \t option file/multi-experiment definition file\n\
00043 \t\t Experiments are separated by @, comments start with #\n\
00044 \t\t Experiment specific parameters are : x0,x1,nms,y0,nobs\n\
00045 \t\t p,spline\n\
00046 \t -p <list> \t initial parameters, e.g. -p 0.5,2.3\n \
00047 \t -y0 <list> \t initial values, e.g. -y0 1,2\n\
00048 \t -nobs <num> \t number of observed quntities\n\
00049 \t -eps <val> \t integration accuracy\n\
00050 \t -nognu \t turns gnuplot animation off\n\
00051 \t -nomeasure \t using evently spaced multiple shooting intervals\n\
00052 \t -doP <list> \t specifies paramaters to be fitted, e.g. -doP 110L\n\
00053 \t -maxit <num> \t maximal number of iterations\n\
00054 \t -wait \t\t wait after each iteration\n\
00055 \t -usesig \t use weigths from input file\n\
00056 \t -int <num> \t integrator: 1-> Runge-Kutta (default)\n \
00057 \t \t 2-> CVODES \n \
00058 \t \t 3-> ODESSA \n \
00059 \t -maxstp \t maximal integration steps\n\
00060 \t -minimp \t minimal improvement for iteration truncation \n\
00061 \t -nowait \t do not wait until keypressed after iteration\n\
00062 \t -elastic <val> a factor (0< <val> <= 1) weakens cont. constraints \n\
00063 \t -reg \t\t regularisation for ill posed problems \n\
00064 \t -epsilon <val> singular value threshold, default: 1e-10\n\
00065 \t -lambda <val> regularisation parameter, default: 1e6\n\
00066 \t -spline <list> specifies spline data for non-autonomous ode's \n\
00067 \t -siminit \t simulate initial state\n\
00068 \t -pert <val>\t perturbate initial try - only in combination with siminit\n\
00069 \t -y0fix <list> fixing initial state\n\
00070 \t -nodamp \t switches damping off\n\
00071 \t -strat <var> optimisation strategy:\n\
00072 \t 1: MS 2: SRES (global) \n\
00073 \t -opt <num> \t selects the minimiser:\n\
00074 \t 1: LSEI 2: NAG (default)\n\
00075 \t -Lconst <list> set local Parameter Constraints, \n\
00076 \t e.g. -Lconst 2,4,3\n";
00077
00078 long get_list(long n,double *arg,char *string,char *desc)
00079 {
00080 long k=1,l=0,ind=0;
00081 char sep=',';
00082 char *dum=new char[strlen(string)+1];
00083 while(string[ind]!='\0')
00084 {
00085 if(k>n)
00086 {
00087 cerr << "Too many arguments in" << desc << "list." << endl;
00088 exit(1);
00089 }
00090 dum[l]=string[ind];
00091
00092 if(dum[l]==sep)
00093 {
00094 dum[l]='\0';
00095 arg[k]=atof(dum);
00096 l=0;
00097 k++;
00098 ind++;
00099 }
00100 else
00101 {
00102 ind++;
00103 l++;
00104 }
00105 }
00106 dum[l]='\0';
00107 arg[k]=atof(dum);
00108
00109 delete dum;
00110 return(k);
00111 }
00112
00120 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr)
00121 {
00122 int longindex,opt;
00123 int parlistSpecified=FALSE;
00124 long k,l,nExp,sp;
00125 char name[MAXLENGTH],line[MAXLENGTH];
00126
00127
00128 static struct option longopts[]={
00129 {"help", 0, 0, 'h'},
00130 {"help", 0, 0, '?'},
00131 {"x0" , 1, 0, 1 },
00132 {"x1" , 1, 0, 2 },
00133 {"nms" ,1, 0, 3 },
00134 {"f" ,1, 0, 4 },
00135 {"p" ,1, 0, 5 },
00136 {"y0" ,1, 0, 6 },
00137 {"nobs" ,1, 0, 7 },
00138 {"eps" ,1, 0, 8 },
00139 {"nognu",0, 0, 9 },
00140 {"nomeasure",0,0,10},
00141 {"doP" ,1, 0, 11},
00142 {"maxit",1, 0, 12},
00143 {"wait" ,0, 0, 13},
00144 {"usesig",0,0, 14},
00145 {"int" ,1,0, 15},
00146 {"maxstp",1,0, 17},
00147 {"minimp",1,0, 18},
00148 {"nowait",0,0, 19},
00149 {"elastic" ,1,0, 21},
00150 {"reg" ,0,0, 22},
00151 {"epsilon" ,1,0, 23},
00152 {"lambda" ,1,0, 24},
00153 {"spline" ,1,0, 25},
00154 {"siminit" ,0,0, 26},
00155 {"pert" ,1,0, 27},
00156 {"y0fix" ,1,0, 28},
00157 {"nodamp" ,0,0, 29},
00158 {"strat" ,1,0, 30},
00159 {"opt" ,1,0, 31},
00160 {"Lconst" ,1,0, 32},
00161 {0, 0, 0, 0}
00162 };
00163
00164 globs->noGnu=FALSE;
00165 globs->eps=1e-6;
00166 globs->npar=NPARAMS;
00167 globs->noMeasurements=FALSE;
00168 globs->doP=ivector(1,globs->npar);
00169 for(k=1;k<=globs->npar;k++)
00170 globs->doP[k]=TRUE;
00171 globs->y0fix=ivector(1,NEQNS);
00172 for(k=1;k<=NEQNS;k++)
00173 globs->y0fix[k]=TRUE;
00174 globs->maxit=1000;
00175 globs->gnuFp=NULL;
00176 globs->wait=FALSE;
00177 globs->usesig=FALSE;
00178 globs->integrator=1;
00179 globs->stiff=TRUE;
00180 globs->maxstp=5000;
00181 globs->minimp=1e-4;
00182 globs->nowait=FALSE;
00183 globs->elastic=1.;
00184 globs->reg=FALSE;
00185 globs->epsilon=1e-10;
00186 globs->lambda=1e6;
00187 globs->simInit=FALSE;
00188 globs->pert=0.;
00189 globs->nodamp=FALSE;
00190 globs->initSpline=TRUE;
00191 globs->wquer=-1;
00192 globs->silent=FALSE;
00193 globs->strategy=1;
00194 globs->minimiser=2;
00195 globs->faktorL=dvector(1,NPARAMS);
00196 globs->faktorLexist=FALSE;
00197 for(k=1;k<=NPARAMS;k++)
00198 {
00199 globs->faktorL[k]=-1;
00200 }
00201 while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00202 {
00203 switch(opt)
00204 {
00205 case 4:
00206 parlistSpecified=TRUE;
00207 strncpy(name,optarg,MAXLENGTH);
00208 break;
00209 }
00210 }
00211 optind=0;
00212
00213 GlobExp *ex;
00214 if(parlistSpecified==FALSE)
00215 {
00216 globs->nrExp=1;
00217 ex=new GlobExp[globs->nrExp+1];
00218
00219 ex[1].fitstart=-1e99;
00220 ex[1].fitend=1e99;
00221 ex[1].nobs=NOBS;
00222 ex[1].nvar=NEQNS;
00223 ex[1].nPoints=2;
00224 ex[1].splinesDefined=FALSE;
00225 ex[1].y0=NULL;
00226
00227 ex[1].par=dvector(1,NPARAMS);
00228 for(k=1;k<=globs->npar;k++)
00229 ex[1].par[k]=DefParameters[k-1];
00230
00231 if(NSPLINES!=0)
00232 {
00233 ex[1].nNodes=lvector(1,NSPLINES);
00234 }
00235 }
00236 else
00237 {
00238 ifstream in; ofstream out;
00239 char str[MAXLENGTH];
00240
00241 in.open(name);
00242
00243 if(!in)
00244 {
00245 cerr << "parameter list " << name << " not found\n";
00246 exit(1);
00247 }
00248
00249
00250 globs->nrExp=0;
00251 sprintf(str,"%s.%d",name,globs->nrExp);
00252 out.open(str);
00253 while(!in.eof())
00254 {
00255 in.getline(str,MAXLENGTH,'\n');
00256 if(str[0]=='@')
00257 {
00258 out.flush();
00259 out.close();
00260 globs->nrExp++;
00261 sprintf(str,"%s.%d",name,globs->nrExp);
00262 out.open(str);
00263 }
00264 else if(str[0]!='#')
00265 {
00266 out << str << endl;
00267 }
00268 }
00269 out.flush();
00270 out.close();
00271 in.close();
00272
00273 if(globs->nrExp==0)
00274 {
00275 cerr << "insuffient number of experiments, @ missing\n";
00276 exit(1);
00277 }
00278
00279 long _argc;
00280 char *_argv[MAXLENGTH];
00281
00282 for(k=0;k<=MAXLENGTH;k++)
00283 _argv[k]=(char *)malloc((size_t) (150*sizeof(char)));
00284
00285 {
00286 ifstream inp;
00287
00288 _argc=0;
00289 sprintf(str,"%s.0",name);
00290 inp.open(str);
00291
00292
00293 while(!inp.eof())
00294 {
00295 inp >> _argv[_argc+1];
00296 _argc++;
00297 }
00298 inp.close();
00299 }
00300
00301
00302
00303 while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00304 {
00305 switch(opt)
00306 {
00307 case 'h':
00308 cout << usage << endl;
00309 exit(1);
00310 break;
00311 case '?':
00312 cerr << usage << endl;
00313 exit(1);
00314 break;
00315 case 8:
00316 globs->eps=fabs(atof(optarg));
00317 break;
00318 case 9:
00319 globs->noGnu=TRUE;
00320 break;
00321 case 10:
00322 globs->noMeasurements=TRUE;
00323 break;
00324 case 11:
00325 if(strlen(optarg)!=globs->npar)
00326 {
00327 cerr << strlen(optarg) << " parameter(s) specified instead of " << globs->npar << ".\n";
00328 exit(1);
00329 }
00330 for(k=0;k < globs->npar;k++)
00331 {
00332 if(optarg[k]=='0')
00333 globs->doP[k+1]=FALSE;
00334 else if (optarg[k]=='L')
00335 globs->doP[k+1]='L';
00336 else
00337 globs->doP[k+1]=TRUE;
00338 }
00339 break;
00340 case 12:
00341 globs->maxit=abs(atol(optarg));
00342 break;
00343 case 13:
00344 globs->wait=TRUE;
00345 break;
00346 case 14:
00347 globs->usesig=TRUE;
00348 break;
00349 case 15:
00350 globs->integrator=abs(atoi(optarg));
00351 break;
00352 case 17:
00353 globs->maxstp=abs(atoi(optarg));
00354 break;
00355 case 18:
00356 globs->minimp=fabs(atof(optarg));
00357 break;
00358 case 19:
00359 globs->nowait=TRUE;
00360 break;
00361 case 21:
00362 globs->elastic=fabs(atof(optarg));
00363 if(globs->elastic > 1. || globs->elastic == 0.)
00364 {
00365 cerr << "Insufficient range of -elast <val> \n";
00366 exit(1);
00367 }
00368 break;
00369 case 22:
00370 globs->reg=TRUE;
00371 break;
00372 case 23:
00373 globs->epsilon=fabs(atof(optarg));
00374 break;
00375 case 24:
00376 globs->lambda=fabs(atof(optarg));
00377 break;
00378 case 26:
00379 globs->simInit=TRUE;
00380 break;
00381 case 27:
00382 globs->pert=fabs(atof(optarg));
00383 break;
00384 case 28:
00385 if(strlen(optarg)!=NEQNS)
00386 {
00387 cerr << strlen(optarg) << " variable(s) specified instead of " << NEQNS << ".\n";
00388 exit(1);
00389 }
00390 for(k=0;k < NEQNS;k++)
00391 {
00392 if(optarg[k]=='0')
00393 globs->y0fix[k+1]=FALSE;
00394 else
00395 globs->y0fix[k+1]=TRUE;
00396 }
00397 break;
00398 case 29:
00399 globs->nodamp=TRUE;
00400 break;
00401 case 30:
00402 globs->strategy=abs(atoi(optarg));
00403 break;
00404 case 31:
00405 globs->minimiser=abs(atoi(optarg));
00406 break;
00407 case 32:
00408 get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00409 globs->faktorLexist=TRUE;
00410 break;
00411 default:
00412 cerr << endl;
00413 cerr << "Parsing parameter list produced errors.\n\n";
00414 exit(1);
00415 }
00416 }
00417
00418 optind=0;
00419 ex=new GlobExp[globs->nrExp+1];
00420
00421
00422 for(nExp=1;nExp<=globs->nrExp;nExp++)
00423 {
00424
00425 ex[nExp].fitstart=-1e99;
00426 ex[nExp].fitend=1e99;
00427 ex[nExp].nobs=NOBS;
00428 ex[nExp].nvar=NEQNS;
00429 ex[nExp].nPoints=2;
00430 ex[nExp].splinesDefined=FALSE;
00431 ex[nExp].y0=NULL;
00432
00433 ex[nExp].par=dvector(1,NPARAMS);
00434 for(k=1;k<=globs->npar;k++)
00435 ex[nExp].par[k]=DefParameters[k-1];
00436
00437 if(NSPLINES!=0)
00438 {
00439 ex[nExp].nNodes=lvector(1,NSPLINES);
00440 }
00441
00442 _argc=0;
00443
00444 ifstream inp;
00445 sprintf(str,"%s.%d",name,nExp);
00446 inp.open(str);
00447
00448 while(!inp.eof())
00449 {
00450 inp >> _argv[_argc+1];
00451 _argc++;
00452 }
00453
00454 inp.close();
00455 while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00456 {
00457 switch(opt)
00458 {
00459 case 1:
00460 ex[nExp].fitstart=atof(optarg);
00461 break;
00462 case 2:
00463 ex[nExp].fitend=atof(optarg);
00464 break;
00465 case 3:
00466 ex[nExp].nPoints=abs(atol(optarg));
00467 ex[nExp].nPoints+=1;
00468 break;
00469 case 5:
00470 get_list(globs->npar,ex[nExp].par,optarg," parameter ");
00471 break;
00472 case 6:
00473 ex[nExp].y0=dvector(1,ex[nExp].nvar);
00474 for(k=1;k<=ex[nExp].nvar;k++)
00475 ex[nExp].y0[k]=0.;
00476 k=get_list(ex[nExp].nvar,ex[nExp].y0,optarg," initial values ");
00477 if(k!=ex[nExp].nvar)
00478 {
00479 cerr << ex[nExp].nvar << " initial values required.\n";
00480 exit(1);
00481 }
00482 break;
00483 case 7:
00484 ex[nExp].nobs=atol(optarg);
00485 break;
00486 case 25:
00487 k=0;
00488 l=1;
00489 ex[nExp].splinesDefined=TRUE;
00490 while(optarg[k]!='\0')
00491 {
00492 if(optarg[k]==',')
00493 l++;
00494 k++;
00495 }
00496 if(l!=NSPLINES)
00497 {
00498 cerr << "Incompatible number of splines.\n";
00499 exit(1);
00500 }
00501 ex[nExp].splineFile=new string[l+1];
00502 l=1;
00503 sp=0;
00504 k=0;
00505 while(optarg[k]!='\0')
00506 {
00507 if(optarg[k]==',')
00508 {
00509 line[sp]='\0';
00510 ex[nExp].splineFile[l]=(string)line;
00511 sp=0;
00512 l++;
00513 }
00514 else
00515 {
00516 line[sp]=optarg[k];
00517 sp++;
00518 }
00519 k++;
00520 }
00521 line[sp]='\0';
00522 ex[nExp].splineFile[l]=(string)line;
00523 break;
00524 default:
00525 cerr << endl;
00526 cerr << "Parsing parameter list produced errors.\n\n";
00527 exit(1);
00528 }
00529 }
00530
00531
00532 if((_argc-optind)!=1)
00533 {
00534 cerr << "No/Too many datafile(s) specified for experiment " << nExp << " . \n\n";
00535 exit(1);
00536 }
00537 else
00538 {
00539 ex[nExp].fileName=new char[strlen(_argv[optind])+1];
00540 strcpy(ex[nExp].fileName,_argv[optind]);
00541 }
00542
00543 if(ex[nExp].fitstart >= ex[nExp].fitend)
00544 {
00545 cerr << "x0 must be smaller than x1.\n";
00546 exit(1);
00547 }
00548 optind=0;
00549 }
00550
00551
00552
00553
00554
00555 sprintf(str,"rm -f %s.*",name);
00556 system(str);
00557 }
00558
00559 optind=0;
00560
00561 while((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00562 {
00563 switch(opt)
00564 {
00565 case 'h':
00566 cout << usage << endl;
00567 exit(1);
00568 break;
00569 case '?':
00570 cerr << usage << endl;
00571 exit(1);
00572 break;
00573 case 1:
00574 ex[1].fitstart=atof(optarg);
00575 break;
00576 case 2:
00577 ex[1].fitend=atof(optarg);
00578 break;
00579 case 3:
00580 ex[1].nPoints=abs(atol(optarg));
00581 ex[1].nPoints+=1;
00582 break;
00583 case 4:
00584 break;
00585 case 5:
00586 get_list(globs->npar,ex[1].par,optarg," parameter ");
00587 break;
00588 case 6:
00589 ex[1].y0=dvector(1,ex[1].nvar);
00590 for(k=1;k<=ex[1].nvar;k++)
00591 ex[1].y0[k]=0.;
00592 k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
00593 if(k!=ex[1].nvar)
00594 {
00595 cerr << ex[1].nvar << " initial values required.\n";
00596 exit(1);
00597 }
00598 break;
00599 case 7:
00600 ex[1].nobs=atol(optarg);
00601 break;
00602 case 8:
00603 globs->eps=fabs(atof(optarg));
00604 break;
00605 case 9:
00606 globs->noGnu=TRUE;
00607 break;
00608 case 10:
00609 globs->noMeasurements=TRUE;
00610 break;
00611 case 11:
00612 if(strlen(optarg)!=globs->npar)
00613 {
00614 cerr << strlen(optarg) << " parameter(s) specified instead of " << globs->npar << ".\n";
00615 exit(1);
00616 }
00617 for(k=0;k < globs->npar;k++)
00618 {
00619 if(optarg[k]=='0')
00620 globs->doP[k+1]=FALSE;
00621 else if (optarg[k]=='L')
00622 globs->doP[k+1]='L';
00623 else
00624 globs->doP[k+1]=TRUE;
00625 }
00626 break;
00627 case 12:
00628 globs->maxit=abs(atol(optarg));
00629 break;
00630 case 13:
00631 globs->wait=TRUE;
00632 break;
00633 case 14:
00634 globs->usesig=TRUE;
00635 break;
00636 case 15:
00637 globs->integrator=abs(atoi(optarg));
00638 break;
00639 case 17:
00640 globs->maxstp=abs(atoi(optarg));
00641 break;
00642 case 18:
00643 globs->minimp=fabs(atof(optarg));
00644 break;
00645 case 19:
00646 globs->nowait=TRUE;
00647 break;
00648 case 21:
00649 globs->elastic=fabs(atof(optarg));
00650 if(globs->elastic > 1. || globs->elastic == 0.)
00651 {
00652 cerr << "Insufficient range of -elast <val> \n";
00653 exit(1);
00654 }
00655 break;
00656 case 22:
00657 globs->reg=TRUE;
00658 break;
00659 case 23:
00660 globs->epsilon=fabs(atof(optarg));
00661 break;
00662 case 24:
00663 globs->lambda=fabs(atof(optarg));
00664 break;
00665 case 25:
00666 k=0;
00667 l=1;
00668 ex[1].splinesDefined=TRUE;
00669 while(optarg[k]!='\0')
00670 {
00671 if(optarg[k]==',')
00672 l++;
00673 k++;
00674 }
00675 if(l!=NSPLINES)
00676 {
00677 cerr << "Incompatible number of splines.\n";
00678 exit(1);
00679 }
00680 ex[1].splineFile=new string[l+1];
00681 l=1;
00682 sp=0;
00683 k=0;
00684 while(optarg[k]!='\0')
00685 {
00686 if(optarg[k]==',')
00687 {
00688 line[sp]='\0';
00689 ex[1].splineFile[l]=(string)line;
00690 sp=0;
00691 l++;
00692 }
00693 else
00694 {
00695 line[sp]=optarg[k];
00696 sp++;
00697 }
00698 k++;
00699 }
00700 line[sp]='\0';
00701 ex[1].splineFile[l]=(string)line;
00702 break;
00703 case 26:
00704 globs->simInit=TRUE;
00705 break;
00706 case 27:
00707 globs->pert=fabs(atof(optarg));
00708 break;
00709 case 28:
00710 if(strlen(optarg)!=NEQNS)
00711 {
00712 cerr << strlen(optarg) << " variable(s) specified instead of " << NEQNS << ".\n";
00713 exit(1);
00714 }
00715 for(k=0;k < NEQNS;k++)
00716 {
00717 if(optarg[k]=='0')
00718 globs->y0fix[k+1]=FALSE;
00719 else
00720 globs->y0fix[k+1]=TRUE;
00721 }
00722 break;
00723 case 29:
00724 globs->nodamp=TRUE;
00725 break;
00726 case 30:
00727 globs->strategy=abs(atoi(optarg));
00728 break;
00729 case 31:
00730 globs->minimiser=abs(atoi(optarg));
00731 break;
00732 case 32:
00733 get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00734 globs->faktorLexist=TRUE;
00735 break;
00736 default:
00737 cerr << endl;
00738 cerr << "Parsing command line options produced errors \n\n";
00739 exit(1);
00740 }
00741 }
00742
00743 if(ex[1].fitstart >= ex[1].fitend)
00744 {
00745 cerr << "x0 must be smaller than x1.\n";
00746 exit(1);
00747 }
00748
00749
00750 if((argc-optind)!=1 && parlistSpecified==FALSE)
00751 {
00752 cerr << "No/Too many datafile(s) specified. \n\n";
00753 exit(1);
00754 }
00755 else if(parlistSpecified==FALSE)
00756 {
00757 ex[1].fileName=new char[strlen(argv[optind])+1];
00758 strcpy(ex[1].fileName,argv[optind]);
00759 }
00760
00761 if(NSPLINES!=0)
00762 {
00763 for(k=1;k<=globs->nrExp;k++)
00764 {
00765 if(ex[k].splinesDefined==FALSE)
00766 {
00767 cerr << "Please specify spline data in experiment " << k << ".\n";
00768 exit(1);
00769 }
00770 }
00771 }
00772
00773 return(ex);
00774 }
00775