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

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