- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/fitIt.cc

Go to the documentation of this file.

    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 //#define PRINTINITVALUES
    00013
    00014 //Begin definition of module prototypes
    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 //End definition of module prototypes
    00026
    00027
    00028
    00029
    00030 //      Main routine
    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   // GGN iteration
    00061   // do the actual fitting
    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       // all experiments
    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           // set up internal data (residues, etc.)
    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       // solve linearized minimisation problem
    00150       solvLin(globs,ex,FALSE);
    00151
    00152       //save old values and new updates
    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       //damping
    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       // next guess
    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           // set up new initial values for next iteration
    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                       // positiveness of initial cond --> test
    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       // else:
    00238       // keep old values; parCovar is covariance matrix at solution point
    00239       // (update is too small to be worth iterating once more)
    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   //to obtain the covariance matrix
    00270   *dbg << "last call of solvLin to obtain the covariance matrix\n";
    00271   if(globs->reg=TRUE)
    00272     {
    00273       globs->minimiser=1;
    00274       //globs->reg=FALSE;
    00275       solvLin(globs,ex,TRUE);
    00276     }
    00277   else
    00278     {
    00279       globs->minimiser=1;
    00280       //globs->reg=FALSE;
    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 }