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

#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "../def.h"
#include "../model.h"
#include "../nr.h"

Go to the source code of this file.

Functions

void intODE (GlobExp *ex, Glob *globs, int doDerivs, int doPlot, long expNr)
double computeRight (Glob *globs, GlobExp *ex)
void solvLin (Glob *globs, GlobExp *ex, int computeCovar)
double Norm (Glob *globs, GlobExp *ex, double ***dS, double **dP, double uS)
double subIt (Glob *globs, GlobExp *ex, double t, double normdx, double **p0, double ***s0, double **dp, double ***ds, double uS)
double dampIt (Glob *globs, GlobExp *ex, double **P0, double ***S0, double **dP, double ***dS, double *uS)


Function Documentation

double computeRight Glob globs,
GlobExp ex
 

Definition at line 39 of file computeRight.cc.

References GlobExp::h, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, and GlobExp::nvar.

00040 {
00041     long i,j,k;
00042 
00043     // calculate discrepancies, residues and constraints
00044     long nPoints=ex->nPoints, nvar=ex->nvar;
00045     long nMeasure=ex->nMeasure, nobs=ex->nobs;
00046     long nP=globs->npar;
00047 
00048     // no continuity constraints for first and last mesh point
00049     for (j=1;j<=nvar;++j) {
00050         ex->h[1][j]=0;
00051         ex->h[nPoints][j]=0;
00052     }
00053     for (i=2;i<nPoints;++i)
00054       {
00055         for (j=1;j<=nvar;++j)
00056           ex->h[i][j]=ex->yComp[i][j]-ex->yTry[i][j];
00057       }
00058     for (i=1;i<=nMeasure;++i)
00059       {
00060         for (j=1;j<=nobs;++j)
00061           ex->residues[i][j]=ex->yPred[i][j]-ex->yMeasure[i][j];
00062       }
00063 
00064     //extra equality and inequality constraints
00065     R2(globs,ex, TRUE);
00066     R3(globs,ex ,TRUE);
00067 
00068 #ifdef PRINTDERIVS
00069     for (i=1; i<=nPoints; ++i) {
00070         *dbg << "dyds[" << i << "] = (";
00071         for (j=1; j<=nvar; ++j) {
00072             for (k=1; k<=nvar; ++k)
00073                 *dbg << " " << ex->dyds[i][j][k];
00074             *dbg << ", ";
00075         }
00076         *dbg << ")\ndydp[" << i << "] = (";
00077         for (j=1; j<=nvar; ++j) {
00078             for (k=1; k<=nP; ++k)
00079                 *dbg << " " << ex->dydp[i][j][k];
00080             *dbg << ", ";
00081         }
00082         *dbg << ")\n";
00083     }
00084     for (i=1; i<=nMeasure; ++i) {
00085         *dbg << "dmds[" << i << "] = (";
00086         for (j=1; j<=nobs; ++j) {
00087             for (k=1; k<=nvar; ++k)
00088                 *dbg << " " << ex->dmds[i][j][k];
00089             *dbg << ", ";
00090         }
00091         *dbg << ")\ndmdp[" << i << "] = (";
00092         for (j=1; j<=nobs; ++j) {
00093             for (k=1; k<=nP; ++k)
00094                 *dbg << " " << ex->dmdp[i][j][k];
00095             *dbg << ", ";
00096         }
00097         *dbg << ")\n";
00098     }
00099     for (i=1; i<=ex->me; ++i) {
00100         *dbg << "dR2ds[" << i << "] = (";
00101         for (j=1; j<=nPoints; ++j) {
00102             for (k=1; k<=nvar; ++k)
00103                 *dbg << " " << ex->dR2ds[i][j][k];
00104             *dbg << ", ";
00105         }
00106         *dbg << ")\n";
00107         *dbg << "dR2dp[" << i << "] = (";
00108         for (j=1; j<=nP; ++j)
00109             *dbg << " " << ex->dR2dp[i][j];
00110         *dbg << ")\n";
00111     }
00112     for (i=1; i<=ex->mg; ++i) {
00113         *dbg << "dR3ds[" << i << "] = (";
00114         for (j=1; j<=nPoints; ++j) {
00115             for (k=1; k<=nvar; ++k)
00116                 *dbg << " " << ex->dR3ds[i][j][k];
00117             *dbg << ", ";
00118         }
00119         *dbg << ")\n";
00120         *dbg << "dR3dp[" << i << "] = (";
00121         for (j=1; j<=nP; ++j)
00122             *dbg << " " << ex->dR3dp[i][j];
00123         *dbg << ")\n";
00124     }
00125     dbg->flush();
00126 #endif
00127 
00128     return(objectiveF(globs,ex));
00129 }

double dampIt Glob globs,
GlobExp ex,
double **  P0,
double ***  S0,
double **  dP,
double ***  dS,
double *  uS
 

Definition at line 155 of file dampIt.cc.

References Glob::chisq, FALSE, Norm(), Glob::silent, subIt(), TRUE, and Glob::wquer.

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 }

