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 > > ## Check measures for a SISe model > model <- SISe(u0 = data.frame(S = 99, I = 1), + tspan = 0:1000, + events = NULL, + phi = 1, + upsilon = 1, + gamma = 0.1, + alpha = 1, + 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(trajectory(model, compartments = "S", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > res <- assertError(trajectory(model, compartments = "I", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > result <- run(model) > result Model: SISe Number of nodes: 1 Number of transitions: 2 Number of scheduled events: 0 Global data ----------- Parameter Value upsilon 1.0 gamma 0.1 alpha 1.0 beta_t1 1.0 beta_t2 1.0 beta_t3 1.0 beta_t4 1.0 epsilon 0.0 Local data ---------- Parameter Value end_t1 91 end_t2 182 end_t3 273 end_t4 365 Continuous state variables -------------------------- Min. 1st Qu. Median Mean 3rd Qu. Max. phi 0.600 0.880 0.900 0.898 0.920 0.990 Compartments ------------ Min. 1st Qu. Median Mean 3rd Qu. Max. S 1.0 8.0 10.0 10.2 12.0 40.0 I 60.0 88.0 90.0 89.8 92.0 99.0 > > res <- assertError(trajectory(result, I ~ S, format = "matrix")) > check_error(res, "Invalid formula specification of 'compartments'.") > > stopifnot(identical( + length(trajectory(result, compartments = "S", format = "matrix")), + 1001L)) > stopifnot(identical( + length(trajectory(result, compartments = "I", format = "matrix")), + 1001L)) > > p <- prevalence(result, I ~ S + I, format = "matrix") > stopifnot(identical(dim(p), c(1L, 1001L))) > > p <- prevalence(result, I ~ S + I, level = 3, format = "matrix") > stopifnot(identical(dim(p), c(1L, 1001L))) > > if (SimInf:::have_openmp() && max_threads > 1) { + stopifnot(identical(set_num_threads(2), 1L)) + result <- run(model) + set_num_threads(1) + result + + stopifnot(identical( + length(trajectory(result, compartments = "S", format = "matrix")), + 1001L)) + stopifnot(identical( + length(trajectory(result, compartments = "I", format = "matrix")), + 1001L)) + + p <- prevalence(result, I ~ S + I, format = "matrix") + stopifnot(identical(dim(p), c(1L, 1001L))) + + p <- prevalence(result, I ~ S + I, level = 3, format = "matrix") + stopifnot(identical(dim(p), c(1L, 1001L))) + } > > ## Check measures for a SISe_sp model > model <- SISe_sp(u0 = data.frame(S = 99, I = 1), + tspan = 0:1000, + events = NULL, + phi = 1, + upsilon = 1, + gamma = 0.1, + alpha = 1, + beta_t1 = 1, + beta_t2 = 1, + beta_t3 = 1, + beta_t4 = 1, + end_t1 = 91, + end_t2 = 182, + end_t3 = 273, + end_t4 = 365, + coupling = 0, + distance = distance_matrix(1, 1, 1)) > > res <- assertError(trajectory(model, compartments = "S", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > res <- assertError(trajectory(model, compartments = "I", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > result <- run(model) > result Model: SISe_sp Number of nodes: 1 Number of transitions: 2 Number of scheduled events: 0 Global data ----------- Parameter Value upsilon 1.0 gamma 0.1 alpha 1.0 beta_t1 1.0 beta_t2 1.0 beta_t3 1.0 beta_t4 1.0 coupling 0.0 Continuous state variables -------------------------- Min. 1st Qu. Median Mean 3rd Qu. Max. phi 0.530 0.870 0.900 0.895 0.920 0.990 Compartments ------------ Min. 1st Qu. Median Mean 3rd Qu. Max. S 1.0 8.0 10.0 10.5 13.0 47.0 I 53.0 87.0 90.0 89.5 92.0 99.0 > > stopifnot(identical( + length(trajectory(result, compartments = "S", format = "matrix")), + 1001L)) > stopifnot(identical( + length(trajectory(result, compartments = "I", format = "matrix")), + 1001L)) > > p <- prevalence(result, I ~ S + I, format = "matrix") > stopifnot(identical(dim(p), c(1L, 1001L))) > > p <- prevalence(result, I ~ S + I, level = 3, format = "matrix") > stopifnot(identical(dim(p), c(1L, 1001L))) > > if (SimInf:::have_openmp() && max_threads > 1) { + set_num_threads(2) + result <- run(model) + set_num_threads(1) + result + + stopifnot(identical( + length(trajectory(result, compartments = "S", format = "matrix")), + 1001L)) + stopifnot(identical( + length(trajectory(result, compartments = "I", format = "matrix")), + 1001L)) + + p <- prevalence(result, I ~ S + I, format = "matrix") + stopifnot(identical(dim(p), c(1L, 1001L))) + + p <- prevalence(result, I ~ S + I, level = 3, format = "matrix") + stopifnot(identical(dim(p), c(1L, 1001L))) + } > > ## Check 'susceptible' and 'infected' methods for a SISe3 model > u0 <- data.frame(S_1 = rep(10, 10), I_1 = rep(0, 10), + S_2 = rep(20, 10), I_2 = rep(0, 10), + S_3 = rep(70, 10), I_3 = rep(0, 10)) > model <- SISe3(u0 = u0, + tspan = seq_len(1000) - 1, + events = NULL, + phi = rep(1, 10), + upsilon_1 = 0.0357, + upsilon_2 = 0.0357, + upsilon_3 = 0.00935, + gamma_1 = 0.1, + gamma_2 = 0.1, + gamma_3 = 0.1, + alpha = 1.0, + beta_t1 = 0.19, + beta_t2 = 0.085, + beta_t3 = 0.075, + beta_t4 = 0.185, + end_t1 = 91, + end_t2 = 182, + end_t3 = 273, + end_t4 = 365, + epsilon = 0.000011) > > res <- assertError(trajectory(model, compartments = "S_1", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > res <- assertError(trajectory(model, compartments = "I_1", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > result <- run(model) > result Model: SISe3 Number of nodes: 10 Number of transitions: 6 Number of scheduled events: 0 Global data ----------- Parameter Value upsilon_1 0.035700 upsilon_2 0.035700 upsilon_3 0.009350 gamma_1 0.100000 gamma_2 0.100000 gamma_3 0.100000 alpha 1.000000 beta_t1 0.190000 beta_t2 0.085000 beta_t3 0.075000 beta_t4 0.185000 epsilon 0.000011 Local data ---------- Parameter Value end_t1 91 end_t2 182 end_t3 273 end_t4 365 Continuous state variables -------------------------- Min. 1st Qu. Median Mean 3rd Qu. Max. phi 5.79e-05 6.17e-05 1.35e-04 3.54e-01 1.86e-02 6.13e+00 Compartments ------------ Min. 1st Qu. Median Mean 3rd Qu. Max. S_1 0.00 10.00 10.00 9.37 10.00 10.00 I_1 0.00 0.00 0.00 0.63 0.00 10.00 S_2 1.00 20.00 20.00 18.74 20.00 20.00 I_2 0.00 0.00 0.00 1.26 0.00 19.00 S_3 39.00 70.00 70.00 68.29 70.00 70.00 I_3 0.00 0.00 0.00 1.71 0.00 31.00 > > stopifnot(identical( + length(trajectory(result, compartments = "S_1", format = "matrix")), + 10000L)) > stopifnot(identical( + length(trajectory(result, compartments = "I_1", format = "matrix")), + 10000L)) > > p <- prevalence(result, I_1 + I_2 + I_3 ~ ., format = "matrix") > stopifnot(identical(dim(p), c(1L, 1000L))) > > p <- prevalence(result, I_1 + I_2 + I_3 ~ ., level = 3, format = "matrix") > stopifnot(identical(dim(p), c(10L, 1000L))) > > if (SimInf:::have_openmp() && max_threads > 1) { + set_num_threads(2) + result <- run(model) + set_num_threads(1) + result + + stopifnot(identical( + length(trajectory(result, compartments = "S_1", format = "matrix")), + 10000L)) + stopifnot(identical( + length(trajectory(result, compartments = "I_1", format = "matrix")), + 10000L)) + + p <- prevalence(result, I_1 + I_2 + I_3 ~ ., format = "matrix") + stopifnot(identical(dim(p), c(1L, 1000L))) + + p <- prevalence(result, I_1 + I_2 + I_3 ~ ., + level = 3, format = "matrix") + stopifnot(identical(dim(p), c(10L, 1000L))) + } > > ## Check measures with a SISe3_sp model > u0 <- data.frame(S_1 = rep(10, 10), I_1 = rep(0, 10), + S_2 = rep(20, 10), I_2 = rep(0, 10), + S_3 = rep(70, 10), I_3 = rep(0, 10)) > > model <- SISe3_sp(u0 = u0, + tspan = seq_len(1000) - 1, + events = NULL, + phi = rep(1, 10), + upsilon_1 = 0.0357, + upsilon_2 = 0.0357, + upsilon_3 = 0.00935, + gamma_1 = 0.1, + gamma_2 = 0.1, + gamma_3 = 0.1, + alpha = 1.0, + beta_t1 = 0.19, + beta_t2 = 0.085, + beta_t3 = 0.075, + beta_t4 = 0.185, + end_t1 = 91, + end_t2 = 182, + end_t3 = 273, + end_t4 = 365, + coupling = 0, + distance = distance_matrix(1:10, 1:10, 1)) > > res <- assertError(trajectory(model, compartments = "S_1", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > res <- assertError(trajectory(model, compartments = "I_1", format = "matrix")) > check_error(res, "Please run the model first, the trajectory is empty.") > > result <- run(model) > result Model: SISe3_sp Number of nodes: 10 Number of transitions: 6 Number of scheduled events: 0 Global data ----------- Parameter Value upsilon_1 0.03570 upsilon_2 0.03570 upsilon_3 0.00935 gamma_1 0.10000 gamma_2 0.10000 gamma_3 0.10000 alpha 1.00000 beta_t1 0.19000 beta_t2 0.08500 beta_t3 0.07500 beta_t4 0.18500 coupling 0.00000 Continuous state variables -------------------------- Min. 1st Qu. Median Mean 3rd Qu. Max. phi 4.83e-59 2.99e-25 4.03e-04 7.87e-01 5.78e-01 6.69e+00 Compartments ------------ Min. 1st Qu. Median Mean 3rd Qu. Max. S_1 0.00 8.00 10.00 8.64 10.00 10.00 I_1 0.00 0.00 0.00 1.36 2.00 10.00 S_2 1.00 16.00 20.00 17.39 20.00 20.00 I_2 0.00 0.00 0.00 2.61 4.00 19.00 S_3 32.00 66.00 70.00 66.12 70.00 70.00 I_3 0.00 0.00 0.00 3.88 4.00 38.00 > > stopifnot(identical( + length(trajectory(result, compartments = "S_1", format = "matrix")), + 10000L)) > stopifnot(identical( + length(trajectory(result, compartments = "I_1", format = "matrix")), + 10000L)) > > p <- prevalence(result, I_1 + I_2 + I_3 ~ ., format = "matrix") > stopifnot(identical(dim(p), c(1L, 1000L))) > > p <- prevalence(result, I_1 + I_2 + I_3 ~ ., level = 3, format = "matrix") > stopifnot(identical(dim(p), c(10L, 1000L))) > > if (SimInf:::have_openmp() && max_threads > 1) { + set_num_threads(2) + result <- run(model) + set_num_threads(1) + result + + stopifnot(identical( + length(trajectory(result, compartments = "S_1", format = "matrix")), + 10000L)) + stopifnot(identical( + length(trajectory(result, compartments = "I_1", format = "matrix")), + 10000L)) + + p <- prevalence(result, I_1 + I_2 + I_3 ~ ., format = "matrix") + stopifnot(identical(dim(p), c(1L, 1000L))) + + p <- prevalence(result, I_1 + I_2 + I_3 ~ ., + level = 3, format = "matrix") + stopifnot(identical(dim(p), c(10L, 1000L))) + } > > proc.time() user system elapsed 1.46 0.23 1.73