diff --git a/.github/workflows/coverage.yaml b/.github/workflows/coverage.yaml index 8b99c08e..6fc769e3 100644 --- a/.github/workflows/coverage.yaml +++ b/.github/workflows/coverage.yaml @@ -15,7 +15,7 @@ jobs: - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 - + - uses: r-lib/actions/setup-r@v2.9.0 - uses: r-lib/actions/setup-r-dependencies@v2 @@ -36,4 +36,3 @@ jobs: shell: Rscript {0} env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} - diff --git a/.github/workflows/pre-commit.yaml b/.github/workflows/pre-commit.yaml new file mode 100644 index 00000000..2efb1421 --- /dev/null +++ b/.github/workflows/pre-commit.yaml @@ -0,0 +1,23 @@ +name: Pre-commit + +on: + push: + branches: + - main + pull_request: + branches: + - main + +jobs: + pre-commit: + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v3 + + - name: Run pre-commit hooks + uses: pre-commit/action@v3.0.1 diff --git a/.github/workflows/r.yml b/.github/workflows/r.yml index 72f9c523..bf69e222 100644 --- a/.github/workflows/r.yml +++ b/.github/workflows/r.yml @@ -12,7 +12,7 @@ on: name: R-CMD-check-final -env: +env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} jobs: @@ -82,7 +82,7 @@ jobs: name: epiworldR-built-package-${{ matrix.config.os }}-${{ matrix.config.r }} path: epiworldR_*.tar.gz retention-days: 7 - + epiworldShiny: runs-on: ubuntu-latest container: rocker/tidyverse:4.4.0 @@ -98,7 +98,7 @@ jobs: run: | install2.r -n 2 shinyjs shinydashboard DT shinycssloaders plotly installGithub.r UofUEpiBio/epiworldR@${{ github.sha }} - + - name: Check the package run: | R CMD build . @@ -118,11 +118,7 @@ jobs: run: | install2.r -n 2 shinyjs shinydashboard DT shinycssloaders plotly installGithub.r UofUEpiBio/epiworldR@${{ github.sha }} - + - name: Check the package run: | R CMD check --no-manual epiworldRShiny_*tar.gz - - - - diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..93c401eb --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,31 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v3.2.0 + hooks: + - id: trailing-whitespace + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + - id: end-of-file-fixer + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + - id: check-yaml + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + - id: check-added-large-files + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + - id: detect-private-key + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + +- repo: https://github.com/lorenzwalthert/precommit + rev: v0.4.2 + hooks: + - id: style-files + args: + [ + '--ignore-start="^# styler: off$"', + '--ignore-stop="^# styler: on$"', + '--strict=FALSE' + ] + exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) + - id: readme-rmd-rendered + # - id: lintr + # args: [--warn_only] + # verbose: true + # exclude: ^(inst/|man/|src/cpp11.cpp|playgroud/|R/cpp11.R) diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index eb3d24ec..2b240b0c 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -17,4 +17,4 @@ } ], "version": 4 -} \ No newline at end of file +} diff --git a/.vscode/settings.json b/.vscode/settings.json index a09d2f0d..b4549178 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -62,4 +62,4 @@ "editor.insertSpaces": true, "editor.detectIndentation": false } -} \ No newline at end of file +} diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md index 0ed3b818..a8187564 100644 --- a/CODE_OF_CONDUCT.md +++ b/CODE_OF_CONDUCT.md @@ -59,7 +59,7 @@ representative at an online or offline event. ## Enforcement Instances of abusive, harassing, or otherwise unacceptable behavior may be -reported to the community leaders responsible for enforcement at derekmeyer37@gmail.com. +reported to the community leaders responsible for enforcement at derekmeyer37@gmail.com. All complaints will be reviewed and investigated promptly and fairly. All community leaders are obligated to respect the privacy and security of the diff --git a/DESCRIPTION b/DESCRIPTION index d14e78bd..e98c661e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -5,7 +5,7 @@ Version: 0.3-2 Authors@R: c( person(given="George", family="Vega Yon", role=c("aut","cre"), email="g.vegayon@gmail.com", comment = c(ORCID = "0000-0002-3171-0844")), - person(given="Derek", family="Meyer", role=c("aut"), + person(given="Derek", family="Meyer", role=c("aut"), email="derekmeyer37@gmail.com", comment = c(ORCID = "0009-0005-1350-6988")), person(given="Susan", family="Holmes", role = "rev", comment = c(what = "JOSS reviewer", ORCID="0000-0002-2208-8168")), @@ -17,20 +17,20 @@ Authors@R: c( Description: A flexible framework for Agent-Based Models (ABM), the 'epiworldR' package provides methods for prototyping disease outbreaks and transmission models using a 'C++' backend, making it very fast. It supports multiple epidemiological models, including the Susceptible-Infected-Susceptible (SIS), Susceptible-Infected-Removed (SIR), Susceptible-Exposed-Infected-Removed (SEIR), and others, involving arbitrary mitigation policies and multiple-disease models. Users can specify infectiousness/susceptibility rates as a function of agents' features, providing great complexity for the model dynamics. Furthermore, 'epiworldR' is ideal for simulation studies featuring large populations. URL: https://github.com/UofUEpiBio/epiworldR, https://uofuepibio.github.io/epiworldR/, - https://uofuepibio.github.io/epiworldR-workshop/ + https://uofuepibio.github.io/epiworldR-workshop/ BugReports: https://github.com/UofUEpiBio/epiworldR/issues License: MIT + file LICENSE RoxygenNote: 7.3.2 Encoding: UTF-8 Roxygen: list(markdown = TRUE) LinkingTo: cpp11 -Suggests: +Suggests: knitr, rmarkdown, tinytest, netplot, igraph, data.table -Imports: +Imports: utils VignetteBuilder: knitr diff --git a/Makefile b/Makefile index 7effd62f..7e0f0bbd 100644 --- a/Makefile +++ b/Makefile @@ -2,11 +2,11 @@ VERSION := $(shell grep Version DESCRIPTION | sed -e 's/Version: //') -build: +build: cd .. && R CMD build epiworldR debug: clean - docker run --rm -ti -w/mnt -v $(PWD):/mnt uofuepibio/epiworldr:debug make docker-debug + docker run --rm -ti -w/mnt -v $(PWD):/mnt uofuepibio/epiworldr:debug make docker-debug docker-debug: EPI_CONFIG="-DEPI_DEBUG -Wall -pedantic -g" R CMD INSTALL \ @@ -15,33 +15,33 @@ docker-debug: install-dev: clean sed -i -E 's/@useDynLib\s+[a-zA-Z]+/@useDynLib epiworldRdev/g' R/epiworldR-package.R sed -i -E 's/useDynLib\(+[a-zA-Z]+/useDynLib(epiworldRdev/g' NAMESPACE - sed -i -E 's/^Package:.+/Package: epiworldRdev/g' DESCRIPTION + sed -i -E 's/^Package:.+/Package: epiworldRdev/g' DESCRIPTION sed -i -E 's/^library\(epiworldR\)/library(epiworldRdev)/g' README.* Rscript --vanilla -e 'roxygen2::roxygenize()' EPI_DEV=yes R CMD INSTALL .& $(MAKE) clean -install: +install: cd .. && \ R CMD INSTALL epiworldR_$(VERSION).tar.gz - + README.md: README.Rmd Rscript --vanilla -e 'rmarkdown::render("README.Rmd")' # update: # wget https://raw.githubusercontent.com/UofUEpiBio/epiworld/master/epiworld.hpp && \ -# mv epiworld.hpp inst/include/epiworld.hpp +# mv epiworld.hpp inst/include/epiworld.hpp local-update: rsync -avz ../epiworld/include/epiworld inst/include/. check: build cd .. && R CMD check epiworldR_*.tar.gz -clean: +clean: rm -f src/*.dll src/*.so src/*.o sed -i -E 's/@useDynLib\s+[a-zA-Z]+/@useDynLib epiworldR/g' R/epiworldR-package.R sed -i -E 's/useDynLib\(+[a-zA-Z]+/useDynLib(epiworldR/g' NAMESPACE - sed -i -E 's/^Package:.+/Package: epiworldR/g' DESCRIPTION + sed -i -E 's/^Package:.+/Package: epiworldR/g' DESCRIPTION # sed -i -E 's/^\\(name|alias|title)\{[a-zA-Z]+/\\\1{epiworldR/g' man/epiworldR-package.Rd sed -i -E 's/^library\(epiworldRdev\)/library(epiworldR)/g' README.* @@ -52,4 +52,3 @@ docs: checkv: build R CMD check --as-cran --use-valgrind epiworldR*.tar.gz - diff --git a/NEWS.md b/NEWS.md index 344e6951..d0be8221 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ * Ports the `Entity` class. Entities are used to group agents within a model. * Refactors `add_tool`, `add_virus`, and `add_entity` simplifying syntax. Now, - these functions only receive the model and object. Prevalence is + these functions only receive the model and object. Prevalence is specified in the object itself. `add_tool_n` and `add_virus_n` are now deprecated. diff --git a/R/ModelDiffNet.R b/R/ModelDiffNet.R index ac96e68e..3a7faee2 100644 --- a/R/ModelDiffNet.R +++ b/R/ModelDiffNet.R @@ -1,10 +1,10 @@ #' Network Diffusion Model -#' -#' The network diffusion model is a simple model that assumes that -#' the probability of adoption of a behavior is proportional to the +#' +#' The network diffusion model is a simple model that assumes that +#' the probability of adoption of a behavior is proportional to the #' number of adopters in the network. -#' -#' @export +#' +#' @export #' @param name Name of the model. #' @param prevalence Prevalence of the disease. #' @param prob_adopt Probability of adoption. @@ -22,22 +22,22 @@ #' \deqn{ #' P(adopt) = \mbox{Logit}^{-1}(prob\_adopt + params * data + exposure) #' } -#' Where exposure is the number of adopters in the agent's network. -#' +#' Where exposure is the number of adopters in the agent's network. +#' #' Another important difference is that the transmission network is not #' necesary useful since adoption in this model is not from a particular #' neighbor. -#' +#' #' @examples #' set.seed(2223) #' n <- 10000 -#' +#' #' # Generating synthetic data on a matrix with 2 columns. #' X <- cbind( -#' age = sample(1:100, n, replace = TRUE), -#' female = sample.int(2, n, replace = TRUE) - 1 +#' age = sample(1:100, n, replace = TRUE), +#' female = sample.int(2, n, replace = TRUE) - 1 #' ) -#' +#' #' adopt_chatgpt <- ModelDiffNet( #' "ChatGPT", #' prevalence = .01, @@ -45,29 +45,29 @@ #' data = X, #' params = c(1, 4) #' ) -#' +#' #' # Simulating a population from smallworld #' agents_smallworld(adopt_chatgpt, n, 8, FALSE, .01) -#' +#' #' # Running the model for 50 steps #' run(adopt_chatgpt, 50) -#' +#' #' # Plotting the model #' plot(adopt_chatgpt) #' @aliases epiworld_diffnet ModelDiffNet <- function( - name, - prevalence, - prob_adopt, - normalize_exposure = TRUE, - data = matrix(nrow = 0, ncol = 0), - data_cols = 1L:ncol(data), - params = vector("double") - ) { + name, + prevalence, + prob_adopt, + normalize_exposure = TRUE, + data = matrix(nrow = 0, ncol = 0), + data_cols = 1L:ncol(data), + params = vector("double") + ) { if (length(data) == 0L) data_cols <- vector("integer") - else + else data_cols <- as.integer(data_cols) - 1L structure( @@ -80,7 +80,7 @@ ModelDiffNet <- function( ncol(data), data_cols, params - ), + ), class = c("epiworld_diffnet", "epiworld_model") ) @@ -94,4 +94,3 @@ ModelDiffNet <- function( plot.epiworld_diffnet <- function(x, main = get_name(x), ...) { plot_epi(x, main = main, ...) } - diff --git a/R/ModelSEIR.R b/R/ModelSEIR.R index a081109e..505d20fa 100644 --- a/R/ModelSEIR.R +++ b/R/ModelSEIR.R @@ -2,27 +2,27 @@ #' #' @param name String. Name of the virus. #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param incubation_days Numeric scalar greater than 0. Average number of +#' @param incubation_days Numeric scalar greater than 0. Average number of #' incubation days. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. -#' @param x Object of class SEIR. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. +#' @param x Object of class SEIR. #' @param ... Currently ignore. #' @export #' @family Models #' @aliases epiworld_seir -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. The user must provide a vector of proportions indicating the following #' values: (1) Proportion of non-infected agents who are removed, and (2) #' Proportion of exposed agents to be set as infected. #' @returns #' - The `ModelSEIR`function returns a model of class [epiworld_model]. -#' @examples -#' model_seir <- ModelSEIR(name = "COVID-19", prevalence = 0.01, -#' transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4) -#' +#' @examples +#' model_seir <- ModelSEIR(name = "COVID-19", prevalence = 0.01, +#' transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4) +#' #' # Adding a small world population #' agents_smallworld( #' model_seir, @@ -30,12 +30,12 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_seir, ndays = 100, seed = 1912) #' model_seir -#' +#' #' plot(model_seir, main = "SEIR Model") #' @seealso epiworld-methods ModelSEIR <- function( @@ -44,20 +44,20 @@ ModelSEIR <- function( transmission_rate, incubation_days, recovery_rate -) { - + ) { + structure( ModelSEIR_cpp(name, prevalence, transmission_rate, incubation_days, recovery_rate), class = c("epiworld_seir", "epiworld_model") ) - + } #' @rdname ModelSEIR #' @param main Title of the plot -#' @returns The `plot` function returns a plot of the SEIR model of class +#' @returns The `plot` function returns a plot of the SEIR model of class #' [epiworld_model]. #' @export plot.epiworld_seir <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSEIRCONN.R b/R/ModelSEIRCONN.R index 3e01bff1..d55a6d8c 100644 --- a/R/ModelSEIRCONN.R +++ b/R/ModelSEIRCONN.R @@ -1,71 +1,71 @@ #' Susceptible Exposed Infected Removed model (SEIR connected) -#' +#' #' The SEIR connected model implements a model where all agents are connected. #' This is equivalent to a compartmental model ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SEIR_model)). -#' +#' #' @param name String. Name of the virus. #' @param n Integer greater than zero. Population size. #' @param prevalence Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. -#' @param incubation_days Numeric scalar greater than 0. Average number of +#' @param incubation_days Numeric scalar greater than 0. Average number of #' incubation days. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery_rate. -#' @param x Object of class SEIRCONN. -#' @param ... Currently ignore. +#' @param x Object of class SEIRCONN. +#' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export #' @family Models #' @aliases epiworld_seirconn #' @returns #' - The `ModelSEIRCONN`function returns a model of class [epiworld_model]. -#' @examples +#' @examples #' # An example with COVID-19 #' model_seirconn <- ModelSEIRCONN( #' name = "COVID-19", -#' prevalence = 0.01, +#' prevalence = 0.01, #' n = 10000, -#' contact_rate = 2, -#' incubation_days = 7, +#' contact_rate = 2, +#' incubation_days = 7, #' transmission_rate = 0.5, #' recovery_rate = 0.3 #' ) -#' +#' #' # Running and printing #' run(model_seirconn, ndays = 100, seed = 1912) #' model_seirconn -#' +#' #' plot(model_seirconn) -#' +#' #' # Adding the flu -#' flu <- virus("Flu", .9, 1/7, prevalence = 0.001, as_proportion = TRUE) +#' flu <- virus("Flu", .9, 1 / 7, prevalence = 0.001, as_proportion = TRUE) #' add_virus(model_seirconn, flu) -#' +#' #' #' # Running and printing #' run(model_seirconn, ndays = 100, seed = 1912) #' model_seirconn -#' +#' #' plot(model_seirconn) #' @seealso epiworld-methods ModelSEIRCONN <- function( - name, n, prevalence, contact_rate, transmission_rate, + name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate -) { - + ) { + structure( - ModelSEIRCONN_cpp(name, n, prevalence, contact_rate, - transmission_rate, incubation_days, recovery_rate), + ModelSEIRCONN_cpp(name, n, prevalence, contact_rate, + transmission_rate, incubation_days, recovery_rate), class = c("epiworld_seirconn", "epiworld_model") ) - + } #' @rdname ModelSEIRCONN #' @export -#' @returns The `plot` function returns a plot of the SEIRCONN model of class +#' @returns The `plot` function returns a plot of the SEIRCONN model of class #' [epiworld_model]. #' @param main Title of the plot. plot.epiworld_seirconn <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSEIRD.R b/R/ModelSEIRD.R index fc1394e4..25a45a1a 100644 --- a/R/ModelSEIRD.R +++ b/R/ModelSEIRD.R @@ -2,16 +2,16 @@ #' #' @param name String. Name of the virus. #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param incubation_days Numeric scalar greater than 0. Average number of +#' @param incubation_days Numeric scalar greater than 0. Average number of #' incubation days. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. #' @param death_rate Numeric scalar between 0 and 1. Rate of death from virus. -#' @param x Object of class SEIRD. +#' @param x Object of class SEIRD. #' @param ... Currently ignore. #' @export -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. The user must provide a vector of proportions indicating the following #' values: (1) Proportion of exposed agents who are infected, (2) @@ -21,11 +21,11 @@ #' @aliases epiworld_seird #' @returns #' - The `ModelSEIRD`function returns a model of class [epiworld_model]. -#' @examples -#' model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, -#' transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4, -#' death_rate = 0.01) -#' +#' @examples +#' model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, +#' transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4, +#' death_rate = 0.01) +#' #' # Adding a small world population #' agents_smallworld( #' model_seird, @@ -33,12 +33,12 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_seird, ndays = 100, seed = 1912) #' model_seird -#' +#' #' plot(model_seird, main = "SEIRD Model") #' @seealso epiworld-methods ModelSEIRD <- function( @@ -48,8 +48,8 @@ ModelSEIRD <- function( incubation_days, recovery_rate, death_rate -) { - + ) { + structure( ModelSEIRD_cpp( name = name, @@ -58,17 +58,17 @@ ModelSEIRD <- function( incubation_days = incubation_days, recovery_rate = recovery_rate, death_rate = death_rate - ), + ), class = c("epiworld_seird", "epiworld_model") ) - + } #' @rdname ModelSEIRD #' @param main Title of the plot -#' @returns The `plot` function returns a plot of the SEIRD model of class +#' @returns The `plot` function returns a plot of the SEIRD model of class #' [epiworld_model]. #' @export plot.epiworld_seird <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) -} + plot_epi(x, main = main, ...) +} diff --git a/R/ModelSEIRDCONN.R b/R/ModelSEIRDCONN.R index 42290d5a..f019e5a0 100644 --- a/R/ModelSEIRDCONN.R +++ b/R/ModelSEIRDCONN.R @@ -1,23 +1,23 @@ #' Susceptible Exposed Infected Removed Deceased model (SEIRD connected) -#' +#' #' The SEIRD connected model implements a model where all agents are connected. #' This is equivalent to a compartmental model ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SEIR_model)). -#' +#' #' @param name String. Name of the virus. #' @param n Integer greater than zero. Population size. #' @param prevalence Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. -#' @param incubation_days Numeric scalar greater than 0. Average number of +#' @param incubation_days Numeric scalar greater than 0. Average number of #' incubation days. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery_rate. #' @param death_rate Numeric scalar between 0 and 1. Probability of death. -#' @param x Object of class SEIRCONN. -#' @param ... Currently ignore. +#' @param x Object of class SEIRCONN. +#' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. The user must provide a vector of proportions indicating the following #' values: (1) Proportion of exposed agents who are infected, (2) @@ -27,58 +27,58 @@ #' @aliases epiworld_seirdconn #' @returns #' - The `ModelSEIRDCONN`function returns a model of class [epiworld_model]. -#' @examples +#' @examples #' # An example with COVID-19 #' model_seirdconn <- ModelSEIRDCONN( #' name = "COVID-19", -#' prevalence = 0.01, +#' prevalence = 0.01, #' n = 10000, -#' contact_rate = 2, -#' incubation_days = 7, +#' contact_rate = 2, +#' incubation_days = 7, #' transmission_rate = 0.5, #' recovery_rate = 0.3, #' death_rate = 0.01 #' ) -#' +#' #' # Running and printing #' run(model_seirdconn, ndays = 100, seed = 1912) #' model_seirdconn -#' +#' #' plot(model_seirdconn) -#' +#' #' # Adding the flu #' flu <- virus( -#' "Flu", prob_infecting = .3, recovery_rate = 1/7, +#' "Flu", prob_infecting = .3, recovery_rate = 1 / 7, #' prob_death = 0.001, #' prevalence = 0.001, as_proportion = TRUE #' ) #' add_virus(model = model_seirdconn, virus = flu) -#' +#' #' #' # Running and printing #' run(model_seirdconn, ndays = 100, seed = 1912) #' model_seirdconn -#' +#' #' plot(model_seirdconn) #' @seealso epiworld-methods ModelSEIRDCONN <- function( - name, n, prevalence, contact_rate, transmission_rate, + name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate, death_rate -) { - + ) { + structure( - ModelSEIRDCONN_cpp(name, n, prevalence, contact_rate, - transmission_rate, incubation_days, recovery_rate, - death_rate), + ModelSEIRDCONN_cpp(name, n, prevalence, contact_rate, + transmission_rate, incubation_days, recovery_rate, + death_rate), class = c("epiworld_seirdconn", "epiworld_model") ) - + } #' @rdname ModelSEIRDCONN #' @export -#' @returns The `plot` function returns a plot of the SEIRDCONN model of class +#' @returns The `plot` function returns a plot of the SEIRDCONN model of class #' [epiworld_model]. #' @param main Title of the plot. plot.epiworld_seirdconn <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSEIRMixing.R b/R/ModelSEIRMixing.R index 1dee4eac..161803e7 100644 --- a/R/ModelSEIRMixing.R +++ b/R/ModelSEIRMixing.R @@ -2,46 +2,46 @@ #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. #' @param incubation_days Numeric scalar. Average number of days in the #' incubation period. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery. -#' @param x Object of class SIRCONN. +#' @param x Object of class SIRCONN. #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @param contact_matrix Matrix of contact rates between individuals. #' @export #' @family Models -#' @details +#' @details #' The `contact_matrix` is a matrix of contact rates between entities. The -#' matrix should be of size `n x n`, where `n` is the number of entities. +#' matrix should be of size `n x n`, where `n` is the number of entities. #' This is a row-stochastic matrix, i.e., the sum of each row should be 1. -#' +#' #' The [initial_states] function allows the user to set the initial state of the #' model. In particular, the user can specify how many of the non-infected #' agents have been removed at the beginning of the simulation. #' @returns #' - The `ModelSEIRMixing`function returns a model of class [epiworld_model]. #' @aliases epiworld_seirmixing -#' -#' @examples -#' +#' +#' @examples +#' #' # Start off creating three entities. #' # Individuals will be distribured randomly between the three. #' e1 <- entity("Population 1", 3e3, as_proportion = FALSE) #' e2 <- entity("Population 2", 3e3, as_proportion = FALSE) #' e3 <- entity("Population 3", 3e3, as_proportion = FALSE) -#' +#' #' # Row-stochastic matrix (rowsums 1) #' cmatrix <- c( #' c(0.9, 0.05, 0.05), #' c(0.1, 0.8, 0.1), #' c(0.1, 0.2, 0.7) #' ) |> matrix(byrow = TRUE, nrow = 3) -#' +#' #' N <- 9e3 -#' +#' #' flu_model <- ModelSEIRMixing( #' name = "Flu", #' n = N, @@ -52,40 +52,40 @@ #' incubation_days = 7, #' contact_matrix = cmatrix #' ) -#' +#' #' # Adding the entities to the model #' flu_model |> #' add_entity(e1) |> #' add_entity(e2) |> #' add_entity(e3) -#' +#' #' set.seed(331) #' run(flu_model, ndays = 100) #' summary(flu_model) #' plot_incidence(flu_model) -#' +#' #' @seealso epiworld-methods ModelSEIRMixing <- function( - name, n, prevalence, contact_rate, transmission_rate, + name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate, contact_matrix -) { - + ) { + structure( ModelSEIRMixing_cpp( - name, n, prevalence, contact_rate, - transmission_rate, incubation_days, + name, n, prevalence, contact_rate, + transmission_rate, incubation_days, recovery_rate, as.vector(contact_matrix) - ), + ), class = c("epiworld_seirmixing", "epiworld_model") ) - + } #' @rdname ModelSEIRMixing #' @export -#' @returns The `plot` function returns a plot of the SEIRMixing model of class +#' @returns The `plot` function returns a plot of the SEIRMixing model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_seirmixing <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSIR.R b/R/ModelSIR.R index 76a38297..0fce74d5 100644 --- a/R/ModelSIR.R +++ b/R/ModelSIR.R @@ -1,29 +1,29 @@ #' SIR model #' @param name String. Name of the virus. -#' +#' #' Susceptible-Infected-Recovered model ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SIR_model)). -#' +#' #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. -#' @param x Object of class SIR. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. +#' @param x Object of class SIR. #' @param ... Additional arguments passed to [graphics::plot]. #' @export #' @family Models -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. In particular, the user can specify how many of the non-infected #' agents have been removed at the beginning of the simulation. #' @aliases epiworld_sir #' @returns #' - The `ModelSIR` function returns a model of class [epiworld_model]. -#' @examples -#' model_sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, -#' transmission_rate = 0.9, recovery_rate = 0.1) -#' +#' @examples +#' model_sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, +#' transmission_rate = 0.9, recovery_rate = 0.1) +#' #' # Adding a small world population #' agents_smallworld( #' model_sir, @@ -31,33 +31,32 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_sir, ndays = 100, seed = 1912) #' model_sir -#' +#' #' # Plotting #' plot(model_sir) #' @seealso epiworld-methods ModelSIR <- function( name, prevalence, transmission_rate, recovery_rate -) { - + ) { + structure( ModelSIR_cpp(name, prevalence, transmission_rate, recovery_rate), class = c("epiworld_sir", "epiworld_model") ) - + } #' @rdname ModelSIR #' @export -#' @returns -#' - The `plot` function returns a plot of the SIR model of class +#' @returns +#' - The `plot` function returns a plot of the SIR model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_sir <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } - diff --git a/R/ModelSIRCONN.R b/R/ModelSIRCONN.R index 6a92b8a2..cd6af26f 100644 --- a/R/ModelSIRCONN.R +++ b/R/ModelSIRCONN.R @@ -2,23 +2,23 @@ #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery. -#' @param x Object of class SIRCONN. +#' @param x Object of class SIRCONN. #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export #' @family Models -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. In particular, the user can specify how many of the non-infected #' agents have been removed at the beginning of the simulation. #' @returns #' - The `ModelSIRCONN`function returns a model of class [epiworld_model]. #' @aliases epiworld_sirconn -#' -#' @examples +#' +#' @examples #' model_sirconn <- ModelSIRCONN( #' name = "COVID-19", #' n = 10000, @@ -27,30 +27,30 @@ #' transmission_rate = 0.4, #' recovery_rate = 0.95 #' ) -#' +#' #' # Running and printing #' run(model_sirconn, ndays = 100, seed = 1912) #' model_sirconn -#' +#' #' plot(model_sirconn, main = "SIRCONN Model") #' @seealso epiworld-methods ModelSIRCONN <- function( name, n, prevalence, contact_rate, transmission_rate, recovery_rate -) { - + ) { + structure( - ModelSIRCONN_cpp(name, n, prevalence, contact_rate, - transmission_rate, recovery_rate), + ModelSIRCONN_cpp(name, n, prevalence, contact_rate, + transmission_rate, recovery_rate), class = c("epiworld_sirconn", "epiworld_model") ) - + } #' @rdname ModelSIRCONN #' @export -#' @returns The `plot` function returns a plot of the SIRCONN model of class +#' @returns The `plot` function returns a plot of the SIRCONN model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_sirconn <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSIRD.R b/R/ModelSIRD.R index 08272f17..20f0297f 100644 --- a/R/ModelSIRD.R +++ b/R/ModelSIRD.R @@ -1,16 +1,16 @@ #' SIRD model #' @param name String. Name of the virus. -#' +#' #' Susceptible-Infected-Recovered-Deceased model ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1166587407#The_SIRD_model)). -#' +#' #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery_rate from virus. #' @param death_rate Numeric scalar between 0 and 1. Rate of death from virus. -#' @param x Object of class SIR. +#' @param x Object of class SIR. #' @param ... Additional arguments passed to [graphics::plot]. #' @export #' @details @@ -22,15 +22,15 @@ #' @aliases epiworld_sird #' @returns #' - The `ModelSIRD` function returns a model of class [epiworld_model]. -#' @examples +#' @examples #' model_sird <- ModelSIRD( -#' name = "COVID-19", -#' prevalence = 0.01, -#' transmission_rate = 0.9, -#' recovery_rate = 0.1, -#' death_rate = 0.01 +#' name = "COVID-19", +#' prevalence = 0.01, +#' transmission_rate = 0.9, +#' recovery_rate = 0.1, +#' death_rate = 0.01 #' ) -#' +#' #' # Adding a small world population #' agents_smallworld( #' model_sird, @@ -38,33 +38,32 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_sird, ndays = 100, seed = 1912) #' model_sird -#' +#' #' # Plotting #' plot(model_sird) #' @seealso epiworld-methods ModelSIRD <- function( name, prevalence, transmission_rate, recovery_rate, death_rate -) { - + ) { + structure( ModelSIRD_cpp(name, prevalence, transmission_rate, recovery_rate, death_rate), class = c("epiworld_sird", "epiworld_model") ) - + } #' @rdname ModelSIRD #' @export -#' @returns -#' - The `plot` function returns a plot of the SIRD model of class +#' @returns +#' - The `plot` function returns a plot of the SIRD model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_sird <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } - diff --git a/R/ModelSIRDCONN.R b/R/ModelSIRDCONN.R index b4b3fb13..8a78314f 100644 --- a/R/ModelSIRDCONN.R +++ b/R/ModelSIRDCONN.R @@ -2,15 +2,15 @@ #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery. #' @param death_rate Numeric scalar between 0 and 1. Probability of death. -#' @param x Object of class SIRDCONN. +#' @param x Object of class SIRDCONN. #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export -#' @details +#' @details #' The [initial_states] function allows the user to set the initial state of the #' model. The user must provide a vector of proportions indicating the following #' values: (1) proportion of non-infected agents already removed, and (2) proportion of @@ -19,8 +19,8 @@ #' @returns #' - The `ModelSIRDCONN`function returns a model of class [epiworld_model]. #' @aliases epiworld_sirdconn -#' -#' @examples +#' +#' @examples #' model_sirdconn <- ModelSIRDCONN( #' name = "COVID-19", #' n = 100000, @@ -30,31 +30,31 @@ #' recovery_rate = 0.5, #' death_rate = 0.1 #' ) -#' +#' #' # Running and printing #' run(model_sirdconn, ndays = 100, seed = 1912) #' model_sirdconn -#' +#' #' plot(model_sirdconn, main = "SIRDCONN Model") #' @seealso epiworld-methods ModelSIRDCONN <- function( name, n, prevalence, contact_rate, transmission_rate, recovery_rate, death_rate -) { - + ) { + structure( - ModelSIRDCONN_cpp(name, n, prevalence, contact_rate, - transmission_rate, recovery_rate, death_rate), + ModelSIRDCONN_cpp(name, n, prevalence, contact_rate, + transmission_rate, recovery_rate, death_rate), class = c("epiworld_sirdconn", "epiworld_model") ) - + } #' @rdname ModelSIRDCONN #' @export -#' @returns The `plot` function returns a plot of the SIRDCONN model of class +#' @returns The `plot` function returns a plot of the SIRDCONN model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_sirdconn <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSIRLogit.R b/R/ModelSIRLogit.R index dea4168b..9d01c8f2 100644 --- a/R/ModelSIRLogit.R +++ b/R/ModelSIRLogit.R @@ -1,4 +1,4 @@ -#' SIR Logistic model +#' SIR Logistic model #' @param vname Name of the virus. #' @param data A numeric matrix with `n` rows. #' @param coefs_infect Numeric vector. Coefficients associated to infect. @@ -10,69 +10,69 @@ #' @param prevalence Numeric scalar. Prevalence (initial state) in proportion. #' #' @export -#' @returns +#' @returns #' - The `ModelSIRLogit` function returns a model of class [epiworld_model]. -#' @examples +#' @examples #' set.seed(2223) #' n <- 100000 -#' -#' # Creating the data to use for the "ModelSIRLogit" function. It contains -#' # information on the sex of each agent and will be used to determine -#' # differences in disease progression between males and females. Note that -#' # the number of rows in these data are identical to n (100000). +#' +#' # Creating the data to use for the "ModelSIRLogit" function. It contains +#' # information on the sex of each agent and will be used to determine +#' # differences in disease progression between males and females. Note that +#' # the number of rows in these data are identical to n (100000). #' X <- cbind( #' Intercept = 1, #' Female = sample.int(2, n, replace = TRUE) - 1 #' ) -#' +#' #' # Declare coefficients for each sex regarding transmission_rate and recovery. #' coef_infect <- c(.1, -2, 2) #' coef_recover <- rnorm(2) -#' -#' # Feed all above information into the "ModelSIRLogit" function. +#' +#' # Feed all above information into the "ModelSIRLogit" function. #' model_logit <- ModelSIRLogit( #' "covid2", #' data = X, #' coefs_infect = coef_infect, -#' coefs_recover = coef_recover, +#' coefs_recover = coef_recover, #' coef_infect_cols = 1L:ncol(X), #' coef_recover_cols = 1L:ncol(X), #' prob_infection = .8, #' recovery_rate = .3, #' prevalence = .01 #' ) -#' +#' #' agents_smallworld(model_logit, n, 8, FALSE, .01) -#' +#' #' run(model_logit, 50) -#' +#' #' plot(model_logit) -#' +#' #' # Females are supposed to be more likely to become infected. #' rn <- get_reproductive_number(model_logit) -#' -#' # Probability of infection for males and females. +#' +#' # Probability of infection for males and females. #' (table( #' X[, "Female"], #' (1:n %in% rn$source) -#' ) |> prop.table())[,2] -#' +#' ) |> prop.table())[, 2] +#' #' # Looking into the individual agents. #' get_agents(model_logit) #' @family Models ModelSIRLogit <- function( - vname, - data, - coefs_infect, - coefs_recover, - coef_infect_cols, - coef_recover_cols, - prob_infection, - recovery_rate, - prevalence -) { - + vname, + data, + coefs_infect, + coefs_recover, + coef_infect_cols, + coef_recover_cols, + prob_infection, + recovery_rate, + prevalence + ) { + structure( ModelSIRLogit_cpp( vname, @@ -88,5 +88,5 @@ ModelSIRLogit <- function( ), class = c("epiworld_sir", "epiworld_model") ) - + } diff --git a/R/ModelSIRMixing.R b/R/ModelSIRMixing.R index d10b5c29..5feb2562 100644 --- a/R/ModelSIRMixing.R +++ b/R/ModelSIRMixing.R @@ -2,45 +2,45 @@ #' @param name String. Name of the virus #' @param prevalence Double. Initial proportion of individuals with the virus. #' @param contact_rate Numeric scalar. Average number of contacts per step. -#' @param transmission_rate Numeric scalar between 0 and 1. Probability of +#' @param transmission_rate Numeric scalar between 0 and 1. Probability of #' transmission. #' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery. -#' @param x Object of class SIRCONN. +#' @param x Object of class SIRCONN. #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @param contact_matrix Matrix of contact rates between individuals. #' @export #' @family Models -#' @details +#' @details #' The `contact_matrix` is a matrix of contact rates between entities. The -#' matrix should be of size `n x n`, where `n` is the number of entities. +#' matrix should be of size `n x n`, where `n` is the number of entities. #' This is a row-stochastic matrix, i.e., the sum of each row should be 1. -#' +#' #' The [initial_states] function allows the user to set the initial state of the #' model. In particular, the user can specify how many of the non-infected #' agents have been removed at the beginning of the simulation. #' @returns #' - The `ModelSIRMixing`function returns a model of class [epiworld_model]. #' @aliases epiworld_sirmixing -#' -#' @examples +#' +#' @examples #' # From the vignette -#' +#' #' # Start off creating three entities. #' # Individuals will be distribured randomly between the three. #' e1 <- entity("Population 1", 3e3, as_proportion = FALSE) #' e2 <- entity("Population 2", 3e3, as_proportion = FALSE) #' e3 <- entity("Population 3", 3e3, as_proportion = FALSE) -#' +#' #' # Row-stochastic matrix (rowsums 1) #' cmatrix <- c( #' c(0.9, 0.05, 0.05), #' c(0.1, 0.8, 0.1), #' c(0.1, 0.2, 0.7) #' ) |> matrix(byrow = TRUE, nrow = 3) -#' +#' #' N <- 9e3 -#' +#' #' flu_model <- ModelSIRMixing( #' name = "Flu", #' n = N, @@ -50,39 +50,39 @@ #' recovery_rate = 1 / 7, #' contact_matrix = cmatrix #' ) -#' +#' #' # Adding the entities to the model #' flu_model |> #' add_entity(e1) |> #' add_entity(e2) |> #' add_entity(e3) -#' +#' #' set.seed(331) #' run(flu_model, ndays = 100) #' summary(flu_model) #' plot_incidence(flu_model) -#' +#' #' @seealso epiworld-methods ModelSIRMixing <- function( name, n, prevalence, contact_rate, transmission_rate, recovery_rate, contact_matrix -) { - + ) { + structure( ModelSIRMixing_cpp( - name, n, prevalence, contact_rate, + name, n, prevalence, contact_rate, transmission_rate, recovery_rate, as.vector(contact_matrix) - ), + ), class = c("epiworld_sirmixing", "epiworld_model") ) - + } #' @rdname ModelSIRMixing #' @export -#' @returns The `plot` function returns a plot of the SIRMixing model of class +#' @returns The `plot` function returns a plot of the SIRMixing model of class #' [epiworld_model]. #' @param main Title of the plot plot.epiworld_sirmixing <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/ModelSIS.R b/R/ModelSIS.R index 075035b7..a874a3ee 100644 --- a/R/ModelSIS.R +++ b/R/ModelSIS.R @@ -1,23 +1,23 @@ #' SIS model -#' +#' #' Susceptible-Infected-Susceptible model (SIS) ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SIS_model)) #' #' @param name String. Name of the virus. #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery from virus. -#' @param x Object of class SIS. -#' @param ... Currently ignore. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery from virus. +#' @param x Object of class SIS. +#' @param ... Currently ignore. #' @export #' @family Models -#' @returns +#' @returns #' - The `ModelSIS` function returns a model of class [epiworld_model]. #' @aliases epiworld_sis -#' @examples -#' model_sis <- ModelSIS(name = "COVID-19", prevalence = 0.01, -#' transmission_rate = 0.9, recovery_rate = 0.1) -#' +#' @examples +#' model_sis <- ModelSIS(name = "COVID-19", prevalence = 0.01, +#' transmission_rate = 0.9, recovery_rate = 0.1) +#' #' # Adding a small world population #' agents_smallworld( #' model_sis, @@ -25,33 +25,33 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_sis, ndays = 100, seed = 1912) #' model_sis -#' +#' #' # Plotting #' plot(model_sis, main = "SIS Model") -#' +#' #' @seealso epiworld-methods ModelSIS <- function( name, prevalence, transmission_rate, recovery_rate) { - + structure( ModelSIS_cpp(name, prevalence, transmission_rate, recovery_rate), class = c("epiworld_sis", "epiworld_model") ) - + } #' @rdname ModelSIS #' @export -#' @returns -#' - The `plot` function returns a plot of the SIS model of class +#' @returns +#' - The `plot` function returns a plot of the SIS model of class #' [epiworld_model]. #' @param main Title of the plot. -plot.epiworld_sis <- function(x, main = get_name(x),...) { # col = NULL - plot_epi(x, main = main, ...) +plot.epiworld_sis <- function(x, main = get_name(x), ...) { # col = NULL + plot_epi(x, main = main, ...) } diff --git a/R/ModelSISD.R b/R/ModelSISD.R index 30abb6b5..bea44a08 100644 --- a/R/ModelSISD.R +++ b/R/ModelSISD.R @@ -1,29 +1,29 @@ #' SISD model -#' +#' #' Susceptible-Infected-Susceptible-Deceased model (SISD) ([wiki](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SIS_model)) #' #' @param name String. Name of the virus. #' @param prevalence Double. Initial proportion of individuals with the virus. -#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of +#' @param transmission_rate Numeric scalar between 0 and 1. Virus's rate of #' infection. -#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery from virus. +#' @param recovery_rate Numeric scalar between 0 and 1. Rate of recovery from virus. #' @param death_rate Numeric scalar between 0 and 1. Rate of death from virus. -#' @param x Object of class SISD. -#' @param ... Currently ignore. +#' @param x Object of class SISD. +#' @param ... Currently ignore. #' @export #' @family Models -#' @returns +#' @returns #' - The `ModelSISD` function returns a model of class [epiworld_model]. #' @aliases epiworld_sisd -#' @examples +#' @examples #' model_sisd <- ModelSISD( -#' name = "COVID-19", -#' prevalence = 0.01, -#' transmission_rate = 0.9, -#' recovery_rate = 0.1, -#' death_rate = 0.01 -#' ) -#' +#' name = "COVID-19", +#' prevalence = 0.01, +#' transmission_rate = 0.9, +#' recovery_rate = 0.1, +#' death_rate = 0.01 +#' ) +#' #' # Adding a small world population #' agents_smallworld( #' model_sisd, @@ -31,33 +31,33 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_sisd, ndays = 100, seed = 1912) #' model_sisd -#' +#' #' # Plotting #' plot(model_sisd, main = "SISD Model") -#' +#' #' @seealso epiworld-methods ModelSISD <- function( name, prevalence, transmission_rate, recovery_rate, death_rate) { - + structure( ModelSISD_cpp(name, prevalence, transmission_rate, recovery_rate, death_rate), class = c("epiworld_sisd", "epiworld_model") ) - + } #' @rdname ModelSISD #' @export -#' @returns -#' - The `plot` function returns a plot of the SISD model of class +#' @returns +#' - The `plot` function returns a plot of the SISD model of class #' [epiworld_model]. #' @param main Title of the plot. -plot.epiworld_sisd <- function(x, main = get_name(x),...) { # col = NULL - plot_epi(x, main = main, ...) +plot.epiworld_sisd <- function(x, main = get_name(x), ...) { # col = NULL + plot_epi(x, main = main, ...) } diff --git a/R/ModelSURV.R b/R/ModelSURV.R index b9312b72..33899c4e 100644 --- a/R/ModelSURV.R +++ b/R/ModelSURV.R @@ -1,5 +1,5 @@ #' SURV model -#' +#' #' #' @param name String. Name of the virus. #' @param prevalence Initial number of individuals with the virus. @@ -18,11 +18,11 @@ #' the chances of becoming infected. #' @param surveillance_prob Double. Probability of testing an agent. #' @param transmission_rate Double. Raw transmission probability. -#' @param prob_death Double. Raw probability of death for symptomatic +#' @param prob_death Double. Raw probability of death for symptomatic #' individuals. #' @param prob_noreinfect Double. Probability of no re-infection. -#' @param x Object of class SURV. -#' @param ... Currently ignore. +#' @param x Object of class SURV. +#' @param ... Currently ignore. #' @export #' @family Models #' @aliases epiworld_surv @@ -52,37 +52,37 @@ #' k = 5, #' d = FALSE, #' p = .01 -#' ) -#' +#' ) +#' #' # Running and printing #' run(model_surv, ndays = 100, seed = 1912) -#' model_surv -#' +#' model_surv +#' #' # Plotting #' plot(model_surv, main = "SURV Model") -#' +#' #' @seealso epiworld-methods ModelSURV <- function( - name, prevalence, efficacy_vax, latent_period, infect_period, prob_symptoms, - prop_vaccinated, prop_vax_redux_transm, prop_vax_redux_infect, - surveillance_prob, transmission_rate, prob_death, prob_noreinfect -) { - + name, prevalence, efficacy_vax, latent_period, infect_period, prob_symptoms, + prop_vaccinated, prop_vax_redux_transm, prop_vax_redux_infect, + surveillance_prob, transmission_rate, prob_death, prob_noreinfect + ) { + structure( - ModelSURV_cpp(name, prevalence, efficacy_vax, latent_period, infect_period, - prob_symptoms, prop_vaccinated, prop_vax_redux_transm, - prop_vax_redux_infect, surveillance_prob, transmission_rate, - prob_death, prob_noreinfect), + ModelSURV_cpp(name, prevalence, efficacy_vax, latent_period, infect_period, + prob_symptoms, prop_vaccinated, prop_vax_redux_transm, + prop_vax_redux_infect, surveillance_prob, transmission_rate, + prob_death, prob_noreinfect), class = c("epiworld_surv", "epiworld_model") ) - + } #' @rdname ModelSURV #' @export -#' @returns The `plot` function returns a plot of the SURV model of class +#' @returns The `plot` function returns a plot of the SURV model of class #' [epiworld_model]. #' @param main Title of the plot. plot.epiworld_surv <- function(x, main = get_name(x), ...) { # col = NULL - plot_epi(x, main = main, ...) + plot_epi(x, main = main, ...) } diff --git a/R/agents-methods.R b/R/agents-methods.R index 2c1bd5fb..81968896 100644 --- a/R/agents-methods.R +++ b/R/agents-methods.R @@ -1,57 +1,57 @@ #' Agents in epiworldR -#' +#' #' These functions provide read-access to the agents of the model. The #' `get_agents` function returns an object of class [epiworld_agents] which #' contains all the information about the agents in the model. The #' `get_agent` function returns the information of a single agent. #' And the `get_state` function returns the state of a single agent. -#' +#' #' @param model An object of class [epiworld_model]. #' @param x An object of class [epiworld_agents]. #' @seealso agents #' @export #' @aliases epiworld_agents -#' @return +#' @return #' - The `get_agents` function returns an object of class [epiworld_agents]. #' @examples -#' +#' #' model_sirconn <- ModelSIRCONN( -#' name = "COVID-19", -#' n = 10000, -#' prevalence = 0.01, -#' contact_rate = 5, -#' transmission_rate = 0.4, -#' recovery_rate = 0.95 +#' name = "COVID-19", +#' n = 10000, +#' prevalence = 0.01, +#' contact_rate = 5, +#' transmission_rate = 0.4, +#' recovery_rate = 0.95 #' ) -#' +#' #' run(model_sirconn, ndays = 100, seed = 1912) -#' -#' x <- get_agents(model_sirconn) # Storing all agent information into object of -#' # class epiworld_agents -#' -#' print(x, compressed = FALSE, max_print = 5) # Displaying detailed information of -#' # the first 5 agents using -#' # compressed=F. Using compressed=T -#' # results in less-detailed -#' # information about each agent. -#' -#' x[0] # Print information about the first agent. Substitute the agent of -#' # interest's position where '0' is. +#' +#' x <- get_agents(model_sirconn) # Storing all agent information into object of +#' # class epiworld_agents +#' +#' print(x, compressed = FALSE, max_print = 5) # Displaying detailed information of +#' # the first 5 agents using +#' # compressed=F. Using compressed=T +#' # results in less-detailed +#' # information about each agent. +#' +#' x[0] # Print information about the first agent. Substitute the agent of +#' # interest's position where '0' is. #' @name agents get_agents <- function(model, ...) UseMethod("get_agents") -#' @export +#' @export #' @rdname agents get_agents.epiworld_model <- function(model, ...) { - + res <- get_agents_cpp(model) - + structure( res, class = "epiworld_agents", model = model ) - + } #' @param x An object of class [epiworld_agents] @@ -62,61 +62,61 @@ get_agents.epiworld_model <- function(model, ...) { #' - The `[` method returns an object of class [epiworld_agent]. #' @aliases epiworld_agent `[.epiworld_agents` <- function(x, i) { - + structure( get_agent_cpp(x, i), class = "epiworld_agent", model = attr(x, "model") ) - + } `[.epiworld_agents` <- function(x, i) { - + structure( get_agent_cpp(x, i), class = "epiworld_agent", model = attr(x, "model") ) - + } #' @export #' @param compressed Logical scalar. When FALSE, it prints detailed information #' about the agent. #' @param ... Ignored -#' @returns -#' - The `print` function returns information about each individual agent of +#' @returns +#' - The `print` function returns information about each individual agent of #' class [epiworld_agent]. #' @rdname agents print.epiworld_agent <- function(x, compressed = FALSE, ...) { - + invisible(print_agent_cpp(x, attr(x, "model"), compressed)) - + } #' @export #' @param max_print Integer scalar. Maximum number of agents to print. #' @rdname agents print.epiworld_agents <- function(x, compressed = TRUE, max_print = 10, ...) { - + model <- attr(x, "model") cat(sprintf("Agents from the model \"%s\":\n", get_name(model))) n <- size(model) for (i in 1L:min(max_print, n)) { - + print(x[i - 1L], compressed) - + } - + if (n > max_print) cat(sprintf("... %i more agents ...\n", n - max_print)) - + invisible(x) - + } #' @export -#' @returns +#' @returns #' - The `get_state` function returns the state of the [epiworld_agents] object. #' @rdname agents get_state <- function(x) { diff --git a/R/agents.R b/R/agents.R index f4b4f8a1..34e775ae 100644 --- a/R/agents.R +++ b/R/agents.R @@ -4,13 +4,13 @@ stopifnot_agent <- function(x) { } #' Load agents to a model -#' +#' #' These functions provide access to the network of the model. The network is #' represented by an edgelist. The `agents_smallworld` function generates a -#' small world network with the Watts-Strogatz algorithm. The +#' small world network with the Watts-Strogatz algorithm. The #' `agents_from_edgelist` function loads a network from an edgelist. #' The `get_network` function returns the edgelist of the network. -#' +#' #' @param model Model object of class [epiworld_model]. #' @param source,target Integer vectors describing the source and target of #' in the edgelist. @@ -20,48 +20,48 @@ stopifnot_agent <- function(x) { #' @param p Probability of rewiring. #' @export #' @return -#' - The 'agents_smallworld' function returns a model with the agents +#' - The 'agents_smallworld' function returns a model with the agents #' loaded. #' @examples -#' +#' #' # Initializing SIR model with agents_smallworld -#' sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, -#' recovery_rate = 0.1) +#' sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, +#' recovery_rate = 0.1) #' agents_smallworld( -#' sir, -#' n = 1000, -#' k = 5, -#' d = FALSE, -#' p = .01 -#' ) +#' sir, +#' n = 1000, +#' k = 5, +#' d = FALSE, +#' p = .01 +#' ) #' run(sir, ndays = 100, seed = 1912) #' sir -#' +#' #' # We can also retrieve the network #' net <- get_network(sir) #' head(net) -#' +#' #' # Simulating a bernoulli graph #' set.seed(333) #' n <- 1000 -#' g <- matrix(runif(n ^ 2) < .01, nrow = n) +#' g <- matrix(runif(n^2) < .01, nrow = n) #' diag(g) <- FALSE #' el <- which(g, arr.ind = TRUE) - 1L -#' -#' +#' +#' #' # Generating an empty model #' sir <- ModelSIR("COVID-19", .01, .8, .3) #' agents_from_edgelist( #' sir, -#' source = el[,1], -#' target = el[,2], +#' source = el[, 1], +#' target = el[, 2], #' size = n, #' directed = TRUE #' ) -#' +#' #' # Running the simulation #' run(sir, 50) -#' +#' #' plot(sir) agents_smallworld <- function(model, n, k, d, p) UseMethod("agents_smallworld") @@ -75,17 +75,17 @@ agents_smallworld.epiworld_model <- function(model, n, k, d, p) { #' @export #' @return #' - The `agents_from_edgelist` function returns an empty model of class -#' `epiworld_model`. +#' `epiworld_model`. #' @rdname agents_smallworld agents_from_edgelist <- function( - model, source, target, size, directed - ) UseMethod("agents_from_edgelist") + model, source, target, size, directed + ) UseMethod("agents_from_edgelist") #' @export agents_from_edgelist.epiworld_model <- function( - model, source, target, size, directed - ) { - + model, source, target, size, directed + ) { + agents_from_edgelist_cpp( model, source, @@ -93,12 +93,12 @@ agents_from_edgelist.epiworld_model <- function( size, directed ) - + invisible(model) - + } -#' @export +#' @export #' @rdname agents_smallworld #' @aliases network #' @return @@ -109,8 +109,8 @@ get_network <- function(model) { get_network_cpp(model) } -#' @export -#' @return +#' @export +#' @return #' - `get_agents_states` returns an character vector with the states of the #' agents by the end of the simulation. #' @rdname agents_smallworld @@ -119,7 +119,7 @@ get_agents_states <- function(model) { get_agents_states_cpp(model) } -#' @export +#' @export #' @param agent Agent object of class `epiworld_agent`. #' @param virus Virus object of class `epiworld_virus`. #' @param state_new Integer scalar. New state of the agent after the action is executed. @@ -132,8 +132,8 @@ get_agents_states <- function(model) { #' - The function `add_virus_agent` adds a virus to an agent and #' returns the agent invisibly. add_virus_agent <- function( - agent, model, virus, state_new = -99, queue = -99 -) { + agent, model, virus, state_new = -99, queue = -99 + ) { stopifnot_model(model) stopifnot_virus(virus) @@ -152,8 +152,8 @@ add_virus_agent <- function( #' - The function `add_tool_agent` adds a tool to an agent and #' returns the agent invisibly. add_tool_agent <- function( - agent, model, tool, state_new = -99, queue = -99 -) { + agent, model, tool, state_new = -99, queue = -99 + ) { stopifnot_model(model) stopifnot_tool(tool) diff --git a/R/data.R b/R/data.R index 0f0f1508..20983a89 100644 --- a/R/data.R +++ b/R/data.R @@ -1,11 +1,11 @@ #' Accessing the database of epiworld -#' +#' #' Models in `epiworld` are stored in a database. This database can be accessed #' using the functions described in this manual page. Some elements of the -#' database are: the transition matrix, the incidence matrix, the reproductive +#' database are: the transition matrix, the incidence matrix, the reproductive #' number, the generation time, and daily incidence at the virus and tool level. -#' -#' +#' +#' #' @param x An object of class [`epiworld_sir`], [`epiworld_seir`], etc. #' any model. #' @param skip_zeros Logical scalar. When `FALSE` it will return all the @@ -20,75 +20,75 @@ #' name = "Disease", #' n = 10000, #' prevalence = 0.1, -#' contact_rate = 2.0, +#' contact_rate = 2.0, #' transmission_rate = 0.8, #' incubation_days = 7.0, #' recovery_rate = 0.3 #' ) -#' +#' #' # Running the simulation for 50 steps (days) #' set.seed(937) #' run(seirconn, 50) #' #' # Retrieving the transition probability -#' get_transition_probability(seirconn) -#' -#' # Retrieving date, state, and counts dataframe including any added tools +#' get_transition_probability(seirconn) +#' +#' # Retrieving date, state, and counts dataframe including any added tools #' get_hist_tool(seirconn) -#' -#' # Retrieving overall date, state, and counts dataframe +#' +#' # Retrieving overall date, state, and counts dataframe #' head(get_hist_total(seirconn)) -#' +#' #' # Retrieving date, state, and counts dataframe by variant #' head(get_hist_virus(seirconn)) -#' +#' #' # Retrieving (and plotting) the reproductive number #' rp <- get_reproductive_number(seirconn) #' plot(rp) # Also equivalent to plot_reproductive_number(seirconn) -#' +#' #' # We can go further and get all the history #' t_hist <- get_hist_transition_matrix(seirconn) -#' +#' #' head(t_hist) -#' +#' #' # And turn it into an array -#' as.array(t_hist)[,,1:3] -#' +#' as.array(t_hist)[, , 1:3] +#' #' # We cam also get (and plot) the incidence, as well as #' # the generation time #' inci <- plot_incidence(seirconn) #' gent <- plot_generation_time(seirconn) -#' +#' NULL #' @export -#' @returns -#' - The `get_hist_total` function returns an object of class +#' @returns +#' - The `get_hist_total` function returns an object of class #' [epiworld_hist_total]. #' @rdname epiworld-data #' @aliases epiworld_hist_total get_hist_total <- function(x) UseMethod("get_hist_total") #' @export -get_hist_total.epiworld_model <- function(x) { - +get_hist_total.epiworld_model <- function(x) { + res <- get_hist_total_cpp(x) structure( res, class = c("epiworld_hist_total", "epiworld_hist", "data.frame"), states = sort(unique(res$state)) ) - + } -#' @export +#' @export #' @rdname epiworld-data #' @returns #' - The `get_today_total` function returns a named vector with the #' total number of individuals in each state at the end of the simulation. get_today_total <- function(x) UseMethod("get_today_total") -#' @export +#' @export get_today_total.epiworld_model <- function(x) { get_today_total_cpp(x) } @@ -99,17 +99,17 @@ plot.epiworld_hist <- function(x, y, ...) { } #' @export -#' @returns -#' - The `get_hist_virus` function returns an object of class +#' @returns +#' - The `get_hist_virus` function returns an object of class #' [epiworld_hist_virus]. #' @rdname epiworld-data #' @aliases epiworld_hist_variant epiworld_hist_virus get_hist_virus <- function(x) UseMethod("get_hist_virus") #' @export -get_hist_virus.epiworld_model <- function(x) { +get_hist_virus.epiworld_model <- function(x) { res <- get_hist_virus_cpp(x) - + structure( res, class = c("epiworld_hist_virus", "epiworld_hist", "data.frame"), @@ -118,14 +118,14 @@ get_hist_virus.epiworld_model <- function(x) { } #' @export -#' @returns +#' @returns #' - The `get_hist_tool` function returns an object of [epiworld_hist_virus]. #' @rdname epiworld-data #' @aliases epiworld_hist_tool get_hist_tool <- function(x) UseMethod("get_hist_tool") #' @export -get_hist_tool.epiworld_model <- function(x) { +get_hist_tool.epiworld_model <- function(x) { res <- get_hist_tool_cpp(x) structure( res, @@ -135,26 +135,26 @@ get_hist_tool.epiworld_model <- function(x) { } #' @export -#' @returns -#' - The `get_transition_probability` function returns an object of class -#' `matrix`. +#' @returns +#' - The `get_transition_probability` function returns an object of class +#' `matrix`. #' @rdname epiworld-data get_transition_probability <- function(x) { UseMethod("get_transition_probability") } #' @export -get_transition_probability.epiworld_model <- function(x) { +get_transition_probability.epiworld_model <- function(x) { res <- get_transition_probability_cpp(x) s <- get_states(x) - + ns <- length(s) - + matrix(res, nrow = ns, ncol = ns, dimnames = list(s, s)) } #' @export -#' @returns +#' @returns #' - The `get_reproductive_number` function returns an object of class #' [epiworld_repnum]. #' @rdname epiworld-data @@ -173,8 +173,8 @@ get_reproductive_number.epiworld_model <- function(x) { #' @param plot Logical scalar. If `TRUE` (default), the function will the #' desired statistic. #' @param ylab,xlab,main,type Further parameters passed to [graphics::plot()] -#' @returns -#' - The `plot` function returns a plot of the reproductive number over time. +#' @returns +#' - The `plot` function returns a plot of the reproductive number over time. #' @export #' @importFrom stats sd quantile aggregate plot.epiworld_repnum <- function( @@ -186,7 +186,7 @@ plot.epiworld_repnum <- function( type = "b", plot = TRUE, ...) { - + if (nrow(x) == 0) { repnum <- data.frame( @@ -200,7 +200,6 @@ plot.epiworld_repnum <- function( ub = numeric() ) } else { - # Computing stats # Compute the mean and 95% CI of rt by virus and source_exposure_date using the repnum data.frame with the tapply function @@ -212,7 +211,7 @@ plot.epiworld_repnum <- function( by = list( virus_comb = x[["virus_comb"]], date = x[["source_exposure_date"]] - ), + ), FUN = function(x) { ci <- stats::quantile(x, c(0.025, 0.975), na.rm = TRUE) @@ -222,10 +221,10 @@ plot.epiworld_repnum <- function( sd = stats::sd(x, na.rm = TRUE), lb = ci[1], ub = ci[2] - ) + ) }, simplify = FALSE - ) + ) repnum <- cbind(repnum[, -3, drop = FALSE], do.call(rbind, repnum[, 3])) repnum <- repnum[order(repnum[["virus_comb"]], repnum[["date"]]), , drop = FALSE] @@ -237,15 +236,15 @@ plot.epiworld_repnum <- function( by = "virus_comb", all.x = TRUE, all.y = FALSE - ) - + ) + rownames(repnum) <- NULL # Reordering columns repnum <- repnum[, c( "virus_id", "virus", "date", "avg", "n", "sd", "lb", "ub", "virus_comb" - )] + )] } @@ -253,19 +252,19 @@ plot.epiworld_repnum <- function( # Nviruses vlabs <- sort(unique(x[, "virus_comb"])) nviruses <- length(vlabs) - + # # Figuring out the range yran <- range(repnum[["avg"]], na.rm = TRUE) xran <- range(repnum[["date"]], na.rm = TRUE) - + # Plotting ------------------------------------------------------------------- if (plot) { for (i in seq_along(vlabs)) { tmp <- repnum[repnum[["virus_comb"]] == vlabs[i], ] - + if (i == 1L) { - + graphics::plot( x = tmp[["date"]], y = tmp[["avg"]], @@ -283,7 +282,7 @@ plot.epiworld_repnum <- function( ) next } - + graphics::lines( x = tmp[["date"]], y = tmp[["avg"]], @@ -294,11 +293,11 @@ plot.epiworld_repnum <- function( type = type, ... ) - + } if (nviruses > 1L) { - + graphics::legend( "topright", legend = vlabs, @@ -308,7 +307,7 @@ plot.epiworld_repnum <- function( title = "Virus", bty = "n" ) - + } } @@ -317,7 +316,7 @@ plot.epiworld_repnum <- function( repnum[["virus_comb"]] <- NULL invisible(repnum) - + } #' @export @@ -332,7 +331,7 @@ plot_reproductive_number <- function(x, ...) { #' @export #' @rdname epiworld-data -#' @returns +#' @returns #' - `get_hist_transition_matrix` returns a [data.frame] with four columns: #' "state_from", "state_to", "date", and "counts." get_hist_transition_matrix <- function(x, skip_zeros = FALSE) @@ -340,47 +339,47 @@ get_hist_transition_matrix <- function(x, skip_zeros = FALSE) #' @export get_hist_transition_matrix.epiworld_model <- function(x, skip_zeros = FALSE) { - + res <- get_hist_transition_matrix_cpp(x, skip_zeros) class(res) <- c("epiworld_hist_transition", class(res)) - + attr(res, "states") <- get_states(x) attr(res, "nsteps") <- get_ndays(x) - + res - + } #' @export -#' @returns -#' - The `as.array` method for `epiworld_hist_transition` objects turns the +#' @returns +#' - The `as.array` method for `epiworld_hist_transition` objects turns the #' `data.frame` returned by `get_hist_transition_matrix` into an array of #' `nstates x nstates x (ndays + 1)` #' entries, where the first entry is the initial state. #' @rdname epiworld-data as.array.epiworld_hist_transition <- function(x, ...) { - + states <- attr(x, "states") n_states <- length(states) n_steps <- attr(x, "nsteps") - + res <- array( 0L, dim = c(n_states, n_states, n_steps + 1), # Includes the baseline dimnames = list(states, states, 0:n_steps) - ) - - res[cbind(x[,1], x[,2], x[,3])] <- x[,4] - + ) + + res[cbind(x[, 1], x[, 2], x[, 3])] <- x[, 4] + res } #' @export #' @rdname epiworld-data -#' @returns -#' - The `plot_incidence` function returns a plot originating from the object +#' @returns +#' - The `plot_incidence` function returns a plot originating from the object #' `get_hist_transition_matrix`. #' @details The `plot_incidence` function is a wrapper between #' [get_hist_transition_matrix] and it's plot method. @@ -388,28 +387,28 @@ plot_incidence <- function(x, ...) { plot(get_hist_transition_matrix(x), ...) } -#' @export +#' @export #' @returns -#' - The `plot` function returns a plot which originates from the +#' - The `plot` function returns a plot which originates from the #' `epiworld_hist_transition` object. #' @rdname epiworld-data #' @details The plot method for the `epiworld_hist_transition` class plots the #' daily incidence of each state. The function returns the data frame used for #' plotting. plot.epiworld_hist_transition <- function( - x, - type = "b", - xlab = "Day (step)", - ylab = "Counts", - main = "Daily incidence", - plot = TRUE, - ... - ) { - + x, + type = "b", + xlab = "Day (step)", + ylab = "Counts", + main = "Daily incidence", + plot = TRUE, + ... + ) { + if (!inherits(x, "epiworld_hist_transition")) { stop("The object must be of class 'epiworld_hist_transition'") } - + states <- attr(x, "states") n_states <- length(states) n_steps <- attr(x, "nsteps") @@ -431,27 +430,27 @@ plot.epiworld_hist_transition <- function( yran <- range(res) xran <- range(0:n_steps) for (i in is_not_zero) { - + col <- states[i] if (i == 1L) { - - graphics::plot( - x = as.integer(rownames(res)), - y = res[[col]], - col = i, - lwd = 2, - lty = i, - type = type, - xlab = xlab, - ylab = ylab, - main = main, - ylim = yran, - xlim = xran, - ... - ) - - next + + graphics::plot( + x = as.integer(rownames(res)), + y = res[[col]], + col = i, + lwd = 2, + lty = i, + type = type, + xlab = xlab, + ylab = ylab, + main = main, + ylim = yran, + xlim = xran, + ... + ) + + next } @@ -469,7 +468,7 @@ plot.epiworld_hist_transition <- function( # Creating a legend if (n_states > 1L) { - + graphics::legend( "topright", legend = states, @@ -479,10 +478,10 @@ plot.epiworld_hist_transition <- function( title = "States", bty = "n" ) - + } } - + invisible(res) } @@ -502,13 +501,13 @@ get_transmissions.epiworld_diffnet <- function(x) { #' @export get_transmissions.epiworld_model <- function(x) { - + res <- get_transmissions_cpp(x) structure( res, class = c("epiworld_transmissions", class(res)) ) - + } #' @export @@ -517,33 +516,33 @@ get_transmissions.epiworld_model <- function(x) { #' - The function `get_generation_time` returns a `data.frame` with #' the following columns: "agent", "virus_id", "virus", "date", and "gentime". get_generation_time <- function(x) { - - stopifnot_model(x) - res <- get_generation_time_cpp(x) - # Replacing -1 with NAs - res[["gentime"]][res[["gentime"]] == -1] <- NA_integer_ + stopifnot_model(x) + res <- get_generation_time_cpp(x) + + # Replacing -1 with NAs + res[["gentime"]][res[["gentime"]] == -1] <- NA_integer_ + + structure( + res, + class = c("epiworld_generation_time", class(res)), + n_steps = get_ndays(x) + ) - structure( - res, - class = c("epiworld_generation_time", class(res)), - n_steps = get_ndays(x) - ) - } #' @export #' @rdname epiworld-data plot.epiworld_generation_time <- function( - x, - type = "b", - xlab = "Day (step)", - ylab = "Avg. Generation Time", - main = "Generation Time", - plot = TRUE, - ... - ) { - + x, + type = "b", + xlab = "Day (step)", + ylab = "Avg. Generation Time", + main = "Generation Time", + plot = TRUE, + ... + ) { + if (!inherits(x, "epiworld_generation_time")) { stop("The object must be of class 'epiworld_generation_time'") } @@ -559,12 +558,12 @@ plot.epiworld_generation_time <- function( x[["gentime"]], by = list( date = x[["date"]], virus_comb = x[["virus_comb"]] - ), + ), FUN = function(x) { ci <- stats::quantile( x, probs = c(0.025, 0.975), na.rm = TRUE - ) - + ) + data.frame( avg = mean(x, na.rm = TRUE), n = sum(!is.na(x)), @@ -574,11 +573,11 @@ plot.epiworld_generation_time <- function( ) }, simplify = FALSE - ) + ) gt <- cbind(gt[, -3, drop = FALSE], do.call(rbind, gt[, 3])) gt <- gt[order(gt[["virus_comb"]], gt[["date"]]), , drop = FALSE] - + # Merging the virus and virus_id column of x to repnum gt <- merge( gt, @@ -586,15 +585,15 @@ plot.epiworld_generation_time <- function( by = "virus_comb", all.x = TRUE, all.y = FALSE - ) - + ) + rownames(gt) <- NULL # Replacing NaNs with NAs gt <- as.data.frame(lapply(gt, function(x) { x[is.nan(x)] <- NA x - })) + })) if (plot) { # Number of viruses @@ -604,9 +603,9 @@ plot.epiworld_generation_time <- function( for (i in 1L:n_viruses) { gt_i <- gt[gt[["virus_comb"]] == viruses[i], , drop = FALSE] - + if (i == 1L) { - + graphics::plot( x = gt_i[["date"]], y = gt_i[["avg"]], @@ -621,11 +620,11 @@ plot.epiworld_generation_time <- function( xlim = range(gt[["date"]], na.rm = TRUE), ... ) - + next - + } - + graphics::points( x = gt_i[["date"]], y = gt_i[["avg"]], @@ -639,7 +638,7 @@ plot.epiworld_generation_time <- function( # Creating a legend if (n_viruses > 1L) { - + graphics::legend( "topright", legend = viruses, @@ -649,7 +648,7 @@ plot.epiworld_generation_time <- function( title = "Virus", bty = "n" ) - + } } @@ -658,16 +657,15 @@ plot.epiworld_generation_time <- function( gt[["virus_comb"]] <- NULL invisible(gt) - - + + } -#' @export +#' @export #' @rdname epiworld-data -#' @return +#' @return #' - The function `plot_generation_time` is a wrapper for [plot] and #' [get_generation_time]. plot_generation_time <- function(x, ...) { plot(get_generation_time(x), ...) } - diff --git a/R/entity.R b/R/entity.R index e10c4227..52403368 100644 --- a/R/entity.R +++ b/R/entity.R @@ -12,37 +12,37 @@ stopifnot_entity_distfun <- function(distfun) { } #' Get entities -#' +#' #' Entities in `epiworld` are objects that can contain agents. #' @param model Model object of class `epiworld_model`. -#' -#' @details +#' +#' @details #' Epiworld entities are especially useful for mixing models, particularly #' [ModelSIRMixing] and [ModelSEIRMixing]. -#' +#' #' @name entities -#' @export -#' @examples +#' @export +#' @examples #' # Creating a mixing model #' mymodel <- ModelSIRMixing( -#' name = "My model", -#' n = 10000, -#' prevalence = .001, -#' contact_rate = 10, -#' transmission_rate = .1, -#' recovery_rate = 1/7, -#' contact_matrix = matrix(c(.9, .1, .1, .9), 2, 2) +#' name = "My model", +#' n = 10000, +#' prevalence = .001, +#' contact_rate = 10, +#' transmission_rate = .1, +#' recovery_rate = 1 / 7, +#' contact_matrix = matrix(c(.9, .1, .1, .9), 2, 2) #' ) -#' +#' #' ent1 <- entity("First", 5000, FALSE) #' ent2 <- entity("Second", 5000, FALSE) -#' +#' #' mymodel |> #' add_entity(ent1) |> #' add_entity(ent2) -#' +#' #' run(mymodel, ndays = 100, seed = 1912) -#' +#' #' summary(mymodel) get_entities <- function(model) { @@ -67,14 +67,14 @@ print.epiworld_entities <- function(x, ...) { invisible(x) } -#' @export +#' @export #' @rdname entities #' @param x Object of class `epiworld_entities`. #' @param i Integer index. `[.epiworld_entities` <- function(x, i) { stopifnot_entity(x) - + if (i > get_entity_size(x)) { stop("Index out of bounds.") } @@ -87,14 +87,14 @@ print.epiworld_entities <- function(x, ...) { } -#' @export +#' @export #' @param name Character scalar. Name of the entity. #' @param prevalence Numeric scalar. Prevalence of the entity. #' @param as_proportion Logical scalar. If `TRUE`, `prevalence` is interpreted #' as a proportion. #' @param to_unassigned Logical scalar. If `TRUE`, the entity is added to the #' unassigned pool. -#' @return +#' @return #' - The function `entity` creates an entity object. #' @rdname entities entity <- function(name, prevalence, as_proportion, to_unassigned = TRUE) { @@ -105,13 +105,13 @@ entity <- function(name, prevalence, as_proportion, to_unassigned = TRUE) { as.double(prevalence), as.logical(as_proportion), as.logical(to_unassigned) - ), + ), class = "epiworld_entity" ) } -#' @export +#' @export #' @rdname entities #' @param entity Entity object of class `epiworld_entity`. #' @return @@ -123,7 +123,7 @@ get_entity_size <- function(entity) { #' @export #' @rdname entities -#' @return +#' @return #' - The function `get_entity_name` returns the name of the entity. get_entity_name <- function(entity) { stopifnot_entity(entity) @@ -133,13 +133,13 @@ get_entity_name <- function(entity) { #' @export #' @rdname entities #' @param agent Agent object of class `epiworld_agent`. -#' @return +#' @return #' - The function `entity_add_agent` adds an agent to the entity. entity_add_agent <- function( - entity, - agent, - model = attr(entity, "model") - ) { + entity, + agent, + model = attr(entity, "model") + ) { stopifnot_entity(entity) stopifnot_agent(agent) @@ -152,10 +152,10 @@ entity_add_agent <- function( #' @export #' @rdname entities #' @param id Integer scalar. Entity id to remove (starting from zero). -#' @return +#' @return #' - The function `rm_entity` removes an entity from the model. rm_entity <- function(model, id) { - + stopifnot_model(model) rm_entity_cpp(model, entity) @@ -165,32 +165,32 @@ rm_entity <- function(model, id) { #' @export #' @rdname entities add_entity <- function( - model, - entity -) { + model, + entity + ) { stopifnot_model(model) stopifnot_entity(entity) add_entity_cpp( model, entity - ) + ) invisible(model) } -#' @export +#' @export #' @rdname entities -#' @param agents_id Integer vector. -#' @param entities_id Integer vector. -#' @return +#' @param agents_id Integer vector. +#' @param entities_id Integer vector. +#' @return #' - The function `load_agents_entities_ties` loads agents into entities. load_agents_entities_ties <- function( - model, - agents_id, - entities_id -) { + model, + agents_id, + entities_id + ) { stopifnot_model(model) if (!inherits(agents_id, "integer")) { @@ -207,9 +207,9 @@ load_agents_entities_ties <- function( } -#' @export +#' @export #' @rdname entities -#' @return +#' @return #' - The function `entity_get_agents` returns an integer vector with the agents #' in the entity (ids). entity_get_agents <- function(entity) { @@ -219,7 +219,7 @@ entity_get_agents <- function(entity) { } -#' @export +#' @export print.epiworld_entity <- function(x, ...) { print_entity_cpp(x) invisible(x) @@ -231,13 +231,13 @@ print.epiworld_entity <- function(x, ...) { #' as a proportion. #' @rdname entities distribute_entity_randomly <- function( - prevalence, - as_proportion, - to_unassigned = TRUE -) { + prevalence, + as_proportion, + to_unassigned = TRUE + ) { structure( - distribute_entity_randomly_cpp( + distribute_entity_randomly_cpp( as.double(prevalence), as.logical(as_proportion), as.logical(to_unassigned) @@ -251,9 +251,9 @@ distribute_entity_randomly <- function( #' @param agents_ids Integer vector. Ids of the agents to distribute. #' @rdname entities distribute_entity_to_set <- function( - agents_ids -) { - + agents_ids + ) { + structure( distribute_entity_to_set_cpp( as.integer(agents_ids) @@ -263,13 +263,13 @@ distribute_entity_to_set <- function( } -#' @export +#' @export #' @rdname entities #' @param distfun Distribution function object of class `epiworld_distribution_entity`. set_distribution_entity <- function( - entity, - distfun -) { + entity, + distfun + ) { stopifnot_entity(entity) stopifnot_entity_distfun(distfun) @@ -277,4 +277,4 @@ set_distribution_entity <- function( invisible(entity) -} \ No newline at end of file +} diff --git a/R/epiworldR-package.R b/R/epiworldR-package.R index d1a5744d..cfa3ff5c 100644 --- a/R/epiworldR-package.R +++ b/R/epiworldR-package.R @@ -2,4 +2,3 @@ #' @useDynLib epiworldR, .registration = TRUE #' @importFrom graphics boxplot plot "_PACKAGE" - diff --git a/R/epiworldR-package.R.in b/R/epiworldR-package.R.in index 21204cbb..415ae58f 100644 --- a/R/epiworldR-package.R.in +++ b/R/epiworldR-package.R.in @@ -2,4 +2,3 @@ #' @useDynLib @EPIWORLD_NAME@, .registration = TRUE #' @importFrom graphics boxplot plot "_PACKAGE" - diff --git a/R/functions-renamed.R b/R/functions-renamed.R index 61556bed..d6094f5a 100644 --- a/R/functions-renamed.R +++ b/R/functions-renamed.R @@ -2,7 +2,7 @@ #' @description #' Starting version 0.0-4, epiworld changed how it refered to "actions." #' Following more traditional ABMs, actions are now called "events." -#' +#' #' @param ... Arguments to be passed to the new function. #' @param model Model object of class `epiworld_model`. #' @param tool Tool object of class `epiworld_tool`. @@ -10,8 +10,8 @@ #' @name epiworldR-deprecated NULL -#' @param n Deprecated. -#' @export +#' @param n Deprecated. +#' @export #' @rdname epiworldR-deprecated add_tool_n <- function(model, tool, n) { @@ -24,12 +24,12 @@ add_tool_n <- function(model, tool, n) { as_proportion = TRUE ) ) - + add_tool(model, tool) } -#' @export +#' @export #' @rdname epiworldR-deprecated add_virus_n <- function(model, virus, n) { @@ -40,9 +40,9 @@ add_virus_n <- function(model, virus, n) { distfun = distribute_virus_randomly( prevalence = n, as_proportion = TRUE - ) + ) ) add_virus(model, virus) -} \ No newline at end of file +} diff --git a/R/global-actions.R b/R/global-actions.R index 2504e204..e305a298 100644 --- a/R/global-actions.R +++ b/R/global-actions.R @@ -1,15 +1,15 @@ #' Global Actions -#' +#' #' Global actions are functions that are executed at each time step of the #' simulation. They are useful for implementing interventions, such as #' vaccination, isolation, and social distancing by means of tools. -#' +#' #' @export #' @param prob Numeric scalar. A probability between 0 and 1. #' @param tool An object of class [tool]. #' @name global-actions -#' @examples +#' @examples #' # Simple model #' model_sirconn <- ModelSIRCONN( #' name = "COVID-19", @@ -19,7 +19,7 @@ #' transmission_rate = 0.4, #' recovery_rate = 0.95 #' ) -#' +#' #' # Creating a tool #' epitool <- tool( #' name = "Vaccine", @@ -27,20 +27,20 @@ #' as_proportion = FALSE, #' susceptibility_reduction = .9, #' transmission_reduction = .5, -#' recovery_enhancer = .5, +#' recovery_enhancer = .5, #' death_reduction = .9 #' ) -#' -#' +#' +#' #' # Adding a global action #' vaccine_day_20 <- globalevent_tool(epitool, .2, day = 20) #' add_globalevent(model_sirconn, vaccine_day_20) -#' +#' #' # Running and printing #' run(model_sirconn, ndays = 40, seed = 1912) #' model_sirconn #' plot_incidence(model_sirconn) -#' +#' #' # Example 2: Changing the contact rate ------------------------------------- #' model_sirconn2 <- ModelSIRCONN( #' name = "COVID-19", @@ -50,10 +50,10 @@ #' transmission_rate = 0.4, #' recovery_rate = 0.95 #' ) -#' +#' #' closure_day_10 <- globalevent_set_params("Contact rate", 0, day = 10) #' add_globalevent(model_sirconn2, closure_day_10) -#' +#' #' # Running and printing #' run(model_sirconn2, ndays = 40, seed = 1912) #' model_sirconn2 @@ -61,23 +61,23 @@ #' @returns #' - The `globalevent_set_params` function returns an object of class #' [epiworld_globalevent_set_param] and [epiworld_globalevent]. -#' -#' - `globalevent_tool` returns an object of class +#' +#' - `globalevent_tool` returns an object of class #' [epiworld_globalevent_tool] and [epiworld_globalevent]. -#' +#' #' - `globalevent_tool_logit` returns an object of class #' [epiworld_globalevent_tool_logit] and [epiworld_globalevent]. -#' @aliases +#' @aliases #' epiworld_globalevent_set_param #' epiworld_globalevent_tool #' epiworld_globalevent_tool_logit #' epiworld_globalevent #' actions -#' +#' globalevent_tool <- function( - tool, prob, - name = get_name_tool(tool), day = -99 - ) { + tool, prob, + name = get_name_tool(tool), day = -99 + ) { structure( globalevent_tool_cpp(tool, prob, name, day), @@ -93,7 +93,7 @@ globalaction_tool <- function(...) { .Defunct( new = "globalevent_tool" - ) + ) } @@ -107,9 +107,9 @@ globalaction_tool <- function(...) { #' `vars` is an integer vector indicating the position of the variables in the #' model. globalevent_tool_logit <- function( - tool, vars, coefs, - name = get_name_tool(tool), day = -99 - ) { + tool, vars, coefs, + name = get_name_tool(tool), day = -99 + ) { stopifnot_tool(tool) @@ -120,7 +120,7 @@ globalevent_tool_logit <- function( as.double(coefs), name, as.integer(day) - ), + ), class = c("epiworld_globalevent_tool_logit", "epiworld_globalevent"), tool = tool, call = match.call() @@ -134,13 +134,13 @@ globalaction_tool_logit <- function(...) { .Defunct( new = "globalevent_tool_logit" - ) + ) globalevent_tool_logit(...) - + } -#' @export +#' @export #' @param param Character scalar. The name of the parameter to be set. #' @param value Numeric scalar. The value of the parameter. #' @rdname global-actions @@ -148,9 +148,9 @@ globalaction_tool_logit <- function(...) { #' the model. The parameter is specified by its name `param` and the value by #' `value`. globalevent_set_params <- function( - param, value, - name = paste0("Set ", param, " to ", value), day = -99 - ) { + param, value, + name = paste0("Set ", param, " to ", value), day = -99 + ) { structure( globalevent_set_param_cpp( @@ -158,7 +158,7 @@ globalevent_set_params <- function( as.double(value), name, as.integer(day) - ), + ), class = c("epiworld_globalevent_set_param", "epiworld_globalevent"), param = param, value = as.double(value), @@ -173,22 +173,22 @@ globalaction_set_params <- function(...) { .Defunct( new = "globalevent_set_params" - ) + ) globalevent_set_params(...) - + } -#' @export +#' @export #' @rdname global-actions #' @param fun Function. The function to be executed. #' @details The function `globalevent_fun` allows to specify a function to be #' executed at a given day. The function object must receive an object of class #' [epiworld_model] as only argument. -#' @examples +#' @examples #' # Example using `globalevent_fun` to record the state of the #' # agents at each time step. -#' +#' #' # We start by creating an SIR connected model #' model <- ModelSIRCONN( #' name = "SIR with Global Saver", @@ -197,28 +197,28 @@ globalaction_set_params <- function(...) { #' contact_rate = 5, #' transmission_rate = 0.4, #' recovery_rate = 0.3 -#' ) -#' +#' ) +#' #' # We create the object where the history of the agents will be stored #' agents_history <- NULL -#' +#' #' # This function prints the total number of agents in each state #' # and stores the history of the agents in the object `agents_history` #' hist_saver <- function(m) { -#' +#' #' message("Today's totals are: ", paste(get_today_total(m), collapse = ", ")) -#' +#' #' # We use the `<<-` operator to assign the value to the global variable #' # `agents_history` (see ?"<<-") #' agents_history <<- cbind( #' agents_history, #' get_agents_states(m) -#' ) -#' +#' ) +#' #' } globalevent_fun <- function( - fun, name = deparse(substitute(fun)), day = -99 - ) { + fun, name = deparse(substitute(fun)), day = -99 + ) { structure( globalevent_fun_cpp(fun, name, as.integer(day)), @@ -235,10 +235,10 @@ globalaction_fun <- function(...) { .Defunct( new = "globalevent_fun" - ) + ) globalevent_fun(...) - + } #' @export @@ -273,11 +273,10 @@ print.epiworld_globalevent <- function(x, ...) { #' - The function `add_globalevent` returns the model with the added #' action. add_globalevent <- function(model, action) { - + if (length(attr(action, "tool"))) add_tool(model, attr(action, "tool")) invisible(add_globalevent_cpp(model, action)) } - diff --git a/R/make_saver.R b/R/make_saver.R index 74ecf634..72064ae7 100644 --- a/R/make_saver.R +++ b/R/make_saver.R @@ -1,9 +1,9 @@ #' Run multiple simulations at once -#' +#' #' The `run_multiple` function allows running multiple simulations at once. #' When available, users can take advantage of parallel computing to speed up #' the process. -#' +#' #' @param m,ndays,seed See [run]. #' @param saver An object of class [epiworld_saver]. #' @param nsims Integer. Number of replicats @@ -12,10 +12,10 @@ #' @param nthreads Integer. Number of threads (parallel computing.) #' @param reset When `TRUE` (default,) resets the simulation. #' @param verbose When `TRUE` (default,) prints a progress bar. -#' +#' #' @details #' Currently, the following elements can be saved: -#' +#' #' - `total_hist` History of the model (total numbers per time). #' - `virus_info` Information about `viruses`. #' - `virus_hist` Changes in `viruses`. @@ -25,7 +25,7 @@ #' - `transition` Transition matrices. #' - `reproductive` Reproductive number. #' - `generation` Estimation of generation time. -#' +#' #' @returns #' - In the case of `make_saver`, an list of class `epiworld_saver`. #' @examples @@ -35,28 +35,28 @@ #' n = 1000, #' contact_rate = 2, #' transmission_rate = 0.9, recovery_rate = 0.1 -#' ) -#' +#' ) +#' #' # Generating a saver #' saver <- make_saver("total_hist", "reproductive") -#' +#' #' # Running and printing #' run_multiple(model_sir, ndays = 100, nsims = 50, saver = saver, nthreads = 2) -#' +#' #' # Retrieving the results #' ans <- run_multiple_get_results(model_sir) -#' +#' #' head(ans$total_hist) #' head(ans$reproductive) -#' +#' #' # Plotting #' multi_sir <- run_multiple_get_results(model_sir)$total_hist -#' multi_sir <- multi_sir[multi_sir$date <= 20,] +#' multi_sir <- multi_sir[multi_sir$date <= 20, ] #' plot(multi_sir) -#' +#' #' @export -#' @returns -#' - The `run_multiple` function runs a specified number of simulations and +#' @returns +#' - The `run_multiple` function runs a specified number of simulations and #' returns a model object of class [epiworld_model]. run_multiple <- function( m, ndays, nsims, @@ -65,7 +65,7 @@ run_multiple <- function( reset = TRUE, verbose = TRUE, nthreads = 1L -) UseMethod("run_multiple") + ) UseMethod("run_multiple") #' @export run_multiple.epiworld_model <- function( @@ -75,8 +75,8 @@ run_multiple.epiworld_model <- function( reset = TRUE, verbose = TRUE, nthreads = 1L -) { - + ) { + if (!inherits(saver, "epiworld_saver")) stop("-saver- should be of class \"epiworld_saver\"") @@ -85,12 +85,12 @@ run_multiple.epiworld_model <- function( fnames <- list.files( path = dirname(saver$fn), full.names = TRUE - ) + ) if (length(fnames)) { unlink(fnames, expand = FALSE) } - + run_multiple_cpp( m, ndays, @@ -101,45 +101,44 @@ run_multiple.epiworld_model <- function( verbose, nthreads ) - + attr(m, "saver") <- saver - + invisible(m) - + } #' @export #' @rdname run_multiple -#' @returns +#' @returns #' - The `run_multiple_get_results` function returns a named list with the #' data specified by `make_saver`. #' @importFrom utils read.table run_multiple_get_results <- function(m) { - + if (!inherits(m, "epiworld_model")) stop("-m- must be of class `epiworld_model`.") - + # Get the filepath saver <- attr(m, "saver") - - if (!length(saver)) + + if (!length(saver)) stop("No -saver- found. -run_multiple_get_results- can only be used after using -run_multiple-.") - + output <- vector("list", length(saver$what)) names(output) <- saver$what - + for (i in saver$what) { - # Listing the files fnames <- list.files( path = dirname(saver$fn), pattern = sprintf("%s\\.csv", i), full.names = TRUE ) - + # Reading the files output[[i]] <- lapply(fnames, utils::read.table, sep = " ", header = TRUE) - + # Getting number of simulation output[[i]] <- lapply(seq_along(fnames), function(j) { if (nrow(output[[i]][[j]]) > 0) @@ -147,14 +146,16 @@ run_multiple_get_results <- function(m) { else NULL }) - + # Putting all together output[[i]] <- do.call(rbind, output[[i]]) - + # If there are no observations, then - err_msg <- tryCatch({ - class(output[[i]]) <- c("epiworld_multiple_save_i", class(output[[i]])) - }, error = function(e) e + err_msg <- tryCatch( + { + class(output[[i]]) <- c("epiworld_multiple_save_i", class(output[[i]])) + }, + error = function(e) e ) if (inherits(err_msg, "error")) { @@ -162,26 +163,25 @@ run_multiple_get_results <- function(m) { warning( "When retrieving the saved results, for the case of ", i, ", there were no observations." - ) + ) class(output[[i]]) <- c( "epiworld_multiple_save_i", class(output[[i]]) - ) + ) } attr(output[[i]], "what") <- i - + } - + structure(output, class = c("epiworld_multiple_save", class(output))) - + } #' @export plot.epiworld_multiple_save <- function(x, y = NULL, ...) { - # what <- attr(x, "what") lapply(x, plot) @@ -196,61 +196,59 @@ plot.epiworld_multiple_save_i <- function(x, y = NULL, ...) { warning( "When plotting the saved results, for the case of ", what, ", there were no observations." - ) + ) return(NULL) } - + # If it is not reproductive number, then... if (what != "reproductive") { - + oldpar <- graphics::par( - mfrow = c(2, floor(length(unique(x$state))/2)) - ) + mfrow = c(2, floor(length(unique(x$state)) / 2)) + ) on.exit(graphics::par(oldpar)) - + for (what in unique(x$state)) { graphics::boxplot( counts ~ date, - data = x[x$state == what,,drop=FALSE], + data = x[x$state == what, , drop = FALSE], main = what, xlab = "Date", ylab = "Counts", border = "black", las = 2 - ) - + ) + } - + } else { - + plot.epiworld_multiple_save_reproductive_number(x, ...) } - - + + } #' @export plot.epiworld_multiple_save_reproductive_number <- function(x, y = NULL, ...) { - # Identifying sims sims <- sort(unique(x[["sim_num"]])) - + totals <- NULL for (s in sims) { - # Subsetting the data - x_tmp <- x[x[["sim_num"]] == s,, drop = FALSE] - + x_tmp <- x[x[["sim_num"]] == s, , drop = FALSE] + # Computing daily values totals <- rbind( totals, plot.epiworld_repnum(x_tmp, plot = FALSE) ) - + } - + graphics::boxplot( avg ~ date, data = totals, @@ -259,10 +257,10 @@ plot.epiworld_multiple_save_reproductive_number <- function(x, y = NULL, ...) { ylab = "rt", border = "black", las = 2 - ) - + ) + invisible(totals) - + } #' @export @@ -272,9 +270,9 @@ make_saver <- function( ..., fn = "" ) { - + what <- list(...) - + # Any missmatch? available <- c( "total_hist", @@ -287,22 +285,22 @@ make_saver <- function( "reproductive", "generation" ) - + not_in_available <- which(!(what %in% available)) if (length(not_in_available)) { stop( "The following elements in -what- are not supported: \"", paste(what[not_in_available], collapse = "\" , \""), "\"" - ) + ) } - + what_bool <- as.list(available %in% what) names(what_bool) <- available - + # Checking the filename file_output <- TRUE - + # Using tempfile to generate directories id <- basename(tempfile("epiworldR-")) @@ -316,9 +314,9 @@ make_saver <- function( } else if (!dir.exists(dirname(fn))) { stop("The directory \"", dirname(fn), "\" does not exists.") } - + what_bool$fn <- fn - + # Generating the saver structure( list( @@ -327,22 +325,21 @@ make_saver <- function( file_output = file_output, what = available[which(available %in% what)], id = id - ), + ), class = "epiworld_saver" ) - + } #' @export print.epiworld_saver <- function(x, ...) { - + cat("A saver for -run_multiple-\n") cat("Saves the following:", paste(x$what, sep = ", "), "\n") cat("To file :", ifelse(x$file_output, "yes", "no"), "\n") if (x$file_output) cat("Saver pattern :", x$fn) - + invisible(x) } - diff --git a/R/model-methods.R b/R/model-methods.R index 61210d0d..464927b2 100644 --- a/R/model-methods.R +++ b/R/model-methods.R @@ -2,19 +2,19 @@ stopifnot_model <- function(model) { if (!inherits(model, "epiworld_model")) { stop( "The -model- object must be of class \"epiworld_model\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(model), collapse = ", ") ) } } #' Methods for epiworldR objects -#' +#' #' The functions described in this section are methods for objects of class #' `epiworld_model`. Besides of printing and plotting, other methods provide #' access to manipulate model parameters, getting information about the model #' and running the simulation. -#' +#' #' @param x An object of class `epiworld_model`. #' @param ndays Number of days (steps) of the simulation. #' @param seed Seed to set for initializing random number generator. @@ -23,83 +23,83 @@ stopifnot_model <- function(model) { #' @name epiworld-methods #' @aliases epiworld_model #' @examples -#' +#' #' model_sirconn <- ModelSIRCONN( -#' name = "COVID-19", -#' n = 10000, -#' prevalence = 0.01, -#' contact_rate = 5, -#' transmission_rate = 0.4, -#' recovery_rate = 0.95 +#' name = "COVID-19", +#' n = 10000, +#' prevalence = 0.01, +#' contact_rate = 5, +#' transmission_rate = 0.4, +#' recovery_rate = 0.95 #' ) -#' -#' # Queuing - If you wish to implement the queuing function, declare whether -#' # you would like it "on" or "off", if any. +#' +#' # Queuing - If you wish to implement the queuing function, declare whether +#' # you would like it "on" or "off", if any. #' queuing_on(model_sirconn) #' queuing_off(model_sirconn) #' run(model_sirconn, ndays = 100, seed = 1912) -#' -#' # Verbose - "on" prints the progress bar on the screen while "off" -#' # deactivates the progress bar. Declare which function you want to implement, -#' # if any. +#' +#' # Verbose - "on" prints the progress bar on the screen while "off" +#' # deactivates the progress bar. Declare which function you want to implement, +#' # if any. #' verbose_on(model_sirconn) #' verbose_off(model_sirconn) #' run(model_sirconn, ndays = 100, seed = 1912) -#' +#' #' get_states(model_sirconn) # Returns all unique states found within the model. -#' -#' get_param(model_sirconn, 'Contact rate') # Returns the value of the selected -#' # parameter within the model object. -#' # In order to view the parameters, -#' # run the model object and find the -#' # "Model parameters" section. -#' -#' set_param(model_sirconn, 'Contact rate', 2) # Allows for adjustment of model -#' # parameters within the model -#' # object. In this example, the -#' # Contact rate parameter is -#' # changed to 2. You can now rerun -#' # the model to observe any -#' # differences. -#' -#' set_name(model_sirconn, 'My Epi-Model') # This function allows for setting -#' # a name for the model. Running the -#' # model object, the name of the model -#' # is now reflected next to "Name of -#' # the model". -#' -#' get_name(model_sirconn) # Returns the set name of the model. -#' -#' get_n_viruses(model_sirconn) # Returns the number of viruses in the model. -#' # In this case, there is only one virus: -#' # "COVID-19". -#' -#' get_n_tools(model_sirconn) # Returns the number of tools in the model. In -#' # this case, there are zero tools. -#' +#' +#' get_param(model_sirconn, "Contact rate") # Returns the value of the selected +#' # parameter within the model object. +#' # In order to view the parameters, +#' # run the model object and find the +#' # "Model parameters" section. +#' +#' set_param(model_sirconn, "Contact rate", 2) # Allows for adjustment of model +#' # parameters within the model +#' # object. In this example, the +#' # Contact rate parameter is +#' # changed to 2. You can now rerun +#' # the model to observe any +#' # differences. +#' +#' set_name(model_sirconn, "My Epi-Model") # This function allows for setting +#' # a name for the model. Running the +#' # model object, the name of the model +#' # is now reflected next to "Name of +#' # the model". +#' +#' get_name(model_sirconn) # Returns the set name of the model. +#' +#' get_n_viruses(model_sirconn) # Returns the number of viruses in the model. +#' # In this case, there is only one virus: +#' # "COVID-19". +#' +#' get_n_tools(model_sirconn) # Returns the number of tools in the model. In +#' # this case, there are zero tools. +#' #' get_ndays(model_sirconn) # Returns the length of the simulation in days. This -#' # will match "ndays" within the "run" function. -#' -#' get_n_replicates(model_sirconn) # Returns the number of replicates of the -#' # model. -#' -#' size(model_sirconn) # Returns the population size in the model. In this case, -#' # there are 10,000 agents in the model. +#' # will match "ndays" within the "run" function. +#' +#' get_n_replicates(model_sirconn) # Returns the number of replicates of the +#' # model. +#' +#' size(model_sirconn) # Returns the population size in the model. In this case, +#' # there are 10,000 agents in the model. #' # Set Agents Data -#' # First, your data matrix must have the same number of rows as agents in the -#' # model. Below is a generated matrix which will be passed into the -#' # "set_agents_data" function. -#' data <- matrix(data=runif(20000, min=0, max=100), nrow=10000, ncol=2) +#' # First, your data matrix must have the same number of rows as agents in the +#' # model. Below is a generated matrix which will be passed into the +#' # "set_agents_data" function. +#' data <- matrix(data = runif(20000, min = 0, max = 100), nrow = 10000, ncol = 2) #' set_agents_data(model_sirconn, data) -#' get_agents_data_ncols(model_sirconn) # Returns number of columns -#' -#' get_virus(model_sirconn, 0) # Returns information about the first virus in -#' # the model (index begins at 0). -#' -#' add_tool(model_sirconn, tool("Vaccine", .9, .9, .5, 1, prevalence = 0.5, as_prop = TRUE)) -#' get_tool(model_sirconn, 0) # Returns information about the first tool in the -#' # model. In this case, there are no tools so an -#' # error message will occur. +#' get_agents_data_ncols(model_sirconn) # Returns number of columns +#' +#' get_virus(model_sirconn, 0) # Returns information about the first virus in +#' # the model (index begins at 0). +#' +#' add_tool(model_sirconn, tool("Vaccine", .9, .9, .5, 1, prevalence = 0.5, as_prop = TRUE)) +#' get_tool(model_sirconn, 0) # Returns information about the first tool in the +#' # model. In this case, there are no tools so an +#' # error message will occur. queuing_on <- function(x) UseMethod("queuing_on") #' @export @@ -130,9 +130,9 @@ queuing_off.epiworld_model <- function(x) { #' @name epiworld-methods #' @export -#' @returns -#' - The `verbose_on` and `verbose_off` functions return the same model, however -#' `verbose_off` returns the model with no progress bar. +#' @returns +#' - The `verbose_on` and `verbose_off` functions return the same model, however +#' `verbose_off` returns the model with no progress bar. #' @details #' The `verbose_on` and `verbose_off` functions activate and deactivate printing #' progress on screen, respectively. Both functions return the model (`x`) invisibly. @@ -154,8 +154,8 @@ verbose_on.epiworld_model <- function(x) { } #' @export -#' @returns -#' - The `run` function returns the simulated model of class `epiworld_model`. +#' @returns +#' - The `run` function returns the simulated model of class `epiworld_model`. #' @rdname epiworld-methods run <- function(model, ndays, seed = sample.int(1e4, 1)) UseMethod("run") @@ -171,12 +171,12 @@ print.epiworld_model <- function(x, ...) { invisible(x) } -#' @export +#' @export #' @returns #' - The `summary` function prints a more detailed view of the model, and returns the same model invisibly. #' @rdname epiworld-methods #' @param object Object of class `epiworld_model`. -#' @param ... Additional arguments. +#' @param ... Additional arguments. summary.epiworld_model <- function(object, ...) { print_cpp(object, lite = FALSE) invisible(object) @@ -184,7 +184,7 @@ summary.epiworld_model <- function(object, ...) { #' @export #' @returns -#' - The `get_states` function returns the unique states found in a model. +#' - The `get_states` function returns the unique states found in a model. #' @rdname epiworld-methods get_states <- function(x) UseMethod("get_states") @@ -193,8 +193,8 @@ get_states.epiworld_model <- function(x) get_states_cpp(x) #' @export #' @param pname String. Name of the parameter. -#' @returns -#' - The `get_param` function returns a selected parameter from the model object +#' @returns +#' - The `get_param` function returns a selected parameter from the model object #' of class `epiworld_model`. #' @rdname epiworld-methods get_param <- function(x, pname) UseMethod("get_param") @@ -207,9 +207,9 @@ get_param.epiworld_model <- function(x, pname) { #' @export #' @param pval Numeric. Value of the parameter. -#' @returns +#' @returns #' - The `set_param` function does not return a value but instead alters a -#' parameter value. +#' parameter value. #' @rdname epiworld-methods set_param <- function(x, pname, pval) UseMethod("set_param") @@ -221,8 +221,8 @@ set_param.epiworld_model <- function(x, pname, pval) { #' @export #' @param mname String. Name of the model. -#' @returns -#' - The `set_name` function does not return a value but instead alters an object +#' @returns +#' - The `set_name` function does not return a value but instead alters an object #' of `epiworld_model`. #' @rdname epiworld-methods set_name <- function(x, mname) UseMethod("set_name") @@ -244,7 +244,7 @@ get_name.epiworld_model <- function(x) { get_name_cpp(x) } -#' @export +#' @export #' @rdname epiworld-methods #' @returns #' - `get_n_viruses` returns the number of viruses of the model. @@ -254,7 +254,7 @@ get_n_viruses <- function(x) UseMethod("get_n_viruses") get_n_viruses.epiworld_model <- function(x) get_n_viruses_cpp(x) -#' @export +#' @export #' @rdname epiworld-methods #' @returns #' - `get_n_tools` returns the number of tools of the model. @@ -264,7 +264,7 @@ get_n_tools <- function(x) UseMethod("get_n_tools") get_n_tools.epiworld_model <- function(x) get_n_tools_cpp(x) -#' @export +#' @export #' @rdname epiworld-methods #' @returns #' - `get_ndays` returns the number of days of the model. @@ -274,9 +274,9 @@ get_ndays <- function(x) UseMethod("get_ndays") get_ndays.epiworld_model <- function(x) get_ndays_cpp(x) -#' @export +#' @export #' @rdname epiworld-methods -#' @returns +#' @returns #' - `get_n_replicates` returns the number of replicates of the model. get_n_replicates <- function(x) UseMethod("get_n_replicates") @@ -286,9 +286,9 @@ get_n_replicates.epiworld_model <- function(x) get_n_replicates_cpp(x) #' @export #' @rdname epiworld-methods -#' @returns +#' @returns #' - `size.epiworld_model` returns the number of agents in the model. -#' +#' size <- function(x) UseMethod("size") #' @export @@ -298,38 +298,38 @@ size.epiworld_model <- function(x) size_cpp(x) #' @export #' @param data A numeric matrix. #' @returns -#' - The 'set_agents_data' function returns an object of class DataFrame. +#' - The 'set_agents_data' function returns an object of class DataFrame. #' @rdname epiworld-methods set_agents_data <- function(model, data) { - + if (!inherits(data, "matrix") | mode(data) != "numeric") stop("-data- must be a numeric (mode) matrix (class).") - + if (size(model) != nrow(data)) stop( "The number of rows in -data- (", nrow(data), ") doesn't match the number of agents in the model (", size(model), ")." - ) - + ) + invisible(set_agents_data_cpp(model = model, data = data, ncols = ncol(data))) - + } #' @export -#' @returns -#' - 'get_agents_data_ncols' returns the number of columns in the model dataframe. +#' @returns +#' - 'get_agents_data_ncols' returns the number of columns in the model dataframe. #' @rdname epiworld-methods get_agents_data_ncols <- function(model) { - + get_agents_data_ncols_cpp(model) - + } #' @export #' @param virus_pos Integer. Relative location (starting from 0) of the virus #' in the model -#' @returns +#' @returns #' - 'get_virus' returns a [virus]. #' @rdname epiworld-methods get_virus <- function(model, virus_pos) { @@ -342,7 +342,7 @@ get_virus <- function(model, virus_pos) { #' @export #' @param tool_pos Integer. Relative location (starting from 0) of the tool #' in the model -#' @returns +#' @returns #' - `get_tool` returns a [tool]. #' @rdname epiworld-methods get_tool <- function(model, tool_pos) { @@ -352,10 +352,10 @@ get_tool <- function(model, tool_pos) { ) } -#' @export +#' @export #' @param proportions Numeric vector. Proportions in which agents will be #' distributed (see details). -#' @return +#' @return #' - `inital_states` returns the model with an updated initial state. #' @rdname epiworld-methods initial_states <- function(model, proportions) { @@ -379,4 +379,3 @@ clone_model <- function(model) { class = class(model) ) } - diff --git a/R/plot_epi.R b/R/plot_epi.R index 36234322..b263a931 100644 --- a/R/plot_epi.R +++ b/R/plot_epi.R @@ -1,8 +1,8 @@ - #------------------------------------------------------------------------------ +#------------------------------------------------------------------------------ # BUILDING AND INITIALIZING SEIR MODEL # library(epiworldR) -# sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, +# sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, # recovery_rate = 0.1) @@ -37,15 +37,15 @@ plot_epi.epiworld_model <- function( x, main = "", counts_scale, ... -) { - + ) { + plot_epi( x = get_hist_total(x), main = main, counts_scale = counts_scale, ... ) - + } #' @export @@ -53,15 +53,15 @@ plot_epi.epiworld_hist_virus <- function( x, main = "", counts_scale, ... -) { - - res <- lapply(sort(unique(x$id)), function(i) x[x$id == i,]) - + ) { + + res <- lapply(sort(unique(x$id)), function(i) x[x$id == i, ]) + lapply(res, function(r) { plot_epi.epiworld_hist(r, main = paste0("Variant id ", r$id[1])) - }) + }) invisible(x) - + } #' @export @@ -69,84 +69,83 @@ plot_epi.epiworld_hist <- function( x, main = "", counts_scale, ... - ) { - + ) { + curves <- x state_names <- attr(curves, "states") - + # If the user didn't say what scale if (missing(counts_scale)) counts_scale <- find_scale(max(curves$counts)) - - curves$counts <- curves$counts/counts_scale - + + curves$counts <- curves$counts / counts_scale + # Initialize date vector of size length for state names - date_candidates <- integer(length = length(state_names)) + date_candidates <- integer(length = length(state_names)) # Identify max date when the counts stop significantly changing by state - - benchmark_value <- diff(range(curves$counts))/200 # 0.5% of range - + + benchmark_value <- diff(range(curves$counts)) / 200 # 0.5% of range + for (i in 1L:length(state_names)) { date_candidates[i] <- with( - curves[curves$state == state_names[i],], - sum(abs(diff(counts)) > benchmark_value ) - ) + curves[curves$state == state_names[i], ], + sum(abs(diff(counts)) > benchmark_value) + ) } - # Round the maximum date up to the nearest 10th + # Round the maximum date up to the nearest 10th max_date <- min( diff(range(curves$date)), max(ceiling(max(date_candidates) / 10L) * 10L, 10L) ) - + # Defining range of x values by max date as the max - curves <- curves[curves$date < max_date,] - # Defining range of y values + curves <- curves[curves$date < max_date, ] + # Defining range of y values counts_range <- range(curves$counts) # Plot the first state with( - curves[curves$state == state_names[1L],], + curves[curves$state == state_names[1L], ], graphics::plot( x = date, - y = counts, - type = 'l', + y = counts, + type = "l", col = 1, - ylim = counts_range, - xlab = "Day (step)", + ylim = counts_range, + xlab = "Day (step)", ylab = ifelse( counts_scale == 1L, "Population", paste("Population (", counts_scale, "'s)", sep = "") - ), + ), main = main ) ) - + # Plot the remaining states for (i in 2L:length(state_names)) { - + with( - curves[curves$state == state_names[i],], + curves[curves$state == state_names[i], ], graphics::lines( x = date, y = counts, - type = 'l', + type = "l", col = i - ) + ) ) - + } - + # Legend graphics::legend( "right", legend = state_names, col = 1L:length(state_names), - lty = 1L, + lty = 1L, lwd = 2L, bty = "n" - ) + ) } # plot_epi(sir, main = "SIR Model") - diff --git a/R/tool.R b/R/tool.R index d311891e..ff86257d 100644 --- a/R/tool.R +++ b/R/tool.R @@ -1,9 +1,9 @@ #' Tools in epiworld -#' +#' #' Tools are functions that affect how agents react to the virus. They can be #' used to simulate the effects of vaccination, isolation, and social #' distancing. -#' +#' #' @param model Model #' @param name Name of the tool #' @param susceptibility_reduction Numeric. Proportion it reduces susceptibility. @@ -11,7 +11,7 @@ #' @param recovery_enhancer Numeric. Proportion it improves recovery. #' @param death_reduction Numeric. Proportion it reduces probability of death.e #' @param tool_pos Positive integer. Index of the tool's position in the model. -#' @examples +#' @examples #' # Simple model #' model_sirconn <- ModelSIRCONN( #' name = "COVID-19", @@ -21,49 +21,49 @@ #' transmission_rate = 0.4, #' recovery_rate = 0.95 #' ) -#' +#' #' # Running and printing #' run(model_sirconn, ndays = 100, seed = 1912) #' plot(model_sirconn) -#' +#' #' epitool <- tool( #' name = "Vaccine", #' prevalence = 0.5, #' as_proportion = TRUE, #' susceptibility_reduction = .9, #' transmission_reduction = .5, -#' recovery_enhancer = .5, +#' recovery_enhancer = .5, #' death_reduction = .9 #' ) -#' +#' #' epitool -#' -#' set_name_tool(epitool, 'Pfizer') # Assigning name to the tool +#' +#' set_name_tool(epitool, "Pfizer") # Assigning name to the tool #' get_name_tool(epitool) # Returning the name of the tool #' add_tool(model_sirconn, epitool) #' run(model_sirconn, ndays = 100, seed = 1912) #' model_sirconn #' plot(model_sirconn) -#' +#' #' # To declare a certain number of individuals with the tool #' rm_tool(model_sirconn, 0) # Removing epitool from the model #' # Setting prevalence to 0.1 #' set_distribution_tool(epitool, distribute_tool_randomly(0.1, TRUE)) #' add_tool(model_sirconn, epitool) #' run(model_sirconn, ndays = 100, seed = 1912) -#' +#' #' # Adjusting probabilities due to tool -#' set_susceptibility_reduction(epitool, 0.1) # Susceptibility reduction +#' set_susceptibility_reduction(epitool, 0.1) # Susceptibility reduction #' set_transmission_reduction(epitool, 0.2) # Transmission reduction #' set_recovery_enhancer(epitool, 0.15) # Probability increase of recovery #' set_death_reduction(epitool, 0.05) # Probability reduction of death -#' -#' rm_tool(model_sirconn, 0) +#' +#' rm_tool(model_sirconn, 0) #' add_tool(model_sirconn, epitool) #' run(model_sirconn, ndays = 100, seed = 1912) # Run model to view changes -#' +#' #' @export -#' @returns +#' @returns #' - The `tool` function creates a tool of class [epiworld_tool]. #' @aliases epiworld_tool tool <- function( @@ -74,7 +74,7 @@ tool <- function( transmission_reduction, recovery_enhancer, death_reduction -) { + ) { uses_deprecated <- FALSE if (missing(prevalence)) { @@ -82,7 +82,7 @@ tool <- function( warning( "Starting version 0.3-0, the 'prevalence' argument is required.", " It will be set to be 0.5. Next versions will fail with an error." - ) + ) prevalence <- 0.5 as_proportion <- TRUE @@ -107,7 +107,7 @@ tool <- function( as_proportion = as_proportion ) ) - + } #' @export @@ -120,7 +120,7 @@ stopifnot_tool <- function(tool) { if (!inherits(tool, "epiworld_tool")) { stop( "The -tool- object must be of class \"epiworld_tool\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(tool), collapse = ", ") ) } @@ -130,7 +130,7 @@ stopifnot_tfun <- function(tfun) { if (!inherits(tfun, "epiworld_tool_fun")) { stop( "The -tfun- object must be of class \"epiworld_tool_fun\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(tfun), collapse = ", ") ) } @@ -140,7 +140,7 @@ stopifnot_tool_distfun <- function(tool_distfun) { if (!inherits(tool_distfun, "epiworld_tool_distfun")) { stop( "The -tool_distfun- object must be of class \"epiworld_tool_distfun\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(tool_distfun), collapse = ", ") ) } @@ -150,8 +150,8 @@ stopifnot_tool_distfun <- function(tool_distfun) { #' @details #' The name of the `epiworld_tool` object can be manipulated with the functions #' [set_name_tool()] and [get_name_tool()]. -#' @returns -#' - The `set_name_tool` function assigns a name to the tool of class +#' @returns +#' - The `set_name_tool` function assigns a name to the tool of class #' [epiworld_tool] and returns the tool. #' @rdname tool set_name_tool <- function(tool, name) { @@ -161,7 +161,7 @@ set_name_tool <- function(tool, name) { #' @returns -#' - The `get_name_tool` function returns the name of the tool of class +#' - The `get_name_tool` function returns the name of the tool of class #' [epiworld_tool]. #' @rdname tool #' @export @@ -173,18 +173,18 @@ get_name_tool <- function(tool) { #' @export #' @param tool An object of class `epiworld_tool` #' @param proportion Deprecated. -#' @details -#' The `add_tool` function adds the specified tool to the model of class +#' @details +#' The `add_tool` function adds the specified tool to the model of class #' [epiworld_model] with specified proportion. #' @rdname tool add_tool <- function(model, tool, proportion) { - + if (!missing(proportion)) { warning( "The 'proportion' argument is deprecated. ", "Use 'set_distribution_tool' instead." - ) + ) set_distribution_tool( tool = tool, @@ -217,7 +217,7 @@ add_tool.epiworld_model <- function(model, tool, proportion) { } #' @export -#' @returns +#' @returns #' - The `rm_tool` function removes the specified tool from a model. #' @rdname tool rm_tool <- function(model, tool_pos) { @@ -236,13 +236,13 @@ rm_tool <- function(model, tool_pos) { #' coefficients associated to the logit probability. #' @rdname tool #' @examples -#' +#' #' # Using the logit function -------------- #' sir <- ModelSIR( -#' name = "COVID-19", prevalence = 0.01, +#' name = "COVID-19", prevalence = 0.01, #' transmission_rate = 0.9, recovery_rate = 0.1 -#' ) -#' +#' ) +#' #' # Adding a small world population #' agents_smallworld( #' sir, @@ -251,7 +251,7 @@ rm_tool <- function(model, tool_pos) { #' d = FALSE, #' p = .01 #' ) -#' +#' #' # Creating a tool #' mask_wearing <- tool( #' name = "Mask", @@ -262,52 +262,53 @@ rm_tool <- function(model, tool_pos) { #' recovery_enhancer = 0.0, #' death_reduction = 0.0 #' ) -#' +#' #' add_tool(sir, mask_wearing) -#' +#' #' run(sir, ndays = 50, seed = 11) #' hist_0 <- get_hist_total(sir) -#' +#' #' # And adding features #' dat <- cbind( #' female = sample.int(2, 10000, replace = TRUE) - 1, #' x = rnorm(10000) #' ) -#' +#' #' set_agents_data(sir, dat) -#' +#' #' # Creating the logit function #' tfun <- tool_fun_logit( #' vars = c(0L, 1L), #' coefs = c(-1, 1), #' model = sir #' ) -#' +#' #' # The infection prob is lower -#' hist(plogis(dat %*% rbind(.5,1))) -#' +#' hist(plogis(dat %*% rbind(.5, 1))) +#' #' tfun # printing -#' -#' +#' +#' #' set_susceptibility_reduction_fun( #' tool = get_tool(sir, 0), #' model = sir, #' tfun = tfun -#' ) -#' +#' ) +#' #' run(sir, ndays = 50, seed = 11) #' hist_1 <- get_hist_total(sir) -#' +#' #' op <- par(mfrow = c(1, 2)) -#' plot(hist_0); abline(v = 30) -#' plot(hist_1); abline(v = 30) +#' plot(hist_0) +#' abline(v = 30) +#' plot(hist_1) +#' abline(v = 30) #' par(op) -#' -#' +#' tool_fun_logit <- function(vars, coefs, model) { - + stopifnot_model(model) - + structure( tool_fun_logit_cpp(as.integer(vars), as.double(coefs), model), class = "epiworld_tool_fun", @@ -316,12 +317,12 @@ tool_fun_logit <- function(vars, coefs, model) { coefs = coefs, model = model ) - + } #' @export print.epiworld_tool_fun <- function(x, ...) { - + cat("An epiworld_tool_function object.\n") cat("(model: ", get_name(attr(x, "model")), ")\n", sep = "") cat("This function was built using -tool_fun_logit()-. and it features ") @@ -334,9 +335,9 @@ print.epiworld_tool_fun <- function(x, ...) { ), collapse = "\n" ), "\n" ) - + invisible(x) - + } # Susceptibility reduction ----------------------------------------------------- @@ -345,75 +346,75 @@ print.epiworld_tool_fun <- function(x, ...) { #' @export #' @param prob Numeric scalar. A probability (between zero and one). #' @returns -#' - The `set_susceptibility_reduction` function assigns a probability reduction -#' to the specified tool of class [epiworld_tool]. +#' - The `set_susceptibility_reduction` function assigns a probability reduction +#' to the specified tool of class [epiworld_tool]. #' @rdname tool set_susceptibility_reduction <- function(tool, prob) { - + stopifnot_tool(tool) set_susceptibility_reduction_cpp(tool, as.double(prob)) - + } #' @export #' @param param Character scalar. Name of the parameter featured in `model` that #' will be added to the tool (see details). #' @details -#' In the case of `set_susceptibility_reduction_ptr`, `set_transmission_reduction_ptr`, +#' In the case of `set_susceptibility_reduction_ptr`, `set_transmission_reduction_ptr`, #' `set_recovery_enhancer`, and #' `set_death_reduction_ptr`, the corresponding parameters are passed as a pointer to #' the tool. The implication of using pointers is that the values will be #' read directly from the `model` object, so changes will be reflected. -#' +#' #' @rdname tool set_susceptibility_reduction_ptr <- function(tool, model, param) { - + stopifnot_tool(tool) stopifnot_model(model) invisible(set_susceptibility_reduction_ptr_cpp(tool, model, param)) - + } #' @export #' @param tfun An object of class `epiworld_tool_fun`. #' @rdname tool set_susceptibility_reduction_fun <- function(tool, model, tfun) { - + stopifnot_tool(tool) stopifnot_model(model) stopifnot_tfun(tfun) invisible(set_susceptibility_reduction_fun_cpp(tool, model, tfun)) - + } # Transmission reduction ------------------------------------------------------- #' @export #' @returns -#' - The `set_transmission_reduction` function assigns a probability reduction -#' to the specified tool of class [epiworld_tool]. +#' - The `set_transmission_reduction` function assigns a probability reduction +#' to the specified tool of class [epiworld_tool]. #' @rdname tool set_transmission_reduction <- function(tool, prob) { - + stopifnot_tool(tool) invisible(set_transmission_reduction_cpp(tool, as.double(prob))) - + } #' @export #' @rdname tool set_transmission_reduction_ptr <- function(tool, model, param) { - + stopifnot_tool(tool) stopifnot_model(model) invisible(set_transmission_reduction_ptr_cpp(tool, model, param)) - + } #' @export #' @rdname tool set_transmission_reduction_fun <- function(tool, model, tfun) { - + stopifnot_tool(tool) stopifnot_model(model) stopifnot_tfun(tfun) @@ -424,111 +425,111 @@ set_transmission_reduction_fun <- function(tool, model, tfun) { #' @export #' @returns -#' - The `set_recovery_enhancer` function assigns a probability increase -#' to the specified tool of class [epiworld_tool]. +#' - The `set_recovery_enhancer` function assigns a probability increase +#' to the specified tool of class [epiworld_tool]. #' @rdname tool set_recovery_enhancer <- function(tool, prob) { - + stopifnot_tool(tool) invisible(set_recovery_enhancer_cpp(tool, as.double(prob))) - + } #' @export #' @rdname tool set_recovery_enhancer_ptr <- function(tool, model, param) { - + stopifnot_tool(tool) stopifnot_model(model) invisible(set_recovery_enhancer_ptr_cpp(tool, model, param)) - + } #' @export #' @rdname tool set_recovery_enhancer_fun <- function(tool, model, tfun) { - + stopifnot_tool(tool) stopifnot_model(model) stopifnot_tfun(tfun) invisible(set_recovery_enhancer_fun_cpp(tool, model, tfun)) - + } # Death reduction -------------------------------------------------------------- #' @export #' @returns -#' - The `set_death_reduction` function assigns a probability decrease -#' to the specified tool of class [epiworld_tool]. +#' - The `set_death_reduction` function assigns a probability decrease +#' to the specified tool of class [epiworld_tool]. #' @rdname tool set_death_reduction <- function(tool, prob) { - + stopifnot_tool(tool) invisible(set_death_reduction_cpp(tool, as.double(prob))) - + } #' @export #' @rdname tool set_death_reduction_ptr <- function(tool, model, param) { - + stopifnot_tool(tool) stopifnot_model(model) invisible(set_death_reduction_ptr_cpp(tool, model, param)) - + } #' @export #' @rdname tool set_death_reduction_fun <- function(tool, model, tfun) { - + stopifnot_tool(tool) stopifnot_model(model) stopifnot_tfun(tfun) invisible(set_death_reduction_fun_cpp(tool, model, tfun)) - + } #' @export #' @rdname agents_smallworld -#' @returns +#' @returns #' - `get_agents_tools` returns a list of class `epiworld_agents_tools` #' with `epiworld_tools` (list of lists). get_agents_tools <- function(model) { - + stopifnot_model(model) - + res <- lapply( get_agents_tools_cpp(model), `class<-`, "epiworld_tools" ) - + structure(res, class = c("epiworld_agents_tools", class(res))) - + } -#' @export +#' @export #' @rdname tool #' @param max_print Numeric scalar. Maximum number of tools to print. #' @param ... Currently ignored. #' @param x An object of class `epiworld_agents_tools`. print.epiworld_agents_tools <- function(x, max_print = 10, ...) { - + for (i in 1:min(max_print, length(x))) { print_agent_tools_cpp(x[[i]]) } - + if (length(x) > max_print) { cat(sprintf("Showing first %s of %s tools.\n", max_print, length(x))) } - + invisible(x) - + } -#' @export +#' @export #' @details #' The `set_distribution_tool` function assigns a distribution function to the #' specified tool of class [epiworld_tool]. The distribution function can be @@ -546,19 +547,19 @@ set_distribution_tool <- function(tool, distfun) { #' @export #' @rdname tool -#' @details +#' @details #' The `distribute_tool_randomly` function creates a distribution function that #' randomly assigns the tool to a proportion of the population. #' @param as_proportion Logical scalar. If `TRUE`, `prevalence` is interpreted #' as a proportion of the total number of agents in the model. #' @param prevalence Numeric scalar. Prevalence of the tool. -#' @return +#' @return #' - The `distribute_tool_randomly` function returns a distribution function of #' class `epiworld_tool_distfun`. distribute_tool_randomly <- function( - prevalence, - as_proportion -) { + prevalence, + as_proportion + ) { structure( distribute_tool_randomly_cpp( @@ -572,7 +573,7 @@ distribute_tool_randomly <- function( #' @export #' @rdname tool -#' @details +#' @details #' The `distribute_tool_to_set` function creates a distribution function that #' assigns the tool to a set of agents. #' @param agents_ids Integer vector. Indices of the agents to which the tool @@ -581,8 +582,8 @@ distribute_tool_randomly <- function( #' - The `distribute_tool_to_set` function returns a distribution function of #' class `epiworld_tool_distfun`. distribute_tool_to_set <- function( - agents_ids -) { + agents_ids + ) { structure( distribute_tool_to_set_cpp( @@ -591,4 +592,4 @@ distribute_tool_to_set <- function( class = "epiworld_tool_distfun" ) -} \ No newline at end of file +} diff --git a/R/virus.R b/R/virus.R index b93b13d4..6ffa7c0e 100644 --- a/R/virus.R +++ b/R/virus.R @@ -1,9 +1,9 @@ #' Virus design -#' +#' #' Viruses can be considered to be anything that can be transmitted (e.g., -#' diseases, as well as ideas.) Most models in epiworldR can feature multiple +#' diseases, as well as ideas.) Most models in epiworldR can feature multiple #' viruses. -#' +#' #' @param name of the virus #' @param post_immunity Numeric scalar. Post immunity (prob of re-infection). #' @param prob_infecting Numeric scalar. Probability of infection (transmission). @@ -14,44 +14,44 @@ #' @details #' The [virus()] function can be used to initialize a virus. Virus features can #' then be modified using the functions `set_prob_*`. -#' +#' #' The function [virus_fun_logit()] creates a "virus function" that can be #' evaluated for transmission, recovery, and death. As the name sugests, it #' computes those probabilities using a logit function (see examples). -#' -#' @examples +#' +#' @examples #' mseirconn <- ModelSEIRCONN( #' name = "COVID-19", -#' prevalence = 0.01, +#' prevalence = 0.01, #' n = 10000, -#' contact_rate = 4, -#' incubation_days = 7, +#' contact_rate = 4, +#' incubation_days = 7, #' transmission_rate = 0.5, #' recovery_rate = 0.99 #' ) -#' +#' #' delta <- virus( #' "Delta Variant", 0, .5, .2, .01, prevalence = 0.3, as_proportion = TRUE #' ) -#' +#' #' # Adding virus and setting/getting virus name #' add_virus(mseirconn, delta) #' set_name_virus(delta, "COVID-19 Strain") #' get_name_virus(delta) -#' +#' #' run(mseirconn, ndays = 100, seed = 992) #' mseirconn -#' +#' #' rm_virus(mseirconn, 0) # Removing the first virus from the model object #' set_distribution_virus(delta, distribute_virus_randomly(100, as_proportion = FALSE)) -#' add_virus(mseirconn, delta) -#' +#' add_virus(mseirconn, delta) +#' #' # Setting parameters for the delta virus manually #' set_prob_infecting(delta, 0.5) #' set_prob_recovery(delta, 0.9) #' set_prob_death(delta, 0.01) #' run(mseirconn, ndays = 100, seed = 992) # Run the model to observe changes -#' +#' #' # If the states were (for example): #' # 1: Infected #' # 2: Recovered @@ -63,15 +63,15 @@ #' @export #' @aliases epiworld_virus virus <- function( - name, - prevalence, - as_proportion, - prob_infecting, - recovery_rate = 0.5, - prob_death = 0.0, - post_immunity = -1.0, - incubation = 7.0 - ) { + name, + prevalence, + as_proportion, + prob_infecting, + recovery_rate = 0.5, + prob_death = 0.0, + post_immunity = -1.0, + incubation = 7.0 + ) { uses_deprecated <- FALSE if (missing(prevalence)) { @@ -79,14 +79,14 @@ virus <- function( warning( "Starting version 0.3-0, the 'prevalence' argument is required.", " It will be set to be 0.5. Next versions will fail with an error." - ) + ) prevalence <- 0.5 as_proportion <- TRUE uses_deprecated <- TRUE } - + structure( virus_cpp( name, @@ -97,7 +97,7 @@ virus <- function( prob_death, post_immunity, incubation - ), + ), class = "epiworld_virus", uses_deprecated = uses_deprecated, deprecated_args = list( @@ -105,7 +105,7 @@ virus <- function( as_proportion = as_proportion ) ) - + } #' @export @@ -117,7 +117,7 @@ stopifnot_virus <- function(virus) { if (!inherits(virus, "epiworld_virus")) { stop( "The -virus- object must be of class \"epiworld_virus\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(virus), collapse = ", ") ) } @@ -127,7 +127,7 @@ stopifnot_vfun <- function(vfun) { if (!inherits(vfun, "epiworld_virus_fun")) { stop( "The -vfun- object must be of class \"epiworld_virus_fun\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(vfun), collapse = ", ") ) } @@ -137,7 +137,7 @@ stopifnot_virus_distfun <- function(virus_distfun) { if (!inherits(virus_distfun, "epiworld_virus_distfun")) { stop( "The -virus_distfun- object must be of class \"epiworld_virus_distfun\". ", - "The object passed to the function is of class(es): ", + "The object passed to the function is of class(es): ", paste(class(virus_distfun), collapse = ", ") ) } @@ -148,9 +148,9 @@ stopifnot_virus_distfun <- function(virus_distfun) { #' @details #' The name of the `epiworld_virus` object can be manipulated with the functions #' [set_name_virus()] and [get_name_virus()]. -#' @returns -#' - The `set_name_virus` function does not return a value, but merely assigns -#' a name to the virus of choice. +#' @returns +#' - The `set_name_virus` function does not return a value, but merely assigns +#' a name to the virus of choice. #' @rdname virus set_name_virus <- function(virus, name) { stopifnot_virus(virus) @@ -158,24 +158,24 @@ set_name_virus <- function(virus, name) { } #' @export -#' @returns -#' - The `get_name_virus` function returns the name of the virus of class +#' @returns +#' - The `get_name_virus` function returns the name of the virus of class #' [epiworld_virus]. #' @rdname virus get_name_virus <- function(virus) { stopifnot_virus(virus) get_name_virus_cpp(virus) } - + # Virus add -------------------------------------------------------------------- #' @export #' @rdname virus #' @param model An object of class `epiworld_model`. #' @param virus An object of class `epiworld_virus` -#' @param proportion Deprecated. -#' @returns -#' - The `add_virus` function does not return a value, instead it adds the +#' @param proportion Deprecated. +#' @returns +#' - The `add_virus` function does not return a value, instead it adds the #' virus of choice to the model object of class [epiworld_model]. add_virus <- function(model, virus, proportion) { @@ -184,11 +184,11 @@ add_virus <- function(model, virus, proportion) { warning( "The argument 'proportion' is deprecated and will be removed in ", "the next version." - ) + ) set_distribution_virus( - virus=virus, - distfun=distribute_virus_randomly(proportion, as_proportion = TRUE) + virus = virus, + distfun = distribute_virus_randomly(proportion, as_proportion = TRUE) ) } else if (isTRUE(attr(tool, "uses_deprecated"))) { @@ -209,80 +209,80 @@ add_virus <- function(model, virus, proportion) { #' @export add_virus.epiworld_model <- function(model, virus, proportion) { - + stopifnot_virus(virus) add_virus_cpp(model, virus) invisible(model) - + } #' @export add_virus.epiworld_sir <- function(model, virus, proportion) { - + stopifnot_virus(virus) virus_set_state(virus, init = 1, end = 2, removed = 2) invisible(add_virus_cpp(model, virus)) - + } #' @export add_virus.epiworld_sird <- function(model, virus, proportion) { - + stopifnot_virus(virus) virus_set_state(virus, init = 1, end = 2, removed = 3) invisible(add_virus_cpp(model, virus)) - + } #' @export add_virus.epiworld_sirconn <- function(model, virus, proportion) { - + stopifnot_virus(virus) add_virus.epiworld_sir(model, virus) - + } #' @export add_virus.epiworld_sirdconn <- function(model, virus, proportion) { - + stopifnot_virus(virus) add_virus.epiworld_sird(model, virus) - + } #' @export add_virus.epiworld_seir <- function(model, virus, proportion) { - + stopifnot_virus(virus) virus_set_state(virus, init = 1, end = 3, removed = 3) invisible(add_virus_cpp(model, virus)) - + } #' @export add_virus.epiworld_seird <- function(model, virus, proportion) { - + stopifnot_virus(virus) virus_set_state(virus, init = 1, end = 3, removed = 4) invisible(add_virus_cpp(model, virus)) - + } #' @export add_virus.epiworld_seirconn <- function(model, virus, proportion) { - + stopifnot_virus(virus) add_virus.epiworld_seir(model, virus) - + } #' @export add_virus.epiworld_seirdconn <- function(model, virus, proportion) { - + stopifnot_virus(virus) add_virus.epiworld_seird(model, virus) - + } # Virus MISC ------------------------------------------------------------------- @@ -291,26 +291,26 @@ add_virus.epiworld_seirdconn <- function(model, virus, proportion) { #' @rdname virus #' @param init,end,removed states after acquiring a virus, removing a virus, #' and removing the agent as a result of the virus, respectively. -#' @returns -#' - The `virus_set_state` function does not return a value but assigns +#' @returns +#' - The `virus_set_state` function does not return a value but assigns #' epidemiological properties to the specified virus of class [epiworld_virus]. virus_set_state <- function(virus, init, end, removed) { - + stopifnot_virus(virus) invisible(virus_set_state_cpp(virus, init, end, removed)) - + } #' @export -#' @returns -#' - The `rm_virus` function does not return a value, but instead removes +#' @returns +#' - The `rm_virus` function does not return a value, but instead removes #' a specified virus from the model of class [epiworld_model]. #' @rdname virus rm_virus <- function(model, virus_pos) { - + stopifnot_model(model) invisible(rm_virus_cpp(model, virus_pos)) - + } # Virus functions -------------------------------------------------------------- @@ -324,10 +324,10 @@ rm_virus <- function(model, virus_pos) { #' @examples #' # Using the logit function -------------- #' sir <- ModelSIR( -#' name = "COVID-19", prevalence = 0.01, +#' name = "COVID-19", prevalence = 0.01, #' transmission_rate = 0.9, recovery = 0.1 -#' ) -#' +#' ) +#' #' # Adding a small world population #' agents_smallworld( #' sir, @@ -336,44 +336,43 @@ rm_virus <- function(model, virus_pos) { #' d = FALSE, #' p = .01 #' ) -#' +#' #' run(sir, ndays = 50, seed = 11) #' plot(sir) -#' +#' #' # And adding features #' dat <- cbind( #' female = sample.int(2, 10000, replace = TRUE) - 1, #' x = rnorm(10000) #' ) -#' +#' #' set_agents_data(sir, dat) -#' +#' #' # Creating the logit function #' vfun <- virus_fun_logit( #' vars = c(0L, 1L), #' coefs = c(-1, 1), #' model = sir #' ) -#' +#' #' # The infection prob is lower -#' hist(plogis(dat %*% rbind(-1,1))) -#' +#' hist(plogis(dat %*% rbind(-1, 1))) +#' #' vfun # printing -#' +#' #' set_prob_infecting_fun( #' virus = get_virus(sir, 0), #' model = sir, #' vfun = vfun -#' ) -#' +#' ) +#' #' run(sir, ndays = 50, seed = 11) #' plot(sir) -#' -#' +#' virus_fun_logit <- function(vars, coefs, model) { - + stopifnot_model(model) - + structure( virus_fun_logit_cpp(vars, coefs, model), class = "epiworld_virus_fun", @@ -382,12 +381,12 @@ virus_fun_logit <- function(vars, coefs, model) { coefs = coefs, model = model ) - + } #' @export print.epiworld_virus_fun <- function(x, ...) { - + cat("An epiworld_virus_function object.\n") cat("(model: ", get_name(attr(x, "model")), ")\n", sep = "") cat("This function was built using -virus_fun_logit()-. and it features ") @@ -397,27 +396,27 @@ print.epiworld_virus_fun <- function(x, ...) { " % 2i: %5.2f", attr(x, "vars"), attr(x, "coefs") - ), collapse = "\n" + ), collapse = "\n" ), "\n" ) - + invisible(x) - + } #' @export #' @param prob Numeric scalar. A probability (between zero and one). -#' @returns -#' - The `set_prob_infecting` function does not return a value, but instead -#' assigns a probability to infection for the specified virus of class +#' @returns +#' - The `set_prob_infecting` function does not return a value, but instead +#' assigns a probability to infection for the specified virus of class #' [epiworld_virus]. #' @rdname virus set_prob_infecting <- function(virus, prob) { - + stopifnot_virus(virus) invisible(set_prob_infecting_cpp(virus, as.numeric(prob))) - + } #' @export @@ -430,52 +429,52 @@ set_prob_infecting <- function(virus, prob) { #' read directly from the `model` object, so changes will be reflected. #' @rdname virus set_prob_infecting_ptr <- function(virus, model, param) { - + stopifnot_virus(virus) stopifnot_model(model) invisible(set_prob_infecting_ptr_cpp(virus, model, param)) - + } #' @export #' @param vfun An object of class `epiworld_virus_fun`. #' @rdname virus set_prob_infecting_fun <- function(virus, model, vfun) { - + stopifnot_virus(virus) stopifnot_model(model) stopifnot_vfun(vfun) invisible(set_prob_infecting_fun_cpp(virus, model, vfun)) - + } #' @export -#' @returns -#' - The `set_prob_recovery` function does not return a value, but instead -#' assigns a probability to recovery for the specified virus of class +#' @returns +#' - The `set_prob_recovery` function does not return a value, but instead +#' assigns a probability to recovery for the specified virus of class #' [epiworld_virus]. #' @rdname virus set_prob_recovery <- function(virus, prob) { - + stopifnot_virus(virus) invisible(set_prob_recovery_cpp(virus, as.numeric(prob))) - + } #' @export #' @rdname virus set_prob_recovery_ptr <- function(virus, model, param) { - + stopifnot_virus(virus) stopifnot_model(model) invisible(set_prob_recovery_ptr_cpp(virus, model, param)) - + } #' @export #' @rdname virus set_prob_recovery_fun <- function(virus, model, vfun) { - + stopifnot_virus(virus) stopifnot_model(model) stopifnot_vfun(vfun) @@ -484,81 +483,81 @@ set_prob_recovery_fun <- function(virus, model, vfun) { } #' @export -#' @returns -#' - The `set_prob_death` function does not return a value, but instead -#' assigns a probability to death for the specified virus of class +#' @returns +#' - The `set_prob_death` function does not return a value, but instead +#' assigns a probability to death for the specified virus of class #' [epiworld_virus]. #' @rdname virus set_prob_death <- function(virus, prob) { - + stopifnot_virus(virus) invisible(set_prob_death_cpp(virus, as.numeric(prob))) - + } #' @export #' @rdname virus set_prob_death_ptr <- function(virus, model, param) { - + stopifnot_virus(virus) stopifnot_model(model) invisible(set_prob_death_ptr_cpp(virus, model, param)) - + } #' @export #' @rdname virus set_prob_death_fun <- function(virus, model, vfun) { - + stopifnot_virus(virus) stopifnot_model(model) stopifnot_vfun(vfun) invisible(set_prob_death_fun_cpp(virus, model, vfun)) - + } -#' @export -#' @return +#' @export +#' @return #' - The `set_incubation` function does not return a value, but instead #' assigns an incubation period to the specified virus of class [epiworld_virus]. #' @rdname virus set_incubation <- function(virus, incubation) { - + stopifnot_virus(virus) invisible(set_incubation_cpp(virus, as.numeric(incubation))) - + } #' @export #' @rdname virus set_incubation_ptr <- function(virus, model, param) { - + stopifnot_virus(virus) stopifnot_model(model) invisible(set_incubation_ptr_cpp(virus, model, param)) - + } #' @export #' @rdname virus set_incubation_fun <- function(virus, model, vfun) { - + stopifnot_virus(virus) stopifnot_model(model) stopifnot_vfun(vfun) invisible(set_incubation_fun_cpp(virus, model, vfun)) - + } #' @export #' @rdname virus #' @param distfun An object of class `epiworld_distribution_virus`. set_distribution_virus <- function(virus, distfun) { - + stopifnot_virus(virus) stopifnot_virus_distfun(distfun) invisible(set_distribution_virus_cpp(virus, distfun)) - + } #' @export @@ -570,13 +569,13 @@ set_distribution_virus <- function(virus, distfun) { #' @param prevalence Numeric scalar. Prevalence of the virus. #' @param as_proportion Logical scalar. If `TRUE`, the prevalence is set as a #' proportion of the total number of agents in the model. -#' @return +#' @return #' - The `distribute_virus_randomly` function returns a function that can be #' used to distribute the virus in the model. distribute_virus_randomly <- function( - prevalence, - as_proportion -) { + prevalence, + as_proportion + ) { structure( distribute_virus_randomly_cpp( @@ -593,10 +592,10 @@ distribute_virus_randomly <- function( #' @param agents_ids Integer vector. Indices of the agents that will receive the #' virus. distribute_virus_set <- function(agents_ids) { - + structure( distribute_virus_to_set_cpp(as.vector(agents_ids)), class = "epiworld_virus_distfun" ) - -} \ No newline at end of file + +} diff --git a/README.Rmd b/README.Rmd index ba6d4105..2d336627 100644 --- a/README.Rmd +++ b/README.Rmd @@ -13,7 +13,7 @@ knitr::opts_chunk$set( ) ``` -# epiworldR +# epiworldR [![CRAN status](https://www.r-pkg.org/badges/version/epiworldR)](https://CRAN.R-project.org/package=epiworldR) @@ -39,8 +39,8 @@ From the package's description: Current available models: ```{r print-models, echo=FALSE, results='asis'} -models <- list.files(path="R/", pattern = "Model.*.R", full.names = FALSE) |> - gsub(pattern = "(Model.*)\\.R", replacement = "\\1") +models <- list.files(path = "R/", pattern = "Model.*.R", full.names = FALSE) |> + gsub(pattern = "(Model.*)\\.R", replacement = "\\1") sprintf("%i. `%s`\n", 1:length(models), models) |> cat() @@ -80,8 +80,8 @@ sir <- ModelSIR( prevalence = .01, transmission_rate = .7, recovery = .3 - ) |> - # Adding a Small world population +) |> + # Adding a Small world population agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |> # Running the model for 50 days run(ndays = 50, seed = 1912) @@ -104,22 +104,22 @@ The SEIR model is similar to the SIR model but includes an exposed state. Here, ```{r seir-conn} model_seirconn <- ModelSEIRCONN( name = "COVID-19", - prevalence = 0.01, + prevalence = 0.01, n = 10000, - contact_rate = 10, - incubation_days = 7, + contact_rate = 10, + incubation_days = 7, transmission_rate = 0.1, - recovery_rate = 1/7 + recovery_rate = 1 / 7 ) |> add_virus( - virus( - name = "COVID-19", - prevalence = 0.01, - as_proportion = TRUE, - prob_infecting = 0.01, - recovery_rate = 0.6, - prob_death = 0.5, - incubation = 7 - )) + virus( + name = "COVID-19", + prevalence = 0.01, + as_proportion = TRUE, + prob_infecting = 0.01, + recovery_rate = 0.6, + prob_death = 0.5, + incubation = 7 +)) set.seed(132) run(model_seirconn, ndays = 100) @@ -141,9 +141,9 @@ head(plot_generation_time(model_seirconn)) ## SIR Logit -This model provides a more complex transmission and recovery pattern based +This model provides a more complex transmission and recovery pattern based on agents' features. With it, we can reflect co-morbidities that could change -the probability of infection and recovery. Here, we simulate a population +the probability of infection and recovery. Here, we simulate a population including a dataset with two features: an intercept and a binary variable `Female`. The probability of infection and recovery are functions of the intercept and the `Female` variables. The following code simulates a population of 100,000 agents @@ -160,7 +160,7 @@ n <- 100000 X <- cbind( Intercept = 1, Female = sample.int(2, n, replace = TRUE) - 1 - ) +) coef_infect <- c(.1, -2, 2) coef_recover <- rnorm(2) @@ -170,7 +170,7 @@ model_logit <- ModelSIRLogit( "covid2", data = X, coefs_infect = coef_infect, - coefs_recover = coef_recover, + coefs_recover = coef_recover, coef_infect_cols = 1L:ncol(X), coef_recover_cols = 1L:ncol(X), prob_infection = .8, @@ -191,7 +191,7 @@ rn <- get_reproductive_number(model_logit) (table( X[, "Female"], (1:n %in% rn$source) -) |> prop.table())[,2] +) |> prop.table())[, 2] # Looking into the agents get_agents(model_logit) @@ -208,11 +208,11 @@ sir <- ModelSIR( prevalence = .01, transmission_rate = .5, recovery = .5 - ) |> - # Adding a Small world population - agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |> - # Running the model for 50 days - run(ndays = 50, seed = 1912) +) |> + # Adding a Small world population + agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |> + # Running the model for 50 days + run(ndays = 50, seed = 1912) # Transmission network net <- get_transmissions(sir) @@ -221,10 +221,10 @@ net <- get_transmissions(sir) library(epiworldR) library(netplot) x <- igraph::graph_from_edgelist( - as.matrix(net[,2:3]) + 1 - ) + as.matrix(net[, 2:3]) + 1 +) -nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE) +nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex = TRUE) ``` ## Multiple simulations @@ -238,7 +238,7 @@ model_sir <- ModelSIRCONN( n = 1000, contact_rate = 2, transmission_rate = 0.9, recovery_rate = 0.1 - ) +) # Generating a saver saver <- make_saver("total_hist", "reproductive") @@ -295,5 +295,5 @@ You may want to check out other R packages for agent-based modeling: [`ABM`](htt [`RNetLogo`](https://cran.r-project.org/package=RNetLogo){target="_blank"}. ## Code of Conduct - + The epiworldR project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/1/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms. diff --git a/README.md b/README.md index 8c19d6c4..b2bc2132 100644 --- a/README.md +++ b/README.md @@ -105,8 +105,8 @@ sir <- ModelSIR( prevalence = .01, transmission_rate = .7, recovery = .3 - ) |> - # Adding a Small world population +) |> + # Adding a Small world population agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |> # Running the model for 50 days run(ndays = 50, seed = 1912) @@ -114,6 +114,9 @@ sir <- ModelSIR( #> |Running the model... #> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. #> | done. +``` + +``` r sir #> ________________________________________________________________________________ @@ -130,39 +133,42 @@ summary(sir) #> ________________________________________________________________________________ #> ________________________________________________________________________________ #> SIMULATION STUDY -#> +#> #> Name of the model : Susceptible-Infected-Recovered (SIR) #> Population size : 100000 #> Agents' data : (none) #> Number of entities : 0 #> Days (duration) : 50 (of 50) #> Number of viruses : 1 -#> Last run elapsed t : 134.00ms -#> Last run speed : 37.07 million agents x day / second +#> Last run elapsed t : 62.00ms +#> Last run speed : 79.47 million agents x day / second #> Rewiring : off -#> +#> #> Global events: #> (none) -#> +#> #> Virus(es): #> - COVID-19 -#> +#> #> Tool(s): #> (none) -#> +#> #> Model parameters: #> - Recovery rate : 0.3000 #> - Transmission rate : 0.7000 -#> +#> #> Distribution of the population at time 50: #> - (0) Susceptible : 99000 -> 822 #> - (1) Infected : 1000 -> 415 #> - (2) Recovered : 0 -> 98763 -#> +#> #> Transition Probabilities: #> - Susceptible 0.91 0.09 0.00 #> - Infected 0.00 0.70 0.30 #> - Recovered 0.00 0.00 1.00 +``` + +``` r plot(sir) ``` @@ -185,22 +191,22 @@ agents can transmit the disease to any other agent: ``` r model_seirconn <- ModelSEIRCONN( name = "COVID-19", - prevalence = 0.01, + prevalence = 0.01, n = 10000, - contact_rate = 10, - incubation_days = 7, + contact_rate = 10, + incubation_days = 7, transmission_rate = 0.1, - recovery_rate = 1/7 + recovery_rate = 1 / 7 ) |> add_virus( - virus( - name = "COVID-19", - prevalence = 0.01, - as_proportion = TRUE, - prob_infecting = 0.01, - recovery_rate = 0.6, - prob_death = 0.5, - incubation = 7 - )) + virus( + name = "COVID-19", + prevalence = 0.01, + as_proportion = TRUE, + prob_infecting = 0.01, + recovery_rate = 0.6, + prob_death = 0.5, + incubation = 7 +)) set.seed(132) run(model_seirconn, ndays = 100) @@ -208,47 +214,50 @@ run(model_seirconn, ndays = 100) #> Running the model... #> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. #> done. +``` + +``` r summary(model_seirconn) #> ________________________________________________________________________________ #> ________________________________________________________________________________ #> SIMULATION STUDY -#> +#> #> Name of the model : Susceptible-Exposed-Infected-Removed (SEIR) (connected) #> Population size : 10000 #> Agents' data : (none) #> Number of entities : 0 #> Days (duration) : 100 (of 100) #> Number of viruses : 2 -#> Last run elapsed t : 81.00ms -#> Last run speed : 12.21 million agents x day / second +#> Last run elapsed t : 15.00ms +#> Last run speed : 64.01 million agents x day / second #> Rewiring : off -#> +#> #> Global events: #> - Update infected individuals (runs daily) -#> +#> #> Virus(es): #> - COVID-19 #> - COVID-19 -#> +#> #> Tool(s): #> (none) -#> +#> #> Model parameters: #> - Avg. Incubation days : 7.0000 #> - Contact rate : 10.0000 #> - Prob. Recovery : 0.1429 #> - Prob. Transmission : 0.1000 -#> +#> #> Distribution of the population at time 100: -#> - (0) Susceptible : 9800 -> 13 +#> - (0) Susceptible : 9800 -> 11 #> - (1) Exposed : 200 -> 0 -#> - (2) Infected : 0 -> 1 +#> - (2) Infected : 0 -> 3 #> - (3) Recovered : 0 -> 9986 -#> +#> #> Transition Probabilities: #> - Susceptible 0.94 0.06 0.00 0.00 -#> - Exposed 0.00 0.86 0.14 0.00 -#> - Infected 0.00 0.00 0.85 0.15 +#> - Exposed 0.00 0.85 0.15 0.00 +#> - Infected 0.00 0.00 0.86 0.14 #> - Recovered 0.00 0.00 0.00 1.00 ``` @@ -269,25 +278,28 @@ head(plot(repnum)) - #> virus_id virus date avg n sd lb ub - #> 1 0 COVID-19 0 5.615385 91 4.832228 1.0 17.0 - #> 2 0 COVID-19 2 5.000000 9 3.605551 0.2 10.4 - #> 3 0 COVID-19 3 6.000000 13 5.049752 0.0 13.0 - #> 4 0 COVID-19 4 4.592593 27 3.885469 0.0 12.7 - #> 5 0 COVID-19 5 4.846154 26 4.920913 0.0 14.5 - #> 6 0 COVID-19 6 4.236842 38 3.241906 0.0 12.0 + #> virus_id virus date avg n sd lb ub + #> 1 0 COVID-19 0 5.769231 91 5.455022 1.000 20.750 + #> 2 0 COVID-19 2 6.400000 10 4.880801 0.450 14.875 + #> 3 0 COVID-19 3 5.166667 18 4.422536 0.425 13.000 + #> 4 0 COVID-19 4 4.659091 44 3.784566 0.000 12.850 + #> 5 0 COVID-19 5 5.205882 34 3.273210 0.000 12.175 + #> 6 0 COVID-19 6 3.137255 51 2.713077 0.000 8.750 - head(plot_generation_time(model_seirconn)) +``` r + +head(plot_generation_time(model_seirconn)) +``` #> date avg n sd ci_lower ci_upper virus virus_id - #> 1 2 7.125000 8 2.474874 2.7 9.825 COVID-19 0 - #> 2 3 8.090909 11 7.203534 2.0 23.750 COVID-19 0 - #> 3 4 6.708333 24 4.338695 2.0 16.425 COVID-19 0 - #> 4 5 7.428571 21 4.738897 2.0 15.500 COVID-19 0 - #> 5 6 7.628571 35 4.173345 2.0 15.300 COVID-19 0 - #> 6 7 6.921053 38 4.675304 2.0 16.150 COVID-19 0 + #> 1 2 4.444444 9 2.185813 2.2 8.000 COVID-19 0 + #> 2 3 7.411765 17 3.922034 2.4 15.000 COVID-19 0 + #> 3 4 8.538462 39 7.100208 2.0 22.000 COVID-19 0 + #> 4 5 6.312500 32 3.905641 2.0 13.225 COVID-19 0 + #> 5 6 7.200000 40 4.052223 2.0 15.100 COVID-19 0 + #> 6 7 7.660000 50 4.461216 2.0 17.000 COVID-19 0 ## SIR Logit @@ -311,7 +323,7 @@ n <- 100000 X <- cbind( Intercept = 1, Female = sample.int(2, n, replace = TRUE) - 1 - ) +) coef_infect <- c(.1, -2, 2) coef_recover <- rnorm(2) @@ -321,7 +333,7 @@ model_logit <- ModelSIRLogit( "covid2", data = X, coefs_infect = coef_infect, - coefs_recover = coef_recover, + coefs_recover = coef_recover, coef_infect_cols = 1L:ncol(X), coef_recover_cols = 1L:ncol(X), prob_infection = .8, @@ -338,6 +350,9 @@ run(model_logit, 50) #> |Running the model... #> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. #> | done. +``` + +``` r plot(model_logit) ``` @@ -351,9 +366,12 @@ rn <- get_reproductive_number(model_logit) (table( X[, "Female"], (1:n %in% rn$source) -) |> prop.table())[,2] -#> 0 1 +) |> prop.table())[, 2] +#> 0 1 #> 0.12984 0.14201 +``` + +``` r # Looking into the agents get_agents(model_logit) @@ -386,15 +404,18 @@ sir <- ModelSIR( prevalence = .01, transmission_rate = .5, recovery = .5 - ) |> - # Adding a Small world population - agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |> - # Running the model for 50 days - run(ndays = 50, seed = 1912) +) |> + # Adding a Small world population + agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |> + # Running the model for 50 days + run(ndays = 50, seed = 1912) #> _________________________________________________________________________ #> |Running the model... #> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. #> | done. +``` + +``` r # Transmission network net <- get_transmissions(sir) @@ -403,11 +424,14 @@ net <- get_transmissions(sir) library(epiworldR) library(netplot) #> Loading required package: grid +``` + +``` r x <- igraph::graph_from_edgelist( - as.matrix(net[,2:3]) + 1 - ) + as.matrix(net[, 2:3]) + 1 +) -nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE) +nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex = TRUE) ``` @@ -428,7 +452,7 @@ model_sir <- ModelSIRCONN( n = 1000, contact_rate = 2, transmission_rate = 0.9, recovery_rate = 0.1 - ) +) # Generating a saver saver <- make_saver("total_hist", "reproductive") @@ -436,11 +460,14 @@ saver <- make_saver("total_hist", "reproductive") # Running and printing # Notice the use of nthread = 2 to run the simulations in parallel run_multiple(model_sir, ndays = 100, nsims = 50, saver = saver, nthread = 2) -#> Starting multiple runs (50) using 2 thread(s) +#> Starting multiple runs (50) #> _________________________________________________________________________ #> _________________________________________________________________________ #> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. #> done. +``` + +``` r # Retrieving the results ans <- run_multiple_get_results(model_sir) @@ -450,17 +477,23 @@ head(ans$total_hist) #> 1 1 0 1 Susceptible 990 #> 2 1 0 1 Infected 10 #> 3 1 0 1 Recovered 0 -#> 4 1 1 1 Susceptible 977 -#> 5 1 1 1 Infected 22 -#> 6 1 1 1 Recovered 1 +#> 4 1 1 1 Susceptible 971 +#> 5 1 1 1 Infected 29 +#> 6 1 1 1 Recovered 0 +``` + +``` r head(ans$reproductive) #> sim_num virus_id virus source source_exposure_date rt -#> 1 1 0 COVID-19 976 9 0 -#> 2 1 0 COVID-19 644 9 0 -#> 3 1 0 COVID-19 608 9 0 -#> 4 1 0 COVID-19 314 9 0 -#> 5 1 0 COVID-19 41 9 0 -#> 6 1 0 COVID-19 32 9 0 +#> 1 1 0 COVID-19 943 9 0 +#> 2 1 0 COVID-19 40 9 0 +#> 3 1 0 COVID-19 6 9 0 +#> 4 1 0 COVID-19 973 8 0 +#> 5 1 0 COVID-19 495 9 0 +#> 6 1 0 COVID-19 480 8 0 +``` + +``` r plot(ans$reproductive) ``` @@ -482,16 +515,16 @@ If you use `epiworldR` in your research, please cite it as follows: ``` r citation("epiworldR") #> To cite epiworldR in publications use: -#> +#> #> Meyer, Derek and Vega Yon, George (2023). epiworldR: Fast Agent-Based #> Epi Models. Journal of Open Source Software, 8(90), 5781, #> https://doi.org/10.21105/joss.05781 -#> +#> #> And the actual R package: -#> -#> Meyer D, Vega Yon G (2024). _epiworldR: Fast Agent-Based Epi Models_. +#> +#> Meyer D, Vega Yon G (????). _epiworldR: Fast Agent-Based Epi Models_. #> R package version 0.3-2, . -#> +#> #> To see these entries in BibTeX format, use 'print(, #> bibtex=TRUE)', 'toBibtex(.)', or set #> 'options(citation.bibtex.max=999)'. diff --git a/_pkgdown.yml b/_pkgdown.yml index d71acfb9..b7517c2d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,4 +1,3 @@ url: ~ template: bootstrap: 5 - diff --git a/configure b/configure index 65317b69..9f0e850a 100755 --- a/configure +++ b/configure @@ -4535,4 +4535,3 @@ if test -n "$ac_unrecognized_opts" && test "$enable_option_checking" != no; then { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: unrecognized options: $ac_unrecognized_opts" >&5 printf "%s\n" "$as_me: WARNING: unrecognized options: $ac_unrecognized_opts" >&2;} fi - diff --git a/configure.ac b/configure.ac index 97048765..f06e2b13 100644 --- a/configure.ac +++ b/configure.ac @@ -184,4 +184,4 @@ AC_CONFIG_FILES([src/Makevars R/epiworldR-package.R]) # Use sed to replace the line useDynLib\([a-zA-Z]+ in the NAMESPACE # file with useDynLib(${epiworldname} -AC_OUTPUT \ No newline at end of file +AC_OUTPUT diff --git a/docker/Dockerfile b/docker/Dockerfile index e82145ec..d186bc94 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -16,4 +16,4 @@ RUN mkdir ~/.R && \ echo "CXX11FLAGS=-g -O0" >> ~/.R/Makevars && \ echo "SAFE_FLAGS=-g -O0" >> ~/.R/Makevars -CMD ["bash"] \ No newline at end of file +CMD ["bash"] diff --git a/docker/Makefile b/docker/Makefile index b161464e..3259f7f6 100644 --- a/docker/Makefile +++ b/docker/Makefile @@ -4,4 +4,4 @@ all: podman: podman build -t uofuepibio/epiworldr:debug -f Dockerfile . -.PHONY: all podman \ No newline at end of file +.PHONY: all podman diff --git a/man/figures/README-logit-model-1.png b/man/figures/README-logit-model-1.png index 3cd4a67d..d4e76bc4 100644 Binary files a/man/figures/README-logit-model-1.png and b/man/figures/README-logit-model-1.png differ diff --git a/man/figures/README-multiple-example-1.png b/man/figures/README-multiple-example-1.png index 819f2583..1cd860cd 100644 Binary files a/man/figures/README-multiple-example-1.png and b/man/figures/README-multiple-example-1.png differ diff --git a/man/figures/README-seir-conn-figures-1.png b/man/figures/README-seir-conn-figures-1.png index fe7d1ac4..e2ddcad3 100644 Binary files a/man/figures/README-seir-conn-figures-1.png and b/man/figures/README-seir-conn-figures-1.png differ diff --git a/man/figures/README-seir-conn-figures-2.png b/man/figures/README-seir-conn-figures-2.png index dcdf7d07..fecf53df 100644 Binary files a/man/figures/README-seir-conn-figures-2.png and b/man/figures/README-seir-conn-figures-2.png differ diff --git a/man/figures/README-seir-conn-figures-3.png b/man/figures/README-seir-conn-figures-3.png index 516d9201..8e525065 100644 Binary files a/man/figures/README-seir-conn-figures-3.png and b/man/figures/README-seir-conn-figures-3.png differ diff --git a/man/figures/README-sir-figures-1.png b/man/figures/README-sir-figures-1.png index 7203b10e..20a8417c 100644 Binary files a/man/figures/README-sir-figures-1.png and b/man/figures/README-sir-figures-1.png differ diff --git a/man/figures/README-sir-figures-2.png b/man/figures/README-sir-figures-2.png index 1963afc5..b6f289d9 100644 Binary files a/man/figures/README-sir-figures-2.png and b/man/figures/README-sir-figures-2.png differ diff --git a/man/figures/README-transmission-net-1.png b/man/figures/README-transmission-net-1.png index 6277139c..b3e48c13 100644 Binary files a/man/figures/README-transmission-net-1.png and b/man/figures/README-transmission-net-1.png differ diff --git a/paper.bib b/paper.bib index 28310667..413563e2 100644 --- a/paper.bib +++ b/paper.bib @@ -13,7 +13,7 @@ @Manual{abmR note = {R package version 1.0.9}, url = {https://CRAN.R-project.org/package=abmR}, } - + @article{gochanour2021abmr, title={abmR: An R package for agent-based model analysis of large-scale movements across taxa}, author={Gochanour, Benjamin and Fern{\'a}ndez-L{\'o}pez, Javier and Contina, Andrea}, @@ -30,7 +30,7 @@ @Manual{cystiSim note = {R package version 0.1.0}, url = {https://CRAN.R-project.org/package=cystiSim}, } - + @Manual{villager, title = {villager: A Framework for Designing and Running Agent Based Models}, author = {Thomas Thelen and Gerardo Aldana and Marcus Thomson and Toni Gonzalez}, @@ -38,7 +38,7 @@ @Manual{villager note = {R package version 1.1.1}, url = {https://CRAN.R-project.org/package=villager}, } - + @Article{RNetLogo, title = {{RNetLogo}: An {R} Package for Running and Exploring Individual-Based Models Implemented in {NetLogo}}, @@ -50,4 +50,4 @@ @Article{RNetLogo pages = {480--483}, url = {http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2011.00180.x/abstract}, -} \ No newline at end of file +} diff --git a/paper.md b/paper.md index 4ed5fc8e..3726b609 100644 --- a/paper.md +++ b/paper.md @@ -58,4 +58,4 @@ The `epiworldR` package addresses the need for sophisticated epidemiological res The development of the `epiworldR` package has ushered in a new era of agent-based modeling in the field of social science and epidemiology. By harnessing the power of C++ and the flexibility of R, this comprehensive package offers a multitude of features that enhance the modeling and analysis of complex infectious disease dynamics. The package's ability to handle multiple viruses and tools, its diverse set of epidemiological models, its capability to run simulations multiple times, and the inclusion of global actions capability empower researchers to explore a wide range of scenarios and make informed decisions regarding public health interventions. `epiworldR` serves as a valuable resource for the social science and epidemiological communities, enabling the study of real-world phenomena, prediction of outcomes, and policy analysis. As the field of epidemiology continues to evolve, `epiworldR` stands at the forefront, providing researchers and practitioners with a powerful tool to navigate the complexities of infectious diseases and contribute to the improvement of global health outcomes. -# References +# References diff --git a/playground/benchmark-seir.R b/playground/benchmark-seir.R index 0c860d18..78fd8f1b 100644 --- a/playground/benchmark-seir.R +++ b/playground/benchmark-seir.R @@ -10,29 +10,29 @@ for (n in ns) { sir <- epiworldR::ModelSEIR( name = "COVID-19", - prevalence = 0.01, - incubation_days = 7, + prevalence = 0.01, + incubation_days = 7, transmission_rate = 0.6, recovery_rate = 0.5 - ) |> + ) |> epiworldR::agents_smallworld(n = n, k = 20, p = 0.0, d = FALSE) |> epiworldR::add_virus( epiworldR::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 - ) |> + ) |> epiworldR::verbose_off() sirfast <- epiworldRdev::ModelSEIR( name = "COVID-19", - prevalence = 0.01, - incubation_days = 7, + prevalence = 0.01, + incubation_days = 7, transmission_rate = 0.6, recovery_rate = 0.5 - ) |> + ) |> epiworldRdev::agents_smallworld(n = n, k = 20, p = 0, d = FALSE) |> epiworldRdev::add_virus( epiworldRdev::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 - ) |> + ) |> epiworldRdev::verbose_off() @@ -47,4 +47,3 @@ for (n in ns) { } saveRDS(ans, "playground/benchmark-seir.rds") - diff --git a/playground/benchmark-seirconn.R b/playground/benchmark-seirconn.R index 6e6a947b..785e8b24 100644 --- a/playground/benchmark-seirconn.R +++ b/playground/benchmark-seirconn.R @@ -10,31 +10,31 @@ for (n in ns) { sir <- epiworldR::ModelSEIRCONN( name = "COVID-19", - prevalence = 0.01, + prevalence = 0.01, n = n, - contact_rate = 4, - incubation_days = 7, + contact_rate = 4, + incubation_days = 7, transmission_rate = 0.6, recovery_rate = 0.5 - ) |> + ) |> epiworldR::add_virus( epiworldR::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 - ) |> + ) |> epiworldR::verbose_off() sirfast <- epiworldRfaster::ModelSEIRCONN( name = "COVID-19", - prevalence = 0.01, + prevalence = 0.01, n = n, - contact_rate = 4, - incubation_days = 7, + contact_rate = 4, + incubation_days = 7, transmission_rate = 0.6, recovery_rate = 0.5 - ) |> + ) |> epiworldRfaster::add_virus( epiworldRfaster::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 - ) |> + ) |> epiworldRfaster::verbose_off() @@ -49,4 +49,3 @@ for (n in ns) { } saveRDS(ans, "playground/benchmark-seirconn.rds") - diff --git a/playground/epiworld_presentation.html b/playground/epiworld_presentation.html index 9f816de5..81b26249 100644 --- a/playground/epiworld_presentation.html +++ b/playground/epiworld_presentation.html @@ -33,11 +33,11 @@ .callout { margin-top: 1em; - margin-bottom: 1em; + margin-bottom: 1em; border-radius: .25rem; } - .callout.callout-style-simple { + .callout.callout-style-simple { padding: 0em 0.5em; border-left: solid #acacac .3rem; border-right: solid 1px silver; @@ -88,7 +88,7 @@ margin-top: 0.5em; margin-bottom: 0.5em; } - + .callout.callout-captioned.callout-style-simple .callout-content p { margin-top: 0; } @@ -139,7 +139,7 @@ .callout-caption { display: flex } - + .callout-icon::before { margin-top: 1rem; padding-right: .5rem; @@ -251,12 +251,12 @@ } .reveal .footnotes ol { counter-reset: ol; - list-style-type: none; + list-style-type: none; margin-left: 0; } .reveal .footnotes ol li:before { counter-increment: ol; - content: counter(ol) ". "; + content: counter(ol) ". "; } .reveal .footnotes ol li > p:first-child { display: inline-block; @@ -298,19 +298,19 @@ .reveal .slide > img.r-stretch.quarto-figure-center { display: block; margin-left: auto; - margin-right: auto; + margin-right: auto; } .reveal .slide > img.stretch.quarto-figure-left, .reveal .slide > img.r-stretch.quarto-figure-left { display: block; margin-left: 0; - margin-right: auto; + margin-right: auto; } .reveal .slide > img.stretch.quarto-figure-right, .reveal .slide > img.r-stretch.quarto-figure-right { display: block; margin-left: auto; - margin-right: 0; + margin-right: 0; } @@ -458,7 +458,7 @@

