- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/simit_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<math.h>
    00006
    00007 #include "def.h"
    00008 #include "model.h"
    00009 #include "nr.h"
    00010
    00011 using namespace std;
    00012
    00013 //Begin definition of module prototypes
    00014
    00015 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr);
    00016
    00017 void initialise(GlobExp ex[],Glob *globs,int simit);
    00018
    00019 void tabulateValues (Glob *globs,GlobExp *ex,double t0, double t1, double *t, long n, double *state,
    00020                      double **y, double ***dmds, double ***dmdp, double **dyds,double **dydp);
    00021
    00022 void outSimit(GlobExp ex[],Glob *globs,double *t,long n,double **y);
    00023 void freeMem(GlobExp *ex,Glob *globs,int simit);
    00024
    00025 //End definition of module prototypes
    00026
    00027
    00028 // DEBUG stream
    00029 ofstream *dbg;
    00030 // Input
    00031 #define IN   prhs[0]
    00032
    00033 // Output
    00034
    00035 #define N    plhs[0]
    00036 #define YOUT plhs[1]
    00037 #define DMDS plhs[2]
    00038 #define DMDP plhs[3]
    00039
    00040
    00041 void mexFunction( int nlhs,
    00042                   mxArray *plhs[],
    00043                   int nrhs,
    00044                   const mxArray *prhs[]
    00045                   )
    00046 {
    00047   long k,i,j,n,ldum;
    00048   char outstr[100];
    00049   Glob globs;
    00050   GlobExp *ex;
    00051   double time,**y,*state;
    00052   double t0; //init. time
    00053   double *t,*y0,*p;
    00054   double **yout,*yout_,**ydum;
    00055   double ***dmds,***dmdp,*dmds_,*dmdp_;
    00056   double **dyds,**dydp;
    00057   double *n_,*spline_;
    00058   double *tdum=dvector(1,1);
    00059   mxArray *splinedat;
    00060   mxArray *tmp;
    00061   mxArray *SPLINE;
    00062   double *tmp_;
    00063   int hidden;
    00064
    00065   // input -> initial time
    00066   tmp=mxGetField(IN,0,"t");
    00067   tmp_=mxGetPr(tmp);
    00068   n=mxGetNumberOfElements(tmp);
    00069   t=dvector(1,n);
    00070   for(i=1;i<=n;i++)
    00071     t[i]=tmp_[i-1];
    00072
    00073   //allocation of init. state and parameters
    00074   y0=dvector(1,NEQNS);
    00075   p=dvector(1,NPARAMS);
    00076
    00077   // input -> initial time
    00078   tmp=mxGetField(IN,0,"t0");
    00079   tmp_=mxGetPr(tmp);
    00080   t0=tmp_[0];
    00081
    00082   // input -> y0
    00083   tmp=mxGetField(IN,0,"y0");
    00084   tmp_=mxGetPr(tmp);
    00085   if(mxGetNumberOfElements(tmp)!=NEQNS)
    00086     {
    00087       if(mxGetNumberOfElements(tmp)!=0)
    00088         mexWarnMsgTxt("Insufficient number of initial values.\nUse predefined.\n");
    00089       for(i=1;i<=NEQNS;i++)
    00090         y0[i]=DefYValues[i-1];
    00091     }
    00092   else
    00093     {
    00094       for(i=1;i<=NEQNS;i++)
    00095         y0[i]=tmp_[i-1];
    00096     }
    00097
    00098   // input -> p
    00099   tmp=mxGetField(IN,0,"p");
    00100   tmp_=mxGetPr(tmp);
    00101   if(mxGetNumberOfElements(tmp)!=NPARAMS)
    00102     {
    00103       if(mxGetNumberOfElements(tmp)!=0)
    00104         mexWarnMsgTxt("Insufficient number of parameters.\nUse predefined.\n");
    00105       for(i=1;i<=NPARAMS;i++)
    00106         p[i]=DefParameters[i-1];
    00107     }
    00108   else
    00109     {
    00110       for(i=1;i<=NPARAMS;i++)
    00111         p[i]=tmp_[i-1];
    00112     }
    00113
    00114   // input -> eps
    00115   tmp=mxGetField(IN,0,"eps");
    00116   tmp_=mxGetPr(tmp);
    00117   globs.eps=fabs(tmp_[0]);
    00118
    00119   // input -> integrator
    00120   tmp=mxGetField(IN,0,"int");
    00121   tmp_=mxGetPr(tmp);
    00122   globs.integrator=abs((int)tmp_[0]);
    00123   globs.stiff=TRUE;
    00124   // input -> spline
    00125   SPLINE=mxGetField(IN,0,"spline");
    00126
    00127   // input -> hidden
    00128   tmp=mxGetField(IN,0,"hidden");
    00129   tmp_=mxGetPr(tmp);
    00130   hidden=(int)tmp_[0];
    00131
    00132   //allocation of trajectory vect. for output
    00133   if(hidden==TRUE)
    00134     YOUT=mxCreateDoubleMatrix(n,NOBS+NEQNS,mxREAL);
    00135   else
    00136     YOUT=mxCreateDoubleMatrix(n,NOBS,mxREAL);
    00137   yout_=mxGetPr(YOUT);
    00138   yout=dmatrix(1,n,1,NOBS+NEQNS);
    00139   ydum=dmatrix(1,2,1,NOBS+NEQNS);
    00140
    00141   //initialise some global parameters in globs
    00142   globs.noGnu=TRUE;
    00143   globs.npar=NPARAMS;
    00144   globs.noMeasurements=FALSE;
    00145   globs.doP=ivector(1,globs.npar);
    00146   for(k=1;k<=globs.npar;k++)
    00147     globs.doP[k]=TRUE;
    00148   globs.maxit=1000;
    00149   globs.gnuFp=NULL;
    00150   globs.wait=FALSE;
    00151   globs.usesig=FALSE;
    00152   globs.stiff=TRUE;
    00153   globs.maxstp=5000;
    00154   globs.minimp=0.05;
    00155   globs.nowait=FALSE;
    00156   globs.elastic=1.;
    00157   globs.reg=FALSE;
    00158   globs.epsilon=1e-10;
    00159   globs.lambda=1e6;
    00160   globs.dt=0.1;
    00161   globs.sig=0.;
    00162   globs.nrExp=1;
    00163
    00164   ex=new GlobExp[globs.nrExp+1];
    00165   //initialise some global parameters
    00166
    00167   ex[1].nobs=NOBS;
    00168   ex[1].nvar=NEQNS;
    00169   ex[1].par=dvector(1,NPARAMS);
    00170   ex[1].y0=dvector(1,NEQNS);
    00171
    00172
    00173   for(k=1;k<=NPARAMS;k++)
    00174     ex[1].par[k]=p[k];
    00175
    00176
    00177   ex[1].fitstart=t[1];
    00178   ex[1].fitend=t[n];
    00179   //initial values
    00180
    00181   state=dvector(1,ex[1].nvar);
    00182   for(k=1;k<=NEQNS;k++)
    00183     state[k]=y0[k];
    00184
    00185
    00186   N=mxCreateDoubleMatrix(3,1,mxREAL);
    00187   n_=mxGetPr(N);
    00188
    00189   n_[0]=NEQNS;
    00190   n_[1]=NOBS;
    00191   n_[2]=NPARAMS;
    00192
    00193   // spline section
    00194   if(NSPLINES==0)
    00195     {
    00196       globs.initSpline=TRUE;
    00197       initialise(ex,&globs,TRUE);
    00198     }
    00199   else
    00200     {
    00201       globs.initSpline=FALSE;
    00202       initialise(ex,&globs,TRUE);
    00203       ex[1].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
    00204       ex[1].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
    00205       ex[1].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*));
    00206       ex[1].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*));
    00207
    00208       for(i=1;i<=NSPLINES;i++)
    00209         {
    00210           if(mxGetNumberOfElements(SPLINE) < NSPLINES)
    00211             {
    00212               cerr << "Please specify " << NSPLINES << " spline(s).\n";
    00213               return;
    00214             }
    00215
    00216           splinedat=mxGetCell(SPLINE,i-1);
    00217           spline_=mxGetPr(splinedat);
    00218           ex[1].nNodes[i]=mxGetM(splinedat);
    00219           if(mxGetN(splinedat)!=3)
    00220             {
    00221               cerr << "Invalid data for spline " << i << ".\n";
    00222               return;
    00223             }
    00224           ex[1].splineNodes[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
    00225           ex[1].splineY[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
    00226           ex[1].splineGam[i]=(double *) malloc((size_t)(ex[1].nNodes[i]+1)*sizeof(double));
    00227
    00228           //copy spline data
    00229           for(j=1;j<=ex[1].nNodes[i];j++)
    00230              {
    00231                ex[1].splineNodes[i][j]=spline_[j-1];
    00232                ex[1].splineY[i][j]=spline_[ex[1].nNodes[i]+j-1];
    00233                ex[1].splineGam[i][j]=spline_[2*ex[1].nNodes[i]+j-1];
    00234              }
    00235         }
    00236     }
    00237
    00238   //... and integrate
    00239   try
    00240     {
    00241       if(nlhs!=2)
    00242         {
    00243           dmds=d3tensor(1,n,1,NEQNS,1,NEQNS);
    00244           dmdp=d3tensor(1,n,1,NEQNS,1,NPARAMS);
    00245           dyds=dmatrix(1,NEQNS,1,NEQNS);
    00246           dydp=dmatrix(1,NEQNS,1,NPARAMS);
    00247
    00248           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,dmds,dmdp,dyds,dydp);
    00249
    00250           DMDS=mxCreateDoubleMatrix(NOBS*NEQNS,n,mxREAL);
    00251           dmds_=mxGetPr(DMDS);
    00252           DMDP=mxCreateDoubleMatrix(NOBS*NPARAMS,n,mxREAL);
    00253           dmdp_=mxGetPr(DMDP);
    00254
    00255
    00256           //copying
    00257
    00258           for(k=1;k<=NEQNS;k++)
    00259             for(i=1;i<=NOBS;i++)
    00260               for(j=1;j<=n;j++)
    00261                 {
    00262                   dmds_[(j-1)*NOBS*NEQNS+(k-1)*NOBS+i-1]=dmds[j][i][k];
    00263                 }
    00264
    00265
    00266           for(k=1;k<=NPARAMS;k++)
    00267             for(i=1;i<=NOBS;i++)
    00268               for(j=1;j<=n;j++)
    00269                 {
    00270                   dmdp_[(j-1)*NOBS*NPARAMS+(k-1)*NOBS+i-1]=dmdp[j][i][k];
    00271                 }
    00272
    00273           free_d3tensor(dmds,1,n,1,NEQNS,1,NEQNS);
    00274           free_d3tensor(dmdp,1,n,1,NEQNS,1,NPARAMS);
    00275           free_dmatrix(dyds,1,NEQNS,1,NEQNS);
    00276           free_dmatrix(dydp,1,NEQNS,1,NPARAMS);
    00277         }
    00278       else
    00279         {
    00280           tabulateValues(&globs,&ex[1],t0,t[n],t,n,state,yout,NULL,NULL,NULL,NULL);
    00281
    00282         }
    00283
    00284       if(hidden==TRUE)
    00285         {
    00286           ldum=NOBS+NEQNS;
    00287           for(i=1;i<=NEQNS;i++)
    00288             state[i]=y0[i];
    00289           tdum[1]=t[1];
    00290           if(t0!=t[1])
    00291             tabulateValues(&globs,&ex[1],t0,t[1],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
    00292           for(j=1;j<=NEQNS;j++)
    00293             yout[1][NOBS+j]=state[j];
    00294
    00295           for(i=2;i<=n;i++)
    00296             {
    00297               tdum[1]=t[i];
    00298               tabulateValues(&globs,&ex[1],t[i-1],t[i],tdum,1,state,ydum,NULL,NULL,NULL,NULL);
    00299               for(j=1;j<=NEQNS;j++)
    00300                 yout[i][NOBS+j]=state[j];
    00301
    00302             }
    00303         }
    00304       else
    00305         ldum=NOBS;
    00306
    00307       for(i=1;i<=ldum;i++)
    00308         for(j=1;j<=n;j++)
    00309           yout_[(i-1)*n+j-1]=yout[j][i];
    00310     }
    00311   catch (int i)
    00312     {
    00313       return;
    00314     }
    00315   //free mem
    00316   freeMem(ex,&globs,TRUE);
    00317   free_dvector(state,1,ex[1].nvar);
    00318   free_dmatrix(yout,1,n,1,NEQNS+NOBS);
    00319   free_dvector(t,1,n);
    00320   free_dvector(y0,1,NEQNS);
    00321   free_dvector(p,1,NPARAMS);
    00322   free_dvector(ex[1].par,1,NPARAMS);
    00323   free_dvector(ex[1].y0,1,NEQNS);
    00324   free_dvector(tdum,1,1);
    00325   free_dmatrix(ydum,1,2,1,NOBS+NEQNS);
    00326
    00327   delete ex;
    00328 }