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
00054
00055 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr);
00056
00057 void readData(GlobExp *ex,Glob *globs,long expNr);
00058
00059 void setMesh(GlobExp *ex,Glob *globs,long expNr);
00060
00061 void outFit(GlobExp ex[],Glob *globs);
00062
00063 void freeMem(GlobExp *ex,Glob *globs,int simit);
00064
00065 void initialise(GlobExp ex[],Glob *globs,int simit);
00066
00067 void simInit(GlobExp *ex,Glob *globs);
00068
00069
00070
00071
00072 void fitIt(GlobExp ex[],Glob *globs);
00073
00074 void globOpt(GlobExp ex[],Glob *globs);
00075
00076
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
00089 ex=parseopts(argc,argv,&globs,outstr);
00090
00091
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);
00099 dbg->unsetf(ios::floatfield);
00100
00101
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
00122 readData(ex,&globs,i);
00123
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
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
00164 if(globs.simInit==TRUE)
00165 {
00166 for (i=1;i<=globs.nrExp;++i)
00167 simInit(&ex[i],&globs);
00168 }
00169
00170
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)
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
00196 outFit(ex,&globs);
00197 if(!globs.noGnu)
00198 system("rm -f gnuout.dat");
00199
00200 freeMem(ex,&globs,FALSE);
00201 }