[Pomp-commits] r1044 - in pkg/pompExamples: . inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 7 00:07:54 CET 2015
Author: kingaa
Date: 2015-01-07 00:07:54 +0100 (Wed, 07 Jan 2015)
New Revision: 1044
Removed:
pkg/pompExamples/src/ebola.c
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/inst/examples/ebola.R
pkg/pompExamples/tests/bbp.Rout.save
pkg/pompExamples/tests/budmoth.Rout.save
pkg/pompExamples/tests/ebola.Rout.save
pkg/pompExamples/tests/parus.Rout.save
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- pull Ebola example back to using on-the-fly compiled codes
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/DESCRIPTION 2015-01-06 23:07:54 UTC (rev 1044)
@@ -1,8 +1,8 @@
Package: pompExamples
Type: Package
Title: Additional pomp examples
-Version: 0.25-3
-Date: 2015-01-02
+Version: 0.25-4
+Date: 2015-01-05
Maintainer: Aaron A. King <kingaa at umich.edu>
Authors at R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),
email="kingaa at umich.edu"),
Modified: pkg/pompExamples/inst/examples/ebola.R
===================================================================
--- pkg/pompExamples/inst/examples/ebola.R 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/inst/examples/ebola.R 2015-01-06 23:07:54 UTC (rev 1044)
@@ -53,6 +53,117 @@
dat <- melt(dat,id="week",variable.name="country",value.name="cases")
mutate(dat,deaths=NA) -> dat
+
+paruntrans <- Csnippet('
+ double *IC = &S_0;
+ double *TIC = &TS_0;
+ TN = log(N);
+ TR0 = log(R0);
+ Talpha = log(alpha);
+ Tgamma = log(gamma);
+ Trho = logit(rho);
+ Tk = log(k);
+ Tcfr = logit(cfr);
+ to_log_barycentric(TIC,IC,4);
+')
+
+partrans <- Csnippet('
+ double *IC = &S_0;
+ double *TIC = &TS_0;
+ TN = exp(N);
+ TR0 = exp(R0);
+ Talpha = exp(alpha);
+ Tgamma = exp(gamma);
+ Trho = expit(rho);
+ Tk = exp(k);
+ Tcfr = expit(cfr);
+ from_log_barycentric(TIC,IC,4);
+')
+
+## Observation model: hierarchical model for cases and deaths
+## p(R_t, D_t| C_t) = p(R_t | C_t) * p(D_t | C_t, R_t)
+## p(R_t | C_t): Negative binomial with mean rho * C_t and dispersion parameter 1 / k
+## p(D_t | C_t, R_t): Binomial B(R_t, cfr)
+
+dObs <- Csnippet('
+ double f;
+ if (k > 0.0)
+ f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
+ else
+ f = dpois(nearbyint(cases),rho*N_EI,1);
+ lik = (give_log) ? f : exp(f);
+')
+
+dObsLS <- Csnippet('
+ double f;
+ f = dnorm(cases,rho*N_EI,k,1);
+ lik = (give_log) ? f : exp(f);
+')
+
+rObs <- Csnippet('
+ if (k > 0) {
+ cases = rnbinom_mu(1.0/k,rho*N_EI);
+ deaths = rnbinom_mu(1.0/k,rho*cfr*N_IR);
+ } else {
+ cases = rpois(rho*N_EI);
+ deaths = rpois(rho*cfr*N_IR);
+ }')
+
+rObsLS <- Csnippet('
+ cases = rnorm(rho*N_EI,k);
+ deaths = NA_REAL;
+')
+
+rSim <- Csnippet('
+ double lambda, beta;
+ double *E = &E1;
+ beta = R0 * gamma; // Transmission rate
+ lambda = beta * I / N; // Force of infection
+ int i;
+
+ // Transitions
+ // From class S
+ double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
+ // From class E
+ double transE[nstageE]; // No of transitions between classes E
+ for(i = 0; i < nstageE; i++){
+ transE[i] = rbinom(E[i], 1.0 - exp(- nstageE * alpha * dt));
+ }
+ // From class I
+ double transI = rbinom(I, 1.0 - exp(- gamma * dt)); // No of transitions I->R
+
+ // Balance the equations
+ S -= transS;
+ E[0] += transS - transE[0];
+ for(i=1; i < nstageE; i++) {
+ E[i] += transE[i-1] - transE[i];
+ }
+ I += transE[nstageE - 1] - transI;
+ R += transI;
+ N_EI += transE[nstageE - 1]; // No of transitions from E to I
+ N_IR += transI; // No of transitions from I to R
+')
+
+skel <- Csnippet('
+ double lambda, beta;
+ double *E = &E1;
+ double *DE = &DE1;
+ beta = R0 * gamma; // Transmission rate
+ lambda = beta * I / N; // Force of infection
+ int i;
+
+ // Balance the equations
+ DS = - lambda * S;
+ DE[0] = lambda * S - nstageE * alpha * E[0];
+ for (i=1; i < nstageE; i++)
+ DE[i] = nstageE * alpha * (E[i-1]-E[i]);
+ DI = nstageE * alpha * E[nstageE-1] - gamma * I;
+ DR = gamma * I;
+ DN_EI = nstageE * alpha * E[nstageE-1];
+ DN_IR = gamma * I;
+')
+
+
ebolaModel <- function (country=c("Guinea", "SierraLeone", "Liberia", "WestAfrica"),
data = NULL,
timestep = 0.01, nstageE = 3L,
@@ -69,7 +180,10 @@
infectious_period <- 7/7
index_case <- 10/pop
dt <- timestep
+ nstageE <- as.integer(nstageE)
+ globs <- paste0("static int nstageE = ",nstageE,";");
+
theta <- c(N=pop,R0=1.4,
alpha=-1/(nstageE*dt)*log(1-nstageE*dt/incubation_period),
gamma=-log(1-dt/infectious_period)/dt,
@@ -89,7 +203,7 @@
} else {
dat <- data
}
-
+
if (na.rm) {
dat <- mutate(subset(dat,!is.na(cases)),week=week-min(week)+1)
}
@@ -103,20 +217,20 @@
times="week",
t0=0,
params=theta,
+ globals=globs,
obsnames=c("cases","deaths"),
statenames=c("S","E1","I","R","N_EI","N_IR"),
zeronames=if (type=="raw") c("N_EI","N_IR") else character(0),
paramnames=c("N","R0","alpha","gamma","rho","k","cfr",
"S_0","E_0","I_0","R_0"),
- nstageE=as.integer(nstageE),
- PACKAGE="pompExamples",
- dmeasure=if (least.sq) "_ebola_dObsLS" else "_ebola_dObs",
- rmeasure=if (least.sq) "_ebola_rObsLS" else "_ebola_rObs",
- rprocess=discrete.time.sim(step.fun="_ebola_rSim",delta.t=timestep),
- skeleton="_ebola_skel",
+ nstageE=nstageE,
+ dmeasure=if (least.sq) dObsLS else dObs,
+ rmeasure=if (least.sq) rObsLS else rObs,
+ rprocess=discrete.time.sim(step.fun=rSim,delta.t=timestep),
+ skeleton=skel,
skeleton.type="vectorfield",
- parameter.transform="_ebola_par_trans",
- parameter.inv.transform="_ebola_par_untrans",
+ parameter.transform=partrans,
+ parameter.inv.transform=paruntrans,
initializer=function (params, t0, nstageE, ...) {
all.state.names <- c("S",paste0("E",1:nstageE),"I","R","N_EI","N_IR")
comp.names <- c("S",paste0("E",1:nstageE),"I","R")
Deleted: pkg/pompExamples/src/ebola.c
===================================================================
--- pkg/pompExamples/src/ebola.c 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/src/ebola.c 2015-01-06 23:07:54 UTC (rev 1044)
@@ -1,187 +0,0 @@
-// SEIR Ebola model
-
-#include <R.h>
-#include <Rmath.h>
-#include <R_ext/Rdynload.h>
-#include <pomp.h>
-
-// State variables
-#define S (x[stateindex[0]]) // Susceptible
-#define E(J) (x[stateindex[1] + (J)]) // Exposed
-#define I (x[stateindex[2]]) // Infected
-#define R (x[stateindex[3]]) // Removed
-#define N_EI (x[stateindex[4]]) // Number of transitions from E to I
-#define N_IR (x[stateindex[5]]) // Number of transitions from I to R
-
-// Variations
-#define DS (f[stateindex[0]]) // Susceptible
-#define DE(J) (f[stateindex[1] + (J)]) // Exposed
-#define DI (f[stateindex[2]]) // Infected
-#define DR (f[stateindex[3]]) // Removed
-#define DN_EI (f[stateindex[4]]) // Number of transitions from E to I
-#define DN_IR (f[stateindex[5]]) // Number of transitions from I to R
-
-// Parameters on the natural scale (all rates are per day)
-#define N (p[parindex[0]]) // Population size
-#define R0 (p[parindex[1]]) // Basic reproduction number
-#define alpha (p[parindex[2]]) // Inverse of latency period
-#define gamma (p[parindex[3]]) // Inverse of duration of infection
-#define rho (p[parindex[4]]) // Reporting probability
-#define k (p[parindex[5]]) // Reporting overdispersion
-#define cfr (p[parindex[6]]) // Case-fatality ratio
-#define IC(J) (p[parindex[7] + (J)]) // Initial conditions
-
-// Parameters on the transformed scale (all rates are per day)
-#define TN (pt[parindex[0]]) // Population size
-#define TR0 (pt[parindex[1]]) // Basic reproduction number
-#define Talpha (pt[parindex[2]]) // Inverse of latency period
-#define Tgamma (pt[parindex[3]]) // Inverse of duration of infection
-#define Trho (pt[parindex[4]]) // Reporting probability
-#define Tk (pt[parindex[5]]) // Reporting overdispersion
-#define Tcfr (pt[parindex[6]]) // Case-fatality ratio
-#define TIC(J) (pt[parindex[7] + (J)]) // Initial conditions
-
-// Observations
-#define cases (y[obsindex[0]]) // Number of reported cases
-#define deaths (y[obsindex[1]]) // Number of reported deaths
-
-// Transforms the parameters to the transformed scale
-void _ebola_par_untrans (double *pt, double *p, int *parindex){
- TN = log(N);
- TR0 = log(R0);
- Talpha = log(alpha);
- Tgamma = log(gamma);
- Trho = logit(rho);
- Tk = log(k);
- Tcfr = logit(cfr);
- to_log_barycentric(&(TIC(0)),&(IC(0)),4);
-}
-
-// Transforms the parameters to the natural scale
-void _ebola_par_trans (double *pt, double *p, int *parindex){
- TN = exp(N);
- TR0 = exp(R0);
- Talpha = exp(alpha);
- Tgamma = exp(gamma);
- Trho = expit(rho);
- Tk = exp(k);
- Tcfr = expit(cfr);
- from_log_barycentric(&(TIC(0)),&(IC(0)),4);
-}
-
-// Observation model: hierarchical model for cases and deaths
-// p(R_t, D_t| C_t) = p(R_t | C_t) * p(D_t | C_t, R_t)
-// p(R_t | C_t): Negative binomial with mean rho * C_t and dispersion parameter 1 / k
-// p(D_t | C_t, R_t): Binomial B(R_t, cfr)
-void _ebola_dObs (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
- double f;
- if (k > 0.0)
- f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
- // f += dnbinom_mu(nearbyint(deaths), 1.0 / k, rho * cfr * N_IR, 1);
- else
- f = dpois(nearbyint(cases),rho*N_EI,1);
- *lik = (give_log) ? f : exp(f);
-}
-
-// For least-squares trajectory-matching:
-void _ebola_dObsLS (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
- double f;
- f = dnorm(cases,rho*N_EI,k,1);
- *lik = (give_log) ? f : exp(f);
-}
-
-void _ebola_rObs (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t)
-{
- if (k > 0) {
- cases = rnbinom_mu(1.0/k,rho*N_EI);
- deaths = rnbinom_mu(1.0/k,rho*cfr*N_IR);
- } else {
- cases = rpois(rho*N_EI);
- deaths = rpois(rho*cfr*N_IR);
- }
-}
-
-// For least-squares trajectory-matching:
-void _ebola_rObsLS (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t)
-{
- cases = rnorm(rho*N_EI,k);
- deaths = NA_REAL;
-}
-
-
-// Process model
-void _ebola_rSim (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covars,
- double t, double dt) {
-
- // Retrieve user data in the pomp object
- int *(*get_pomp_userdata_int)(const char *);
- get_pomp_userdata_int = (int *(*)(const char *)) R_GetCCallable("pomp","get_pomp_userdata_int");
- int nstageE = *(get_pomp_userdata_int("nstageE")); // Number of stages in the E class
-
- // Other parameters
- double lambda, beta;
- beta = R0 * gamma; // Transmission rate
- lambda = beta * I / N; // Force of infection
- int i;
-
- // Transitions
-
- // From class S
- double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
-
- // From class E
- double transE[nstageE]; // No of transitions between classes E
- for(i = 0; i < nstageE; i++){
- transE[i] = rbinom(E(i), 1.0 - exp(- nstageE * alpha * dt));
- }
-
- // From class I
- double transI = rbinom(I, 1.0 - exp(- gamma * dt)); // No of transitions I->R
-
- // Balance the equations
- S -= transS;
- E(0) += transS - transE[0];
- for(i=1; i < nstageE; i++) {
- E(i) += transE[i-1] - transE[i];
- }
- I += transE[nstageE - 1] - transI;
- R += transI;
- N_EI += transE[nstageE - 1]; // No of transitions from E to I
- N_IR += transI; // No of transitions from I to R
-}
-
-// Continuous-time deterministic skeleton
-void _ebola_skel (double *f, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
-
- // Retrieve user data in the pomp object
- int *(*get_pomp_userdata_int)(const char *);
- get_pomp_userdata_int = (int *(*)(const char *)) R_GetCCallable("pomp","get_pomp_userdata_int");
- int nstageE = *(get_pomp_userdata_int("nstageE")); // Number of stages in the E class
-
- // Other parameters
- double lambda, beta;
- beta = R0 * gamma; // Transmission rate
- lambda = beta * I / N; // Force of infection
- int i;
-
- // Balance the equations
- DS = - lambda * S;
- DE(0) = lambda * S - nstageE * alpha * E(0);
- for(i=1; i < nstageE; i++) DE(i) = nstageE * alpha * (E(i - 1) - E(i));
- DI = nstageE * alpha * E(nstageE - 1) - gamma * I;
- DR = gamma * I;
- DN_EI = nstageE * alpha * E(nstageE - 1);
- DN_IR = gamma * I;
-}
Modified: pkg/pompExamples/tests/bbp.Rout.save
===================================================================
--- pkg/pompExamples/tests/bbp.Rout.save 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/tests/bbp.Rout.save 2015-01-06 23:07:54 UTC (rev 1044)
@@ -23,11 +23,11 @@
> set.seed(47575684L)
>
> pompExample(bbp)
-newly created pomp object(s):
+newly created object(s):
bbp
> pf <- pfilter(simulate(bbp),Np=100,max.fail=Inf)
> tj <- trajectory(bbp)
>
> proc.time()
user system elapsed
- 0.573 0.038 0.597
+ 0.728 0.028 0.786
Modified: pkg/pompExamples/tests/budmoth.Rout.save
===================================================================
--- pkg/pompExamples/tests/budmoth.Rout.save 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/tests/budmoth.Rout.save 2015-01-06 23:07:54 UTC (rev 1044)
@@ -127,4 +127,4 @@
>
> proc.time()
user system elapsed
- 0.614 0.056 0.652
+ 0.732 0.060 0.818
Modified: pkg/pompExamples/tests/ebola.Rout.save
===================================================================
--- pkg/pompExamples/tests/ebola.Rout.save 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/tests/ebola.Rout.save 2015-01-06 23:07:54 UTC (rev 1044)
@@ -25,19 +25,43 @@
> pompExample(ebola)
Loading required package: plyr
Loading required package: reshape2
-newly created pomp object(s):
+newly created object(s):
ebolaModel
> ebolaModel(country="Guinea") -> po
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+gcc -std=gnu99 -I/usr/local/apps/R/R-3.1.2/lib64/R/include -DNDEBUG -I/usr/local/include -I/usr/local/apps/R/R-3.1.2/lib64/R/library/pomp/include -fpic -g -O2 -c /tmp/RtmpU0NeMt/pomp2583FF7C2A2A.c -o /tmp/RtmpU0NeMt/pomp2583FF7C2A2A.o
+gcc -std=gnu99 -shared -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -L/usr/local/lib64 -o /tmp/RtmpU0NeMt/pomp2583FF7C2A2A.so /tmp/RtmpU0NeMt/pomp2583FF7C2A2A.o -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -lR
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+model codes written to '/tmp/RtmpU0NeMt/pomp2583FF7C2A2A.c'
+link to shared-object library '/tmp/RtmpU0NeMt/pomp2583FF7C2A2A.so'
> pf <- pfilter(simulate(po),Np=100)
> tj <- trajectory(po)
>
> ebolaModel(country="SierraLeone",na.rm=TRUE,type='cum') -> po
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+gcc -std=gnu99 -I/usr/local/apps/R/R-3.1.2/lib64/R/include -DNDEBUG -I/usr/local/include -I/usr/local/apps/R/R-3.1.2/lib64/R/library/pomp/include -fpic -g -O2 -c /tmp/RtmpU0NeMt/pompDA1B7D4C88B1.c -o /tmp/RtmpU0NeMt/pompDA1B7D4C88B1.o
+gcc -std=gnu99 -shared -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -L/usr/local/lib64 -o /tmp/RtmpU0NeMt/pompDA1B7D4C88B1.so /tmp/RtmpU0NeMt/pompDA1B7D4C88B1.o -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -lR
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+model codes written to '/tmp/RtmpU0NeMt/pompDA1B7D4C88B1.c'
+link to shared-object library '/tmp/RtmpU0NeMt/pompDA1B7D4C88B1.so'
> pf <- pfilter(simulate(po),Np=100)
> tj <- trajectory(po)
> dd <- simulate(po,as.data.frame=TRUE,obs=TRUE)
> dd$week <- dd$time
> po <- ebolaModel(data=subset(dd,select=c(week,cases,deaths)))
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+gcc -std=gnu99 -I/usr/local/apps/R/R-3.1.2/lib64/R/include -DNDEBUG -I/usr/local/include -I/usr/local/apps/R/R-3.1.2/lib64/R/library/pomp/include -fpic -g -O2 -c /tmp/RtmpU0NeMt/pompBF9800D29D72.c -o /tmp/RtmpU0NeMt/pompBF9800D29D72.o
+gcc -std=gnu99 -shared -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -L/usr/local/lib64 -o /tmp/RtmpU0NeMt/pompBF9800D29D72.so /tmp/RtmpU0NeMt/pompBF9800D29D72.o -L/usr/local/apps/R/R-3.1.2/lib64/R/lib -lR
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Entering directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+make[1]: Leaving directory `/userdata/kingaa/projects/pomp/pkg/pompExamples.Rcheck/tests'
+model codes written to '/tmp/RtmpU0NeMt/pompBF9800D29D72.c'
+link to shared-object library '/tmp/RtmpU0NeMt/pompBF9800D29D72.so'
>
> proc.time()
user system elapsed
- 1.100 0.034 1.119
+ 2.032 0.228 3.798
Modified: pkg/pompExamples/tests/parus.Rout.save
===================================================================
--- pkg/pompExamples/tests/parus.Rout.save 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/tests/parus.Rout.save 2015-01-06 23:07:54 UTC (rev 1044)
@@ -41,4 +41,4 @@
>
> proc.time()
user system elapsed
- 0.572 0.037 0.596
+ 0.716 0.060 0.806
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2015-01-02 14:04:49 UTC (rev 1043)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2015-01-06 23:07:54 UTC (rev 1044)
@@ -168,4 +168,4 @@
>
> proc.time()
user system elapsed
- 23.697 0.033 23.738
+ 18.417 0.068 18.567
More information about the pomp-commits
mailing list