Nothing
tests/sample_select.R
In SimInf: A Framework for Data-Driven Stochastic Disease Spread Simulations
## This file is part of SimInf, a framework for stochastic
## disease spread simulations.
##
## Copyright (C) 2015 Pavol Bauer
## Copyright (C) 2017 -- 2019 Robin Eriksson
## Copyright (C) 2015 -- 2019 Stefan Engblom
## Copyright (C) 2015 -- 2021 Stefan Widgren
##
## SimInf 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.
##
## SimInf 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 <https://www.gnu.org/licenses/>.
library(SimInf)
library(tools)
source("util/check.R")
## Specify the number of threads to use.
max_threads <- set_num_threads(1)
## For debugging
sessionInfo()
## Test that sample_select in events.c works
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, two individuals are moved to node = 1. This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 2,
proportion = 1,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(run(model))
set_num_threads(1)
check_error(res, "Unable to sample individuals for event.")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, -1 individuals are moved to node = 1. This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = -1,
proportion = 1,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(run(model))
set_num_threads(1)
check_error(res, "Unable to sample individuals for event.")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 10 individuals are moved to node = 1.
## This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 10,
select = 4,
shift = 0)
## We should not be able to create model with prop = 10
res <- assertError(SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0))
check_error(res, "prop must be in the range 0 <= prop <= 1", FALSE)
## Replace proportion = 10 to proportion = 1
events$proportion <- 1
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Replace proportion = 1 with proportion = 10
model@events@proportion <- 10
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
set_num_threads(1)
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of -1 individuals are moved to node = 1.
## This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = -1,
select = 4,
shift = 0)
## We should not be able to create model with proportion = -1
res <- assertError(SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0))
check_error(res, "prop must be in the range 0 <= prop <= 1", FALSE)
## Replace proportion = -1 with proportion = 0
events$proportion <- 0
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Replace proportion = 0 with proportion = -1
model@events@proportion <- -1
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
set_num_threads(1)
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 0 individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, proportion of all (1) individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, Nkind = 1, one individual is moved to node = 1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state and 8 individuals in
## infected state in node = 2, with a zero probability of becoming
## infected.
##
## At t = 1, Nkind = 2, one individual is moved to node = 1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 8),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 0,
gamma_2 = 0,
gamma_3 = 0,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 8L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 2L, 7L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 2L, 7L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
set.seed(1)
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set.seed(1)
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, one individual is moved to node = 1 using a select matrix
## that can only select a susceptible from S_1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 1,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U <- structure(
c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L),
.Dim = c(12L, 11L))
res <- run(model)
stopifnot(identical(res@U, U))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- run(model)
set_num_threads(1)
stopifnot(identical(res@U, U))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## 6 individuals start in the susceptible compartments in node = 2,
## with a zero probability of becoming infected.
##
## At t = 1, three individuals are moved to node = 1 using a select
## matrix that can select from the susceptible S_1, S_2 and S_3.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 2),
I_2 = c(0, 0),
S_3 = c(0, 2),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 3,
proportion = 0,
select = 1,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Add S_2 and S_3 to first column in select matrix
model@events@E[3, 1] <- 1
model@events@E[5, 1] <- 1
S_expected <- c(0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)
res <- run(model)
S_observed <- colSums(trajectory(res, compartments = c("S_1", "S_2", "S_3"),
index = 1, format = "matrix"))
stopifnot(identical(S_observed, S_expected))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- run(model)
set_num_threads(1)
S_observed <- colSums(trajectory(res, compartments = c("S_1", "S_2", "S_3"),
index = 1, format = "matrix"))
stopifnot(identical(S_observed, S_expected))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## 10 individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 0.5 individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 10),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 0.5,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## In this example, a proportion of 0.5 of 10 animals are scheduled to
## move and this results in 6 being moved because the number of
## animals moved is now sampled from a binomial distribution:
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 6L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 6L))
set_num_threads(1)
}
## With a very small proportion, most often, 0 animals are moved:
model@events@proportion <- 0.01
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 0L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 0L))
set_num_threads(1)
}
## In some cases >0 are moved when 0.01 of 10 animals are scheduled to
## move:
set.seed(17)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 1L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(17)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 1L))
set_num_threads(1)
}
## 1 Node in a SIR model
##
## One individual starts in S, one in I and one in R with a zero
## probability of becoming infected or recover.
##
## At t = 1, one individual exit, with a higher weight for being
## sampled from the R compartment.
model <- SIR(u0 = data.frame(S = 1, I = 1, R = 1),
tspan = 1,
events = data.frame(event = 0,
time = 1,
node = 1,
dest = 0,
n = 1,
proportion = 0,
select = 4,
shift = 0),
beta = 0,
gamma = 0)
## With equal weight the individual is sampled from the S compartment.
set.seed(2)
stopifnot(identical(run(model)@U, structure(c(0L, 1L, 1L), .Dim = c(3L, 1L))))
## With non-equal weight the individual is sampled from the R
## compartment (using the same seed).
model@events@E[3, 4] <- 1000
set.seed(2)
stopifnot(identical(run(model)@U, structure(c(1L, 1L, 0L), .Dim = c(3L, 1L))))
## Check that an error is raised if there are no non-zero elements in
## the selected column in E.
model <- SIR(u0 = data.frame(S = 1, I = 1, R = 1),
tspan = 1,
events = data.frame(event = 0,
time = 1,
node = 1,
dest = 0,
n = 1,
proportion = 0,
select = 4,
shift = 0),
beta = 0,
gamma = 0)
model@events@E <- as(matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0),
nrow = 3, ncol = 4,
dimnames = list(c("S", "I", "R"),
c("1", "2", "3", "4"))),
"dgCMatrix")
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to an enter event
model@events@event <- 1L
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to enter event and n = -1
model@events@n <- -1L
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to enter event and proportion = -1
model@events@E <- as(matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1),
nrow = 3, ncol = 4,
dimnames = list(c("S", "I", "R"),
c("1", "2", "3", "4"))),
"dgCMatrix")
model@events@n <- 0L
model@events@proportion <- -1
res <- assertError(.Call(SimInf:::SIR_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
## Change to enter event and proportion = 1.1
model@events@n <- 0L
model@events@proportion <- 1.1
res <- assertError(.Call(SimInf:::SIR_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
## Check that proportion works for an enter event.
model@events@n <- 0L
model@events@proportion <- 0.1
model@events@select <- 1L
set.seed(3)
stopifnot(identical(run(model)@U, structure(c(2L, 1L, 1L), .Dim = c(3L, 1L))))
## Check that an enter event fails with shift without N.
model@events@n <- 1L
model@events@shift <- 1L
model@events@proportion <- 0
res <- assertError(run(model))
check_error(res, "'N' is invalid.")
## Check that shift fails for an enter event when the shift is out of
## bounds.
model@events@N <- matrix(c(3L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
res <- assertError(run(model))
check_error(res, "'shift' is out of bounds.")
model@events@N <- matrix(c(-1L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
res <- assertError(run(model))
check_error(res, "'shift' is out of bounds.")
## Check that shift works for an enter event.
model@events@N <- matrix(c(1L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
stopifnot(identical(run(model)@U, structure(c(1L, 2L, 1L), .Dim = c(3L, 1L))))
Try the SimInf package in your browser
Any scripts or data that you put into this service are public.
SimInf documentation built on Jan. 23, 2023, 5:43 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## This file is part of SimInf, a framework for stochastic
## disease spread simulations.
##
## Copyright (C) 2015 Pavol Bauer
## Copyright (C) 2017 -- 2019 Robin Eriksson
## Copyright (C) 2015 -- 2019 Stefan Engblom
## Copyright (C) 2015 -- 2021 Stefan Widgren
##
## SimInf 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.
##
## SimInf 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 <https://www.gnu.org/licenses/>.
library(SimInf)
library(tools)
source("util/check.R")
## Specify the number of threads to use.
max_threads <- set_num_threads(1)
## For debugging
sessionInfo()
## Test that sample_select in events.c works
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, two individuals are moved to node = 1. This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 2,
proportion = 1,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(run(model))
set_num_threads(1)
check_error(res, "Unable to sample individuals for event.")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, -1 individuals are moved to node = 1. This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = -1,
proportion = 1,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(run(model))
set_num_threads(1)
check_error(res, "Unable to sample individuals for event.")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 10 individuals are moved to node = 1.
## This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 10,
select = 4,
shift = 0)
## We should not be able to create model with prop = 10
res <- assertError(SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0))
check_error(res, "prop must be in the range 0 <= prop <= 1", FALSE)
## Replace proportion = 10 to proportion = 1
events$proportion <- 1
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Replace proportion = 1 with proportion = 10
model@events@proportion <- 10
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
set_num_threads(1)
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of -1 individuals are moved to node = 1.
## This should fail.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = -1,
select = 4,
shift = 0)
## We should not be able to create model with proportion = -1
res <- assertError(SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0))
check_error(res, "prop must be in the range 0 <= prop <= 1", FALSE)
## Replace proportion = -1 with proportion = 0
events$proportion <- 0
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Replace proportion = 0 with proportion = -1
model@events@proportion <- -1
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- assertError(.Call(SimInf:::SISe3_run, model, NULL))
set_num_threads(1)
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 0 individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## One individual start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, proportion of all (1) individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 1),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, Nkind = 1, one individual is moved to node = 1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state and 8 individuals in
## infected state in node = 2, with a zero probability of becoming
## infected.
##
## At t = 1, Nkind = 2, one individual is moved to node = 1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 8),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:2,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 0,
gamma_2 = 0,
gamma_3 = 0,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 8L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 2L, 7L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 2L, 7L, 0L, 0L, 0L, 0L),
.Dim = c(12L, 3L))
set.seed(1)
result <- run(model)
stopifnot(identical(model@G, result@G))
stopifnot(identical(model@S, result@S))
stopifnot(identical(trajectory(result, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result@ldata))
stopifnot(identical(model@tspan, result@tspan))
stopifnot(identical(model@u0, result@u0))
stopifnot(identical(model@events, result@events))
if (SimInf:::have_openmp() && max_threads > 1) {
set.seed(1)
set_num_threads(2)
result_omp <- run(model)
set_num_threads(1)
stopifnot(identical(model@G, result_omp@G))
stopifnot(identical(model@S, result_omp@S))
stopifnot(identical(trajectory(result_omp, format = "matrix"), U_expected))
stopifnot(identical(model@ldata, result_omp@ldata))
stopifnot(identical(model@tspan, result_omp@tspan))
stopifnot(identical(model@u0, result_omp@u0))
stopifnot(identical(model@events, result_omp@events))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## Two individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, one individual is moved to node = 1 using a select matrix
## that can only select a susceptible from S_1.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 1,
proportion = 0,
select = 1,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
U <- structure(
c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L),
.Dim = c(12L, 11L))
res <- run(model)
stopifnot(identical(res@U, U))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- run(model)
set_num_threads(1)
stopifnot(identical(res@U, U))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## 6 individuals start in the susceptible compartments in node = 2,
## with a zero probability of becoming infected.
##
## At t = 1, three individuals are moved to node = 1 using a select
## matrix that can select from the susceptible S_1, S_2 and S_3.
u0 <- data.frame(S_1 = c(0, 2),
I_1 = c(0, 0),
S_2 = c(0, 2),
I_2 = c(0, 0),
S_3 = c(0, 2),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 3,
proportion = 0,
select = 1,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## Add S_2 and S_3 to first column in select matrix
model@events@E[3, 1] <- 1
model@events@E[5, 1] <- 1
S_expected <- c(0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)
res <- run(model)
S_observed <- colSums(trajectory(res, compartments = c("S_1", "S_2", "S_3"),
index = 1, format = "matrix"))
stopifnot(identical(S_observed, S_expected))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
res <- run(model)
set_num_threads(1)
S_observed <- colSums(trajectory(res, compartments = c("S_1", "S_2", "S_3"),
index = 1, format = "matrix"))
stopifnot(identical(S_observed, S_expected))
}
## 2 Nodes
## 3 Age categories
## 2 Disease-states: Susceptible & Infected
##
## 10 individuals start in susceptible state in node = 2, with a zero
## probability of becoming infected.
##
## At t = 1, a proportion of 0.5 individuals are moved to node = 1.
u0 <- data.frame(S_1 = c(0, 10),
I_1 = c(0, 0),
S_2 = c(0, 0),
I_2 = c(0, 0),
S_3 = c(0, 0),
I_3 = c(0, 0))
events <- data.frame(event = 3,
time = 1,
node = 2,
dest = 1,
n = 0,
proportion = 0.5,
select = 4,
shift = 0)
model <- SISe3(u0 = u0,
tspan = 0:10,
events = events,
phi = rep(0, 2),
upsilon_1 = 0,
upsilon_2 = 0,
upsilon_3 = 0,
gamma_1 = 1,
gamma_2 = 1,
gamma_3 = 1,
alpha = 0,
beta_t1 = 1,
beta_t2 = 1,
beta_t3 = 1,
beta_t4 = 1,
end_t1 = 91,
end_t2 = 182,
end_t3 = 273,
end_t4 = 365,
epsilon = 0)
## In this example, a proportion of 0.5 of 10 animals are scheduled to
## move and this results in 6 being moved because the number of
## animals moved is now sampled from a binomial distribution:
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 6L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 6L))
set_num_threads(1)
}
## With a very small proportion, most often, 0 animals are moved:
model@events@proportion <- 0.01
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 0L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(42)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 0L))
set_num_threads(1)
}
## In some cases >0 are moved when 0.01 of 10 animals are scheduled to
## move:
set.seed(17)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 1L))
if (SimInf:::have_openmp() && max_threads > 1) {
set_num_threads(2)
set.seed(17)
res <- .Call(SimInf:::SISe3_run, model, NULL)
stopifnot(identical(res@U[1, 2], 1L))
set_num_threads(1)
}
## 1 Node in a SIR model
##
## One individual starts in S, one in I and one in R with a zero
## probability of becoming infected or recover.
##
## At t = 1, one individual exit, with a higher weight for being
## sampled from the R compartment.
model <- SIR(u0 = data.frame(S = 1, I = 1, R = 1),
tspan = 1,
events = data.frame(event = 0,
time = 1,
node = 1,
dest = 0,
n = 1,
proportion = 0,
select = 4,
shift = 0),
beta = 0,
gamma = 0)
## With equal weight the individual is sampled from the S compartment.
set.seed(2)
stopifnot(identical(run(model)@U, structure(c(0L, 1L, 1L), .Dim = c(3L, 1L))))
## With non-equal weight the individual is sampled from the R
## compartment (using the same seed).
model@events@E[3, 4] <- 1000
set.seed(2)
stopifnot(identical(run(model)@U, structure(c(1L, 1L, 0L), .Dim = c(3L, 1L))))
## Check that an error is raised if there are no non-zero elements in
## the selected column in E.
model <- SIR(u0 = data.frame(S = 1, I = 1, R = 1),
tspan = 1,
events = data.frame(event = 0,
time = 1,
node = 1,
dest = 0,
n = 1,
proportion = 0,
select = 4,
shift = 0),
beta = 0,
gamma = 0)
model@events@E <- as(matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0),
nrow = 3, ncol = 4,
dimnames = list(c("S", "I", "R"),
c("1", "2", "3", "4"))),
"dgCMatrix")
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to an enter event
model@events@event <- 1L
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to enter event and n = -1
model@events@n <- -1L
res <- assertError(run(model))
check_error(res, "Unable to sample individuals for event.")
## Change to enter event and proportion = -1
model@events@E <- as(matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1),
nrow = 3, ncol = 4,
dimnames = list(c("S", "I", "R"),
c("1", "2", "3", "4"))),
"dgCMatrix")
model@events@n <- 0L
model@events@proportion <- -1
res <- assertError(.Call(SimInf:::SIR_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
## Change to enter event and proportion = 1.1
model@events@n <- 0L
model@events@proportion <- 1.1
res <- assertError(.Call(SimInf:::SIR_run, model, NULL))
check_error(res, "Invalid proportion detected (< 0.0 or > 1.0).")
## Check that proportion works for an enter event.
model@events@n <- 0L
model@events@proportion <- 0.1
model@events@select <- 1L
set.seed(3)
stopifnot(identical(run(model)@U, structure(c(2L, 1L, 1L), .Dim = c(3L, 1L))))
## Check that an enter event fails with shift without N.
model@events@n <- 1L
model@events@shift <- 1L
model@events@proportion <- 0
res <- assertError(run(model))
check_error(res, "'N' is invalid.")
## Check that shift fails for an enter event when the shift is out of
## bounds.
model@events@N <- matrix(c(3L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
res <- assertError(run(model))
check_error(res, "'shift' is out of bounds.")
model@events@N <- matrix(c(-1L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
res <- assertError(run(model))
check_error(res, "'shift' is out of bounds.")
## Check that shift works for an enter event.
model@events@N <- matrix(c(1L, 0L, 0L),
nrow = 3, ncol = 1,
dimnames = list(c("S", "I", "R"), "1"))
stopifnot(identical(run(model)@U, structure(c(1L, 2L, 1L), .Dim = c(3L, 1L))))
Try the SimInf package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.