00001 #include<string.h>
00002 #include<stdlib.h>
00003 #include<stddef.h>
00004 #include<iostream>
00005 #include<fstream>
00006 #include<getopt.h>
00007 #include<math.h>
00008
00009 #include "../nr.h"
00010 #include "../def.h"
00011 #include "../model.h"
00012
00013 #define MAXLENGTH 1000
00014
00015 using namespace std;
00016
00017 char usage[]="\n\t --- simit : Simulating differential equations ---\n \
00018 simit, v 2.3 12/Apr./2006 Univ. of Freiburg, Author: Martin Peifer\n \n \
00019 SYNOPSIS \n \t simit [Options] <File>\n\n \
00020 DESCRIPTION\n \
00021 \t -h -? -help \t help\n \
00022 \t -f <File> \t parameter file for simulation \n \
00023 \t -x0 <val> \t integration starts at time <val>\n \
00024 \t -x1 <val> \t integration stops at time <val>\n \
00025 \t -dt <val> \t integration step\n \
00026 \t -p <list> \t parameters, e.g. -p 0.5,2.3\n \
00027 \t -y0 <list> \t initial values, e.g. -y0 1,2\n \
00028 \t -nobs <num> \t number of observed quantities\n \
00029 \t -eps <val> \t integration accuracy\n \
00030 \t -sig <val> \t noise level\n \
00031 \t -int <num> \t integrator: 1-> Runge-Kutta (default)\n \
00032 \t \t 2-> CVODES \n \
00033 \t \t 3-> ODESSA \n \
00034 \t -maxstp \t maximal integration steps\n \
00035 \t -spline <list> specifies spline data for non-autonomous ode's \n";
00036
00037 long get_list(long n,double *arg,char *string,char *desc)
00038 {
00039 long k=1,l=0,ind=0;
00040 char sep=',';
00041 char *dum=new char[strlen(string)+1];
00042 while(string[ind]!='\0')
00043 {
00044 if(k>n)
00045 {
00046 cerr << "Too many arguments in" << desc << "list." << endl;
00047 exit(1);
00048 }
00049 dum[l]=string[ind];
00050
00051 if(dum[l]==sep)
00052 {
00053 dum[l]='\0';
00054 arg[k]=atof(dum);
00055 l=0;
00056 k++;
00057 ind++;
00058 }
00059 else
00060 {
00061 ind++;
00062 l++;
00063 }
00064 }
00065 dum[l]='\0';
00066 arg[k]=atof(dum);
00067
00068 delete dum;
00069 return(k);
00070 }
00071
00072
00073 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr)
00074 {
00075 int longindex,opt;
00076 int parlistSpecified=FALSE;
00077 long k,l,nExp,sp;
00078 char name[MAXLENGTH],line[MAXLENGTH];
00079
00080
00081 static struct option longopts[]={
00082 {"help", 0, 0, 'h'},
00083 {"help", 0, 0, '?'},
00084 {"x0" , 1, 0, 1 },
00085 {"x1" , 1, 0, 2 },
00086 {"dt" , 1, 0, 3 },
00087 {"f" ,1, 0, 4 },
00088 {"p" ,1, 0, 5 },
00089 {"y0" ,1, 0, 6 },
00090 {"nobs" ,1, 0, 7 },
00091 {"eps" ,1, 0, 8 },
00092 {"sig" ,1, 0, 9 },
00093 {"maxit",1, 0, 12 },
00094 {"int" ,1, 0, 15 },
00095 {"maxstp",1,0, 17 },
00096 {"spline",1,0, 21 },
00097 {0, 0, 0, 0}
00098 };
00099
00100 globs->noGnu=FALSE;
00101 globs->eps=1e-6;
00102 globs->npar=NPARAMS;
00103 globs->noMeasurements=FALSE;
00104 globs->doP=ivector(1,globs->npar);
00105 for(k=1;k<=globs->npar;k++)
00106 globs->doP[k]=TRUE;
00107 globs->maxit=1000;
00108 globs->gnuFp=NULL;
00109 globs->wait=FALSE;
00110 globs->usesig=FALSE;
00111 globs->integrator=1;
00112 globs->stiff=TRUE;
00113 globs->maxstp=5000;
00114 globs->minimp=0.05;
00115 globs->nowait=FALSE;
00116 globs->elastic=1.;
00117 globs->reg=FALSE;
00118 globs->epsilon=1e-10;
00119 globs->lambda=1e6;
00120 globs->dt=0.1;
00121 globs->sig=0.;
00122 globs->nrExp=1;
00123 globs->initSpline=TRUE;
00124
00125 while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00126 {
00127 switch(opt)
00128 {
00129 case 4:
00130 parlistSpecified=TRUE;
00131 strncpy(name,optarg,MAXLENGTH);
00132 break;
00133 }
00134 }
00135 optind=0;
00136
00137 GlobExp *ex;
00138 if(parlistSpecified==FALSE)
00139 {
00140 globs->nrExp=1;
00141 ex=new GlobExp[globs->nrExp+1];
00142
00143 ex[1].fitstart=0;
00144 ex[1].fitend=1;
00145 ex[1].nobs=NOBS;
00146 ex[1].nvar=NEQNS;
00147 ex[1].y0=NULL;
00148 ex[1].par=dvector(1,NPARAMS);
00149 for(k=1;k<=globs->npar;k++)
00150 ex[1].par[k]=DefParameters[k-1];
00151 }
00152 else
00153 {
00154 ifstream in;
00155 ofstream out;
00156 char str[MAXLENGTH];
00157
00158 in.open(name);
00159
00160 if(!in)
00161 {
00162 cerr << "parameter list " << name << " not found\n";
00163 exit(1);
00164 }
00165
00166
00167 globs->nrExp=1;
00168 ex=new GlobExp[globs->nrExp+1];
00169
00170 ex[1].fitstart=0;
00171 ex[1].fitend=1;
00172 ex[1].nobs=NOBS;
00173 ex[1].nvar=NEQNS;
00174 ex[1].par=dvector(1,NPARAMS);
00175 for(k=1;k<=globs->npar;k++)
00176 ex[1].par[k]=DefParameters[k-1];
00177
00178 sprintf(str,"%s.%d",name,globs->nrExp);
00179 out.open(str);
00180 while(!in.eof())
00181 {
00182 in.getline(str,MAXLENGTH,'\n');
00183 if(str[0]!='#')
00184 {
00185 out << str << endl;
00186 }
00187 }
00188 out.flush();
00189 out.close();
00190 in.close();
00191
00192 long _argc;
00193 char *_argv[MAXLENGTH];
00194
00195 for(k=0;k<=MAXLENGTH;k++)
00196 _argv[k]=(char *)malloc((size_t) (500*sizeof(char)));
00197
00198
00199 ifstream inp;
00200
00201 _argc=0;
00202 sprintf(str,"%s.1",name);
00203 inp.open(str);
00204
00205
00206 while(!inp.eof())
00207 {
00208 inp >> _argv[_argc+1];
00209 _argc++;
00210 }
00211 inp.close();
00212
00213
00214
00215 optind=0;
00216 while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00217 {
00218 switch(opt)
00219 {
00220 case 'h':
00221 cout << usage << endl;
00222 exit(1);
00223 break;
00224 case '?':
00225 cerr << usage << endl;
00226 exit(1);
00227 break;
00228 case 1:
00229 ex[1].fitstart=atof(optarg);
00230 break;
00231 case 2:
00232 ex[1].fitend=atof(optarg);
00233 break;
00234 case 3:
00235 globs->dt=fabs(atof(optarg));
00236 break;
00237 case 4:
00238 break;
00239 case 5:
00240 get_list(globs->npar,ex[1].par,optarg," parameter ");
00241 break;
00242 case 6:
00243 ex[1].y0=dvector(1,ex[1].nvar);
00244 for(k=1;k<=ex[1].nvar;k++)
00245 ex[1].y0[k]=0.;
00246 k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
00247 if(k!=ex[1].nvar)
00248 {
00249 cerr << ex[1].nvar << " initial values required.\n";
00250 exit(1);
00251 }
00252 break;
00253 case 7:
00254 ex[1].nobs=atol(optarg);
00255 break;
00256 case 8:
00257 globs->eps=fabs(atof(optarg));
00258 break;
00259 case 9:
00260 globs->sig=fabs(atof(optarg));
00261 break;
00262 case 11:
00263 if(strlen(optarg)!=globs->npar)
00264 {
00265 cerr << strlen(optarg) << " parameter(s) specified instead of " << globs->npar << ".\n";
00266 exit(1);
00267 }
00268 for(k=0;k < globs->npar;k++)
00269 {
00270 if(optarg[k]=='0')
00271 globs->doP[k+1]=FALSE;
00272 else
00273 globs->doP[k+1]=TRUE;
00274 }
00275 break;
00276 case 12:
00277 globs->maxit=abs(atol(optarg));
00278 break;
00279 case 15:
00280 globs->integrator=abs(atoi(optarg));
00281 break;
00282 case 17:
00283 globs->maxstp=abs(atoi(optarg));
00284 break;
00285 case 21:
00286 k=0;
00287 l=1;
00288 ex[1].splinesDefined=TRUE;
00289 while(optarg[k]!='\0')
00290 {
00291 if(optarg[k]==',')
00292 l++;
00293 k++;
00294 }
00295 if(l!=NSPLINES)
00296 {
00297 cerr << "Incompatible number of splines.\n";
00298 exit(1);
00299 }
00300 ex[1].splineFile=new string[l+1];
00301 l=1;
00302 sp=0;
00303 k=0;
00304 while(optarg[k]!='\0')
00305 {
00306 if(optarg[k]==',')
00307 {
00308 line[sp]='\0';
00309 ex[1].splineFile[l]=(string)line;
00310 sp=0;
00311 l++;
00312 }
00313 else
00314 {
00315 line[sp]=optarg[k];
00316 sp++;
00317 }
00318 k++;
00319 }
00320 line[sp]='\0';
00321 ex[1].splineFile[l]=(string)line;
00322 break;
00323 default:
00324 cerr << endl;
00325 cerr << "Parsing command line options produced errors \n\n";
00326 exit(1);
00327 }
00328 }
00329 nExp=1;
00330
00331 if((_argc-optind)!=1)
00332 {
00333 cerr << "No/Too many datafile(s) specified for experiment " << nExp << " . \n\n";
00334 exit(1);
00335 }
00336 else
00337 {
00338 ex[nExp].fileName=new char[strlen(_argv[optind])+1];
00339 strcpy(ex[nExp].fileName,_argv[optind]);
00340 }
00341
00342 if(ex[nExp].fitstart >= ex[nExp].fitend)
00343 {
00344 cerr << "x0 must be smaller than x1.\n";
00345 exit(1);
00346 }
00347
00348 sprintf(str,"rm -f %s.*",name);
00349 system(str);
00350 }
00351
00352 optind=0;
00353
00354 while((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00355 {
00356 switch(opt)
00357 {
00358 case 'h':
00359 cout << usage << endl;
00360 exit(1);
00361 break;
00362 case '?':
00363 cerr << usage << endl;
00364 exit(1);
00365 break;
00366 case 1:
00367 ex[1].fitstart=atof(optarg);
00368 break;
00369 case 2:
00370 ex[1].fitend=atof(optarg);
00371 break;
00372 case 3:
00373 globs->dt=fabs(atof(optarg));
00374 break;
00375 case 4:
00376 break;
00377 case 5:
00378 get_list(globs->npar,ex[1].par,optarg," parameter ");
00379 break;
00380 case 6:
00381 ex[1].y0=dvector(1,ex[1].nvar);
00382 for(k=1;k<=ex[1].nvar;k++)
00383 ex[1].y0[k]=0.;
00384 k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
00385 if(k!=ex[1].nvar)
00386 {
00387 cerr << ex[1].nvar << " initial values required.\n";
00388 exit(1);
00389 }
00390 break;
00391 case 7:
00392 ex[1].nobs=atol(optarg);
00393 break;
00394 case 8:
00395 globs->eps=fabs(atof(optarg));
00396 break;
00397 case 9:
00398 globs->sig=fabs(atof(optarg));
00399 break;
00400 case 11:
00401 if(strlen(optarg)!=globs->npar)
00402 {
00403 cerr << strlen(optarg) << " parameter(s) specified instead of " << globs->npar << ".\n";
00404 exit(1);
00405 }
00406 for(k=0;k < globs->npar;k++)
00407 {
00408 if(optarg[k]=='0')
00409 globs->doP[k+1]=FALSE;
00410 else
00411 globs->doP[k+1]=TRUE;
00412 }
00413 break;
00414 case 12:
00415 globs->maxit=abs(atol(optarg));
00416 break;
00417 case 15:
00418 globs->integrator=abs(atol(optarg));
00419 break;
00420 case 17:
00421 globs->maxstp=abs(atoi(optarg));
00422 break;
00423 case 21:
00424 k=0;
00425 l=1;
00426 ex[1].splinesDefined=TRUE;
00427 while(optarg[k]!='\0')
00428 {
00429 if(optarg[k]==',')
00430 l++;
00431 k++;
00432 }
00433 if(l!=NSPLINES)
00434 {
00435 cerr << "Incompatible number of splines.\n";
00436 exit(1);
00437 }
00438 ex[1].splineFile=new string[l+1];
00439 l=1;
00440 sp=0;
00441 k=0;
00442 while(optarg[k]!='\0')
00443 {
00444 if(optarg[k]==',')
00445 {
00446 line[sp]='\0';
00447 ex[1].splineFile[l]=(string)line;
00448 sp=0;
00449 l++;
00450 }
00451 else
00452 {
00453 line[sp]=optarg[k];
00454 sp++;
00455 }
00456 k++;
00457 }
00458 line[sp]='\0';
00459 ex[1].splineFile[l]=(string)line;
00460 break;
00461 default:
00462 cerr << endl;
00463 cerr << "Parsing command line options produced errors \n\n";
00464 exit(1);
00465 }
00466 }
00467
00468 if(ex[1].fitstart >= ex[1].fitend)
00469 {
00470 cerr << "x0 must be smaller than x1.\n";
00471 exit(1);
00472 }
00473
00474
00475 if((argc-optind)!=1 && parlistSpecified==FALSE)
00476 {
00477 cerr << "No/Too many datafile(s) specified. \n\n";
00478 exit(1);
00479 }
00480 else if(parlistSpecified==FALSE)
00481 {
00482 ex[1].fileName=new char[strlen(argv[optind])+1];
00483 strcpy(ex[1].fileName,argv[optind]);
00484 }
00485
00486
00487 if(NSPLINES!=0)
00488 {
00489 if(ex[1].splinesDefined==FALSE)
00490 {
00491 cerr << "Please specify spline data.\n";
00492 exit(1);
00493 }
00494 }
00495
00496
00497 return(ex);
00498 }
00499