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

Root program for diffit (contains function main). More...

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include "def.h"
#include "model.h"
#include "nr.h"

Go to the source code of this file.

Functions

GlobExpparseopts (int argc, char *argv[], Glob *globs, char *outstr)
void readData (GlobExp *ex, Glob *globs, long expNr)
void setMesh (GlobExp *ex, Glob *globs, long expNr)
void outFit (GlobExp ex[], Glob *globs)
void freeMem (GlobExp *ex, Glob *globs, int simit)
void initialise (GlobExp ex[], Glob *globs, int simit)
void simInit (GlobExp *ex, Glob *globs)
void fitIt (GlobExp ex[], Glob *globs)
void globOpt (GlobExp ex[], Glob *globs)
 main (int argc, char *argv[])

Variables

ofstream * dbg
 debug stream (C++)


Detailed Description

Root program for diffit (contains function main).

Diffit is has a modular structure. This file can be considered as root module for the entire C stand-alone version. A matlab interface is also available under: diffit_mex.cc. The structure of the matlab interface is similar than this file, apart from parsing the data and the program parameters. Note that all global variables are capsuled into two C structures named Glob and GlobExp. The first structure Glob stores all experiment unspecific parameters such as the used integrator, or the integration accuracy, etc. Experiment specific program variables are held in GlobExp, these are, e.g., the number of multiple shooting intervals. Modules called by function main is in order of occurrence:

  1. parseopts(): Parsing all program parameters either from the commant line or a file.
  2. readData(): Reading the data.
  3. setMesh(): This module sets up the multiple shooting intervals for each observation and experiment.
  4. initialise(): Allocates memory and initializes vectors and matrices.
  5. simInit(): If some state variables are not directly obseved (non trivial observation function), this function integrates a trajectory using the initial guess to set up the initial guess for the state variables at multiple shooting intervals. This is activated by the command line argument "-siminit". If not used, the state variables are initialized with 1. Since the integrated trajectory for the initial guess is continuous, rather than discontinuous the "-pert <value>" perturbs this guess using mean zero Gaussian noise. The standard deviation of this noise can be chosen by the value of the -pert argument.
  6. globOpt(): Experimental !! Global optimizer instead of multiple shooting. Currently not working properly.
  7. fitIt(): General fitting function, all the numerics is done in here.
  8. outFit(): This routine writes the output such as parameter estimates and its standard error.
  9. freeMem(): Deallocation of vectors and matrices.

Definition in file diffit.cc.


Function Documentation

void fitIt GlobExp  ex[],
Glob globs
 

Definition at line 33 of file fitIt.cc.

