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
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
00026
00027
00028 void fitIt(GlobExp ex[],Glob *globs);
00029
00030
00031 ofstream *dbg;
00032
00033
00034 #define IN prhs[0]
00035 #define INEX prhs[1]
00036
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
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
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
00079 tmp=mxGetField(IN,0,"eps");
00080 tmp_=mxGetPr(tmp);
00081 globs.eps=fabs(tmp_[0]);
00082
00083 tmp=mxGetField(IN,0,"nognu");
00084 tmp_=mxGetPr(tmp);
00085 globs.noGnu=(int)tmp_[0];
00086
00087 tmp=mxGetField(IN,0,"nomeasure");
00088 tmp_=mxGetPr(tmp);
00089 globs.noMeasurements=(int)tmp_[0];
00090
00091 tmp=mxGetField(IN,0,"maxit");
00092 tmp_=mxGetPr(tmp);
00093 globs.maxit=abs((int)tmp_[0]);
00094
00095 tmp=mxGetField(IN,0,"int");
00096 tmp_=mxGetPr(tmp);
00097 globs.integrator=abs((int)tmp_[0]);
00098 globs.stiff=TRUE;
00099
00100 tmp=mxGetField(IN,0,"minimp");
00101 tmp_=mxGetPr(tmp);
00102 globs.minimp=fabs(tmp_[0]);
00103
00104 tmp=mxGetField(IN,0,"reg");
00105 tmp_=mxGetPr(tmp);
00106 globs.reg=(int)tmp_[0];
00107
00108 tmp=mxGetField(IN,0,"epsilon");
00109 tmp_=mxGetPr(tmp);
00110 globs.epsilon=fabs(tmp_[0]);
00111
00112 tmp=mxGetField(IN,0,"lambda");
00113 tmp_=mxGetPr(tmp);
00114 globs.lambda=fabs(tmp_[0]);
00115
00116 tmp=mxGetField(IN,0,"siminit");
00117 tmp_=mxGetPr(tmp);
00118 globs.simInit=(int)tmp_[0];
00119
00120 tmp=mxGetField(IN,0,"pert");
00121 tmp_=mxGetPr(tmp);
00122 globs.pert=fabs(tmp_[0]);
00123
00124 tmp=mxGetField(IN,0,"nodamp");
00125 tmp_=mxGetPr(tmp);
00126 globs.nodamp=(int)tmp_[0];
00127
00128 tmp=mxGetField(IN,0,"nexp");
00129 tmp_=mxGetPr(tmp);
00130 globs.nrExp=(long)tmp_[0];
00131
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
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
00190 tmp=mxGetField(IN,0,"wquer");
00191 tmp_=mxGetPr(tmp);
00192 globs.wquer=tmp_[0];
00193
00194 tmp=mxGetField(IN,0,"init");
00195 tmp_=mxGetPr(tmp);
00196 init=(int)tmp_[0];
00197
00198 tmp=mxGetField(IN,0,"silent");
00199 tmp_=mxGetPr(tmp);
00200 globs.silent=(int)tmp_[0];
00201
00202 tmp=mxGetField(IN,0,"opt");
00203 tmp_=mxGetPr(tmp);
00204 globs.minimiser=(int)tmp_[0];
00205
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
00227 ndat=lvector(1,globs.nrExp);
00228 ex=new GlobExp[globs.nrExp+1];
00229
00230 for(nExp=1;nExp<=globs.nrExp;nExp++)
00231 {
00232
00233 ex[nExp].nvar=NEQNS;
00234 ex[nExp].splinesDefined=FALSE;
00235
00236 tmp=mxGetField(INEX,nExp-1,"t0");
00237 tmp_=mxGetPr(tmp);
00238 ex[nExp].fitstart=tmp_[0];
00239
00240
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
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
00259 tmp=mxGetField(INEX,nExp-1,"nms");
00260 tmp_=mxGetPr(tmp);
00261 ex[nExp].nPoints=(long)tmp_[0]+1;
00262
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
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
00306 tmp=mxGetField(INEX,nExp-1,"n");
00307 tmp_=mxGetPr(tmp);
00308 ndat[nExp]=(long)tmp_[0];
00309
00310 }
00311
00312
00313 dbg=new ofstream("diffit.dbg");
00314 if (!dbg)
00315 {
00316 cerr << "Error opening DEBUG file.\n";
00317 return;
00318 }
00319 dbg->precision(4);
00320 dbg->unsetf(ios::floatfield);
00321
00322
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
00344
00345
00346 tmp=mxGetField(INEX,i-1,"data");
00347 tmp_=mxGetPr(tmp);
00348 sig=mxGetField(INEX,i-1,"sig");
00349 sig_=mxGetPr(sig);
00350
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
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
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
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
00393 ex[i].firstMeasure=1;
00394 ex[i].lastMeasure=ex[i].nMeasure;
00395
00396
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
00416 for(j=1;j<=ex[i].nPoints-2;j++)
00417 {
00418 ex[i].mesh[j+1]=tmp_[j-1];
00419 }
00420
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
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 }
00471
00472
00473
00474
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
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
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
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)
00565 {
00566 dims[0]=1;
00567 dims[1]=1;
00568 OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT);
00569
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
00592 if(init!=1)
00593 {
00594 outFit(ex,&globs);
00595 if(!globs.noGnu)
00596 system("rm gnuout.dat");
00597 }
00598
00599
00600
00601 dims[0]=1;
00602 dims[1]=1;
00603
00604
00605 OUT = mxCreateStructArray(2,dims,sizeof(fieldNamesOUT)/sizeof(*fieldNamesOUT),fieldNamesOUT);
00606
00607
00608 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00609 tmp_= mxGetPr(tmp);
00610 tmp_[0]= globs.wquer;
00611 mxSetField(OUT,0,"wquer",tmp);
00612
00613 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00614 tmp_= mxGetPr(tmp);
00615 tmp_[0]= globs.fitConverged;
00616 mxSetField(OUT,0,"converged",tmp);
00617
00618 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00619 tmp_= mxGetPr(tmp);
00620 tmp_[0]= globs.chisq;
00621 mxSetField(OUT,0,"chisq",tmp);
00622
00623 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00624 tmp_= mxGetPr(tmp);
00625 tmp_[0]= 0;
00626 mxSetField(OUT,0,"except",tmp);
00627
00628 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00629 tmp_= mxGetPr(tmp);
00630 tmp_[0]= globs.Lambda;
00631 mxSetField(OUT,0,"Lambda",tmp);
00632
00633 tmp = mxCreateDoubleMatrix(1,1,mxREAL);
00634 tmp_= mxGetPr(tmp);
00635 tmp_[0]= globs.cond;
00636 mxSetField(OUT,0,"cond",tmp);
00637
00638
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
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
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
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
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
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 }