Skip to content

Commit

Permalink
Merge pull request #2 from stewid/SIR
Browse files Browse the repository at this point in the history
Add 'SIR' compartment model
  • Loading branch information
stewid committed Nov 23, 2016
2 parents 8febf95 + ee326bd commit 83ef763
Show file tree
Hide file tree
Showing 19 changed files with 809 additions and 41 deletions.
6 changes: 4 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@ Description: Livestock movements are important for the spread of many
can be extended with user defined models.
Acknowledgements: This work was financially supported by the Swedish
Research Council within the UPMARC Linnaeus center of Excellence
(Pavol Bauer and Stefan Engblom), The Swedish Research Council
Formas (Stefan Engblom and Stefan Widgren), and The Swedish Board
(Pavol Bauer and Stefan Engblom), the Swedish Research Council
Formas (Stefan Engblom and Stefan Widgren), and the Swedish Board
of Agriculture (Stefan Widgren).
License: GPL-3
URL: https://github.com/stewid/SimInf
Expand All @@ -33,12 +33,14 @@ Depends:
Imports:
graphics,
stats,
utils,
Matrix
Collate:
'AllGenerics.R'
'check_arguments.R'
'scheduled_events.R'
'siminf_model.R'
'SIR.R'
'SISe.R'
'SISe3.R'
'SISe3_sp.R'
Expand Down
7 changes: 7 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,16 +1,21 @@
# Generated by roxygen2: do not edit by hand

