[Pomp-commits] r66 - pkg/inst/examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 27 15:50:10 CET 2009
Author: kingaa
Date: 2009-01-27 15:50:10 +0100 (Tue, 27 Jan 2009)
New Revision: 66
Added:
pkg/inst/examples/euler_sir.R
pkg/inst/examples/euler_sir.c
Removed:
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
Log:
moved 'sir.*' examples to 'euler_sir.*'
Copied: pkg/inst/examples/euler_sir.R (from rev 65, pkg/inst/examples/sir.R)
===================================================================
--- pkg/inst/examples/euler_sir.R (rev 0)
+++ pkg/inst/examples/euler_sir.R 2009-01-27 14:50:10 UTC (rev 66)
@@ -0,0 +1,190 @@
+require(pomp)
+
+## basis functions for the seasonality
+tbasis <- seq(0,25,by=1/52)
+basis <- periodic.bspline.basis(tbasis,nbasis=3)
+colnames(basis) <- paste("seas",1:3,sep='')
+
+## some parameters
+params <- c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
+ )
+
+## set up the pomp object
+po <- pomp(
+ times=1/52*seq.int(length=4*52),
+ data=rbind(measles=numeric(52*4)),
+ t0=0,
+ tcovar=tbasis,
+ covar=basis,
+ delta.t=1/52/20,
+ zeronames=c("cases"),
+ step.fun=function(t,x,params,covars,delta.t,...) {
+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
+ {
+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ beta.var <- beta.sd^2
+ dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
+ foi <- (iota+beta*I*dW/delta.t)/pop
+ trans <- c(
+ rpois(n=1,lambda=mu*pop*delta.t),
+ reulermultinom(n=1,size=S,rate=c(foi,mu),dt=delta.t),
+ reulermultinom(n=1,size=I,rate=c(gamma,mu),dt=delta.t),
+ reulermultinom(n=1,size=R,rate=c(mu),dt=delta.t)
+ )
+ c(
+ S=S+trans[1]-trans[2]-trans[3],
+ I=I+trans[2]-trans[4]-trans[5],
+ R=R+trans[4]-trans[6],
+ cases=cases+trans[4],
+ W=if (beta.sd>0) W+(dW-delta.t)/beta.sd else W,
+ B=trans[1],
+ SI=trans[2],
+ SD=trans[3],
+ IR=trans[4],
+ ID=trans[5],
+ RD=trans[6],
+ dW=dW
+ )
+ }
+ )
+ },
+ dens.fun=function(t1,t2,params,x1,x2,covars,...) {
+ params <- exp(params)
+ with(
+ as.list(params),
+ {
+ dt <- t2-t1
+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ beta.var <- beta.sd^2
+ dW <- x2['dW']
+ foi <- (iota+beta*x1["I"]*dW/dt)/pop
+ probs <- c(
+ dpois(x=x2["B"],lambda=mu*pop*dt,log=T),
+ deulermultinom(x=x2[c("SI","SD")],size=x1["S"],rate=c(foi,mu),dt=dt,log=T),
+ deulermultinom(x=x2[c("IR","ID")],size=x1["I"],rate=c(gamma,mu),dt=dt,log=T),
+ deulermultinom(x=x2["RD"],size=x1["R"],rate=c(mu),dt=dt,log=T),
+ dgamma(x=dW,shape=dt/beta.var,scale=beta.var,log=T)
+ )
+ sum(probs)
+ }
+ )
+ },
+ skeleton.vectorfield=function(x,t,params,covars,...) {
+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
+ {
+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ foi <- (iota+beta*I)/pop
+ terms <- c(
+ mu*pop,
+ foi*S,
+ mu*S,
+ gamma*I,
+ mu*I,
+ mu*R
+ )
+ xdot <- c(
+ terms[1]-terms[2]-terms[3],
+ terms[2]-terms[4]-terms[5],
+ terms[4]-terms[6],
+ terms[4]
+ )
+ ifelse(x>0,xdot,0)
+ }
+ )
+ },
+ rprocess=euler.simulate,
+ dprocess=euler.density,
+ measurement.model=measles~binom(size=cases,prob=exp(rho)),
+ initializer=function(params,t0,...){
+ p <- exp(params)
+ with(
+ as.list(p),
+ {
+ fracs <- c(S.0,I.0,R.0)
+ x0 <- c(
+ round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
+ rep(0,9) # zeros for 'cases', 'W', and the transition numbers
+ )
+ names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
+ x0
+ }
+ )
+ }
+ )
+
+## alternatively, one can define the computationally intensive bits using native routines:
+## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory (file "euler_sir.c")
+
+if (Sys.info()['sysname']=='Linux') {
+
+ modelfile <- system.file("examples/euler_sir.c",package="pomp")
+ includedir <- system.file("include",package="pomp")
+ lib <- system.file("libs/pomp.so",package="pomp")
+
+ ## compile the model into a shared-object library
+ system(paste("cp",modelfile,"."))
+ system(paste("cp ",includedir,"/pomp.h .",sep=""))
+ system(paste("R CMD SHLIB -o ./euler_sir.so euler_sir.c",lib))
+
+ po <- pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
+ t0=0,
+ tcovar=tbasis,
+ 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"),
+ covarnames=c("seas1"),
+ zeronames=c("cases"),
+ step.fun="sir_euler_simulator",
+ rprocess=euler.simulate,
+ dens.fun="sir_euler_density",
+ dprocess=euler.density,
+ skeleton.vectorfield="sir_ODE",
+ 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(
+ as.list(p),
+ {
+ fracs <- c(S.0,I.0,R.0)
+ x0 <- c(
+ round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
+ rep(0,9) # zeros for 'cases', 'W', and the transition numbers
+ )
+ names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
+ x0
+ }
+ )
+ }
+ )
+
+ dyn.load("euler_sir.so") ## load the shared-object library
+
+ ## compute a trajectory of the deterministic skeleton
+ tic <- Sys.time()
+ X <- trajectory(po,params=log(params),hmax=1/52)
+ toc <- Sys.time()
+ print(toc-tic)
+
+ ## simulate from the model
+ tic <- Sys.time()
+ x <- simulate(po,params=log(params),nsim=3)
+ toc <- Sys.time()
+ print(toc-tic)
+
+ dyn.unload("euler_sir.so")
+
+}
Copied: pkg/inst/examples/euler_sir.c (from rev 65, pkg/inst/examples/sir.c)
===================================================================
--- pkg/inst/examples/euler_sir.c (rev 0)
+++ pkg/inst/examples/euler_sir.c 2009-01-27 14:50:10 UTC (rev 66)
@@ -0,0 +1,265 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+static double expit (double x) {
+ return 1.0/(1.0 + exp(-x));
+}
+
+static double logit (double x) {
+ return log(x/(1-x));
+}
+
+static double term_time (double t, double b0, double b1)
+{
+ static double correction = 0.52328767123287671233;
+ double day = 365.0 * (t - floor(t));
+ int tt = (day >= 7.0 && day <= 100.0)
+ || (day >= 115.0 && day <= 200.0)
+ || (day >= 252.0 && day <= 300.0)
+ || (day >= 308.0 && day <= 356.0);
+ return b0 * (1.0 + b1 * (-1.0 + 2.0 * ((double) tt))) / (1.0 + correction * b1);
+}
+
+#define LOGGAMMA (p[parindex[0]]) // recovery rate
+#define LOGMU (p[parindex[1]]) // baseline birth and death rate
+#define LOGIOTA (p[parindex[2]]) // import rate
+#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 SUSC (x[stateindex[0]]) // number of susceptibles
+#define INFD (x[stateindex[1]]) // number of infectives
+#define RCVD (x[stateindex[2]]) // number of recovereds
+#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x[stateindex[4]]) // integrated white noise
+#define BIRTHS (x[stateindex[5]]) // births
+#define dW (x[stateindex[6]]) // white noise process
+
+#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
+
+// SIR model with Euler multinomial step
+// forced transmission (basis functions passed as covariates)
+// constant population size as a parameter
+// environmental stochasticity on transmission
+void sir_euler_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double *trans; // transition numbers
+ double gamma, mu, iota, beta_sd, beta_var, popsize;
+ double beta;
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ beta_sd = exp(LOGBETA_SD);
+ popsize = exp(LOGPOPSIZE);
+ beta_var = beta_sd*beta_sd;
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+
+ // test to make sure the parameters and state variable values are sane
+ if (!(R_FINITE(beta)) ||
+ !(R_FINITE(gamma)) ||
+ !(R_FINITE(mu)) ||
+ !(R_FINITE(beta_sd)) ||
+ !(R_FINITE(iota)) ||
+ !(R_FINITE(popsize)) ||
+ !(R_FINITE(SUSC)) ||
+ !(R_FINITE(INFD)) ||
+ !(R_FINITE(RCVD)) ||
+ !(R_FINITE(CASE)) ||
+ !(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;
+ } else { // environmental noise is OFF
+ dW = dt;
+ }
+
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the transition numbers
+ trans = &BIRTHS;
+ trans[0] = rpois(rate[0]*dt); // births are Poisson
+ reulermultinom(2,SUSC,&rate[1],dt,&trans[1]);
+ 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];
+ RCVD += trans[3]-trans[5];
+ CASE += trans[3]; // cases are cumulative recoveries
+ if (beta_sd > 0.0) {
+ W += (dW-dt)/beta_sd; // mean zero, variance = dt
+ }
+
+}
+
+#undef BIRTHS
+#undef W
+#undef dW
+
+#define DSDT (f[stateindex[0]])
+#define DIDT (f[stateindex[1]])
+#define DRDT (f[stateindex[2]])
+#define DCDT (f[stateindex[3]])
+
+void sir_ODE (double *f, double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar, double t)
+{
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double term[nrate]; // terms in the equations
+ double gamma, mu, iota, popsize;
+ double beta;
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ popsize = exp(LOGPOPSIZE);
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*INFD)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the several terms
+ term[0] = rate[0];
+ term[1] = rate[1]*SUSC;
+ term[2] = rate[2]*SUSC;
+ term[3] = rate[3]*INFD;
+ term[4] = rate[4]*INFD;
+ term[5] = rate[5]*RCVD;
+
+ // balance the equations
+ DSDT = term[0]-term[1]-term[2];
+ DIDT = term[1]-term[3]-term[4];
+ DRDT = term[3]-term[5];
+ DCDT = term[3]; // cases are cumulative recoveries
+
+}
+
+#undef DSDT
+#undef DIDT
+#undef DRDT
+#undef DCDT
+
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+
+#define SUSC (x1[stateindex[0]]) // number of susceptibles
+#define INFD (x1[stateindex[1]]) // number of infectives
+#define RCVD (x1[stateindex[2]]) // number of recovereds
+#define CASE (x1[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x1[stateindex[4]]) // integrated white noise
+#define BIRTHS (x2[stateindex[5]]) // births
+#define dW (x2[stateindex[6]]) // white noise process
+
+// SIR model with Euler multinomial step
+// forced transmission (basis functions passed as covariates)
+// constant population size as a parameter
+// environmental stochasticity on transmission
+void sir_euler_density (double *f, double *x1, double *x2, double t1, double t2, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar)
+{
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double *trans; // transition numbers
+ double gamma, mu, iota, beta_sd, popsize;
+ double beta;
+ double dt = t2-t1;
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ beta_sd = exp(LOGBETA_SD);
+ popsize = exp(LOGPOPSIZE);
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+
+ // test to make sure the parameters and state variable values are sane
+ if (!(R_FINITE(beta)) ||
+ !(R_FINITE(gamma)) ||
+ !(R_FINITE(mu)) ||
+ !(R_FINITE(beta_sd)) ||
+ !(R_FINITE(iota)) ||
+ !(R_FINITE(popsize)) ||
+ !(R_FINITE(SUSC)) ||
+ !(R_FINITE(INFD)) ||
+ !(R_FINITE(RCVD)) ||
+ !(R_FINITE(CASE)) ||
+ !(R_FINITE(W))) {
+ *f = R_NaN;
+ return;
+ }
+
+ // compute the transition rates
+ trans = &BIRTHS;
+ if (beta_sd > 0.0) { // environmental noise is ON
+ double beta_var = beta_sd*beta_sd;
+ *f = dgamma(dW,dt/beta_var,beta_var,1);
+ } else { // environmental noise is OFF
+ *f = 0; // THIS ASSUMES THAT dw = dt !!!
+ }
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the transition numbers
+ *f += dpois(trans[0],rate[0]*dt,1); // births are Poisson
+ *f += deulermultinom(2,SUSC,&rate[1],dt,&trans[1],1);
+ *f += deulermultinom(2,INFD,&rate[3],dt,&trans[3],1);
+ *f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
+}
+
+#undef GAMMA
+#undef MU
+#undef IOTA
+#undef BETA
+#undef BETA_SD
+#undef POPSIZE
+
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
+#undef BIRTHS
+#undef DW
+
+#undef SEASBASIS
Deleted: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2009-01-27 14:35:56 UTC (rev 65)
+++ pkg/inst/examples/sir.R 2009-01-27 14:50:10 UTC (rev 66)
@@ -1,190 +0,0 @@
-require(pomp)
-
-## basis functions for the seasonality
-tbasis <- seq(0,25,by=1/52)
-basis <- periodic.bspline.basis(tbasis,nbasis=3)
-colnames(basis) <- paste("seas",1:3,sep='')
-
-## some parameters
-params <- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
-
-## set up the pomp object
-po <- pomp(
- times=1/52*seq.int(length=4*52),
- data=rbind(measles=numeric(52*4)),
- t0=0,
- tcovar=tbasis,
- covar=basis,
- delta.t=1/52/20,
- zeronames=c("cases"),
- step.fun=function(t,x,params,covars,delta.t,...) {
- params <- exp(params)
- with(
- as.list(c(x,params)),
- {
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
- beta.var <- beta.sd^2
- dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
- foi <- (iota+beta*I*dW/delta.t)/pop
- trans <- c(
- rpois(n=1,lambda=mu*pop*delta.t),
- reulermultinom(n=1,size=S,rate=c(foi,mu),dt=delta.t),
- reulermultinom(n=1,size=I,rate=c(gamma,mu),dt=delta.t),
- reulermultinom(n=1,size=R,rate=c(mu),dt=delta.t)
- )
- c(
- S=S+trans[1]-trans[2]-trans[3],
- I=I+trans[2]-trans[4]-trans[5],
- R=R+trans[4]-trans[6],
- cases=cases+trans[4],
- W=if (beta.sd>0) W+(dW-delta.t)/beta.sd else W,
- B=trans[1],
- SI=trans[2],
- SD=trans[3],
- IR=trans[4],
- ID=trans[5],
- RD=trans[6],
- dW=dW
- )
- }
- )
- },
- dens.fun=function(t1,t2,params,x1,x2,covars,...) {
- params <- exp(params)
- with(
- as.list(params),
- {
- dt <- t2-t1
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
- beta.var <- beta.sd^2
- dW <- x2['dW']
- foi <- (iota+beta*x1["I"]*dW/dt)/pop
- probs <- c(
- dpois(x=x2["B"],lambda=mu*pop*dt,log=T),
- deulermultinom(x=x2[c("SI","SD")],size=x1["S"],rate=c(foi,mu),dt=dt,log=T),
- deulermultinom(x=x2[c("IR","ID")],size=x1["I"],rate=c(gamma,mu),dt=dt,log=T),
- deulermultinom(x=x2["RD"],size=x1["R"],rate=c(mu),dt=dt,log=T),
- dgamma(x=dW,shape=dt/beta.var,scale=beta.var,log=T)
- )
- sum(probs)
- }
- )
- },
- skeleton.vectorfield=function(x,t,params,covars,...) {
- params <- exp(params)
- with(
- as.list(c(x,params)),
- {
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
- foi <- (iota+beta*I)/pop
- terms <- c(
- mu*pop,
- foi*S,
- mu*S,
- gamma*I,
- mu*I,
- mu*R
- )
- xdot <- c(
- terms[1]-terms[2]-terms[3],
- terms[2]-terms[4]-terms[5],
- terms[4]-terms[6],
- terms[4]
- )
- ifelse(x>0,xdot,0)
- }
- )
- },
- rprocess=euler.simulate,
- dprocess=euler.density,
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
- initializer=function(params,t0,...){
- p <- exp(params)
- with(
- as.list(p),
- {
- fracs <- c(S.0,I.0,R.0)
- x0 <- c(
- round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
- rep(0,9) # zeros for 'cases', 'W', and the transition numbers
- )
- names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
- x0
- }
- )
- }
- )
-
-## alternatively, one can define the computationally intensive bits using native routines:
-## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory (file "sir.c")
-
-if (Sys.info()['sysname']=='Linux') {
-
- modelfile <- system.file("examples/sir.c",package="pomp")
- includedir <- system.file("include",package="pomp")
- lib <- system.file("libs/pomp.so",package="pomp")
-
- ## compile the model into a shared-object library
- system(paste("cp",modelfile,"."))
- system(paste("cp ",includedir,"/pomp.h .",sep=""))
- system(paste("R CMD SHLIB -o ./sir_example.so sir.c",lib))
-
- po <- pomp(
- times=seq(1/52,4,by=1/52),
- data=rbind(measles=numeric(52*4)),
- t0=0,
- tcovar=tbasis,
- 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"),
- covarnames=c("seas1"),
- zeronames=c("cases"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
- dens.fun="sir_euler_density",
- dprocess=euler.density,
- skeleton.vectorfield="sir_ODE",
- PACKAGE="sir_example", ## name of the shared-object library
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
- initializer=function(params,t0,...){
- p <- exp(params)
- with(
- as.list(p),
- {
- fracs <- c(S.0,I.0,R.0)
- x0 <- c(
- round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
- rep(0,9) # zeros for 'cases', 'W', and the transition numbers
- )
- names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
- x0
- }
- )
- }
- )
-
- dyn.load("sir_example.so") ## load the shared-object library
-
- ## compute a trajectory of the deterministic skeleton
- tic <- Sys.time()
- X <- trajectory(po,params=log(params),hmax=1/52)
- toc <- Sys.time()
- print(toc-tic)
-
- ## simulate from the model
- tic <- Sys.time()
- x <- simulate(po,params=log(params),nsim=3)
- toc <- Sys.time()
- print(toc-tic)
-
- dyn.unload("sir_example.so")
-
-}
Deleted: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2009-01-27 14:35:56 UTC (rev 65)
+++ pkg/inst/examples/sir.c 2009-01-27 14:50:10 UTC (rev 66)
@@ -1,265 +0,0 @@
-// dear emacs, please treat this as -*- C++ -*-
-
-#include <Rmath.h>
-
-#include "pomp.h"
-
-static double expit (double x) {
- return 1.0/(1.0 + exp(-x));
-}
-
-static double logit (double x) {
- return log(x/(1-x));
-}
-
-static double term_time (double t, double b0, double b1)
-{
- static double correction = 0.52328767123287671233;
- double day = 365.0 * (t - floor(t));
- int tt = (day >= 7.0 && day <= 100.0)
- || (day >= 115.0 && day <= 200.0)
- || (day >= 252.0 && day <= 300.0)
- || (day >= 308.0 && day <= 356.0);
- return b0 * (1.0 + b1 * (-1.0 + 2.0 * ((double) tt))) / (1.0 + correction * b1);
-}
-
-#define LOGGAMMA (p[parindex[0]]) // recovery rate
-#define LOGMU (p[parindex[1]]) // baseline birth and death rate
-#define LOGIOTA (p[parindex[2]]) // import rate
-#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 SUSC (x[stateindex[0]]) // number of susceptibles
-#define INFD (x[stateindex[1]]) // number of infectives
-#define RCVD (x[stateindex[2]]) // number of recovereds
-#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
-#define W (x[stateindex[4]]) // integrated white noise
-#define BIRTHS (x[stateindex[5]]) // births
-#define dW (x[stateindex[6]]) // white noise process
-
-#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
-
-// SIR model with Euler multinomial step
-// forced transmission (basis functions passed as covariates)
-// constant population size as a parameter
-// environmental stochasticity on transmission
-void sir_euler_simulator (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar,
- double t, double dt)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double *trans; // transition numbers
- double gamma, mu, iota, beta_sd, beta_var, popsize;
- double beta;
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
- beta_var = beta_sd*beta_sd;
-
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
-
- // test to make sure the parameters and state variable values are sane
- if (!(R_FINITE(beta)) ||
- !(R_FINITE(gamma)) ||
- !(R_FINITE(mu)) ||
- !(R_FINITE(beta_sd)) ||
- !(R_FINITE(iota)) ||
- !(R_FINITE(popsize)) ||
- !(R_FINITE(SUSC)) ||
- !(R_FINITE(INFD)) ||
- !(R_FINITE(RCVD)) ||
- !(R_FINITE(CASE)) ||
- !(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;
- } else { // environmental noise is OFF
- dW = dt;
- }
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
-
- // compute the transition numbers
- trans = &BIRTHS;
- trans[0] = rpois(rate[0]*dt); // births are Poisson
- reulermultinom(2,SUSC,&rate[1],dt,&trans[1]);
- 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];
- RCVD += trans[3]-trans[5];
- CASE += trans[3]; // cases are cumulative recoveries
- if (beta_sd > 0.0) {
- W += (dW-dt)/beta_sd; // mean zero, variance = dt
- }
-
-}
-
-#undef BIRTHS
-#undef W
-#undef dW
-
-#define DSDT (f[stateindex[0]])
-#define DIDT (f[stateindex[1]])
-#define DRDT (f[stateindex[2]])
-#define DCDT (f[stateindex[3]])
-
-void sir_ODE (double *f, double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar, double t)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double term[nrate]; // terms in the equations
- double gamma, mu, iota, popsize;
- double beta;
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- popsize = exp(LOGPOPSIZE);
-
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD)/popsize; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
-
- // compute the several terms
- term[0] = rate[0];
- term[1] = rate[1]*SUSC;
- term[2] = rate[2]*SUSC;
- term[3] = rate[3]*INFD;
- term[4] = rate[4]*INFD;
- term[5] = rate[5]*RCVD;
-
- // balance the equations
- DSDT = term[0]-term[1]-term[2];
- DIDT = term[1]-term[3]-term[4];
- DRDT = term[3]-term[5];
- DCDT = term[3]; // cases are cumulative recoveries
-
-}
-
-#undef DSDT
-#undef DIDT
-#undef DRDT
-#undef DCDT
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-
-#define SUSC (x1[stateindex[0]]) // number of susceptibles
-#define INFD (x1[stateindex[1]]) // number of infectives
-#define RCVD (x1[stateindex[2]]) // number of recovereds
-#define CASE (x1[stateindex[3]]) // number of cases (accumulated per reporting period)
-#define W (x1[stateindex[4]]) // integrated white noise
-#define BIRTHS (x2[stateindex[5]]) // births
-#define dW (x2[stateindex[6]]) // white noise process
-
-// SIR model with Euler multinomial step
-// forced transmission (basis functions passed as covariates)
-// constant population size as a parameter
-// environmental stochasticity on transmission
-void sir_euler_density (double *f, double *x1, double *x2, double t1, double t2, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double *trans; // transition numbers
- double gamma, mu, iota, beta_sd, popsize;
- double beta;
- double dt = t2-t1;
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
-
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
-
- // test to make sure the parameters and state variable values are sane
- if (!(R_FINITE(beta)) ||
- !(R_FINITE(gamma)) ||
- !(R_FINITE(mu)) ||
- !(R_FINITE(beta_sd)) ||
- !(R_FINITE(iota)) ||
- !(R_FINITE(popsize)) ||
- !(R_FINITE(SUSC)) ||
- !(R_FINITE(INFD)) ||
- !(R_FINITE(RCVD)) ||
- !(R_FINITE(CASE)) ||
- !(R_FINITE(W))) {
- *f = R_NaN;
- return;
- }
-
- // compute the transition rates
- trans = &BIRTHS;
- if (beta_sd > 0.0) { // environmental noise is ON
- double beta_var = beta_sd*beta_sd;
- *f = dgamma(dW,dt/beta_var,beta_var,1);
- } else { // environmental noise is OFF
- *f = 0; // THIS ASSUMES THAT dw = dt !!!
- }
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
-
- // compute the transition numbers
- *f += dpois(trans[0],rate[0]*dt,1); // births are Poisson
- *f += deulermultinom(2,SUSC,&rate[1],dt,&trans[1],1);
- *f += deulermultinom(2,INFD,&rate[3],dt,&trans[3],1);
- *f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
-}
-
-#undef GAMMA
-#undef MU
-#undef IOTA
-#undef BETA
-#undef BETA_SD
-#undef POPSIZE
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
-#undef BIRTHS
-#undef DW
-
-#undef SEASBASIS
More information about the pomp-commits
mailing list