/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 }

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