R version 4.4.0 RC (2024-04-16 r86468 ucrt) -- "Puppy Cup" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## 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 . > > library(SimInf) > library(tools) > source("util/check.R") > > ## Specify the number of threads to use. > max_threads <- set_num_threads(1) > > ## For debugging > sessionInfo() R version 4.4.0 RC (2024-04-16 r86468 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 [3] LC_MONETARY=C LC_NUMERIC=C [5] LC_TIME=C time zone: Europe/Berlin tzcode source: internal attached base packages: [1] tools stats graphics grDevices utils datasets methods [8] base other attached packages: [1] SimInf_9.7.0 loaded via a namespace (and not attached): [1] MASS_7.3-60.2 compiler_4.4.0 Matrix_1.7-0 grid_4.4.0 digest_0.6.35 [6] lattice_0.22-6 > > ## 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)) Cannot sample 2 for event from 1 individuals. Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 2 proportion: 1 select: 4 shift: 0 > 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.") + } Cannot sample 2 for event from 1 individuals. Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 2 proportion: 1 select: 4 shift: 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, -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)) Cannot sample -1 individuals for event. Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: -1 proportion: 1 select: 4 shift: 0 > 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.") + } Cannot sample -1 individuals for event. Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: -1 proportion: 1 select: 4 shift: 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 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)) Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 0 proportion: 10 select: 4 shift: 0 > 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).") + } Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 0 proportion: 10 select: 4 shift: 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)) Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 0 proportion: -1 select: 4 shift: 0 > 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).") + } Current state in node --------------------- {1, 0, 0, 0, 0, 0} Current state in dest --------------------- {0, 0, 0, 0, 0, 0} Scheduled event --------------- event: 3 (external transfer event) time: 1 node: 2 dest: 1 n: 0 proportion: -1 select: 4 shift: 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)) No states to sample from. Cannot sample 1 for event from 0 individuals. Current state in node --------------------- {1, 1, 1} Scheduled event --------------- event: 0 (exit event) time: 1 node: 1 dest: 0 n: 1 proportion: 0 select: 4 shift: 0 > check_error(res, "Unable to sample individuals for event.") > > ## Change to an enter event > model@events@event <- 1L > res <- assertError(run(model)) No states to sample from. Cannot sample 1 for event from 0 individuals. Current state in node --------------------- {1, 1, 1} Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 1 proportion: 0 select: 4 shift: 0 > check_error(res, "Unable to sample individuals for event.") > > ## Change to enter event and n = -1 > model@events@n <- -1L > res <- assertError(run(model)) No states to sample from. Cannot sample -1 individuals for event. Current state in node --------------------- {1, 1, 1} Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: -1 proportion: 0 select: 4 shift: 0 > 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)) Current state in node --------------------- {1, 1, 1} Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 0 proportion: -1 select: 4 shift: 0 > 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)) Current state in node --------------------- {1, 1, 1} Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 0 proportion: 1.1 select: 4 shift: 0 > 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)) Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 1 proportion: 0 select: 1 shift: 1 > 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)) Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 1 proportion: 0 select: 1 shift: 1 > 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)) Scheduled event --------------- event: 1 (enter event) time: 1 node: 1 dest: 0 n: 1 proportion: 0 select: 1 shift: 1 > 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)))) > > proc.time() user system elapsed 1.10 0.09 1.20