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 -- 2020 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. > 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 invalid u0 > res <- assertError(SEIR(u0 = "u0")) > check_error(res, "Missing columns in u0.") > > u0 <- data.frame(S = c(0, 1, 2, 3, 4, 5), + E = c(0, 0, 0, 0, 0, 0), + I = c(0, 0, 0, 0, 0, 0), + R = c(0, 0, 0, 0, 0, 0)) > > ## Check missing columns in u0 > res <- assertError(SEIR(u0 = u0[, c("E", "I", "R"), drop = FALSE])) > check_error(res, "Missing columns in u0.") > > res <- assertError(SEIR(u0 = u0[, c("S", "I", "R"), drop = FALSE])) > check_error(res, "Missing columns in u0.") > > res <- assertError(SEIR(u0 = u0[, c("S", "E", "R"), drop = FALSE])) > check_error(res, "Missing columns in u0.") > > res <- assertError(SEIR(u0 = u0[, c("S", "E", "I"), drop = FALSE])) > check_error(res, "Missing columns in u0.") > > ## Check missing beta > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + epsilon = 0.5, + gamma = 0.5)) > check_error(res, "'beta' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check missing epsilon > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + gamma = 0.5)) > check_error(res, "'epsilon' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check missing gamma > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + epsilon = 0.5)) > check_error(res, "'gamma' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check non-numeric beta > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = "0.5", + epsilon = 0.3, + gamma = 0.1)) > check_error(res, "'beta' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check non-numeric epsilon > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + epsilon = "0.3", + gamma = 0.1)) > check_error(res, "'epsilon' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check non-numeric gamma > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + epsilon = 0.3, + gamma = "0.1")) > check_error(res, "'gamma' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check that length of beta equals 1 > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = c(0.5, 0.5), + epsilon = 0.3, + gamma = 0.1)) > check_error(res, "'beta' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check that length of epsilon equals 1 > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + epsilon = c(0.3, 0.3), + gamma = 0.1)) > check_error(res, "'epsilon' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check that length of gamma equals 1 > res <- assertError(SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0.5, + epsilon = 0.3, + gamma = c(0.1, 0.1))) > check_error(res, "'gamma' must be numeric of length 1 or 'nrow(u0)'.") > > ## Check extraction of data from 'suscpetible', 'infected' and > ## 'recovered' compartments > model <- SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0, + epsilon = 0, + gamma = 0) > > result <- run(model) > > S_expected <- structure(c(0L, 1L, 2L, 3L, 4L, 5L, 0L, 1L, 2L, 3L, 4L, 5L, 0L, + 1L, 2L, 3L, 4L, 5L, 0L, 1L, 2L, 3L, 4L, 5L, 0L, 1L, + 2L, 3L, 4L, 5L, 0L, 1L, 2L, 3L, 4L, 5L, 0L, 1L, 2L, + 3L, 4L, 5L, 0L, 1L, 2L, 3L, 4L, 5L, 0L, 1L, 2L, 3L, + 4L, 5L, 0L, 1L, 2L, 3L, 4L, 5L), + .Dim = c(6L, 10L)) > > S_observed <- trajectory(result, compartments = "S", format = "matrix") > stopifnot(identical(S_observed, S_expected)) > > I_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), + .Dim = c(6L, 10L)) > > I_observed <- trajectory(result, compartments = "I", format = "matrix") > stopifnot(identical(I_observed, I_expected)) > > R_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), + .Dim = c(6L, 10L)) > > R_observed <- trajectory(result, compartments = "R", format = "matrix") > stopifnot(identical(R_observed, R_expected)) > > R_expected <- data.frame( + node = c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, + 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, + 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, + 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, + 5L, 6L), + time = c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, + 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, + 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, + 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, + 9L), + R = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)) > R_observed <- trajectory(result, compartments = "R") > stopifnot(identical(R_observed, R_expected)) > > ## Extract the number of recovered individuals in the first node after > ## each time step in the simulation > R_expected <- data.frame( + node = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), + time = 0:9L, + R = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)) > R_observed <- trajectory(result, compartments = "R", index = 1) > stopifnot(identical(R_observed, R_expected)) > > R_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), + .Dim = c(1L, 10L)) > R_observed <- trajectory(result, compartments = "R", + index = 1, format = "matrix") > stopifnot(identical(R_observed, R_expected)) > > ## Extract the number of recovered individuals in the first and third > ## node after each time step in the simulation > R_expected <- data.frame( + node = c(1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, + 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L), + time = c(0L, 0L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, + 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L), + R = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)) > R_observed <- trajectory(result, compartments = "R", index = c(1, 3)) > stopifnot(identical(R_observed, R_expected)) > > R_expected <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), + .Dim = c(2L, 10L)) > R_observed <- trajectory(result, compartments = "R", index = c(1, 3), + format = "matrix") > stopifnot(identical(R_observed, R_expected)) > > ## A more complex test to extract data from U from a trajectory of 6 > ## nodes and 10 time steps (0:9). Generate events that adds one > ## individual per node in each compartment each day. To achieve this, > ## using the template SEIR model, first add them to the 'S' > ## compartment, and then modify the select matrix E and update the > ## select index. > events <- do.call("rbind", lapply(seq_len(9), function(time) { + do.call("rbind", lapply(seq_len(6), function(node) { + data.frame(event = 1, time = time, node = node, dest = 0, + n = 1, proportion = 0, select = rep(1, 4), shift = 0) + })) + })) > model <- SEIR(u0 = data.frame(S = c(110, 210, 310, 410, 510, 610), + E = c(120, 220, 320, 420, 520, 620), + I = c(130, 230, 330, 430, 530, 630), + R = c(140, 240, 340, 440, 540, 640)), + tspan = 0:9, events = events, beta = 0, + epsilon = 0, gamma = 0) > model@events@E <- as(diag(4), "dgCMatrix") > model@events@select <- rep(1:4, length.out = length(model@events@select)) > > # Check that this fails because rownames (compartments) are missing > res <- assertError(run(model)) > check_error( + res, + "'S' and 'E' must have rownames matching the compartments.", + FALSE) > > rownames(model@events@E) <- c("S", "E", "I", "R") > result <- run(model) > > U_expected <- data.frame( + node = c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, + 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, + 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, + 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, + 5L, 6L), + time = c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, + 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, + 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, + 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, + 9L), + S = c(110L, 210L, 310L, 410L, 510L, 610L, 111L, 211L, 311L, + 411L, 511L, 611L, 112L, 212L, 312L, 412L, 512L, 612L, 113L, 213L, + 313L, 413L, 513L, 613L, 114L, 214L, 314L, 414L, 514L, 614L, 115L, + 215L, 315L, 415L, 515L, 615L, 116L, 216L, 316L, 416L, 516L, 616L, + 117L, 217L, 317L, 417L, 517L, 617L, 118L, 218L, 318L, 418L, 518L, + 618L, 119L, 219L, 319L, 419L, 519L, 619L), + E = c(120L, 220L, + 320L, 420L, 520L, 620L, 121L, 221L, 321L, 421L, 521L, 621L, 122L, + 222L, 322L, 422L, 522L, 622L, 123L, 223L, 323L, 423L, 523L, 623L, + 124L, 224L, 324L, 424L, 524L, 624L, 125L, 225L, 325L, 425L, 525L, + 625L, 126L, 226L, 326L, 426L, 526L, 626L, 127L, 227L, 327L, 427L, + 527L, 627L, 128L, 228L, 328L, 428L, 528L, 628L, 129L, 229L, 329L, + 429L, 529L, 629L), + I = c(130L, 230L, 330L, 430L, 530L, 630L, + 131L, 231L, 331L, 431L, 531L, 631L, 132L, 232L, 332L, 432L, 532L, + 632L, 133L, 233L, 333L, 433L, 533L, 633L, 134L, 234L, 334L, 434L, + 534L, 634L, 135L, 235L, 335L, 435L, 535L, 635L, 136L, 236L, 336L, + 436L, 536L, 636L, 137L, 237L, 337L, 437L, 537L, 637L, 138L, 238L, + 338L, 438L, 538L, 638L, 139L, 239L, 339L, 439L, 539L, 639L), + R = c(140L, 240L, 340L, 440L, 540L, 640L, 141L, 241L, 341L, + 441L, 541L, 641L, 142L, 242L, 342L, 442L, 542L, 642L, 143L, + 243L, 343L, 443L, 543L, 643L, 144L, 244L, 344L, 444L, 544L, + 644L, 145L, 245L, 345L, 445L, 545L, 645L, 146L, 246L, 346L, + 446L, 546L, 646L, 147L, 247L, 347L, 447L, 547L, 647L, 148L, + 248L, 348L, 448L, 548L, 648L, 149L, 249L, 349L, 449L, 549L, + 649L)) > > U_observed <- trajectory(result) > stopifnot(identical(U_observed, U_expected)) > U_observed <- trajectory(result, compartments = c("S", "E", "I", "R"), + index = 1:6) > stopifnot(identical(U_observed, U_expected)) > > U_expected <- data.frame( + node = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), + time = 0:9, S = 110:119, I = 130:139) > U_observed <- trajectory(result, compartments = c("S", "I"), index = 1) > stopifnot(identical(U_observed, U_expected)) > > U_expected <- data.frame( + node = c(2L, 4L, 2L, 4L, 2L, 4L, 2L, 4L, 2L, 4L, + 2L, 4L, 2L, 4L, 2L, 4L, 2L, 4L, 2L, 4L), + time = c(0L, 0L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, + 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L), + E = c(220L, 420L, 221L, 421L, 222L, 422L, 223L, 423L, 224L, 424L, + 225L, 425L, 226L, 426L, 227L, 427L, 228L, 428L, 229L, 429L), + R = c(240L, 440L, 241L, 441L, 242L, 442L, 243L, 443L, 244L, 444L, + 245L, 445L, 246L, 446L, 247L, 447L, 248L, 448L, 249L, 449L)) > U_observed <- trajectory(result, compartments = c("E", "R"), index = c(2, 4)) > stopifnot(identical(U_observed, U_expected)) > > U_expected <- structure( + c(220L, 240L, 420L, 440L, 221L, 241L, 421L, 441L, 222L, + 242L, 422L, 442L, 223L, 243L, 423L, 443L, 224L, 244L, 424L, 444L, + 225L, 245L, 425L, 445L, 226L, 246L, 426L, 446L, 227L, 247L, 427L, + 447L, 228L, 248L, 428L, 448L, 229L, 249L, 429L, 449L), + .Dim = c(4L, 10L)) > U_observed <- trajectory(result, compartments = c("E", "R"), + index = c(2, 4), format = "matrix") > stopifnot(identical(U_observed, U_expected)) > > ## Check prevalence > > ## Define a tolerance > tol <- 1e-8 > > p_expected <- data.frame( + time = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), + prevalence = c(0.253333333333333, 0.253324468085106, 0.253315649867374, + 0.253306878306878, 0.253298153034301, 0.253289473684211, + 0.253280839895013, 0.253272251308901, 0.253263707571802, + 0.253255208333333)) > p_observed <- prevalence(result, I ~ S + E + I + R) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > p_observed <- prevalence(result, I ~ .) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > > p_expected <- structure( + list(time = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), + prevalence = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), + .Names = c("time", "prevalence"), + row.names = c(NA, -10L), + class = "data.frame") > p_observed <- prevalence(result, I ~ S + E + I + R, level = 2) > stopifnot(identical(p_observed, p_expected)) > p_observed <- prevalence(result, I ~ ., level = 2) > stopifnot(identical(p_observed, p_expected)) > > p_expected <- data.frame( + node = c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, + 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, + 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, + 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L), + time = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, + 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, + 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9), + prevalence = c(0.26, 0.255555555555556, 0.253846153846154, + 0.252941176470588, 0.252380952380952, 0.252, + 0.259920634920635, 0.255530973451327, 0.253834355828221, + 0.252934272300469, 0.252376425855513, 0.251996805111821, + 0.259842519685039, 0.255506607929515, 0.253822629969419, + 0.252927400468384, 0.252371916508539, 0.251993620414673, + 0.259765625, 0.255482456140351, 0.253810975609756, + 0.252920560747664, 0.252367424242424, 0.251990445859873, + 0.25968992248062, 0.255458515283843, 0.253799392097264, + 0.252913752913753, 0.252362948960302, 0.251987281399046, + 0.259615384615385, 0.255434782608696, 0.253787878787879, + 0.252906976744186, 0.252358490566038, 0.251984126984127, + 0.259541984732824, 0.255411255411255, 0.253776435045317, + 0.252900232018561, 0.252354048964218, 0.251980982567353, + 0.259469696969697, 0.255387931034483, 0.253765060240964, + 0.252893518518519, 0.25234962406015, 0.251977848101266, + 0.259398496240602, 0.255364806866953, 0.253753753753754, + 0.252886836027714, 0.25234521575985, 0.251974723538705, + 0.259328358208955, 0.25534188034188, 0.25374251497006, + 0.252880184331797, 0.252340823970037, 0.251971608832808)) > p_observed <- prevalence(result, I ~ ., level = 3) > stopifnot(identical(p_observed$node, p_expected$node)) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > p_observed <- prevalence(result, I ~ S + E + I + R, level = 3) > stopifnot(identical(p_observed$node, p_expected$node)) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > > p_expected <- data.frame( + node = c(2L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, + 2L, 3L, 2L, 3L), + time = c(0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9), + prevalence = c(0.255555555555556, 0.253846153846154, 0.255530973451327, + 0.253834355828221, 0.255506607929515, 0.253822629969419, + 0.255482456140351, 0.253810975609756, 0.255458515283843, + 0.253799392097264, 0.255434782608696, 0.253787878787879, + 0.255411255411255, 0.253776435045317, 0.255387931034483, + 0.253765060240964, 0.255364806866953, 0.253753753753754, + 0.25534188034188, 0.25374251497006)) > p_observed <- prevalence(result, I ~ ., level = 3, i = 2:3) > stopifnot(identical(p_observed$node, p_expected$node)) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > p_observed <- prevalence(result, I ~ S + E + I + R, level = 3, i = 2:3) > stopifnot(identical(p_observed$node, p_expected$node)) > stopifnot(identical(p_observed$time, p_expected$time)) > stopifnot(all(abs(p_observed$prevalence - p_expected$prevalence) < tol)) > > ## Check 'V' > res <- assertError(trajectory(result, "phi")) > check_error(res, "Non-existing compartment(s) in model: 'phi'.") > > ## Check data > stopifnot(identical(nrow(events_SEIR()), 466692L)) > stopifnot(identical(nrow(u0_SEIR()), 1600L)) > > ## Try to plot non-extisting compartment. > res <- assertError(plot(result, "X")) > check_error(res, "Non-existing compartment(s) in model: 'X'.") > > ## Try to plot with invalid range argument. > res <- assertError(plot(result, range = 1.2)) > check_error(res, "'range' must be FALSE or a value between 0 and 1.") > > ## Check SEIR plot method > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > plot(result) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check SEIR plot method with range = FALSE > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > plot(result, "S", lty = 1, range = FALSE) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check SEIR boxplot method > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > boxplot(result) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check SEIR pairs plot method > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > pairs(result) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check SEIR events plot with no events > model <- SEIR(u0 = u0, + tspan = seq_len(10) - 1, + events = NULL, + beta = 0, + epsilon = 0, + gamma = 0) > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > plot(model@events) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check SEIR events plot method > model <- SEIR(u0 = u0_SEIR(), + tspan = seq_len(365 * 4), + events = events_SEIR(), + beta = 0, + epsilon = 0, + gamma = 0) > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > plot(model@events) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check that C SEIR run function fails for a misspecified SEIR model > res <- assertError(.Call(SimInf:::SEIR_run, NULL, NULL)) > check_error(res, "Invalid model.") > > res <- assertError(.Call(SimInf:::SEIR_run, "SEIR", NULL)) > check_error(res, "Invalid model.") > > ## Check that an invalid rate error is raised during the simulation. > model <- SEIR(u0 = data.frame(S = rep(100, 2), E = 0, I = 10, R = 0), + tspan = 1:100, + beta = 0.16, + epsilon = -0.3, + gamma = 0.077) > set.seed(1) > res <- assertError(run(model, solver = "ssm")) Status: ------- Time: 1.34252 Node: 1 Current state in node: u(length: 4) = {99, 1, 10, 0} v(length: 0) = {} ldata(length: 3) = {0.16, -0.3, 0.077} Transition: 2 Rate: -0.3 > check_error(res, "Invalid rate detected (non-finite or < 0.0).") > res <- assertError(run(model, solver = "aem")) Status: ------- Time: 1.36216 Node: 1 Current state in node: u(length: 4) = {99, 1, 10, 0} v(length: 0) = {} ldata(length: 3) = {0.16, -0.3, 0.077} Transition: 2 Rate: -0.3 Status: ------- Time: 1.59832 Node: 1 Current state in node: u(length: 4) = {98, 1, 11, 0} v(length: 0) = {} ldata(length: 3) = {0.16, -0.3, 0.077} Transition: 2 Rate: -0.3 > check_error(res, "Invalid rate detected (non-finite or < 0.0).") > > proc.time() user system elapsed 2.03 0.37 2.39