Skip to content

Commit

Permalink
adding general SIR model
Browse files Browse the repository at this point in the history
  • Loading branch information
Lam Ho authored and Lam Ho committed Mar 21, 2018
1 parent 6fdcc9b commit 564ba03
Show file tree
Hide file tree
Showing 7 changed files with 82 additions and 30 deletions.
4 changes: 2 additions & 2 deletions R/RcppExports.R
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ SEIR_Cpp <- function(t, alpha, beta, kappa, S0, E0, I0, Ap1, Bp1, Cp1, direction
.Call('_MultiBD_SEIR_Cpp', PACKAGE = 'MultiBD', t, alpha, beta, kappa, S0, E0, I0, Ap1, Bp1, Cp1, direction, nblocks, tol, Lmax, computeMode, nThreads)
}

SIR_Cpp <- function(t, alpha, beta, S0, I0, Ap1, Bp1, direction, nblocks, tol, Lmax, computeMode, nThreads) {
.Call('_MultiBD_SIR_Cpp', PACKAGE = 'MultiBD', t, alpha, beta, S0, I0, Ap1, Bp1, direction, nblocks, tol, Lmax, computeMode, nThreads)
SIR_Cpp <- function(t, alpha, beta, S0, I0, Ap1, Bp1, direction, powS, powI_inf, powI_rem, nblocks, tol, Lmax, computeMode, nThreads) {
.Call('_MultiBD_SIR_Cpp', PACKAGE = 'MultiBD', t, alpha, beta, S0, I0, Ap1, Bp1, direction, powS, powI_inf, powI_rem, nblocks, tol, Lmax, computeMode, nThreads)
}

tb_lt_invert_Cpp <- function(t, lambda1, lambda2, lambda3, Ap1, Bp1, Cp1, direction, nblocks, tol, Lmax, computeMode, nThreads) {
Expand Down
21 changes: 18 additions & 3 deletions R/SIR_prob.R
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,16 @@
#' @param nSI number of infection events
#' @param nIR number of removal events
#' @param direction direction of the transition probabilities (either \code{Forward} or \code{Backward})
#' @param power the power of the general SIR model (see Note for more details)
#' @param nblocks number of blocks
#' @param tol tolerance
#' @param computeMode computation mode
#' @param nThreads number of threads
#' @return a matrix of the transition probabilities
#' @note The infection rate and the removal rate of a general SIR model
#' are \code{beta*S^{powS}*I^{powI_inf}} and \code{alpha*I^{powI_rem}} respectively.
#' The parameter \code{power} is a list of \code{powS, powI_inf, powI_rem}.
#' Their default values are \code{powS = powI_inf = pwoI_rem = 1}, which correspond to the classic SIR model.
#' @examples
#' data(Eyam)
#'
Expand Down Expand Up @@ -49,11 +54,20 @@
#' @export

SIR_prob <- function(t, alpha, beta, S0, I0, nSI, nIR, direction = c("Forward","Backward"),
nblocks = 20, tol = 1e-12, computeMode = 0, nThreads = 4) {
power = NULL, nblocks = 20, tol = 1e-12, computeMode = 0, nThreads = 4) {

direction <- match.arg(direction)
dir = 0
if (direction == "Backward") dir = 1
if (is.null(power)) {
powS = 1
powI_inf = 1
powI_rem =1
} else {
powS = power$powS
powI_inf = power$powI_inf
powI_rem = power$powI_rem
}

################################
### t is too small
Expand All @@ -69,8 +83,9 @@ SIR_prob <- function(t, alpha, beta, S0, I0, nSI, nIR, direction = c("Forward","
}

Lmax = nblocks # initialize Lmax
res = matrix(SIR_Cpp(t, alpha, beta, S0, I0, nSI + 1, nIR + 1, dir, nblocks, tol,
Lmax, computeMode, nThreads),
res = matrix(SIR_Cpp(t, alpha, beta, S0, I0, nSI + 1, nIR + 1, dir,
powS, powI_inf, powI_rem,
nblocks, tol, Lmax, computeMode, nThreads),
nrow = nSI + 1, byrow = T)

rownames(res) = 0:nSI # Infection events
Expand Down
11 changes: 10 additions & 1 deletion man/SIR_prob.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

11 changes: 7 additions & 4 deletions src/RcppExports.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,8 @@ BEGIN_RCPP
END_RCPP
}
// SIR_Cpp
std::vector<double> SIR_Cpp(const double t, const double alpha, const double beta, const long int S0, const long int I0, const int Ap1, const int Bp1, const int direction, const int nblocks, const double tol, int& Lmax, const int computeMode, const int nThreads);
RcppExport SEXP _MultiBD_SIR_Cpp(SEXP tSEXP, SEXP alphaSEXP, SEXP betaSEXP, SEXP S0SEXP, SEXP I0SEXP, SEXP Ap1SEXP, SEXP Bp1SEXP, SEXP directionSEXP, SEXP nblocksSEXP, SEXP tolSEXP, SEXP LmaxSEXP, SEXP computeModeSEXP, SEXP nThreadsSEXP) {
std::vector<double> SIR_Cpp(const double t, const double alpha, const double beta, const long int S0, const long int I0, const int Ap1, const int Bp1, const int direction, const double powS, const double powI_inf, const double powI_rem, const int nblocks, const double tol, int& Lmax, const int computeMode, const int nThreads);
RcppExport SEXP _MultiBD_SIR_Cpp(SEXP tSEXP, SEXP alphaSEXP, SEXP betaSEXP, SEXP S0SEXP, SEXP I0SEXP, SEXP Ap1SEXP, SEXP Bp1SEXP, SEXP directionSEXP, SEXP powSSEXP, SEXP powI_infSEXP, SEXP powI_remSEXP, SEXP nblocksSEXP, SEXP tolSEXP, SEXP LmaxSEXP, SEXP computeModeSEXP, SEXP nThreadsSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Expand All @@ -92,12 +92,15 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const int >::type Ap1(Ap1SEXP);
Rcpp::traits::input_parameter< const int >::type Bp1(Bp1SEXP);
Rcpp::traits::input_parameter< const int >::type direction(directionSEXP);
Rcpp::traits::input_parameter< const double >::type powS(powSSEXP);
Rcpp::traits::input_parameter< const double >::type powI_inf(powI_infSEXP);
Rcpp::traits::input_parameter< const double >::type powI_rem(powI_remSEXP);
Rcpp::traits::input_parameter< const int >::type nblocks(nblocksSEXP);
Rcpp::traits::input_parameter< const double >::type tol(tolSEXP);
Rcpp::traits::input_parameter< int& >::type Lmax(LmaxSEXP);
Rcpp::traits::input_parameter< const int >::type computeMode(computeModeSEXP);
Rcpp::traits::input_parameter< const int >::type nThreads(nThreadsSEXP);
rcpp_result_gen = Rcpp::wrap(SIR_Cpp(t, alpha, beta, S0, I0, Ap1, Bp1, direction, nblocks, tol, Lmax, computeMode, nThreads));
rcpp_result_gen = Rcpp::wrap(SIR_Cpp(t, alpha, beta, S0, I0, Ap1, Bp1, direction, powS, powI_inf, powI_rem, nblocks, tol, Lmax, computeMode, nThreads));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -129,7 +132,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_MultiBD_bb_lt_invert_Cpp", (DL_FUNC) &_MultiBD_bb_lt_invert_Cpp, 11},
{"_MultiBD_bbd_lt_invert_Cpp", (DL_FUNC) &_MultiBD_bbd_lt_invert_Cpp, 16},
{"_MultiBD_SEIR_Cpp", (DL_FUNC) &_MultiBD_SEIR_Cpp, 16},
{"_MultiBD_SIR_Cpp", (DL_FUNC) &_MultiBD_SIR_Cpp, 13},
{"_MultiBD_SIR_Cpp", (DL_FUNC) &_MultiBD_SIR_Cpp, 16},
{"_MultiBD_tb_lt_invert_Cpp", (DL_FUNC) &_MultiBD_tb_lt_invert_Cpp, 13},
{NULL, NULL, 0}
};
Expand Down
7 changes: 5 additions & 2 deletions src/SIR_Cpp.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
// [[Rcpp::export]]
std::vector<double> SIR_Cpp(const double t, const double alpha, const double beta,
const long int S0, const long int I0, const int Ap1, const int Bp1, const int direction,
const double powS, const double powI_inf, const double powI_rem,
const int nblocks, const double tol, int& Lmax, const int computeMode, const int nThreads) {

const int matsize = Ap1*Bp1;
Expand All @@ -12,8 +13,10 @@ std::vector<double> SIR_Cpp(const double t, const double alpha, const double bet
for (int b=0; b<Bp1; ++b) {
double Spop = std::max(0.0, S0-a+0.0);
double Ipop = std::max(0.0, a+I0-b+0.0);
lambda1[a + b*Ap1] = beta*Spop*Ipop; // Infection rate is beta*S*I
lambda2[a + b*Ap1] = alpha*Ipop;
lambda1[a + b*Ap1] = beta*pow(Spop, powS)*pow(Ipop, powI_inf);
// Infection rate is beta*S^{powS}*I^{powI_inf}
lambda2[a + b*Ap1] = alpha*pow(Ipop, powI_rem);
// Removal rate is alpha*I^{powI_rem}
}
}

Expand Down
1 change: 1 addition & 0 deletions src/bbd.h
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,7 @@ std::vector<double> bbd_lt_invert_Cpp(double t, const int a0, const int b0,

std::vector<double> SIR_Cpp(const double t, const double alpha, const double beta,
const long int S0, const long int I0, const int Ap1, const int Bp1, const int direction,
const double powS, const double powI_inf, const double powI_rem,
const int nblocks, const double tol, int& Lmax, const int computeMode, const int nThreads);

std::vector<double> bb_lt_invert_Cpp(double t, const std::vector<double>& lambda1,
Expand Down
57 changes: 39 additions & 18 deletions tests/testthat/test_SIR.R
Original file line number Diff line number Diff line change
@@ -1,47 +1,68 @@
library(testthat)

test_that("Transition probabilities of SIR model", {

data(Eyam)
data <- Eyam

alpha = 3.204
beta = 0.019
k = 1
N <- data$S[1] + data$I[1] + data$R[1]

power = list(powS = 0.9, powI_inf = 0.8, powI_rem = 0.9)

p1 <- bbd_prob(
t = data$time[k + 1] - data$time[k], # Time increment
a0 = data$R[k], b0 = data$I[k], # From: R(t_k), I(t_k)
brates1 <- function(a,b){ 0 },
brates2 <- function(a,b){ beta * max(N - a - b, 0) * b },
drates2 <- function(a,b){ 0 },
a0 = data$R[k], b0 = data$I[k], # From: R(t_k), I(t_k)
brates1 <- function(a,b){ 0 },
brates2 <- function(a,b){ beta * max(N - a - b, 0) * b },
drates2 <- function(a,b){ 0 },
trans21 <- function(a,b){ alpha * b },
A = data$R[k + 1], B = data$R[k + 1] + data$I[k] - data$R[k],
computeMode = 4, nblocks = 80 # Compute using 4 threads
)[data$R[k + 1] - data$R[k] + 1, data$I[k + 1] + 1]

p2 <- dbd_prob( # Compute the transition probability matrix
t = data$time[k + 1] - data$time[k], # Time increment
a0 = data$S[k], b0 = data$I[k], # From: S(t_k), I(t_k)
drates1 <- function(a, b) { 0 },
brates2 <- function(a, b) { 0 },
drates2 <- function(a, b) { alpha * b },
a0 = data$S[k], b0 = data$I[k], # From: S(t_k), I(t_k)
drates1 <- function(a, b) { 0 },
brates2 <- function(a, b) { 0 },
drates2 <- function(a, b) { alpha * b },
trans12 <- function(a, b) { beta * a * b },
a = data$S[k + 1], B = data$S[k] + data$I[k] - data$S[k + 1],
computeMode = 4, nblocks = 80 # Compute using 4 threads
)[1, data$I[k + 1] + 1]
p3 <- SIR_prob(

p3 <- SIR_prob(
t = data$time[k + 1] - data$time[k], # Time increment
alpha = alpha, beta = beta,
S0 = data$S[k], I0 = data$I[k], # From: R(t_k), I(t_k)
S0 = data$S[k], I0 = data$I[k], # From: R(t_k), I(t_k)
nSI = data$S[k] - data$S[k + 1], nIR = data$R[k + 1] - data$R[k],
direction = "Forward",
computeMode = 4, nblocks = 80 # Compute using 4 threads
)[data$S[k] - data$S[k + 1] + 1, data$R[k + 1] - data$R[k] + 1]


p4 <- dbd_prob( # Compute the transition probability matrix
t = data$time[k + 1] - data$time[k], # Time increment
a0 = data$S[k], b0 = data$I[k], # From: S(t_k), I(t_k)
drates1 <- function(a, b) { 0 },
brates2 <- function(a, b) { 0 },
drates2 <- function(a, b) { alpha * b^0.9 },
trans12 <- function(a, b) { beta * a^0.9 * b^0.8 },
a = data$S[k + 1], B = data$S[k] + data$I[k] - data$S[k + 1],
computeMode = 4, nblocks = 80 # Compute using 4 threads
)[1, data$I[k + 1] + 1]

p5 <- SIR_prob(
t = data$time[k + 1] - data$time[k], # Time increment
alpha = alpha, beta = beta,
S0 = data$S[k], I0 = data$I[k], # From: R(t_k), I(t_k)
nSI = data$S[k] - data$S[k + 1], nIR = data$R[k + 1] - data$R[k],
direction = "Forward", power = power,
computeMode = 4, nblocks = 80 # Compute using 4 threads
)[data$S[k] - data$S[k + 1] + 1, data$R[k + 1] - data$R[k] + 1]

expect_equal(0.0, abs(p1-p2), 1E-5)
expect_equal(0.0, abs(p2-p3), 1E-8)
})
expect_equal(0.0, abs(p4-p5), 1E-8)
})

0 comments on commit 564ba03

Please sign in to comment.