[Pomp-commits] r98 - in pkg: . data inst inst/doc inst/examples inst/include man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 17 17:22:34 CEST 2009
Author: kingaa
Date: 2009-04-17 17:22:33 +0200 (Fri, 17 Apr 2009)
New Revision: 98
Modified:
pkg/DESCRIPTION
pkg/data/ou2.rda
pkg/inst/ChangeLog
pkg/inst/doc/compiled_code_in_pomp.Rnw
pkg/inst/doc/compiled_code_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/examples/euler_sir.R
pkg/inst/examples/euler_sir.c
pkg/inst/include/pomp.h
pkg/man/euler.Rd
pkg/src/euler.c
pkg/src/euler_sir.c
pkg/src/ou2.c
pkg/src/pomp.h
pkg/src/rmeasure.c
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
The 'euler' and 'rmeasure' evaluation codes now call GetRNGstate() and PutRNGstate() when the user-supplied function they will evaluate is a native routine. This means that potentially many fewer calls to these expensive functions can be made in, e.g., particle filtering or simulation of a large number of replicates and consequently significant performance gains.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/DESCRIPTION 2009-04-17 15:22:33 UTC (rev 98)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.23-1
-Date: 2009-04-14
+Version: 0.23-2
+Date: 2009-04-17
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/inst/ChangeLog 2009-04-17 15:22:33 UTC (rev 98)
@@ -1,5 +1,15 @@
+2009-04-16 kingaa
+
+ * [r97] R/mif-methods.R, R/mif.R: use rowMeans and colMeans in
+ place of apply where appropriate
+ * [r96] man/mif-class.Rd, man/mif-methods.Rd, man/mif.Rd: update
+ references for 'mif'
+ * [r95] man/mif.Rd: fix order of arguments in documentation for
+ 'mif'
+
2009-04-14 kingaa
+ * [r94] inst/ChangeLog:
* [r93] DESCRIPTION, R/mif.R, data/ou2.rda, inst/ChangeLog,
inst/doc/compiled_code_in_pomp.pdf, inst/doc/intro_to_pomp.Rnw,
inst/doc/intro_to_pomp.pdf, man/mif.Rd, tests/ou2-mif.R,
Modified: pkg/inst/doc/compiled_code_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-04-17 15:22:33 UTC (rev 98)
@@ -64,6 +64,7 @@
likewise, we let $\tau\tau^T$ be that of $\varepsilon_{t}$.
Since many of the methods we will use require us to simulate the process and/or measurement models many times, it is a good idea to use native (compiled) codes for the computational heavy lifting.
+This can result in many-fold speedup.
The package includes some C codes that were written to implement the OU example.
Read the source (file `ou2.c') for details.
Modified: pkg/inst/doc/compiled_code_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/euler_sir.R
===================================================================
--- pkg/inst/examples/euler_sir.R 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/inst/examples/euler_sir.R 2009-04-17 15:22:33 UTC (rev 98)
@@ -144,7 +144,7 @@
covar=basis,
delta.t=1/52/20,
statenames=c("S","I","R","cases","W","B","dW"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
covarnames=c("seas1"),
zeronames=c("cases"),
step.fun="sir_euler_simulator",
@@ -152,8 +152,9 @@
dens.fun="sir_euler_density",
dprocess=euler.density,
skeleton.vectorfield="sir_ODE",
+ rmeasure="binom_rmeasure",
+ dmeasure="binom_dmeasure",
PACKAGE="euler_sir", ## name of the shared-object library
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
p <- exp(params)
with(
Modified: pkg/inst/examples/euler_sir.c
===================================================================
--- pkg/inst/examples/euler_sir.c 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/inst/examples/euler_sir.c 2009-04-17 15:22:33 UTC (rev 98)
@@ -29,6 +29,7 @@
#define LOGBETA (p[parindex[3]]) // transmission rate
#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
#define LOGPOPSIZE (p[parindex[5]]) // population size
+#define LOGRHO (p[parindex[6]]) // reporting probability
#define SUSC (x[stateindex[0]]) // number of susceptibles
#define INFD (x[stateindex[1]]) // number of infectives
@@ -40,6 +41,22 @@
#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
+#define MEASLES (y[0])
+
+void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dbinom(MEASLES,CASE,exp(LOGRHO),give_log);
+}
+
+void binom_rmeasure (double *y, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ MEASLES = rbinom(CASE,exp(LOGRHO));
+}
+
+#undef MEASLES
+
// SIR model with Euler multinomial step
// forced transmission (basis functions passed as covariates)
// constant population size as a parameter
@@ -79,8 +96,6 @@
!(R_FINITE(W)))
return;
- GetRNGstate(); // initialize R's random number generator
-
if (beta_sd > 0.0) { // environmental noise is ON
dW = rgamma(dt/beta_var,beta_var); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
if (!(R_FINITE(dW))) return;
@@ -103,8 +118,6 @@
reulermultinom(2,INFD,&rate[3],dt,&trans[3]);
reulermultinom(1,RCVD,&rate[5],dt,&trans[5]);
- PutRNGstate(); // finished with R's random number generator
-
// balance the equations
SUSC += trans[0]-trans[1]-trans[2];
INFD += trans[1]-trans[3]-trans[4];
@@ -116,10 +129,6 @@
}
-#undef BIRTHS
-#undef W
-#undef dW
-
#define DSDT (f[stateindex[0]])
#define DIDT (f[stateindex[1]])
#define DRDT (f[stateindex[2]])
@@ -172,10 +181,13 @@
#undef DRDT
#undef DCDT
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
+#undef BIRTHS
+#undef dW
#define SUSC (x1[stateindex[0]]) // number of susceptibles
#define INFD (x1[stateindex[1]]) // number of infectives
@@ -247,19 +259,20 @@
*f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
}
-#undef GAMMA
-#undef MU
-#undef IOTA
-#undef BETA
-#undef BETA_SD
-#undef POPSIZE
+#undef SEASBASIS
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
#undef BIRTHS
-#undef DW
+#undef dW
-#undef SEASBASIS
+#undef LOGGAMMA
+#undef LOGMU
+#undef LOGIOTA
+#undef LOGBETA
+#undef LOGBETA_SD
+#undef LOGPOPSIZE
+#undef LOGRHO
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/inst/include/pomp.h 2009-04-17 15:22:33 UTC (rev 98)
@@ -27,7 +27,12 @@
// dt = size (duration) of the Euler step
// on output:
// x = contains the new state vector (i.e., at time t+dt)
+//
+// NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function.
+// The RNG is initialized before any call to this function, and the RNG state is written afterward.
+// Inclusion of these calls in the user-defined function may result in significant slowdown.
+
// Prototype for one-step Euler PDF, as used by "euler.density":
typedef void euler_step_pdf(double *f,
double *x1, double *x2, double t1, double t2, const double *p,
@@ -99,6 +104,10 @@
// t = time at the beginning of the Euler step
// on output:
// y = pointer to vector containing simulated observations (length = nobs = nrow(data))
+//
+// NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function.
+// The RNG is initialized before any call to this function, and the RNG state is written afterward.
+// Inclusion of these calls in the user-defined function may result in significant slowdown.
// Prototype for measurement model density evaluator
Modified: pkg/man/euler.Rd
===================================================================
--- pkg/man/euler.Rd 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/man/euler.Rd 2009-04-17 15:22:33 UTC (rev 98)
@@ -132,7 +132,7 @@
covar=basis,
delta.t=1/52/20,
statenames=c("S","I","R","cases","W","B","dW"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
covarnames=c("seas1"),
zeronames=c("cases"),
step.fun="sir_euler_simulator",
@@ -140,7 +140,8 @@
skeleton="sir_ODE",
rprocess=euler.simulate,
dprocess=euler.density,
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
+ dmeasure="binom_dmeasure",
+ rmeasure="binom_rmeasure",
initializer=function(params,t0,...){
p <- exp(params)
with(
Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/src/euler.c 2009-04-17 15:22:33 UTC (rev 98)
@@ -23,7 +23,7 @@
double dt, tol;
struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
-
+
tol = sqrt(DOUBLE_EPS); // relative tolerance in choosing Euler stepsize
// copy the start values into the result array
@@ -179,6 +179,7 @@
SEXP tcovar, SEXP covar, SEXP args)
{
int nprotect = 0;
+ int use_native = 0;
int *dim, xdim[3], ndim[7];
int nvar, npar, nrep, ntimes;
int covlen, covdim;
@@ -201,9 +202,18 @@
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
if (inherits(func,"NativeSymbol")) {
+ use_native = 1;
+ } else if (isFunction(func)) {
+ use_native = 0;
+ } else {
+ UNPROTECT(nprotect);
+ error("illegal input: supplied function must be either an R function or a compiled native function");
+ }
+
+ if (use_native) {
ff = (euler_step_sim *) R_ExternalPtrAddr(func);
VINDEX = 0;
- } else if (isFunction(func)) {
+ } else {
PROTECT(fn = func); nprotect++;
PROTECT(RHO = (CLOENV(fn))); nprotect++;
NVAR = nvar; // for internal use
@@ -231,9 +241,6 @@
ff = (euler_step_sim *) default_euler_step_fn;
VINDEX = (int *) Calloc(nvar,int);
FIRST = 1;
- } else {
- UNPROTECT(nprotect);
- error("illegal input: supplied function must be either an R function or a compiled native function");
}
xdim[0] = nvar; xdim[1] = nrep; xdim[2] = ntimes;
@@ -268,10 +275,14 @@
ndim[0] = nvar; ndim[1] = npar; ndim[2] = nrep; ndim[3] = ntimes;
ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nzeros;
+ if (use_native) GetRNGstate();
+
euler_simulator(ff,REAL(X),REAL(xstart),REAL(times),REAL(params),
ndim,REAL(dt),sidx,pidx,cidx,zidx,
REAL(tcovar),REAL(covar));
+ if (use_native) PutRNGstate();
+
if (VINDEX != 0) Free(VINDEX);
VINDEX = 0;
Modified: pkg/src/euler_sir.c
===================================================================
--- pkg/src/euler_sir.c 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/src/euler_sir.c 2009-04-17 15:22:33 UTC (rev 98)
@@ -29,6 +29,7 @@
#define LOGBETA (p[parindex[3]]) // transmission rate
#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
#define LOGPOPSIZE (p[parindex[5]]) // population size
+#define LOGRHO (p[parindex[6]]) // reporting probability
#define SUSC (x[stateindex[0]]) // number of susceptibles
#define INFD (x[stateindex[1]]) // number of infectives
@@ -40,6 +41,22 @@
#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
+#define MEASLES (y[0])
+
+void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dbinom(MEASLES,CASE,exp(LOGRHO),give_log);
+}
+
+void binom_rmeasure (double *y, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ MEASLES = rbinom(CASE,exp(LOGRHO));
+}
+
+#undef MEASLES
+
// SIR model with Euler multinomial step
// forced transmission (basis functions passed as covariates)
// constant population size as a parameter
@@ -79,8 +96,6 @@
!(R_FINITE(W)))
return;
- GetRNGstate(); // initialize R's random number generator
-
if (beta_sd > 0.0) { // environmental noise is ON
dW = rgamma(dt/beta_var,beta_var); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
if (!(R_FINITE(dW))) return;
@@ -103,8 +118,6 @@
reulermultinom(2,INFD,&rate[3],dt,&trans[3]);
reulermultinom(1,RCVD,&rate[5],dt,&trans[5]);
- PutRNGstate(); // finished with R's random number generator
-
// balance the equations
SUSC += trans[0]-trans[1]-trans[2];
INFD += trans[1]-trans[3]-trans[4];
@@ -116,10 +129,6 @@
}
-#undef BIRTHS
-#undef W
-#undef dW
-
#define DSDT (f[stateindex[0]])
#define DIDT (f[stateindex[1]])
#define DRDT (f[stateindex[2]])
@@ -172,10 +181,13 @@
#undef DRDT
#undef DCDT
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
+#undef BIRTHS
+#undef dW
#define SUSC (x1[stateindex[0]]) // number of susceptibles
#define INFD (x1[stateindex[1]]) // number of infectives
@@ -247,19 +259,20 @@
*f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
}
-#undef GAMMA
-#undef MU
-#undef IOTA
-#undef BETA
-#undef BETA_SD
-#undef POPSIZE
+#undef SEASBASIS
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
#undef BIRTHS
-#undef DW
+#undef dW
-#undef SEASBASIS
+#undef LOGGAMMA
+#undef LOGMU
+#undef LOGIOTA
+#undef LOGBETA
+#undef LOGBETA_SD
+#undef LOGPOPSIZE
+#undef LOGRHO
Modified: pkg/src/ou2.c
===================================================================
--- pkg/src/ou2.c 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/src/ou2.c 2009-04-17 15:22:33 UTC (rev 98)
@@ -107,10 +107,8 @@
double t)
{
double sd = fabs(TAU);
- GetRNGstate();
Y1 = rnorm(X1,sd);
Y2 = rnorm(X2,sd);
- PutRNGstate();
}
#undef X1
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/src/pomp.h 2009-04-17 15:22:33 UTC (rev 98)
@@ -27,7 +27,12 @@
// dt = size (duration) of the Euler step
// on output:
// x = contains the new state vector (i.e., at time t+dt)
+//
+// NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function.
+// The RNG is initialized before any call to this function, and the RNG state is written afterward.
+// Inclusion of these calls in the user-defined function may result in significant slowdown.
+
// Prototype for one-step Euler PDF, as used by "euler.density":
typedef void euler_step_pdf(double *f,
double *x1, double *x2, double t1, double t2, const double *p,
@@ -99,6 +104,10 @@
// t = time at the beginning of the Euler step
// on output:
// y = pointer to vector containing simulated observations (length = nobs = nrow(data))
+//
+// NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function.
+// The RNG is initialized before any call to this function, and the RNG state is written afterward.
+// Inclusion of these calls in the user-defined function may result in significant slowdown.
// Prototype for measurement model density evaluator
Modified: pkg/src/rmeasure.c
===================================================================
--- pkg/src/rmeasure.c 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/src/rmeasure.c 2009-04-17 15:22:33 UTC (rev 98)
@@ -239,10 +239,14 @@
ndim[0] = nvars; ndim[1] = npars; ndim[2] = nreps; ndim[3] = ntimes;
ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nobs;
+ if (use_native) GetRNGstate();
+
simul_meas(ff,REAL(Y),REAL(x),REAL(times),REAL(params),ndim,
sidx,pidx,cidx,
REAL(tcovar),REAL(covar));
+ if (use_native) PutRNGstate();
+
if (OIDX != 0) Free(OIDX);
OIDX = 0;
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/tests/sir.R 2009-04-17 15:22:33 UTC (rev 98)
@@ -192,7 +192,7 @@
covar=basis,
delta.t=1/52/20,
statenames=c("S","I","R","cases","W","B","dW"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
covarnames=c("seas1"),
zeronames=c("cases"),
step.fun="sir_euler_simulator",
@@ -200,8 +200,9 @@
dens.fun="sir_euler_density",
dprocess=euler.density,
skeleton.vectorfield="sir_ODE",
+ rmeasure="binom_rmeasure",
+ dmeasure="binom_dmeasure",
PACKAGE="pomp",
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
p <- exp(params)
with(
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-04-16 20:11:12 UTC (rev 97)
+++ pkg/tests/sir.Rout.save 2009-04-17 15:22:33 UTC (rev 98)
@@ -148,7 +148,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.595228 secs
+Time difference of 4.017606 secs
>
> pdf(file='sir.pdf')
>
@@ -165,7 +165,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 5.991744 secs
+Time difference of 9.691413 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -223,7 +223,7 @@
+ covar=basis,
+ delta.t=1/52/20,
+ statenames=c("S","I","R","cases","W","B","dW"),
-+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
++ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
+ covarnames=c("seas1"),
+ zeronames=c("cases"),
+ step.fun="sir_euler_simulator",
@@ -231,8 +231,9 @@
+ dens.fun="sir_euler_density",
+ dprocess=euler.density,
+ skeleton.vectorfield="sir_ODE",
++ rmeasure="binom_rmeasure",
++ dmeasure="binom_dmeasure",
+ PACKAGE="pomp",
-+ measurement.model=measles~binom(size=cases,prob=exp(rho)),
+ initializer=function(params,t0,...){
+ p <- exp(params)
+ with(
@@ -256,7 +257,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.2616198 secs
+Time difference of 0.2553811 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -264,7 +265,7 @@
> X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.824738 secs
+Time difference of 7.702386 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
More information about the pomp-commits
mailing list