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

#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++)


Define Documentation

#define IN   prhs[0]
 

Definition at line 34 of file diffit_mex.cc.

Referenced by mexFunction().

#define INEX   prhs[1]
 

Definition at line 35 of file diffit_mex.cc.

#define OUT   plhs[0]
 

Definition at line 37 of file diffit_mex.cc.

#define OUTEX   plhs[1]
 

Definition at line 38 of file diffit_mex.cc.


Function Documentation

void fitIt GlobExp  ex[],
Glob globs
 

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 }

void freeMem GlobExp ex,
Glob globs,
int  simit
 

Definition at line 11 of file freeMem.cc.

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

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

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 }

void mexFunction int  nlhs,
mxArray *  plhs[],
int  nrhs,
const mxArray *  prhs[]
 

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 }

void outFit GlobExp  ex[],
Glob globs
 

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 }

void setMesh GlobExp ex,
Glob globs,
long  expNr
 

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 }

void simInit GlobExp ex,
Glob globs
 

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 }


Variable Documentation

ofstream* dbg
 

debug stream (C++)

Definition at line 31 of file diffit_mex.cc.


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