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
00016
00017
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
00026
00027
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 }
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 }
00058
00059
00060
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);
00078 double ***newdS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
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
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
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
00128
00129
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 }
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;
00164
00165 double normdx=Norm(globs,ex,dS,dP,*uS);
00166
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 }