void intODE GlobExp ex,
Glob globs,
int  doDerivs,
int  doPlot,
long  expNr
 

Definition at line 313 of file intODE.cc.

References d3tensor(), dmatrix(), GlobExp::dmdp, GlobExp::dmds, dvector(), GlobExp::dydp, GlobExp::dyds, maxTableLen(), GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, GlobExp::nvar, and GlobExp::xMeasure.

00314 {
00315   long i,j,k,ksav,l,m,n,nok,nbad;
00316   long nobs=ex->nobs;
00317   long nPoints=ex->nPoints, nvar=ex->nvar, nMeasure=ex->nMeasure;
00318   double *mesh=ex->mesh, *xMeasure=ex->xMeasure;
00319   long tableLen=maxTableLen(mesh, nPoints, xMeasure, nMeasure)+1; 
00320   double temp, h,ymin,ymax;
00321   
00322   // allocate memory
00323   double *ys=dvector(1,nvar);
00324   double *xTab=dvector(1,tableLen);
00325   double **yTab=dmatrix(1,tableLen,1,nvar);
00326   double ***dTabds, ***dTabdp;
00327   char name[50];
00328   ofstream gnuout;
00329   ofstream outout;
00330 
00331 
00332   //gnuout.open("gnuout.dat");
00333   
00334   // compute derivatives?
00335   // if not, don't write to graphics either
00336   // (no genuine integration, but only estimation of omega)
00337   if (doDerivs) 
00338     {
00339       dTabds=d3tensor(1,tableLen,1,nvar,1,nvar);
00340       dTabdp=d3tensor(1,tableLen,1,nvar,1,globs->npar);
00341     } 
00342   else 
00343     {
00344       dTabds=NULL; 
00345       dTabdp=NULL;
00346     }
00347   
00348 
00349   // first mesh point
00350   if (doDerivs) {
00351     for (i=1; i<=nvar; ++i) 
00352       {
00353         for (j=1; j<=nvar; ++j)
00354           { 
00355             ex->dyds[1][i][j] = (i==j) ? 1 : 0;
00356             ex->dmds[1][i][j] = (i==j) ? 1 : 0;
00357           }
00358         for (j=1; j<=globs->npar; ++j)  
00359           {
00360             ex->dydp[1][i][j] = 0;
00361             ex->dmdp[1][i][j] = 0;
00362           }
00363       }
00364   }
00365   for (j=1; j<=nvar; ++j) 
00366     ex->yComp[1][j]=ex->yTry[1][j];
00367 
00368 
00369   if (doDerivs && mesh[1]==xMeasure[1]) 
00370     {
00371       observation(globs,ex,mesh[1],ex->yTry[1],ex->yPred[1],ex->par,ex->dmds[1],ex->dmdp[1]);
00372     }
00373   
00374   // integrate
00375   k=1;
00376   ksav=1;
00377   if (mesh[1]==xMeasure[1]) 
00378     {
00379       k=2; 
00380       ksav=2;
00381     }
00382   for (i=1; i< nPoints; ++i) 
00383     {   // given starting values
00384       for (j=1; j<=nvar; ++j) 
00385         ys[j]=ex->yTry[i][j];  
00386       if (xMeasure[k]<mesh[i]) 
00387         {
00388           cerr << "fatal: intODE: xMeasure table corrupted (mesh: "
00389                << mesh[i] << ", xMeasure: " << xMeasure[k] << ")\n";
00390           throw 1;
00391         }
00392       // build table of measuring points
00393       n=1;
00394       while (k <= nMeasure && xMeasure[k]<=mesh[i+1]) 
00395         xTab[n++]=xMeasure[k++];
00396       xTab[n]=mesh[i+1];
00397 #ifdef PRINTMESH
00398       *dbg << "Mesh point " << mesh[i] << ": xTab = (";
00399       for (j=1; j<=n; ++j) 
00400         *dbg << '\t' << xTab[j];
00401       *dbg << ")\n";
00402       dbg->flush();
00403 #endif      
00404       //integration 
00405       tabulateValues(globs,ex,mesh[i],mesh[i+1],xTab,n,ys,yTab, dTabds, dTabdp,ex->dyds[i+1],ex->dydp[i+1]);
00406 
00407 #ifdef PRINTINTEGRATE
00408       *dbg << "Integration returned\t  (";
00409       for (j=1; j<=n; ++j) 
00410         {
00411           for (l=1; l<=nvar; ++l) 
00412             *dbg << " " << yTab[j][l];
00413           *dbg << ",";
00414         }
00415       *dbg << ")\n";
00416       dbg->flush();
00417 #endif      
00418       for (j=1; j<=ex->nvar; ++j)
00419         ex->yComp[i+1][j]=ys[j];   // store value at next mesh point
00420 
00421       for (j=ksav; j<k; ++j)          // store value at measuring points
00422         {
00423           for (l=1; l<=ex->nobs; ++l) 
00424             {
00425               ex->yPred[j][l]=yTab[j-ksav+1][l];
00426               if (doDerivs) 
00427                 {
00428                   for (m=1; m<=ex->nvar; ++m)
00429                     ex->dmds[j][l][m]=dTabds[j-ksav+1][l][m];
00430                   for (m=1; m<=globs->npar; ++m)
00431                     ex->dmdp[j][l][m]=dTabdp[j-ksav+1][l][m];
00432                 }
00433             }
00434         }
00435       ksav=k;
00436     }
00437 
00438       //GnuPlotting
00439   long mind=2;
00440   if(doDerivs && (globs->noGnu==FALSE) && doPlot==TRUE)
00441     { 
00442       for(l=1;l<=ex->nobs;l++)
00443         {
00444           if(globs->gnuindex <= globs->ngnu)
00445             {
00446               mind=2;
00447               ofstream gnuout;
00448               gnuout.open("gnuout.dat");
00449               ymax=ex->yPred[1][l];
00450               ymin=ex->yPred[1][l];           
00451               for(j=1;j<=ex->nMeasure;j++)
00452                 {
00453                   gnuout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
00454                   if(mesh[mind] <= ex->xMeasure[j])
00455                     {
00456                       observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
00457                       gnuout << mesh[mind] << "\t" << ys[l] << endl << endl;
00458                       observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
00459                       gnuout << mesh[mind] << "\t" << ys[l] << endl;
00460                       mind++;
00461                     }
00462                   //determine y-range
00463                   if(ymax <= ex->yPred[j][l])
00464                     ymax=ex->yPred[j][l];
00465                   if(ymin >= ex->yPred[j][l])
00466                     ymin=ex->yPred[j][l];
00467                   if(ymax <= ex->yMeasure[j][l])
00468                     ymax=ex->yMeasure[j][l];
00469                   if(ymin >= ex->yMeasure[j][l])
00470                     ymin=ex->yMeasure[j][l];
00471                 }
00472               gnuout.flush();
00473               gnuout.close();
00474               
00475               fprintf(globs->gnuFp[globs->gnuindex],"reset\n");
00476               fprintf(globs->gnuFp[globs->gnuindex],"set yrange[%f:%f]\n",ymin-0.01*ymin,ymax+0.01*ymax);
00477               fprintf(globs->gnuFp[globs->gnuindex],"set xrange[%f:%f]\n",ex->xMeasure[1],ex->xMeasure[ex->nMeasure]);        
00478               fflush(globs->gnuFp[globs->gnuindex]);
00479               fprintf(globs->gnuFp[globs->gnuindex],"plot \"gnuout.dat\" title \"\"w l 3,\"gnuout.dat\" u 1:3 title \"\" 1 \n");
00480               fflush(globs->gnuFp[globs->gnuindex]);
00481               //DELAY
00482               for(long dd=1;dd<=1000000;dd++)
00483                1+1;
00484             }
00485           sprintf(name,"Gnuout_Exp%d_Obs%d.dat",expNr,l);
00486           outout.open(name);
00487           
00488           globs->gnuindex++;
00489           mind=2;
00490           for(j=1;j<=ex->nMeasure;j++)
00491             {
00492               
00493               outout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
00494               if(mesh[mind] <= ex->xMeasure[j])
00495                 {
00496                   observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
00497                   outout << mesh[mind] << "\t" << ys[l] << endl << endl;
00498                   observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
00499                   outout << mesh[mind] << "\t" << ys[l] << endl;
00500                   mind++;
00501                 }
00502             }
00503           outout.flush();
00504           outout.close();
00505           for(long dd=1;dd<=1000000;dd++)
00506             1+1; 
00507         }
00508     }
00509   
00510   if (doDerivs) {
00511     free_d3tensor(dTabdp,1,tableLen,1,nvar,1,globs->npar);
00512     free_d3tensor(dTabds,1,tableLen,1,nvar,1,nvar);
00513   }
00514   free_dmatrix(yTab,1,tableLen,1,nvar);
00515   free_dvector(xTab,1,tableLen);
00516   free_dvector(ys,1,nvar);
00517 }

