[Pomp-commits] r688 - in pkg/pomp: . R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 30 20:49:33 CEST 2012
Author: kingaa
Date: 2012-04-30 20:49:32 +0200 (Mon, 30 Apr 2012)
New Revision: 688
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/inst/ChangeLog
pkg/pomp/man/builder.Rd
pkg/pomp/src/trajectory.c
Log:
- rework the guts of 'iterate_map' for much faster trajectory calculation
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-29 14:24:52 UTC (rev 687)
+++ pkg/pomp/DESCRIPTION 2012-04-30 18:49:32 UTC (rev 688)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.42-1
-Date: 2012-04-28
+Date: 2012-04-30
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/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-04-29 14:24:52 UTC (rev 687)
+++ pkg/pomp/R/builder.R 2012-04-30 18:49:32 UTC (rev 688)
@@ -20,9 +20,9 @@
}
pompBuilder <- function (data, times, t0, name,
- statenames, paramnames, zeronames,
+ statenames, paramnames,
rmeasure, dmeasure, step.fn, step.fn.delta.t,
- skeleton, skeleton.type, skelmap.delta.t = 1,
+ skeleton, skeleton.type, skelmap.delta.t = 1, ...,
link = TRUE) {
obsnames <- names(data)
obsnames <- setdiff(obsnames,times)
@@ -37,7 +37,6 @@
skeleton=skeleton
)
if (link) pompLink(name)
- if (missing(zeronames)) zeronames <- character(0)
pomp(
data=data,times=times,t0=t0,
rprocess=euler.sim(
@@ -54,7 +53,7 @@
obsnames=obsnames,
statenames=statenames,
paramnames=paramnames,
- zeronames=zeronames
+ ...
)
}
Modified: pkg/pomp/inst/ChangeLog
===================================================================
--- pkg/pomp/inst/ChangeLog 2012-04-29 14:24:52 UTC (rev 687)
+++ pkg/pomp/inst/ChangeLog 2012-04-30 18:49:32 UTC (rev 688)
@@ -1,5 +1,35 @@
+2012-04-29 kingaa
+
+ * [r687] src/dmeasure.c, src/skeleton.c: - more work on the guts
+
+2012-04-28 kingaa
+
+ * [r686] DESCRIPTION, src/euler.c, src/rmeasure.c: - more changes
+ to the guts
+
+2012-04-27 kingaa
+
+ * [r685] DESCRIPTION, R/pomp-fun.R, R/pomp.R, data/bbs.rda,
+ data/blowflies.rda, data/dacca.rda, data/euler.sir.rda,
+ data/gillespie.sir.rda, data/gompertz.rda, data/ou2.rda,
+ data/ricker.rda, data/rw2.rda, data/verhulst.rda,
+ inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/bsmc-ricker-flat-prior.rda,
+ inst/doc/gompertz-multi-mif.rda, inst/doc/gompertz-trajmatch.rda,
+ inst/doc/intro_to_pomp.pdf, inst/doc/nlf-block-boot.rda,
+ inst/doc/nlf-boot.rda, inst/doc/nlf-fit-from-truth.rda,
+ inst/doc/nlf-fits.rda, inst/doc/nlf-lag-tests.rda,
+ inst/doc/nlf-multi-short.rda, inst/doc/ricker-mif.rda,
+ inst/doc/ricker-probe-match.rda, src/dmeasure.c, src/euler.c,
+ src/lookup_table.c, src/partrans.c, src/pomp_fun.c,
+ src/pomp_internal.h, src/rmeasure.c, src/skeleton.c: - changes to
+ the guts of 'rmeasure', 'dmeasure', 'partrans', and the Euler
+ rprocess plugins
+
2012-04-26 kingaa
+ * [r684] R/builder.R, demo/gompertz.R, inst/ChangeLog: - zeronames
+ argument in 'pompBuilder' is optional
* [r683] DESCRIPTION, NAMESPACE, R/builder.R, data/bbs.rda,
data/blowflies.rda, data/dacca.rda, data/euler.sir.rda,
data/gillespie.sir.rda, data/gompertz.rda, data/ou2.rda,
Modified: pkg/pomp/man/builder.Rd
===================================================================
--- pkg/pomp/man/builder.Rd 2012-04-29 14:24:52 UTC (rev 687)
+++ pkg/pomp/man/builder.Rd 2012-04-30 18:49:32 UTC (rev 688)
@@ -5,9 +5,10 @@
\code{pompBuilder} is an EXPERIMENTAL facility for producing compiled \code{pomp} objects.
}
\usage{
-pompBuilder(data, times, t0, name, statenames, paramnames, zeronames,
+pompBuilder(data, times, t0, name, statenames, paramnames,
rmeasure, dmeasure, step.fn, step.fn.delta.t,
- skeleton, skeleton.type, skelmap.delta.t = 1, link = TRUE)
+ skeleton, skeleton.type, skelmap.delta.t = 1,
+ \dots, link = TRUE)
}
\arguments{
\item{data, times, t0}{
@@ -18,8 +19,8 @@
\item{name}{
character; the stem of the name for the files that will be produced.
}
- \item{statenames, paramnames, zeronames}{
- names of state-variables, parameters, and accumulator variables, respectively
+ \item{statenames, paramnames}{
+ names of state-variables and parameters, respectively
}
\item{rmeasure, dmeasure}{
C codes implementing the measurement model
@@ -33,6 +34,9 @@
As in \code{pomp}, \code{skeleton.type} indicates whether the skeleton is a map (discrete-time) or vectorfield (continuous-time).
If the former, \code{skelmap.delta.t} is the time-step of the map.
}
+ \item{\dots}{
+ additional arguments are passed to \code{\link{pomp}}
+ }
\item{link}{
logical; if TRUE, the resulting code will be linked after compilation.
}
Modified: pkg/pomp/src/trajectory.c
===================================================================
--- pkg/pomp/src/trajectory.c 2012-04-29 14:24:52 UTC (rev 687)
+++ pkg/pomp/src/trajectory.c 2012-04-30 18:49:32 UTC (rev 688)
@@ -5,101 +5,252 @@
#include "pomp_internal.h"
-// copy t(x[-1,-1]) -> y[,rep,]
-SEXP traj_transp_and_copy (SEXP y, SEXP x, SEXP rep) {
+SEXP iterate_map (SEXP object, SEXP times, SEXP t0, SEXP x0, SEXP params)
+{
int nprotect = 0;
- SEXP ans = R_NilValue;
- int nvars, nreps, ntimes;
- int i, j, k, m, n;
- double *xp, *yp;
+ SEXP fn, fcall, rho, ans, nm;
+ SEXP X;
+ SEXP Snames, Pnames, Cnames;
+ SEXP zeronames;
+ SEXP cvec, tvec, xvec, pvec;
+ int nvars, npars, nreps, nrepx, nrepp, ntimes, ncovars, nzeros, nsteps;
+ int mode = -1;
+ int *dim;
+ int *zidx = 0, *sidx = 0, *pidx = 0, *cidx = 0;
+ pomp_skeleton *ff = NULL;
+ struct lookup_table covariate_table;
- j = INTEGER(rep)[0]-1;
- nvars = INTEGER(GET_DIM(y))[0];
- nreps = INTEGER(GET_DIM(y))[1];
- ntimes = INTEGER(GET_DIM(y))[2];
- n = INTEGER(GET_DIM(x))[0];
- m = nvars*nreps;
+ PROTECT(x0 = as_matrix(x0)); nprotect++;
+ dim = INTEGER(GET_DIM(x0));
+ nvars = dim[0]; nrepx = dim[1];
- for (i = 0, xp = REAL(x)+n+1; i < nvars; i++, xp += n) {
- for (k = 0, yp = REAL(y)+i+nvars*j; k < ntimes; k++, yp += m) {
- *yp = xp[k];
- }
- }
+ PROTECT(params = as_matrix(params)); nprotect++;
+ dim = INTEGER(GET_DIM(params));
+ npars = dim[0]; nrepp = dim[1];
- UNPROTECT(nprotect);
- return ans;
-}
+ // 2nd dimension of 'x0' and 'params' need not agree
+ nreps = (nrepp > nrepx) ? nrepp : nrepx;
+ if ((nreps % nrepp != 0) || (nreps % nrepx != 0))
+ error("trajectory error: 2nd dimensions of 'x0' and 'params' are incompatible");
-SEXP iterate_map (SEXP object, SEXP times, SEXP t0, SEXP x0, SEXP params)
-{
- int nprotect = 0;
- SEXP ans;
- SEXP x, skel, time, zeronames;
- int nvars, nreps, ntimes;
- int nzeros;
- int nsteps;
- int h, i, j, k;
- int ndim[3];
- double *tm, *tp, *xp, *fp, *ap;
- int *zidx;
- double deltat;
-
- PROTECT(x = as_state_array(duplicate(AS_NUMERIC(x0)))); nprotect++;
- xp = REAL(x);
-
PROTECT(times = AS_NUMERIC(times)); nprotect++;
ntimes = LENGTH(times);
- tp = REAL(times);
- nvars = INTEGER(GET_DIM(x0))[0];
- nreps = INTEGER(GET_DIM(x0))[1];
- if (nreps != INTEGER(GET_DIM(params))[1])
- error("mismatch in dimensions of 'x0' and 'params'");
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x0))); nprotect++;
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(GET_SLOT(object,install("covar"))))); nprotect++;
+
+ // set up the covariate table
+ covariate_table = make_covariate_table(object,&ncovars);
- PROTECT(time = duplicate(AS_NUMERIC(t0))); nprotect++;
- tm = REAL(time);
+ // vector for interpolated covariates
+ PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
+ SET_NAMES(cvec,Cnames);
- ndim[0] = nvars; ndim[1] = nreps; ndim[2] = ntimes;
- PROTECT(ans = makearray(3,ndim)); nprotect++;
- setrownames(ans,GET_ROWNAMES(GET_DIMNAMES(x0)),3);
- ap = REAL(ans);
+ // extract user-defined function
+ PROTECT(fn = pomp_fun_handler(GET_SLOT(object,install("skeleton")),&mode)); nprotect++;
- PROTECT(skel = get_pomp_fun(GET_SLOT(object,install("skeleton")))); nprotect++;
- deltat = *(REAL(GET_SLOT(object,install("skelmap.delta.t"))));
+ // extract 'userdata' as pairlist
+ PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ // get names and indices of accumulator variables
PROTECT(zeronames = GET_SLOT(object,install("zeronames"))); nprotect++;
nzeros = LENGTH(zeronames);
- if (nzeros>0) {
- zidx = INTEGER(PROTECT(MATCHROWNAMES(x0,zeronames))); nprotect++;
- } else {
- zidx = 0;
+ if (nzeros > 0) {
+ zidx = INTEGER(PROTECT(matchnames(Snames,zeronames))); nprotect++;
}
- if (nzeros>0) {
- for (j = 0; j < nreps; j++)
- for (i = 0; i < nzeros; i++) xp[zidx[i]+nvars*j] = 0.0;
+ // set up the computations
+ switch (mode) {
+
+ case 0: // R function
+
+ PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
+ SET_NAMES(xvec,Snames);
+ SET_NAMES(pvec,Pnames);
+
+ // set up the function call
+ 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++;
+
+ // get function's environment
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
+
+ break;
+
+ case 1: // native skeleton
+
+ // construct state, parameter, covariate, observable indices
+ 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++;
+
+ ff = (pomp_skeleton *) R_ExternalPtrAddr(fn);
+
+ break;
+
+ default:
+ error("unrecognized 'mode' in 'iterate_map'");
+ break;
}
- for (k = 0; k < ntimes; k++, tp++) {
+ // create array to store results
+ {
+ int dim[3] = {nvars, nreps, ntimes};
+ PROTECT(X = makearray(3,dim)); nprotect++;
+ setrownames(X,Snames,3);
+ }
- nsteps = num_map_steps(*tm,*tp,deltat);
+ switch (mode) {
- for (h = 0; h < nsteps; h++) {
- fp = REAL(do_skeleton(object,x,time,params,skel));
+ case 0: // R function
+
+ {
+ int first = 1;
+ int use_names;
+ double *time = REAL(times);
+ double *x0s = REAL(x0);
+ double *Xt = REAL(X);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *tp = REAL(tvec);
+ double *xp = REAL(xvec);
+ double *pp = REAL(pvec);
+ double *fs, *xm;
+ double deltat = *(REAL(GET_SLOT(object,install("skelmap.delta.t"))));
+ double t = *(REAL(AS_NUMERIC(t0)));
+ int *posn = 0;
+ int h, i, j, k;
+
for (j = 0; j < nreps; j++)
- for (i = 0; i < nvars; i++) xp[i+nvars*j] = fp[i+nvars*j];
- *tm += deltat;
+ for (i = 0; i < nvars; i++)
+ Xt[i+nvars*j] = x0s[i+nvars*(j%nrepx)];
+
+ for (k = 0; k < ntimes; k++, time++, Xt += nvars*nreps) {
+
+ R_CheckUserInterrupt();
+
+ nsteps = num_map_steps(t,*time,deltat);
+
+ for (i = 0; i < nzeros; i++)
+ for (j = 0, xm = &Xt[zidx[i]]; j < nreps; j++, xm += nvars)
+ *xm = 0.0;
+
+ for (h = 0; h < nsteps; h++) {
+
+ table_lookup(&covariate_table,t,cp,0);
+
+ for (j = 0, xm = Xt; j < nreps; j++, xm += nvars) {
+
+ *tp = t;
+ memcpy(xp,xm,nvars*sizeof(double));
+ memcpy(pp,&ps[npars*(j%nrepp)],npars*sizeof(double));
+
+ 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);
+
+ // get name information to fix possible alignment problems
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) {
+ posn = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
+ }
+
+ fs = REAL(AS_NUMERIC(ans));
+
+ first = 0;
+
+ } else {
+
+ fs = REAL(AS_NUMERIC(eval(fcall,rho)));
+
+ }
+
+ if (use_names)
+ for (i = 0; i < nvars; i++) xm[posn[i]] = fs[i];
+ else
+ for (i = 0; i < nvars; i++) xm[i] = fs[i];
+
+ }
+
+ t += deltat;
+
+ }
+
+ if (k+1 < ntimes)
+ memcpy(Xt+nvars*nreps,Xt,nvars*nreps*sizeof(double));
+
+ }
}
+ break;
- 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;
+ case 1: // native code
+
+ {
+ double *time = REAL(times);
+ double *Xt = REAL(X);
+ double *ps = REAL(params);
+ double *x0s = REAL(x0);
+ double *cp = REAL(cvec);
+ double deltat = *(REAL(GET_SLOT(object,install("skelmap.delta.t"))));
+ double t = *(REAL(AS_NUMERIC(t0)));
+ double *xs = (double *) R_alloc(nvars*nreps,sizeof(double)); // array to temporarily hold states
+ double *fp, *xp;
+ int h, i, j, k;
+
+ // initialize the state matrix
+ for (j = 0; j < nreps; j++)
+ for (i = 0; i < nvars; i++)
+ xs[i+nvars*j] = x0s[i+nvars*(j%nrepx)];
+
+ for (k = 0; k < ntimes; k++, time++, Xt += nvars*nreps) {
+
+ R_CheckUserInterrupt();
+
+ for (i = 0; i < nzeros; i++)
+ for (j = 0, xp = &Xt[zidx[i]]; j < nreps; j++, xp += nvars)
+ *xp = 0.0;
+
+ nsteps = num_map_steps(t,*time,deltat);
+
+ for (h = 0; h < nsteps; h++) {
+
+ table_lookup(&covariate_table,t,cp,0);
+
+ for (j = 0, fp = Xt, xp = xs; j < nreps; j++, fp += nvars, xp += nvars) {
+ (*ff)(fp,xp,&ps[npars*(j%nrepp)],sidx,pidx,cidx,ncovars,cp,t);
+ memcpy(xp,fp,nvars*sizeof(double));
+ }
+
+ t += deltat;
+
+ }
+
+ memcpy(Xt,xs,nvars*nreps*sizeof(double));
+
+ }
}
+ break;
+ default:
+ error("unrecognized 'mode' in 'iterate_map'");
+ break;
}
UNPROTECT(nprotect);
- return ans;
+ return X;
}
static struct {
@@ -149,3 +300,29 @@
}
#undef COMMON
+
+// copy t(x[-1,-1]) -> y[,rep,]
+SEXP traj_transp_and_copy (SEXP y, SEXP x, SEXP rep) {
+ int nprotect = 0;
+ SEXP ans = R_NilValue;
+ int nvars, nreps, ntimes;
+ int i, j, k, m, n;
+ double *xp, *yp;
+
+ j = INTEGER(rep)[0]-1;
+ nvars = INTEGER(GET_DIM(y))[0];
+ nreps = INTEGER(GET_DIM(y))[1];
+ ntimes = INTEGER(GET_DIM(y))[2];
+ n = INTEGER(GET_DIM(x))[0];
+ m = nvars*nreps;
+
+ for (i = 0, xp = REAL(x)+n+1; i < nvars; i++, xp += n) {
+ for (k = 0, yp = REAL(y)+i+nvars*j; k < ntimes; k++, yp += m) {
+ *yp = xp[k];
+ }
+ }
+
+ UNPROTECT(nprotect);
+ return ans;
+}
+
More information about the pomp-commits
mailing list