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

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