export(SIR)
export(SISe)
export(SISe3)
export(SISe3_sp)
export(SISe_sp)
export(demo_model)
export(distance_matrix)
export(events_SIR)
export(infected)
export(prevalence)
export(recovered)
export(scheduled_events)
export(siminf_model)
export(susceptible)
export(u0_SIR)
exportClasses(SIR)
exportClasses(SISe)
exportClasses(SISe3)
exportClasses(SISe3_sp)
Expand All @@ -20,6 +25,7 @@ exportClasses(siminf_model)
exportMethods(infected)
exportMethods(plot)
exportMethods(prevalence)
exportMethods(recovered)
exportMethods(run)
exportMethods(run_outer)
exportMethods(show)
Expand All @@ -36,4 +42,5 @@ importFrom(graphics,plot)
importFrom(graphics,title)
importFrom(stats,terms)
importFrom(stats,xtabs)
importFrom(utils,data)
useDynLib(SimInf, .registration=TRUE)
46 changes: 43 additions & 3 deletions R/AllGenerics.R
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
## siminf, a framework for stochastic disease spread simulations
## SimInf, a framework for stochastic disease spread simulations
## Copyright (C) 2015 Pavol Bauer
## Copyright (C) 2015 Stefan Engblom
## Copyright (C) 2015 Stefan Widgren
## Copyright (C) 2015 - 2016 Stefan Engblom
## Copyright (C) 2015 - 2016 Stefan Widgren
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
Expand Down Expand Up @@ -172,6 +172,46 @@ setGeneric("susceptible",
setGeneric("infected",
function(model, ...) standardGeneric("infected"))

##' Recovered
##'
##' Extracts the number of recovered
##' @rdname recovered-methods
##' @docType methods
##' @param model The \code{model} to extract the recovered from
##' @param ... Additional arguments affecting the measure
##' @param i Indices specifying the nodes to include when extracting
##' the number of recovered. Default is NULL, which includes all nodes.
##' @param by The number to increment the sequence of time points
##' starting from 1. Default is 1, which gives the number of
##' recovered at every time point.
##' @keywords methods
##' @export
##' @examples
##' ## Create a 'SIR' demo model with 5 nodes and initialize
##' ## it to run over 10 days.
##' model <- demo_model(nodes = 5, days = 10, model = "SIR")
##'
##' ## Run the model and save the result
##' result <- run(model)
##'
##' ## Extract the number of recovered individuals in each
##' ## node after each time step in the simulation
##' recovered(result)
##'
##' ## Extract the number of recovered individuals in the
##' ## first node after each time step in the simulation
##' recovered(result, i = 1)
##'
##' ## Extract the number of recovered individuals in the
##' ## first and third node after each time step in the simulation
##' recovered(result, i = c(1, 3))
##'
##' ## Extract the number of recovered individuals in the first
##' ## and third node after every other time step in the simulation
##' recovered(result, i = c(1, 3), by = 2)
setGeneric("recovered",
function(model, ...) standardGeneric("recovered"))

##' Prevalence
##'
##' Calculate the proportion infected individuals
Expand Down
221 changes: 221 additions & 0 deletions R/SIR.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,221 @@
## SimInf, a framework for stochastic disease spread simulations
## Copyright (C) 2015 - 2016 Stefan Widgren
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.

##' Class \code{"SIR"}
##'
##' Class to handle the SIR \code{\link{siminf_model}}.
##' @include siminf_model.R
##' @include AllGenerics.R
##' @export
setClass("SIR", contains = c("siminf_model"))

##' Create a SIR model
##'
##' Create a SIR model to be used by the simulation framework.
##'
##'
##' The argument \code{u0} must be a \code{data.frame} with one row for
##' each node with the following columns:
##' \describe{
##' \item{S}{The number of sucsceptible in each node}
##' \item{I}{The number of infected in each node}
##' \item{R}{The number of recovered in each node}
##' }
##'
##' @param u0 A \code{data.frame} with the initial state in each node,
##' see details.
##' @param tspan An increasing sequence of points in time where the
##' state of the system is to be returned.
##' @param events a \code{data.frame} with the scheduled events, see
##' \code{\link{siminf_model}}.
##' @param beta The transmission rate from susceptible to infected.
##' @param gamma The recovery rate from infected to recovered.
##' @return \code{SIR}
##' @include check_arguments.R
##' @export
SIR <- function(u0,
tspan,
events = NULL,
beta = NULL,
gamma = NULL)
{
compartments <- c("S", "I", "R")

## Check arguments.

## Check u0
if (!is.data.frame(u0))
stop("'u0' must be a data.frame")
if (!all(compartments %in% names(u0)))
stop("Missing columns in u0")
u0 <- u0[, compartments]

## Check for non-numeric parameters
check_gdata_arg(beta, gamma)

## Arguments seems ok...go on

E <- Matrix(c(1, 1,
0, 1,
0, 1),
nrow = 3,
ncol = 2,
byrow = TRUE,
sparse = TRUE)
E <- as(E, "dgCMatrix")
colnames(E) <- as.character(1:2)
rownames(E) <- compartments

N <- matrix(integer(0), nrow = 0, ncol = 0)

G <- Matrix(c(1, 1,
1, 1),
nrow = 2,
ncol = 2,
byrow = TRUE,
sparse = TRUE)
G <- as(G, "dgCMatrix")
colnames(G) <- as.character(1:2)
rownames(G) <- c("S -> I", "I -> R")

S <- Matrix(c(-1, 0,
1, -1,
0, 1),
nrow = 3,
ncol = 2,
byrow = TRUE,
sparse = TRUE)
S <- as(S, "dgCMatrix")
colnames(S) <- as.character(1:2)
rownames(S) <- compartments

v0 <- matrix(numeric(0), nrow = 0, ncol = nrow(u0))
storage.mode(v0) <- "double"

ldata <- matrix(numeric(0), nrow = 0, ncol = nrow(u0))
storage.mode(ldata) <- "double"

gdata <- c(beta, gamma)
storage.mode(gdata) <- "double"
names(gdata) <- c("beta", "gamma")

model <- siminf_model(G = G,
S = S,
E = E,
N = N,
tspan = tspan,
events = events,
ldata = ldata,
gdata = gdata,
u0 = u0,
v0 = v0)

return(as(model, "SIR"))
}

##' @rdname susceptible-methods
##' @export
setMethod("susceptible",
signature("SIR"),
function(model, i = NULL, by = 1, ...) {
if (identical(dim(model@U), c(0L, 0L)))
stop("Please run the model first, the 'U' matrix is empty")

ii <- seq(from = 1, to = dim(model@U)[1], by = 3)
if (!is.null(i))
ii <- ii[i]
j <- seq(from = 1, to = dim(model@U)[2], by = by)
as.matrix(model@U[ii, j, drop = FALSE])
}
)

##' @rdname infected-methods
##' @export
setMethod("infected",
signature("SIR"),
function(model, i = NULL, by = 1, ...) {
if (identical(dim(model@U), c(0L, 0L)))
stop("Please run the model first, the 'U' matrix is empty")

ii <- seq(from = 2, to = dim(model@U)[1], by = 3)
if (!is.null(i))
ii <- ii[i]
j <- seq(from = 1, to = dim(model@U)[2], by = by)
as.matrix(model@U[ii, j, drop = FALSE])
}
)

##' @rdname recovered-methods
##' @export
setMethod("recovered",
signature("SIR"),
function(model, i = NULL, by = 1, ...) {
if (identical(dim(model@U), c(0L, 0L)))
stop("Please run the model first, the 'U' matrix is empty")

ii <- seq(from = 3, to = dim(model@U)[1], by = 3)
if (!is.null(i))
ii <- ii[i]
j <- seq(from = 1, to = dim(model@U)[2], by = by)
as.matrix(model@U[ii, j, drop = FALSE])
}
)

##' @name plot-methods
##' @aliases plot plot-methods plot,SIR-method
##' @importFrom graphics plot
##' @export
setMethod("plot",
signature(x = "SIR"),
function(x, t0 = NULL, ...)
{
callNextMethod(x, t0 = t0, legend = c("S", "I", "R"), ...)
}
)

##' Scheduled events example data for the \code{SIR} model
##'
##' Synthetic scheduled events data to demonstrate the \code{SIR}
##' model. The data contains 466692 events for 1600 nodes over 365 * 4
##' days.
##' @return A \code{data.frame}
##' @keywords methods
##' @importFrom utils data
##' @export
events_SIR <- function() {
utils::data(events_SISe3, envir = environment())
events_SISe3$select[events_SISe3$event == 0] <- 2
events_SISe3$select[events_SISe3$event == 1] <- 1
events_SISe3 <- events_SISe3[events_SISe3$event != 2,]
events_SISe3$select[events_SISe3$event == 3] <- 2
events_SISe3
}

##' Example data to initialize the \code{SIR} model
##'
##' Synthetic init data for 1600 nodes to demonstrate the \code{SIR}
##' model.
##' @return A \code{data.frame}
##' @keywords methods
##' @importFrom utils data
##' @export
u0_SIR <- function() {
utils::data(u0_SISe3, envir = environment())
u0_SISe3$S <- u0_SISe3$S_1 + u0_SISe3$S_2 + u0_SISe3$S_3
u0_SISe3$I <- u0_SISe3$I_1 + u0_SISe3$I_2 + u0_SISe3$I_3
u0_SISe3$R <- 0
u0_SISe3[, c("S", "I", "R")]
}
12 changes: 11 additions & 1 deletion R/demo_model.R
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
##' plot(result)
demo_model <- function(nodes = 1,
days = 1000,
model = c("SISe", "SISe3", "SISe_sp"))
model = c("SISe", "SISe3", "SISe_sp", "SIR"))
{
## Check 'nodes' argument
if (!is.numeric(nodes))
Expand Down Expand Up @@ -137,6 +137,16 @@ demo_model <- function(nodes = 1,
end_t4 = 365,
distance = distance,
coupling = 0.002)
} else if (identical(model, "SIR")) {
u0 <- data.frame(S = rep(392, nodes),
I = rep(8, nodes),
R = rep(0, nodes))

model <- SIR(u0 = u0,
tspan = seq_len(days) - 1,
events = NULL,
beta = 0.16,
gamma = 0.077)
}

model
Expand Down
Loading

0 comments on commit 83ef763

Please sign in to comment.