- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /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
GlobExp * | parseopts (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.