- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/modules/parseSimit.cc

Go to the documentation of this file.

    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   //option definition
    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   //initialise some global parameters in globs
    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       //initialise some global parameters
    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 //if parlist specified
    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       //preprocessing of data list
    00167       globs->nrExp=1;
    00168       ex=new GlobExp[globs->nrExp+1];
    00169       //initialise some global parameters
    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       //PARSING LIST
    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       //Parsing file name
    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       //checking if x0 < x1
    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   //checking if x0 < x1
    00468   if(ex[1].fitstart >= ex[1].fitend)
    00469     {
    00470       cerr << "x0 must be smaller than x1.\n";
    00471       exit(1);
    00472     }
    00473
    00474   //Parsing file name for a single experiment
    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