/home/peifer/diffit/diffit.cc

Go to the documentation of this file.
00001 
00043 #include<iostream>
00044 #include<fstream>
00045 #include<stdlib.h>
00046 
00047 #include "def.h"
00048 #include "model.h"
00049 #include "nr.h"
00050 
00051 using namespace std;
00052 
00053 //Begin definition of module prototypes
00054 //Module modules/parse.cc
00055 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr);
00056 //Module modules/readData.cc
00057 void readData(GlobExp *ex,Glob *globs,long expNr);
00058 //Module modules/setMesh.cc
00059 void setMesh(GlobExp *ex,Glob *globs,long expNr);
00060 //Module modules/outFit.cc
00061 void outFit(GlobExp ex[],Glob *globs);
00062 //Module modules/freeMem.cc
00063 void freeMem(GlobExp *ex,Glob *globs,int simit);
00064 //Module modules/initialise.cc
00065 void initialise(GlobExp ex[],Glob *globs,int simit);
00066 //Module modules/simInit.cc
00067 void simInit(GlobExp *ex,Glob *globs);
00068 
00069 //End definition of module prototypes
00070 
00071 //The numerics subprogram
00072 void fitIt(GlobExp ex[],Glob *globs);
00073 //Global optimiser stuff
00074 void globOpt(GlobExp ex[],Glob *globs);
00075 
00076 // DEBUG stream
00078 ofstream *dbg;
00079 
00081 main(int argc, char *argv[])
00082 {
00083   long k,i,j;
00084   char outstr[100];
00085   Glob globs;
00086   GlobExp *ex;
00087   
00088   //  parseopt
00089   ex=parseopts(argc,argv,&globs,outstr);
00090 
00091   //open DEBUG stream
00092   dbg=new ofstream("diffit.dbg");
00093   if (!dbg) 
00094     {
00095       cerr << "Error opening DEBUG file.\n";
00096       exit(1);
00097     }
00098   dbg->precision(4);     // set output format for floating point data
00099   dbg->unsetf(ios::floatfield);
00100 
00101   //print debugging information
00102   *dbg << DefModelDescription << "\n";
00103   if (globs.nrExp==1) 
00104     *dbg << "1 experiment\n";
00105   else
00106     *dbg << globs.nrExp << " experiments\n";
00107    *dbg << "\n";
00108    for(i=1;i<=globs.nrExp;i++)
00109      {
00110        *dbg << "Experiment: " << i << "\n";
00111        *dbg << NPARAMS << " parameter(s):";
00112        for (k=1; k<=NPARAMS; ++k) 
00113          *dbg << " " << ex[i].par[k];
00114        *dbg << "\n";
00115      }
00116    
00117   dbg->flush();
00118  
00119   for(i=1;i<=globs.nrExp;i++)
00120     {
00121       // readData()
00122       readData(ex,&globs,i);
00123       //set mesh
00124       setMesh(ex,&globs,i);
00125 
00126 #ifdef PRINTDATA
00127       *dbg << "\nData: \n";
00128       for (j=1; j<=ex[i].nMeasure; j++) {
00129         *dbg << ex[i].xMeasure[j];
00130         for (k=1; k<=ex[i].nobs; k++)
00131           *dbg << "\t" << ex[i].yMeasure[j][k] << "\t" << ex[i].sigma[j][k];
00132         *dbg << "\n";
00133       }
00134       *dbg << "\n";
00135 #endif
00136       *dbg << "Mesh:";
00137       for (k=1;k<=ex[i].nPoints;++k) 
00138         *dbg << " " << ex[i].mesh[k];
00139       *dbg << "\n";
00140     }
00141 
00142   //initial values
00143   for(i=1;i<=globs.nrExp;i++)
00144     {
00145       *dbg << "\nExperiment " << i << ":\n";
00146       if (ex[i].y0) 
00147         {
00148           *dbg << NEQNS << " starting value(s):";
00149           for (k=1; k<=NEQNS; ++k) 
00150             *dbg << ' ' << ex[i].y0[k];
00151           *dbg << "\n";
00152         } 
00153       else 
00154         {
00155           *dbg << "No starting values specified\n";
00156         }
00157       *dbg << "\n";
00158       dbg->flush();
00159     }
00160   
00161    initialise(ex,&globs,FALSE);
00162    
00163    //simulate initial state
00164    if(globs.simInit==TRUE)
00165      {
00166        for (i=1;i<=globs.nrExp;++i) 
00167          simInit(&ex[i],&globs);
00168      }
00169    
00170    // starting of the numerics
00171    //*************************
00172       
00173    if(globs.strategy==2)
00174      globOpt(ex,&globs);
00175    else
00176      {
00177        try 
00178          {            
00179             fitIt(ex,&globs);
00180          }
00181        catch(int i)  //Exeption handling
00182          {
00183            if(globs.noGnu==FALSE)
00184              {
00185                for(j=1;j<=globs.ngnu;j++)
00186                  pclose(globs.gnuFp[j]);
00187              }       
00188            exit(1);
00189          }
00190      }
00191 
00192    
00193    //*************************
00194 
00195   //Output after convergence
00196   outFit(ex,&globs);
00197   if(!globs.noGnu)
00198     system("rm -f gnuout.dat");
00199 
00200   freeMem(ex,&globs,FALSE);
00201 }

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