[Pomp-commits] r648 - in pkg: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 10 14:45:52 CEST 2012
Author: kingaa
Date: 2012-04-10 14:45:52 +0200 (Tue, 10 Apr 2012)
New Revision: 648
Modified:
pkg/DESCRIPTION
pkg/R/trajectory-pomp.R
pkg/src/lookup_table.c
pkg/src/skeleton.c
pkg/src/trajectory.c
Log:
- put vectorfield evaluation routines (needed for trajetory computation) into C
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-07 20:47:07 UTC (rev 647)
+++ pkg/DESCRIPTION 2012-04-10 12:45:52 UTC (rev 648)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-3
-Date: 2012-04-07
+Version: 0.41-4
+Date: 2012-04-10
Revision: $Rev$
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>
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2012-04-07 20:47:07 UTC (rev 647)
+++ pkg/R/trajectory-pomp.R 2012-04-10 12:45:52 UTC (rev 648)
@@ -12,8 +12,8 @@
as.data.frame <- as.logical(as.data.frame)
-if (length(times)==0)
- stop("if ",sQuote("times")," is empty, there is no work to do",call.=FALSE)
+ if (length(times)==0)
+ stop(sQuote("times")," is empty, there is no work to do",call.=FALSE)
if (any(diff(times)<0))
stop(sQuote("times")," must be a nondecreasing sequence of times",call.=FALSE)
@@ -23,7 +23,7 @@
else
t0 <- as.numeric(t0)
- if (t0>times[1])
+ if (t0>times[1])
stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=FALSE)
ntimes <- length(times)
@@ -33,22 +33,11 @@
stop("trajectory error: ",sQuote("params")," must be supplied",call.=FALSE)
}
}
- nrep <- NCOL(params)
- if (is.null(dim(params))) {
- params <- matrix(
- params,
- nrow=length(params),
- ncol=nrep,
- dimnames=list(
- names(params),
- NULL
- )
- )
- }
+ params <- as.matrix(params)
+ nrep <- ncol(params)
paramnames <- rownames(params)
if (is.null(paramnames))
stop("trajectory error: ",sQuote("params")," must have rownames",call.=FALSE)
- params <- as.matrix(params)
x0 <- init.state(object,params=params,t0=t0)
nvar <- nrow(x0)
@@ -62,7 +51,7 @@
type <- object at skeleton.type # map or vectorfield?
if (is.na(type))
- stop("trajectory error: no skeleton specified",call.=FALSE)
+ stop("trajectory error: 'skeleton.type' unspecified",call.=FALSE)
if (type=="map") {
Modified: pkg/src/lookup_table.c
===================================================================
--- pkg/src/lookup_table.c 2012-04-07 20:47:07 UTC (rev 647)
+++ pkg/src/lookup_table.c 2012-04-10 12:45:52 UTC (rev 648)
@@ -44,6 +44,7 @@
int flag = 0;
int j, k;
double e;
+ if ((tab == 0) || (tab->length < 1) || (tab->width < 1)) return;
tab->index = findInterval(tab->x,tab->length,x,TRUE,TRUE,tab->index,&flag);
if (flag != 0) // we are extrapolating
warning("table_lookup: extrapolating (flag %d) at %le", flag, x);
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2012-04-07 20:47:07 UTC (rev 647)
+++ pkg/src/skeleton.c 2012-04-10 12:45:52 UTC (rev 648)
@@ -7,307 +7,215 @@
#include <Rinternals.h>
#include <R_ext/Rdynload.h>
-static void eval_skel (pomp_skeleton *vf,
- double *f,
- double *x, double *times, double *params,
- int *ndim,
- int *stateindex, int *parindex, int *covindex,
- double *time_table, double *covar_table)
-{
- double t, *xp, *pp, *fp;
- int nvar = ndim[0];
- int npar = ndim[1];
- int nrepp = ndim[2];
- int nrepx = ndim[3];
- int ntimes = ndim[4];
- int covlen = ndim[5];
- int covdim = ndim[6];
- double covar_fn[covdim];
- int nrep;
- int k, p;
-
- // set up the covariate table
- struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
-
- nrep = (nrepp > nrepx) ? nrepp : nrepx;
-
- for (k = 0; k < ntimes; k++) { // loop over times
-
- R_CheckUserInterrupt(); // check for user interrupt
-
- t = times[k];
-
- // interpolate the covar functions for the covariates
- if (covdim > 0)
- table_lookup(&covariate_table,t,covar_fn,0);
-
- for (p = 0; p < nrep; p++) { // loop over replicates
-
- fp = &f[nvar*(p+nrep*k)];
- xp = &x[nvar*((p%nrepx)+nrepx*k)];
- pp = ¶ms[npar*(p%nrepp)];
-
- (*vf)(fp,xp,pp,stateindex,parindex,covindex,covdim,covar_fn,t);
-
- }
- }
-}
-
-// these global objects will pass the needed information to the user-defined function (see 'default_skel_fn')
-// each of these is allocated once, globally, and refilled many times
-static SEXP _pomp_skel_Xvec; // state variable vector
-static SEXP _pomp_skel_Pvec; // parameter vector
-static SEXP _pomp_skel_Cvec; // covariate vector
-static SEXP _pomp_skel_time; // time
-static int _pomp_skel_nvar; // number of state variables
-static int _pomp_skel_npar; // number of parameters
-static SEXP _pomp_skel_envir; // function's environment
-static SEXP _pomp_skel_fcall; // function call
-static SEXP _pomp_skel_vnames; // names of state variables
-static int *_pomp_skel_vindex; // indices of state variables
-static int _pomp_skel_first; // first evaluation?
-
-#define XVEC (_pomp_skel_Xvec)
-#define PVEC (_pomp_skel_Pvec)
-#define CVEC (_pomp_skel_Cvec)
-#define TIME (_pomp_skel_time)
-#define NVAR (_pomp_skel_nvar)
-#define NPAR (_pomp_skel_npar)
-#define RHO (_pomp_skel_envir)
-#define FCALL (_pomp_skel_fcall)
-#define VNAMES (_pomp_skel_vnames)
-#define VIDX (_pomp_skel_vindex)
-#define FIRST (_pomp_skel_first)
-
-// this is the vectorfield 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_skel_fn (double *f, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t)
-{
- int nprotect = 0;
- int k;
- double *xp;
- SEXP ans, nm, idx;
- int use_names = 0, *op;
-
- xp = REAL(XVEC);
- for (k = 0; k < NVAR; k++) xp[k] = x[k];
- xp = REAL(PVEC);
- for (k = 0; k < NPAR; k++) xp[k] = p[k];
- xp = REAL(CVEC);
- for (k = 0; k < covdim; k++) xp[k] = covar[k];
- xp = REAL(TIME);
- xp[0] = t;
-
- PROTECT(ans = eval(FCALL,RHO)); nprotect++; // evaluate the call
-
- if (FIRST) {
- if (LENGTH(ans)!=NVAR)
- error("user 'skeleton' returns a vector of %d state variables but %d are expected",LENGTH(ans),NVAR);
- PROTECT(nm = GET_NAMES(ans)); nprotect++;
- use_names = !isNull(nm);
- if (use_names) { // match names against known names of states
- PROTECT(idx = matchnames(VNAMES,nm)); nprotect++;
- op = INTEGER(idx);
- for (k = 0; k < NVAR; k++) VIDX[k] = op[k];
- } else {
- VIDX = 0;
- }
- FIRST = 0;
- }
-
- xp = REAL(AS_NUMERIC(ans));
- if (VIDX == 0) {
- for (k = 0; k < NVAR; k++) f[k] = xp[k];
- } else {
- for (k = 0; k < NVAR; k++) f[VIDX[k]] = xp[k];
- }
-
- UNPROTECT(nprotect);
-}
-
SEXP do_skeleton (SEXP object, SEXP x, SEXP t, SEXP params, SEXP fun)
{
int nprotect = 0;
- int *dim, nvars, npars, nrepsp, nrepsx, nreps, ntimes, covlen, covdim;
+ int nvars, npars, nrepp, nrepx, nreps, ntimes, covlen, covdim;
int use_native;
- int nstates, nparams, ncovars;
- double *xp;
- SEXP dimP, dimX, fn, F;
- SEXP tcovar, covar;
+ int fdim[3];
+ SEXP tcovar, covar, fn, Snames, F;
+ double *xs, *ts, *ps, *ft;
+ int *dim;
+
+ SEXP rho, fcall, ans, nm;
+ SEXP tvec, xvec, pvec, cvec;
+ int first = 1;
+ int use_names = 0;
+ int *op;
+
+ int *sidx, *pidx, *cidx;
SEXP statenames, paramnames, covarnames;
- SEXP sindex, pindex, cindex;
- SEXP Pnames, Cnames;
- pomp_skeleton *ff = NULL;
- int k, len;
+ pomp_skeleton *vf;
+ double *tp, *xp, *pp, *cp, *fs;
+ int i, j, k;
+
PROTECT(t = AS_NUMERIC(t)); nprotect++;
ntimes = LENGTH(t);
+ ts = REAL(t);
PROTECT(x = as_state_array(x)); nprotect++;
dim = INTEGER(GET_DIM(x));
- nvars = dim[0]; nrepsx = dim[1];
+ nvars = dim[0]; nrepx = dim[1];
if (ntimes != dim[2])
error("skeleton error: length of 't' and 3rd dimension of 'x' do not agree");
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ xs = REAL(x);
PROTECT(params = as_matrix(params)); nprotect++;
dim = INTEGER(GET_DIM(params));
- npars = dim[0]; nrepsp = dim[1];
+ npars = dim[0]; nrepp = dim[1];
+ ps = REAL(params);
- nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
+ // 2nd dimension of 'x' and 'params' need not agree
+ nreps = (nrepp > nrepx) ? nrepp : nrepx;
+ if ((nreps % nrepp != 0) || (nreps % nrepx != 0))
+ error("skeleton error: 2nd dimensions of 'x' and 'params' are incompatible");
- if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
- error("skeleton error: larger number of replicates is not a multiple of smaller");
-
+ // extract the covariates
PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
-
- PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
- PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
- PROTECT(covarnames = GET_SLOT(object,install("covarnames"))); nprotect++;
- nstates = LENGTH(statenames);
- nparams = LENGTH(paramnames);
- ncovars = LENGTH(covarnames);
-
dim = INTEGER(GET_DIM(covar));
covlen = dim[0]; covdim = dim[1];
- PROTECT(VNAMES = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
- PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
+ // set up the covariate table
+ struct lookup_table covariate_table = {covlen, covdim, 0, REAL(tcovar), REAL(covar)};
+
+ // set up the array to be returned
+ fdim[0] = nvars; fdim[1] = nreps; fdim[2] = ntimes;
+ PROTECT(F = makearray(3,fdim)); nprotect++;
+ setrownames(F,Snames,3);
+
+ // extract the user-defined function
PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
use_native = INTEGER(VECTOR_ELT(fun,1))[0];
+ PROTECT(cvec = NEW_NUMERIC(covdim)); nprotect++;
+ SET_NAMES(cvec,GET_COLNAMES(GET_DIMNAMES(covar)));
+ cp = REAL(cvec);
+
+ // first do setup
switch (use_native) {
- case 0: // R skeleton
- ff = (pomp_skeleton *) default_skel_fn;
- PROTECT(RHO = (CLOENV(fn))); nprotect++;
- NVAR = nvars; // for internal use
- NPAR = npars; // for internal use
- PROTECT(TIME = NEW_NUMERIC(1)); nprotect++; // for internal use
- PROTECT(XVEC = NEW_NUMERIC(nvars)); nprotect++; // for internal use
- xp = REAL(XVEC);
- for (k = 0; k < nvars; k++) xp[k] = 0.0;
- PROTECT(PVEC = NEW_NUMERIC(npars)); nprotect++; // for internal use
- PROTECT(CVEC = NEW_NUMERIC(covdim)); nprotect++; // for internal use
- SET_NAMES(XVEC,VNAMES); // 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
+ case 0: // R skeleton
+
+ 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(pvec,GET_ROWNAMES(GET_DIMNAMES(params)));
+ pp = REAL(pvec);
+
// set up the function call
- PROTECT(FCALL = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
- PROTECT(FCALL = LCONS(CVEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("covars"));
- PROTECT(FCALL = LCONS(PVEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("params"));
- PROTECT(FCALL = LCONS(TIME,FCALL)); nprotect++;
- SET_TAG(FCALL,install("t"));
- PROTECT(FCALL = LCONS(XVEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("x"));
- PROTECT(FCALL = LCONS(fn,FCALL)); nprotect++;
- VIDX = (int *) R_alloc(nvars,sizeof(int));
- FIRST = 1;
+ PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ PROTECT(fcall = LCONS(cvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("covars"));
+ PROTECT(fcall = LCONS(pvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("params"));
+ PROTECT(fcall = LCONS(tvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("t"));
+ PROTECT(fcall = LCONS(xvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("x"));
+ PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
+
break;
+
case 1: // native skeleton
- ff = (pomp_skeleton *) R_ExternalPtrAddr(fn);
- VIDX = 0;
+
+ vf = (pomp_skeleton *) R_ExternalPtrAddr(fn);
+
+ PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
+ if (LENGTH(statenames) > 0) {
+ sidx = INTEGER(PROTECT(MATCHROWNAMES(x,statenames))); nprotect++;
+ } else {
+ sidx = 0;
+ }
+
+ PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
+ if (LENGTH(paramnames) > 0) {
+ pidx = INTEGER(PROTECT(MATCHROWNAMES(params,paramnames))); nprotect++;
+ } else {
+ pidx = 0;
+ }
+
+ PROTECT(covarnames = GET_SLOT(object,install("covarnames"))); nprotect++;
+ if (LENGTH(covarnames) > 0) {
+ cidx = INTEGER(PROTECT(MATCHCOLNAMES(covar,covarnames))); nprotect++;
+ } else {
+ cidx = 0;
+ }
+
break;
+
default:
error("unrecognized 'use' slot in 'skeleton'");
break;
}
- {
- int fdim[3];
- fdim[0] = nvars; fdim[1] = nreps; fdim[2] = ntimes;
- PROTECT(F = makearray(3,fdim)); nprotect++;
- setrownames(F,VNAMES,3);
- xp = REAL(F);
- // for (k = 0, len = nvars*nreps*ntimes; k < len; k++) xp[k] = 0.0;
- }
- if (nstates > 0) {
- PROTECT(sindex = MATCHROWNAMES(x,statenames)); nprotect++;
- } else {
- PROTECT(sindex = NEW_INTEGER(0)); nprotect++;
- }
- if (nparams > 0) {
- PROTECT(pindex = MATCHROWNAMES(params,paramnames)); nprotect++;
- } else {
- PROTECT(pindex = NEW_INTEGER(0)); nprotect++;
- }
- if (ncovars > 0) {
- PROTECT(cindex = MATCHCOLNAMES(covar,covarnames)); nprotect++;
- } else {
- PROTECT(cindex = NEW_INTEGER(0)); nprotect++;
- }
+ // now do computations
+ switch (use_native) {
+ case 0: // R skeleton
- {
- int ndim[7];
- ndim[0] = nvars; ndim[1] = npars; ndim[2] = nrepsp; ndim[3] = nrepsx;
- ndim[4] = ntimes; ndim[5] = covlen; ndim[6] = covdim;
+ for (k = 0, ft = REAL(F); k < ntimes; k++, ts++) { // loop over times
- eval_skel(ff,REAL(F),REAL(x),REAL(t),REAL(params),
- ndim,INTEGER(sindex),INTEGER(pindex),INTEGER(cindex),
- REAL(tcovar),REAL(covar));
- }
+ R_CheckUserInterrupt(); // check for user interrupt
+
+ *tp = *ts; // copy the time
- UNPROTECT(nprotect);
- return F;
-}
+ // interpolate the covar functions for the covariates
+ if (covdim > 0) table_lookup(&covariate_table,*ts,cp,0);
+
+ for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
+
+ for (i = 0; i < nvars; i++) xp[i] = xs[i+nvars*((j%nrepx)+nrepx*k)];
+ for (i = 0; i < npars; i++) pp[i] = ps[i+npars*(j%nrepp)];
+
+ if (first) {
+
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (LENGTH(ans)!=nvars)
+ error("user 'skeleton' returns a vector of %d state variables but %d are expected",LENGTH(ans),nvars);
-#undef XVEC
-#undef PVEC
-#undef CVEC
-#undef TIME
-#undef NVAR
-#undef NPAR
-#undef RHO
-#undef FCALL
-#undef VNAMES
+ // get name information to fix possible alignment problems
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) {
+ op = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
+ } else {
+ op = 0;
+ }
+
+ fs = REAL(AS_NUMERIC(ans));
+
+ first = 0;
+
+ } else {
+
+ fs = REAL(AS_NUMERIC(eval(fcall,rho)));
+ }
+
+ if (use_names)
+ for (i = 0; i < nvars; i++) ft[op[i]] = fs[i];
+ else
+ for (i = 0; i < nvars; i++) ft[i] = fs[i];
+
+ }
+ }
-static SEXP *_pomp_vf_eval_object;
-static SEXP *_pomp_vf_eval_params;
-static SEXP *_pomp_vf_eval_skelfun;
-static SEXP *_pomp_vf_eval_xnames;
-static int _pomp_vf_eval_xdim[3];
-#define OBJECT (_pomp_vf_eval_object)
-#define PARAMS (_pomp_vf_eval_params)
-#define SKELFUN (_pomp_vf_eval_skelfun)
-#define XNAMES (_pomp_vf_eval_xnames)
-#define XDIM (_pomp_vf_eval_xdim)
+ break;
-void pomp_desolve_init (SEXP object, SEXP params, SEXP fun, SEXP statenames, SEXP nvar, SEXP nrep) {
- OBJECT = &object;
- PARAMS = ¶ms;
- SKELFUN = &fun;
- XNAMES = &statenames;
- XDIM[0] = INTEGER(AS_INTEGER(nvar))[0];
- XDIM[1] = INTEGER(AS_INTEGER(nrep))[0];
- XDIM[2] = 1;
-}
+ case 1: // native skeleton
+
+ for (k = 0, ft = REAL(F); k < ntimes; k++, ts++) { // loop over times
+ R_CheckUserInterrupt(); // check for user interrupt
+
+ // interpolate the covar functions for the covariates
+ if (covdim > 0) table_lookup(&covariate_table,*ts,cp,0);
+
+ for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
+
+ xp = &xs[nvars*((j%nrepx)+nrepx*k)];
+ pp = &ps[npars*(j%nrepp)];
+
+ (*vf)(ft,xp,pp,sidx,pidx,cidx,covdim,cp,*ts);
+
+ }
+ }
-void pomp_vf_eval (int *neq, double *t, double *y, double *ydot, double *yout, int *ip)
-{
- SEXP T, X, dXdt;
- int dim[3];
-
- PROTECT(T = NEW_NUMERIC(1));
- PROTECT(X = makearray(3,XDIM));
- setrownames(X,*(XNAMES),3);
- REAL(T)[0] = *t;
- memcpy(REAL(X),y,(*neq)*sizeof(double));
- PROTECT(dXdt = do_skeleton(*(OBJECT),X,T,*(PARAMS),*(SKELFUN)));
- memcpy(ydot,REAL(dXdt),(*neq)*sizeof(double));
-
- UNPROTECT(3);
-}
+ break;
-#undef XDIM
-#undef XNAMES
-#undef SKELFUN
-#undef PARAMS
-#undef OBJECT
+ default:
+ error("unrecognized 'use' slot in 'skeleton'");
+ break;
+ }
+
+ UNPROTECT(nprotect);
+ return F;
+}
Modified: pkg/src/trajectory.c
===================================================================
--- pkg/src/trajectory.c 2012-04-07 20:47:07 UTC (rev 647)
+++ pkg/src/trajectory.c 2012-04-10 12:45:52 UTC (rev 648)
@@ -34,7 +34,7 @@
{
int nprotect = 0;
SEXP ans;
- SEXP x, f, skel, time, zeronames, zindex, dt;
+ SEXP x, skel, time, zeronames;
int nvars, nreps, ntimes;
int nzeros;
int nsteps;
@@ -44,8 +44,7 @@
int *zidx;
double deltat;
- PROTECT(t0 = AS_NUMERIC(t0)); nprotect++;
- PROTECT(x = duplicate(AS_NUMERIC(x0))); nprotect++;
+ PROTECT(x = as_state_array(duplicate(AS_NUMERIC(x0)))); nprotect++;
xp = REAL(x);
PROTECT(times = AS_NUMERIC(times)); nprotect++;
@@ -57,9 +56,8 @@
if (nreps != INTEGER(GET_DIM(params))[1])
error("mismatch in dimensions of 'x0' and 'params'");
- PROTECT(time = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(time = duplicate(AS_NUMERIC(t0))); nprotect++;
tm = REAL(time);
- *tm = REAL(t0)[0];
ndim[0] = nvars; ndim[1] = nreps; ndim[2] = ntimes;
PROTECT(ans = makearray(3,ndim)); nprotect++;
@@ -67,14 +65,12 @@
ap = REAL(ans);
PROTECT(skel = get_pomp_fun(GET_SLOT(object,install("skeleton")))); nprotect++;
- PROTECT(dt = GET_SLOT(object,install("skelmap.delta.t"))); nprotect++;
- deltat = REAL(dt)[0];
+ deltat = *(REAL(GET_SLOT(object,install("skelmap.delta.t"))));
PROTECT(zeronames = GET_SLOT(object,install("zeronames"))); nprotect++;
nzeros = LENGTH(zeronames);
if (nzeros>0) {
- PROTECT(zindex = MATCHROWNAMES(x0,zeronames)); nprotect++;
- zidx = INTEGER(zindex);
+ zidx = INTEGER(PROTECT(MATCHROWNAMES(x0,zeronames))); nprotect++;
} else {
zidx = 0;
}
@@ -92,26 +88,58 @@
nsteps = num_map_steps(*tm,*tp,deltat);
for (h = 0; h < nsteps; h++) {
- PROTECT(f = do_skeleton(object,x,time,params,skel));
- fp = REAL(f);
- for (j = 0; j < nreps; j++) {
- for (i = 0; i < nvars; i++) {
- xp[i+nvars*j] = fp[i+nvars*j];
- }
- }
- UNPROTECT(1);
+ fp = REAL(do_skeleton(object,x,time,params,skel));
+ for (j = 0; j < nreps; j++)
+ for (i = 0; i < nvars; i++) xp[i+nvars*j] = fp[i+nvars*j];
*tm += deltat;
}
+
for (j = 0; j < nreps; j++) {
- for (i = 0; i < nvars; i++) {
- ap[i+nvars*(j+nreps*k)] = xp[i+nvars*j];
- }
- for (i = 0; i < nzeros; i++) {
- xp[zidx[i]+nvars*j] = 0.0;
- }
+ for (i = 0; i < nvars; i++) ap[i+nvars*(j+nreps*k)] = xp[i+nvars*j];
+ for (i = 0; i < nzeros; i++) xp[zidx[i]+nvars*j] = 0.0;
}
+
}
UNPROTECT(nprotect);
return ans;
}
+
+static struct {
+ SEXP *object;
+ SEXP *params;
+ SEXP *skelfun;
+ SEXP *xnames;
+ int xdim[3];
+} _pomp_vf_eval_common;
+
+
+#define COMMON(X) (_pomp_vf_eval_common.X)
+
+void pomp_desolve_init (SEXP object, SEXP params, SEXP fun, SEXP statenames, SEXP nvar, SEXP nrep) {
+ COMMON(object) = &object;
+ COMMON(params) = ¶ms;
+ COMMON(skelfun) = &fun;
+ COMMON(xnames) = &statenames;
+ COMMON(xdim)[0] = INTEGER(AS_INTEGER(nvar))[0];
+ COMMON(xdim)[1] = INTEGER(AS_INTEGER(nrep))[0];
+ COMMON(xdim)[2] = 1;
+}
+
+
+void pomp_vf_eval (int *neq, double *t, double *y, double *ydot, double *yout, int *ip)
+{
+ SEXP T, X, dXdt;
+
+ PROTECT(T = NEW_NUMERIC(1));
+ PROTECT(X = makearray(3,COMMON(xdim)));
+ setrownames(X,*(COMMON(xnames)),3);
+ REAL(T)[0] = *t;
+ memcpy(REAL(X),y,(*neq)*sizeof(double));
+ PROTECT(dXdt = do_skeleton(*(COMMON(object)),X,T,*(COMMON(params)),*(COMMON(skelfun))));
+ memcpy(ydot,REAL(dXdt),(*neq)*sizeof(double));
+
+ UNPROTECT(3);
+}
+
+#undef COMMON
More information about the pomp-commits
mailing list