[Pomp-commits] r686 - in pkg/pomp: . src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 28 20:01:23 CEST 2012
Author: kingaa
Date: 2012-04-28 20:01:22 +0200 (Sat, 28 Apr 2012)
New Revision: 686
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/src/euler.c
pkg/pomp/src/rmeasure.c
Log:
- more changes to the guts
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-27 21:51:50 UTC (rev 685)
+++ pkg/pomp/DESCRIPTION 2012-04-28 18:01:22 UTC (rev 686)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.42-1
-Date: 2012-04-27
+Date: 2012-04-28
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/pomp/src/euler.c
===================================================================
--- pkg/pomp/src/euler.c 2012-04-27 21:51:50 UTC (rev 685)
+++ pkg/pomp/src/euler.c 2012-04-28 18:01:22 UTC (rev 686)
@@ -10,28 +10,23 @@
SEXP tcovar, SEXP covar, SEXP args)
{
int nprotect = 0;
- int first = 1;
int mode = -1;
- int meth = 0;
- int use_names;
- int *dim, *posn, xdim[3];
int nstep, nvars, npars, nreps, ntimes, nzeros, ncovars, covlen;
- pomp_onestep_sim *ff = NULL;
SEXP X;
- SEXP fn, fcall, ans, rho, nm;
+ SEXP fn, fcall, rho, ans, nm;
SEXP Snames, Pnames, Cnames;
SEXP tvec, xvec, pvec, cvec, dtvec;
- int *sidx, *pidx, *cidx, *zidx;
- double *tp, *xp, *pp, *cp, *dtp, *xt, *ps;
- double *time;
- double *Xt;
- double t, dt;
- int i, j, k, step;
+ int *pidx = 0, *sidx = 0, *cidx = 0, *zidx = 0;
+ pomp_onestep_sim *ff = NULL;
+ int meth = *(INTEGER(AS_INTEGER(method))); // 0 = Euler, 1 = one-step, 2 = fixed step
- dim = INTEGER(GET_DIM(xstart)); nvars = dim[0]; nreps = dim[1];
- dim = INTEGER(GET_DIM(params)); npars = dim[0];
- dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; ncovars = dim[1];
- ntimes = LENGTH(times);
+ {
+ int *dim;
+ dim = INTEGER(GET_DIM(xstart)); nvars = dim[0]; nreps = dim[1];
+ dim = INTEGER(GET_DIM(params)); npars = dim[0];
+ dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; ncovars = dim[1];
+ ntimes = LENGTH(times);
+ }
PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(xstart))); nprotect++;
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
@@ -43,7 +38,6 @@
// vector for interpolated covariates
PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
SET_NAMES(cvec,Cnames);
- cp = REAL(cvec);
// indices of accumulator variables
nzeros = LENGTH(zeronames);
@@ -55,32 +49,14 @@
// set up
switch (mode) {
- case 1: // native code
-
- // construct state, parameter, covariate, observable indices
- sidx = INTEGER(PROTECT(matchnames(Snames,statenames))); nprotect++;
- pidx = INTEGER(PROTECT(matchnames(Pnames,paramnames))); nprotect++;
- cidx = INTEGER(PROTECT(matchnames(Cnames,covarnames))); nprotect++;
-
- ff = (pomp_onestep_sim *) R_ExternalPtrAddr(fn);
-
- break;
-
case 0: // R function
- // get function's environment
- PROTECT(rho = (CLOENV(fn))); nprotect++;
-
PROTECT(dtvec = NEW_NUMERIC(1)); nprotect++;
- dtp = REAL(dtvec);
PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
- tp = REAL(tvec);
PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
- SET_NAMES(xvec,Snames);
- xp = REAL(xvec);
PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
+ SET_NAMES(xvec,Snames);
SET_NAMES(pvec,Pnames);
- pp = REAL(pvec);
// set up the function call
PROTECT(fcall = LCONS(cvec,args)); nprotect++;
@@ -95,68 +71,20 @@
SET_TAG(fcall,install("x"));
PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
- // to hold indices of variables that must be rearranged
- posn = (int *) R_alloc(nvars,sizeof(int));
+ // get function's environment
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
- // this is the euler step function that is evaluated when the user supplies an R function
- // (and not a native routine)
- // Note that stateindex, parindex, covindex are ignored.
- void R_step_fn (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int ncovar, const double *covar,
- double t, double dt)
- {
- int nprotect = 0;
- int *op, i;
- double *xs;
- SEXP ans, nm;
+ break;
- for (i = 0; i < nvars; i++) xp[i] = x[i];
- for (i = 0; i < npars; i++) pp[i] = p[i];
- *tp = t;
- *dtp = dt;
-
- if (first) {
+ case 1: // native code
- PROTECT(ans = eval(fcall,rho)); nprotect++; // evaluate the call
-
- if (LENGTH(ans) != nvars) {
- error("user 'step.fun' returns a vector of %d states but %d are expected: compare initial conditions?",
- LENGTH(ans),nvars);
- }
+ // construct state, parameter, covariate indices
+ sidx = INTEGER(PROTECT(matchnames(Snames,statenames))); nprotect++;
+ pidx = INTEGER(PROTECT(matchnames(Pnames,paramnames))); nprotect++;
+ cidx = INTEGER(PROTECT(matchnames(Cnames,covarnames))); nprotect++;
- PROTECT(nm = GET_NAMES(ans)); nprotect++;
- use_names = !isNull(nm);
- if (use_names) {
- op = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
- for (i = 0; i < nvars; i++) posn[i] = op[i];
- }
+ ff = (pomp_onestep_sim *) R_ExternalPtrAddr(fn);
- xs = REAL(AS_NUMERIC(ans));
-
- first = 0;
-
- } else {
-
- xs = REAL(AS_NUMERIC(eval(fcall,rho)));
-
- }
-
- if (use_names) {
- for (i = 0; i < nvars; i++) x[posn[i]] = xs[i];
- } else {
- for (i = 0; i < nvars; i++) x[i] = xs[i];
- }
-
- UNPROTECT(nprotect);
- }
-
- sidx = 0;
- pidx = 0;
- cidx = 0;
-
- ff = (pomp_onestep_sim *) R_step_fn;
-
break;
default:
@@ -164,83 +92,145 @@
break;
}
+ // create array to hold results
+ {
+ int dim[3] = {nvars, nreps, ntimes};
+ PROTECT(X = makearray(3,dim)); nprotect++;
+ setrownames(X,Snames,3);
+ }
+
+ // copy the start values into the result array
+ memcpy(REAL(X),REAL(xstart),nvars*nreps*sizeof(double));
+
if (mode==1) {
set_pomp_userdata(args);
GetRNGstate();
}
- // create array to hold results
- xdim[0] = nvars; xdim[1] = nreps; xdim[2] = ntimes;
- PROTECT(X = makearray(3,xdim)); nprotect++;
- setrownames(X,Snames,3);
- Xt = REAL(X);
+ // now do computations
+ {
+ int first = 1;
+ int use_names = 0;
+ int *posn;
+ double *time = REAL(times);
+ double *xs = REAL(X);
+ double *xt = REAL(X)+nvars*nreps;
+ double *cp = REAL(cvec);
+ double *ps = REAL(params);
+ double t = time[0];
+ double dt;
+ double *pm, *xm;
+ int i, j, k, step;
- // copy the start values into the result array
- xt = REAL(xstart);
- for (j = 0; j < nreps; j++)
- for (i = 0; i < nvars; i++)
- Xt[i+nvars*j] = xt[i+nvars*j];
+ for (step = 1; step < ntimes; step++, xs = xt, xt += nvars*nreps) {
- meth = *(INTEGER(AS_INTEGER(method))); // 0 = Euler, 1 = one-step, 2 = fixed step
+ R_CheckUserInterrupt();
+
+ if (t > time[step]) {
+ error("'times' is not an increasing sequence");
+ }
- // now do computations
- // loop over times
- time = REAL(times);
- t = time[0];
+ memcpy(xt,xs,nreps*nvars*sizeof(double));
+
+ // set accumulator variables to zero
+ for (j = 0; j < nreps; j++)
+ for (i = 0; i < nzeros; i++)
+ xt[zidx[i]+nvars*j] = 0.0;
- for (step = 1; step < ntimes; step++) {
+ switch (meth) {
+ case 0: // Euler method
+ dt = *(REAL(deltat));
+ nstep = num_euler_steps(t,time[step],&dt);
+ break;
+ case 1: // one step
+ dt = time[step]-t;
+ nstep = (dt > 0) ? 1 : 0;
+ break;
+ case 2: // fixed step
+ dt = *(REAL(deltat));
+ nstep = num_map_steps(t,time[step],dt);
+ break;
+ default:
+ error("unrecognized 'method' in 'euler_model_simulator'");
+ break;
+ }
- R_CheckUserInterrupt();
+ for (k = 0; k < nstep; k++) { // loop over Euler steps
- if (t > time[step]) {
- error("'times' is not an increasing sequence");
- }
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,t,cp,0);
- switch (meth) {
- case 0: // Euler method
- dt = *(REAL(deltat));
- nstep = num_euler_steps(t,time[step],&dt);
- break;
- case 1: // one step
- dt = time[step]-t;
- nstep = (dt > 0) ? 1 : 0;
- break;
- case 2: // fixed step
- dt = *(REAL(deltat));
- nstep = num_map_steps(t,time[step],dt);
- break;
- default:
- error("unrecognized 'method' in 'stepwise_simulator'");
- }
+ for (j = 0, pm = ps, xm = xt; j < nreps; j++, pm += npars, xm += nvars) { // loop over replicates
+
+ switch (mode) {
- for (j = 0; j < nreps; j++) {
- xt = &Xt[nvars*(j+nreps*step)];
- // copy in the previous values of the state variables
- for (i = 0; i < nvars; i++) xt[i] = Xt[i+nvars*(j+nreps*(step-1))];
- // set some variables to zero
- for (i = 0; i < nzeros; i++) xt[zidx[i]] = 0.0;
- }
+ case 0: // R function
- for (k = 0; k < nstep; k++) { // loop over Euler steps
+ {
+ double *xp = REAL(xvec);
+ double *pp = REAL(pvec);
+ double *tp = REAL(tvec);
+ double *dtp = REAL(dtvec);
+ double *ap;
+
+ *tp = t;
+ *dtp = dt;
+ memcpy(xp,xm,nvars*sizeof(double));
+ memcpy(pp,pm,npars*sizeof(double));
+
+ if (first) {
- // interpolate the covar functions for the covariates
- table_lookup(&covariate_table,t,cp,0);
+ PROTECT(ans = eval(fcall,rho)); nprotect++; // evaluate the call
+ if (LENGTH(ans) != nvars) {
+ error("user 'step.fun' returns a vector of %d states but %d are expected: compare initial conditions?",
+ LENGTH(ans),nvars);
+ }
+
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) {
+ posn = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
+ }
- for (j = 0, ps = REAL(params); j < nreps; j++, ps += npars) { // loop over replicates
-
- xt = &Xt[nvars*(j+nreps*step)];
+ ap = REAL(AS_NUMERIC(ans));
+
+ first = 0;
- (*ff)(xt,ps,sidx,pidx,cidx,ncovars,cp,t,dt);
+ } else {
+
+ ap = REAL(AS_NUMERIC(eval(fcall,rho)));
- }
+ }
+
+ if (use_names) {
+ for (i = 0; i < nvars; i++) xm[posn[i]] = ap[i];
+ } else {
+ for (i = 0; i < nvars; i++) xm[i] = ap[i];
+ }
- t += dt;
+ }
- if ((method == 0) && (k == nstep-2)) { // penultimate step
- dt = time[step]-t;
- t = time[step]-dt;
- }
+ break;
+
+ case 1: // native code
+ (*ff)(xm,pm,sidx,pidx,cidx,ncovars,cp,t,dt);
+ break;
+
+ default:
+ error("unrecognized 'mode' in 'euler_simulator'");
+ break;
+ }
+
+ }
+
+ t += dt;
+
+ if ((method == 0) && (k == nstep-2)) { // penultimate step
+ dt = time[step]-t;
+ t = time[step]-dt;
+ }
+ }
}
}
@@ -248,225 +238,194 @@
PutRNGstate();
unset_pomp_userdata();
}
-
+
UNPROTECT(nprotect);
return X;
}
// compute pdf of a sequence of Euler steps
-static void euler_densities (pomp_onestep_pdf *estep,
- double *f,
- double *x, double *times, double *params,
- int *ndim,
- int *stateindex, int *parindex, int *covindex,
- double *time_table, double *covar_table,
- int *give_log)
+SEXP euler_model_density (SEXP func,
+ SEXP x, SEXP times, SEXP params,
+ SEXP statenames, SEXP paramnames, SEXP covarnames,
+ SEXP tcovar, SEXP covar, SEXP log, SEXP args)
{
- double *x1p, *x2p, *pp, *fp, t1, t2;
- int nvar = ndim[0];
- int npar = ndim[1];
- int nrep = ndim[2];
- int ntimes = ndim[3];
- int covlen = ndim[4];
- int covdim = ndim[5];
- double covar_fn[covdim];
- int p, step;
+ int nprotect = 0;
+ int mode;
+ int give_log;
+ int nvars, npars, nreps, ntimes, ncovars, covlen;
+ pomp_onestep_pdf *ff = NULL;
+ SEXP t1vec, t2vec, x1vec, x2vec, pvec, cvec;
+ SEXP Snames, Pnames, Cnames;
+ SEXP ans, rho, fcall, fn;
+ SEXP F;
+ int *pidx = 0, *sidx = 0, *cidx = 0;
- // set up the covariate table
- struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
-
- for (step = 0; step < ntimes-1; step++) { // loop over times
+ {
+ int *dim;
+ dim = INTEGER(GET_DIM(x)); nvars = dim[0]; nreps = dim[1];
+ dim = INTEGER(GET_DIM(params)); npars = dim[0];
+ dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; ncovars = dim[1];
+ ntimes = LENGTH(times);
+ }
- R_CheckUserInterrupt(); // check for user interrupt
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
- t1 = times[step];
- t2 = times[step+1];
+ // set up the covariate table
+ struct lookup_table covariate_table = {covlen, ncovars, 0, REAL(tcovar), REAL(covar)};
- // interpolate the covariates at time t1
- if (covdim > 0)
- table_lookup(&covariate_table,t1,covar_fn,0);
-
- for (p = 0; p < nrep; p++) { // loop over replicates
-
- fp = &f[p+nrep*step];
- x1p = &x[nvar*(p+nrep*step)];
- x2p = &x[nvar*(p+nrep*(step+1))];
- pp = ¶ms[npar*p];
-
- (*estep)(fp,x1p,x2p,t1,t2,pp,
- stateindex,parindex,covindex,
- covdim,covar_fn);
-
- if (!(*give_log)) *fp = exp(*fp);
-
- }
- }
-}
+ // vector for interpolated covariates
+ PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
+ SET_NAMES(cvec,Cnames);
-// these global objects will pass the needed information to the user-defined function (see 'default_onestep_dens_fn')
-// each of these is allocated once, globally, and refilled many times
-static SEXP euler_dens_Xvec1; // state variable vector
-static SEXP euler_dens_Xvec2; // state variable vector
-static SEXP euler_dens_Pvec; // parameter vector
-static SEXP euler_dens_Cvec; // covariate vector
-static SEXP euler_dens_time1; // time1
-static SEXP euler_dens_time2; // time2
-static int euler_dens_nvar; // number of state variables
-static int euler_dens_npar; // number of parameters
-static SEXP euler_dens_envir; // function's environment
-static SEXP euler_dens_fcall; // function call
+ PROTECT(fn = pomp_fun_handler(func,&mode)); nprotect++;
-#define X1VEC (euler_dens_Xvec1)
-#define X2VEC (euler_dens_Xvec2)
-#define PVEC (euler_dens_Pvec)
-#define CVEC (euler_dens_Cvec)
-#define TIME1 (euler_dens_time1)
-#define TIME2 (euler_dens_time2)
-#define NVAR (euler_dens_nvar)
-#define NPAR (euler_dens_npar)
-#define RHO (euler_dens_envir)
-#define FCALL (euler_dens_fcall)
+ give_log = *(INTEGER(log));
-// this is the euler dens function that is evaluated when the user supplies an R function
-// (and not a native routine)
-// Note that stateindex, parindex, covindex are ignored.
-static void default_onestep_dens_fn (double *f, double *x1, double *x2, double t1, double t2, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int ncovar, const double *covar)
-{
- int nprotect = 0;
- int k;
- double *xp;
- SEXP ans;
- xp = REAL(X1VEC);
- for (k = 0; k < NVAR; k++) xp[k] = x1[k];
- xp = REAL(X2VEC);
- for (k = 0; k < NVAR; k++) xp[k] = x2[k];
- xp = REAL(PVEC);
- for (k = 0; k < NPAR; k++) xp[k] = p[k];
- xp = REAL(CVEC);
- for (k = 0; k < ncovar; k++) xp[k] = covar[k];
- xp = REAL(TIME1); xp[0] = t1;
- xp = REAL(TIME2); xp[0] = t2;
+ switch (mode) {
- PROTECT(ans = eval(FCALL,RHO)); nprotect++; // evaluate the call
+ case 0: // R function
- xp = REAL(AS_NUMERIC(ans));
- f[0] = xp[0];
- UNPROTECT(nprotect);
-}
+ PROTECT(t1vec = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(t2vec = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(x1vec = NEW_NUMERIC(nvars)); nprotect++;
+ SET_NAMES(x1vec,Snames);
+ PROTECT(x2vec = NEW_NUMERIC(nvars)); nprotect++;
+ SET_NAMES(x2vec,Snames);
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
+ SET_NAMES(pvec,Pnames);
+ // set up the function call
+ PROTECT(fcall = LCONS(cvec,args)); nprotect++;
+ SET_TAG(fcall,install("covars"));
+ PROTECT(fcall = LCONS(pvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("params"));
+ PROTECT(fcall = LCONS(t2vec,fcall)); nprotect++;
+ SET_TAG(fcall,install("t2"));
+ PROTECT(fcall = LCONS(t1vec,fcall)); nprotect++;
+ SET_TAG(fcall,install("t1"));
+ PROTECT(fcall = LCONS(x2vec,fcall)); nprotect++;
+ SET_TAG(fcall,install("x2"));
+ PROTECT(fcall = LCONS(x1vec,fcall)); nprotect++;
+ SET_TAG(fcall,install("x1"));
+ PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
-SEXP euler_model_density (SEXP func,
- SEXP x, SEXP times, SEXP params,
- SEXP statenames, SEXP paramnames, SEXP covarnames,
- SEXP tcovar, SEXP covar, SEXP log, SEXP args)
-{
- int nprotect = 0;
- int *dim, fdim[2], ndim[6];
- int nvar, npar, nrep, ntimes;
- int covlen, covdim;
- int nstates = LENGTH(statenames);
- int nparams = LENGTH(paramnames);
- int ncovars = LENGTH(covarnames);
- pomp_onestep_pdf *ff = NULL;
- SEXP F, pindex, sindex, cindex;
- int *pidx, *sidx, *cidx;
- SEXP fn, Xnames, Pnames, Cnames;
- int use_native;
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
- dim = INTEGER(GET_DIM(x)); nvar = dim[0]; nrep = dim[1];
- dim = INTEGER(GET_DIM(params)); npar = dim[0];
- dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; covdim = dim[1];
- ntimes = LENGTH(times);
+ break;
- PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
- PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
+ case 1: // native code
- PROTECT(fn = pomp_fun_handler(func,&use_native)); nprotect++;
+ // construct state, parameter, covariate indices
+ sidx = INTEGER(PROTECT(matchnames(Snames,statenames))); nprotect++;
+ pidx = INTEGER(PROTECT(matchnames(Pnames,paramnames))); nprotect++;
+ cidx = INTEGER(PROTECT(matchnames(Cnames,covarnames))); nprotect++;
- if (use_native) {
ff = (pomp_onestep_pdf *) R_ExternalPtrAddr(fn);
- } else {
- PROTECT(RHO = (CLOENV(fn))); nprotect++;
- NVAR = nvar; // for internal use
- NPAR = npar; // for internal use
- PROTECT(TIME1 = NEW_NUMERIC(1)); nprotect++; // for internal use
- PROTECT(TIME2 = NEW_NUMERIC(1)); nprotect++; // for internal use
- PROTECT(X1VEC = NEW_NUMERIC(nvar)); nprotect++; // for internal use
- PROTECT(X2VEC = NEW_NUMERIC(nvar)); nprotect++; // for internal use
- PROTECT(PVEC = NEW_NUMERIC(npar)); nprotect++; // for internal use
- PROTECT(CVEC = NEW_NUMERIC(covdim)); nprotect++; // for internal use
- SET_NAMES(X1VEC,Xnames); // make sure the names attribute is copied
- SET_NAMES(X2VEC,Xnames); // make sure the names attribute is copied
- SET_NAMES(PVEC,Pnames); // make sure the names attribute is copied
- SET_NAMES(CVEC,Cnames); // make sure the names attribute is copied
- // set up the function call
- PROTECT(FCALL = LCONS(CVEC,args)); nprotect++;
- SET_TAG(FCALL,install("covars"));
- PROTECT(FCALL = LCONS(PVEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("params"));
- PROTECT(FCALL = LCONS(TIME2,FCALL)); nprotect++;
- SET_TAG(FCALL,install("t2"));
- PROTECT(FCALL = LCONS(TIME1,FCALL)); nprotect++;
- SET_TAG(FCALL,install("t1"));
- PROTECT(FCALL = LCONS(X2VEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("x2"));
- PROTECT(FCALL = LCONS(X1VEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("x1"));
- PROTECT(FCALL = LCONS(fn,FCALL)); nprotect++;
- ff = (pomp_onestep_pdf *) default_onestep_dens_fn;
- }
- fdim[0] = nrep; fdim[1] = ntimes-1;
- PROTECT(F = makearray(2,fdim)); nprotect++;
+ break;
- if (nstates>0) {
- PROTECT(sindex = MATCHROWNAMES(x,statenames)); nprotect++;
- sidx = INTEGER(sindex);
- } else {
- sidx = 0;
+ default:
+ error("unrecognized 'mode' in 'euler_model_density'");
+ break;
}
- if (nparams>0) {
- PROTECT(pindex = MATCHROWNAMES(params,paramnames)); nprotect++;
- pidx = INTEGER(pindex);
- } else {
- pidx = 0;
+
+ // create array to hold results
+ {
+ int dim[2] = {nreps, ntimes-1};
+ PROTECT(F = makearray(2,dim)); nprotect++;
}
- if (ncovars>0) {
- PROTECT(cindex = MATCHCOLNAMES(covar,covarnames)); nprotect++;
- cidx = INTEGER(cindex);
- } else {
- cidx = 0;
- }
- ndim[0] = nvar; ndim[1] = npar; ndim[2] = nrep; ndim[3] = ntimes;
- ndim[4] = covlen; ndim[5] = covdim;
+ switch (mode) {
- if (use_native) set_pomp_userdata(args);
+ case 0: // R function
- euler_densities(ff,REAL(F),REAL(x),REAL(times),REAL(params),
- ndim,sidx,pidx,cidx,
- REAL(tcovar),REAL(covar),INTEGER(log));
+ {
+ double *cp = REAL(cvec);
+ double *t1p = REAL(t1vec);
+ double *t2p = REAL(t2vec);
+ double *x1p = REAL(x1vec);
+ double *x2p = REAL(x2vec);
+ double *pp = REAL(pvec);
+ double *t1s = REAL(times);
+ double *t2s = t1s+1;
+ double *x1s = REAL(x);
+ double *x2s = x1s + nvars*nreps;
+ double *ps;
+ double *fs = REAL(F);
+ int j, k;
- if (use_native) unset_pomp_userdata();
+ for (k = 0; k < ntimes-1; k++, t1s++, t2s++) { // loop over times
+ R_CheckUserInterrupt();
+
+ *t1p = *t1s; *t2p = *t2s;
+
+ // interpolate the covariates at time t1, store the results in cvec
+ table_lookup(&covariate_table,*t1p,cp,0);
+
+ for (j = 0, ps = REAL(params); j < nreps; j++, fs++, x1s += nvars, x2s += nvars, ps += npars) { // loop over replicates
+
+ memcpy(x1p,x1s,nvars*sizeof(double));
+ memcpy(x2p,x2s,nvars*sizeof(double));
+ memcpy(pp,ps,npars*sizeof(double));
+
+ *fs = *(REAL(AS_NUMERIC(eval(fcall,rho))));
+
+ if (!give_log) *fs = exp(*fs);
+
+ }
+ }
+ }
+
+ break;
+
+ case 1: // native code
+
+ set_pomp_userdata(args);
+
+ {
+ double *t1s = REAL(times);
+ double *t2s = t1s+1;
+ double *x1s = REAL(x);
+ double *x2s = x1s + nvars*nreps;
+ double *fs = REAL(F);
+ double *cp = REAL(cvec);
+ double *ps;
+ int j, k;
+
+ for (k = 0; k < ntimes-1; k++, t1s++, t2s++) { // loop over times
+
+ R_CheckUserInterrupt();
+
+ // interpolate the covariates at time t1, store the results in cvec
+ table_lookup(&covariate_table,*t1s,cp,0);
+
+ for (j = 0, ps = REAL(params); j < nreps; j++, fs++, x1s += nvars, x2s += nvars, ps += npars) { // loop over replicates
+
+ (*ff)(fs,x1s,x2s,*t1s,*t2s,ps,sidx,pidx,cidx,ncovars,cp);
+
+ if (!give_log) *fs = exp(*fs);
+
+ }
+ }
+ }
+
+ unset_pomp_userdata();
+
+ break;
+
+ default:
+ error("unrecognized 'mode' in 'euler_model_density'");
+ break;
+
+ }
+
UNPROTECT(nprotect);
return F;
}
-#undef X1VEC
-#undef X2VEC
-#undef PVEC
-#undef CVEC
-#undef TIME1
-#undef TIME2
-#undef NVAR
-#undef NPAR
-#undef RHO
-#undef FCALL
-
int num_euler_steps (double t1, double t2, double *dt) {
double tol = sqrt(DOUBLE_EPS);
int nstep;
Modified: pkg/pomp/src/rmeasure.c
===================================================================
--- pkg/pomp/src/rmeasure.c 2012-04-27 21:51:50 UTC (rev 685)
+++ pkg/pomp/src/rmeasure.c 2012-04-28 18:01:22 UTC (rev 686)
@@ -11,20 +11,15 @@
SEXP do_rmeasure (SEXP object, SEXP x, SEXP times, SEXP params, SEXP fun)
{
int nprotect = 0;
- int first = 1;
int mode = -1;
- int use_names;
int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp, nobs;
SEXP Snames, Pnames, Cnames, Onames;
SEXP statenames, paramnames, covarnames, obsnames;
SEXP tvec, xvec, pvec, cvec;
SEXP fn, fcall, rho, ans, nm;
SEXP Y;
- double *ts, *xs, *ps, *ys;
- double *cp, *xp, *pp, *tp, *yt;
int *dim, ndim[3], *op;
int *sidx, *pidx, *cidx, *oidx;
- int i, j, k;
struct lookup_table covariate_table;
pomp_measure_model_simulator *ff = NULL;
@@ -36,7 +31,6 @@
PROTECT(x = as_state_array(x)); nprotect++;
dim = INTEGER(GET_DIM(x));
nvars = dim[0]; nrepsx = dim[1];
- xs = REAL(x);
if (ntimes != dim[2])
error("rmeasure error: length of 'times' and 3rd dimension of 'x' do not agree");
@@ -44,7 +38,6 @@
PROTECT(params = as_matrix(params)); nprotect++;
dim = INTEGER(GET_DIM(params));
npars = dim[0]; nrepsp = dim[1];
- ps = REAL(params);
nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
@@ -65,7 +58,6 @@
// vector for interpolated covariates
PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
SET_NAMES(cvec,Cnames);
- cp = REAL(cvec);
ndim[0] = nobs; ndim[1] = nreps; ndim[2] = ntimes;
PROTECT(Y = makearray(3,ndim)); nprotect++;
@@ -82,15 +74,10 @@
case 0: // use R function
PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
- tp = REAL(tvec);
-
PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
SET_NAMES(xvec,Snames);
- xp = REAL(xvec);
-
- PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
SET_NAMES(pvec,Pnames);
- pp = REAL(pvec);
// set up the function call
PROTECT(fcall = LCONS(cvec,fcall)); nprotect++;
@@ -133,53 +120,68 @@
case 0: // R function
- for (k = 0, yt = REAL(Y), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
+ {
+ int first = 1;
+ int use_names = 0;
+ double *yt = REAL(Y);
+ double *time = REAL(times);
+ double *tp = REAL(tvec);
+ double *cp = REAL(cvec);
+ double *xp = REAL(xvec);
+ double *pp = REAL(pvec);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *ys;
+ int *posn;
+ int i, j, k;
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++) { // loop over times
- *tp = *ts; // copy the time
- table_lookup(&covariate_table,*tp,cp,0); // interpolate the covariates
+ R_CheckUserInterrupt(); // check for user interrupt
+
+ *tp = *time; // copy the time
+ table_lookup(&covariate_table,*tp,cp,0); // interpolate the covariates
- for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
+ for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
- // copy the states and parameters into place
- for (i = 0; i < nvars; i++) xp[i] = xs[i+nvars*((j%nrepsx)+nrepsx*k)];
- for (i = 0; i < npars; i++) pp[i] = ps[i+npars*(j%nrepsp)];
+ // copy the states and parameters into place
+ for (i = 0; i < nvars; i++) xp[i] = xs[i+nvars*((j%nrepsx)+nrepsx*k)];
+ for (i = 0; i < npars; i++) pp[i] = ps[i+npars*(j%nrepsp)];
- if (first) {
- // evaluate the call
- PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (first) {
+ // evaluate the call
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (LENGTH(ans) != nobs) {
+ error("user 'rmeasure' returns a vector of %d observables but %d are expected: compare 'data' slot?",
+ LENGTH(ans),nobs);
+ }
- if (LENGTH(ans) != nobs) {
- error("user 'rmeasure' returns a vector of %d observables but %d are expected: compare 'data' slot?",
- LENGTH(ans),nobs);
- }
+ // get name information to fix potential alignment problems
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) { // match names against names from data slot
+ posn = INTEGER(PROTECT(matchnames(Onames,nm))); nprotect++;
+ } else {
+ posn = 0;
+ }
- // get name information to fix potential alignment problems
- PROTECT(nm = GET_NAMES(ans)); nprotect++;
- use_names = !isNull(nm);
- if (use_names) { // match names against names from data slot
- op = INTEGER(PROTECT(matchnames(Onames,nm))); nprotect++;
- } else {
- op = 0;
- }
+ ys = REAL(AS_NUMERIC(ans));
- ys = REAL(AS_NUMERIC(ans));
+ first = 0;
- first = 0;
+ } else {
- } else {
+ ys = REAL(AS_NUMERIC(eval(fcall,rho)));
- ys = REAL(AS_NUMERIC(eval(fcall,rho)));
+ }
+ if (use_names) {
+ for (i = 0; i < nobs; i++) yt[posn[i]] = ys[i];
+ } else {
+ for (i = 0; i < nobs; i++) yt[i] = ys[i];
+ }
+
}
-
- if (use_names) {
- for (i = 0; i < nobs; i++) yt[op[i]] = ys[i];
- } else {
- for (i = 0; i < nobs; i++) yt[i] = ys[i];
- }
-
}
}
@@ -187,29 +189,39 @@
case 1: // native routine
- set_pomp_userdata(fcall);
- GetRNGstate();
+ {
+ double *yt = REAL(Y);
+ double *time = REAL(times);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *xp, *pp;
+ int j, k;
- for (k = 0, yt = REAL(Y), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
+ set_pomp_userdata(fcall);
+ GetRNGstate();
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++) { // loop over times
- // interpolate the covar functions for the covariates
- table_lookup(&covariate_table,*ts,cp,0);
+ R_CheckUserInterrupt(); // check for user interrupt
+
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,*time,cp,0);
- for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
+ for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
- xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
- pp = &ps[npars*(j%nrepsp)];
+ xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
+ pp = &ps[npars*(j%nrepsp)];
- (*ff)(yt,xp,pp,oidx,sidx,pidx,cidx,ncovars,cp,*ts);
+ (*ff)(yt,xp,pp,oidx,sidx,pidx,cidx,ncovars,cp,*time);
+ }
}
+
+ PutRNGstate();
+ unset_pomp_userdata();
}
-
- PutRNGstate();
- unset_pomp_userdata();
-
+
break;
default:
More information about the pomp-commits
mailing list