Visualization

- + @@ -480,7 +480,7 @@

Visualization

'autoAnimateUnmatched': true, 'menu': {"side":"left","useTextContentForMissingTitles":true,"markers":false,"loadIcons":false,"custom":[{"title":"Tools","icon":"","content":""}],"openButton":true}, 'smaller': false, - + // Display controls in the bottom right corner controls: false, @@ -658,7 +658,7 @@

Visualization

] }); - + - - \ No newline at end of file + + diff --git a/playground/epiworld_presentation.qmd b/playground/epiworld_presentation.qmd index e67c29be..cb3a76ad 100644 --- a/playground/epiworld_presentation.qmd +++ b/playground/epiworld_presentation.qmd @@ -81,49 +81,48 @@ Schools as entities: tbd library(epiworldR) seir <- ModelSEIR(name = "COVID-19", prevalence = 0.01, infectiousness = 0.9, incubation_days = 4, recovery = 0.3) -# Adding a Small world population +# Adding a Small world population agents_smallworld( seir, n = 200000, k = 5, d = FALSE, p = .01 - ) +) # Running and printing run(seir, ndays = 100, seed = 1912) seir - ``` ## Visualization ```{r} x <- get_hist_total(seir) -x$counts <- x$counts/1000 -x <- x[x$date < 50,] +x$counts <- x$counts / 1000 +x <- x[x$date < 50, ] with( - x[x$status == "Susceptible",], + x[x$status == "Susceptible", ], plot( x = date, y = counts, type = "l", col = "blue", ylim = range(x$counts), ylab = "Population (thousands)", xlab = "days", main = "SEIR model") - ) +) with( - x[x$status == "Exposed",], + x[x$status == "Exposed", ], lines(x = date, y = counts, col = "purple") - ) +) with( - x[x$status == "Infected",], + x[x$status == "Infected", ], lines(x = date, y = counts, col = "red") - ) +) with( - x[x$status == "Removed",], + x[x$status == "Removed", ], lines(x = date, y = counts, col = "darkgreen") - ) +) legend( "right", @@ -132,5 +131,5 @@ legend( lty = 1, lwd = 2, bty = "n" - ) +) ``` diff --git a/playground/transition_matrix_checker.Rmd b/playground/transition_matrix_checker.Rmd index 817920e0..d89c6a52 100644 --- a/playground/transition_matrix_checker.Rmd +++ b/playground/transition_matrix_checker.Rmd @@ -9,9 +9,9 @@ output: html_document ```{r eval=FALSE} library(epiworldR) -model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, -transmission_rate = 0.9, recovery_rate = 0.2, incubation_days = 4, -death_rate = 0.1) +model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, + transmission_rate = 0.9, recovery_rate = 0.2, incubation_days = 4, + death_rate = 0.1) # Adding a small world population agents_smallworld( @@ -20,8 +20,8 @@ agents_smallworld( k = 5, d = FALSE, p = .01 - ) - +) + # Running and printing run(model_seird, ndays = 100, seed = 1912) model_seird @@ -38,7 +38,7 @@ precovery <- 0.2 pnone <- (1 - pdie) * (1 - precovery) patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) -p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.08163265 +p_die_given_infected <- pdie * (1 - precovery) / (pnone + patmost_one) # 0.08163265 p_recovery_given_infected <- precovery * (1 - pdie) / (pnone + patmost_one) # 0.1836735 p_neither <- (1 - pdie) * (1 - precovery) / (pnone + patmost_one) # 0.7346939 @@ -51,15 +51,15 @@ $P(Dying | Infected) = \frac{p_{die} * (1-p_{recovery})}{(1-p_{die})*(1-p_{recov # An example with COVID-19 model_seirdconn <- ModelSEIRDCONN( name = "COVID-19", - prevalence = 0.01, + prevalence = 0.01, n = 100000, - contact_rate = 2, - incubation_days = 7, + contact_rate = 2, + incubation_days = 7, transmission_rate = 0.5, recovery_rate = 0.3, death_rate = 0.01 ) - + # Running and printing run(model_seirdconn, ndays = 100, seed = 1912) model_seirdconn @@ -76,11 +76,11 @@ get_transition_probability(model_seirdconn) pdie <- 0.01 precovery <- 0.3 -pnone <- (1-pdie)*(1-precovery) -patmost_one <- pdie*(1-precovery) + precovery*(1-pdie) +pnone <- (1 - pdie) * (1 - precovery) +patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) -p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.007021063 -p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.2978937 +p_die_given_infected <- pdie * (1 - precovery) / (pnone + patmost_one) # 0.007021063 +p_recovery_given_infected <- precovery * (1 - pdie) / (pnone + patmost_one) # 0.2978937 ``` # SIRDCONN Example @@ -94,7 +94,7 @@ p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.29 #' recovery_rate = 0.5, #' death_rate = 0.1 #' ) -#' +#' #' # Running and printing #' run(model_sirdconn, ndays = 100, seed = 1912) #' model_sirdconn @@ -107,21 +107,21 @@ p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.29 pdie <- 0.1 precovery <- 0.5 -pnone <- (1-pdie)*(1-precovery) -patmost_one <- pdie*(1-precovery) + precovery*(1-pdie) +pnone <- (1 - pdie) * (1 - precovery) +patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) -p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.05263158 -p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.4736842 +p_die_given_infected <- pdie * (1 - precovery) / (pnone + patmost_one) # 0.05263158 +p_recovery_given_infected <- precovery * (1 - pdie) / (pnone + patmost_one) # 0.4736842 ``` # SIRD Example ```{r eval = FALSE} model_sird <- ModelSIRD( - name = "COVID-19", - prevalence = 0.01, - transmission_rate = 0.9, - recovery_rate = 0.2, - death_rate = 0.01 + name = "COVID-19", + prevalence = 0.01, + transmission_rate = 0.9, + recovery_rate = 0.2, + death_rate = 0.01 ) # Adding a small world population @@ -131,8 +131,8 @@ agents_smallworld( k = 5, d = FALSE, p = .01 - ) - +) + # Running and printing run(model_sird, ndays = 100, seed = 1912) model_sird @@ -145,22 +145,22 @@ get_transition_probability(model_sird) pdie <- 0.01 precovery <- 0.2 -pnone <- (1-pdie)*(1-precovery) -patmost_one <- pdie*(1-precovery) + precovery*(1-pdie) +pnone <- (1 - pdie) * (1 - precovery) +patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) -p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.008016032 -p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.1983968 +p_die_given_infected <- pdie * (1 - precovery) / (pnone + patmost_one) # 0.008016032 +p_recovery_given_infected <- precovery * (1 - pdie) / (pnone + patmost_one) # 0.1983968 ``` # SISD Example ```{r eval = FALSE} model_sisd <- ModelSISD( - name = "COVID-19", - prevalence = 0.01, - transmission_rate = 0.9, - recovery_rate = 0.2, - death_rate = 0.01 - ) + name = "COVID-19", + prevalence = 0.01, + transmission_rate = 0.9, + recovery_rate = 0.2, + death_rate = 0.01 +) # Adding a small world population agents_smallworld( @@ -169,8 +169,8 @@ agents_smallworld( k = 5, d = FALSE, p = .01 - ) - +) + # Running and printing run(model_sisd, ndays = 100, seed = 1912) model_sisd @@ -183,9 +183,8 @@ get_transition_probability(model_sisd) pdie <- 0.01 precovery <- 0.2 -pnone <- (1-pdie)*(1-precovery) -patmost_one <- pdie*(1-precovery) + precovery*(1-pdie) +pnone <- (1 - pdie) * (1 - precovery) +patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) -p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.008016032 +p_die_given_infected <- pdie * (1 - precovery) / (pnone + patmost_one) # 0.008016032 ``` - diff --git a/pre-commit.yaml b/pre-commit.yaml new file mode 100644 index 00000000..2efb1421 --- /dev/null +++ b/pre-commit.yaml @@ -0,0 +1,23 @@ +name: Pre-commit + +on: + push: + branches: + - main + pull_request: + branches: + - main + +jobs: + pre-commit: + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v3 + + - name: Run pre-commit hooks + uses: pre-commit/action@v3.0.1 diff --git a/src/actions.cpp b/src/actions.cpp index 3868f090..2dbc7cf2 100644 --- a/src/actions.cpp +++ b/src/actions.cpp @@ -28,7 +28,7 @@ SEXP globalevent_tool_logit_cpp( cpp11::external_pointer> ptr( new GlobalEvent(action, name, day) ); - + return ptr; } @@ -45,11 +45,11 @@ SEXP globalevent_tool_cpp( *cpp11::external_pointer>(tool), prob )); - + cpp11::external_pointer> ptr( new GlobalEvent(action, name, day) ); - + return ptr; } @@ -66,11 +66,11 @@ SEXP globalevent_set_param_cpp( param, value )); - + cpp11::external_pointer> ptr( new GlobalEvent(action, name, day) ); - + return ptr; } @@ -79,13 +79,13 @@ SEXP globalevent_set_param_cpp( SEXP print_global_action_cpp( SEXP action ) { - + external_pointer> actionptr(action); - + actionptr->print(); - + return action; - + } @@ -94,14 +94,14 @@ SEXP add_globalevent_cpp( SEXP model, SEXP action ) { - + external_pointer> modelptr(model); external_pointer> actionptr(action); - + modelptr->add_globalevent(*actionptr); - + return model; - + } [[cpp11::register]] @@ -109,13 +109,13 @@ SEXP rm_globalevent_cpp( SEXP model, std::string name ) { - + external_pointer> modelptr(model); - + modelptr->rm_globalevent(name); - + return model; - + } [[cpp11::register]] @@ -124,24 +124,23 @@ SEXP globalevent_fun_cpp( std::string name, int day ) { - + GlobalFun fun_call = [fun](Model * model) -> void { - + cpp11::external_pointer> modelptr(model, false); sexp modelptrs(modelptr); modelptrs.attr("class") = "epiworld_model"; fun(modelptr); - + return; }; - + return external_pointer>( new GlobalEvent(fun_call, name, day) ); - - -} + +} diff --git a/src/agents.cpp b/src/agents.cpp index 04c6b8d1..cc22775a 100644 --- a/src/agents.cpp +++ b/src/agents.cpp @@ -12,17 +12,17 @@ using namespace cpp11; SEXP get_agents_cpp( SEXP model ) { - + // Making some room - + cpp11::external_pointer> ptr(model); cpp11::external_pointer >> agents( &ptr->get_agents(), false ); - + return agents; - + } [[cpp11::register]] @@ -30,12 +30,12 @@ SEXP get_agent_cpp( SEXP agents, size_t i ) { - + cpp11::external_pointer>> ptr(agents); - + if (i >= ptr->size()) stop("The agent index %lu is out of range.\n", i); - + return cpp11::external_pointer< Agent<> >( new Agent<>(ptr->operator[](i)) ); @@ -48,14 +48,14 @@ SEXP print_agent_cpp( SEXP model, bool compressed ) { - + cpp11::external_pointer> ptr(agent); cpp11::external_pointer> ptr_model(model); - + ptr->print(&(*ptr_model), compressed); - + return agent; - + } [[cpp11::register]] @@ -66,82 +66,80 @@ int get_state_agent_cpp(SEXP agent) { // Function to get agent's states using get_agents_states() [[cpp11::register]] std::vector get_agents_states_cpp(SEXP model) { - + cpp11::external_pointer> ptr(model); - + std::vector states; states.reserve(ptr->size()); - + auto states_uint = ptr->get_agents_states(); - + // Getting the model's states auto model_states = ptr->get_states(); // Copying the data to states for (auto i : states_uint) states.push_back(model_states[i]); - + return states; - + } [[cpp11::register]] SEXP add_virus_agent_cpp(SEXP agent, SEXP model, SEXP virus, int state_new, int queue) { - + cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_model(model); cpp11::external_pointer> ptr_virus(virus); - + ptr_agent->set_virus(*ptr_virus, &(*ptr_model)); - + return agent; - + } [[cpp11::register]] SEXP add_tool_agent_cpp(SEXP agent, SEXP model, SEXP tool, int state_new, int queue) { - + cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_model(model); cpp11::external_pointer> ptr_tool(tool); - + ptr_agent->add_tool(*ptr_tool, &(*ptr_model)); - + return agent; - + } [[cpp11::register]] bool has_virus_cpp(SEXP agent, SEXP virus) { - + cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_virus(virus); - + return ptr_agent->has_virus(*ptr_virus); - + } [[cpp11::register]] bool has_tool_cpp(SEXP agent, SEXP tool) { - + cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_tool(tool); - + return ptr_agent->has_tool(*ptr_tool); - + } [[cpp11::register]] SEXP change_state_cpp(SEXP agent, SEXP model, int new_state, int queue) { - + cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_model(model); - + ptr_agent->change_state(&(*ptr_model), new_state, queue); - - return agent; - -} + return agent; +} diff --git a/src/db.cpp b/src/db.cpp index 3731caa7..1f7aa397 100644 --- a/src/db.cpp +++ b/src/db.cpp @@ -13,23 +13,23 @@ using namespace cpp11; cpp11::data_frame get_hist_total_cpp( SEXP model ) { - + // Making some room std::vector< int > date; std::vector< std::string > state; std::vector< int > counts; - + cpp11::external_pointer> ptr(model); ptr->get_db().get_hist_total(&date, &state, &counts); - + // Preparing the output cpp11::writable::data_frame res({ "date"_nm = date, "state"_nm = state, "counts"_nm = counts }); - - + + return res; } @@ -37,23 +37,23 @@ cpp11::data_frame get_hist_total_cpp( cpp11::data_frame get_hist_virus_cpp( SEXP model ) { - + cpp11::external_pointer> ptr(model); - + std::vector date; std::vector id; std::vector state; std::vector counts; - + ptr->get_db().get_hist_virus( date, id, state, counts ); - + // Mapping the id to the name std::vector< std::string > viruses; for (auto i : ptr->get_viruses()) viruses.push_back(i->get_name()); - + // Mapping using std::transform std::vector< std::string > vnames(id.size()); std::transform( @@ -62,27 +62,27 @@ cpp11::data_frame get_hist_virus_cpp( ); return cpp11::writable::data_frame({ - "date"_nm = date, + "date"_nm = date, "virus_id"_nm = id, "virus"_nm = vnames, "state"_nm = state, "counts"_nm = counts, }); - + } [[cpp11::register]] cpp11::data_frame get_hist_tool_cpp( SEXP model ) { - + cpp11::external_pointer> ptr(model); - + std::vector date; std::vector id; std::vector state; std::vector counts; - + ptr->get_db().get_hist_tool( date, id, state, counts ); @@ -97,25 +97,25 @@ cpp11::data_frame get_hist_tool_cpp( id.begin(), id.end(), tnames.begin(), [&tools](int i) { return tools[i]; } ); - + return cpp11::writable::data_frame({ - "date"_nm = date, + "date"_nm = date, "tool_id"_nm = id, "tool"_nm = tnames, "state"_nm = state, "counts"_nm = counts, }); - + } [[cpp11::register]] doubles get_transition_probability_cpp( SEXP model ) { - + cpp11::external_pointer> ptr(model); return cpp11::writable::doubles(ptr->get_db().transition_probability(false)); - + } [[cpp11::register]] @@ -123,45 +123,45 @@ cpp11::data_frame get_hist_transition_matrix_cpp( SEXP model, bool skip_zeros ) { - + cpp11::external_pointer> ptr(model); - + std::vector< std::string > state_from; std::vector< std::string > state_to; std::vector< int > date; std::vector< int > counts; - + ptr->get_db().get_hist_transition_matrix( state_from, state_to, date, counts, skip_zeros ); - + return cpp11::writable::data_frame({ - "state_from"_nm = state_from, + "state_from"_nm = state_from, "state_to"_nm = state_to, "date"_nm = date, "counts"_nm = counts, }); - + } - + [[cpp11::register]] cpp11::data_frame get_reproductive_number_cpp( SEXP model ) { - + // Making some room std::vector< int > virus; std::vector< int > source; std::vector< int > source_exposure_date; std::vector< int > counts; - + // Getting the right class cpp11::external_pointer> ptr(model); std::unordered_map< std::vector< int >, int, epiworld::vecHasher> rn = ptr->get_db().reproductive_number(); - - for (const auto & m : rn) + + for (const auto & m : rn) { virus.push_back(m.first[0u]); source.push_back(m.first[1u]); @@ -179,7 +179,7 @@ cpp11::data_frame get_reproductive_number_cpp( virus.begin(), virus.end(), vnames.begin(), [&viruses](int i) { return viruses[i]; } ); - + return cpp11::writable::data_frame({ "virus_id"_nm = virus, "virus"_nm = vnames, @@ -187,22 +187,22 @@ cpp11::data_frame get_reproductive_number_cpp( "source_exposure_date"_nm = source_exposure_date, "rt"_nm = counts }); - + } [[cpp11::register]] cpp11::data_frame get_transmissions_cpp( SEXP model ) { - + cpp11::external_pointer> ptr(model); - + std::vector date; std::vector source; std::vector target; std::vector virus; std::vector source_exposure_date; - + ptr->get_db().get_transmissions( date, source, @@ -221,7 +221,7 @@ cpp11::data_frame get_transmissions_cpp( virus.begin(), virus.end(), vnames.begin(), [&viruses](int i) { return viruses[i]; } ); - + return cpp11::writable::data_frame({ "date"_nm = date, "source"_nm = source, @@ -230,21 +230,21 @@ cpp11::data_frame get_transmissions_cpp( "virus"_nm = vnames, "source_exposure_date"_nm = source_exposure_date, }); - + } [[cpp11::register]] cpp11::data_frame get_generation_time_cpp( SEXP model ) { - + cpp11::external_pointer> ptr(model); - + std::vector agent_id; std::vector virus_id; std::vector date; std::vector gentime; - + ptr->get_db().generation_time( agent_id, virus_id, @@ -262,7 +262,7 @@ cpp11::data_frame get_generation_time_cpp( virus_id.begin(), virus_id.end(), vnames.begin(), [&viruses](int i) { return viruses[i]; } ); - + return cpp11::writable::data_frame({ "agent"_nm = agent_id, "virus_id"_nm = virus_id, @@ -270,7 +270,7 @@ cpp11::data_frame get_generation_time_cpp( "date"_nm = date, "gentime"_nm = gentime }); - + } [[cpp11::register]] @@ -282,7 +282,7 @@ cpp11::writable::doubles get_today_total_cpp(SEXP model) { std::vector< std::string > names; ptr->get_db().get_today_total(&names, &totals); - cpp11::writable::doubles totals_r(totals.begin(), totals.end()); + cpp11::writable::doubles totals_r(totals.begin(), totals.end()); cpp11::writable::strings names_r(names.begin(), names.end()); totals_r.names() = names_r; diff --git a/src/entities.cpp b/src/entities.cpp index f20534c6..33d1622e 100644 --- a/src/entities.cpp +++ b/src/entities.cpp @@ -12,9 +12,9 @@ using namespace epiworld; SEXP get_entities_cpp( SEXP model ) { - + // Making some room - + cpp11::external_pointer> ptr(model); cpp11::writable::list res; @@ -23,9 +23,9 @@ SEXP get_entities_cpp( cpp11::external_pointer> entity(&i, false); res.push_back(entity); } - + return res; - + } [[cpp11::register]] @@ -33,15 +33,15 @@ SEXP get_entity_cpp( SEXP entities, int idx ) { - + cpp11::external_pointer> > ptr(entities); - + cpp11::external_pointer> entity( &ptr->at(static_cast(idx)), false ); - + return entity; - + } @@ -52,7 +52,7 @@ SEXP entity_cpp( bool as_proportion, bool to_unassigned ) { - + cpp11::external_pointer> ptr( new Entity<>( name, @@ -63,9 +63,9 @@ SEXP entity_cpp( ) ) ); - + return ptr; - + } [[cpp11::register]] @@ -76,13 +76,13 @@ int get_entity_size_cpp(SEXP entity) { [[cpp11::register]] int entity_add_agent_cpp(SEXP entity, SEXP agent, SEXP model) { - + cpp11::external_pointer> ptr(entity); cpp11::external_pointer> ptr_agent(agent); cpp11::external_pointer> ptr_model(model); - + ptr->add_agent(&(*ptr_agent), &(*ptr_model)); - + return 0; } @@ -96,24 +96,24 @@ int add_entity_cpp( SEXP model, SEXP entity ) { - + cpp11::external_pointer> ptr_model(model); cpp11::external_pointer> ptr_entity(entity); - + ptr_model->add_entity(*ptr_entity); - + return 0; } [[cpp11::register]] int rm_entity_cpp(SEXP model, int entity_pos) { - + cpp11::external_pointer> ptr_model(model); - + ptr_model->rm_entity( static_cast(entity_pos) ); - + return 0; } @@ -142,12 +142,12 @@ int load_agents_entities_ties_cpp( [[cpp11::register]] cpp11::data_frame entity_get_agents_cpp(SEXP entity) { - + cpp11::external_pointer> ptr(entity); - + cpp11::writable::integers agent; cpp11::writable::integers entity_id; - + int id = static_cast(ptr->get_id()); auto agents_ids = ptr->get_agents(); size_t nagents = ptr->size(); @@ -160,7 +160,7 @@ cpp11::data_frame entity_get_agents_cpp(SEXP entity) { "agent"_nm = agent, "entity"_nm = entity_id }); - + } [[cpp11::register]] @@ -234,4 +234,3 @@ SEXP distribute_entity_to_set_cpp( // cpp11::external_pointer>(entity)->set_name(name); // return 0; // } - diff --git a/src/epimodels.cpp b/src/epimodels.cpp index bdf41cd8..388acfdd 100644 --- a/src/epimodels.cpp +++ b/src/epimodels.cpp @@ -18,18 +18,18 @@ SEXP ModelSURV_cpp( std::string name, double prevalence, double efficacy_vax, - double latent_period, - double prob_symptoms, - double prop_vaccinated, - double prop_vax_redux_transm, - double infect_period, - double prop_vax_redux_infect, - double surveillance_prob, - double transmission_rate, - double prob_death, - double prob_noreinfect + double latent_period, + double prob_symptoms, + double prop_vaccinated, + double prop_vax_redux_transm, + double infect_period, + double prop_vax_redux_infect, + double surveillance_prob, + double transmission_rate, + double prob_death, + double prob_noreinfect ) { - + // Creating a pointer to a ModelSIR model WrapSURV(ptr)( @@ -37,16 +37,16 @@ SEXP ModelSURV_cpp( name, prevalence, efficacy_vax, - latent_period, - infect_period, - prob_symptoms, - prop_vaccinated, - prop_vax_redux_transm, - prop_vax_redux_infect, - surveillance_prob, - transmission_rate, - prob_death, - prob_noreinfect + latent_period, + infect_period, + prob_symptoms, + prop_vaccinated, + prop_vax_redux_transm, + prop_vax_redux_infect, + surveillance_prob, + transmission_rate, + prob_death, + prob_noreinfect ) ); @@ -65,9 +65,9 @@ SEXP ModelSURV_cpp( double transmission_rate, double incubation_days, double recovery_rate - + ) { - + // Creating a pointer to a ModelSIR model WrapSEIR(ptr)( new epiworld::epimodels::ModelSEIR<>( @@ -78,12 +78,12 @@ SEXP ModelSURV_cpp( recovery_rate ) ); - + return ptr; } - + #undef WrapSEIR - + #define WrapSIS(a) \ cpp11::external_pointer> (a) @@ -94,7 +94,7 @@ SEXP ModelSIS_cpp( double transmission_rate, double recovery_rate ) { - + // Creating a pointer to a ModelSIR model WrapSIS(ptr)( new epiworld::epimodels::ModelSIS<>( @@ -104,15 +104,15 @@ SEXP ModelSIS_cpp( recovery_rate ) ); - - + + return ptr; } - - + + #undef WrapSIS - - + + #define WrapSIRCONN(a) \ cpp11::external_pointer> (a) @@ -122,10 +122,10 @@ SEXP ModelSIRCONN_cpp( unsigned int n, double prevalence, double contact_rate, - double transmission_rate, + double transmission_rate, double recovery_rate ) { - + // Creating a pointer to a ModelSIR model WrapSIRCONN(ptr)( new epiworld::epimodels::ModelSIRCONN<>( @@ -137,13 +137,13 @@ SEXP ModelSIRCONN_cpp( recovery_rate ) ); - + return ptr; } - + #undef WrapSIRCONN - + #define WrapSIR(a) \ cpp11::external_pointer> (a) @@ -154,7 +154,7 @@ SEXP ModelSIR_cpp( double transmission_rate, double recovery_rate ) { - + // Creating a pointer to a ModelSIR model WrapSIR(ptr)( new epiworld::epimodels::ModelSIR<>( @@ -164,10 +164,10 @@ SEXP ModelSIR_cpp( recovery_rate ) ); - + return ptr; } - + #undef WrapSIR #define WrapSIRD(a) \ @@ -181,7 +181,7 @@ SEXP ModelSIRD_cpp( double recovery_rate, double death_rate ) { - + // Creating a pointer to a ModelSIRD model WrapSIRD(ptr)( new epiworld::epimodels::ModelSIRD<>( @@ -192,12 +192,12 @@ SEXP ModelSIRD_cpp( death_rate ) ); - + return ptr; } - + #undef WrapSIRD - + #define WrapSEIRD(a) \ cpp11::external_pointer> (a) @@ -209,9 +209,9 @@ SEXP ModelSIRD_cpp( double incubation_days, double recovery_rate, double death_rate - + ) { - + // Creating a pointer to a ModelSEIRD model WrapSEIRD(ptr)( new epiworld::epimodels::ModelSEIRD<>( @@ -223,12 +223,12 @@ SEXP ModelSIRD_cpp( death_rate ) ); - + return ptr; } - + #undef WrapSEIRD - + #define WrapSISD(a) \ cpp11::external_pointer> (a) @@ -240,7 +240,7 @@ SEXP ModelSISD_cpp( double recovery_rate, double death_rate ) { - + // Creating a pointer to a ModelSISD model WrapSISD(ptr)( new epiworld::epimodels::ModelSISD<>( @@ -251,10 +251,10 @@ SEXP ModelSISD_cpp( death_rate ) ); - + return ptr; } - + #undef WrapSISD #define WrapSIRDCONN(a) \ @@ -266,11 +266,11 @@ SEXP ModelSIRDCONN_cpp( unsigned int n, double prevalence, double contact_rate, - double transmission_rate, - double recovery_rate, + double transmission_rate, + double recovery_rate, double death_rate ) { - + // Creating a pointer to a ModelSIR model WrapSIRDCONN(ptr)( new epiworld::epimodels::ModelSIRDCONN<>( @@ -283,12 +283,12 @@ SEXP ModelSIRDCONN_cpp( death_rate ) ); - + return ptr; } - + #undef WrapSIRDCONN - + #define WrapSEIRDCONN(a) \ cpp11::external_pointer> (a) @@ -303,7 +303,7 @@ SEXP ModelSEIRDCONN_cpp( double recovery_rate, double death_rate ) { - + // Creating a pointer to a ModelSIR model WrapSEIRDCONN(ptr)( new epiworld::epimodels::ModelSEIRDCONN<>( @@ -317,13 +317,13 @@ SEXP ModelSEIRDCONN_cpp( death_rate ) ); - + return ptr; } - - + + #undef WrapSEIRDCONN - + #define WrapSEIRCONN(a) \ cpp11::external_pointer> (a) @@ -337,7 +337,7 @@ SEXP ModelSEIRCONN_cpp( double incubation_days, double recovery_rate ) { - + // Creating a pointer to a ModelSIR model WrapSEIRCONN(ptr)( new epiworld::epimodels::ModelSEIRCONN<>( @@ -350,13 +350,13 @@ SEXP ModelSEIRCONN_cpp( recovery_rate ) ); - + return ptr; } - - + + #undef WrapSEIRCONN - + [[cpp11::register]] SEXP ModelSIRLogit_cpp( std::string vname, @@ -370,16 +370,16 @@ SEXP ModelSIRLogit_cpp( double recovery_rate, double prevalence ) { - + std::vector< size_t > cinfect; std::vector< size_t > crecover; - + for (auto i : coef_infect_cols) cinfect.push_back(static_cast(i)); - + for (auto i : coef_recover_cols) crecover.push_back(static_cast(i)); - + cpp11::external_pointer> ptr( new epiworld::epimodels::ModelSIRLogit<>( vname, @@ -394,9 +394,9 @@ SEXP ModelSIRLogit_cpp( prevalence ) ); - + return ptr; - + } @@ -416,7 +416,7 @@ SEXP ModelDiffNet_cpp( std::vector< size_t > data_cols_s; for (auto i : data_cols) data_cols_s.push_back(static_cast(i)); - + cpp11::external_pointer> ptr( new epiworld::epimodels::ModelDiffNet<>( name, @@ -429,9 +429,9 @@ SEXP ModelDiffNet_cpp( params ) ); - + return ptr; - + } [[cpp11::register]] @@ -444,7 +444,7 @@ SEXP ModelSIRMixing_cpp( double recovery_rate, std::vector< double > contact_matrix ) { - + // Creating a pointer to a ModelSIRMixing model cpp11::external_pointer> ptr( new epiworld::epimodels::ModelSIRMixing<>( @@ -457,7 +457,7 @@ SEXP ModelSIRMixing_cpp( contact_matrix ) ); - + return ptr; } @@ -473,7 +473,7 @@ SEXP ModelSEIRMixing_cpp( double recovery_rate, std::vector< double > contact_matrix ) { - + // Creating a pointer to a ModelSIRMixing model cpp11::external_pointer> ptr( new epiworld::epimodels::ModelSEIRMixing<>( @@ -487,7 +487,7 @@ SEXP ModelSEIRMixing_cpp( contact_matrix ) ); - + return ptr; - -} \ No newline at end of file + +} diff --git a/src/epiworld-common.h b/src/epiworld-common.h index 6c5892ac..bd6317fb 100644 --- a/src/epiworld-common.h +++ b/src/epiworld-common.h @@ -3,4 +3,4 @@ #include "epiworld/epiworld.hpp" -#endif \ No newline at end of file +#endif diff --git a/src/model.cpp b/src/model.cpp index 3a054d04..44d0c8a5 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -7,12 +7,12 @@ using namespace cpp11; [[cpp11::register]] SEXP print_cpp(SEXP m, bool lite) { - + external_pointer> ptr(m); ptr->print(lite); - + return m; - + } [[cpp11::register]] @@ -24,12 +24,12 @@ SEXP agents_smallworld_cpp( double p = .01 ) { - + external_pointer> ptr(m); ptr->agents_smallworld(n, k, d, p); - + return m; - + } [[cpp11::register]] @@ -40,22 +40,22 @@ SEXP agents_from_edgelist_cpp( int size, bool directed ) { - + external_pointer> ptr(m); - ptr->agents_from_edgelist(source, target, size, directed); - + ptr->agents_from_edgelist(source, target, size, directed); + return m; - + } [[cpp11::register]] SEXP run_cpp(SEXP m, int ndays, int seed) { - + external_pointer> ptr(m); ptr->run(ndays, seed); - + return m; - + } typedef std::function*)> funptr; @@ -73,7 +73,7 @@ SEXP make_saver_cpp( bool reproductive, bool generation ) { - + funptr* saver = new funptr(make_save_run( fn, total_hist, @@ -83,14 +83,14 @@ SEXP make_saver_cpp( tool_hist, transmission, transition, - reproductive, + reproductive, generation )); - + external_pointer sav_ptr(saver); - + return sav_ptr; - + } [[cpp11::register]] @@ -104,10 +104,10 @@ SEXP run_multiple_cpp( bool verbose, int nthreads ) { - + external_pointer> ptr(m); external_pointer sav_ptr(saver); - + ptr->run_multiple( static_cast< epiworld_fast_uint >(ndays), static_cast< epiworld_fast_uint >(nsims), @@ -117,31 +117,31 @@ SEXP run_multiple_cpp( verbose, nthreads ); - + return m; - + } [[cpp11::register]] SEXP queuing_on_cpp( SEXP model ) { - + external_pointer> ptr(model); ptr->queuing_on(); return model; - + } [[cpp11::register]] SEXP queuing_off_cpp( SEXP model ) { - + external_pointer> ptr(model); ptr->queuing_off(); return model; - + } [[cpp11::register]] @@ -152,10 +152,10 @@ double get_param_cpp(SEXP model, std::string pname) { [[cpp11::register]] SEXP set_param_cpp(SEXP model, std::string pname, double val) { - + external_pointer> ptr(model); ptr->operator()(pname) = val; - + return model; } @@ -176,60 +176,60 @@ std::string get_name_cpp(SEXP model) { strings get_states_cpp( SEXP model ) { - + external_pointer> ptr(model); return writable::strings(ptr->get_states()); - + } [[cpp11::register]] SEXP verbose_on_cpp(SEXP model) { - + external_pointer> ptr(model); - ptr->verbose_on(); + ptr->verbose_on(); return model; - + } [[cpp11::register]] SEXP verbose_off_cpp(SEXP model) { - + external_pointer> ptr(model); - ptr->verbose_off(); + ptr->verbose_off(); return model; - + } [[cpp11::register]] int get_n_viruses_cpp(SEXP model) { - + external_pointer> ptr(model); return static_cast(ptr->get_n_viruses()); - + } [[cpp11::register]] int get_n_tools_cpp(SEXP model) { - + external_pointer> ptr(model); return static_cast(ptr->get_n_tools()); - + } [[cpp11::register]] int get_ndays_cpp(SEXP model) { - + external_pointer> ptr(model); return static_cast(ptr->get_ndays()); - + } [[cpp11::register]] int get_n_replicates_cpp(SEXP model) { - + external_pointer> ptr(model); return static_cast(ptr->get_n_replicates()); - + } [[cpp11::register]] @@ -245,52 +245,52 @@ SEXP set_agents_data_cpp(SEXP model, SEXP data, int ncols) { REAL(data), ncols ); - + return model; - + } [[cpp11::register]] int get_agents_data_ncols_cpp(SEXP model) { - + return static_cast( external_pointer>(model)->get_agents_data_ncols() ); - + } [[cpp11::register]] SEXP get_virus_model_cpp(SEXP model, int virus_pos) { external_pointer> modelptr(model); - + external_pointer> res( &modelptr->get_virus(static_cast(virus_pos)), false ); - + return res; - + } [[cpp11::register]] SEXP get_tool_model_cpp(SEXP model, int tool_pos) { external_pointer> modelptr(model); - + external_pointer> res( &modelptr->get_tool(static_cast(tool_pos)), false ); - + return res; - + } [[cpp11::register]] cpp11::data_frame get_network_cpp(SEXP model) { - + external_pointer> modelptr(model); - + std::vector from; std::vector to; @@ -300,7 +300,7 @@ cpp11::data_frame get_network_cpp(SEXP model) { "from"_nm = from, "to"_nm = to }); - + } [[cpp11::register]] @@ -311,20 +311,19 @@ SEXP initial_states_cpp(SEXP model, cpp11::doubles proportions) { std::vector< double > states_vec(proportions.begin(), proportions.end()); modelptr->initial_states(states_vec, std::vector< int >({})); - + return model; - + } // Function for cloning a model [[cpp11::register]] SEXP clone_model_cpp(const SEXP & model) { - + external_pointer> modelptr(model); - + return external_pointer>( new Model<>(*modelptr) ); - -} +} diff --git a/src/tool.cpp b/src/tool.cpp index d39f8ab6..a7dd9efd 100644 --- a/src/tool.cpp +++ b/src/tool.cpp @@ -22,40 +22,40 @@ SEXP tool_cpp( double recovery_enhancer, double death_reduction ) { - + WrapTool(tool)(new epiworld::Tool( name, prevalence, as_proportion )); - + if (susceptibility_reduction > 0) tool->set_susceptibility_reduction(susceptibility_reduction); - + if (transmission_reduction > 0) tool->set_transmission_reduction(transmission_reduction); - + if (recovery_enhancer > 0) tool->set_recovery_enhancer(recovery_enhancer); - + if (death_reduction > 0) tool->set_death_reduction(death_reduction); - + return tool; - + } - + [[cpp11::register]] int add_tool_cpp(SEXP m, SEXP t) { - + cpp11::external_pointer>(m)->add_tool( *cpp11::external_pointer>(t) ); - + return 0; } - + [[cpp11::register]] SEXP rm_tool_cpp(SEXP m, size_t tool_pos) { cpp11::external_pointer>(m)->rm_tool(tool_pos); @@ -69,9 +69,9 @@ SEXP tool_fun_logit_cpp( doubles coefs, SEXP model ) { - - external_pointer> mptr(model); - + + external_pointer> mptr(model); + external_pointer> res( new ToolFun<>( tool_fun_logit( @@ -81,157 +81,157 @@ SEXP tool_fun_logit_cpp( ) ) ); - + return res; - + } // Probability of transmission ------------------------------------------------- [[cpp11::register]] SEXP set_transmission_reduction_cpp(SEXP tool, double prob) { - + WrapTool(toolptr)(tool); toolptr->set_transmission_reduction(prob); return tool; - + } [[cpp11::register]] SEXP set_transmission_reduction_ptr_cpp(SEXP tool, SEXP model, std::string param) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); - + toolptr->set_transmission_reduction( &(mptr->operator()(param)) ); - + return tool; - + } [[cpp11::register]] SEXP set_transmission_reduction_fun_cpp(SEXP tool, SEXP model, SEXP tfun) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); external_pointer> tfunptr(tfun); - + toolptr->set_transmission_reduction_fun(*tfunptr); - + return tool; - + } // Probability of recovery ----------------------------------------------------- [[cpp11::register]] SEXP set_susceptibility_reduction_cpp(SEXP tool, double prob) { - + WrapTool(toolptr)(tool); toolptr->set_susceptibility_reduction(prob); return tool; - + } [[cpp11::register]] SEXP set_susceptibility_reduction_ptr_cpp(SEXP tool, SEXP model, std::string param) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); - + toolptr->set_susceptibility_reduction( &(mptr->operator()(param)) ); - + return tool; - + } [[cpp11::register]] SEXP set_susceptibility_reduction_fun_cpp(SEXP tool, SEXP model, SEXP tfun) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); external_pointer> tfunptr(tfun); - + toolptr->set_susceptibility_reduction_fun(*tfunptr); - + return tool; - + } // Recovery enhancer ----------------------------------------------------------- [[cpp11::register]] SEXP set_recovery_enhancer_cpp(SEXP tool, double prob) { - + WrapTool(toolptr)(tool); toolptr->set_recovery_enhancer(prob); return tool; - + } [[cpp11::register]] SEXP set_recovery_enhancer_ptr_cpp(SEXP tool, SEXP model, std::string param) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); - + toolptr->set_recovery_enhancer( &(mptr->operator()(param)) ); - + return tool; - + } [[cpp11::register]] SEXP set_recovery_enhancer_fun_cpp(SEXP tool, SEXP model, SEXP tfun) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); external_pointer> tfunptr(tfun); - + toolptr->set_recovery_enhancer_fun(*tfunptr); - + return tool; - + } // Death reduction ------------------------------------------------------------- [[cpp11::register]] SEXP set_death_reduction_cpp(SEXP tool, double prob) { - + WrapTool(toolptr)(tool); toolptr->set_death_reduction(prob); return tool; - + } [[cpp11::register]] SEXP set_death_reduction_ptr_cpp(SEXP tool, SEXP model, std::string param) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); - + toolptr->set_death_reduction( &(mptr->operator()(param)) ); - + return tool; - + } [[cpp11::register]] SEXP set_death_reduction_fun_cpp(SEXP tool, SEXP model, SEXP tfun) { - + WrapTool(toolptr)(tool); external_pointer> mptr(model); external_pointer> tfunptr(tfun); - + toolptr->set_death_reduction_fun(*tfunptr); - + return tool; - + } @@ -248,30 +248,30 @@ SEXP set_name_tool_cpp(SEXP tool, std::string name) { [[cpp11::register]] SEXP print_tool_cpp(SEXP t) { - + WrapTool(tptr)(t); tptr->print(); return t; - + } // Function to get agent's viruses using get_agents_viruses() [[cpp11::register]] cpp11::writable::list get_agents_tools_cpp(SEXP model) { - + cpp11::external_pointer> ptr(model); - + cpp11::writable::list tools; - + for (auto & agent : ptr->get_agents()) tools.push_back( cpp11::external_pointer< Tools<> >( new Tools<>(agent.get_tools()) ) ); - + return tools; - + } [[cpp11::register]] @@ -286,14 +286,14 @@ SEXP set_distribution_tool_cpp( SEXP tool, SEXP distfun ) { - + WrapTool(toolptr)(tool); external_pointer> tfunptr(distfun); - + toolptr->set_distribution(*tfunptr); - + return tool; - + } [[cpp11::register]] @@ -310,9 +310,9 @@ SEXP distribute_tool_randomly_cpp( ) ) ); - + return res; - + } [[cpp11::register]] @@ -335,9 +335,9 @@ SEXP distribute_tool_to_set_cpp( distribute_tool_to_set(ids) ) ); - + return res; - + } diff --git a/src/virus.cpp b/src/virus.cpp index eed27cb1..8dd8950e 100644 --- a/src/virus.cpp +++ b/src/virus.cpp @@ -23,22 +23,22 @@ SEXP virus_cpp( double post_immunity, double incubation ) { - + WrapVirus(virus)(new epiworld::Virus( name, prevalence, as_proportion )); - + virus->set_prob_infecting(prob_infecting); virus->set_prob_recovery(prob_recovery); virus->set_prob_death(prob_death); - + if (post_immunity > 0.0) virus->set_post_immunity(post_immunity); - + virus->set_incubation(incubation); return virus; - + } [[cpp11::register]] @@ -48,41 +48,41 @@ SEXP virus_set_state_cpp( size_t end, size_t removed ) { - + WrapVirus(vptr)(v); vptr->set_state( init, end, removed ); - + return v; - + } [[cpp11::register]] SEXP add_virus_cpp(SEXP m, SEXP v) { - + external_pointer>(m)->add_virus( *external_pointer>(v) ); - + return m; } [[cpp11::register]] SEXP rm_virus_cpp(SEXP m, size_t virus_pos) { - + external_pointer>(m)->rm_virus(virus_pos); return m; - + } [[cpp11::register]] SEXP print_virus_cpp(SEXP v) { - + WrapVirus(vptr)(v); vptr->print(); return v; - + } // Virus function -------------------------------------------------------------- @@ -92,9 +92,9 @@ SEXP virus_fun_logit_cpp( doubles coefs, SEXP model ) { - - external_pointer> mptr(model); - + + external_pointer> mptr(model); + external_pointer> res( new VirusFun<>( virus_fun_logit( @@ -104,157 +104,157 @@ SEXP virus_fun_logit_cpp( ) ) ); - + return res; - + } - + // Probability of infection ---------------------------------------------------- [[cpp11::register]] SEXP set_prob_infecting_cpp(SEXP virus, double prob) { - + WrapVirus(vptr)(virus); vptr->set_prob_infecting(prob); return virus; - + } - + [[cpp11::register]] SEXP set_prob_infecting_ptr_cpp(SEXP virus, SEXP model, std::string param) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); - + vptr->set_prob_infecting( &(mptr->operator()(param)) ); - + return virus; - + } - + [[cpp11::register]] SEXP set_prob_infecting_fun_cpp(SEXP virus, SEXP model, SEXP vfun) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); external_pointer> vfunptr(vfun); - + vptr->set_prob_infecting_fun(*vfunptr); - + return virus; - + } - + // Probability of recovery ----------------------------------------------------- [[cpp11::register]] SEXP set_prob_recovery_cpp(SEXP virus, double prob) { - + WrapVirus(vptr)(virus); vptr->set_prob_recovery(prob); return virus; - + } [[cpp11::register]] SEXP set_prob_recovery_ptr_cpp(SEXP virus, SEXP model, std::string param) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); - + vptr->set_prob_recovery( &(mptr->operator()(param)) ); - + return virus; - + } [[cpp11::register]] SEXP set_prob_recovery_fun_cpp(SEXP virus, SEXP model, SEXP vfun) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); external_pointer> vfunptr(vfun); - + vptr->set_prob_recovery_fun(*vfunptr); - + return virus; - + } - + // Probability of death -------------------------------------------------------- [[cpp11::register]] SEXP set_prob_death_cpp(SEXP virus, double prob) { - + WrapVirus(vptr)(virus); vptr->set_prob_death(prob); return virus; - + } [[cpp11::register]] SEXP set_prob_death_ptr_cpp(SEXP virus, SEXP model, std::string param) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); - + vptr->set_prob_death( &(mptr->operator()(param)) ); - + return virus; - + } [[cpp11::register]] SEXP set_prob_death_fun_cpp(SEXP virus, SEXP model, SEXP vfun) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); external_pointer> vfunptr(vfun); - + vptr->set_prob_death_fun(*vfunptr); return virus; - + } // Incubation period ---------------------------------------------------------- [[cpp11::register]] SEXP set_incubation_cpp(SEXP virus, double prob) { - + WrapVirus(vptr)(virus); vptr->set_incubation(prob); return virus; - + } [[cpp11::register]] SEXP set_incubation_ptr_cpp(SEXP virus, SEXP model, std::string param) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); - + vptr->set_incubation( &(mptr->operator()(param)) ); - + return virus; - + } [[cpp11::register]] SEXP set_incubation_fun_cpp(SEXP virus, SEXP model, SEXP vfun) { - + WrapVirus(vptr)(virus); external_pointer> mptr(model); external_pointer> vfunptr(vfun); - + vptr->set_incubation_fun(*vfunptr); - + return virus; - + } [[cpp11::register]] @@ -324,4 +324,4 @@ SEXP distribute_virus_to_set_cpp( } -#undef WrapVirus \ No newline at end of file +#undef WrapVirus diff --git a/tests/tinytest.R b/tests/tinytest.R index 3d8963c3..0e78dc56 100644 --- a/tests/tinytest.R +++ b/tests/tinytest.R @@ -1,5 +1,4 @@ -if ( requireNamespace("tinytest", quietly=TRUE) ){ +if (requireNamespace("tinytest", quietly = TRUE)) { tinytest::test_package("epiworldR") } - diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd index 90d39fb2..2ae8b68e 100644 --- a/vignettes/getting-started.Rmd +++ b/vignettes/getting-started.Rmd @@ -14,7 +14,7 @@ vignette: > ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, - comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5, + comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5, fig.align = "center" ) ``` @@ -23,33 +23,33 @@ knitr::opts_chunk$set( epiworldR is an R package that provides a fast (C++ backend) and highly- customizable framework for building network-based transmission/diffusion agent- -based models [ABM]. Some key features of epiworldR are the ability to construct -multi-virus models (e.g., models of competing multi-pathogens/multi-rumor,) -design mutating pathogens, architect population-level interventions, and build -models with an arbitrary number of compartments/states (beyond SIR/SEIR.) +based models [ABM]. Some key features of epiworldR are the ability to construct +multi-virus models (e.g., models of competing multi-pathogens/multi-rumor,) +design mutating pathogens, architect population-level interventions, and build +models with an arbitrary number of compartments/states (beyond SIR/SEIR.) # Example 1: Simulating an SIR model ## Setup and running the model -This example implements a social network with parameters listed within the -`ModelSIRCONN` function. The virus name is specified (COVID-19), 50000 agents -are initialized, the virus prevalence of 0.001 is declared, each agent will -contact two others (contact_rate), the transmission rate for +This example implements a social network with parameters listed within the +`ModelSIRCONN` function. The virus name is specified (COVID-19), 50000 agents +are initialized, the virus prevalence of 0.001 is declared, each agent will +contact two others (contact_rate), the transmission rate for any given agent is 0.3, and the recovery rate is set to $\frac{1}{3}$. -To create this model on epiworldR, simply use the `ModelSIRCONN` function. +To create this model on epiworldR, simply use the `ModelSIRCONN` function. From here, the example will take you through the basic features of epiworldR. ```{r sirconn-setup} library(epiworldR) model_sir <- ModelSIRCONN( name = "COVID-19", - n = 50000, - prevalence = 0.0001, + n = 50000, + prevalence = 0.0001, contact_rate = 2, transmission_rate = 0.5, - recovery_rate = 1/3 - ) + recovery_rate = 1 / 3 +) # Printing the model model_sir @@ -61,65 +61,65 @@ Printing the model shows us some information. Nevertheless, we can extract detai summary(model_sir) ``` -First, the name of the model, population size, number of entities (think of these as public spaces in which agents can make social contact with one another), the duration in days, number of viruses, amount of time the last replicate took to run (last run elapsed t), and rewiring -status (on or off). The model also includes a list of global actions (interventions) that are called during the model run. Next, you will see a list of the -viruses used in the model. In this case, COVID-19 was the only virus used. -Note that epiworldR can include more than one virus in a +First, the name of the model, population size, number of entities (think of these as public spaces in which agents can make social contact with one another), the duration in days, number of viruses, amount of time the last replicate took to run (last run elapsed t), and rewiring +status (on or off). The model also includes a list of global actions (interventions) that are called during the model run. Next, you will see a list of the +viruses used in the model. In this case, COVID-19 was the only virus used. +Note that epiworldR can include more than one virus in a model. Tool(s) lists agents' tools to fight the virus. Examples of -this may include masking, vaccines, social distancing, etc. In this model, no -tools are specified. Lastly, model parameters are listed. +this may include masking, vaccines, social distancing, etc. In this model, no +tools are specified. Lastly, model parameters are listed. + +To execute the model, use the run function with the SIR model object, the number of +simulation days, and an optional seed for reproducibility. Next, print out the +results from the simulated model using model_sir. -To execute the model, use the run function with the SIR model object, the number of -simulation days, and an optional seed for reproducibility. Next, print out the -results from the simulated model using model_sir. - -```{r} +```{r} run(model_sir, ndays = 50, seed = 1912) summary(model_sir) ``` ```{r getting-totals, echo=FALSE} -initials <- get_hist_total(model_sir)[1:3,]$counts |> prettyNum(big.mark = ",") +initials <- get_hist_total(model_sir)[1:3, ]$counts |> prettyNum(big.mark = ",") finals <- get_today_total(model_sir) |> prettyNum(big.mark = ",") tmat <- get_transition_probability(model_sir) tmat <- round(tmat, digits = 2) ``` -There are two additional sections included in the summary after running the model. First, we see the distribution of the population at time 50. This section describes the flow of agents from each state (SIR) after 50 days. In the +There are two additional sections included in the summary after running the model. First, we see the distribution of the population at time 50. This section describes the flow of agents from each state (SIR) after 50 days. In the example, you'll see the number of agents in the susceptible state decreased -from `r initials[1]` to `r finals[1]`, the number of agents in the infected state increased from `r initials[2]` to `r finals[2]`, and recovered agents increased to `r finals[3]` after 50 days. -The counts for these states will change based on model parameters or +from `r initials[1]` to `r finals[1]`, the number of agents in the infected state increased from `r initials[2]` to `r finals[2]`, and recovered agents increased to `r finals[3]` after 50 days. +The counts for these states will change based on model parameters or simulation run-time. The transmission probabilities section outputs a 3x3 matrix -that describes the probability of moving from one state to another. For example, -in the susceptible row, each agent has a `r tmat[1]` probability of remaining in the susceptible state with a `r tmat[1,2]` probability of moving from the susceptible state to the infected state. Notice there is no chance of skipping states. In other words, an agent can't jump from the -susceptible state to the recovered state; that agent must pass through the +that describes the probability of moving from one state to another. For example, +in the susceptible row, each agent has a `r tmat[1]` probability of remaining in the susceptible state with a `r tmat[1,2]` probability of moving from the susceptible state to the infected state. Notice there is no chance of skipping states. In other words, an agent can't jump from the +susceptible state to the recovered state; that agent must pass through the infected state to progress to the recovered state. The same logic -applies to moving backward; an agent cannot become susceptible again after -infection. - +applies to moving backward; an agent cannot become susceptible again after +infection. + ## Plot ## Extracting information -After running the epiworldR model, below is a list of all the functions that can -be called using the epiworld model object. +After running the epiworldR model, below is a list of all the functions that can +be called using the epiworld model object. ```{r showing-methods} methods(class = "epiworld_model") ``` -To demonstrate, start with the basic plot and get_hist_total functions. +To demonstrate, start with the basic plot and get_hist_total functions. ```{r} plot(model_sir) ``` -As evident from the above plot, the SIR model constructed from epiworldR -displays the changes in susceptible, infected, and recovered case counts over -time (days). Notice after a certain amount of time, the curves flatten. Below, -a table representation of the above plot is printed, complete with each state -within the SIR model, date, and agent counts. +As evident from the above plot, the SIR model constructed from epiworldR +displays the changes in susceptible, infected, and recovered case counts over +time (days). Notice after a certain amount of time, the curves flatten. Below, +a table representation of the above plot is printed, complete with each state +within the SIR model, date, and agent counts. ```{r get-hist-total} head(get_hist_total(model_sir)) @@ -132,8 +132,8 @@ repnum <- get_reproductive_number(model_sir) head(repnum) ``` -epiworldR has a method to plot the reproductive number automatically. -The function takes the average of values in the above table for each date and +epiworldR has a method to plot the reproductive number automatically. +The function takes the average of values in the above table for each date and repeats until all data have been accounted for. ```{r} @@ -151,12 +151,12 @@ plot_incidence(model_sir) ## Adding more viruses/viruses -epiworldR supports multi-virus models. The below code gives instructions on -how to implement this. Using the `virus` function, give a name to the new -virus/virus with its corresponding probability of infecting any given agent. -In this example, `prob_infecting` is set to 1.0, making it highly contagious. -To officially add this new virus to the model, use the `add_virus` -function by calling the original epiworldR model object, the new virus, and +epiworldR supports multi-virus models. The below code gives instructions on +how to implement this. Using the `virus` function, give a name to the new +virus/virus with its corresponding probability of infecting any given agent. +In this example, `prob_infecting` is set to 1.0, making it highly contagious. +To officially add this new virus to the model, use the `add_virus` +function by calling the original epiworldR model object, the new virus, and the new virus' prevalence (which is set to 0.01 in this example). ```{r design-and-add} @@ -164,47 +164,47 @@ the new virus' prevalence (which is set to 0.01 in this example). flu <- virus( name = "Flu", prob_infecting = .3, prevalence = .0001, as_proportion = TRUE - ) +) # Adding the virus to the model add_virus(model_sir, flu) ``` After running the updated model with the new virus included for 50 days, the -output below describes the simulation. To confirm that the flu is included, +output below describes the simulation. To confirm that the flu is included, notice the presence of "Flu" in the Virus(es) section of the output. All other -output is interpretable as specified in previous sections. +output is interpretable as specified in previous sections. ```{r} run(model_sir, ndays = 50, seed = 1912) model_sir ``` -Plotting the previous model (including the flu) yields the following. Notice -the presence of two reproductive numbers plotted over time. Variant 0 refers -to COVID-19, and virus 1 refers to the flu. +Plotting the previous model (including the flu) yields the following. Notice +the presence of two reproductive numbers plotted over time. Variant 0 refers +to COVID-19, and virus 1 refers to the flu. ```{r, fig.height=10} repnum2 <- get_reproductive_number(model_sir) -op <- par(mfrow = c(2,1)) +op <- par(mfrow = c(2, 1)) plot(model_sir) -plot(repnum2, type="b") +plot(repnum2, type = "b") par(op) ``` ## Tools -Now, the implementation of tools to combat any viruses and viruses in the model -will be demonstrated. First, for the sake of simplicity, remove the flu virus -from the SIR model object (keep in mind the index for the flu virus in the +Now, the implementation of tools to combat any viruses and viruses in the model +will be demonstrated. First, for the sake of simplicity, remove the flu virus +from the SIR model object (keep in mind the index for the flu virus in the model object is 1). Next, provide parameters for the new tool using the `tool` -function. These parameters include the name of the tool, any reduction in -probabilities for the SIR model parameters, and increased probability of +function. These parameters include the name of the tool, any reduction in +probabilities for the SIR model parameters, and increased probability of recovery option. In order to add the tool to the SIR model, use the `add_tool` -function with the SIR model object, new tool, and prevalence of the tool. -In this example, assume that 85% of the population will have received the -vaccination. +function with the SIR model object, new tool, and prevalence of the tool. +In this example, assume that 85% of the population will have received the +vaccination. ```{r} # Removing the flu virus from the model @@ -216,7 +216,7 @@ vaccine <- tool( as_proportion = TRUE, susceptibility_reduction = .9, transmission_reduction = .5, - recovery_enhancer = .5, + recovery_enhancer = .5, death_reduction = .9 ) @@ -227,9 +227,8 @@ run(model_sir, ndays = 50, seed = 1231) ```{r curves-including-vaccine, fig.height=10} repnum3 <- get_reproductive_number(model_sir) -op <- par(mfrow = c(2,1)) +op <- par(mfrow = c(2, 1)) plot_incidence(model_sir) -plot(repnum3, type="b") +plot(repnum3, type = "b") par(op) ``` - diff --git a/vignettes/implementation.Rmd b/vignettes/implementation.Rmd index a72a83ca..dae17362 100644 --- a/vignettes/implementation.Rmd +++ b/vignettes/implementation.Rmd @@ -13,9 +13,9 @@ vignette: > # Introduction -The following vignette provides detailed information about the implementation of `epiworldR`. The package is a wrapper of the C++ package `epiworld`, a framework for building agent-based models.[^where-cpp] +The following vignette provides detailed information about the implementation of `epiworldR`. The package is a wrapper of the C++ package `epiworld`, a framework for building agent-based models.[^where-cpp] -[^where-cpp]: The C++ package is available at [https://github.com/UofUEpiBio/epiworld](https://github.com/UofUEpiBio/epiworld). +[^where-cpp]: The C++ package is available at [https://github.com/UofUEpiBio/epiworld](https://github.com/UofUEpiBio/epiworld). # General flow of the models @@ -31,12 +31,12 @@ The core function of `epiworldR` is the `run()` function. This function executes a. The state of each agent is updated. States are updated according to their corresponding `update_state` function. - + Since the model is discrete-time, **state changes are stored as promises**, meaning that agents' states are not updated immediately. Instead, the state is updated at the end of the updates. This is done to avoid updating the state of an agent and then using the updated state to update the state of another agent. For example, if agent $i$ infects agent $j$, then agent $j$ should not be able to infect agent $i$ in the same step. - + b. Once the update schedule is laid out, the changes are made effective, so, for instance, individuals who became infected during the update will start the next step in the infected state. - + c. Global actions are executed. These could also change agents' states, so just like in the previous step, these changes are stored as promises and made effective once all actions have been evaluated. d. The model's state is recorded in the database, and the `current_date` is incremented by 1. @@ -69,9 +69,9 @@ Viruses and tools provide a way to adjust how agents move between states. Viruses in `epiworldR` contain various baseline probabilities used across models, including transmission, recovery, and death. On the other hand, tools alter these probabilities by reducing/increasing -them. Furthermore, tools alter agents' susceptibility, infectiousness, +them. Furthermore, tools alter agents' susceptibility, infectiousness, recovery, and death probabilities. Currently, tools alter these probabilities -by a constant factor, +by a constant factor, $$ p_{ij} = p_{v} \times \left(1 - factor_{host}\right) \times \left(1 - factor_{target}\right) @@ -85,8 +85,8 @@ $$ factor_{agent} = 1 - \prod_{t\in tools_{agent}}\left(1 - factor_{t}\right) $$ -Therefore, for example, a vaccinated agent wearing a mask would have a factor of $1 - (1 - 0.30) \times (1 - 0.5) = 0.65$. -The adjusted probabilities principle also applies to recovery rates in the SIR and SEIR models. +Therefore, for example, a vaccinated agent wearing a mask would have a factor of $1 - (1 - 0.30) \times (1 - 0.5) = 0.65$. +The adjusted probabilities principle also applies to recovery rates in the SIR and SEIR models. ## Transmission in connected models[^models] @@ -98,4 +98,4 @@ The "connected" models provide a version where agents live in a fully connected 2. Then, $c$ agents are randomly selected from the population. Transmission can then occur from any of these agents to the susceptible agent. -3. The probability of transmission from each of the $c$ agents is calculated as described in the previous section. \ No newline at end of file +3. The probability of transmission from each of the $c$ agents is calculated as described in the previous section. diff --git a/vignettes/mixing.Rmd b/vignettes/mixing.Rmd index 7fb29c8f..ccdb36b8 100644 --- a/vignettes/mixing.Rmd +++ b/vignettes/mixing.Rmd @@ -40,7 +40,7 @@ This matrix represents the probability of an agent from population $i$ interacti We will build this model using the `entity` class in epiworld. The following code chunk instantiates three entities and the contact matrix. -```{r entity-matrix-setup} +```{r entity-matrix-setup} library(epiworldR) e1 <- entity("Population 1", 3e3, as_proportion = FALSE) @@ -71,7 +71,7 @@ flu_model <- ModelSEIRMixing( contact_matrix = cmatrix ) -# Adding the entities +# Adding the entities flu_model |> add_entity(e1) |> add_entity(e2) |> @@ -145,4 +145,4 @@ legend( col = c("black", "red", "blue"), lty = 1 ) -``` \ No newline at end of file +``` diff --git a/vignettes/run-multiple.Rmd b/vignettes/run-multiple.Rmd index a080d22b..34efebee 100644 --- a/vignettes/run-multiple.Rmd +++ b/vignettes/run-multiple.Rmd @@ -13,88 +13,85 @@ vignette: > ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, - comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5, + comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5, fig.align = "center" ) ``` # Introduction -The purpose of the "run_multiple" function is to run a specified number of +The purpose of the "run_multiple" function is to run a specified number of simulations using the same model object. That is, this function makes it possible to compare model results across several separate and repeated -simulations. +simulations. -# Example: Simulating a SEIRCONN Model 50 Times +# Example: Simulating a SEIRCONN Model 50 Times ## Setup and Running Model To use the "run_multiple" function in epiworld, create your epimodel of choice; -in this case, the example uses a SEIRCONN model for COVID-19, 100000 people, -an initial prevalence of 0.0001 (0.01%), a contact rate of 2, probability of -transmission 0.5, a total of 7 incubation days, and probability of recovery -$\frac{1}{3}$. +in this case, the example uses a SEIRCONN model for COVID-19, 100000 people, +an initial prevalence of 0.0001 (0.01%), a contact rate of 2, probability of +transmission 0.5, a total of 7 incubation days, and probability of recovery +$\frac{1}{3}$. ```{r sirconn-setup} library(epiworldR) model_seirconn <- ModelSEIRCONN( name = "COVID-19", - n = 10000, - prevalence = 0.0001, + n = 10000, + prevalence = 0.0001, contact_rate = 2, transmission_rate = 0.5, incubation_days = 7, - recovery_rate = 1/3 - ) - + recovery_rate = 1 / 3 +) ``` ## Generating a Saver -Next, generate a saver for the purpose of extracting the "total_hist" and -"reproductive" information from the model object. Now, use the -"run_multiple" function with the model object, number of desired days to run -the simulation, number of simulations to run, and number of threads for -parallel computing. +Next, generate a saver for the purpose of extracting the "total_hist" and +"reproductive" information from the model object. Now, use the +"run_multiple" function with the model object, number of desired days to run +the simulation, number of simulations to run, and number of threads for +parallel computing. ```{r saver-generation} # Generating a saver saver <- make_saver("total_hist", "reproductive") # Running and printing run_multiple(model_seirconn, ndays = 50, nsims = 50, saver = saver, nthreads = 2) - ``` -Using the "run_multiple_get_results" function, extract the results from the -model object that was simulated 50 times for comparison across simulations. +Using the "run_multiple_get_results" function, extract the results from the +model object that was simulated 50 times for comparison across simulations. ```{r retrieving results} # Retrieving the results ans <- run_multiple_get_results(model_seirconn) head(ans$total_hist) head(ans$reproductive) - ``` ## Plotting -To plot the epicurves and reproductive numbers over time using boxplots, extract -the results from the model object using "run_multiple_get_results". For this -example, the dates are filtered to be less than or equal to 20 to observe the -epicurves in the first 20 days. Notice each boxplot in the below table -represents the observed values from each of the 50 simulations for each date. +To plot the epicurves and reproductive numbers over time using boxplots, extract +the results from the model object using "run_multiple_get_results". For this +example, the dates are filtered to be less than or equal to 20 to observe the +epicurves in the first 20 days. Notice each boxplot in the below table +represents the observed values from each of the 50 simulations for each date. ```{r plotting seirconn epicurves} seirconn_50 <- run_multiple_get_results(model_seirconn)$total_hist -seirconn_50 <- seirconn_50[seirconn_50$date <= 20,] +seirconn_50 <- seirconn_50[seirconn_50$date <= 20, ] plot(seirconn_50) ``` - -To view the a plot of the reproductive number over all 50 days for each of the -50 simulations, store the reproductive results to a new object using -"run_multiple_get_results", then plot using the "boxplot" function. Notice + +To view the a plot of the reproductive number over all 50 days for each of the +50 simulations, store the reproductive results to a new object using +"run_multiple_get_results", then plot using the "boxplot" function. Notice each source exposure date displays a boxplot representing the distribution of reproductive numbers across all 50 simulations. As expected, the reproductive -number on average, decreases over time. +number on average, decreases over time. ```{r reproductive number plot} seirconn_50_r <- run_multiple_get_results(model_seirconn)$reproductive @@ -105,7 +102,4 @@ plot(seirconn_50_r) # ylab = "rt", # border = "black", # las = 2) - ``` - -