/home/peifer/diffit/modules/dampIt.cc

Go to the documentation of this file.
00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<string.h>
00005 #include<stdio.h>
00006 
00007 #include "../def.h"
00008 #include "../model.h"
00009 #include "../nr.h"
00010 
00011 using namespace std;
00012 
00013 
00014 
00015 //Begin definition of module prototypes
00016 
00017 //Module modules/setInitialValues.cc
00018 
00019 void intODE(GlobExp *ex,Glob *globs,int doDerivs,int doPlot,long expNr);
00020 
00021 double computeRight(Glob *globs,GlobExp *ex);
00022 
00023 void solvLin(Glob *globs,GlobExp *ex,int computeCovar);
00024 
00025 //End definition of module prototypes
00026 
00027 //norm of dx
00028 
00029 double Norm(Glob *globs, GlobExp *ex,double ***dS, double **dP,double uS) 
00030 {
00031   double norm=0;
00032   if (!dS) 
00033     {
00034     cerr <<"dS==NULL in Norm\n";
00035     throw 1;
00036     }
00037   long i,ivar,iexp,ip;
00038   for (iexp=1;iexp<=globs->nrExp;++iexp)
00039     {
00040       for (i=1; i< ex[iexp].nPoints; ++i)
00041         for (ivar=1; ivar<=ex[iexp].nvar; ++ivar)
00042           norm+=pow(uS*dS[iexp][i][ivar],2);
00043 
00044       for (ip=1;ip<=globs->npar; ++ip)
00045         {
00046           if(globs->doP[ip]=='L')
00047             norm+=pow(dP[iexp][ip],2);
00048         }
00049     } // end of: for (iexp=1;iexp<=globs->nrExp;++iexp)
00050   for (ip=1;ip<=globs->npar; ++ip)
00051     {
00052       if(globs->doP[ip]==TRUE)
00053         norm+=pow(dP[1][ip],2);
00054     }
00055   return sqrt(norm);
00056   
00057 } // Norm
00058 
00059 
00060 // perform one damping subiteration
00061 
00062 double subIt(Glob *globs, GlobExp *ex, double t, double normdx,
00063              double **p0, double ***s0, double **dp, double ***ds,double uS) 
00064 {
00065   long iexp,i,j;
00066   long maxPoints=0,maxVar=0;
00067   long nP=globs->npar;
00068 
00069   for (i=1;i<=globs->nrExp;++i) 
00070     {
00071       if (ex[i].nPoints>maxPoints) 
00072         maxPoints=ex[i].nPoints;
00073       if (ex[i].nvar>maxVar)       
00074         maxVar=ex[i].nvar;
00075     }
00076   
00077   double **newdP=dmatrix(1,globs->nrExp,1,nP);                      // parameter updates
00078   double ***newdS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar); // yTry updates
00079   
00080   
00081   for(iexp=1;iexp<=globs->nrExp;iexp++) 
00082     {
00083       GlobExp *expi=&ex[iexp];
00084       for(i=1; i< expi->nPoints; i++)
00085         {
00086           for(j=1; j<=expi->nvar; j++)
00087             {
00088               expi->yTry[i][j] = s0[iexp][i][j]+t*uS*ds[iexp][i][j];
00089               // positiveness of initial values --> test
00090               if(expi->yTry[i][j] < 0.0)
00091                 expi->yTry[i][j] = 0.0;
00092               expi->h[i][j]*=uS;
00093             }
00094         }
00095       for(i=1; i<=nP; i++)
00096         {
00097           expi->par[i] = p0[iexp][i]+t*dp[iexp][i];
00098         }
00099     }
00100 
00101   //integrate and compute value of objective function
00102   globs->chisq=0.;
00103   for (iexp=1;iexp<=globs->nrExp;iexp++) 
00104     {
00105       intODE(&ex[iexp],globs,FALSE,FALSE,iexp);
00106       globs->chisq+=computeRight(globs,&ex[iexp]);
00107       
00108     }
00109   *dbg  <<"subtotal chisq=" << globs->chisq << endl;
00110   solvLin(globs,ex,FALSE);
00111 
00112   for (iexp=1;iexp<=globs->nrExp;iexp++) 
00113     {      
00114       GlobExp *expi=&ex[iexp];
00115       for (i=1; i< expi->nPoints; i++)
00116         {
00117           for (j=1; j<=expi->nvar; j++)
00118             {
00119               newdS[iexp][i][j]=ex[iexp].dS[i][j];
00120             }
00121         }
00122       for(i=1; i<=nP; i++)
00123         {
00124         newdP[iexp][i]=ex[iexp].dP[i];
00125         }
00126     }
00127   //new increments in newdP and newdS
00128 
00129   //w-estimator
00130   double sum=0;
00131   for (iexp=1;iexp<=globs->nrExp;++iexp) 
00132     {
00133       GlobExp *expi=&ex[iexp];
00134       for (i=1; i< expi->nPoints; ++i)
00135         for (j=1; j<=expi->nvar; ++j)
00136           sum+=pow(newdS[iexp][i][j]-(1-t)*uS*ds[iexp][i][j],2);
00137       for (i=1; i<=nP; ++i)
00138         {
00139           if(globs->doP[i]=='L')
00140             sum+=pow(newdP[iexp][i]-(1-t)*dp[iexp][i],2);
00141         }
00142     }
00143   for (i=1; i<=nP; ++i)
00144     {
00145       if(globs->doP[i]==TRUE)
00146         sum+=pow(newdP[1][i]-(1-t)*dp[1][i],2);
00147     }
00148         
00149   free_d3tensor(newdS,1,globs->nrExp,1,maxPoints,1,maxVar);
00150   free_dmatrix(newdP,1,globs->nrExp,1,nP);
00151   return 2*sqrt(sum)/(pow(t*normdx,2));
00152 } // subIt
00153 
00154 
00155 double dampIt(Glob *globs,GlobExp *ex,double **P0, double ***S0, 
00156               double **dP, double ***dS,double *uS)
00157 {
00158   int etaOk=FALSE,i=1;
00159   double tau=0.95;
00160   double t;
00161   double tau_min=0.1;
00162   double eta0=1.;
00163   double eta2=1.8; //should be parsed
00164 
00165   double normdx=Norm(globs,ex,dS,dP,*uS);
00166   // damping algorithm
00167 
00168   if(globs->wquer<0)
00169     {
00170       globs->wquer=1.*eta0/normdx;
00171     }
00172   
00173   else
00174     t=eta0/(globs->wquer*normdx);
00175   
00176   if (t > tau) 
00177     {
00178       t=1;
00179     }
00180   else if(t <= tau_min)
00181     {
00182       t=tau_min;
00183     }
00184  
00185   globs->wquer=subIt(globs, ex, t,normdx,P0, S0, dP, dS,*uS);
00186   if(globs->silent!=TRUE)
00187     {
00188       cout << "Start damping:\n";
00189       cout << "#" << i << " (Predicted)  t = " << t << "\t";
00190       cout << "chisq=" << globs->chisq;
00191       cout << endl;
00192     }
00193 
00194   etaOk=((globs->wquer*t*normdx)<eta2);
00195   i++;
00196 
00197   while (!etaOk) 
00198     {
00199       t=eta0/(normdx*globs->wquer);
00200         
00201       if (t > tau) 
00202         {
00203           t=1;
00204           etaOk=TRUE;
00205           continue;
00206         }
00207       else if(t < tau_min)
00208         {
00209               t=tau_min;
00210               etaOk=TRUE;
00211               if(globs->silent!=TRUE)
00212                 {
00213                   cout << "#" << i <<  " (Corrected)  t = " << t << "\t";
00214                   cout << "\nMinimal damping factor tau_min= " << tau_min << " reached.\n";
00215                 }
00216               continue;
00217         }
00218       globs->wquer=subIt(globs, ex, t, normdx,P0, S0, dP, dS,*uS);
00219       if(globs->silent!=TRUE)
00220         {
00221           cout << "#" << i <<  " (Corrected)  t = " << t << "\t";
00222           cout << "chisq=" << globs->chisq;
00223           cout << endl;
00224         }
00225       etaOk=((globs->wquer*t*normdx)<eta2);
00226       i++;
00227       if(i>=10)
00228         {
00229           cerr << "Too many corrections.\n";
00230           throw 1;
00231         }
00232     }
00233   return(t);
00234 }

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