- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/diffit_mex.cc

Go to the documentation of this file.

    00001 #include "mex.h"
    00002 #include<iostream>
    00003 #include<fstream>
    00004 #include<stdlib.h>
    00005 #include<string>
    00006 #include<math.h>
    00007
    00008
    00009 #include "def.h"
    00010 #include "model.h"
    00011 #include "nr.h"
    00012
    00013 using namespace std;
    00014
    00015
    00016 //Begin definition of module prototypes
    00017
    00018 void setMesh(GlobExp *ex,Glob *globs,long expNr);
    00019 void outFit(GlobExp ex[],Glob *globs);
    00020 void freeMem(GlobExp *ex,Glob *globs,int simit);
    00021 void initialise(GlobExp ex[],Glob *globs,int simit);
    00022 void simInit(GlobExp *ex,Glob *globs);
    00023
    00024
    00025 //End definition of module prototypes
    00026
    00027 //The numerics subprogram
    00028 void fitIt(GlobExp ex[],Glob *globs);
    00029
    00030 // DEBUG stream
    00031 ofstream *dbg;
    00032
    00033 //input
    00034 #define IN prhs[0]
    00035 #define INEX prhs[1]
    00036 //output
    00037 #define OUT plhs[0]
    00038 #define OUTEX plhs[1]
    00039
    00040
    00041 void mexFunction( int nlhs,
    00042                   mxArray *plhs[],
    00043                   int nrhs,
    00044                   const mxArray *prhs[])
    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 }