double Norm Glob globs,
GlobExp ex,
double ***  dS,
double **  dP,
double  uS
 

Definition at line 29 of file dampIt.cc.

References Glob::doP, Glob::npar, GlobExp::nPoints, Glob::nrExp, and GlobExp::nvar.

Referenced by dampIt().

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

void solvLin Glob globs,
GlobExp ex,
int  computeCovar
 

Definition at line 313 of file solvLin.cc.

References condense(), FALSE, Glob::npar, Glob::nrExp, GlobExp::nvar, and Glob::y0fix.

00314 {
00315   
00316   long i,j,k,l,xind=0,yind=0,nL;
00317   long nvarFit=0;  // # der zu fittenden Anfangswerte
00318   long nLpara=0;   // # lokale Parameter
00319   long nparL=0;    // nLpara*nrExp
00320   long npar=globs->npar;
00321   long nvar=ex[1].nvar;
00322   long nrExp=globs->nrExp;
00323   long nparG=0;
00324   long fitDim=0;
00325   long aDim=0,eDim=0,gDim=0;
00326 
00327   //big system (condensed)
00328   double **Ma,*Ra;
00329   double **Me,*Re;
00330   double **Mg,*Rg;
00331 
00332   //solution
00333   double *dX;
00334   //condense experiments
00335   condense(globs,ex);
00336 
00337   //determine dimension of big system
00338   for(i=1;i<=nvar;i++)
00339     {
00340       if(globs->y0fix[i]!=FALSE)
00341         nvarFit++;
00342     }
00343   nvarFit*=nrExp;
00344   for(i=1;i<=npar;i++)
00345     {
00346       if(globs->doP[i]=='L')
00347         nparL++;
00348       else if(globs->doP[i]!=FALSE)
00349         nparG++;
00350     }
00351   nLpara=nparL;
00352   nparL*=nrExp;
00353   fitDim=nvarFit+nparL+nparG;
00354   for(i=1;i<=nrExp;i++)
00355     {
00356       aDim+=ex[i].nMeasure*ex[i].nobs;
00357       eDim+=ex[i].me;
00358       gDim+=ex[i].mg;
00359     }
00360   if(fitDim==0)
00361     {
00362       cerr << "Nothing to fit.\n";
00363       throw 1;
00364     }
00365 #ifdef PRINTDIMENSION
00366   //print dimension of bigsys;
00367   *dbg << "Dimension of Big-System\n";
00368   *dbg << "-----------------------\n";
00369   *dbg << "Least-squares:\n";
00370   *dbg << "#Rows : " << aDim  << "  #Columns : " << fitDim << endl;
00371   *dbg << "Equality constraints:\n";
00372   *dbg << "#Rows : " << eDim  << "  #Columns : " << fitDim << endl;  
00373   *dbg << "Inequality constraints:\n";
00374   *dbg << "#Rows : " << gDim  << "  #Columns : " << fitDim;  
00375   *dbg << "\n\n";
00376   dbg->flush();
00377 #endif
00378 
00379 
00380   long aDimold=aDim;
00381   long eDimold=eDim;
00382   long fitDimold=fitDim;
00383 
00384   // Beginn Local_Parameter_Constraints
00385   if(globs->faktorLexist)
00386   {
00387         // Erweiterung der Dimensionen um nLpara*nExp
00388         aDim=aDim+nparL;
00389         eDim=eDim+nparL;
00390         fitDim=fitDim+nparL;
00391   }
00392   // Ende Local_Parameter_Constraints
00393   
00394   //allocate memory
00395   Ma=dmatrix(0,aDim,1,fitDim);
00396   Me=dmatrix(0,eDim,1,fitDim);
00397   Mg=dmatrix(0,gDim,1,fitDim);
00398   Ra=dvector(0,aDim);
00399   Re=dvector(0,eDim);
00400   Rg=dvector(0,gDim);
00401 
00402   //initialise 
00403   for(i=1;i<=aDim;i++)
00404     {
00405       Ra[i]=0.;
00406       for(j=1;j<=fitDim;j++)
00407         Ma[i][j]=0.;
00408     }
00409    for(i=1;i<=eDim;i++)
00410     {
00411       Re[i]=0.;
00412       for(j=1;j<=fitDim;j++)
00413         Me[i][j]=0.;
00414     } 
00415    for(i=1;i<=gDim;i++)
00416      {
00417       Rg[i]=0.;
00418       for(j=1;j<=fitDim;j++)
00419         Mg[i][j]=0.;
00420     } 
00421 
00422   //least-squares
00423   for(i=1;i<=nrExp;i++)
00424     {
00425       for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)    
00426         Ra[yind+k]=-ex[i].ua[k];
00427       //initial values
00428       for(j=1;j<=nvar;j++)
00429         {
00430           if(globs->y0fix[j]!=FALSE)
00431             {
00432               xind++;
00433               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)    
00434                 Ma[yind+k][xind]=ex[i].Ea[k][j];
00435             }
00436           
00437         }
00438       //parameters
00439       nL=0;
00440       for(j=1;j<=npar;j++)
00441         {
00442           if(globs->doP[j]=='L')
00443             {
00444               xind++;
00445               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
00446                 Ma[yind+k][xind]=ex[i].Pa[k][j];
00447             }
00448           else if(globs->doP[j]!=FALSE)
00449             {
00450               nL++;
00451               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
00452                 Ma[yind+k][nvarFit+nparL+nL]=ex[i].Pa[k][j];
00453             }
00454         }
00455       yind=yind+(ex[i].nMeasure*ex[i].nobs);
00456     }//end loop over experiments
00457   //-------------> end lest squares
00458   
00459   //equality constraints
00460   yind=0;
00461   xind=0;
00462   for(i=1;i<=nrExp;i++)
00463     {
00464       for(k=1;k<=ex[i].me;k++)    
00465         Re[yind+k]=-ex[i].ue[k];
00466       //initial values
00467       for(j=1;j<=nvar;j++)
00468         {
00469           if(globs->y0fix[j]!=FALSE)
00470             {
00471               xind++;
00472               for(k=1;k<=ex[i].me;k++)    
00473                 Me[yind+k][xind]=ex[i].Ee[k][xind];
00474             }
00475         }
00476       //parameters
00477       nL=0;
00478       for(j=1;j<=npar;j++)
00479         {
00480           if(globs->doP[j]=='L')
00481             {
00482               xind++;
00483               for(k=1;k<=ex[i].me;k++)
00484                 Me[yind+k][xind]=ex[i].Pe[k][j];
00485             }
00486           else if(globs->doP[j]!=FALSE)
00487             {
00488               nL++;
00489               for(k=1;k<=ex[i].me;k++)
00490                 Me[yind+k][nvarFit+nparL+nL]=ex[i].Pe[k][j];
00491             }
00492         }
00493       yind+=ex[i].me;
00494     }//end loop over experiments
00495   //-------------> end equality  constraints
00496 
00497   //inequality constraints
00498   yind=0;
00499   xind=0;
00500   for(i=1;i<=nrExp;i++)
00501     {
00502       for(k=1;k<=ex[i].mg;k++)    
00503         Rg[yind+k]=-ex[i].ug[k];
00504       //initial values
00505       for(j=1;j<=nvar;j++)
00506         {
00507           if(globs->y0fix[j]!=FALSE)
00508             {
00509               xind++;
00510               for(k=1;k<=ex[i].mg;k++)  
00511                 Mg[yind+k][xind]=ex[i].Eg[k][xind];
00512             }
00513         }
00514       //parameters
00515       nL=0;
00516       for(j=1;j<=npar;j++)
00517         {
00518           if(globs->doP[j]=='L')
00519             {
00520               xind++;
00521               for(k=1;k<=ex[i].mg;k++)
00522                 Mg[yind+k][xind]=ex[i].Pg[k][j];
00523             }
00524           else if(globs->doP[j]!=FALSE)
00525             {
00526               nL++;
00527               for(k=1;k<=ex[i].mg;k++)
00528                 Mg[yind+k][nvarFit+nparL+nL]=ex[i].Pg[k][j];
00529             }
00530         }
00531       yind+=ex[i].mg;
00532     }//end loop over experiments
00533   //-------------> end equality  constraints
00534   
00535   //Big System Ready to Solve
00536   //-------------------------
00537   //
00538   //  Ma*dX  ~  Ra
00539   //  Me*dX  =  Re
00540   //  Mg*dX >=  Rg
00541   //
00542   //-------------------------
00543   
00544   //****************************************************************************
00545   // Beginn Local_Parameter_Constraints
00546   // Erweiterung der Systems auf Local_Parameter_Constraints:
00547   //
00548   // pLmean und pLsig bestimmen
00549   if(globs->faktorLexist)
00550   {
00551         double *pLmean, *pLsig;
00552         int index;
00553         int index_local_parameter;
00554         pLmean=dvector(1,nLpara);
00555         pLsig=dvector(1,nLpara);
00556   
00557 //      for(i=1;i<=nLpara;i++)
00558 //      {
00559 //         cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
00560 //         if(globs->faktorL[i]<0)
00561 //         {
00562 //            //i;
00563 //            break;
00564 //         }
00565 //      }
00566 //      cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
00567 //      if(i<=nLpara)
00568 //      {
00569 //         cerr<<"number of local parameter constraints too small"<<endl;
00570 //         throw(1);
00571 //      }
00572         for(i=1;i<=nLpara;i++)
00573         {
00574                 pLmean[i]=0;
00575                 index=0;
00576                 index_local_parameter=0;
00577                 while(index_local_parameter!=i)
00578                 {
00579                         index++;
00580 //                      cerr<<"doP[index]="<<globs->doP[index]<<" index="<<index<<" i="<<i<<endl;
00581                         if(globs->doP[index]=='L')
00582                         {
00583                                 index_local_parameter++;
00584 //                              cerr<<"index="<<index<<" i="<<i<<endl;
00585                         }
00586                         
00587                 }               
00588                 for(j=1;j<=nrExp;j++)
00589                 {
00590                         pLmean[i]=pLmean[i]+ex[j].par[index]; 
00591 //                      cerr<<"ex[j].par[index]="<<ex[j].par[index]<<endl;
00592                 }
00593                 pLmean[i]=pLmean[i]/double(nrExp);
00594           //    cerr<<"pLmean["<<i<<"]="<<pLmean[i]<<" globs->faktorL["<<i<<"]"<<globs->faktorL[i]<<endl;
00595         }
00596         // Berechnung pLsig
00597         for(i=1;i<=nLpara;i++)
00598         {
00599                 pLsig[i]=(globs->faktorL[i]*pLmean[i]-pLmean[i])/(1+globs->faktorL[i]);
00600                 if(pLsig[i]<0)
00601                 {
00602                 cerr<<"ERROR - Negative varianz bei den Lokalen Parameter Constraints..."<<endl;
00603                 }
00604         }       
00605         // Ra,Re
00606         k=1;
00607         for(i=1;i<=nrExp;i++)
00608         {
00609                 for(j=1;j<=nLpara;j++)
00610                 {       
00611                         index=0;
00612                         index_local_parameter=0;
00613                         while(index_local_parameter!=j)
00614                         {
00615                                 index++;
00616                                 if(globs->doP[index]=='L')
00617                                 {
00618                                         index_local_parameter++;
00619                                 }
00620                         
00621                         }  
00622                         Ra[aDimold+k]=(ex[i].par[index]-pLmean[j])/pLsig[j];
00623                         Re[eDimold+k]=0;
00624                         double pLsig_real=0;
00625                         for(int z=1;z<=nrExp;z++)
00626                         {
00627                                 pLsig_real=pLsig_real+(ex[z].par[index]-pLmean[j])*(ex[z].par[index]-pLmean[j]);
00628                         }
00629                         pLsig_real=sqrt(pLsig_real/nrExp);
00630                         //*dbg<<"ex["<<i<<"].par["<<index<<"]="<<ex[i].par[index]<<" pLmean["<<j<<"]="<<pLmean[j]<<" pLsig["<<j<<"]="<<pLsig[j]<<" pLsig_real="<<pLsig_real<<endl;
00631                         // Ma
00632                         for(int l=1;l<=nrExp;l++)
00633                         {
00634                                 Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=1/(pLsig[j]*nrExp);
00635                                 if(i == l)
00636                                 {
00637                                         Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]-1/pLsig[j];
00638                                 }
00639                         }
00640                         k++;
00641                 }
00642         }
00643         // Me
00644         for(i=eDimold+1;i<=eDim;i++)
00645         {
00646                 for(j=nvarFit+1;j<=fitDim;j++)   
00647                 {
00648                         if(j==fitDimold+i)
00649                         {
00650                                 Me[i][j]=Me[i][j]+1;
00651                         }
00652                         if(j==nvarFit+i)
00653                         {
00654                                 Me[i][j]=Me[i][j]-1;
00655                         }
00656                 }
00657         }
00658   }
00659   // Ende Local_Parameter_Constraints
00660   //****************************************************************************
00661   
00662   //regularising
00663   if(globs->reg==TRUE)
00664     {
00665       double *diag,**V,**MMa;
00666       double cond,maxx;
00667       double reg=0.;
00668 
00669       V=dmatrix(1,fitDim,1,fitDim);
00670       MMa=dmatrix(1,aDim,1,fitDim);
00671       diag=dvector(1,fitDim);
00672 
00673        for(i=1;i<=aDim;i++)
00674         {
00675           for(j=1;j<=fitDim;j++)
00676             {
00677               MMa[i][j]=Ma[i][j];
00678             }     
00679         }
00680       
00681       //SVD
00682       svdcmp(MMa,aDim,fitDim,diag,V);
00683       
00684       //determine condition number
00685       *dbg << "\ncondition number : ";
00686       maxx=diag[1];
00687       cond=diag[1];
00688       for(i=1;i<=fitDim;i++)
00689         {
00690           if(diag[i]>maxx)
00691             maxx=diag[i];
00692           if(diag[i]<cond)
00693             cond=diag[i];
00694         }
00695 
00696       //Is the matrix approximately singular?
00697       for(i=1;i<=fitDim;i++)
00698         {
00699           if(fabs(diag[i]/maxx)<=globs->epsilon)
00700             {
00701               diag[i]=0.;
00702               reg=globs->lambda*maxx;
00703             }
00704         }
00705       
00706       cond*=1./maxx;
00707       *dbg << cond << endl;
00708       globs->cond=cond;
00709       
00710       //writing regularised system back
00711       for(i=1;i<=aDim;i++)
00712         {
00713           for(j=1;j<=fitDim;j++)
00714             {
00715               Ma[i][j]=0.;
00716               for(k=1;k<=fitDim;k++)
00717                 {
00718                   Ma[i][j]+=MMa[i][k]*diag[k]*V[j][k];
00719                 }
00720               if(i <=fitDim)
00721                 {
00722                   if(diag[i]==0.)
00723                     Ma[i][j]+=reg*V[j][i];
00724                 }
00725             }
00726         }
00727       if(reg!=0.)
00728         {
00729           *dbg << "regularisation active, linear dependency of parameters:\n";
00730           for(i=1;i<=fitDim;i++)
00731               {
00732                 if(diag[i]==0.)
00733                   *dbg << "line " << i << ": ";
00734                 for(j=1;j<=fitDim;j++)
00735                   {
00736                     if(diag[i]==0.)
00737                       *dbg << V[j][i] << " ";
00738                   }
00739                 if(diag[i]==0.)
00740                   *dbg << endl;
00741               }
00742         }
00743 
00744       free_dmatrix(MMa,1,aDim,1,fitDim);
00745       free_dmatrix(V,1,fitDim,1,fitDim);
00746       free_dvector(diag,1,fitDim);
00747     }//end if(globs->reg==TRUE)
00748 
00749 #ifdef PRINTBIGSYS
00750   *dbg << "\nLeast-squares :\n";
00751    for(i=1;i<=aDim;i++)
00752     {
00753       for(j=1;j<=fitDim;j++)
00754         {
00755           *dbg << Ma[i][j] << " ";
00756         }
00757       *dbg << "| " << Ra[i] << endl;
00758     }  
00759    *dbg << "\nEquality constraints :\n";
00760    for(i=1;i<=eDim;i++)
00761     {
00762       for(j=1;j<=fitDim;j++)
00763         {
00764           *dbg << Me[i][j] << " ";
00765         }
00766       *dbg << "| " << Re[i] << endl;
00767     }
00768    *dbg << "\nInquality constraints :\n";
00769    for(i=1;i<=gDim;i++)
00770     {
00771       for(j=1;j<=fitDim;j++)
00772         {
00773           *dbg << Mg[i][j] << " ";
00774         }
00775       *dbg << "| " << Rg[i] << endl;
00776     }
00777    dbg->flush();
00778   #endif
00779   //allocate solution vector
00780   dX=(double*)malloc((fitDim+10)*sizeof(double));
00781   //allocate covariance matrix
00782   globs->fitdim=fitDim;
00783   if(globs->covar==NULL)
00784     globs->covar=dmatrix(1,fitDim,1,fitDim);
00785 
00786   if(globs->minimiser==1)
00787     {
00788       // ###############
00789       // #     LSEI    #
00790       // ###############
00791       
00792       //allocation/initialisation for LSEI
00793       double *w;
00794       double *ws;
00795       double *prgopt;
00796       double rnorme,rnorml;
00797       long mdw=aDim+eDim+gDim;
00798       long mode;
00799       long ipstor=mdw+npar+2*fitDim+2;
00800       long wsstor=3*(mdw+npar)+(mdw+npar+2)*(fitDim+8)+fitDim+npar;
00801       long *ip;
00802       
00803       w=(double*)malloc(((fitDim+1)*mdw)*sizeof(double));
00804       ws=(double*)malloc((wsstor+1)*sizeof(double));
00805       ip=(long*)malloc((ipstor+10)*sizeof(long));
00806       prgopt=(double*)malloc(30*sizeof(double));
00807       
00808       xind=0;
00809       yind=0;
00810       
00811       for(i=1;i<=fitDim;i++)
00812         {
00813           for(j=1;j<=eDim;j++)
00814             {
00815               w[xind*mdw+yind]=Me[j][i];
00816               yind++;
00817             }
00818           for(j=1;j<=aDim;j++)
00819             {
00820               w[xind*mdw+yind]=Ma[j][i];
00821               yind++;
00822             }
00823           for(j=1;j<=gDim;j++)
00824             {
00825               w[xind*mdw+yind]=Mg[j][i];
00826               yind++;
00827             } 
00828           yind=0;
00829           xind++;
00830         }
00831  
00832       for(j=1;j<=eDim;j++)
00833         {
00834           w[(fitDim)*mdw+yind]=Re[j];
00835           yind++;
00836         }
00837       for(j=1;j<=aDim;j++)
00838         {
00839           w[(fitDim)*mdw+yind]=Ra[j];
00840           yind++;
00841         }
00842       for(j=1;j<=gDim;j++)
00843         {
00844           w[(fitDim)*mdw+yind]=Rg[j];
00845           yind++;
00846         }
00847       
00848       ip[0]=wsstor;
00849       ip[1]=ipstor;
00850       
00851       if(computeCovar==TRUE)
00852         {
00853           prgopt[0]=4.0; 
00854           prgopt[1]=1.0; 
00855           prgopt[2]=1.0;
00856           prgopt[3]=7.0;
00857           prgopt[4]=10.0;
00858           prgopt[5]=1.0;
00859         }
00860       else
00861         {
00862           prgopt[0]=4.0; 
00863           prgopt[1]=0.0; 
00864           prgopt[2]=1.0;
00865           prgopt[3]=7.0;
00866           prgopt[4]=10.0;
00867           prgopt[5]=0.0;
00868         }
00869       prgopt[6]=1.0;
00870       
00871       //call LSEI
00872       lsei_(w,mdw,eDim,aDim,gDim,fitDim,prgopt,dX,&rnorme,&rnorml,&mode,ws,ip);
00873       
00874       if(mode>=4)
00875         {
00876           cerr << "Minimisation (LSEI) failed.\n";
00877           throw 1;
00878         }
00879       
00880       if(computeCovar==TRUE)
00881         {
00882           for (i=1; i<=fitDim; ++i)
00883             {
00884               for (j=1;j<=fitDim;++j)
00885                 globs->covar[i][j] = w[(j-1)*mdw+i-1];
00886             } 
00887         }
00888       
00889       free(w);
00890       free(ws);
00891       free(ip);
00892       free(prgopt);
00893       
00894       // #######  END OF LSEI ######
00895     }
00896   else if(globs->minimiser==2)
00897     {
00898       //to solve ----->
00899       //-------------------------
00900       //
00901       //  Ma*dX  ~  Ra
00902       //  Me*dX  =  Re
00903       //  Mg*dX >=  Rg
00904       //
00905       //-------------------------
00906 
00907 
00908       //initialisation of  E04NCF
00909       long NCLIN=eDim+gDim;
00910       double *C=(double*)malloc(((eDim+gDim+1)*fitDim)*sizeof(double));
00911       double *BL=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
00912       double *BU=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
00913       double *CVEC=(double*)malloc((fitDim+1)*sizeof(double));
00914       long *ISTATE=(long*)malloc((fitDim+eDim+gDim+1)*sizeof(long));
00915       long *KX=(long*)malloc((fitDim+1)*sizeof(long));
00916       double *X=(double*)malloc((1+fitDim)*sizeof(double));
00917       double *A=(double*)malloc((1+fitDim*aDim)*sizeof(double));
00918       double *B=(double*)malloc((1+aDim)*sizeof(double));
00919       long ITER;
00920       double OBJ;
00921       double *CLAMDA=(double*)malloc((1+fitDim+eDim+gDim)*sizeof(double));
00922       long *IWORK=(long*)malloc((1+fitDim)*sizeof(long));
00923       long LWORK=2*fitDim*fitDim+9*fitDim+6*(eDim+gDim);
00924       double *WORK=(double*)malloc((1+LWORK)*sizeof(double));
00925       long IFAIL=-1;
00926       double **CTMP=dmatrix(1,1+eDim+gDim,1,1+fitDim);
00927 
00928       //right-hand side
00929       for(i=1;i<=aDim;i++)
00930         B[i-1]=Ra[i];
00931       
00932       //bounds
00933       for(i=1;i<=fitDim;i++)
00934         {
00935           BL[i-1]=-1e25;
00936           BU[i-1]=1e25;
00937           X[i-1]=0.;  //initialisation of X
00938         }
00939       for(i=fitDim+1;i<=fitDim+eDim;i++)
00940         {
00941           BL[i-1]=Re[i-fitDim];
00942           BU[i-1]=Re[i-fitDim];
00943         }
00944       for(i=fitDim+eDim+1;i<=fitDim+eDim+gDim;i++)
00945         {
00946           BL[i-1]=Rg[i-fitDim-eDim];
00947           BU[i-1]=1e25;
00948         }
00949       
00950       for(i=1;i<=fitDim;i++)
00951         {
00952           for(j=1;j<=eDim+gDim;j++)
00953             {
00954               if(j<=eDim)
00955                 CTMP[j][i]=Me[j][i];
00956               else
00957                 CTMP[j][i]=Mg[j-eDim][i];
00958             }
00959         }
00960       
00961       for(i=1;i<=fitDim+eDim+gDim;i++)
00962       {
00963         ISTATE[i]=0;
00964       }
00965       
00966       cmat2fvec(Ma,aDim,fitDim,A);
00967       cmat2fvec(CTMP,eDim+gDim,fitDim,C);
00968       
00969       setopte04ncf_();
00970       e04ncf_(aDim,fitDim,NCLIN,NCLIN,aDim,C,BL,BU,CVEC,ISTATE,
00971               KX,X,A,B,ITER,OBJ,CLAMDA,IWORK,fitDim,WORK,LWORK,IFAIL);
00972       for(i=0;i< fitDim;i++)
00973       {
00974         dX[i]=X[i];
00975       }
00976       
00977       //free memory
00978       free(C);
00979       free(BL);
00980       free(BU);
00981       free(CVEC);
00982       free(ISTATE);
00983       free(KX);
00984       free(X);
00985       free(A);
00986       free(B);
00987       free(CLAMDA);
00988       free(IWORK);
00989       free(WORK);
00990       free_dmatrix(CTMP,1,eDim+gDim,1,fitDim);
00991     }
00992 
00993   //decondense experiments
00994   decondense(globs,ex,dX);
00995 
00996 #ifdef PRINTSTEP
00997   for(i=1;i<=nrExp;i++)
00998     {
00999       *dbg << "\nExperiment #" << i << ":\n\n";
01000       for(k=1;k<=ex[i].nPoints-1;k++)
01001         {
01002           *dbg << "dS(msP=" << k << ") = ";
01003           for(j=1;j<=nvar;j++)
01004             *dbg << ex[i].dS[k][j] << " ";
01005           *dbg << endl;
01006         }
01007       *dbg << "dP = ";
01008       for(j=1;j<=npar;j++)
01009         *dbg << ex[i].dP[j] << " ";
01010     }
01011   *dbg << endl;
01012   dbg->flush();
01013 #endif
01014   
01015 
01016   //free memory
01017   free_dmatrix(Ma,0,aDim,1,fitDim);
01018   free_dmatrix(Me,0,eDim,1,fitDim);
01019   free_dmatrix(Mg,0,gDim,1,fitDim);
01020   free_dvector(Ra,0,aDim);
01021   free_dvector(Re,0,eDim);
01022   free_dvector(Rg,0,gDim);
01023 
01024   dX[0];
01025   free(dX);
01026 }

double subIt Glob globs,
GlobExp ex,
double  t,
double  normdx,
double **  p0,
double ***  s0,
double **  dp,
double ***  ds,
double  uS
 

Definition at line 62 of file dampIt.cc.

References Glob::npar, and Glob::nrExp.

Referenced by dampIt().

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


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