/home/peifer/diffit/modules/initialise.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 //Begin definition of module prototypes
00014 
00015 void setInitialValues(GlobExp *ex, int nP, double *parameters, double *yValues);
00016 
00017 //End definition of module prototypes
00018 
00019 //Tolerance of unsatisfied equality constraints
00020 #define EPSTOL 1e-8
00021 
00022 //Checking if initial values are compartible with constrains
00023 int initValuesOK(Glob *globs,GlobExp *ex)
00024 {
00025   int i,good=TRUE;
00026   
00027   if (ex->me > 0) 
00028     {
00029       R2(globs,ex, FALSE);
00030       for (i=1; i<=ex->me; ++i)
00031         if (good) 
00032           good = fabs(ex->r2[i]) < EPSTOL;
00033     }
00034   
00035   if (ex->mg>0) 
00036     {
00037       R3(globs,ex, FALSE);
00038       for (i=1; i<=ex->mg; ++i) 
00039         good = good && (ex->r3[i]>=0);
00040     }
00041   
00042   return(good);
00043 }
00044 
00045 //Gnuplotting
00046 void setupGnuFp(Glob *globs, GlobExp ex[])
00047 {
00048   long gnu_n=20,ncomp,k;
00049   int fullx=700; //max. geometry
00050   int fully=500;
00051   int nwi,dx=fullx/2,dy,xpos,ypos;
00052   char outstr[100];
00053 
00054 
00055   ncomp=0;
00056   for(k=1;k<=globs->nrExp;k++)
00057     {
00058       ncomp+=ex[k].nobs;
00059     }
00060 
00061   // determine nr of windows in y-direction
00062   gnu_n = LMIN (gnu_n, ncomp);
00063   globs->ngnu=gnu_n;
00064 
00065   // allocate streams
00066   globs->gnuFp = new FILE *[gnu_n+1];
00067   
00068   if((gnu_n % 2)==1)
00069     nwi=(gnu_n+1)/2;
00070   else
00071     nwi=gnu_n/2;
00072   dy=fully/nwi;
00073  
00074   for(k=1;k<=gnu_n;k++)
00075     {
00076       if(k==1)
00077         {
00078           xpos=0;
00079           ypos=0;
00080         }
00081       else if((k % 2)==0)
00082         {
00083           xpos=dx+10;
00084         }
00085       else
00086         {
00087           xpos=0;
00088           ypos+=dy+25;
00089         }
00090       sprintf(outstr,"gnuplot -noraise -geometry %ix%i+%i+%i -title \"%i\"\0",dx,dy,xpos,ypos,k);
00091       globs->gnuFp[k]=popen(outstr,"w");
00092     }
00093 }                               // setupGnuFp
00094 
00095 
00096 //Initialise all vectors, matrices, tensors, etc. 
00097 //used during the iterations
00098 void initialise(GlobExp ex[],Glob *globs,int simit)
00099 {
00100   long i,j,k,nExp;
00101   long nrExp=globs->nrExp;
00102   
00103   ifstream inSplines;
00104 
00105   if(simit==FALSE)
00106     {
00107       for (nExp=1;nExp<=nrExp;++nExp) 
00108         { 
00109           
00110           //non-local parameters are equal for each experiment
00111           for(i=1;i<=globs->npar;i++)
00112             {
00113               if(globs->doP[i]==TRUE)
00114                 ex[nExp].par[i]=ex[1].par[i];
00115             }
00116           
00117           // set ex->me and ex->mg given in model.cc
00118           setNrConstraints(&ex[nExp],globs->npar,ex[nExp].par);
00119           
00120           long nPoints=ex[nExp].nPoints;
00121           long nP=globs->npar;
00122           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
00123           long nMeas=ex[nExp].nMeasure;
00124           long  me=ex[nExp].me, mg=ex[nExp].mg;
00125           
00126 #ifdef PRINTINITVALUES
00127           *dbg << "initialize: exp. #" << nExp << ", nPoints=" << nPoints
00128                << ", nvar=" << nvar << ", nobs=" << nobs << "\nnMeas="
00129                << nMeas << ", nVal=" << nVal << ", me=" << me
00130                << ", mg=" << mg << ", nP=" << nP << "\n";
00131           dbg->flush();
00132 #endif
00133           //std. error -> output
00134           ex[nExp].errP=dvector(1,nP);
00135           ex[nExp].errY0=dvector(1,nvar);         
00136           // initial guesses at mesh points
00137           ex[nExp].yTry=dmatrix(1,nPoints,1,nvar); 
00138           ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar); 
00139           // computed values at mesh points   
00140           ex[nExp].yComp=dmatrix(1,nPoints,1,nvar);
00141           // computed values at measuring points
00142           ex[nExp].yPred=dmatrix(1,nMeas,1,nobs);
00143           
00144           ex[nExp].h=dmatrix(1,nPoints,1,nvar);      // discrepancies
00145           ex[nExp].residues=dmatrix(1,nMeas,1,nobs); // residues
00146           
00147           // derivatives
00148           // d(yComp[i][j]) / d(yTry[i-1][k])
00149           ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar);
00150           // d(yComp[i][j]) / d(p[k])
00151           ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP);
00152           // d(yPred[i][j]) / d(yTry[**][k])
00153           ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar);
00154           // d(yPred[i][j]) / d(p[k])
00155           ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP);
00156           if (me>0) 
00157             {
00158               ex[nExp].r2=dvector(1,me);                 // equal. constr.
00159               // d(yR2[i]) / d(yTry[j][k])
00160               ex[nExp].dR2ds=d3tensor(1,me,1,nPoints,1,nvar);
00161               for(i=1;i<=me;i++)
00162                 for(j=1;j<=nPoints;j++)
00163                   for(k=1;k<=nvar;k++)
00164                     ex[nExp].dR2ds[i][j][k]=0.;
00165               
00166               // d(yR2[i]) / d(p[j])
00167               ex[nExp].dR2dp=dmatrix(1,me,1,nP);
00168               for(i=1;i<=me;i++)
00169                 for(j=1;j<=nP;j++)
00170                   ex[nExp].dR2dp[i][j]=0.;
00171             } 
00172           else 
00173             {
00174               ex[nExp].r2=NULL; ex[nExp].dR2ds=NULL; ex[nExp].dR2dp=NULL;
00175             }
00176           if (mg>0) 
00177             {
00178               ex[nExp].r3=dvector(1,mg);                 // ineq.  constr.
00179               // d(yR3[i]) / d(yTry[j][k])
00180               ex[nExp].dR3ds=d3tensor(1,mg,1,nPoints,1,nvar);
00181               for(i=1;i<=mg;i++)
00182                 for(j=1;j<=nPoints;j++)
00183                   for(k=1;k<=nvar;k++)
00184                     ex[nExp].dR3ds[i][j][k]=0.;
00185               // d(yR3[i]) / d(p[j])
00186               ex[nExp].dR3dp=dmatrix(1,mg,1,nP);
00187               for(i=1;i<=mg;i++)
00188                 for(j=1;j<=nP;j++)
00189                   ex[nExp].dR3dp[i][j]=0.;
00190             } 
00191           else 
00192             {
00193               ex[nExp].r3=NULL; ex[nExp].dR3ds=NULL; ex[nExp].dR3dp=NULL;
00194             }
00195           
00196           // compute yTry from ex (all ex read from file still present)
00197           setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0);
00198           // check if initial values are compatible with constraints
00199           // CURRENTLY DISABLED (1. Dez. 2004)
00200 //        if (!initValuesOK(globs,&ex[nExp])) 
00201 //          {
00202 //            cerr << "Experiment #" << nExp 
00203 //                 << ": initial values are not compatible with constraints\n";
00204 //            cerr << "R2:";
00205 //            for (i=1;i<=ex[nExp].me;++i) 
00206 //              cerr << "\t" << ex[nExp].r2[i];
00207 //            cerr << "\nR3:";
00208 //            for (i=1;i<=ex[nExp].mg;++i) 
00209 //              cerr << "\t" << ex[nExp].r3[i];
00210 //            cerr << "\n";
00211 //            exit(1);
00212 //          }
00213           
00214 #ifdef PRINTEX
00215           *dbg << "Ex used for fitting: \n\n";
00216           for (i=1; i<=ex[nExp].nMeasure; ++i) 
00217             {
00218               *dbg << ex[nExp].xMeasure[i];
00219               for (j=1; j<=ex[nExp].nobs; ++j)
00220                 *dbg << "\t" << ex[nExp].yMeasure[i][j] << "\t" << ex[nExp].sigma[i][j];
00221               *dbg << '\n';
00222             }
00223           *dbg << '\n';
00224           dbg->flush();
00225 #endif
00226           //Allocate memory for condensation
00227           //least squares
00228           ex[nExp].ua=dvector(1,nMeas*nobs);
00229           ex[nExp].Ea=dmatrix(1,nMeas*nobs,1,nvar);
00230           ex[nExp].Pa=dmatrix(1,nMeas*nobs,1,nP);
00231           //equality constraints
00232           ex[nExp].ue=dvector(1,me);
00233           ex[nExp].Ee=dmatrix(1,me,1,nvar);
00234           ex[nExp].Pe=dmatrix(1,me,1,nP);
00235           //inequality constraints
00236           ex[nExp].ug=dvector(1,mg);
00237           ex[nExp].Eg=dmatrix(1,mg,1,nvar);
00238           ex[nExp].Pg=dmatrix(1,mg,1,nP);
00239           
00240           //update steps
00241           ex[nExp].dS=dmatrix(1,nPoints,1,nvar);
00242           ex[nExp].dP=dvector(1,nP);
00243         }// matches: for(nExp=1;nExp<=nrExp;++nExp) loop over experiments
00244       //covar in solvLin allocated
00245       globs->covar=NULL;
00246     } // matches with if(simit==FALSE)...
00247 
00248   //Initialise GNUplot
00249   if(globs->noGnu==FALSE)
00250     setupGnuFp(globs,ex);
00251   
00252   globs->cond=0;
00253   globs->Lambda=1;
00254   
00255   char line[1000];
00256   //Initialise Splines
00257   if(globs->initSpline==TRUE)
00258     {
00259       for (nExp=1;nExp<=nrExp;++nExp) 
00260         {
00261           ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00262           ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00263           ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
00264           ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
00265           
00266           for(i=1;i<=NSPLINES;i++)
00267             {
00268               j=0;
00269               while(ex[nExp].splineFile[i][j]!='\0')
00270                 {
00271                   line[j]=ex[nExp].splineFile[i][j];
00272                   j++;
00273                 }
00274               line[j]='\0';
00275               inSplines.open(line);
00276               if(inSplines==NULL)
00277                 {
00278                   cerr << "Cannot open " << ex[nExp].splineFile[i] << ".\n";
00279                   exit(1);
00280                 }
00281               inSplines >>  ex[nExp].nNodes[i];
00282               ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00283               ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00284               ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double));
00285               for(j=1;j<= ex[nExp].nNodes[i];j++)
00286                 {
00287                   inSplines >>  ex[nExp].splineNodes[i][j];
00288                   inSplines >>  ex[nExp].splineY[i][j];
00289                   inSplines >>  ex[nExp].splineGam[i][j];             
00290                 }
00291               inSplines.close();
00292             }
00293         }
00294     }
00295 }

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