[Pomp-commits] r685 - in pkg/pomp: . R data inst/doc src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 27 23:51:50 CEST 2012
Author: kingaa
Date: 2012-04-27 23:51:50 +0200 (Fri, 27 Apr 2012)
New Revision: 685
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/pomp-fun.R
pkg/pomp/R/pomp.R
pkg/pomp/data/bbs.rda
pkg/pomp/data/blowflies.rda
pkg/pomp/data/dacca.rda
pkg/pomp/data/euler.sir.rda
pkg/pomp/data/gillespie.sir.rda
pkg/pomp/data/gompertz.rda
pkg/pomp/data/ou2.rda
pkg/pomp/data/ricker.rda
pkg/pomp/data/rw2.rda
pkg/pomp/data/verhulst.rda
pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
pkg/pomp/inst/doc/bsmc-ricker-flat-prior.rda
pkg/pomp/inst/doc/gompertz-multi-mif.rda
pkg/pomp/inst/doc/gompertz-trajmatch.rda
pkg/pomp/inst/doc/intro_to_pomp.pdf
pkg/pomp/inst/doc/nlf-block-boot.rda
pkg/pomp/inst/doc/nlf-boot.rda
pkg/pomp/inst/doc/nlf-fit-from-truth.rda
pkg/pomp/inst/doc/nlf-fits.rda
pkg/pomp/inst/doc/nlf-lag-tests.rda
pkg/pomp/inst/doc/nlf-multi-short.rda
pkg/pomp/inst/doc/ricker-mif.rda
pkg/pomp/inst/doc/ricker-probe-match.rda
pkg/pomp/src/dmeasure.c
pkg/pomp/src/euler.c
pkg/pomp/src/lookup_table.c
pkg/pomp/src/partrans.c
pkg/pomp/src/pomp_fun.c
pkg/pomp/src/pomp_internal.h
pkg/pomp/src/rmeasure.c
pkg/pomp/src/skeleton.c
Log:
- changes to the guts of 'rmeasure', 'dmeasure', 'partrans', and the Euler rprocess plugins
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-26 14:43:38 UTC (rev 684)
+++ pkg/pomp/DESCRIPTION 2012-04-27 21:51:50 UTC (rev 685)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.42-1
-Date: 2012-04-26
+Date: 2012-04-27
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/R/pomp-fun.R
===================================================================
--- pkg/pomp/R/pomp-fun.R 2012-04-26 14:43:38 UTC (rev 684)
+++ pkg/pomp/R/pomp-fun.R 2012-04-27 21:51:50 UTC (rev 685)
@@ -5,10 +5,14 @@
R.fun = 'function',
native.fun = 'character',
PACKAGE = 'character',
- use = 'integer'
+ mode = 'integer'
)
)
+pomp_native_code <- 2L
+pomp_R_function <- 1L
+pomp_undef_mode <- -1L
+
## constructor
pomp.fun <- function (f = NULL, PACKAGE, proto = NULL) {
if (missing(PACKAGE)) PACKAGE <- ""
@@ -29,7 +33,7 @@
R.fun=f,
native.fun=character(0),
PACKAGE=PACKAGE,
- use=1L
+ mode=pomp_R_function
)
} else if (is.character(f)) {
retval <- new(
@@ -37,7 +41,7 @@
R.fun=function(...)stop("unreachable error: please report this bug!"),
native.fun=f,
PACKAGE=PACKAGE,
- use=2L
+ mode=pomp_native_code
)
} else {
retval <- new(
@@ -45,7 +49,7 @@
R.fun=function(...)stop(sQuote(fname)," not specified"),
native.fun=character(0),
PACKAGE=PACKAGE,
- use=-1L
+ mode=pomp_undef_mode
)
}
retval
@@ -55,10 +59,10 @@
'show',
'pomp.fun',
function (object) {
- use <- object at use
- if (use==1L) {
+ mode <- object at mode
+ if (mode==pomp_R_function) {
show(object at R.fun)
- } else if (use==2L) {
+ } else if (mode==pomp_native_code) {
cat("native function ",sQuote(object at native.fun),sep="")
if (length(object at PACKAGE)>0)
cat(", dynamically loaded from ",sQuote(object at PACKAGE),sep="")
@@ -75,14 +79,18 @@
function (x, ...) show(x)
)
-get.pomp.fun <- function (object) {
- use <- object at use
- if (use==1L) {
- f <- object at R.fun
- } else if (use==2L) {
- f <- getNativeSymbolInfo(name=object at native.fun,PACKAGE=object at PACKAGE)$address
- } else {
- stop("function not specified")
- }
- list(f,as.integer(use-1))
-}
+get.pomp.fun <- function (object)
+ .Call(get_pomp_fun,object)
+
+
+## get.pomp.fun <- function (object) {
+## mode <- object at mode
+## if (mode==1L) {
+## f <- object at R.fun
+## } else if (mode==2L) {
+## f <- getNativeSymbolInfo(name=object at native.fun,PACKAGE=object at PACKAGE)$address
+## } else {
+## stop("function not specified")
+## }
+## list(f,as.integer(mode-1))
+## }
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2012-04-26 14:43:38 UTC (rev 684)
+++ pkg/pomp/R/pomp.R 2012-04-27 21:51:50 UTC (rev 685)
@@ -209,17 +209,17 @@
if (nrow(covar)>0) {
if (
- (skeleton at use==1)
+ (skeleton at mode==1)
&&!("covars"%in%names(formals(skeleton at R.fun)))
)
warning("a covariate table has been given, yet the ",sQuote("skeleton")," function does not have ",sQuote("covars")," as a formal argument",call.=FALSE)
if (
- (rmeasure at use==1)
+ (rmeasure at mode==1)
&&!("covars"%in%names(formals(rmeasure at R.fun)))
)
warning("a covariate table has been given, yet the ",sQuote("rmeasure")," function does not have ",sQuote("covars")," as a formal argument",call.=FALSE)
if (
- (dmeasure at use==1)
+ (dmeasure at mode==1)
&&!("covars"%in%names(formals(dmeasure at R.fun)))
)
warning("a covariate table has been given, yet the ",sQuote("dmeasure")," function does not have ",sQuote("covars")," as a formal argument",call.=FALSE)
@@ -263,8 +263,8 @@
}
if (
has.trans &&
- par.trans at use<0 &&
- par.untrans at use<0
+ par.trans at mode<0 &&
+ par.untrans at mode<0
)
has.trans <- FALSE
Modified: pkg/pomp/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/bsmc-ricker-flat-prior.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/gompertz-multi-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/gompertz-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-block-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-fit-from-truth.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-fits.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-lag-tests.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/nlf-multi-short.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/ricker-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/ricker-probe-match.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/src/dmeasure.c
===================================================================
--- pkg/pomp/src/dmeasure.c 2012-04-26 14:43:38 UTC (rev 684)
+++ pkg/pomp/src/dmeasure.c 2012-04-27 21:51:50 UTC (rev 685)
@@ -4,133 +4,40 @@
#include <Rmath.h>
#include <Rdefines.h>
#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
#include "pomp_internal.h"
-static void dens_meas (pomp_measure_model_density *f,
- double *lik, double *y,
- double *x, double *times, double *params,
- int *give_log, int *ndim,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- double *time_table, double *covar_table)
-{
- double t, *lp, *xp, *pp, *yp;
- 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];
- int nobs = ndim[7];
- 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];
- yp = &y[nobs*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
-
- lp = &lik[p+nrep*k];
- xp = &x[nvar*((p%nrepx)+nrepx*k)];
- pp = ¶ms[npar*(p%nrepp)];
-
- (*f)(lp,yp,xp,pp,*give_log,obsindex,stateindex,parindex,covindex,covdim,covar_fn,t);
-
- }
- }
-}
-
-// these global objects will pass the needed information to the user-defined function
-// each of these is allocated once, globally, and refilled many times
-static SEXP _pomp_dmeas_Yvec; // observable vector
-static SEXP _pomp_dmeas_Xvec; // state variable vector
-static SEXP _pomp_dmeas_Pvec; // parameter vector
-static SEXP _pomp_dmeas_Cvec; // covariate vector
-static SEXP _pomp_dmeas_time; // time
-static int _pomp_dmeas_nvar; // number of state variables
-static int _pomp_dmeas_npar; // number of parameters
-static int _pomp_dmeas_nobs; // number of observables
-static SEXP _pomp_dmeas_envir; // function's environment
-static SEXP _pomp_dmeas_fcall; // function call
-
-#define YVEC (_pomp_dmeas_Yvec)
-#define XVEC (_pomp_dmeas_Xvec)
-#define PVEC (_pomp_dmeas_Pvec)
-#define CVEC (_pomp_dmeas_Cvec)
-#define TIME (_pomp_dmeas_time)
-#define NVAR (_pomp_dmeas_nvar)
-#define NPAR (_pomp_dmeas_npar)
-#define NOBS (_pomp_dmeas_nobs)
-#define RHO (_pomp_dmeas_envir)
-#define FCALL (_pomp_dmeas_fcall)
-
-// this is the measurement pdf evaluated when the user supplies an R function
-// (and not a native routine)
-// Note that obsindex, stateindex, parindex, covindex are ignored.
-static void default_meas_dens (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t)
-{
- int nprotect = 0;
- int k;
- double *xp;
- SEXP ans;
- xp = REAL(YVEC);
- for (k = 0; k < NOBS; k++) xp[k] = y[k];
- 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 = AS_NUMERIC(eval(FCALL,RHO))); nprotect++; // evaluate the call
- if (LENGTH(ans) != 1)
- error("dmeasure error: user 'dmeasure' must return a scalar");
- xp = REAL(ans);
- *lik = *xp;
- UNPROTECT(nprotect);
-}
-
SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log, SEXP fun)
{
int nprotect = 0;
- int *dim, nvars, npars, nrepsp, nrepsx, nreps, ntimes, covlen, covdim, nobs;
- int ndim[7];
- SEXP F, fn;
- SEXP tcovar, covar;
+ int first = 1;
+ int mode = -1;
+ int give_log;
+ int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp, nobs;
+ SEXP Snames, Pnames, Cnames, Onames;
SEXP statenames, paramnames, covarnames, obsnames;
- SEXP sindex, pindex, cindex, oindex;
+ SEXP tvec, xvec, yvec, pvec, cvec;
+ SEXP fn, fcall, rho, ans;
+ SEXP F;
+ double *ts, *xs, *ps, *ys, *fs;
+ double *cp, *xp, *pp, *tp, *yp;
+ double *ft;
+ int *dim, ndim[2];
int *sidx, *pidx, *cidx, *oidx;
- SEXP Xnames, Ynames, Pnames, Cnames;
- int use_native;
- int nstates, nparams, ncovars, nobsers;
+ int i, j, k;
+ struct lookup_table covariate_table;
pomp_measure_model_density *ff = NULL;
PROTECT(times = AS_NUMERIC(times)); nprotect++;
- ntimes = LENGTH(times);
+ ntimes = length(times);
if (ntimes < 1)
- error("dmeasure error: no work to do");
+ error("dmeasure error: length('times') = 0, no work to do");
PROTECT(y = as_matrix(y)); nprotect++;
dim = INTEGER(GET_DIM(y));
nobs = dim[0];
+ ys = REAL(y);
if (ntimes != dim[1])
error("dmeasure error: length of 'times' and 2nd dimension of 'y' do not agree");
@@ -138,6 +45,7 @@
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("dmeasure error: length of 'times' and 3rd dimension of 'x' do not agree");
@@ -145,129 +53,175 @@
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;
if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
error("dmeasure error: larger number of replicates is not a multiple of smaller");
- PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
- PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
- PROTECT(obsnames = GET_SLOT(object,install("obsnames"))); 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++;
- nobsers = LENGTH(obsnames);
- nstates = LENGTH(statenames);
- nparams = LENGTH(paramnames);
- ncovars = LENGTH(covarnames);
-
- dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; covdim = dim[1];
- PROTECT(Ynames = GET_ROWNAMES(GET_DIMNAMES(y))); nprotect++;
- PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ PROTECT(Onames = GET_ROWNAMES(GET_DIMNAMES(y))); nprotect++;
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
+ PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(GET_SLOT(object,install("covar"))))); nprotect++;
+
+ give_log = *(INTEGER(log));
- PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
- use_native = INTEGER(VECTOR_ELT(fun,1))[0];
+ // set up the covariate table
+ covariate_table = make_covariate_table(object,&ncovars);
- switch (use_native) {
- case 0: // use R function
- ff = (pomp_measure_model_density *) default_meas_dens;
- PROTECT(RHO = (CLOENV(fn))); nprotect++;
- NVAR = nvars; // for internal use
- NPAR = npars; // for internal use
- NOBS = nobs; // for internal use
- PROTECT(TIME = NEW_NUMERIC(1)); nprotect++; // for internal use
- PROTECT(YVEC = NEW_NUMERIC(nobs)); nprotect++; // for internal use
- PROTECT(XVEC = NEW_NUMERIC(nvars)); nprotect++; // for internal use
- PROTECT(PVEC = NEW_NUMERIC(npars)); nprotect++; // for internal use
- PROTECT(CVEC = NEW_NUMERIC(covdim)); nprotect++; // for internal use
- SET_NAMES(YVEC,Ynames); // make sure the names attribute is copied
- SET_NAMES(XVEC,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
+ // vector for interpolated covariates
+ PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
+ SET_NAMES(cvec,Cnames);
+ cp = REAL(cvec);
+
+ ndim[0] = nreps; ndim[1] = ntimes;
+ PROTECT(F = makearray(2,ndim)); nprotect++;
+
+ // extract the user-defined function
+ PROTECT(fn = unpack_pomp_fun(fun,&mode)); nprotect++;
+
+ // extract 'userdata' as pairlist
+ PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+
+ // first do setup
+ switch (mode) {
+
+ case 0: // R function
+
+ PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
+ tp = REAL(tvec);
+
+ PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
+ SET_NAMES(xvec,Snames);
+ xp = REAL(xvec);
+
+ PROTECT(yvec = NEW_NUMERIC(nobs)); nprotect++;
+ SET_NAMES(yvec,Onames);
+ yp = REAL(yvec);
+
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
+ SET_NAMES(pvec,Pnames);
+ 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(AS_LOGICAL(log),FCALL)); nprotect++;
- SET_TAG(FCALL,install("log"));
- 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(YVEC,FCALL)); nprotect++;
- SET_TAG(FCALL,install("y"));
- PROTECT(FCALL = LCONS(fn,FCALL)); nprotect++;
+ PROTECT(fcall = LCONS(cvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("covars"));
+ PROTECT(fcall = LCONS(AS_LOGICAL(log),fcall)); nprotect++;
+ SET_TAG(fcall,install("log"));
+ 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(yvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("y"));
+ PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
+
+ // get the function's environment
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
+
break;
- case 1: // use native routine
+
+ case 1: // native code
+
+ // construct state, parameter, covariate, observable indices
+ oidx = INTEGER(PROTECT(name_index(Onames,object,"obsnames"))); nprotect++;
+ sidx = INTEGER(PROTECT(name_index(Snames,object,"statenames"))); nprotect++;
+ pidx = INTEGER(PROTECT(name_index(Pnames,object,"paramnames"))); nprotect++;
+ cidx = INTEGER(PROTECT(name_index(Cnames,object,"covarnames"))); nprotect++;
+
+ // address of native routine
ff = (pomp_measure_model_density *) R_ExternalPtrAddr(fn);
+
break;
+
default:
- error("unrecognized 'use' slot in 'dmeasure'");
+
+ error("unrecognized 'mode' slot in 'dmeasure'");
break;
+
}
- ndim[0] = nreps; ndim[1] = ntimes;
- PROTECT(F = makearray(2,ndim)); nprotect++;
+ // now do computations
+ switch (mode) {
- if (nobsers > 0) {
- PROTECT(oindex = matchnames(Ynames,obsnames)); nprotect++;
- oidx = INTEGER(oindex);
- } else {
- oidx = 0;
- }
+ case 0: // R function
- if (nstates > 0) {
- PROTECT(sindex = matchnames(Xnames,statenames)); nprotect++;
- sidx = INTEGER(sindex);
- } else {
- sidx = 0;
- }
+ for (k = 0, ft = REAL(F), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
- if (nparams > 0) {
- PROTECT(pindex = matchnames(Pnames,paramnames)); nprotect++;
- pidx = INTEGER(pindex);
- } else {
- pidx = 0;
- }
-
- if (ncovars > 0) {
- PROTECT(cindex = matchnames(Cnames,covarnames)); nprotect++;
- cidx = INTEGER(cindex);
- } else {
- cidx = 0;
- }
-
- ndim[0] = nvars; ndim[1] = npars; ndim[2] = nrepsp; ndim[3] = nrepsx; ndim[4] = ntimes;
- ndim[5] = covlen; ndim[6] = covdim; ndim[7] = nobs;
+ R_CheckUserInterrupt(); // check for user interrupt
- if (use_native) {
- PROTECT(FCALL = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
- set_pomp_userdata(FCALL);
- }
+ *tp = *ts; // copy the time
+ table_lookup(&covariate_table,*tp,cp,0); // interpolate the covariates
- dens_meas(ff,REAL(F),REAL(y),REAL(x),REAL(times),REAL(params),INTEGER(log),ndim,
- oidx,sidx,pidx,cidx,REAL(tcovar),REAL(covar));
-
- if (use_native) {
+ for (i = 0; i < nobs; i++) yp[i] = ys[i+nobs*k];
+
+ for (j = 0; j < nreps; j++, ft++) { // 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)];
+
+ if (first) {
+ // evaluate the call
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+
+ if (LENGTH(ans) != 1)
+ error("user 'dmeasure' return a vector of length %d when it should return a scalar",LENGTH(ans));
+
+ fs = REAL(AS_NUMERIC(ans));
+
+ first = 0;
+
+ } else {
+
+ fs = REAL(AS_NUMERIC(eval(fcall,rho)));
+
+ }
+
+ *ft = *fs;
+
+ }
+ }
+
+ break;
+
+ case 1: // native code
+
+ set_pomp_userdata(fcall);
+
+ for (k = 0, ft = REAL(F), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
+
+ R_CheckUserInterrupt(); // check for user interrupt
+
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,*ts,cp,0);
+
+ yp = &ys[nobs*k];
+
+ for (j = 0; j < nreps; j++, ft++) { // loop over replicates
+
+ xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
+ pp = &ps[npars*(j%nrepsp)];
+
+ (*ff)(ft,yp,xp,pp,give_log,oidx,sidx,pidx,cidx,ncovars,cp,*ts);
+
+ }
+ }
+
unset_pomp_userdata();
+
+ break;
+
+ default:
+
+ error("unrecognized 'mode' slot in 'dmeasure'");
+ break;
+
}
UNPROTECT(nprotect);
return F;
}
-
-#undef YVEC
-#undef XVEC
-#undef PVEC
-#undef CVEC
-#undef TIME
-#undef NVAR
-#undef NPAR
-#undef NOBS
-#undef RHO
-#undef FCALL
Modified: pkg/pomp/src/euler.c
===================================================================
--- pkg/pomp/src/euler.c 2012-04-26 14:43:38 UTC (rev 684)
+++ pkg/pomp/src/euler.c 2012-04-27 21:51:50 UTC (rev 685)
@@ -3,465 +3,256 @@
#include "pomp_internal.h"
#include <R_ext/Constants.h>
-int num_euler_steps (double t1, double t2, double *dt) {
- double tol = sqrt(DOUBLE_EPS);
- int nstep;
- // nstep will be the number of Euler steps to take in going from t1 to t2.
- // note also that the stepsize changes.
- // this choice is meant to be conservative
- // (i.e., so that the actual dt does not exceed the specified dt
- // by more than the relative tolerance 'tol')
- // and to counteract roundoff error.
- // It seems to work well, but is not guaranteed:
- // suggestions would be appreciated.
-
- if (t1 >= t2) {
- *dt = 0.0;
- nstep = 0;
- } else if (t1+*dt >= t2) {
- *dt = t2-t1;
- nstep = 1;
- } else {
- nstep = (int) ceil((t2-t1)/(*dt)/(1+tol));
- *dt = (t2-t1)/((double) nstep);
- }
- return nstep;
-}
-
-int num_map_steps (double t1, double t2, double dt) {
- double tol = sqrt(DOUBLE_EPS);
- int nstep;
- // nstep will be the number of discrete-time steps to take in going from t1 to t2.
- nstep = (int) floor((t2-t1)/dt/(1-tol));
- return (nstep > 0) ? nstep : 0;
-}
-
-// take Euler-Poisson steps of size at most deltat from t1 to t2
-static void euler_simulator (pomp_onestep_sim *estep,
- double *x, double *xstart, double *times, double *params,
- int *ndim, double *deltat,
- int *stateindex, int *parindex, int *covindex, int *zeroindex,
- double *time_table, double *covar_table)
+SEXP euler_model_simulator (SEXP func,
+ SEXP xstart, SEXP times, SEXP params,
+ SEXP deltat, SEXP method,
+ SEXP statenames, SEXP paramnames, SEXP covarnames, SEXP zeronames,
+ SEXP tcovar, SEXP covar, SEXP args)
{
- double t, *xp, *pp;
- int nvar = ndim[0];
- int npar = ndim[1];
- int nrep = ndim[2];
- int ntimes = ndim[3];
- int covlen = ndim[4];
- int covdim = ndim[5];
- int nzero = ndim[6];
- double covar_fn[covdim];
- int j, k, p, step, nstep = 0;
- double dt;
+ 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 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;
- struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
+ 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);
- // copy the start values into the result array
- for (p = 0; p < nrep; p++)
- for (k = 0; k < nvar; k++)
- x[k+nvar*p] = xstart[k+nvar*p];
-
- // loop over times
- for (step = 1; step < ntimes; step++) {
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(xstart))); nprotect++;
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
- R_CheckUserInterrupt();
+ // set up the covariate table
+ struct lookup_table covariate_table = {covlen, ncovars, 0, REAL(tcovar), REAL(covar)};
- t = times[step-1];
- dt = *deltat;
+ // vector for interpolated covariates
+ PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
+ SET_NAMES(cvec,Cnames);
+ cp = REAL(cvec);
- if (t > times[step]) {
- error("'times' is not an increasing sequence");
- }
+ // indices of accumulator variables
+ nzeros = LENGTH(zeronames);
+ zidx = INTEGER(PROTECT(matchnames(Snames,zeronames))); nprotect++;
- nstep = num_euler_steps(t,times[step],&dt);
+ // extract user function
+ PROTECT(fn = pomp_fun_handler(func,&mode)); nprotect++;
+
+ // set up
+ switch (mode) {
- for (p = 0; p < nrep; p++) {
- xp = &x[nvar*(p+nrep*step)];
- // copy in the previous values of the state variables
- for (k = 0; k < nvar; k++)
- xp[k] = x[k+nvar*(p+nrep*(step-1))];
- // set some variables to zero
- for (k = 0; k < nzero; k++)
- xp[zeroindex[k]] = 0.0;
- }
+ case 1: // native code
- for (j = 0; j < nstep; j++) { // loop over Euler steps
+ // 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++;
- // interpolate the covar functions for the covariates
- if (covdim > 0)
- table_lookup(&covariate_table,t,covar_fn,0);
+ ff = (pomp_onestep_sim *) R_ExternalPtrAddr(fn);
- for (p = 0; p < nrep; p++) { // loop over replicates
-
- pp = ¶ms[npar*p];
- xp = &x[nvar*(p+nrep*step)];
+ break;
- (*estep)(xp,pp,stateindex,parindex,covindex,covdim,covar_fn,t,dt);
+ case 0: // R function
- }
+ // get function's environment
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
- t += dt;
+ 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(pvec,Pnames);
+ pp = REAL(pvec);
- if (j == nstep-2) { // penultimate step
- dt = times[step]-t;
- t = times[step]-dt;
- }
+ // set up the function call
+ PROTECT(fcall = LCONS(cvec,args)); nprotect++;
+ SET_TAG(fcall,install("covars"));
+ PROTECT(fcall = LCONS(dtvec,fcall)); nprotect++;
+ SET_TAG(fcall,install("delta.t"));
+ 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++;
- }
- }
-}
+ // to hold indices of variables that must be rearranged
+ posn = (int *) R_alloc(nvars,sizeof(int));
-// take discrete steps of size deltat from t1 to t2
-static void discrete_time_simulator (pomp_onestep_sim *estep,
- double *x, double *xstart, double *times, double *params,
- int *ndim, double *deltat,
- int *stateindex, int *parindex, int *covindex, int *zeroindex,
- double *time_table, double *covar_table)
-{
- double t, *xp, *pp;
- int nvar = ndim[0];
- int npar = ndim[1];
- int nrep = ndim[2];
- int ntimes = ndim[3];
- int covlen = ndim[4];
- int covdim = ndim[5];
- int nzero = ndim[6];
- double covar_fn[covdim];
- int j, k, p, step, nstep = 0;
- double dt;
+ // 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;
- struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
+ 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) {
- // copy the start values into the result array
- for (p = 0; p < nrep; p++)
- for (k = 0; k < nvar; k++)
- x[k+nvar*p] = xstart[k+nvar*p];
-
- dt = *deltat;
- t = times[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);
+ }
- // loop over times
- for (step = 1; step < ntimes; step++) {
+ 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];
+ }
- R_CheckUserInterrupt();
+ xs = REAL(AS_NUMERIC(ans));
- nstep = num_map_steps(t,times[step],dt);
+ first = 0;
- for (p = 0; p < nrep; p++) {
- xp = &x[nvar*(p+nrep*step)];
- // copy in the previous values of the state variables
- for (k = 0; k < nvar; k++)
- xp[k] = x[k+nvar*(p+nrep*(step-1))];
- // set some variables to zero
- for (k = 0; k < nzero; k++)
- xp[zeroindex[k]] = 0.0;
- }
+ } else {
+
+ xs = REAL(AS_NUMERIC(eval(fcall,rho)));
- for (j = 0; j < nstep; j++) { // loop over steps
+ }
- // 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
+ 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];
+ }
- pp = ¶ms[npar*p];
- xp = &x[nvar*(p+nrep*step)];
+ UNPROTECT(nprotect);
+ }
- (*estep)(xp,pp,stateindex,parindex,covindex,covdim,covar_fn,t,dt);
+ sidx = 0;
+ pidx = 0;
+ cidx = 0;
- }
+ ff = (pomp_onestep_sim *) R_step_fn;
- t += dt;
+ break;
- }
+ default:
+ error("unrecognized 'mode' in 'euler_simulator'");
+ break;
}
-}
-// take one step from t1 to t2
-static void onestep_simulator (pomp_onestep_sim *estep,
- double *x, double *xstart, double *times, double *params,
- int *ndim,
- int *stateindex, int *parindex, int *covindex, int *zeroindex,
- double *time_table, double *covar_table)
-{
- double t, *xp, *pp;
- int nvar = ndim[0];
- int npar = ndim[1];
- int nrep = ndim[2];
- int ntimes = ndim[3];
- int covlen = ndim[4];
- int covdim = ndim[5];
- int nzero = ndim[6];
- double covar_fn[covdim];
- int k, p, step;
- double dt;
+ if (mode==1) {
+ set_pomp_userdata(args);
+ GetRNGstate();
+ }
- struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
+ // 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);
// copy the start values into the result array
- for (p = 0; p < nrep; p++)
- for (k = 0; k < nvar; k++)
- x[k+nvar*p] = xstart[k+nvar*p];
-
+ xt = REAL(xstart);
+ for (j = 0; j < nreps; j++)
+ for (i = 0; i < nvars; i++)
+ Xt[i+nvars*j] = xt[i+nvars*j];
+
+ meth = *(INTEGER(AS_INTEGER(method))); // 0 = Euler, 1 = one-step, 2 = fixed step
+
+ // now do computations
// loop over times
+ time = REAL(times);
+ t = time[0];
+
for (step = 1; step < ntimes; step++) {
R_CheckUserInterrupt();
- t = times[step-1];
- dt = times[step]-t;
+ if (t > time[step]) {
+ error("'times' is not an increasing sequence");
+ }
- // interpolate the covar functions for the covariates
- if (covdim > 0)
- table_lookup(&covariate_table,t,covar_fn,0);
-
- for (p = 0; p < nrep; p++) {
- xp = &x[nvar*(p+nrep*step)];
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 685
More information about the pomp-commits
mailing list