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
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
00026
00027
00028
00029 ofstream *dbg;
00030
00031 #define IN prhs[0]
00032
00033
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;
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
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
00074 y0=dvector(1,NEQNS);
00075 p=dvector(1,NPARAMS);
00076
00077
00078 tmp=mxGetField(IN,0,"t0");
00079 tmp_=mxGetPr(tmp);
00080 t0=tmp_[0];
00081
00082
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
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
00115 tmp=mxGetField(IN,0,"eps");
00116 tmp_=mxGetPr(tmp);
00117 globs.eps=fabs(tmp_[0]);
00118
00119
00120 tmp=mxGetField(IN,0,"int");
00121 tmp_=mxGetPr(tmp);
00122 globs.integrator=abs((int)tmp_[0]);
00123 globs.stiff=TRUE;
00124
00125 SPLINE=mxGetField(IN,0,"spline");
00126
00127
00128 tmp=mxGetField(IN,0,"hidden");
00129 tmp_=mxGetPr(tmp);
00130 hidden=(int)tmp_[0];
00131
00132
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
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
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
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
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
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
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
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
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 }