[Pomp-commits] r308 - pkg/inst/examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 26 22:32:11 CEST 2010
Author: kingaa
Date: 2010-08-26 22:32:11 +0200 (Thu, 26 Aug 2010)
New Revision: 308
Modified:
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
Log:
- simplify the SIR example using compiled code
Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2010-08-26 20:06:50 UTC (rev 307)
+++ pkg/inst/examples/sir.R 2010-08-26 20:32:11 UTC (rev 308)
@@ -19,14 +19,12 @@
model <- "sir"
pkg <- "pomp"
modelfile <- paste(model,".c",sep="")
- cppflags <- paste("PKG_CPPFLAGS=-I",system.file("include",package=pkg),sep="")
- pkglib <- system.file(paste("libs/",pkg,.Platform$dynlib.ext,sep=""),package=pkg)
solib <- paste(model,.Platform$dynlib.ext,sep="")
## compile the model into a shared-object library
if (!file.copy(from=system.file(paste("examples/",modelfile,sep=""),package=pkg),to=getwd()))
stop("cannot copy source code ",modelfile," to ",getwd())
- cmd <- paste(cppflags,R.home("bin/R"),"CMD SHLIB -o",solib,modelfile,pkglib)
+ cmd <- paste(R.home("bin/R"),"CMD SHLIB -o",solib,modelfile)
rv <- system(cmd)
if (rv!=0)
stop("cannot compile shared-object library ",solib)
@@ -39,12 +37,12 @@
times="time",
t0=0,
rprocess=euler.sim(
- step.fun="_sir_euler_simulator", # native routine for the simulation step
+ step.fun="sir_euler_simulator", # native routine for the simulation step
delta.t=1/52/20
),
- skeleton.vectorfield="_sir_ODE", # native routine for the skeleton
- rmeasure="_sir_binom_rmeasure", # binomial measurement model
- dmeasure="_sir_binom_dmeasure", # binomial measurement model
+ skeleton.vectorfield="sir_ODE", # native routine for the skeleton
+ rmeasure="binomial_rmeasure", # binomial measurement model
+ dmeasure="binomial_dmeasure", # binomial measurement model
PACKAGE="sir", ## name of the shared-object library
## the order of the observables assumed in the native routines:
obsnames="reports",
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2010-08-26 20:06:50 UTC (rev 307)
+++ pkg/inst/examples/sir.c 2010-08-26 20:32:11 UTC (rev 308)
@@ -1,9 +1,9 @@
// dear emacs, please treat this as -*- C++ -*-
+#include <R.h>
#include <Rmath.h>
+#include <R_ext/Rdynload.h>
-#include "pomp.h"
-
static double expit (double x) {
return 1.0/(1.0 + exp(-x));
}
@@ -12,18 +12,6 @@
return log(x/(1-x));
}
-static double term_time (double t, double b0, double b1)
-{
- static double correction = 0.4958904;
- double day = 365.0 * (t - floor(t));
- double interm;
- interm = ((day >= 7.0 && day < 100.0)
- || (day >= 116.0 && day < 200.0)
- || (day >= 252.0 && day < 300.0)
- || (day >= 308.0 && day < 356.0)) ? 1.0 : -1.0;
- return b0*(1.0+b1*interm)/(1.0+b1*correction);
-}
-
#define LOGGAMMA (p[parindex[0]]) // recovery rate
#define LOGMU (p[parindex[1]]) // baseline birth and death rate
#define LOGIOTA (p[parindex[2]]) // import rate
@@ -43,13 +31,13 @@
#define REPORTS (y[obsindex[0]])
-void _sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+void binomial_dmeasure (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) {
*lik = dbinom(REPORTS,CASE,exp(LOGRHO),give_log);
}
-void _sir_binom_rmeasure (double *y, double *x, double *p,
+void binomial_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
REPORTS = rbinom(CASE,exp(LOGRHO));
@@ -61,10 +49,10 @@
// 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)
+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
@@ -75,6 +63,9 @@
int nseas = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
double seasonality[nseas];
+ void (*eval_basis)(double,double,int,int,double*);
+ void (*reulmult)(int,double,double*,double,double*);
+ double (*dot)(int,const double *, const double *);
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -84,9 +75,13 @@
popsize = exp(LOGPOPSIZE);
beta_var = beta_sd*beta_sd;
+ eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
+ reulmult = (void (*)(int,double,double*,double,double*)) R_GetCCallable("pomp","reulermultinom");
+ dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
+
if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+ (*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
+ beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
// test to make sure the parameters and state variable values are sane
if (!(R_FINITE(beta)) ||
@@ -119,9 +114,9 @@
// compute the transition numbers
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]);
+ (*reulmult)(2,SUSC,&rate[1],dt,&trans[1]);
+ (*reulmult)(2,INFD,&rate[3],dt,&trans[3]);
+ (*reulmult)(1,RCVD,&rate[5],dt,&trans[5]);
// balance the equations
SUSC += trans[0]-trans[1]-trans[2];
@@ -139,9 +134,9 @@
#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)
+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
@@ -151,6 +146,8 @@
int nseas = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
double seasonality[nseas];
+ void (*eval_basis)(double,double,int,int,double*);
+ double (*dot)(int,const double *, const double *);
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -158,9 +155,12 @@
iota = exp(LOGIOTA);
popsize = exp(LOGPOPSIZE);
+ eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
+ dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
+
if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+ (*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
+ beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
// compute the transition rates
rate[0] = mu*popsize; // birth into susceptible class
@@ -207,175 +207,3 @@
#undef NBASIS
#undef DEGREE
#undef PERIOD
-
-#define LOGGAMMA (p[parindex[0]]) // recovery rate
-#define LOGMU (p[parindex[1]]) // death rate
-#define LOGIOTA (p[parindex[2]]) // import rate
-#define LOGBETA (p[parindex[3]]) // transmission rate
-#define LOGNU (p[parindex[4]]) // birth rate
-#define NBASIS (p[parindex[5]]) // number of periodic B-spline basis functions
-#define DEGREE (p[parindex[6]]) // degree of periodic B-spline basis functions
-#define PERIOD (p[parindex[7]]) // period of B-spline basis functions
-
-#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 POPN (x[stateindex[3]]) // population size
-#define CASE (x[stateindex[4]]) // accumulated cases
-
-double _sir_rates (int j, double t, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovar, double *covar) {
- double beta;
- double rate = 0.0;
- int nseas = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
-
- switch (j) {
- case 1: // birth
- rate = exp(LOGNU)*POPN;
- break;
- case 2: // susceptible death
- rate = exp(LOGMU)*SUSC;
- break;
- case 3: // infection
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
- rate = (beta*INFD+exp(LOGIOTA))*SUSC/POPN;
- break;
- case 4: // infected death
- rate = exp(LOGMU)*INFD;
- break;
- case 5: // recovery
- rate = exp(LOGGAMMA)*INFD;
- break;
- case 6: // recovered death
- rate = exp(LOGMU)*RCVD;
- break;
- default:
- error("unrecognized rate code %d",j);
- break;
- }
- return rate;
-}
-
-#undef LOGGAMMA
-#undef LOGMU
-#undef LOGIOTA
-#undef LOGBETA
-#undef LOGNU
-#undef NBASIS
-#undef DEGREE
-#undef PERIOD
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef POPN
-#undef CASE
-
-#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 LOGRHO (p[parindex[6]]) // reporting probability
-#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
-#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
-#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
-
-#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;
- int nseas = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
-
- periodic_bspline_basis_eval(t1,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&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)) ||
- (nseas <= 0)) {
- *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 SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
-#undef BIRTHS
-#undef dW
-
-#undef LOGGAMMA
-#undef LOGMU
-#undef LOGIOTA
-#undef LOGBETA
-#undef LOGBETA_SD
-#undef LOGPOPSIZE
-#undef LOGRHO
-#undef NBASIS
-#undef DEGREE
-#undef PERIOD
More information about the pomp-commits
mailing list