00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<stdio.h>
00005
00006 #include "def.h"
00007 #include "model.h"
00008 #include "nr.h"
00009
00010 using namespace std;
00011
00012
00013
00014
00015
00016 void intODE(GlobExp *ex,Glob *globs,int doDerivs,int doPlot,long expNr);
00017
00018 double computeRight(Glob *globs,GlobExp *ex);
00019
00020 void solvLin(Glob *globs,GlobExp *ex,int computeCovar);
00021
00022 double dampIt(Glob *globs,GlobExp *ex,double **P0, double ***S0,
00023 double **dP, double ***dS,double *uS);
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 void fitIt(GlobExp ex[],Glob *globs)
00034 {
00035 long nExp,i,j,k,l,n;
00036 long tryOnceMore,ok;
00037 double objF, sum, scale, lambda;
00038 long maxPoints=0,maxVar=0;
00039 double minImprovement=globs->minimp;
00040
00041 double uS=1.;
00042
00043 for (i=1;i<=globs->nrExp;++i)
00044 {
00045 if (ex[i].nPoints>maxPoints)
00046 maxPoints=ex[i].nPoints;
00047 if (ex[i].nvar>maxVar)
00048 maxVar=ex[i].nvar;
00049 }
00050
00051 double **dPSave=dmatrix(1,globs->nrExp,1,globs->npar);
00052 double ***dSSave=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00053 double **P0=dmatrix(1,globs->nrExp,1,globs->npar);
00054 double ***S0=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00055 double **dP=dmatrix(1,globs->nrExp,1,globs->npar);
00056 double ***dS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00057
00058 *dbg << "Minimum norm for Gauss-Newton update: minImprovement = " << minImprovement << "\n\n";
00059
00060
00061
00062
00063 globs->nIter=1;
00064
00065 tryOnceMore=TRUE;
00066 while (globs->nIter < globs->maxit && tryOnceMore)
00067 {
00068 globs->gnuindex=1;
00069 if(globs->silent!=TRUE)
00070 {
00071 cout << "Integration # " << globs->nIter << '\n';
00072 for(j=1;j<=globs->nrExp;j++)
00073 {
00074 cout << "Parameter value(s) of experiment " << j << ":\n";
00075 for (i=1; i<=globs->npar; ++i)
00076 {
00077 cout << ex[j].par[i] << " ";
00078 }
00079 cout << endl;
00080 cout << "Initial value(s) of experiment " << j << ":\n";
00081 for (i=1; i<=ex[j].nvar; ++i)
00082 {
00083 cout << ex[j].yTry[1][i] << " ";
00084 }
00085 cout << endl;
00086 }
00087 cout << "\n";
00088 cout.flush();
00089 }
00090 *dbg << "Integration # " << globs->nIter << '\n';
00091 *dbg << "Parameter value(s):";
00092 for (i=1; i<=globs->npar; ++i)
00093 {
00094 if(globs->doP[i]!='L')
00095 *dbg << '\t' << ex[1].par[i];
00096 }
00097 *dbg << "\n\n";
00098
00099
00100 sum=0;
00101 for (nExp=1;nExp<=globs->nrExp;++nExp)
00102 {
00103 #ifdef PRINTINITVALUES
00104 *dbg << "Experiment #" << nExp << "\n";
00105 *dbg << "Initial values at mesh points:\n";
00106 for (i=1; i<ex[nExp].nPoints; ++i)
00107 {
00108 *dbg << ex[nExp].mesh[i] << ':';
00109 for (j=1; j<=ex[nExp].nvar; ++j) {
00110 *dbg << '\t' << ex[nExp].yTry[i][j];
00111 }
00112 *dbg << '\n';
00113 }
00114 dbg->flush();
00115 #endif
00116 intODE(&ex[nExp],globs,TRUE,TRUE,nExp);
00117
00118
00119 objF=computeRight(globs,&ex[nExp]);
00120 sum += objF;
00121
00122 *dbg << "Experiment #" << nExp << ": objective f ="
00123 << objF << "\n";
00124 *dbg << "Discrepancies: ";
00125 for (i=1; i<ex[nExp].nPoints; ++i)
00126 {
00127 for (j=1; j<=ex[nExp].nvar; ++j)
00128 *dbg << ' ' << ex[nExp].h[i][j];
00129 *dbg << ",";
00130 }
00131 *dbg << "\n\n";
00132 dbg->flush();
00133 }
00134
00135 if(globs->silent!=TRUE)
00136 cout << "Objective function f = " << sum << "\n";
00137 *dbg << "Global objective function f = " << sum << "\n";
00138
00139 for (nExp=1;nExp<=globs->nrExp;++nExp)
00140 {
00141 for (i=1; i<ex[nExp].nPoints; ++i)
00142 {
00143 for (j=1; j<=ex[nExp].nvar; ++j)
00144 {
00145 ex[nExp].yTrySave[i][j] = ex[nExp].yTry[i][j];
00146 }
00147 }
00148 }
00149
00150 solvLin(globs,ex,FALSE);
00151
00152
00153 for (nExp=1;nExp<=globs->nrExp;nExp++)
00154 {
00155 for (i=1; i<ex[nExp].nPoints; i++)
00156 {
00157 for (j=1; j<=ex[nExp].nvar; j++)
00158 {
00159 S0[nExp][i][j] = ex[nExp].yTry[i][j];
00160 dSSave[nExp][i][j]=ex[nExp].dS[i][j];
00161 dS[nExp][i][j]=ex[nExp].dS[i][j];
00162 }
00163 }
00164 for (i=1; i<=globs->npar; i++)
00165 {
00166 P0[nExp][i] = ex[nExp].par[i];
00167 dPSave[nExp][i]=ex[nExp].dP[i];
00168 dP[nExp][i]=ex[nExp].dP[i];
00169 }
00170 }
00171
00172 uS=1.;
00173
00174 if(!globs->nodamp)
00175 lambda=dampIt(globs,ex,P0,S0,dP,dS,&uS);
00176 else
00177 lambda=1.;
00178
00179 globs->Lambda=lambda;
00180
00181 globs->chisq=sum;
00182 *dbg << "Chisq = " << globs->chisq << endl;
00183
00184
00185 *dbg << "using lambda=" << lambda << "\n";
00186
00187 sum=0;
00188 for (nExp=1;nExp<=globs->nrExp;++nExp)
00189 {
00190 for (i=1; i<ex[nExp].nPoints; ++i)
00191 {
00192 for (j=1; j<=ex[nExp].nvar; ++j)
00193 {
00194 sum += (lambda*uS*dSSave[nExp][i][j])*(lambda*uS*dSSave[nExp][i][j]);
00195 }
00196 }
00197 for (i=1; i<=globs->npar; ++i)
00198 {
00199 if(globs->doP[i]=='L')
00200 sum += (lambda*dPSave[nExp][i])*(lambda*dPSave[nExp][i]);
00201 }
00202 }
00203 for (i=1; i<=globs->npar; ++i)
00204 {
00205 if(globs->doP[i]!='L')
00206 sum += (lambda*dPSave[1][i])*(lambda*dPSave[1][i]);
00207 }
00208
00209 tryOnceMore = (sqrt(sum)/lambda > minImprovement);
00210
00211 if(uS==0)
00212 tryOnceMore=TRUE;
00213
00214 if (tryOnceMore)
00215 {
00216
00217 for (nExp=1;nExp<=globs->nrExp;++nExp)
00218 {
00219 for (i=1; i<ex[nExp].nPoints; ++i)
00220 {
00221 for (j=1; j<=ex[nExp].nvar; ++j)
00222 {
00223
00224 ex[nExp].yTry[i][j] = S0[nExp][i][j]+lambda*uS*dSSave[nExp][i][j];
00225
00226 if(ex[nExp].yTry[i][j]<0.0)
00227 ex[nExp].yTry[i][j]=0;
00228 }
00229 }
00230 for (i=1; i<=globs->npar; ++i)
00231 {
00232 ex[nExp].par[i]=P0[nExp][i]+lambda*dPSave[nExp][i];
00233 }
00234 }
00235 }
00236
00237
00238
00239
00240
00241 ++globs->nIter;
00242 if(globs->silent!=TRUE)
00243 {
00244 cout << "Norm of update vector: " << sqrt(sum) <<"\n\n";
00245 cout.flush();
00246 }
00247 *dbg << "Norm of update vector: " << sqrt(sum) <<"\n\n";
00248
00249 dbg->flush();
00250 if(globs->wait)
00251 {
00252 cout << "<Press any key to continue>\n";
00253 getchar();
00254 }
00255 }
00256
00257 if (globs->nIter>=globs->maxit)
00258 {
00259 if(globs->silent!=TRUE)
00260 {
00261 cout << "fitit: no convergence after " << globs->nIter << " iterations\n";
00262 cout << "PRELIMINARY RESULTS! \n";
00263 }
00264 globs->fitConverged=FALSE;
00265
00266 }
00267 else
00268 globs->fitConverged=TRUE;
00269
00270 *dbg << "last call of solvLin to obtain the covariance matrix\n";
00271 if(globs->reg=TRUE)
00272 {
00273 globs->minimiser=1;
00274
00275 solvLin(globs,ex,TRUE);
00276 }
00277 else
00278 {
00279 globs->minimiser=1;
00280
00281 solvLin(globs,ex,TRUE);
00282 }
00283
00284 if(!globs->wait && !globs->nowait)
00285 {
00286 cout << "<Press any key to continue>\n";
00287 getchar();
00288 }
00289
00290
00291 if(globs->noGnu==FALSE)
00292 {
00293 for(j=1;j<=globs->ngnu;j++)
00294 pclose(globs->gnuFp[j]);
00295 }
00296
00297 free_dmatrix(P0,1,globs->nrExp,1,globs->npar);
00298 free_d3tensor(S0,1,globs->nrExp,1,maxPoints,1,maxVar);
00299 free_dmatrix(dPSave,1,globs->nrExp,1,globs->npar);
00300 free_d3tensor(dSSave,1,globs->nrExp,1,maxPoints,1,maxVar);
00301 free_dmatrix(dP,1,globs->nrExp,1,globs->npar);
00302 free_d3tensor(dS,1,globs->nrExp,1,maxPoints,1,maxVar);
00303 }