#include "mex.h"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string>
#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 | INEX prhs[1] |
#define | OUT plhs[0] |
#define | OUTEX plhs[1] |
Functions | |
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 | mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) |
Variables | |
ofstream * | dbg |
debug stream (C++) |
|
Definition at line 34 of file diffit_mex.cc. Referenced by mexFunction(). |
|
Definition at line 35 of file diffit_mex.cc. |
|
Definition at line 37 of file diffit_mex.cc. |
|
Definition at line 38 of file diffit_mex.cc. |
|
Definition at line 33 of file fitIt.cc. References Glob::minimp, and Glob::nrExp. 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 98 of file initialise.cc. 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 diffit_mex.cc. References Glob::doP, Glob::elastic, Glob::eps, Glob::epsilon, FALSE, Glob::gnuFp, IN, Glob::initSpline, Glob::integrator, ivector(), Glob::lambda, Glob::maxit, Glob::maxstp, Glob::minimiser, Glob::minimp, Glob::nodamp, Glob::noGnu, Glob::noMeasurements, Glob::nowait, Glob::npar, Glob::nrExp, Glob::pert, Glob::reg, Glob::silent, Glob::simInit, Glob::stiff, TRUE, Glob::wait, and Glob::wquer. 00045 { 00046 long k,i,j,nExp; 00047 long index,nglob=0; 00048 unsigned long *ndat; 00049 Glob globs; 00050 GlobExp *ex; 00051 mxArray *tmp,*sig; 00052 double *tmp_,*sig_,*spline_; 00053 char tmpstr[1000]; 00054 mxArray *splinedat; 00055 int except=0; 00056 int init=0; 00057 00058 //for output 00059 const char *fieldNamesOUT[]={"wquer","converged","chisq","except","Lambda","cond"}; 00060 const char *fieldNamesOUTEX[]={"p","y0","mesh","errP","errY0"}; 00061 00062 int dims[2]; 00063 00064 //GET DATA FOR GLOBS 00065 00066 globs.npar=NPARAMS; 00067 globs.wait=FALSE; 00068 globs.nowait=TRUE; 00069 globs.elastic=1.; 00070 globs.initSpline=FALSE; 00071 globs.maxstp=5000; 00072 globs.gnuFp=NULL; 00073 globs.initSpline=TRUE; 00074 globs.wquer=100.; 00075 globs.silent=FALSE; 00076 globs.minimiser=2; 00077 00078 //eps 00079 tmp=mxGetField(IN,0,"eps"); 00080 tmp_=mxGetPr(tmp); 00081 globs.eps=fabs(tmp_[0]); 00082 //eps 00083 tmp=mxGetField(IN,0,"nognu"); 00084 tmp_=mxGetPr(tmp); 00085 globs.noGnu=(int)tmp_[0]; 00086 //nomeasure 00087 tmp=mxGetField(IN,0,"nomeasure"); 00088 tmp_=mxGetPr(tmp); 00089 globs.noMeasurements=(int)tmp_[0]; 00090 //maxit 00091 tmp=mxGetField(IN,0,"maxit"); 00092 tmp_=mxGetPr(tmp); 00093 globs.maxit=abs((int)tmp_[0]); 00094 //odessa 00095 tmp=mxGetField(IN,0,"int"); 00096 tmp_=mxGetPr(tmp); 00097 globs.integrator=abs((int)tmp_[0]); 00098 globs.stiff=TRUE; 00099 //minimp 00100 tmp=mxGetField(IN,0,"minimp"); 00101 tmp_=mxGetPr(tmp); 00102 globs.minimp=fabs(tmp_[0]); 00103 //reg 00104 tmp=mxGetField(IN,0,"reg"); 00105 tmp_=mxGetPr(tmp); 00106 globs.reg=(int)tmp_[0]; 00107 //epsilon 00108 tmp=mxGetField(IN,0,"epsilon"); 00109 tmp_=mxGetPr(tmp); 00110 globs.epsilon=fabs(tmp_[0]); 00111 //lambda 00112 tmp=mxGetField(IN,0,"lambda"); 00113 tmp_=mxGetPr(tmp); 00114 globs.lambda=fabs(tmp_[0]); 00115 //siminit 00116 tmp=mxGetField(IN,0,"siminit"); 00117 tmp_=mxGetPr(tmp); 00118 globs.simInit=(int)tmp_[0]; 00119 //pert 00120 tmp=mxGetField(IN,0,"pert"); 00121 tmp_=mxGetPr(tmp); 00122 globs.pert=fabs(tmp_[0]); 00123 //nodamp 00124 tmp=mxGetField(IN,0,"nodamp"); 00125 tmp_=mxGetPr(tmp); 00126 globs.nodamp=(int)tmp_[0]; 00127 //nrExp 00128 tmp=mxGetField(IN,0,"nexp"); 00129 tmp_=mxGetPr(tmp); 00130 globs.nrExp=(long)tmp_[0]; 00131 //doP 00132 globs.doP=ivector(1,globs.npar); 00133 for(k=1;k<=globs.npar;k++) 00134 globs.doP[k]=TRUE; 00135 tmp=mxGetField(IN,0,"doP"); 00136 if(mxGetString(tmp,tmpstr,1000)) 00137 { 00138 cerr << "Parsing doP failed.\n"; 00139 return; 00140 } 00141 if(tmpstr[0]!='N') 00142 { 00143 if(strlen(tmpstr)!=globs.npar) 00144 { 00145 cerr << "doP must be of length " << globs.npar << ".\n"; 00146 return; 00147 } 00148 else 00149 { 00150 for(k=0;k < globs.npar;k++) 00151 { 00152 if(tmpstr[k]=='0') 00153 globs.doP[k+1]=FALSE; 00154 else if (tmpstr[k]=='L') 00155 globs.doP[k+1]='L'; 00156 else 00157 globs.doP[k+1]=TRUE; 00158 } 00159 } 00160 } 00161 //y0fix 00162 globs.y0fix=ivector(1,NEQNS); 00163 for(k=1;k<=NEQNS;k++) 00164 globs.y0fix[k]=TRUE; 00165 tmp=mxGetField(IN,0,"y0fix"); 00166 if(mxGetString(tmp,tmpstr,1000)) 00167 { 00168 cerr << "Parsing y0fix failed.\n"; 00169 return; 00170 } 00171 if(tmpstr[0]!='N') 00172 { 00173 if(strlen(tmpstr)!=NEQNS) 00174 { 00175 cerr << "y0fix must be of length " << NEQNS << ".\n"; 00176 return; 00177 } 00178 else 00179 { 00180 for(k=0;k < NEQNS ;k++) 00181 { 00182 if(tmpstr[k]=='0') 00183 globs.y0fix[k+1]=FALSE; 00184 else 00185 globs.y0fix[k+1]=TRUE; 00186 } 00187 } 00188 } 00189 //wquer 00190 tmp=mxGetField(IN,0,"wquer"); 00191 tmp_=mxGetPr(tmp); 00192 globs.wquer=tmp_[0]; 00193 //init 00194 tmp=mxGetField(IN,0,"init"); 00195 tmp_=mxGetPr(tmp); 00196 init=(int)tmp_[0]; 00197 //silent 00198 tmp=mxGetField(IN,0,"silent"); 00199 tmp_=mxGetPr(tmp); 00200 globs.silent=(int)tmp_[0]; 00201 //opt 00202 tmp=mxGetField(IN,0,"opt"); 00203 tmp_=mxGetPr(tmp); 00204 globs.minimiser=(int)tmp_[0]; 00205 //Lconst 00206 tmp=mxGetField(IN,0,"Lconst"); 00207 tmp_=mxGetPr(tmp); 00208 if(mxIsEmpty(tmp)) 00209 { 00210 globs.faktorLexist=FALSE; 00211 } 00212 else 00213 { 00214 globs.faktorL=dvector(1,NPARAMS); 00215 for(k=1;k<=NPARAMS;k++) 00216 globs.faktorL[k]=-1; 00217 for(k=1;k<=mxGetNumberOfElements(tmp) && k<=NPARAMS;k++) 00218 { 00219 globs.faktorL[k]=tmp_[k-1]; 00220 } 00221 globs.faktorLexist=TRUE; 00222 } 00223 00224 00225 00226 // GET DATA FOR EX 00227 ndat=lvector(1,globs.nrExp); 00228 ex=new GlobExp[globs.nrExp+1]; 00229 00230 for(nExp=1;nExp<=globs.nrExp;nExp++) 00231 { 00232 //general stuff 00233 ex[nExp].nvar=NEQNS; 00234 ex[nExp].splinesDefined=FALSE; 00235 //t0 00236 tmp=mxGetField(INEX,nExp-1,"t0"); 00237 tmp_=mxGetPr(tmp); 00238 ex[nExp].fitstart=tmp_[0]; 00239 00240 //t1 00241 tmp=mxGetField(INEX,nExp-1,"t1"); 00242 tmp_=mxGetPr(tmp); 00243 ex[nExp].fitend=tmp_[0]; 00244 if(ex[nExp].fitstart >= ex[nExp].fitend) 00245 { 00246 cerr << "Fitstart >= fitend in exp. " << nExp << ".\n"; 00247 return; 00248 } 00249 //nobs 00250 tmp=mxGetField(INEX,nExp-1,"nobs"); 00251 tmp_=mxGetPr(tmp); 00252 ex[nExp].nobs=(long)tmp_[0]; 00253 if(ex[nExp].nobs > NOBS) 00254 { 00255 cerr << "More observables than required in exp. " << nExp << ".\n"; 00256 return; 00257 } 00258 //nms 00259 tmp=mxGetField(INEX,nExp-1,"nms"); 00260 tmp_=mxGetPr(tmp); 00261 ex[nExp].nPoints=(long)tmp_[0]+1; 00262 //p 00263 ex[nExp].par=dvector(1,NPARAMS); 00264 tmp=mxGetField(INEX,nExp-1,"p"); 00265 tmp_=mxGetPr(tmp); 00266 if(mxIsEmpty(tmp)) 00267 { 00268 for(k=1;k<=globs.npar;k++) 00269 ex[nExp].par[k]=DefParameters[k-1]; 00270 } 00271 else 00272 { 00273 if(mxGetNumberOfElements(tmp)!=globs.npar) 00274 { 00275 cerr << globs.npar << " parameter(s) required in exp. " << nExp <<".\n"; 00276 return; 00277 } 00278 else 00279 { 00280 for(k=1;k<=globs.npar;k++) 00281 ex[nExp].par[k]= tmp_[k-1]; 00282 } 00283 } 00284 //y0 00285 tmp=mxGetField(INEX,nExp-1,"y0"); 00286 tmp_=mxGetPr(tmp); 00287 if(mxIsEmpty(tmp)) 00288 { 00289 ex[nExp].y0=NULL; 00290 } 00291 else 00292 { 00293 if(mxGetNumberOfElements(tmp)!=ex[nExp].nvar) 00294 { 00295 cerr << ex[nExp].nvar << " initial value(s) required in exp. " << nExp <<".\n"; 00296 return; 00297 } 00298 else 00299 { 00300 ex[nExp].y0=dvector(1,ex[nExp].nvar); 00301 for(k=1;k<=ex[nExp].nvar;k++) 00302 ex[nExp].y0[k]= tmp_[k-1]; 00303 } 00304 } 00305 //n 00306 tmp=mxGetField(INEX,nExp-1,"n"); 00307 tmp_=mxGetPr(tmp); 00308 ndat[nExp]=(long)tmp_[0]; 00309 00310 } //end for(nExp=1;nExp<=globs.nrExp;nExp++) 00311 00312 //open DEBUG stream 00313 dbg=new ofstream("diffit.dbg"); 00314 if (!dbg) 00315 { 00316 cerr << "Error opening DEBUG file.\n"; 00317 return; 00318 } 00319 dbg->precision(4); // set output format for floating point data 00320 dbg->unsetf(ios::floatfield); 00321 00322 //print debugging information 00323 *dbg << DefModelDescription << "\n"; 00324 if (globs.nrExp==1) 00325 *dbg << "1 experiment\n"; 00326 else 00327 *dbg << globs.nrExp << " experiments\n"; 00328 *dbg << "\n"; 00329 for(i=1;i<=globs.nrExp;i++) 00330 { 00331 *dbg << "Experiment: " << i << "\n"; 00332 *dbg << NPARAMS << " parameter(s):"; 00333 for (k=1; k<=NPARAMS; ++k) 00334 *dbg << " " << ex[i].par[k]; 00335 *dbg << "\n"; 00336 } 00337 00338 dbg->flush(); 00339 00340 00341 for(i=1;i<=globs.nrExp;i++) 00342 { 00343 //read data and check for for consistency 00344 00345 //data 00346 tmp=mxGetField(INEX,i-1,"data"); 00347 tmp_=mxGetPr(tmp); 00348 sig=mxGetField(INEX,i-1,"sig"); 00349 sig_=mxGetPr(sig); 00350 //deternine nMeasure 00351 index=1; 00352 for(j=1;j<=ndat[i];j++) 00353 { 00354 if(tmp_[j-1]>= ex[i].fitstart && tmp_[j-1]<= ex[i].fitend) 00355 index++; 00356 } 00357 index--; 00358 if(index==0) 00359 { 00360 cerr << "Nothing to fit in exp " << i << ".\n"; 00361 return; 00362 } 00363 ex[i].nMeasure=index; 00364 //memory allocation 00365 ex[i].xMeasure=dvector(1,ex[i].nMeasure); 00366 ex[i].yMeasure=dmatrix(1,ex[i].nMeasure,1,ex[i].nobs); 00367 ex[i].sigma=dmatrix(1,ex[i].nMeasure,1,ex[i].nobs); 00368 //copy data 00369 index=1; 00370 for(j=1;j<=ndat[i];j++) 00371 { 00372 if(tmp_[j-1]>= ex[i].fitstart && tmp_[j-1]<= ex[i].fitend) 00373 { 00374 ex[i].xMeasure[index]=tmp_[j-1]; 00375 for(k=1;k<=ex[i].nobs;k++) 00376 { 00377 ex[i].yMeasure[index][k]=tmp_[ndat[i]*k+j-1]; 00378 ex[i].sigma[index][k]=sig_[ndat[i]*(k-1)+j-1]; 00379 } 00380 index++; 00381 } 00382 } 00383 //check if time is in ascending order 00384 for(j=2;j<=ex[i].nMeasure;j++) 00385 { 00386 if(ex[i].xMeasure[j-1] >= ex[i].xMeasure[j]) 00387 { 00388 cerr << "Data not in temporal ascending order in exp " << i << ",\n"; 00389 return; 00390 } 00391 } 00392 //setting some variables 00393 ex[i].firstMeasure=1; 00394 ex[i].lastMeasure=ex[i].nMeasure; 00395 00396 //set mesh 00397 tmp=mxGetField(INEX,i-1,"mesh"); 00398 tmp_=mxGetPr(tmp); 00399 if(mxIsEmpty(tmp)) 00400 setMesh(ex,&globs,i); 00401 else 00402 { 00403 ex[i].nPoints=mxGetM(tmp)+2; 00404 if (ex[i].nMeasure < ex[i].nPoints) 00405 { 00406 cerr << "Too few measurements to construct mesh\n"; 00407 return; 00408 } 00409 00410 if(mxGetN(tmp) <= ex[i].nvar+1) 00411 { 00412 ex[i].mesh=dvector(1,ex[i].nPoints); 00413 ex[i].mesh[1]= ex[i].fitstart; 00414 ex[i].mesh[ex[i].nPoints]= ex[i].fitend; 00415 //reading mesh data 00416 for(j=1;j<=ex[i].nPoints-2;j++) 00417 { 00418 ex[i].mesh[j+1]=tmp_[j-1]; 00419 } 00420 //check if mesh is in ascending order 00421 for(j=2;j<=ex[i].nPoints;j++) 00422 { 00423 if(ex[i].mesh[j-1] >= ex[i].mesh[j]) 00424 { 00425 cerr << "Mesh not in ascending order in exp " << i << ",\n"; 00426 return; 00427 } 00428 } 00429 } 00430 else 00431 { 00432 cerr << "Incorrect format of mesh in exp " << i << ",\n"; 00433 return; 00434 } 00435 } 00436 00437 #ifdef PRINTDATA 00438 *dbg << "\nData: \n"; 00439 for (j=1; j<=ex[i].nMeasure; j++) { 00440 *dbg << ex[i].xMeasure[j]; 00441 for (k=1; k<=ex[i].nobs; k++) 00442 *dbg << "\t" << ex[i].yMeasure[j][k] << "\t" << ex[i].sigma[j][k]; 00443 *dbg << "\n"; 00444 } 00445 *dbg << "\n"; 00446 #endif 00447 *dbg << "Mesh:"; 00448 for (k=1;k<=ex[i].nPoints;++k) 00449 *dbg << " " << ex[i].mesh[k]; 00450 *dbg << "\n"; 00451 } 00452 00453 //initial values 00454 for(i=1;i<=globs.nrExp;i++) 00455 { 00456 *dbg << "\nExperiment " << i << ":\n"; 00457 if (ex[i].y0) 00458 { 00459 *dbg << NEQNS << " starting value(s):"; 00460 for (k=1; k<=NEQNS; ++k) 00461 *dbg << ' ' << ex[i].y0[k]; 00462 *dbg << "\n"; 00463 } 00464 else 00465 { 00466 *dbg << "No starting values specified\n"; 00467 } 00468 *dbg << "\n"; 00469 dbg->flush(); 00470 } // end for(i=1;i<=globs.nrExp;i++) 00471 // starting of the numerics 00472 00473 00474 // spline section and initialse 00475 00476 if(NSPLINES==0) 00477 { 00478 globs.initSpline=TRUE; 00479 initialise(ex,&globs,FALSE); 00480 } 00481 else 00482 { 00483 for(nExp=1;nExp<=globs.nrExp;nExp++) 00484 { 00485 globs.initSpline=FALSE; 00486 initialise(ex,&globs,FALSE); 00487 ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00488 ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00489 ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00490 ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*)); 00491 tmp=mxGetField(INEX,nExp-1,"spline"); 00492 for(i=1;i<=NSPLINES;i++) 00493 { 00494 if(mxGetNumberOfElements(tmp) < NSPLINES) 00495 { 00496 cerr << "Please specify " << NSPLINES << " spline(s).\n"; 00497 return; 00498 } 00499 00500 splinedat=mxGetCell(tmp,i-1); 00501 spline_=mxGetPr(splinedat); 00502 ex[nExp].nNodes[i]=mxGetM(splinedat); 00503 if(mxGetN(splinedat)!=3) 00504 { 00505 cerr << "Invalid data for spline " << i << ".\n"; 00506 return; 00507 } 00508 ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00509 ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00510 ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00511 00512 //copy spline data 00513 for(j=1;j<=ex[nExp].nNodes[i];j++) 00514 { 00515 ex[nExp].splineNodes[i][j]=spline_[j-1]; 00516 ex[nExp].splineY[i][j]=spline_[ex[nExp].nNodes[i]+j-1]; 00517 ex[nExp].splineGam[i][j]=spline_[2*ex[nExp].nNodes[i]+j-1]; 00518 } 00519 } 00520 } 00521 } 00522 00523 //simulate initial state 00524 try 00525 { 00526 if(globs.simInit==TRUE) 00527 { 00528 for (i=1;i<=globs.nrExp;++i) 00529 simInit(&ex[i],&globs); 00530 } 00531 } 00532 catch (int i) 00533 { 00534 cerr << "Cannot integrate trajectory, no siminit!\n"; 00535 } 00536 00537 //setting initial values from mesh data 00538 for(nExp=1;nExp<=globs.nrExp;nExp++) 00539 { 00540 tmp=mxGetField(INEX,nExp-1,"mesh"); 00541 tmp_=mxGetPr(tmp); 00542 if(!mxIsEmpty(tmp)) 00543 { 00544 index=mxGetN(tmp); 00545 for(i=1;i<=ex[nExp].nPoints-2;i++) 00546 { 00547 for(j=1;j<=index-1;j++) 00548 { 00549 ex[nExp].yTry[i+1][j]=tmp_[j*(ex[nExp].nPoints-2)+i-1]; 00550 } 00551 } 00552 } 00553 } 00554 00555 //************************* 00556 try 00557 { 00558 00559 if(init!=1) 00560 fitIt(ex,&globs); 00561 else 00562 globs.fitConverged=0; 00563 } 00564 catch(int i) //Exception handling 00565 { 00566 dims[0]=1; 00567 dims[1]=1; 00568 OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT); 00569 //error 00570 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00571 tmp_= mxGetPr(tmp); 00572 tmp_[0]= 1; 00573 mxSetField(OUT,0,"except",tmp); 00574 00575 dims[0]=1; 00576 dims[1]=globs.nrExp; 00577 OUTEX = mxCreateStructArray(2,dims,sizeof(fieldNamesOUTEX)/sizeof(*fieldNamesOUTEX),fieldNamesOUTEX); 00578 if(globs.noGnu==FALSE) 00579 { 00580 for(j=1;j<=globs.ngnu;j++) 00581 pclose(globs.gnuFp[j]); 00582 } 00583 00584 freeMem(ex,&globs,FALSE); 00585 delete ex; 00586 return; 00587 } 00588 00589 //************************* 00590 00591 //Output after convergence 00592 if(init!=1) 00593 { 00594 outFit(ex,&globs); 00595 if(!globs.noGnu) 00596 system("rm gnuout.dat"); 00597 } 00598 00599 //Mex Output 00600 00601 dims[0]=1; 00602 dims[1]=1; 00603 00604 00605 OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT); 00606 00607 //wquer 00608 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00609 tmp_= mxGetPr(tmp); 00610 tmp_[0]= globs.wquer; 00611 mxSetField(OUT,0,"wquer",tmp); 00612 //converged 00613 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00614 tmp_= mxGetPr(tmp); 00615 tmp_[0]= globs.fitConverged; 00616 mxSetField(OUT,0,"converged",tmp); 00617 //chisq 00618 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00619 tmp_= mxGetPr(tmp); 00620 tmp_[0]= globs.chisq; 00621 mxSetField(OUT,0,"chisq",tmp); 00622 //error 00623 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00624 tmp_= mxGetPr(tmp); 00625 tmp_[0]= 0; 00626 mxSetField(OUT,0,"except",tmp); 00627 //Lambda 00628 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00629 tmp_= mxGetPr(tmp); 00630 tmp_[0]= globs.Lambda; 00631 mxSetField(OUT,0,"Lambda",tmp); 00632 //cond 00633 tmp = mxCreateDoubleMatrix(1,1,mxREAL); 00634 tmp_= mxGetPr(tmp); 00635 tmp_[0]= globs.cond; 00636 mxSetField(OUT,0,"cond",tmp); 00637 00638 //Experiment Specific output -> OUTEX 00639 00640 dims[0]=1; 00641 dims[1]=globs.nrExp; 00642 00643 OUTEX = mxCreateStructArray(2,dims,sizeof(fieldNamesOUTEX)/sizeof(*fieldNamesOUTEX),fieldNamesOUTEX); 00644 00645 00646 for(nExp=1;nExp<=globs.nrExp;nExp++) 00647 { 00648 //p 00649 tmp = mxCreateDoubleMatrix(1,globs.npar,mxREAL); 00650 tmp_= mxGetPr(tmp); 00651 for(i=1;i<=globs.npar;i++) 00652 tmp_[i-1]= ex[nExp].par[i]; 00653 mxSetField(OUTEX,nExp-1,"p",tmp); 00654 //errP 00655 tmp = mxCreateDoubleMatrix(1,globs.npar,mxREAL); 00656 tmp_= mxGetPr(tmp); 00657 for(i=1;i<=globs.npar;i++) 00658 tmp_[i-1]= ex[nExp].errP[i]; 00659 mxSetField(OUTEX,nExp-1,"errP",tmp); 00660 //y0 00661 tmp = mxCreateDoubleMatrix(1,ex[nExp].nvar,mxREAL); 00662 tmp_= mxGetPr(tmp); 00663 for(i=1;i<=ex[nExp].nvar;i++) 00664 tmp_[i-1]= ex[nExp].yTry[1][i]; 00665 mxSetField(OUTEX,nExp-1,"y0",tmp); 00666 //errY0 00667 tmp = mxCreateDoubleMatrix(1,ex[nExp].nvar,mxREAL); 00668 tmp_= mxGetPr(tmp); 00669 for(i=1;i<=ex[nExp].nvar;i++) 00670 tmp_[i-1]= ex[nExp].errY0[i]; 00671 mxSetField(OUTEX,nExp-1,"errY0",tmp); 00672 //mesh 00673 tmp = mxCreateDoubleMatrix(ex[nExp].nPoints-2,ex[nExp].nvar+1,mxREAL); 00674 tmp_= mxGetPr(tmp); 00675 for(i=1;i<=ex[nExp].nPoints-2;i++) 00676 { 00677 for(j=1;j<=ex[nExp].nvar+1;j++) 00678 { 00679 if(j==1) 00680 tmp_[i-1]=ex[nExp].mesh[i+1]; 00681 else 00682 tmp_[(j-1)*(ex[nExp].nPoints-2)+i-1]=ex[nExp].yTry[i+1][j-1]; 00683 } 00684 } 00685 mxSetField(OUTEX,nExp-1,"mesh",tmp); 00686 } 00687 00688 freeMem(ex,&globs,FALSE); 00689 free_lvector(ndat,1,globs.nrExp); 00690 delete ex; 00691 00692 }
|
|
Definition at line 11 of file outFit.cc. References Glob::covar, dmatrix(), FALSE, Glob::npar, Glob::nrExp, GlobExp::nvar, Glob::strategy, and Glob::y0fix. 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 }
|
|
Definition at line 11 of file setMesh.cc. References dvector(), FALSE, GlobExp::fitend, GlobExp::fitstart, GlobExp::mesh, Glob::noMeasurements, and GlobExp::nPoints. 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. References dmatrix(), dvector(), GlobExp::nPoints, GlobExp::nvar, GlobExp::y0, and GlobExp::yTry. 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++)
Definition at line 31 of file diffit_mex.cc. |