- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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