00034 {
00035   long nExp,i,j,k,l,n;
00036   long tryOnceMore,ok;
00037   double objF, sum, scale, lambda;
00038   long maxPoints=0,maxVar=0;
00039   double minImprovement=globs->minimp;
00040 
00041   double uS=1.;
00042 
00043   for (i=1;i<=globs->nrExp;++i) 
00044     {
00045       if (ex[i].nPoints>maxPoints) 
00046         maxPoints=ex[i].nPoints;
00047       if (ex[i].nvar>maxVar)       
00048         maxVar=ex[i].nvar;
00049     }
00050   
00051   double **dPSave=dmatrix(1,globs->nrExp,1,globs->npar);      
00052   double ***dSSave=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00053   double **P0=dmatrix(1,globs->nrExp,1,globs->npar);
00054   double ***S0=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00055   double **dP=dmatrix(1,globs->nrExp,1,globs->npar);      
00056   double ***dS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar);
00057   
00058   *dbg << "Minimum norm for Gauss-Newton update: minImprovement = " << minImprovement << "\n\n";
00059   
00060   // GGN iteration
00061   // do the actual fitting
00062   
00063   globs->nIter=1; 
00064 
00065   tryOnceMore=TRUE;
00066   while (globs->nIter < globs->maxit && tryOnceMore) 
00067     {
00068       globs->gnuindex=1;
00069       if(globs->silent!=TRUE)
00070         {
00071           cout << "Integration # " << globs->nIter << '\n';
00072           for(j=1;j<=globs->nrExp;j++)
00073             {
00074               cout << "Parameter value(s) of experiment " << j << ":\n";
00075               for (i=1; i<=globs->npar; ++i) 
00076                 {
00077                   cout << ex[j].par[i] << " ";
00078                 }
00079               cout << endl;
00080               cout << "Initial value(s) of experiment " << j << ":\n";
00081               for (i=1; i<=ex[j].nvar; ++i) 
00082                 {
00083                   cout << ex[j].yTry[1][i] << " ";
00084                 }
00085               cout << endl;
00086             }
00087           cout << "\n";
00088           cout.flush();
00089         }
00090       *dbg << "Integration # " << globs->nIter << '\n';
00091       *dbg << "Parameter value(s):";
00092       for (i=1; i<=globs->npar; ++i) 
00093         {
00094           if(globs->doP[i]!='L')
00095             *dbg << '\t' << ex[1].par[i];
00096         }
00097       *dbg << "\n\n";
00098       
00099       // all experiments
00100       sum=0;
00101       for (nExp=1;nExp<=globs->nrExp;++nExp) 
00102         {
00103 #ifdef PRINTINITVALUES
00104           *dbg << "Experiment #" << nExp << "\n";
00105           *dbg << "Initial values at mesh points:\n";
00106           for (i=1; i<ex[nExp].nPoints; ++i) 
00107             {
00108               *dbg << ex[nExp].mesh[i] << ':';
00109               for (j=1; j<=ex[nExp].nvar; ++j) {
00110                 *dbg << '\t' << ex[nExp].yTry[i][j];
00111               }
00112               *dbg << '\n';
00113             }
00114           dbg->flush();
00115 #endif    
00116           intODE(&ex[nExp],globs,TRUE,TRUE,nExp);
00117           
00118           // set up internal data (residues, etc.)
00119           objF=computeRight(globs,&ex[nExp]);
00120           sum += objF;
00121           
00122           *dbg << "Experiment #" << nExp << ": objective f ="
00123                << objF << "\n";
00124           *dbg << "Discrepancies: ";
00125           for (i=1; i<ex[nExp].nPoints; ++i) 
00126             {
00127               for (j=1; j<=ex[nExp].nvar; ++j)
00128                 *dbg << ' ' << ex[nExp].h[i][j];
00129               *dbg << ",";
00130             }
00131           *dbg << "\n\n";
00132           dbg->flush();
00133         }
00134 
00135       if(globs->silent!=TRUE)
00136         cout << "Objective function f = " << sum << "\n";
00137       *dbg << "Global objective function f = " << sum << "\n";
00138       
00139           for (nExp=1;nExp<=globs->nrExp;++nExp) 
00140             {
00141               for (i=1; i<ex[nExp].nPoints; ++i)
00142                 {
00143                   for (j=1; j<=ex[nExp].nvar; ++j) 
00144                     {
00145                       ex[nExp].yTrySave[i][j] = ex[nExp].yTry[i][j];
00146                     }
00147                 }
00148             }
00149       // solve linearized minimisation problem
00150       solvLin(globs,ex,FALSE);
00151 
00152       //save old values and new updates
00153       for (nExp=1;nExp<=globs->nrExp;nExp++) 
00154         {
00155           for (i=1; i<ex[nExp].nPoints; i++)
00156             {
00157               for (j=1; j<=ex[nExp].nvar; j++)
00158                 {
00159                   S0[nExp][i][j] = ex[nExp].yTry[i][j];
00160                   dSSave[nExp][i][j]=ex[nExp].dS[i][j];
00161                   dS[nExp][i][j]=ex[nExp].dS[i][j];
00162                 }
00163             }
00164           for (i=1; i<=globs->npar; i++)
00165             {
00166               P0[nExp][i] = ex[nExp].par[i];
00167               dPSave[nExp][i]=ex[nExp].dP[i];
00168               dP[nExp][i]=ex[nExp].dP[i];
00169             }
00170         }
00171       
00172       uS=1.;
00173       //damping
00174       if(!globs->nodamp)
00175         lambda=dampIt(globs,ex,P0,S0,dP,dS,&uS);
00176       else
00177         lambda=1.;
00178 
00179       globs->Lambda=lambda;
00180 
00181       globs->chisq=sum;
00182       *dbg << "Chisq = " << globs->chisq << endl;
00183       
00184       // next guess
00185       *dbg << "using lambda=" << lambda << "\n";
00186       
00187       sum=0;
00188       for (nExp=1;nExp<=globs->nrExp;++nExp) 
00189         {
00190           for (i=1; i<ex[nExp].nPoints; ++i)
00191             {
00192               for (j=1; j<=ex[nExp].nvar; ++j) 
00193                 {
00194                   sum += (lambda*uS*dSSave[nExp][i][j])*(lambda*uS*dSSave[nExp][i][j]);
00195                 }
00196             }
00197           for (i=1; i<=globs->npar; ++i) 
00198             {
00199               if(globs->doP[i]=='L')
00200                 sum  += (lambda*dPSave[nExp][i])*(lambda*dPSave[nExp][i]);
00201             }
00202         }
00203       for (i=1; i<=globs->npar; ++i) 
00204         {
00205           if(globs->doP[i]!='L')
00206             sum  += (lambda*dPSave[1][i])*(lambda*dPSave[1][i]);
00207         }
00208       
00209       tryOnceMore = (sqrt(sum)/lambda > minImprovement);
00210       
00211       if(uS==0)
00212         tryOnceMore=TRUE;
00213 
00214       if (tryOnceMore) 
00215         {
00216           // set up new initial values for next iteration
00217           for (nExp=1;nExp<=globs->nrExp;++nExp) 
00218             {
00219               for (i=1; i<ex[nExp].nPoints; ++i)
00220                 {
00221                   for (j=1; j<=ex[nExp].nvar; ++j) 
00222                     {
00223                       
00224                       ex[nExp].yTry[i][j] = S0[nExp][i][j]+lambda*uS*dSSave[nExp][i][j];
00225                       // positiveness of initial cond --> test
00226                       if(ex[nExp].yTry[i][j]<0.0)
00227                         ex[nExp].yTry[i][j]=0;
00228                     }
00229                 }
00230               for (i=1; i<=globs->npar; ++i) 
00231                 {
00232                   ex[nExp].par[i]=P0[nExp][i]+lambda*dPSave[nExp][i];
00233                 }
00234             }
00235         } 
00236 
00237       // else:
00238       // keep old values; parCovar is covariance matrix at solution point
00239       // (update is too small to be worth iterating once more)
00240       
00241         ++globs->nIter;
00242         if(globs->silent!=TRUE)
00243           {
00244             cout << "Norm of update vector: " << sqrt(sum) <<"\n\n";
00245             cout.flush();
00246           }
00247         *dbg << "Norm of update vector: " << sqrt(sum) <<"\n\n";
00248         
00249         dbg->flush();
00250         if(globs->wait)
00251           {
00252             cout << "<Press any key to continue>\n";
00253             getchar();
00254           }
00255     }
00256   
00257   if (globs->nIter>=globs->maxit) 
00258     {
00259       if(globs->silent!=TRUE)
00260         {
00261           cout << "fitit: no convergence after " << globs->nIter << " iterations\n";
00262           cout << "PRELIMINARY RESULTS! \n";
00263         }
00264       globs->fitConverged=FALSE;
00265       
00266     }
00267   else
00268     globs->fitConverged=TRUE;
00269   //to obtain the covariance matrix
00270   *dbg << "last call of solvLin to obtain the covariance matrix\n";
00271   if(globs->reg=TRUE)
00272     {
00273       globs->minimiser=1;
00274       //globs->reg=FALSE;
00275       solvLin(globs,ex,TRUE);
00276     }
00277   else
00278     {
00279       globs->minimiser=1;
00280       //globs->reg=FALSE;
00281       solvLin(globs,ex,TRUE);
00282     }
00283 
00284   if(!globs->wait && !globs->nowait)
00285     {
00286       cout << "<Press any key to continue>\n";
00287       getchar();
00288     }
00289   
00290   
00291   if(globs->noGnu==FALSE)
00292     {
00293       for(j=1;j<=globs->ngnu;j++)
00294         pclose(globs->gnuFp[j]);
00295     }
00296 
00297   free_dmatrix(P0,1,globs->nrExp,1,globs->npar);
00298   free_d3tensor(S0,1,globs->nrExp,1,maxPoints,1,maxVar);
00299   free_dmatrix(dPSave,1,globs->nrExp,1,globs->npar);
00300   free_d3tensor(dSSave,1,globs->nrExp,1,maxPoints,1,maxVar);       
00301   free_dmatrix(dP,1,globs->nrExp,1,globs->npar);
00302   free_d3tensor(dS,1,globs->nrExp,1,maxPoints,1,maxVar);
00303 }

void freeMem GlobExp ex,
Glob globs,
int  simit
 

Definition at line 11 of file freeMem.cc.

