#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++) |
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:
Definition in file diffit.cc.
|
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 }
|
|
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 }
|
|
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, ¶m, 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 }
|
|
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 }
|
|
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 }
|
|
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 }
|
|
Main module routine doing all the pasing.
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 }
|
|
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 }
|
|
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 }
|
|
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 }
|
|
debug stream (C++)
|