/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 

Generated on Mon Jan 29 17:09:13 2007 for Diffit by  doxygen 1.4.6