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 -- 2022 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) > > ## Create a model with parameters in ldata > model <- mparse(transitions = c("S -> beta*S*I/(S+I+R) -> I + Icum", + "I -> gamma*I -> R"), + compartments = c("S", "I", "Icum", "R"), + ldata = data.frame(beta = 1, gamma = 0.5), + u0 = data.frame(S = 9999, I = 0, Icum = 0, R = 0), + events = data.frame(event = 1, time = 25, node = 1, + dest = 0, n = 1, proportion = 0, + select = 1, shift = 0), + E = matrix(c(0, 1, 0, 0), nrow = 4, ncol = 1, + dimnames = list(c("S", "I", "Icum", "R"), + c("1"))), + tspan = 2:75) > > ## Check that a non-numeric distance vector raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + c("1", "2") + }, + tolerance = c(5, 4))) > check_error(res, "The result from the ABC distance function must be numeric.") > > ## Check that a distance vector with the wrong dimension raises an > ## error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:3 + }, + tolerance = c(5, 4))) > check_error( + res, + "Invalid dimension of the result from the ABC distance function.") > > ## Check that a distance vector with NA raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + c(NA, 2:20) + }, + tolerance = c(5, 4))) > check_error( + res, + "The result from the ABC distance function must be non-negative.") > > ## Check that a distance vector with a negative value raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + c(-1L, 2:20) + }, + tolerance = c(5, 4))) > check_error( + res, + "The result from the ABC distance function must be non-negative.") > > ## Check that a non-numeric tolerance raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:20 + }, + tolerance = c("1", "2"))) > check_error(res, "'tolerance' must have non-negative values.") > > ## Check that a tolerance with NA raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:20 + }, + tolerance = c(NA_real_, 2))) > check_error(res, "'tolerance' must have non-negative values.") > > ## Check that a 'tolerance' with a negative value raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:20 + }, + tolerance = c(1, -2))) > check_error(res, "'tolerance' must have non-negative values.") > > ## Check that a tolerance with wrong dimension raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:20 + }, + tolerance = numeric(0))) > check_error(res, "'tolerance' must have columns.") > > ## Check that a non-decreasing tolerance raises an error. > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 2, + distance = function(result, ...) { + 1:20 + }, + tolerance = c(4, 5))) > check_error(res, "'tolerance' must be a decreasing vector.") > > distance_fn_ldata <- function(result, ...) { + ## Extract the time-series for R1 for each node as a data.frame. + sim <- trajectory(result, "Icum") + + ## Split the 'sim' data.frame by node and calculate the sum of the + ## squared distance at each time-point for every node. + tapply(sim$Icum, sim$node, function(Icum) { + ## Observed cases + cases <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 6, 6, 25, 42, 56, + 106, 171, 279, 382, 576, 710, 977, 934, 846, 672, + 585, 430, 346, 221, 192, 172, 122, 66, 48, 57, 26, + 12, 10, 6, 6, 8, 5, 0, 1, 2, 1, 0, 4, 0, 0, 0, 1, + 0, 0, 0, 0, 0, 0) + + ## Simulated cases + sim_cases <- c(0, diff(c(0, Icum))) + + sum((sim_cases - cases)^2) + }) + } > > ## Check invalid npart > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 1, + distance = distance_fn_ldata, + tolerance = c(250000, 225000))) > check_error(res, "'npart' must be an integer > 1.") > > res <- assertError(abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = c(10, 10), + distance = distance_fn_ldata, + tolerance = c(250000, 225000))) > check_error(res, "'npart' must be an integer > 1.") > > set.seed(123) > fit <- abc(model = model, + priors = c(beta ~ uniform(0.5, 1.5), + gamma ~ uniform(0.3, 0.7)), + npart = 10, + distance = distance_fn_ldata, + tolerance = c(250000, 225000), + verbose = TRUE) Generation 1 ... | | | 0% | |===================== | 30% | |========================================== | 60% | |================================================= | 70% | |=============================================================== | 90% | |======================================================================| 100% accrate = 2.00e-02, ESS = 1.00e+01 time = 2.30 secs Distance -------- Min. 1st Qu. Median Mean 3rd Qu. Max. 1: 60655 87160 135192 128036 163555 180489 Parameters ---------- Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.788 0.843 0.978 1.036 1.246 1.317 gamma 0.324 0.392 0.475 0.518 0.676 0.697 Generation 2 ... | | | 0% | |======================================================================| 100% accrate = 1.72e-01, ESS = 8.85e+00 time = 0.11 secs Distance -------- Min. 1st Qu. Median Mean 3rd Qu. Max. 1: 45569 60730 64953 97559 152094 173509 Parameters ---------- Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.873 0.990 1.063 1.079 1.177 1.292 gamma 0.373 0.449 0.556 0.540 0.623 0.685 > fit Number of particles per generation: 10 Number of generations: 2 Generation 2: ------------- Accrate: 1.72e-01 ESS: 8.85e+00 Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.873 0.990 1.063 1.079 1.177 1.292 gamma 0.373 0.449 0.556 0.540 0.623 0.685 > summary(fit) Number of particles per generation: 10 Number of generations: 2 Generation 1: ------------- Accrate: 2.00e-02 ESS: 1.00e+01 Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.788 0.843 0.978 1.036 1.246 1.317 gamma 0.324 0.392 0.475 0.518 0.676 0.697 Generation 2: ------------- Accrate: 1.72e-01 ESS: 8.85e+00 Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.873 0.990 1.063 1.079 1.177 1.292 gamma 0.373 0.449 0.556 0.540 0.623 0.685 > as.data.frame(fit) generation weight beta gamma 1 1 0.10000000 0.9145463 0.4654897 2 1 0.10000000 0.8012289 0.3242882 3 1 0.10000000 1.2862816 0.6919288 4 1 0.10000000 1.0150297 0.4610293 5 1 0.10000000 0.7882393 0.3682581 6 1 0.10000000 1.3166059 0.6928561 7 1 0.10000000 0.8184946 0.3695935 8 1 0.10000000 1.2520375 0.6275822 9 1 0.10000000 1.2263989 0.6966981 10 1 0.10000000 0.9405331 0.4848395 11 2 0.14180235 1.0383278 0.5849378 12 2 0.05941771 1.2043605 0.6356896 13 2 0.11759919 0.9745226 0.4410641 14 2 0.15419004 0.8729579 0.3728507 15 2 0.06247430 1.2917334 0.6747790 16 2 0.06823140 1.0947209 0.5643436 17 2 0.07022759 1.0885134 0.5479743 18 2 0.06710709 1.2356521 0.6848714 19 2 0.11796070 1.0361072 0.4739716 20 2 0.14098963 0.9493243 0.4162528 > > fit <- continue(fit, tolerance = 200000, verbose = TRUE) Generation 3 ... | | | 0% | |========================================== | 60% | |======================================================================| 100% accrate = 7.25e-02, ESS = 8.61e+00 time = 0.19 secs Distance -------- Min. 1st Qu. Median Mean 3rd Qu. Max. 1: 40414 45744 90633 91654 130570 161747 Parameters ---------- Min. 1st Qu. Median Mean 3rd Qu. Max. beta 0.887 0.995 1.041 1.047 1.114 1.214 gamma 0.398 0.459 0.479 0.506 0.561 0.622 > > pdf_file <- tempfile(fileext = ".pdf") > pdf(pdf_file) > plot(fit, xlim = c(0.3, 1.5), ylim = c(0.3, 1.5)) > dev.off() null device 1 > stopifnot(file.exists(pdf_file)) > unlink(pdf_file) > > ## Check that an invalid 'n' is detected. > sigma <- SimInf:::abc_proposal_covariance(SimInf:::abc_particles(fit, 2L)) > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 0L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "'n' must be an integer > 0.") > > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 2, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "'n' must be an integer > 0.") > > ## Check that an invalid 'parameter' is detected. > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + 3L, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "'parameter' must be a character vector.") > > ## Check that an invalid 'distribution' is detected. > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + "a", + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "Unknown distribution.") > > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + "a", + fit@priors$p1, + fit@priors$p2, + 1L, + NULL, + NULL, + NULL)) > check_error(res, "Unknown distribution.") > > ## Check that an invalid weight is detected. > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + numeric(0), + sigma)) > check_error(res, "'w' must have length >= 1 when 'x' is non-null.") > > fit@weight[2, 2] <- -1 > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "Invalid weight detected (non-finite or < 0.0).") > > fit@weight[2, 2] <- NaN > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "Invalid weight detected (non-finite or < 0.0).") > > fit@weight[2, 2] <- NA_real_ > res <- assertError( + .Call(SimInf:::SimInf_abc_proposals, + fit@priors$parameter, + fit@priors$distribution, + fit@priors$p1, + fit@priors$p2, + 1L, + SimInf:::abc_particles(fit, 2L), + fit@weight[, 2], + sigma)) > check_error(res, "Invalid weight detected (non-finite or < 0.0).") > > proc.time() user system elapsed 2.46 0.15 4.17