/home/peifer/diffit/modules/initialise.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.

Defines

#define EPSTOL   1e-8

Functions

void setInitialValues (GlobExp *ex, int nP, double *parameters, double *yValues)
int initValuesOK (Glob *globs, GlobExp *ex)
void setupGnuFp (Glob *globs, GlobExp ex[])
void initialise (GlobExp ex[], Glob *globs, int simit)


Define Documentation

#define EPSTOL   1e-8
 

Definition at line 20 of file initialise.cc.

Referenced by initValuesOK().


Function Documentation

void initialise GlobExp  ex[],
Glob globs,
int  simit
 

Definition at line 98 of file initialise.cc.

References d3tensor(), dbg, dmatrix(), GlobExp::dmdp, GlobExp::dmds, Glob::doP, GlobExp::dR2dp, GlobExp::dR2ds, GlobExp::dR3dp, GlobExp::dR3ds, dvector(), GlobExp::dydp, GlobExp::dyds, GlobExp::errP, GlobExp::errY0, FALSE, GlobExp::h, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, GlobExp::par, GlobExp::r2, GlobExp::residues, TRUE, GlobExp::yComp, GlobExp::yPred, GlobExp::yTry, and GlobExp::yTrySave.

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 }

int initValuesOK Glob globs,
GlobExp ex
 

Definition at line 23 of file initialise.cc.

References EPSTOL, FALSE, GlobExp::me, GlobExp::mg, GlobExp::r2, GlobExp::r3, and TRUE.

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 }

void setInitialValues GlobExp ex,
int  nP,
double *  parameters,
double *  yValues
 

Definition at line 8 of file setInitialValues.cc.

References GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, GlobExp::nPoints, GlobExp::nvar, GlobExp::xMeasure, GlobExp::yMeasure, and GlobExp::yTry.

00009 {
00010   long i,j,k;
00011   double xi, d, w1, w2;
00012   long nobs=ex->nobs, nvar=ex->nvar, nMeasure=ex->nMeasure;
00013   long nPoints=ex->nPoints;
00014   
00015   /* set initial values */
00016   if (yValues) 
00017     {
00018       for (j=1; j<=nvar; ++j) 
00019         ex->yTry[1][j]=yValues[j];
00020       i=2;
00021     } 
00022   else
00023     i=1;
00024   k=1;
00025   
00026   while (i<=nPoints) 
00027     {
00028       xi=ex->mesh[i];
00029       // make educated guess for yTry(xi)
00030       
00031       // could be improved (polynomial interpolation/extrapolation etc.)
00032       
00033       if (ex->xMeasure[k]>=xi) 
00034         {
00035           // before first measurement: take first measured value
00036           for (j=1; j<=nvar; ++j)
00037             if (j<=nobs)
00038               ex->yTry[i][j]=ex->yMeasure[k][j];
00039             else
00040               ex->yTry[i][j]=1;
00041         } else if (k==nMeasure) 
00042           {
00043             // after last measurement: take last measured value
00044             for (j=1; j<=nvar; ++j)
00045               if (j<=nobs)
00046                 ex->yTry[i][j]=ex->yMeasure[k][j];
00047               else
00048                 ex->yTry[i][j]=1;
00049           } 
00050       else 
00051         {
00052           // interpolate linearly
00053           while (k<nMeasure && ex->xMeasure[k]<xi) 
00054             ++k;
00055           d=ex->xMeasure[k]-ex->xMeasure[k-1];
00056           w1=(ex->xMeasure[k]-xi)/d;
00057           w2=(xi-ex->xMeasure[k-1])/d;
00058           for (j=1; j<=nvar; ++j) 
00059             if (j<=nobs)
00060               ex->yTry[i][j]=w1*ex->yMeasure[k-1][j] + w2*ex->yMeasure[k][j];
00061             else
00062               ex->yTry[i][j]=1;
00063         }
00064       ++i;
00065     }  
00066 }

void setupGnuFp Glob globs,
GlobExp  ex[]
 

Definition at line 46 of file initialise.cc.

References Glob::nrExp.

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


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