- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /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.