- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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