- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /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 }