#include "mex.h"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include "def.h"
#include "model.h"
#include "nr.h"
Go to the source code of this file.
Defines | |
#define | IN prhs[0] |
#define | N plhs[0] |
#define | YOUT plhs[1] |
#define | DMDS plhs[2] |
#define | DMDP plhs[3] |
Functions | |
GlobExp * | parseopts (int argc, char *argv[], Glob *globs, char *outstr) |
void | initialise (GlobExp ex[], Glob *globs, int simit) |
void | tabulateValues (Glob *globs, GlobExp *ex, double t0, double t1, double *t, long n, double *state, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp) |
void | outSimit (GlobExp ex[], Glob *globs, double *t, long n, double **y) |
void | freeMem (GlobExp *ex, Glob *globs, int simit) |
void | mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) |
Variables | |
ofstream * | dbg |
debug stream (C++) |
|
Definition at line 38 of file simit_mex.cc. |
|
Definition at line 37 of file simit_mex.cc. |
|
Definition at line 31 of file simit_mex.cc. |
|
Definition at line 35 of file simit_mex.cc. Referenced by mexFunction(). |
|
Definition at line 36 of file simit_mex.cc. Referenced by mexFunction(). |
|
Definition at line 11 of file freeMem.cc. References FALSE, free_d3tensor(), free_dmatrix(), free_dvector(), Glob::gnuFp, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, Glob::nrExp, and GlobExp::nvar. 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 98 of file initialise.cc. References d3tensor(), dbg, dmatrix(), GlobExp::dmdp, GlobExp::dmds, Glob::doP, GlobExp::dR2dp, GlobExp::dR2ds, GlobExp::dR3dp, GlobExp::dR3ds, dvector(), GlobExp::dydp, GlobExp::dyds, GlobExp::errP, GlobExp::errY0, FALSE, GlobExp::h, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, GlobExp::par, GlobExp::r2, GlobExp::residues, TRUE, GlobExp::yComp, GlobExp::yPred, GlobExp::yTry, and GlobExp::yTrySave. 00099 { 00100 long i,j,k,nExp; 00101 long nrExp=globs->nrExp; 00102 00103 ifstream inSplines; 00104 00105 if(simit==FALSE) 00106 { 00107 for (nExp=1;nExp<=nrExp;++nExp) 00108 { 00109 00110 //non-local parameters are equal for each experiment 00111 for(i=1;i<=globs->npar;i++) 00112 { 00113 if(globs->doP[i]==TRUE) 00114 ex[nExp].par[i]=ex[1].par[i]; 00115 } 00116 00117 // set ex->me and ex->mg given in model.cc 00118 setNrConstraints(&ex[nExp],globs->npar,ex[nExp].par); 00119 00120 long nPoints=ex[nExp].nPoints; 00121 long nP=globs->npar; 00122 long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs; 00123 long nMeas=ex[nExp].nMeasure; 00124 long me=ex[nExp].me, mg=ex[nExp].mg; 00125 00126 #ifdef PRINTINITVALUES 00127 *dbg << "initialize: exp. #" << nExp << ", nPoints=" << nPoints 00128 << ", nvar=" << nvar << ", nobs=" << nobs << "\nnMeas=" 00129 << nMeas << ", nVal=" << nVal << ", me=" << me 00130 << ", mg=" << mg << ", nP=" << nP << "\n"; 00131 dbg->flush(); 00132 #endif 00133 //std. error -> output 00134 ex[nExp].errP=dvector(1,nP); 00135 ex[nExp].errY0=dvector(1,nvar); 00136 // initial guesses at mesh points 00137 ex[nExp].yTry=dmatrix(1,nPoints,1,nvar); 00138 ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar); 00139 // computed values at mesh points 00140 ex[nExp].yComp=dmatrix(1,nPoints,1,nvar); 00141 // computed values at measuring points 00142 ex[nExp].yPred=dmatrix(1,nMeas,1,nobs); 00143 00144 ex[nExp].h=dmatrix(1,nPoints,1,nvar); // discrepancies 00145 ex[nExp].residues=dmatrix(1,nMeas,1,nobs); // residues 00146 00147 // derivatives 00148 // d(yComp[i][j]) / d(yTry[i-1][k]) 00149 ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar); 00150 // d(yComp[i][j]) / d(p[k]) 00151 ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP); 00152 // d(yPred[i][j]) / d(yTry[**][k]) 00153 ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar); 00154 // d(yPred[i][j]) / d(p[k]) 00155 ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP); 00156 if (me>0) 00157 { 00158 ex[nExp].r2=dvector(1,me); // equal. constr. 00159 // d(yR2[i]) / d(yTry[j][k]) 00160 ex[nExp].dR2ds=d3tensor(1,me,1,nPoints,1,nvar); 00161 for(i=1;i<=me;i++) 00162 for(j=1;j<=nPoints;j++) 00163 for(k=1;k<=nvar;k++) 00164 ex[nExp].dR2ds[i][j][k]=0.; 00165 00166 // d(yR2[i]) / d(p[j]) 00167 ex[nExp].dR2dp=dmatrix(1,me,1,nP); 00168 for(i=1;i<=me;i++) 00169 for(j=1;j<=nP;j++) 00170 ex[nExp].dR2dp[i][j]=0.; 00171 } 00172 else 00173 { 00174 ex[nExp].r2=NULL; ex[nExp].dR2ds=NULL; ex[nExp].dR2dp=NULL; 00175 } 00176 if (mg>0) 00177 { 00178 ex[nExp].r3=dvector(1,mg); // ineq. constr. 00179 // d(yR3[i]) / d(yTry[j][k]) 00180 ex[nExp].dR3ds=d3tensor(1,mg,1,nPoints,1,nvar); 00181 for(i=1;i<=mg;i++) 00182 for(j=1;j<=nPoints;j++) 00183 for(k=1;k<=nvar;k++) 00184 ex[nExp].dR3ds[i][j][k]=0.; 00185 // d(yR3[i]) / d(p[j]) 00186 ex[nExp].dR3dp=dmatrix(1,mg,1,nP); 00187 for(i=1;i<=mg;i++) 00188 for(j=1;j<=nP;j++) 00189 ex[nExp].dR3dp[i][j]=0.; 00190 } 00191 else 00192 { 00193 ex[nExp].r3=NULL; ex[nExp].dR3ds=NULL; ex[nExp].dR3dp=NULL; 00194 } 00195 00196 // compute yTry from ex (all ex read from file still present) 00197 setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0); 00198 // check if initial values are compatible with constraints 00199 // CURRENTLY DISABLED (1. Dez. 2004) 00200 // if (!initValuesOK(globs,&ex[nExp])) 00201 // { 00202 // cerr << "Experiment #" << nExp 00203 // << ": initial values are not compatible with constraints\n"; 00204 // cerr << "R2:"; 00205 // for (i=1;i<=ex[nExp].me;++i) 00206 // cerr << "\t" << ex[nExp].r2[i]; 00207 // cerr << "\nR3:"; 00208 // for (i=1;i<=ex[nExp].mg;++i) 00209 // cerr << "\t" << ex[nExp].r3[i]; 00210 // cerr << "\n"; 00211 // exit(1); 00212 // } 00213 00214 #ifdef PRINTEX 00215 *dbg << "Ex used for fitting: \n\n"; 00216 for (i=1; i<=ex[nExp].nMeasure; ++i) 00217 { 00218 *dbg << ex[nExp].xMeasure[i]; 00219 for (j=1; j<=ex[nExp].nobs; ++j) 00220 *dbg << "\t" << ex[nExp].yMeasure[i][j] << "\t" << ex[nExp].sigma[i][j]; 00221 *dbg << '\n'; 00222 } 00223 *dbg << '\n'; 00224 dbg->flush(); 00225 #endif 00226 //Allocate memory for condensation 00227 //least squares 00228 ex[nExp].ua=dvector(1,nMeas*nobs); 00229 ex[nExp].Ea=dmatrix(1,nMeas*nobs,1,nvar); 00230 ex[nExp].Pa=dmatrix(1,nMeas*nobs,1,nP); 00231 //equality constraints 00232 ex[nExp].ue=dvector(1,me); 00233 ex[nExp].Ee=dmatrix(1,me,1,nvar); 00234 ex[nExp].Pe=dmatrix(1,me,1,nP); 00235 //inequality constraints 00236 ex[nExp].ug=dvector(1,mg); 00237 ex[nExp].Eg=dmatrix(1,mg,1,nvar); 00238 ex[nExp].Pg=dmatrix(1,mg,1,nP); 00239 00240 //update steps 00241 ex[nExp].dS=dmatrix(1,nPoints,1,nvar); 00242 ex[nExp].dP=dvector(1,nP); 00243 }// matches: for(nExp=1;nExp<=nrExp;++nExp) loop over experiments 00244 //covar in solvLin allocated 00245 globs->covar=NULL; 00246 } // matches with if(simit==FALSE)... 00247 00248 //Initialise GNUplot 00249 if(globs->noGnu==FALSE) 00250 setupGnuFp(globs,ex); 00251 00252 globs->cond=0; 00253 globs->Lambda=1; 00254 00255 char line[1000]; 00256 //Initialise Splines 00257 if(globs->initSpline==TRUE) 00258 { 00259 for (nExp=1;nExp<=nrExp;++nExp) 00260 { 00261 ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00262 ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00263 ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00264 ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*)); 00265 00266 for(i=1;i<=NSPLINES;i++) 00267 { 00268 j=0; 00269 while(ex[nExp].splineFile[i][j]!='\0') 00270 { 00271 line[j]=ex[nExp].splineFile[i][j]; 00272 j++; 00273 } 00274 line[j]='\0'; 00275 inSplines.open(line); 00276 if(inSplines==NULL) 00277 { 00278 cerr << "Cannot open " << ex[nExp].splineFile[i] << ".\n"; 00279 exit(1); 00280 } 00281 inSplines >> ex[nExp].nNodes[i]; 00282 ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00283 ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00284 ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00285 for(j=1;j<= ex[nExp].nNodes[i];j++) 00286 { 00287 inSplines >> ex[nExp].splineNodes[i][j]; 00288 inSplines >> ex[nExp].splineY[i][j]; 00289 inSplines >> ex[nExp].splineGam[i][j]; 00290 } 00291 inSplines.close(); 00292 } 00293 } 00294 } 00295 }
|
|
Definition at line 41 of file simit_mex.cc. References dmatrix(), Glob::doP, Glob::dt, dvector(), Glob::elastic, Glob::eps, Glob::epsilon, FALSE, GlobExp::fitend, GlobExp::fitstart, Glob::gnuFp, IN, initialise(), Glob::initSpline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, Glob::maxstp, Glob::minimp, N, GlobExp::nNodes, GlobExp::nobs, Glob::noGnu, Glob::noMeasurements, Glob::nowait, Glob::npar, Glob::nrExp, GlobExp::nvar, GlobExp::par, Glob::reg, Glob::sig, SPLINE, GlobExp::splineGam, GlobExp::splineNodes, GlobExp::splineY, Glob::stiff, TRUE, Glob::usesig, Glob::wait, GlobExp::y0, and YOUT. 00046 { 00047 long k,i,j,n,ldum; 00048 char outstr[100]; 00049 Glob globs; 00050 GlobExp *ex; 00051 double time,**y,*state; 00052 double t0; //init. time 00053 double *t,*y0,*p; 00054 double **yout,*yout_,**ydum; 00055 double ***dmds,***dmdp,*dmds_,*dmdp_; 00056 double **dyds,**dydp; 00057 double *n_,*spline_; 00058 double *tdum=dvector(1,1); 00059 mxArray *splinedat; 00060 mxArray *tmp; 00061 mxArray *SPLINE; 00062 double *tmp_; 00063 int hidden; 00064 00065 // input -> initial time 00066 tmp=mxGetField(IN,0,"t"); 00067 tmp_=mxGetPr(tmp); 00068 n=mxGetNumberOfElements(tmp); 00069 t=dvector(1,n); 00070 for(i=1;i<=n;i++) 00071 t[i]=tmp_[i-1]; 00072 00073 //allocation of init. state and parameters 00074 y0=dvector(1,NEQNS); 00075 p=dvector(1,NPARAMS); 00076 00077 // input -> initial time 00078 tmp=mxGetField(IN,0,"t0"); 00079 tmp_=mxGetPr(tmp); 00080 t0=tmp_[0]; 00081 00082 // input -> y0 00083 tmp=mxGetField(IN,0,"y0"); 00084 tmp_=mxGetPr(tmp); 00085 if(mxGetNumberOfElements(tmp)!=NEQNS) 00086 { 00087 if(mxGetNumberOfElements(tmp)!=0) 00088 mexWarnMsgTxt("Insufficient number of initial values.\nUse predefined.\n"); 00089 for(i=1;i<=NEQNS;i++) 00090 y0[i]=DefYValues[i-1]; 00091 } 00092 else 00093 { 00094 for(i=1;i<=NEQNS;i++) 00095 y0[i]=tmp_[i-1]; 00096 } 00097 00098 // input -> p 00099 tmp=mxGetField(IN,0,"p"); 00100 tmp_=mxGetPr(tmp); 00101 if(mxGetNumberOfElements(tmp)!=NPARAMS) 00102 { 00103 if(mxGetNumberOfElements(tmp)!=0) 00104 mexWarnMsgTxt("Insufficient number of parameters.\nUse predefined.\n"); 00105 for(i=1;i<=NPARAMS;i++) 00106 p[i]=DefParameters[i-1]; 00107 } 00108 else 00109 { 00110 for(i=1;i<=NPARAMS;i++) 00111 p[i]=tmp_[i-1]; 00112 } 00113 00114 // input -> eps 00115 tmp=mxGetField(IN,0,"eps"); 00116 tmp_=mxGetPr(tmp); 00117 globs.eps=fabs(tmp_[0]); 00118 00119 // input -> integrator 00120 tmp=mxGetField(IN,0,"int"); 00121 tmp_=mxGetPr(tmp); 00122 globs.integrator=abs((int)tmp_[0]); 00123 globs.stiff=TRUE; 00124 // input -> spline 00125 SPLINE=mxGetField(IN,0,"spline"); 00126 00127 // input -> hidden 00128 tmp=mxGetField(IN,0,"hidden"); 00129 tmp_=mxGetPr(tmp); 00130 hidden=(int)tmp_[0]; 00131 00132 //allocation of trajectory vect. for output 00133 if(hidden==TRUE) 00134 YOUT=mxCreateDoubleMatrix(n,NOBS+NEQNS,mxREAL); 00135 else 00136 YOUT=mxCreateDoubleMatrix(n,NOBS,mxREAL); 00137 yout_=mxGetPr(YOUT); 00138 yout=dmatrix(1,n,1,NOBS+NEQNS); 00139 ydum=dmatrix(1,2,1,NOBS+NEQNS); 00140 00141 //initialise some global parameters in globs 00142 globs.noGnu=TRUE; 00143 globs.npar=NPARAMS; 00144 globs.noMeasurements=FALSE; 00145 globs.doP=ivector(1,globs.npar); 00146 for(k=1;k<=globs.npar;k++) 00147 globs.doP[k]=TRUE; 00148 globs.maxit=1000; 00149 globs.gnuFp=NULL; 00150 globs.wait=FALSE; 00151 globs.usesig=FALSE; 00152 globs.stiff=TRUE; 00153 globs.maxstp=5000; 00154 globs.minimp=0.05; 00155 globs.nowait=FALSE; 00156 globs.elastic=1.; 00157 globs.reg=FALSE; 00158 globs.epsilon=1e-10; 00159 globs.lambda=1e6; 00160 globs.dt=0.1; 00161 globs.sig=0.; 00162 globs.nrExp=1; 00163 00164 ex=new GlobExp[globs.nrExp+1]; 00165 //initialise some global parameters 00166 00167 ex[1].nobs=NOBS; 00168 ex[1].nvar=NEQNS; 00169 ex[1].par=dvector(1,NPARAMS); 00170 ex[1].y0=dvector(1,NEQNS); 00171 00172 00173 for(k=1;k<=NPARAMS;k++) 00174 ex[1].par[k]=p[k]; 00175 00176 00177 ex[1].fitstart=t[1]; 00178 ex[1].fitend=t[n]; 00179 //initial values 00180 00181 state=dvector(1,ex[1].nvar); 00182 for(k=1;k<=NEQNS;k++) 00183 state[k]=y0[k]; 00184 00185 00186 N=mxCreateDoubleMatrix(3,1,mxREAL); 00187 n_=mxGetPr(N); 00188 00189 n_[0]=NEQNS; 00190 n_[1]=NOBS; 00191 n_[2]=NPARAMS; 00192 00193 // spline section 00194 if(NSPLINES==0) 00195 { 00196 globs.initSpline=TRUE; 00197 initialise(ex,&globs,TRUE); 00198 } 00199 else 00200 { 00201 globs.initSpline=FALSE; 00202 initialise(ex,&globs,TRUE); 00203 ex[1].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00204 ex[1].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00205 ex[1].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00206 ex[1].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*)); 00207 00208 for(i=1;i<=NSPLINES;i++) 00209 { 00210 if(mxGetNumberOfElements(SPLINE) < NSPLINES) 00211 { 00212 cerr << "Please specify " << NSPLINES << " spline(s).\n"; 00213 return; 00214 } 00215 00216 splinedat=mxGetCell(SPLINE,i-1); 00217 spline_=mxGetPr(splinedat); 00218 ex[1].nNodes[i]=mxGetM(splinedat); 00219 if(mxGetN(splinedat)!=3) 00220 { 00221 cerr << "Invalid data for spline " << i << ".\n"; 00222 return; 00223 } 00224 ex[1].splineNodes[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double)); 00225 ex[1].splineY[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double)); 00226 ex[1].splineGam[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double)); 00227 00228 //copy spline data 00229 for(j=1;j<=ex[1].nNodes[i];j++) 00230 { 00231 ex[1].splineNodes[i][j]=spline_[j-1]; 00232 ex[1].splineY[i][j]=spline_[ex[1].nNodes[i]+j-1]; 00233 ex[1].splineGam[i][j]=spline_[2*ex[1].nNodes[i]+j-1]; 00234 } 00235 } 00236 } 00237 00238 //... and integrate 00239 try 00240 { 00241 if(nlhs!=2) 00242 { 00243 dmds=d3tensor(1,n,1,NEQNS,1,NEQNS); 00244 dmdp=d3tensor(1,n,1,NEQNS,1,NPARAMS); 00245 dyds=dmatrix(1,NEQNS,1,NEQNS); 00246 dydp=dmatrix(1,NEQNS,1,NPARAMS); 00247 00248 tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,dmds,dmdp,dyds,dydp); 00249 00250 DMDS=mxCreateDoubleMatrix(NOBS*NEQNS,n,mxREAL); 00251 dmds_=mxGetPr(DMDS); 00252 DMDP=mxCreateDoubleMatrix(NOBS*NPARAMS,n,mxREAL); 00253 dmdp_=mxGetPr(DMDP); 00254 00255 00256 //copying 00257 00258 for(k=1;k<=NEQNS;k++) 00259 for(i=1;i<=NOBS;i++) 00260 for(j=1;j<=n;j++) 00261 { 00262 dmds_[(j-1)*NOBS*NEQNS+(k-1)*NOBS+i-1]=dmds[j][i][k]; 00263 } 00264 00265 00266 for(k=1;k<=NPARAMS;k++) 00267 for(i=1;i<=NOBS;i++) 00268 for(j=1;j<=n;j++) 00269 { 00270 dmdp_[(j-1)*NOBS*NPARAMS+(k-1)*NOBS+i-1]=dmdp[j][i][k]; 00271 } 00272 00273 free_d3tensor(dmds,1,n,1,NEQNS,1,NEQNS); 00274 free_d3tensor(dmdp,1,n,1,NEQNS,1,NPARAMS); 00275 free_dmatrix(dyds,1,NEQNS,1,NEQNS); 00276 free_dmatrix(dydp,1,NEQNS,1,NPARAMS); 00277 } 00278 else 00279 { 00280 tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,NULL,NULL,NULL,NULL); 00281 00282 } 00283 00284 if(hidden==TRUE) 00285 { 00286 ldum=NOBS+NEQNS; 00287 for(i=1;i<=NEQNS;i++) 00288 state[i]=y0[i]; 00289 tdum[1]=t[1]; 00290 if(t0!=t[1]) 00291 tabulateValues(&globs,&ex[1],t0,t[1],tdum,1,state,ydum,NULL,NULL,NULL,NULL); 00292 for(j=1;j<=NEQNS;j++) 00293 yout[1][NOBS+j]=state[j]; 00294 00295 for(i=2;i<=n;i++) 00296 { 00297 tdum[1]=t[i]; 00298 tabulateValues(&globs,&ex[1],t[i-1],t[i],tdum,1,state,ydum,NULL,NULL,NULL,NULL); 00299 for(j=1;j<=NEQNS;j++) 00300 yout[i][NOBS+j]=state[j]; 00301 00302 } 00303 } 00304 else 00305 ldum=NOBS; 00306 00307 for(i=1;i<=ldum;i++) 00308 for(j=1;j<=n;j++) 00309 yout_[(i-1)*n+j-1]=yout[j][i]; 00310 } 00311 catch (int i) 00312 { 00313 return; 00314 } 00315 //free mem 00316 freeMem(ex,&globs,TRUE); 00317 free_dvector(state,1,ex[1].nvar); 00318 free_dmatrix(yout,1,n,1,NEQNS+NOBS); 00319 free_dvector(t,1,n); 00320 free_dvector(y0,1,NEQNS); 00321 free_dvector(p,1,NPARAMS); 00322 free_dvector(ex[1].par,1,NPARAMS); 00323 free_dvector(ex[1].y0,1,NEQNS); 00324 free_dvector(tdum,1,1); 00325 free_dmatrix(ydum,1,2,1,NOBS+NEQNS); 00326 00327 delete ex; 00328 }
|
|
Definition at line 12 of file outSimit.cc. References GlobExp::nobs, and Glob::sig. 00013 { 00014 long k,l; 00015 long idum=-time(NULL); 00016 long nobs=ex[1].nobs; 00017 ofstream out; 00018 00019 out.open(ex[1].fileName); 00020 00021 if(globs->sig==0) 00022 { 00023 for(k=1;k<=n;k++) 00024 { 00025 out << t[k]; 00026 for(l=1;l<=nobs;l++) 00027 out << " " << y[k][l] << " 1"; 00028 out << endl; 00029 } 00030 00031 } 00032 else 00033 { 00034 double *max=dvector(1,nobs),*min=dvector(1,nobs); 00035 for(l=1;l<=nobs;l++) 00036 { 00037 max[l]=y[1][l]; 00038 min[l]=y[1][l]; 00039 } 00040 for(k=2;k<=n;k++) 00041 { 00042 for(l=1;l<=nobs;l++) 00043 { 00044 if(max[l]<y[k][l]) 00045 max[l]=y[k][l]; 00046 if(min[l]>y[k][l]) 00047 min[l]=y[k][l]; 00048 } 00049 } 00050 for(k=1;k<=n;k++) 00051 { 00052 out << t[k] ; 00053 for(l=1;l<=nobs;l++) 00054 out << " " << y[k][l]+globs->sig*fabs(max[l]-min[l])*gasdev(&idum) << " " << globs->sig*fabs(max[l]-min[l]); 00055 out << endl; 00056 } 00057 } 00058 00059 }
|
|
Main module routine doing all the pasing.
Definition at line 120 of file parse.cc. References Glob::doP, Glob::dt, Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob::gnuFp, Glob::initSpline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, MAXLENGTH, Glob::maxstp, Glob::minimp, Glob::noGnu, Glob::noMeasurements, Glob::nowait, Glob::npar, Glob::nrExp, Glob::reg, Glob::sig, Glob::stiff, TRUE, Glob::usesig, and Glob::wait. 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 108 of file intODE.cc. 00110 { 00111 // integrate from t0 to t1, storing values at t[1..n] 00112 // t[n] may be ==t1 or <t1 00113 // in the latter case we will integrate from t[n] to t1 in the 00114 // final loop k=n+1 without storing any data points 00115 // If dmds!=NULL, compute derivatives as well 00116 // and store them in dmds,dmdp,dyds,dydp 00117 // with respect to observation, see odin.doc: Observation function 00118 00119 //BEGIN von Felix veraendert 00120 // double hmin=max(fabs(t0),fabs(t1))*1e-8; 00121 double hmin = max (fabs (t0), fabs (t1)) * 1e-15; 00122 // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest. 00123 // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern 00124 // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion. 00125 // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht 00126 // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich 00127 // von 1e-18 noch einigermassen weit entfernt bleibt. 00128 // Man koennte das in Zukunft noch in einen externen Parameter stecken. 00129 //END von Felix veraendert 00130 00131 long nPExp=globs->npar; 00132 long nVar=ex->nvar; 00133 long nobs=NOBS; 00134 int generic; 00135 double *p=ex->par; 00136 double eps=globs->eps; 00137 double h1 = (t1 - t0) * .01; 00138 long nok,nbad; 00139 00140 initInt(globs,ex); 00141 if (n > 0 && (t0 > t[1] || t1 < t[n])) 00142 { 00143 cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl; 00144 throw 1; 00145 } 00146 00147 if (globs->integrator==3) 00148 { // said -odessa 00149 #ifdef ODESSA 00150 double &rtol = eps; 00151 //BEGIN von Felix veraendert 00152 // double atol=eps*1e-6; 00153 double atol = eps; 00154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!) 00155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu 00156 // siehe odessa.f: ATOL und RTOL 00157 //END von Felix veraendert 00158 // odeint uses TINY=1e-30 in a similar way as odessa's atol, 00159 // but atol<1e-16 or so results in step size underflow 00160 00161 int odeint_MAXSTP=globs->maxstp; 00162 int d_flag=1; 00163 int stif=globs->stiff; 00164 char doPStr[globs->npar+1]; 00165 00166 for(long i=1;i<=globs->npar;i++) 00167 doPStr[i]=(char)globs->doP[i]; 00168 00169 call_odessa(globs,ex,nVar, nPExp, doPStr, 00170 d_flag, odeint_MAXSTP, stif, 00171 TRUE, n, t, t0, t1, state, y, 00172 dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1); 00173 00174 return; 00175 #else 00176 cerr << "-odessa used but not linked" << endl; 00177 throw 1; 00178 #endif // ODESSA 00179 } // if odessa 00180 double *gy=dvector(1,nobs); 00181 long i, j, k; 00182 double xx1 = t0, xx2; 00183 00184 if (dmds) 00185 { // compute sensitivities 00186 // create extended state vector Y 00187 double **Yt = dmatrix (0, nPExp + nVar, 1, nVar); 00188 initYt (Yt[0], state); 00189 00190 // chain rule for extra observables 00191 double **dgdy, **dgdp; 00192 00193 dgdy = dmatrix (1, nobs, 1, nVar); 00194 dgdp = dmatrix (1, nobs, 1, nPExp); 00195 dfill0 (dgdy[1], nobs * nVar); 00196 dfill0 (dgdp[1], nobs * nPExp); 00197 00198 00199 for (k = 1; k <= n + 1; k++) 00200 { 00201 xx2 = (k > n) ? t1 : // end point 00202 t[k]; // data point 00203 00204 globs->rkqs_ign = (nPExp + nVar) * nVar; 00205 if (xx2 > xx1 && nVar > 0) 00206 { 00207 if(globs->integrator==1) 00208 { 00209 //Runge-Kutta integation 00210 integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2, 00211 eps, h1, hmin, &nok, &nbad, sensDerivs); 00212 } 00213 else if(globs->integrator==2) 00214 { 00215 cerr << "CVODES not optimized skip integrator" << endl; 00216 throw 1; 00217 00218 //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE); 00219 } 00220 } 00221 // write state to y and sensitivities to dm/dp and dm/ds 00222 if (k <= n) 00223 { 00224 dcopy (y[k], Yt[0], nVar); 00225 generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp); 00226 for(i=1;i<=nobs;i++) 00227 y[k][i]=gy[i]; 00228 00229 for (i = 1; i <= nobs; i++) 00230 { 00231 for (j = 1; j <= nPExp; j++) 00232 { 00233 double &dest = dmdp[k][i][j]; 00234 if (!generic) 00235 { 00236 dest = dgdp[i][j]; 00237 for (long l = 1; l <= nVar; l++) 00238 dest += dgdy[i][l] * Yt[j][l]; 00239 } 00240 else 00241 dest = Yt[j][i]; 00242 } // for j 00243 for (j = 1; j <= nVar; j++) 00244 { 00245 double &dest = dmds[k][i][j]; 00246 if (!generic) 00247 { 00248 dest = 0; 00249 for (long l = 1; l <= nVar; l++) 00250 dest += dgdy[i][l] * Yt[nPExp + j][l]; 00251 } 00252 else 00253 dest = Yt[nPExp + j][i]; 00254 } // for j 00255 } // for i 00256 } // if k<=n 00257 xx1 = xx2; 00258 } // for k 00259 00260 // copy state and its derivatives to state, dydp and dyds 00261 dcopy (state, Yt[0], nVar); // original state vector 00262 // write sensitivities to dy/dp and dy/ds 00263 for (i = 1; i <= nVar; i++) 00264 { 00265 for (j = 1; j <= nPExp; j++) 00266 dydp[i][j] = Yt[j][i]; 00267 for (j = 1; j <= nVar; j++) 00268 dyds[i][j] = Yt[nPExp + j][i]; 00269 } 00270 00271 free_dmatrix (dgdp,1 ,nobs , 1, nPExp); 00272 free_dmatrix (dgdy,1 ,nobs, 1, nVar); 00273 free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar); 00274 } 00275 else 00276 { // !dmds 00277 for (k = 1; k <= n + 1; k++) 00278 { 00279 xx2 = (k <= n) ? t[k] : t1; 00280 00281 globs->rkqs_ign = 0; 00282 if (xx2 > xx1 && nVar > 0) 00283 { 00284 if(globs->integrator==1) 00285 { 00286 //Runge-Kutta integration 00287 integrateRK(globs,state, p, nVar, xx1, xx2, 00288 eps, h1, hmin, &nok, &nbad, ode); 00289 } 00290 else if(globs->integrator==2) 00291 { 00292 cerr << "CVODES not optimized skip integrator" << endl; 00293 throw 1; 00294 00295 //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE); 00296 } 00297 } 00298 if (k <= n) 00299 { 00300 dcopy (y[k], state, nVar); 00301 generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL); 00302 for(i=1;i<=nobs;i++) 00303 y[k][i]=gy[i]; 00304 } 00305 xx1 = xx2; 00306 } // for k 00307 } 00308 // if dmds else 00309 free_dvector(gy,1,nobs); 00310 }
|
|
debug stream (C++)
Definition at line 29 of file simit_mex.cc. Referenced by initialise(), main(), and objectiveF(). |