00012 {
00013   long nExp,i;
00014   long nrExp=globs->nrExp;
00015 
00016 
00017   delete globs->gnuFp;
00018   
00019   if(simit==FALSE)
00020     {
00021       for (nExp=1;nExp<=nrExp;++nExp) 
00022         { 
00023           long nPoints=ex[nExp].nPoints;
00024           long nP=globs->npar;
00025           long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs;
00026           long nMeas=ex[nExp].nMeasure;
00027           long  me=ex[nExp].me, mg=ex[nExp].mg;
00028           free_dvector(ex[nExp].errP,1,nP);
00029           free_dvector(ex[nExp].errY0,1,nvar);
00030           free_dmatrix(ex[nExp].yTry,1,nPoints,1,nvar);
00031           free_dmatrix(ex[nExp].yTrySave,1,nPoints,1,nvar);
00032           free_dmatrix(ex[nExp].yComp,1,nPoints,1,nvar);
00033           free_dmatrix(ex[nExp].yPred,1,nMeas,1,nobs);
00034           free_dmatrix(ex[nExp].h,1,nPoints,1,nvar);   
00035           free_dmatrix(ex[nExp].residues,1,nMeas,1,nobs);
00036           free_d3tensor(ex[nExp].dyds,1,nPoints,1,nvar,1,nvar);
00037           free_d3tensor(ex[nExp].dydp,1,nPoints,1,nvar,1,nP);
00038           free_d3tensor(ex[nExp].dmds,1,nMeas,1,nobs,1,nvar);
00039           free_d3tensor(ex[nExp].dmdp,1,nMeas,1,nobs,1,nP);
00040           if (me>0) 
00041             {
00042               free_dvector(ex[nExp].r2,1,me);
00043               free_d3tensor(ex[nExp].dR2ds,1,me,1,nPoints,1,nvar);
00044               free_dmatrix(ex[nExp].dR2dp,1,me,1,nP);
00045             }
00046           if (mg>0) 
00047             {
00048               free_dvector(ex[nExp].r3,1,mg);
00049               free_d3tensor(ex[nExp].dR3ds,1,mg,1,nPoints,1,nvar);
00050               free_dmatrix(ex[nExp].dR3dp,1,mg,1,nP);
00051             }
00052           
00053           free_dvector(ex[nExp].ua,1,nMeas*nobs);
00054           free_dmatrix(ex[nExp].Ea,1,nMeas*nobs,1,nvar);
00055           free_dmatrix(ex[nExp].Pa,1,nMeas*nobs,1,nP);
00056           free_dvector(ex[nExp].ue,1,me);
00057           free_dmatrix(ex[nExp].Ee,1,me,1,nvar);
00058           free_dmatrix(ex[nExp].Pe,1,me,1,nP);
00059           free_dvector(ex[nExp].ug,1,mg);
00060           free_dmatrix(ex[nExp].Eg,1,mg,1,nvar);
00061           free_dmatrix(ex[nExp].Pg,1,mg,1,nP);
00062           free_dmatrix(ex[nExp].dS,1,nPoints,1,nvar);
00063           free_dvector(ex[nExp].dP,1,nP);
00064           
00065         }
00066     }
00067   for (nExp=1;nExp<=nrExp;++nExp) 
00068     {
00069       free(ex[nExp].splineNodes);
00070       free(ex[nExp].splineY);
00071       free(ex[nExp].splineGam);
00072       free(ex[nExp].nNodes);
00073       for(i=1;i<=NSPLINES;i++)
00074         {
00075           free(ex[nExp].splineNodes[i]);
00076           free(ex[nExp].splineY[i]);
00077           free(ex[nExp].splineGam[i]);
00078         }
00079     }
00080   if(simit==FALSE && globs->covar!=NULL)
00081     free_dmatrix(globs->covar,1,globs->fitdim,1,globs->fitdim);
00082 }

void globOpt GlobExp  ex[],
Glob globs
 

Definition at line 89 of file globOpt.cc.

References exPtr, globsPtr, lvector(), Glob::nIter, GlobExp::nPoints, and Glob::nrExp.

00090 {
00091   globsPtr=globs;
00092   exPtr=ex;
00093 
00094 
00095   globsPtr->nIter=0;
00096   long i,j;
00097   //variables for libSRES
00098   ESParameter *param;
00099   ESPopulation *population;
00100   ESStatistics *stats;
00101   ESfcnTrsfm *trsfm;
00102   unsigned int seed;
00103   int es;
00104   int dim=0;
00105   double *ub, *lb;
00106   int miu, lambda, gen;
00107   double gamma, alpha, varphi;
00108   int retry;
00109   double pf; 
00110 
00111   //keeping old values
00112   long unsigned *nms=lvector(1,globsPtr->nrExp);
00113   for(i=1;i<=globsPtr->nrExp;i++)
00114     {
00115       nms[i]=exPtr[i].nPoints;
00116       exPtr[i].nPoints=2;
00117     }
00118 
00119   seed = shareDefSeed;
00120   gamma = esDefGamma;
00121   alpha = esDefAlpha;
00122   varphi = esDefVarphi;
00123   retry = esDefRetry;
00124   pf = essrDefPf;
00125   es = esDefESSlash;
00126   miu = 30;
00127   lambda = 200;
00128   gen = 1750;
00129   trsfm = NULL;
00130 
00131   //determine # of decision variables
00132   for(i=1;i<=globsPtr->nrExp;i++)
00133     {
00134       for(j=1;j<=NEQNS;j++)
00135         {
00136           if(globsPtr->y0fix[j]==TRUE)
00137             dim++;
00138         }
00139       for(j=1;j<=globsPtr->npar;j++)
00140         {
00141           if(globsPtr->doP[j]=='L')
00142             dim++;
00143         }
00144     }      
00145   for(j=1;j<=globsPtr->npar;j++)
00146     {
00147       if(globsPtr->doP[j]==TRUE)
00148         dim++;
00149     }
00150 
00151 
00152   //bounds ...for test only
00153 
00154   lb=dvector(0,dim);
00155   ub=dvector(0,dim);
00156 
00157   for(i=0;i<dim;i++)
00158     {
00159       lb[i]=0;
00160       ub[i]=2;
00161     }
00162 
00163   ESInitial(seed, &param, trsfm, fitness, es,0, dim, ub, lb, \
00164             miu, lambda, gen, gamma, alpha, varphi,  \
00165             retry, &population, &stats);
00166   
00167   for(int i=1;i<=globsPtr->maxit;i++)
00168     { 
00169       globsPtr->nIter++;
00170       ESStep(population, param, stats, pf);
00171       if(stats->curgen >= param->gen)
00172         break;
00173     }
00174   
00175   //copying old stuff back
00176   for(i=1;i<=globsPtr->nrExp;i++)
00177     {
00178       exPtr[i].nPoints=nms[i];
00179     }
00180 
00181   cout << endl;
00182 
00183   free_dvector(lb,0,dim);
00184   free_dvector(ub,0,dim);
00185   free_lvector(nms,1,globsPtr->nrExp);
00186 }

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

Definition at line 98 of file initialise.cc.

Referenced by mexFunction().

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 }

main int  argc,
char *  argv[]
 

function main

Definition at line 81 of file diffit.cc.

References dbg, Glob::nrExp, and parseopts().

00082 {
00083   long k,i,j;
00084   char outstr[100];
00085   Glob globs;
00086   GlobExp *ex;
00087   
00088   //  parseopt
00089   ex=parseopts(argc,argv,&globs,outstr);
00090 
00091   //open DEBUG stream
00092   dbg=new ofstream("diffit.dbg");
00093   if (!dbg) 
00094     {
00095       cerr << "Error opening DEBUG file.\n";
00096       exit(1);
00097     }
00098   dbg->precision(4);     // set output format for floating point data
00099   dbg->unsetf(ios::floatfield);
00100 
00101   //print debugging information
00102   *dbg << DefModelDescription << "\n";
00103   if (globs.nrExp==1) 
00104     *dbg << "1 experiment\n";
00105   else
00106     *dbg << globs.nrExp << " experiments\n";
00107    *dbg << "\n";
00108    for(i=1;i<=globs.nrExp;i++)
00109      {
00110        *dbg << "Experiment: " << i << "\n";
00111        *dbg << NPARAMS << " parameter(s):";
00112        for (k=1; k<=NPARAMS; ++k) 
00113          *dbg << " " << ex[i].par[k];
00114        *dbg << "\n";
00115      }
00116    
00117   dbg->flush();
00118  
00119   for(i=1;i<=globs.nrExp;i++)
00120     {
00121       // readData()
00122       readData(ex,&globs,i);
00123       //set mesh
00124       setMesh(ex,&globs,i);
00125 
00126 #ifdef PRINTDATA
00127       *dbg << "\nData: \n";
00128       for (j=1; j<=ex[i].nMeasure; j++) {
00129         *dbg << ex[i].xMeasure[j];
00130         for (k=1; k<=ex[i].nobs; k++)
00131           *dbg << "\t" << ex[i].yMeasure[j][k] << "\t" << ex[i].sigma[j][k];
00132         *dbg << "\n";
00133       }
00134       *dbg << "\n";
00135 #endif
00136       *dbg << "Mesh:";
00137       for (k=1;k<=ex[i].nPoints;++k) 
00138         *dbg << " " << ex[i].mesh[k];
00139       *dbg << "\n";
00140     }
00141 
00142   //initial values
00143   for(i=1;i<=globs.nrExp;i++)
00144     {
00145       *dbg << "\nExperiment " << i << ":\n";
00146       if (ex[i].y0) 
00147         {
00148           *dbg << NEQNS << " starting value(s):";
00149           for (k=1; k<=NEQNS; ++k) 
00150             *dbg << ' ' << ex[i].y0[k];
00151           *dbg << "\n";
00152         } 
00153       else 
00154         {
00155           *dbg << "No starting values specified\n";
00156         }
00157       *dbg << "\n";
00158       dbg->flush();
00159     }
00160   
00161    initialise(ex,&globs,FALSE);
00162    
00163    //simulate initial state
00164    if(globs.simInit==TRUE)
00165      {
00166        for (i=1;i<=globs.nrExp;++i) 
00167          simInit(&ex[i],&globs);
00168      }
00169    
00170    // starting of the numerics
00171    //*************************
00172       
00173    if(globs.strategy==2)
00174      globOpt(ex,&globs);
00175    else
00176      {
00177        try 
00178          {            
00179             fitIt(ex,&globs);
00180          }
00181        catch(int i)  //Exeption handling
00182          {
00183            if(globs.noGnu==FALSE)
00184              {
00185                for(j=1;j<=globs.ngnu;j++)
00186                  pclose(globs.gnuFp[j]);
00187              }       
00188            exit(1);
00189          }
00190      }
00191 
00192    
00193    //*************************
00194 
00195   //Output after convergence
00196   outFit(ex,&globs);
00197   if(!globs.noGnu)
00198     system("rm -f gnuout.dat");
00199 
00200   freeMem(ex,&globs,FALSE);
00201 }

void outFit GlobExp  ex[],
Glob globs
 

Definition at line 11 of file outFit.cc.

00012 {
00013   long i,j;
00014   long nExp;
00015   long nvar,npar,nglob=0;
00016   long ind=1;
00017 
00018   npar=globs->npar;     
00019   nvar=ex[1].nvar;
00020 
00021   double **errorS=dmatrix(1,globs->nrExp,1,nvar);
00022   double **errorP=dmatrix(1,globs->nrExp,1,npar);
00023   
00024 
00025   if(globs->strategy!=2)
00026     {
00027       //prepare output
00028       for(nExp=1;nExp<=globs->nrExp;++nExp) 
00029         {
00030           for(i=1;i<=nvar;i++)
00031             errorS[nExp][i]=-1;
00032           for(i=1;i<=npar;i++)
00033             errorP[nExp][i]=-1;
00034           
00035           for(i=1;i<=nvar;i++)
00036             {
00037               if(globs->y0fix[i]!=FALSE)
00038                 {
00039                   errorS[nExp][i]=sqrt(globs->covar[ind][ind]);
00040                   ind++;
00041                 }
00042             }
00043           for(i=1;i<=npar;i++)
00044             {
00045               if(globs->doP[i]=='L')
00046                 {
00047                   errorP[nExp][i]=sqrt(globs->covar[ind][ind]);
00048                   ind++;
00049                   nglob++;
00050                 }
00051             }
00052         }
00053       
00054       for(i=1;i<=npar;i++)
00055         {
00056           if(globs->doP[i]==TRUE)
00057             {
00058               for(nExp=1;nExp<=globs->nrExp;++nExp) 
00059                 {
00060                   errorP[nExp][i]=sqrt(globs->covar[ind][ind]);
00061                 }
00062               ind++;
00063             }
00064         } 
00065     }  //end if(globs->strategy..)
00066 
00067   /* print results */
00068   if(globs->silent!=TRUE)
00069     {
00070       cout << "Number of iterations: " << globs->nIter << ", chi^2 = " << globs->chisq << "\n\n";
00071       cout << "\nBest fit parameters +/- standard errors:\n";
00072       cout << "----------------------------------------\n\n"; 
00073       
00074       cout << "Global Parameters:\n";
00075       for(j=1;j<=npar;j++)
00076         { 
00077           if(globs->doP[j]!='L')
00078             {
00079               if(globs->strategy!=2)
00080                 {
00081                   for(nExp=1;nExp<=globs->nrExp;++nExp) 
00082                     ex[nExp].errP[j]=errorP[1][j];
00083                   cout << ParameterNames[j-1];
00084                   if(errorP[1][j]!=-1)
00085                     cout << " = " << ex[1].par[j] <<  " +/- " << errorP[1][j] << endl;
00086                   else
00087                     cout << " (fixed)\n";
00088                 }
00089               else
00090                 cout << ParameterNames[j-1] << " = " << ex[1].par[j] << endl;
00091                 
00092                 
00093             }
00094         }
00095       cout << endl;
00096       for(i=1;i<=globs->nrExp;++i) 
00097         {
00098           cout << "Experiment " << i << ":";
00099           if(nglob!=0)
00100             cout << "\n\nLocal Parameters:\n";
00101           else
00102             cout << endl;
00103           
00104           for(j=1;j<=npar;j++)
00105             {
00106               if(globs->doP[j]=='L')
00107                 {
00108                   if(globs->strategy!=2)
00109                     {
00110                       ex[i].errP[j]=errorP[i][j];
00111                       cout << ParameterNames[j-1] << " = " << ex[i].par[j];
00112                       if(errorP[i][j]!=-1)
00113                         cout <<  " +/- " << errorP[i][j] << endl;
00114                       else
00115                         cout << " (fixed)\n";
00116                     }
00117                   else
00118                     cout << ParameterNames[j-1] << " = " << ex[i].par[j] << endl;
00119                 }
00120             }
00121           cout << "Initial Values" << endl;
00122           for(j=1;j<=nvar;j++)
00123             {
00124               if(globs->strategy!=2)
00125                 {
00126                   ex[i].errY0[j]=errorS[i][j];
00127                   cout << VariableNames[j-1] << " = " << ex[i].yTry[1][j];
00128                   if(errorS[i][j]!=-1)
00129                     cout <<  " +/- " << errorS[i][j] << endl;
00130                   else
00131                     cout << " (fixed)\n";
00132                 }
00133               else
00134                 cout << VariableNames[j-1] << " = " << ex[i].yTry[1][j] << endl;
00135             }
00136           cout << endl;
00137         }
00138     } //end of --> if(silent!=TRUE)
00139 
00140   free_dmatrix(errorS,1,globs->nrExp,1,nvar);
00141   free_dmatrix(errorP,1,globs->nrExp,1,npar);
00142 }

GlobExp * parseopts int  argc,
char *  argv[],
Glob globs,
char *  outstr
 

Main module routine doing all the pasing.

Parameters:
argc The number of command line arguments as passed by main.
*argv[] String of the arguments.
*globs The structure Glob.
*outstr Depreciated, not used.

Definition at line 120 of file parse.cc.

Referenced by main().

00121 {
00122   int longindex,opt;
00123   int parlistSpecified=FALSE;
00124   long k,l,nExp,sp;
00125   char name[MAXLENGTH],line[MAXLENGTH];
00126 
00127   // option definition
00128   static struct option longopts[]={
00129     {"help", 0, 0,  'h'},
00130     {"help", 0, 0,  '?'},
00131     {"x0"  , 1, 0,   1 },
00132     {"x1"  , 1, 0,   2 },
00133     {"nms"  ,1, 0,   3 },
00134     {"f"    ,1, 0,   4 },
00135     {"p"    ,1, 0,   5 },
00136     {"y0"   ,1, 0,   6 },
00137     {"nobs" ,1, 0,   7 },
00138     {"eps"  ,1, 0,   8 },
00139     {"nognu",0, 0,   9 },
00140     {"nomeasure",0,0,10},
00141     {"doP"  ,1, 0,   11},
00142     {"maxit",1, 0,   12},
00143     {"wait" ,0, 0,   13},
00144     {"usesig",0,0,   14},
00145     {"int"   ,1,0,   15},
00146     {"maxstp",1,0,   17},
00147     {"minimp",1,0,   18},
00148     {"nowait",0,0,   19},
00149     {"elastic" ,1,0, 21},
00150     {"reg"     ,0,0, 22},
00151     {"epsilon" ,1,0, 23},
00152     {"lambda"  ,1,0, 24},
00153     {"spline"  ,1,0, 25},
00154     {"siminit" ,0,0, 26},
00155     {"pert"    ,1,0, 27},
00156     {"y0fix"   ,1,0, 28},
00157     {"nodamp"  ,0,0, 29},
00158     {"strat"   ,1,0, 30},
00159     {"opt"     ,1,0, 31},
00160     {"Lconst"  ,1,0, 32},
00161     {0, 0, 0, 0}
00162   };
00163   //initialise some global parameters in globs
00164   globs->noGnu=FALSE;
00165   globs->eps=1e-6;
00166   globs->npar=NPARAMS;
00167   globs->noMeasurements=FALSE;
00168   globs->doP=ivector(1,globs->npar);
00169   for(k=1;k<=globs->npar;k++)
00170     globs->doP[k]=TRUE;
00171   globs->y0fix=ivector(1,NEQNS);
00172   for(k=1;k<=NEQNS;k++)
00173     globs->y0fix[k]=TRUE;
00174   globs->maxit=1000;
00175   globs->gnuFp=NULL;
00176   globs->wait=FALSE;
00177   globs->usesig=FALSE;
00178   globs->integrator=1;
00179   globs->stiff=TRUE;
00180   globs->maxstp=5000;
00181   globs->minimp=1e-4;
00182   globs->nowait=FALSE;
00183   globs->elastic=1.;
00184   globs->reg=FALSE;
00185   globs->epsilon=1e-10;
00186   globs->lambda=1e6;
00187   globs->simInit=FALSE;
00188   globs->pert=0.;
00189   globs->nodamp=FALSE;
00190   globs->initSpline=TRUE;
00191   globs->wquer=-1;
00192   globs->silent=FALSE;
00193   globs->strategy=1;
00194   globs->minimiser=2;
00195   globs->faktorL=dvector(1,NPARAMS);
00196   globs->faktorLexist=FALSE;
00197   for(k=1;k<=NPARAMS;k++)
00198   {
00199         globs->faktorL[k]=-1;
00200   }     
00201   while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00202     {
00203       switch(opt)
00204         {
00205         case 4:
00206           parlistSpecified=TRUE;
00207           strncpy(name,optarg,MAXLENGTH);
00208           break;
00209         }
00210     }
00211   optind=0;
00212   
00213   GlobExp *ex;
00214   if(parlistSpecified==FALSE)
00215     {
00216       globs->nrExp=1;
00217       ex=new GlobExp[globs->nrExp+1];   
00218       //initialise some global parameters
00219       ex[1].fitstart=-1e99;
00220       ex[1].fitend=1e99;
00221       ex[1].nobs=NOBS;
00222       ex[1].nvar=NEQNS;
00223       ex[1].nPoints=2;
00224       ex[1].splinesDefined=FALSE;
00225       ex[1].y0=NULL;
00226       // initialisation of the parameters
00227       ex[1].par=dvector(1,NPARAMS);
00228       for(k=1;k<=globs->npar;k++)
00229         ex[1].par[k]=DefParameters[k-1];
00230       
00231       if(NSPLINES!=0)
00232         {
00233           ex[1].nNodes=lvector(1,NSPLINES);
00234         }
00235     }
00236   else //if parlist specified
00237     {
00238       ifstream in;      ofstream out;
00239       char str[MAXLENGTH];
00240 
00241       in.open(name);
00242 
00243       if(!in)
00244         {
00245           cerr << "parameter list " << name << "  not found\n";
00246           exit(1);
00247         }
00248       
00249       //preprocessing of data list
00250       globs->nrExp=0;
00251       sprintf(str,"%s.%d",name,globs->nrExp);
00252       out.open(str);
00253       while(!in.eof())
00254         {
00255           in.getline(str,MAXLENGTH,'\n');
00256           if(str[0]=='@')
00257             {
00258               out.flush();
00259               out.close();
00260               globs->nrExp++;
00261               sprintf(str,"%s.%d",name,globs->nrExp);
00262               out.open(str);
00263             }
00264           else if(str[0]!='#')
00265             {
00266               out << str << endl;
00267             }
00268         }
00269       out.flush();
00270       out.close();
00271       in.close();
00272 
00273       if(globs->nrExp==0)
00274         {
00275           cerr << "insuffient number of experiments, @ missing\n";
00276           exit(1);
00277         }
00278       
00279       long _argc;
00280       char *_argv[MAXLENGTH];
00281 
00282       for(k=0;k<=MAXLENGTH;k++)
00283         _argv[k]=(char *)malloc((size_t) (150*sizeof(char)));
00284 
00285       {
00286         ifstream inp;
00287         
00288         _argc=0;
00289         sprintf(str,"%s.0",name);
00290         inp.open(str);
00291         
00292         
00293         while(!inp.eof())
00294           {
00295             inp >> _argv[_argc+1];
00296             _argc++;
00297           }
00298         inp.close();
00299       }
00300    
00301       //PARSING LIST
00302 
00303       while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00304         {
00305           switch(opt) 
00306             { 
00307             case 'h':
00308               cout << usage << endl;
00309               exit(1);
00310               break;
00311             case '?':
00312               cerr << usage << endl;
00313               exit(1);  
00314               break;
00315             case 8:
00316               globs->eps=fabs(atof(optarg));
00317               break;
00318             case 9:
00319               globs->noGnu=TRUE;
00320               break;
00321             case 10:
00322               globs->noMeasurements=TRUE;
00323               break;
00324             case 11:
00325               if(strlen(optarg)!=globs->npar)
00326                 {
00327                   cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
00328                   exit(1);
00329                 }
00330               for(k=0;k < globs->npar;k++)
00331                 {
00332                   if(optarg[k]=='0')
00333                     globs->doP[k+1]=FALSE;
00334                   else if (optarg[k]=='L')
00335                     globs->doP[k+1]='L'; 
00336                   else
00337                     globs->doP[k+1]=TRUE;
00338                 }
00339               break;
00340             case 12:
00341               globs->maxit=abs(atol(optarg));
00342               break;
00343             case 13:
00344               globs->wait=TRUE;
00345               break;
00346             case 14:
00347               globs->usesig=TRUE;
00348               break;
00349             case 15:
00350               globs->integrator=abs(atoi(optarg));
00351               break;
00352             case 17:
00353               globs->maxstp=abs(atoi(optarg));
00354               break;
00355             case 18:
00356               globs->minimp=fabs(atof(optarg));
00357               break;
00358             case 19:
00359               globs->nowait=TRUE;
00360               break;
00361             case 21:
00362               globs->elastic=fabs(atof(optarg));
00363               if(globs->elastic > 1. || globs->elastic == 0.)
00364                 {
00365                   cerr << "Insufficient range of -elast <val> \n";
00366                   exit(1);
00367                 }
00368               break;
00369             case 22:
00370               globs->reg=TRUE;
00371               break;
00372             case 23:
00373               globs->epsilon=fabs(atof(optarg));
00374               break;
00375             case 24:
00376               globs->lambda=fabs(atof(optarg));
00377               break;
00378             case 26:
00379               globs->simInit=TRUE;
00380               break;
00381             case 27:
00382               globs->pert=fabs(atof(optarg));
00383               break;
00384             case 28:
00385               if(strlen(optarg)!=NEQNS)
00386                 {
00387                   cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
00388                   exit(1);
00389                 }
00390               for(k=0;k < NEQNS;k++)
00391                 {
00392                   if(optarg[k]=='0')
00393                     globs->y0fix[k+1]=FALSE;
00394                   else
00395                     globs->y0fix[k+1]=TRUE;
00396                 }
00397               break;
00398             case 29:
00399               globs->nodamp=TRUE;
00400               break; 
00401             case 30:
00402               globs->strategy=abs(atoi(optarg));
00403               break;
00404             case 31:
00405               globs->minimiser=abs(atoi(optarg));
00406               break;
00407             case 32:
00408               get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00409               globs->faktorLexist=TRUE;
00410               break;                  
00411             default:
00412               cerr << endl;
00413               cerr << "Parsing parameter list produced errors.\n\n";
00414               exit(1);
00415             }
00416         }
00417 
00418       optind=0;
00419       ex=new GlobExp[globs->nrExp+1];   
00420   
00421       //PARSING experiment specific list
00422       for(nExp=1;nExp<=globs->nrExp;nExp++)
00423         {
00424           //initialise some global parameters
00425           ex[nExp].fitstart=-1e99;
00426           ex[nExp].fitend=1e99;
00427           ex[nExp].nobs=NOBS;
00428           ex[nExp].nvar=NEQNS;
00429           ex[nExp].nPoints=2;
00430           ex[nExp].splinesDefined=FALSE;
00431           ex[nExp].y0=NULL;
00432           // initialisation of the parameters
00433           ex[nExp].par=dvector(1,NPARAMS);
00434           for(k=1;k<=globs->npar;k++)
00435             ex[nExp].par[k]=DefParameters[k-1];
00436           
00437           if(NSPLINES!=0)
00438             {
00439               ex[nExp].nNodes=lvector(1,NSPLINES);
00440             }
00441 
00442           _argc=0;
00443 
00444           ifstream inp;
00445           sprintf(str,"%s.%d",name,nExp);
00446           inp.open(str);
00447 
00448           while(!inp.eof())
00449             {
00450               inp >> _argv[_argc+1];
00451               _argc++;
00452             }
00453 
00454           inp.close();
00455           while((opt=getopt_long_only(_argc,_argv,"h?",longopts,&longindex)) != -1)
00456              {
00457                switch(opt)
00458                  {
00459                  case 1:
00460                    ex[nExp].fitstart=atof(optarg);
00461                    break;
00462                  case 2:
00463                    ex[nExp].fitend=atof(optarg);
00464                    break;
00465                  case 3:
00466                    ex[nExp].nPoints=abs(atol(optarg));
00467                    ex[nExp].nPoints+=1;
00468                    break;
00469                  case 5:
00470                    get_list(globs->npar,ex[nExp].par,optarg," parameter ");
00471                    break;
00472                  case 6:
00473                    ex[nExp].y0=dvector(1,ex[nExp].nvar);
00474                    for(k=1;k<=ex[nExp].nvar;k++)
00475                      ex[nExp].y0[k]=0.;
00476                    k=get_list(ex[nExp].nvar,ex[nExp].y0,optarg," initial values ");
00477                    if(k!=ex[nExp].nvar)
00478                      {
00479                        cerr << ex[nExp].nvar << " initial values required.\n";
00480                        exit(1);
00481                      }
00482                    break;
00483                  case 7:
00484                    ex[nExp].nobs=atol(optarg);
00485                    break;
00486                  case 25:
00487                    k=0;
00488                    l=1;
00489                    ex[nExp].splinesDefined=TRUE;
00490                    while(optarg[k]!='\0')
00491                      {
00492                        if(optarg[k]==',')
00493                          l++;
00494                        k++;
00495                      }
00496                    if(l!=NSPLINES)
00497                      {
00498                        cerr << "Incompatible number of splines.\n";
00499                        exit(1);
00500                      }
00501                    ex[nExp].splineFile=new string[l+1];
00502                    l=1;
00503                    sp=0;
00504                    k=0;
00505                    while(optarg[k]!='\0')
00506                      {
00507                        if(optarg[k]==',')
00508                          {
00509                            line[sp]='\0';
00510                            ex[nExp].splineFile[l]=(string)line;
00511                            sp=0;
00512                            l++;
00513                          }
00514                        else
00515                          {
00516                            line[sp]=optarg[k];
00517                            sp++;
00518                          }
00519                        k++;
00520                      }
00521                    line[sp]='\0';
00522                    ex[nExp].splineFile[l]=(string)line;
00523                    break;
00524                  default:
00525                    cerr << endl;
00526                    cerr << "Parsing parameter list produced errors.\n\n";
00527                    exit(1);
00528                  }
00529              }
00530 
00531            //Parsing file name
00532            if((_argc-optind)!=1)
00533              {
00534                cerr << "No/Too many datafile(s) specified for experiment " << nExp << " .  \n\n";
00535                exit(1);
00536              }
00537            else
00538              {
00539                ex[nExp].fileName=new char[strlen(_argv[optind])+1];
00540                strcpy(ex[nExp].fileName,_argv[optind]);
00541              }
00542            //checking if x0 < x1
00543            if(ex[nExp].fitstart >= ex[nExp].fitend)
00544              {
00545                cerr << "x0 must be smaller than x1.\n";
00546                exit(1);
00547              }
00548            optind=0;
00549         }
00550       
00551      
00552       //for(k=0;k<=MAXLENGTH;k++)
00553       //free(_argv[k]);
00554 
00555       sprintf(str,"rm -f %s.*",name);
00556       system(str);
00557     }
00558 
00559   optind=0;
00560   
00561   while((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
00562     {
00563       switch(opt)
00564         {
00565         case 'h':
00566           cout << usage << endl;
00567           exit(1);
00568           break;
00569         case '?':
00570           cerr << usage << endl;
00571           exit(1);      
00572           break;
00573         case 1:
00574           ex[1].fitstart=atof(optarg);
00575           break;
00576         case 2:
00577           ex[1].fitend=atof(optarg);
00578           break;
00579         case 3:
00580           ex[1].nPoints=abs(atol(optarg));
00581           ex[1].nPoints+=1;
00582           break;
00583         case 4:
00584           break;
00585         case 5:
00586           get_list(globs->npar,ex[1].par,optarg," parameter ");
00587           break;
00588         case 6:
00589           ex[1].y0=dvector(1,ex[1].nvar);
00590           for(k=1;k<=ex[1].nvar;k++)
00591             ex[1].y0[k]=0.;
00592           k=get_list(ex[1].nvar,ex[1].y0,optarg," initial values ");
00593           if(k!=ex[1].nvar)
00594             {
00595               cerr << ex[1].nvar << " initial values required.\n";
00596               exit(1);
00597             }
00598           break;
00599         case 7:
00600           ex[1].nobs=atol(optarg);
00601           break;
00602         case 8:
00603           globs->eps=fabs(atof(optarg));
00604           break;
00605         case 9:
00606           globs->noGnu=TRUE;
00607           break;
00608         case 10:
00609           globs->noMeasurements=TRUE;
00610           break;
00611         case 11:
00612           if(strlen(optarg)!=globs->npar)
00613             {
00614               cerr << strlen(optarg) << " parameter(s) specified instead of " <<  globs->npar << ".\n";
00615               exit(1);
00616             }
00617           for(k=0;k < globs->npar;k++)
00618             {
00619               if(optarg[k]=='0')
00620                 globs->doP[k+1]=FALSE;
00621               else if (optarg[k]=='L')
00622                 globs->doP[k+1]='L'; 
00623               else
00624                 globs->doP[k+1]=TRUE;
00625             }
00626           break;
00627         case 12:
00628           globs->maxit=abs(atol(optarg));
00629           break;
00630         case 13:
00631           globs->wait=TRUE;
00632           break;
00633         case 14:
00634           globs->usesig=TRUE;
00635           break;
00636         case 15:
00637           globs->integrator=abs(atoi(optarg));
00638           break;
00639         case 17:
00640           globs->maxstp=abs(atoi(optarg));
00641           break;
00642         case 18:
00643           globs->minimp=fabs(atof(optarg));
00644           break;
00645         case 19:
00646           globs->nowait=TRUE;
00647           break;
00648         case 21:
00649           globs->elastic=fabs(atof(optarg));
00650           if(globs->elastic > 1. || globs->elastic == 0.)
00651             {
00652               cerr << "Insufficient range of -elast <val> \n";
00653               exit(1);
00654             }
00655           break;
00656         case 22:
00657           globs->reg=TRUE;
00658           break;
00659         case 23:
00660           globs->epsilon=fabs(atof(optarg));
00661           break;
00662         case 24:
00663           globs->lambda=fabs(atof(optarg));
00664           break;
00665         case 25:
00666           k=0;
00667           l=1;
00668           ex[1].splinesDefined=TRUE;
00669           while(optarg[k]!='\0')
00670             {
00671               if(optarg[k]==',')
00672                 l++;
00673               k++;
00674             }
00675           if(l!=NSPLINES)
00676             {
00677               cerr << "Incompatible number of splines.\n";
00678               exit(1);
00679             }
00680           ex[1].splineFile=new string[l+1];
00681           l=1;
00682           sp=0;
00683           k=0;
00684           while(optarg[k]!='\0')
00685             {
00686               if(optarg[k]==',')
00687                 {
00688                   line[sp]='\0';
00689                   ex[1].splineFile[l]=(string)line;
00690                   sp=0;
00691                   l++;
00692                 }
00693               else
00694                 {
00695                   line[sp]=optarg[k];
00696                   sp++;
00697                 }
00698               k++;
00699             }
00700           line[sp]='\0';
00701           ex[1].splineFile[l]=(string)line;
00702           break;           
00703         case 26:
00704           globs->simInit=TRUE;
00705           break;
00706         case 27:
00707           globs->pert=fabs(atof(optarg));
00708           break;
00709         case 28:
00710           if(strlen(optarg)!=NEQNS)
00711             {
00712               cerr << strlen(optarg) << " variable(s) specified instead of " <<  NEQNS << ".\n";
00713               exit(1);
00714             }
00715           for(k=0;k < NEQNS;k++)
00716             {
00717               if(optarg[k]=='0')
00718                 globs->y0fix[k+1]=FALSE;
00719               else
00720                 globs->y0fix[k+1]=TRUE;
00721             }
00722           break;
00723         case 29:
00724           globs->nodamp=TRUE;
00725           break;
00726         case 30:
00727           globs->strategy=abs(atoi(optarg));
00728           break;
00729         case 31:
00730           globs->minimiser=abs(atoi(optarg));
00731           break;
00732         case 32:
00733           get_list(globs->npar,globs->faktorL,optarg," local paramter constraints ");
00734           globs->faktorLexist=TRUE;
00735           break;
00736         default:
00737           cerr << endl;
00738           cerr << "Parsing command line options produced errors \n\n";
00739           exit(1);
00740         }
00741     }
00742   //checking if x0 < x1
00743   if(ex[1].fitstart >= ex[1].fitend)
00744     {
00745       cerr << "x0 must be smaller than x1.\n";
00746       exit(1);
00747     }
00748   
00749   //Parsing file name for a single experiment
00750   if((argc-optind)!=1 && parlistSpecified==FALSE)
00751     {
00752       cerr << "No/Too many datafile(s) specified.  \n\n";
00753       exit(1);
00754     }
00755   else if(parlistSpecified==FALSE)
00756     {
00757       ex[1].fileName=new char[strlen(argv[optind])+1];
00758       strcpy(ex[1].fileName,argv[optind]);
00759     }
00760   
00761   if(NSPLINES!=0)
00762     {
00763       for(k=1;k<=globs->nrExp;k++)
00764         {
00765           if(ex[k].splinesDefined==FALSE)
00766             {
00767               cerr << "Please specify spline data in experiment " << k << ".\n";
00768               exit(1);
00769             }
00770         }
00771     }
00772   
00773   return(ex);
00774 }

void readData GlobExp ex,
Glob globs,
long  expNr
 

Definition at line 17 of file readData.cc.

References GlobExp::fileName, and GlobExp::nobs.

00018 {
00019   int lineNr=0;
00020   long k,l,index;
00021   long nLines=0,maxLeng=0;
00022   long nobsmax=ex[expNr].nobs;
00023   char c,*line,*arg;
00024   double *x;
00025   double **y;
00026   double **sig;
00027 
00028   ifstream in;
00029 
00030   FILE *fp;
00031 
00032   fp=fopen(ex[expNr].fileName,"r");
00033   if (fp==NULL) {
00034     cerr << "Error opening data file " << ex[expNr].fileName << "\n";
00035     exit(1);
00036   }
00037   
00038   //read total number of lines
00039   while((c=fgetc(fp))!=EOF)
00040         {
00041           for(k=1;c!='\n';k++)
00042             {
00043               if(k>maxLeng)
00044                 maxLeng=k;
00045               c=fgetc(fp);
00046             }
00047           nLines++;
00048         }
00049 #ifdef DEBUGREADDATA
00050   *dbg << "Total number of lines in " <<  ex[expNr].fileName << " : " << nLines << endl;
00051   *dbg << "Maximal length of line in " << ex[expNr].fileName << " : " << maxLeng << endl;
00052 #endif
00053   fclose(fp);
00054 
00055   in.open(ex[expNr].fileName);
00056 
00057   //allocate memory
00058   line=new char[maxLeng+2];
00059   arg=new char[maxLeng+2];
00060   x=dvector(1,nLines);
00061   y=dmatrix(1,nLines,1,ex[expNr].nobs);
00062   sig=dmatrix(1,nLines,1,ex[expNr].nobs);
00063 
00064   //extracting the data
00065   while(!in.eof())
00066         {
00067           in.getline(line,maxLeng+2,'\n');
00068           if(isspace(line[0])==FALSE && line[0]!= '#' && line[0] != '\0' )
00069             {
00070               lineNr++;
00071 #ifdef DEBUGREADDATA
00072               *dbg << "Extracting data line no. : " << lineNr << endl;
00073 #endif
00074               l=0;
00075               index=0;
00076               for(k=0;k<=strlen(line);k++)
00077                 {
00078                   if(isspace(line[k])!=FALSE)
00079                     {
00080                       arg[l]='\0';
00081                       l=0;
00082                       if(index==0)
00083                         x[lineNr]=atof(arg);
00084                       else if(index <=2*ex[expNr].nobs && (index % 2)==1)
00085                         y[lineNr][(index+1)/2]=atof(arg);
00086                       else if(index <=2*ex[expNr].nobs && (index % 2)==0)
00087                         sig[lineNr][index/2]=atof(arg);
00088                       index++;
00089                     }
00090                   else
00091                     {
00092                       arg[l]=line[k];
00093                       l++;
00094                     }
00095                 }
00096               arg[l]='\0';      
00097               if(index==0)
00098                 x[lineNr]=atof(arg);
00099               else if(index <=2*ex[expNr].nobs && (index % 2)==1)
00100                 y[lineNr][(index+1)/2]=atof(arg);
00101               else if(index <=2*ex[expNr].nobs && (index % 2)==0)
00102                 sig[lineNr][index/2]=atof(arg);
00103               
00104               //determ. max. complete obs.
00105               if((index % 2)==1)
00106                 index--;
00107               if(index==0)
00108                 {
00109                   cerr << "Data line " << lineNr << " contains no data." << endl;
00110                   exit(1);
00111                 }
00112               if(index/2 <=  nobsmax)
00113                 nobsmax=index/2;
00114             }
00115                     
00116         }
00117   in.close();
00118 
00119   ex[expNr].nobs=nobsmax;
00120   ex[expNr].nMeasure=lineNr;
00121   if(ex[expNr].nMeasure==0)
00122     {
00123       cerr << "File " << ex[expNr].fileName << " contains no data." << endl;
00124       exit(1);
00125     }
00126 
00127   //check if the time is in ascending order
00128     for(k=2;k<=ex[expNr].nMeasure;k++)
00129       {
00130         if(x[k-1] >= x[k])
00131           {
00132             cerr << "Data not in temporal ascending order at dataline " << k << "\n";
00133             exit(1);
00134           }
00135       }
00136     //fit range stuff
00137     if(ex[expNr].fitstart == -1e99)
00138       ex[expNr].fitstart=x[1];
00139     if(ex[expNr].fitend ==  1e99)
00140       ex[expNr].fitend=x[ex[expNr].nMeasure];
00141     
00142     //for memory allocation
00143     index=1;
00144     for(k=1;k<=ex[expNr].nMeasure;k++)
00145       {
00146         if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
00147           { 
00148             index++;
00149           }
00150       }
00151     index--; 
00152   
00153     //the actual amount of data
00154     ex[expNr].nMeasure=index;
00155 
00156     //memory allocation
00157     ex[expNr].xMeasure=dvector(1,ex[expNr].nMeasure);
00158     ex[expNr].yMeasure=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
00159     ex[expNr].sigma=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
00160 
00161     // copy data to data structure
00162     index=1;
00163     for(k=1;k<=lineNr;k++)
00164       {
00165         if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
00166           {
00167             ex[expNr].xMeasure[index]=x[k];
00168             for(l=1;l<=ex[expNr].nobs;l++)
00169               {
00170                 ex[expNr].yMeasure[index][l]=y[k][l];
00171                 ex[expNr].sigma[index][l]=sig[k][l];
00172               }
00173             index++;
00174           }
00175       }
00176     index--;
00177 
00178     //setting some variables
00179     ex[expNr].firstMeasure=1;
00180     ex[expNr].lastMeasure=ex[expNr].nMeasure;
00181     
00182     //free memory
00183     free_dvector(x,1,nLines);
00184     free_dmatrix(y,1,nLines,1,ex[expNr].nobs);
00185     free_dmatrix(sig,1,nLines,1,ex[expNr].nobs);
00186 }

void setMesh GlobExp ex,
Glob globs,
long  expNr
 

Definition at line 11 of file setMesh.cc.

00012 {
00013   long i,nPoints=ex[expNr].nPoints;
00014   
00015   ex[expNr].mesh=dvector(1,nPoints);
00016   ex[expNr].mesh[1]= ex[expNr].fitstart;
00017   ex[expNr].mesh[nPoints]= ex[expNr].fitend;
00018 
00019   if (globs->noMeasurements==FALSE) 
00020     {         
00021       // mesh at data points
00022       if (ex[expNr].nMeasure < nPoints) 
00023         {
00024           cerr << "Too few measurements to construct mesh\n";
00025           exit(1);
00026         }
00027       for (i=2; i < nPoints; ++i)
00028         ex[expNr].mesh[i]=ex[expNr].xMeasure[ex[expNr].firstMeasure+int(double((i-1)*(ex[expNr].nMeasure-1))/double(nPoints-1))];
00029     } 
00030   else 
00031     {        
00032       // evenly spaced mesh
00033       for (i=2; i < nPoints; ++i)
00034         ex[expNr].mesh[i]=ex[expNr].fitstart+(ex[expNr].fitend-ex[expNr].fitstart)*double(i-1)/double(nPoints-1);
00035     }
00036 }

void simInit GlobExp ex,
Glob globs
 

Definition at line 21 of file simInit.cc.

00022 {
00023   long i,j;
00024   long nvar=ex->nvar;
00025   long nPoints=ex->nPoints;
00026   long idum=-42;
00027   double *t=dvector(1,2);
00028   double **y=dmatrix(1,2,1,nvar);
00029   
00030    if (ex->y0) 
00031     {
00032       for (i=1; i<=nvar; ++i) 
00033         ex->yTry[1][i]=ex->y0[i];
00034     } 
00035    else
00036      {
00037         for (i=1; i<=nvar; ++i) 
00038           {
00039             ex->yTry[1][i]=DefYValues[i-1];
00040           }
00041      }
00042 
00043    //integrating 
00044 
00045    for(i=2;i<=nPoints;i++)
00046      {
00047        for(j=1;j<=nvar;j++)
00048          ex->yTry[i][j]=ex->yTry[i-1][j];
00049        t[1]=ex->mesh[i-1];
00050        t[2]=ex->mesh[i];
00051        tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
00052       for(j=1;j<=nvar;j++)
00053         ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
00054      }
00055    
00056 
00057 
00058    free_dvector(t,1,2);
00059    free_dmatrix(y,1,2,1,nvar);
00060 }


Variable Documentation

ofstream* dbg
 

debug stream (C++)

Definition at line 78 of file diffit.cc.


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