diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile
new file mode 100644
index 00000000..6319e186
--- /dev/null
+++ b/.devcontainer/Dockerfile
@@ -0,0 +1,15 @@
+FROM rocker/r-base:latest
+
+RUN \
+ echo 'options(repos=c(CRAN="https://cloud.r-project.org"))' >> ~/.Rprofile && \
+ Rscript --vanilla -e 'getOption("repos")'
+
+# Adding Git
+RUN apt-get update && apt-get install -y --no-install-recommends git
+
+# Adding R packages
+RUN \
+ wget https://github.com/jgm/pandoc/releases/download/3.2.1/pandoc-3.2.1-1-amd64.deb && \
+ dpkg -i pandoc-3.2.1-1-amd64.deb
+
+RUN install2.r cpp11 rmarkdown roxygen2 tinytest data.table netplot
\ No newline at end of file
diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
new file mode 100644
index 00000000..d2974f1d
--- /dev/null
+++ b/.devcontainer/devcontainer.json
@@ -0,0 +1,23 @@
+// For format details, see https://aka.ms/devcontainer.json. For config options, see the
+// README at: https://github.com/devcontainers/templates/tree/main/src/cpp
+{
+ "name": "epiworldR",
+ "build": {
+ "dockerfile": "Dockerfile"
+ },
+
+ // Features to add to the dev container. More info: https://containers.dev/features.
+ // "features": {},
+
+ // Use 'forwardPorts' to make a list of ports inside the container available locally.
+ // "forwardPorts": [],
+
+ // Use 'postCreateCommand' to run commands after the container is created.
+ // "postCreateCommand": "gcc -v",
+
+ // Configure tool-specific properties.
+ // "customizations": {},
+
+ // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root.
+ "remoteUser": "root"
+}
diff --git a/.github/workflows/r.yml b/.github/workflows/r.yml
index 58ac8b27..15336ff2 100644
--- a/.github/workflows/r.yml
+++ b/.github/workflows/r.yml
@@ -36,7 +36,7 @@ jobs:
R_KEEP_PKG_SOURCE: yes
steps:
- - uses: actions/checkout@v2
+ - uses: actions/checkout@v4
- uses: r-lib/actions/setup-pandoc@v2
@@ -69,7 +69,7 @@ jobs:
error-on: '"warning"'
# Upload the built package as an artifact
- - uses: actions/upload-artifact@v2
+ - uses: actions/upload-artifact@v4
if: ${{ matrix.config.os == 'ubuntu-latest' && matrix.config.r == 'release' }}
with:
name: ${{ matrix.config.os }}-pkg
diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json
index a52fc3db..eb3d24ec 100644
--- a/.vscode/c_cpp_properties.json
+++ b/.vscode/c_cpp_properties.json
@@ -7,7 +7,8 @@
"/usr/local/include",
"/usr/lib/R/site-library/cpp11/include",
"/usr/lib/R/site-library/Rcpp/include",
- "/usr/share/R/include"
+ "/usr/share/R/include",
+ "inst/include/epiworld"
],
"intelliSenseMode": "linux-gcc-x64",
"compilerPath": "/usr/bin/gcc",
diff --git a/DESCRIPTION b/DESCRIPTION
index 9c882d09..d67fb3b2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: epiworldR
Type: Package
Title: Fast Agent-Based Epi Models
-Version: 0.1-1
-Date: 2024-04-07
+Version: 0.3-2
+Date: 2024-06-12
Authors@R: c(
person("Derek", "Meyer", role=c("aut","cre"),
email="derekmeyer37@gmail.com", comment = c(ORCID = "0009-0005-1350-6988")),
@@ -21,7 +21,7 @@ URL: https://github.com/UofUEpiBio/epiworldR,
https://uofuepibio.github.io/epiworldR-workshop/
BugReports: https://github.com/UofUEpiBio/epiworldR/issues
License: MIT + file LICENSE
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
LinkingTo: cpp11
diff --git a/Makefile b/Makefile
index 21139d01..7dd0d23e 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,7 @@
# Capture the current value of the version of the package in DESCRIPTION
VERSION := $(shell grep Version DESCRIPTION | sed -e 's/Version: //')
+
build:
cd .. && R CMD build epiworldR
diff --git a/NAMESPACE b/NAMESPACE
index 887ea275..17f23b74 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,8 +1,8 @@
# Generated by roxygen2: do not edit by hand
S3method("[",epiworld_agents)
+S3method("[",epiworld_entities)
S3method(add_tool,epiworld_model)
-S3method(add_tool_n,epiworld_model)
S3method(add_virus,epiworld_model)
S3method(add_virus,epiworld_seir)
S3method(add_virus,epiworld_seirconn)
@@ -12,14 +12,10 @@ S3method(add_virus,epiworld_sir)
S3method(add_virus,epiworld_sirconn)
S3method(add_virus,epiworld_sird)
S3method(add_virus,epiworld_sirdconn)
-S3method(add_virus_n,epiworld_model)
-S3method(add_virus_n,epiworld_seir)
-S3method(add_virus_n,epiworld_seirconn)
-S3method(add_virus_n,epiworld_sir)
-S3method(add_virus_n,epiworld_sirconn)
S3method(agents_from_edgelist,epiworld_model)
S3method(agents_smallworld,epiworld_model)
S3method(as.array,epiworld_hist_transition)
+S3method(get_agents,epiworld_model)
S3method(get_hist_tool,epiworld_model)
S3method(get_hist_total,epiworld_model)
S3method(get_hist_transition_matrix,epiworld_model)
@@ -48,16 +44,23 @@ S3method(plot,epiworld_seir)
S3method(plot,epiworld_seirconn)
S3method(plot,epiworld_seird)
S3method(plot,epiworld_seirdconn)
+S3method(plot,epiworld_seirmixing)
S3method(plot,epiworld_sir)
S3method(plot,epiworld_sirconn)
S3method(plot,epiworld_sird)
S3method(plot,epiworld_sirdconn)
+S3method(plot,epiworld_sirmixing)
S3method(plot,epiworld_sis)
S3method(plot,epiworld_sisd)
S3method(plot,epiworld_surv)
+S3method(plot_epi,epiworld_hist)
+S3method(plot_epi,epiworld_hist_virus)
+S3method(plot_epi,epiworld_model)
S3method(print,epiworld_agent)
S3method(print,epiworld_agents)
S3method(print,epiworld_agents_tools)
+S3method(print,epiworld_entities)
+S3method(print,epiworld_entity)
S3method(print,epiworld_globalevent)
S3method(print,epiworld_model)
S3method(print,epiworld_saver)
@@ -82,14 +85,17 @@ export(ModelSEIR)
export(ModelSEIRCONN)
export(ModelSEIRD)
export(ModelSEIRDCONN)
+export(ModelSEIRMixing)
export(ModelSIR)
export(ModelSIRCONN)
export(ModelSIRD)
export(ModelSIRDCONN)
export(ModelSIRLogit)
+export(ModelSIRMixing)
export(ModelSIS)
export(ModelSISD)
export(ModelSURV)
+export(add_entity)
export(add_globalevent)
export(add_tool)
export(add_tool_agent)
@@ -101,10 +107,22 @@ export(agents_from_edgelist)
export(agents_smallworld)
export(change_state)
export(clone_model)
+export(distribute_entity_randomly)
+export(distribute_entity_to_set)
+export(distribute_tool_randomly)
+export(distribute_tool_to_set)
+export(distribute_virus_randomly)
+export(distribute_virus_set)
+export(entity)
+export(entity_add_agent)
+export(entity_get_agents)
export(get_agents)
export(get_agents_data_ncols)
export(get_agents_states)
export(get_agents_tools)
+export(get_entities)
+export(get_entity_name)
+export(get_entity_size)
export(get_generation_time)
export(get_hist_tool)
export(get_hist_total)
@@ -138,12 +156,14 @@ export(globalevent_tool_logit)
export(has_tool)
export(has_virus)
export(initial_states)
+export(load_agents_entities_ties)
export(make_saver)
export(plot_generation_time)
export(plot_incidence)
export(plot_reproductive_number)
export(queuing_off)
export(queuing_on)
+export(rm_entity)
export(rm_tool)
export(rm_virus)
export(run)
@@ -153,6 +173,9 @@ export(set_agents_data)
export(set_death_reduction)
export(set_death_reduction_fun)
export(set_death_reduction_ptr)
+export(set_distribution_entity)
+export(set_distribution_tool)
+export(set_distribution_virus)
export(set_incubation)
export(set_incubation_fun)
export(set_incubation_ptr)
diff --git a/NEWS.md b/NEWS.md
index c949aca9..2f15b86e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,22 @@
+# epiworldR 0.3-2 (dev)
+
+* Starting version 0.3-0, `epiworldR` is versioned using the same version as the C++ library, `epiworld`.
+
+* Adds the new mixing models `ModelSIRMixing` and `ModelSEIRMixing`.
+
+* 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
+ specified in the object itself. `add_tool_n` and `add_virus_n` are now
+ deprecated.
+
+* `globalaction_*` are now defunct. Use `globalevent_*` instead.
+
+* New functions to specify how viruses, tools, and entities are distributed
+ among agents: `distribute_viruses`, `distribute_tools`, and `distribute_entities`.
+
+
# epiworldR 0.1-0`
* Force model to update agents' states when running a simulation.
diff --git a/R/ModelSEIRCONN.R b/R/ModelSEIRCONN.R
index c5a41cea..3e01bff1 100644
--- a/R/ModelSEIRCONN.R
+++ b/R/ModelSEIRCONN.R
@@ -39,8 +39,8 @@
#' plot(model_seirconn)
#'
#' # Adding the flu
-#' flu <- virus("Flu", .9, 1/7)
-#' add_virus(model_seirconn, flu, .001)
+#' 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)
diff --git a/R/ModelSEIRDCONN.R b/R/ModelSEIRDCONN.R
index 44aed7ba..42290d5a 100644
--- a/R/ModelSEIRDCONN.R
+++ b/R/ModelSEIRDCONN.R
@@ -47,8 +47,12 @@
#' plot(model_seirdconn)
#'
#' # Adding the flu
-#' flu <- virus("Flu", prob_infecting = .3, recovery_rate = 1/7, prob_death = 0.001)
-#' add_virus(model = model_seirdconn, virus = flu, proportion = .001)
+#' flu <- virus(
+#' "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)
diff --git a/R/ModelSEIRMixing.R b/R/ModelSEIRMixing.R
new file mode 100644
index 00000000..1dee4eac
--- /dev/null
+++ b/R/ModelSEIRMixing.R
@@ -0,0 +1,91 @@
+#' Susceptible Exposed Infected Removed model (SEIR) with mixing
+#' @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
+#' 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 ... Currently ignore.
+#' @param n Number of individuals in the population.
+#' @param contact_matrix Matrix of contact rates between individuals.
+#' @export
+#' @family Models
+#' @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.
+#' 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
+#'
+#' # 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,
+#' prevalence = 1 / N,
+#' contact_rate = 20,
+#' transmission_rate = 0.1,
+#' recovery_rate = 1 / 7,
+#' 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,
+ incubation_days, recovery_rate, contact_matrix
+) {
+
+ structure(
+ ModelSEIRMixing_cpp(
+ 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
+#' [epiworld_model].
+#' @param main Title of the plot
+plot.epiworld_seirmixing <- function(x, main = get_name(x), ...) { # col = NULL
+ plot_epi(x, main = main, ...)
+}
diff --git a/R/ModelSIRMixing.R b/R/ModelSIRMixing.R
new file mode 100644
index 00000000..d10b5c29
--- /dev/null
+++ b/R/ModelSIRMixing.R
@@ -0,0 +1,88 @@
+#' Susceptible Infected Removed model (SIR) with mixing
+#' @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
+#' transmission.
+#' @param recovery_rate Numeric scalar between 0 and 1. Probability of recovery.
+#' @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
+#' 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.
+#' 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
+#' # 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,
+#' prevalence = 1 / N,
+#' contact_rate = 20,
+#' transmission_rate = 0.1,
+#' 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,
+ 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
+#' [epiworld_model].
+#' @param main Title of the plot
+plot.epiworld_sirmixing <- function(x, main = get_name(x), ...) { # col = NULL
+ plot_epi(x, main = main, ...)
+}
diff --git a/R/agents-methods.R b/R/agents-methods.R
index d71148c5..2c1bd5fb 100644
--- a/R/agents-methods.R
+++ b/R/agents-methods.R
@@ -37,7 +37,12 @@
#'
#' x[0] # Print information about the first agent. Substitute the agent of
#' # interest's position where '0' is.
-get_agents <- function(model) {
+#' @name agents
+get_agents <- function(model, ...) UseMethod("get_agents")
+
+#' @export
+#' @rdname agents
+get_agents.epiworld_model <- function(model, ...) {
res <- get_agents_cpp(model)
@@ -52,7 +57,7 @@ get_agents <- function(model) {
#' @param x An object of class [epiworld_agents]
#' @param i Index (id) of the agent (from 0 to `n-1`)
#' @export
-#' @rdname get_agents
+#' @rdname agents
#' @return
#' - The `[` method returns an object of class [epiworld_agent].
#' @aliases epiworld_agent
@@ -82,7 +87,7 @@ get_agents <- function(model) {
#' @returns
#' - The `print` function returns information about each individual agent of
#' class [epiworld_agent].
-#' @rdname get_agents
+#' @rdname agents
print.epiworld_agent <- function(x, compressed = FALSE, ...) {
invisible(print_agent_cpp(x, attr(x, "model"), compressed))
@@ -91,7 +96,7 @@ print.epiworld_agent <- function(x, compressed = FALSE, ...) {
#' @export
#' @param max_print Integer scalar. Maximum number of agents to print.
-#' @rdname get_agents
+#' @rdname agents
print.epiworld_agents <- function(x, compressed = TRUE, max_print = 10, ...) {
model <- attr(x, "model")
@@ -113,7 +118,7 @@ print.epiworld_agents <- function(x, compressed = TRUE, max_print = 10, ...) {
#' @export
#' @returns
#' - The `get_state` function returns the state of the [epiworld_agents] object.
-#' @rdname get_agents
+#' @rdname agents
get_state <- function(x) {
get_state_agent_cpp(x)
}
diff --git a/R/agents.R b/R/agents.R
index 9d2e9ca3..f4b4f8a1 100644
--- a/R/agents.R
+++ b/R/agents.R
@@ -19,7 +19,6 @@ stopifnot_agent <- function(x) {
#' @param d,directed Logical scalar. Whether the graph is directed or not.
#' @param p Probability of rewiring.
#' @export
-#' @aliases agents
#' @return
#' - The 'agents_smallworld' function returns a model with the agents
#' loaded.
diff --git a/R/cpp11.R b/R/cpp11.R
index 0b5eb2c2..b1276fdb 100644
--- a/R/cpp11.R
+++ b/R/cpp11.R
@@ -104,6 +104,62 @@ get_today_total_cpp <- function(model) {
.Call(`_epiworldR_get_today_total_cpp`, model)
}
+get_entities_cpp <- function(model) {
+ .Call(`_epiworldR_get_entities_cpp`, model)
+}
+
+get_entity_cpp <- function(entities, idx) {
+ .Call(`_epiworldR_get_entity_cpp`, entities, idx)
+}
+
+entity_cpp <- function(name, preval, as_proportion, to_unassigned) {
+ .Call(`_epiworldR_entity_cpp`, name, preval, as_proportion, to_unassigned)
+}
+
+get_entity_size_cpp <- function(entity) {
+ .Call(`_epiworldR_get_entity_size_cpp`, entity)
+}
+
+entity_add_agent_cpp <- function(entity, agent, model) {
+ .Call(`_epiworldR_entity_add_agent_cpp`, entity, agent, model)
+}
+
+get_entity_name_cpp <- function(entity) {
+ .Call(`_epiworldR_get_entity_name_cpp`, entity)
+}
+
+add_entity_cpp <- function(model, entity) {
+ .Call(`_epiworldR_add_entity_cpp`, model, entity)
+}
+
+rm_entity_cpp <- function(model, entity_pos) {
+ .Call(`_epiworldR_rm_entity_cpp`, model, entity_pos)
+}
+
+load_agents_entities_ties_cpp <- function(model, agents_ids, entities_ids) {
+ .Call(`_epiworldR_load_agents_entities_ties_cpp`, model, agents_ids, entities_ids)
+}
+
+entity_get_agents_cpp <- function(entity) {
+ .Call(`_epiworldR_entity_get_agents_cpp`, entity)
+}
+
+print_entity_cpp <- function(entity) {
+ .Call(`_epiworldR_print_entity_cpp`, entity)
+}
+
+set_distribution_entity_cpp <- function(entity, fun) {
+ .Call(`_epiworldR_set_distribution_entity_cpp`, entity, fun)
+}
+
+distribute_entity_randomly_cpp <- function(prevalence, as_proportion, to_unassigned) {
+ .Call(`_epiworldR_distribute_entity_randomly_cpp`, prevalence, as_proportion, to_unassigned)
+}
+
+distribute_entity_to_set_cpp <- function(agents_ids) {
+ .Call(`_epiworldR_distribute_entity_to_set_cpp`, agents_ids)
+}
+
ModelSURV_cpp <- function(name, prevalence, efficacy_vax, latent_period, prob_symptoms, prop_vaccinated, prop_vax_redux_transm, infect_period, prop_vax_redux_infect, surveillance_prob, transmission_rate, prob_death, prob_noreinfect) {
.Call(`_epiworldR_ModelSURV_cpp`, name, prevalence, efficacy_vax, latent_period, prob_symptoms, prop_vaccinated, prop_vax_redux_transm, infect_period, prop_vax_redux_infect, surveillance_prob, transmission_rate, prob_death, prob_noreinfect)
}
@@ -156,6 +212,14 @@ ModelDiffNet_cpp <- function(name, prevalence, prob_adopt, normalize_exposure, d
.Call(`_epiworldR_ModelDiffNet_cpp`, name, prevalence, prob_adopt, normalize_exposure, data, data_ncols, data_cols, params)
}
+ModelSIRMixing_cpp <- function(name, n, prevalence, contact_rate, transmission_rate, recovery_rate, contact_matrix) {
+ .Call(`_epiworldR_ModelSIRMixing_cpp`, name, n, prevalence, contact_rate, transmission_rate, recovery_rate, contact_matrix)
+}
+
+ModelSEIRMixing_cpp <- function(name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate, contact_matrix) {
+ .Call(`_epiworldR_ModelSEIRMixing_cpp`, name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate, contact_matrix)
+}
+
print_cpp <- function(m, lite) {
.Call(`_epiworldR_print_cpp`, m, lite)
}
@@ -264,16 +328,12 @@ clone_model_cpp <- function(model) {
.Call(`_epiworldR_clone_model_cpp`, model)
}
-tool_cpp <- function(name, susceptibility_reduction, transmission_reduction, recovery_enhancer, death_reduction) {
- .Call(`_epiworldR_tool_cpp`, name, susceptibility_reduction, transmission_reduction, recovery_enhancer, death_reduction)
-}
-
-add_tool_cpp <- function(m, t, preval) {
- .Call(`_epiworldR_add_tool_cpp`, m, t, preval)
+tool_cpp <- function(name, prevalence, as_proportion, susceptibility_reduction, transmission_reduction, recovery_enhancer, death_reduction) {
+ .Call(`_epiworldR_tool_cpp`, name, prevalence, as_proportion, susceptibility_reduction, transmission_reduction, recovery_enhancer, death_reduction)
}
-add_tool_n_cpp <- function(m, t, preval) {
- .Call(`_epiworldR_add_tool_n_cpp`, m, t, preval)
+add_tool_cpp <- function(m, t) {
+ .Call(`_epiworldR_add_tool_cpp`, m, t)
}
rm_tool_cpp <- function(m, tool_pos) {
@@ -352,20 +412,28 @@ print_agent_tools_cpp <- function(tools) {
.Call(`_epiworldR_print_agent_tools_cpp`, tools)
}
-virus_cpp <- function(name, prob_infecting, prob_recovery, prob_death, post_immunity, incubation) {
- .Call(`_epiworldR_virus_cpp`, name, prob_infecting, prob_recovery, prob_death, post_immunity, incubation)
+set_distribution_tool_cpp <- function(tool, distfun) {
+ .Call(`_epiworldR_set_distribution_tool_cpp`, tool, distfun)
}
-virus_set_state_cpp <- function(v, init, end, removed) {
- .Call(`_epiworldR_virus_set_state_cpp`, v, init, end, removed)
+distribute_tool_randomly_cpp <- function(prevalence, as_proportion) {
+ .Call(`_epiworldR_distribute_tool_randomly_cpp`, prevalence, as_proportion)
}
-add_virus_cpp <- function(m, v, preval) {
- .Call(`_epiworldR_add_virus_cpp`, m, v, preval)
+distribute_tool_to_set_cpp <- function(agents_ids) {
+ .Call(`_epiworldR_distribute_tool_to_set_cpp`, agents_ids)
}
-add_virus_n_cpp <- function(m, v, preval) {
- .Call(`_epiworldR_add_virus_n_cpp`, m, v, preval)
+virus_cpp <- function(name, prevalence, as_proportion, prob_infecting, prob_recovery, prob_death, post_immunity, incubation) {
+ .Call(`_epiworldR_virus_cpp`, name, prevalence, as_proportion, prob_infecting, prob_recovery, prob_death, post_immunity, incubation)
+}
+
+virus_set_state_cpp <- function(v, init, end, removed) {
+ .Call(`_epiworldR_virus_set_state_cpp`, v, init, end, removed)
+}
+
+add_virus_cpp <- function(m, v) {
+ .Call(`_epiworldR_add_virus_cpp`, m, v)
}
rm_virus_cpp <- function(m, virus_pos) {
@@ -435,3 +503,15 @@ get_name_virus_cpp <- function(virus) {
set_name_virus_cpp <- function(virus, name) {
.Call(`_epiworldR_set_name_virus_cpp`, virus, name)
}
+
+set_distribution_virus_cpp <- function(virus, dist) {
+ .Call(`_epiworldR_set_distribution_virus_cpp`, virus, dist)
+}
+
+distribute_virus_randomly_cpp <- function(prevalence, as_proportion) {
+ .Call(`_epiworldR_distribute_virus_randomly_cpp`, prevalence, as_proportion)
+}
+
+distribute_virus_to_set_cpp <- function(agents_ids) {
+ .Call(`_epiworldR_distribute_virus_to_set_cpp`, agents_ids)
+}
diff --git a/R/entity.R b/R/entity.R
new file mode 100644
index 00000000..e10c4227
--- /dev/null
+++ b/R/entity.R
@@ -0,0 +1,280 @@
+
+stopifnot_entity <- function(entity) {
+ if (!inherits(entity, "epiworld_entity")) {
+ stop("Argument 'entity' must be an entity object.")
+ }
+}
+
+stopifnot_entity_distfun <- function(distfun) {
+ if (!inherits(distfun, "epiworld_distribution_entity")) {
+ stop("Argument 'distfun' must be a distribution function.")
+ }
+}
+
+#' Get entities
+#'
+#' Entities in `epiworld` are objects that can contain agents.
+#' @param model Model object of class `epiworld_model`.
+#'
+#' @details
+#' Epiworld entities are especially useful for mixing models, particularly
+#' [ModelSIRMixing] and [ModelSEIRMixing].
+#'
+#' @name entities
+#' @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)
+#' )
+#'
+#' 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) {
+
+ stopifnot_model(model)
+ structure(
+ lapply(
+ get_entities_cpp(model), \(e) {
+ structure(
+ e,
+ class = c("epiworld_entity"),
+ model = model
+ )
+ }
+ ),
+ class = c("epiworld_entities")
+ )
+}
+
+#' @export
+print.epiworld_entities <- function(x, ...) {
+ cat("A collection of ", length(x), " entities.\n")
+ invisible(x)
+}
+
+#' @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.")
+ }
+
+ structure(
+ get_entity_cpp(x, i),
+ class = c("epiworld_entity"),
+ model = x$model
+ )
+
+}
+
+#' @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
+#' - The function `entity` creates an entity object.
+#' @rdname entities
+entity <- function(name, prevalence, as_proportion, to_unassigned = TRUE) {
+
+ structure(
+ entity_cpp(
+ name,
+ as.double(prevalence),
+ as.logical(as_proportion),
+ as.logical(to_unassigned)
+ ),
+ class = "epiworld_entity"
+ )
+
+}
+
+#' @export
+#' @rdname entities
+#' @param entity Entity object of class `epiworld_entity`.
+#' @return
+#' - The function `get_entity_size` returns the number of agents in the entity.
+get_entity_size <- function(entity) {
+ stopifnot_entity(entity)
+ get_entity_size_cpp(entity)
+}
+
+#' @export
+#' @rdname entities
+#' @return
+#' - The function `get_entity_name` returns the name of the entity.
+get_entity_name <- function(entity) {
+ stopifnot_entity(entity)
+ get_entity_name_cpp(entity)
+}
+
+#' @export
+#' @rdname entities
+#' @param agent Agent object of class `epiworld_agent`.
+#' @return
+#' - The function `entity_add_agent` adds an agent to the entity.
+entity_add_agent <- function(
+ entity,
+ agent,
+ model = attr(entity, "model")
+ ) {
+
+ stopifnot_entity(entity)
+ stopifnot_agent(agent)
+ entity_add_agent_cpp(entity, agent, model)
+
+ invisible(entity)
+
+}
+
+#' @export
+#' @rdname entities
+#' @param id Integer scalar. Entity id to remove (starting from zero).
+#' @return
+#' - The function `rm_entity` removes an entity from the model.
+rm_entity <- function(model, id) {
+
+ stopifnot_model(model)
+ rm_entity_cpp(model, entity)
+
+ invisible(model)
+}
+
+#' @export
+#' @rdname entities
+add_entity <- function(
+ model,
+ entity
+) {
+
+ stopifnot_model(model)
+ stopifnot_entity(entity)
+ add_entity_cpp(
+ model,
+ entity
+ )
+
+ invisible(model)
+
+}
+
+#' @export
+#' @rdname entities
+#' @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
+) {
+
+ stopifnot_model(model)
+ if (!inherits(agents_id, "integer")) {
+ stop("Argument 'agents_id' must be an integer.")
+ }
+
+ if (!inherits(entities_id, "integer")) {
+ stop("Argument 'entities_id' must be an integer.")
+ }
+
+ load_agents_entities_ties_cpp(model, agents_id, entities_id)
+
+ invisible(model)
+
+}
+
+#' @export
+#' @rdname entities
+#' @return
+#' - The function `entity_get_agents` returns an integer vector with the agents
+#' in the entity (ids).
+entity_get_agents <- function(entity) {
+
+ stopifnot_entity(entity)
+ entity_get_agents_cpp(entity)
+
+}
+
+#' @export
+print.epiworld_entity <- function(x, ...) {
+ print_entity_cpp(x)
+ invisible(x)
+}
+
+#' @export
+#' @param prevalence Numeric scalar. Prevalence of the entity.
+#' @param as_proportion Logical scalar. If `TRUE`, `prevalence` is interpreted
+#' as a proportion.
+#' @rdname entities
+distribute_entity_randomly <- function(
+ prevalence,
+ as_proportion,
+ to_unassigned = TRUE
+) {
+
+ structure(
+ distribute_entity_randomly_cpp(
+ as.double(prevalence),
+ as.logical(as_proportion),
+ as.logical(to_unassigned)
+ ),
+ class = "epiworld_distribution_entity"
+ )
+
+}
+
+#' @export
+#' @param agents_ids Integer vector. Ids of the agents to distribute.
+#' @rdname entities
+distribute_entity_to_set <- function(
+ agents_ids
+) {
+
+ structure(
+ distribute_entity_to_set_cpp(
+ as.integer(agents_ids)
+ ),
+ class = "epiworld_distribution_entity"
+ )
+
+}
+
+#' @export
+#' @rdname entities
+#' @param distfun Distribution function object of class `epiworld_distribution_entity`.
+set_distribution_entity <- function(
+ entity,
+ distfun
+) {
+
+ stopifnot_entity(entity)
+ stopifnot_entity_distfun(distfun)
+ set_distribution_entity_cpp(entity, distfun)
+
+ invisible(entity)
+
+}
\ No newline at end of file
diff --git a/R/functions-renamed.R b/R/functions-renamed.R
index 2726e6a6..61556bed 100644
--- a/R/functions-renamed.R
+++ b/R/functions-renamed.R
@@ -1,7 +1,48 @@
-#' Deprecated functions in epiworldR
+#' Deprecated and removed functions in epiworldR
#' @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`.
+#' @param virus Virus object of class `epiworld_virus`.
#' @name epiworldR-deprecated
-NULL
\ No newline at end of file
+NULL
+
+#' @param n Deprecated.
+#' @export
+#' @rdname epiworldR-deprecated
+add_tool_n <- function(model, tool, n) {
+
+ .Deprecated(new = "add_tool")
+
+ set_distribution_tool(
+ tool,
+ distfun = distribute_tool_randomly(
+ prevalence = n,
+ as_proportion = TRUE
+ )
+ )
+
+ add_tool(model, tool)
+
+}
+
+#' @export
+#' @rdname epiworldR-deprecated
+add_virus_n <- function(model, virus, n) {
+
+ .Deprecated(new = "add_virus")
+
+ set_distribution_virus(
+ virus = virus,
+ 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 163b7fe8..2504e204 100644
--- a/R/global-actions.R
+++ b/R/global-actions.R
@@ -23,6 +23,8 @@
#' # Creating a tool
#' epitool <- tool(
#' name = "Vaccine",
+#' prevalence = 0,
+#' as_proportion = FALSE,
#' susceptibility_reduction = .9,
#' transmission_reduction = .5,
#' recovery_enhancer = .5,
@@ -89,12 +91,10 @@ globalevent_tool <- function(
#' @rdname epiworldR-deprecated
globalaction_tool <- function(...) {
- .Deprecated(
+ .Defunct(
new = "globalevent_tool"
- )
+ )
- globalevent_tool(...)
-
}
#' @export
@@ -132,7 +132,7 @@ globalevent_tool_logit <- function(
#' @rdname epiworldR-deprecated
globalaction_tool_logit <- function(...) {
- .Deprecated(
+ .Defunct(
new = "globalevent_tool_logit"
)
@@ -170,7 +170,8 @@ globalevent_set_params <- function(
#' @rdname epiworldR-deprecated
globalaction_set_params <- function(...) {
- .Deprecated(
+
+ .Defunct(
new = "globalevent_set_params"
)
@@ -232,7 +233,7 @@ globalevent_fun <- function(
#' @rdname epiworldR-deprecated
globalaction_fun <- function(...) {
- .Deprecated(
+ .Defunct(
new = "globalevent_fun"
)
@@ -274,7 +275,7 @@ print.epiworld_globalevent <- function(x, ...) {
add_globalevent <- function(model, action) {
if (length(attr(action, "tool")))
- add_tool_n(model, attr(action, "tool"), 0)
+ add_tool(model, attr(action, "tool"))
invisible(add_globalevent_cpp(model, action))
diff --git a/R/model-methods.R b/R/model-methods.R
index 19991701..61210d0d 100644
--- a/R/model-methods.R
+++ b/R/model-methods.R
@@ -96,7 +96,7 @@ stopifnot_model <- function(model) {
#' 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), proportion = .5)
+#' 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.
diff --git a/R/plot_epi.R b/R/plot_epi.R
index 0922951f..36234322 100644
--- a/R/plot_epi.R
+++ b/R/plot_epi.R
@@ -32,6 +32,7 @@ find_scale <- function(x) {
#' @importFrom graphics legend
plot_epi <- function(x, main = "", counts_scale, ...) UseMethod("plot_epi")
+#' @export
plot_epi.epiworld_model <- function(
x, main = "",
counts_scale,
@@ -47,6 +48,7 @@ plot_epi.epiworld_model <- function(
}
+#' @export
plot_epi.epiworld_hist_virus <- function(
x, main = "",
counts_scale,
@@ -62,6 +64,7 @@ plot_epi.epiworld_hist_virus <- function(
}
+#' @export
plot_epi.epiworld_hist <- function(
x, main = "",
counts_scale,
diff --git a/R/tool.R b/R/tool.R
index 97d0cb8d..3cf628ad 100644
--- a/R/tool.R
+++ b/R/tool.R
@@ -28,6 +28,8 @@
#'
#' epitool <- tool(
#' name = "Vaccine",
+#' prevalence = 0.5,
+#' as_proportion = TRUE,
#' susceptibility_reduction = .9,
#' transmission_reduction = .5,
#' recovery_enhancer = .5,
@@ -38,14 +40,16 @@
#'
#' 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, .5)
+#' 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
-#' add_tool_n(model_sirconn, epitool, 5500)
+#' # 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
@@ -53,6 +57,9 @@
#' 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)
+#' add_tool(model_sirconn, epitool)
#' run(model_sirconn, ndays = 100, seed = 1912) # Run model to view changes
#'
#' @export
@@ -61,6 +68,8 @@
#' @aliases epiworld_tool
tool <- function(
name,
+ prevalence,
+ as_proportion,
susceptibility_reduction,
transmission_reduction,
recovery_enhancer,
@@ -70,6 +79,8 @@ tool <- function(
structure(
tool_cpp(
name,
+ prevalence,
+ as_proportion,
susceptibility_reduction,
transmission_reduction,
recovery_enhancer,
@@ -106,6 +117,16 @@ stopifnot_tfun <- function(tfun) {
}
}
+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): ",
+ paste(class(tool_distfun), collapse = ", ")
+ )
+ }
+}
+
#' @export
#' @details
#' The name of the `epiworld_tool` object can be manipulated with the functions
@@ -132,35 +153,36 @@ get_name_tool <- function(tool) {
#' @export
#' @param tool An object of class `epiworld_tool`
-#' @param proportion In the case of `add_tool`, a proportion, otherwise, an integer.
+#' @param proportion Deprecated.
#' @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) UseMethod("add_tool")
+add_tool <- function(model, tool, proportion) {
+
+ if (!missing(proportion)) {
-#' @export
-add_tool.epiworld_model <- function(model, tool, proportion) {
+ warning(
+ "The 'proportion' argument is deprecated. ",
+ "Use 'set_distribution_tool' instead."
+ )
- stopifnot_tool(tool)
- add_tool_cpp(model, tool, as.double(proportion))
- invisible(model)
+ set_distribution_tool(
+ tool,
+ distribute_tool_randomly(proportion, TRUE)
+ )
-}
+ }
-#' @export
-#' @rdname tool
-#' @returns
-#' - The `add_tool_n` function adds the specified tool to the model of class
-#' [epiworld_model] with specified count n.
-#' @param n A positive integer. Number of agents to initially have the tool.
-add_tool_n <- function(model, tool, n) UseMethod("add_tool_n")
+ UseMethod("add_tool")
+
+}
#' @export
-add_tool_n.epiworld_model <- function(model, tool, n) {
+add_tool.epiworld_model <- function(model, tool, proportion) {
stopifnot_tool(tool)
- add_tool_n_cpp(model, tool, as.integer(n))
+ add_tool_cpp(model, tool)
invisible(model)
}
@@ -204,13 +226,15 @@ rm_tool <- function(model, tool_pos) {
#' # Creating a tool
#' mask_wearing <- tool(
#' name = "Mask",
+#' prevalence = 0.5,
+#' as_proportion = TRUE,
#' susceptibility_reduction = 0.0,
#' transmission_reduction = 0.3, # Only transmission
#' recovery_enhancer = 0.0,
#' death_reduction = 0.0
#' )
#'
-#' add_tool(sir, mask_wearing, .5)
+#' add_tool(sir, mask_wearing)
#'
#' run(sir, ndays = 50, seed = 11)
#' hist_0 <- get_hist_total(sir)
@@ -343,7 +367,7 @@ set_susceptibility_reduction_fun <- function(tool, model, tfun) {
set_transmission_reduction <- function(tool, prob) {
stopifnot_tool(tool)
- set_transmission_reduction_cpp(tool, as.double(prob))
+ invisible(set_transmission_reduction_cpp(tool, as.double(prob)))
}
@@ -353,7 +377,7 @@ set_transmission_reduction_ptr <- function(tool, model, param) {
stopifnot_tool(tool)
stopifnot_model(model)
- set_transmission_reduction_ptr_cpp(tool, model, param)
+ invisible(set_transmission_reduction_ptr_cpp(tool, model, param))
}
@@ -364,7 +388,7 @@ set_transmission_reduction_fun <- function(tool, model, tfun) {
stopifnot_tool(tool)
stopifnot_model(model)
stopifnot_tfun(tfun)
- set_transmission_reduction_fun_cpp(tool, model, tfun)
+ invisible(set_transmission_reduction_fun_cpp(tool, model, tfun))
}
# Recovery enhancer ------------------------------------------------------------
@@ -377,7 +401,7 @@ set_transmission_reduction_fun <- function(tool, model, tfun) {
set_recovery_enhancer <- function(tool, prob) {
stopifnot_tool(tool)
- set_recovery_enhancer_cpp(tool, as.double(prob))
+ invisible(set_recovery_enhancer_cpp(tool, as.double(prob)))
}
@@ -387,7 +411,7 @@ set_recovery_enhancer_ptr <- function(tool, model, param) {
stopifnot_tool(tool)
stopifnot_model(model)
- set_recovery_enhancer_ptr_cpp(tool, model, param)
+ invisible(set_recovery_enhancer_ptr_cpp(tool, model, param))
}
@@ -398,7 +422,7 @@ set_recovery_enhancer_fun <- function(tool, model, tfun) {
stopifnot_tool(tool)
stopifnot_model(model)
stopifnot_tfun(tfun)
- set_recovery_enhancer_fun_cpp(tool, model, tfun)
+ invisible(set_recovery_enhancer_fun_cpp(tool, model, tfun))
}
@@ -412,7 +436,7 @@ set_recovery_enhancer_fun <- function(tool, model, tfun) {
set_death_reduction <- function(tool, prob) {
stopifnot_tool(tool)
- set_death_reduction_cpp(tool, as.double(prob))
+ invisible(set_death_reduction_cpp(tool, as.double(prob)))
}
@@ -422,7 +446,7 @@ set_death_reduction_ptr <- function(tool, model, param) {
stopifnot_tool(tool)
stopifnot_model(model)
- set_death_reduction_ptr_cpp(tool, model, param)
+ invisible(set_death_reduction_ptr_cpp(tool, model, param))
}
@@ -433,7 +457,7 @@ set_death_reduction_fun <- function(tool, model, tfun) {
stopifnot_tool(tool)
stopifnot_model(model)
stopifnot_tfun(tfun)
- set_death_reduction_fun_cpp(tool, model, tfun)
+ invisible(set_death_reduction_fun_cpp(tool, model, tfun))
}
@@ -475,5 +499,67 @@ print.epiworld_agents_tools <- function(x, max_print = 10, ...) {
}
+#' @export
+#' @details
+#' The `set_distribution_tool` function assigns a distribution function to the
+#' specified tool of class [epiworld_tool]. The distribution function can be
+#' created using the functions [distribute_tool_randomly()] and
+#' [distribute_tool_to_set()].
+#' @param distfun An object of class `epiworld_tool_distfun`.
+#' @rdname tool
+set_distribution_tool <- function(tool, distfun) {
+ stopifnot_tool(tool)
+ stopifnot_tool_distfun(distfun)
+ invisible(set_distribution_tool_cpp(tool = tool, distfun = distfun))
+
+}
+
+#' @export
+#' @rdname tool
+#' @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
+#' - The `distribute_tool_randomly` function returns a distribution function of
+#' class `epiworld_tool_distfun`.
+distribute_tool_randomly <- function(
+ prevalence,
+ as_proportion
+) {
+
+ structure(
+ distribute_tool_randomly_cpp(
+ as.double(prevalence),
+ as.logical(as_proportion)
+ ),
+ class = "epiworld_tool_distfun"
+ )
+
+}
+
+#' @export
+#' @rdname tool
+#' @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
+#' will be assigned.
+#' @return
+#' - The `distribute_tool_to_set` function returns a distribution function of
+#' class `epiworld_tool_distfun`.
+distribute_tool_to_set <- function(
+ agents_ids
+) {
+
+ structure(
+ distribute_tool_to_set_cpp(
+ agents_ids
+ ),
+ class = "epiworld_tool_distfun"
+ )
+}
\ No newline at end of file
diff --git a/R/virus.R b/R/virus.R
index e49463db..0bcfd472 100644
--- a/R/virus.R
+++ b/R/virus.R
@@ -30,10 +30,12 @@
#' recovery_rate = 0.99
#' )
#'
-#' delta <- virus("Delta Variant", 0, .5, .2, .01)
+#' 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, .3)
+#' add_virus(mseirconn, delta)
#' set_name_virus(delta, "COVID-19 Strain")
#' get_name_virus(delta)
#'
@@ -41,8 +43,8 @@
#' mseirconn
#'
#' rm_virus(mseirconn, 0) # Removing the first virus from the model object
-#' add_virus_n(mseirconn, delta, 100) # Setting initial count of delta virus
-#' # to n = 100
+#' set_distribution_virus(delta, distribute_virus_randomly(100, as_proportion = FALSE))
+#' add_virus(mseirconn, delta)
#'
#' # Setting parameters for the delta virus manually
#' set_prob_infecting(delta, 0.5)
@@ -54,12 +56,16 @@
#' # 1: Infected
#' # 2: Recovered
#' # 3: Dead
-#' delta2 <- virus("Delta Variant 2", 0, .5, .2, .01)
+#' delta2 <- virus(
+#' "Delta Variant 2", 0, .5, .2, .01, prevalence = 0, as_proportion = TRUE
+#' )
#' virus_set_state(delta2, 1, 2, 3)
#' @export
#' @aliases epiworld_virus
virus <- function(
name,
+ prevalence,
+ as_proportion,
prob_infecting,
recovery_rate = 0.5,
prob_death = 0.0,
@@ -70,6 +76,8 @@ virus <- function(
structure(
virus_cpp(
name,
+ prevalence,
+ as_proportion,
prob_infecting,
recovery_rate,
prob_death,
@@ -106,6 +114,17 @@ stopifnot_vfun <- function(vfun) {
}
}
+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): ",
+ paste(class(virus_distfun), collapse = ", ")
+ )
+ }
+}
+
+
#' @export
#' @details
#' The name of the `epiworld_virus` object can be manipulated with the functions
@@ -135,18 +154,35 @@ get_name_virus <- function(virus) {
#' @rdname virus
#' @param model An object of class `epiworld_model`.
#' @param virus An object of class `epiworld_virus`
-#' @param proportion In the case of `add_virus`, a proportion, otherwise, an integer.
+#' @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) UseMethod("add_virus")
+add_virus <- function(model, virus, proportion) {
+
+ if (!missing(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)
+ )
+
+ }
+
+ UseMethod("add_virus")
+
+}
#' @export
add_virus.epiworld_model <- function(model, virus, proportion) {
stopifnot_virus(virus)
-
- add_virus_cpp(model, virus, proportion)
+ add_virus_cpp(model, virus)
invisible(model)
}
@@ -156,7 +192,7 @@ 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, proportion))
+ invisible(add_virus_cpp(model, virus))
}
@@ -165,7 +201,7 @@ 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, proportion))
+ invisible(add_virus_cpp(model, virus))
}
@@ -173,7 +209,7 @@ add_virus.epiworld_sird <- function(model, virus, proportion) {
add_virus.epiworld_sirconn <- function(model, virus, proportion) {
stopifnot_virus(virus)
- add_virus.epiworld_sir(model, virus, proportion)
+ add_virus.epiworld_sir(model, virus)
}
@@ -181,7 +217,7 @@ add_virus.epiworld_sirconn <- function(model, virus, proportion) {
add_virus.epiworld_sirdconn <- function(model, virus, proportion) {
stopifnot_virus(virus)
- add_virus.epiworld_sird(model, virus, proportion)
+ add_virus.epiworld_sird(model, virus)
}
@@ -191,7 +227,7 @@ 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, proportion))
+ invisible(add_virus_cpp(model, virus))
}
@@ -200,7 +236,7 @@ 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, proportion))
+ invisible(add_virus_cpp(model, virus))
}
@@ -208,7 +244,7 @@ add_virus.epiworld_seird <- function(model, virus, proportion) {
add_virus.epiworld_seirconn <- function(model, virus, proportion) {
stopifnot_virus(virus)
- add_virus.epiworld_seir(model, virus, proportion)
+ add_virus.epiworld_seir(model, virus)
}
@@ -216,58 +252,7 @@ add_virus.epiworld_seirconn <- function(model, virus, proportion) {
add_virus.epiworld_seirdconn <- function(model, virus, proportion) {
stopifnot_virus(virus)
- add_virus.epiworld_seird(model, virus, proportion)
-
-}
-
-#' @export
-#' @rdname virus
-#' @returns
-#' - The `add_virus_n` function does not return a value, but instead adds a
-#' specified number of agents with the virus of choice to the model object
-#' of class [epiworld_model].
-#' @param n A positive integer. Initial count of agents to have the virus.
-add_virus_n <- function(model, virus, n) UseMethod("add_virus_n")
-
-#' @export
-add_virus_n.epiworld_model <- function(model, virus, n) {
-
- stopifnot_virus(virus)
- invisible(add_virus_n_cpp(model, virus, n))
-
-}
-
-#' @export
-add_virus_n.epiworld_sir <- function(model, virus, n) {
-
- stopifnot_virus(virus)
- virus_set_state(virus, init = 1, end = 2, removed = 2)
- invisible(add_virus_n_cpp(model, virus, n))
-
-}
-
-#' @export
-add_virus_n.epiworld_sirconn <- function(model, virus, n) {
-
- stopifnot_virus(virus)
- add_virus_n.epiworld_sir(model, virus, n)
-
-}
-
-#' @export
-add_virus_n.epiworld_seir <- function(model, virus, n) {
-
- stopifnot_virus(virus)
- virus_set_state(virus, init = 1, end = 3, removed = 3)
- invisible(add_virus_n_cpp(model, virus, n))
-
-}
-
-#' @export
-add_virus_n.epiworld_seirconn <- function(model, virus, n) {
-
- stopifnot_virus(virus)
- add_virus_n.epiworld_seir(model, virus, n)
+ add_virus.epiworld_seird(model, virus)
}
@@ -536,5 +521,53 @@ set_incubation_fun <- function(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
+#' @rdname virus
+#' @details The `distribute_virus_randomly` function is a factory function
+#' used to randomly distribute the virus in the model. The prevalence can be set
+#' as a proportion or as a number of agents. The resulting function can then be
+#' passed to `set_distribution_virus`.
+#' @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
+#' - 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
+) {
+
+ structure(
+ distribute_virus_randomly_cpp(
+ as.double(prevalence),
+ as.logical(as_proportion)
+ ),
+ class = "epiworld_virus_distfun"
+ )
+}
+#' @export
+#' @rdname virus
+#' @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 8807652f..ba6d4105 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -20,7 +20,7 @@ knitr::opts_chunk$set(
[](https://github.com/UofUEpiBio/epiworldR/actions/workflows/r.yml)
[](https://cran.r-project.org/package=epiworldR)
[](https://github.com/UofUEpiBio/epiworldR/blob/master/LICENSE.md)
-[](https://codecov.io/gh/UofUEpiBio/epiworldR)
+[](https://app.codecov.io/gh/UofUEpiBio/epiworldR)
This R package is a wrapper of the C++ library [epiworld](https://github.com/UofUEpiBio/epiworld){target="_blank"}. It provides a general framework for modeling disease transmission using [agent-based models](https://en.wikipedia.org/w/index.php?title=Agent-based_model&oldid=1153634802){target="_blank"}. Some of the main features include:
@@ -36,6 +36,16 @@ From the package's 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.
+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")
+
+sprintf("%i. `%s`\n", 1:length(models), models) |>
+ cat()
+```
+
## Installation
You can install the development version of epiworldR from [GitHub](https://github.com/) with:
@@ -52,10 +62,10 @@ install.packages("epiworldR")
# Examples
-This R package includes several popular epidemiological models including
+This R package includes several popular epidemiological models, including
[SIS](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#Variations_on_the_basic_SIR_model){target="_blank"},
[SIR](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SIR_model){target="_blank"}, and
-[SEIR](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SEIR_model){target="_blank"} using either a fully connected graph (similar to a compartmental model) or a user-defined network. Here are some examples:
+[SEIR](https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SEIR_model){target="_blank"} using either a fully connected graph (similar to a compartmental model) or a user-defined network.
## SIR model using a random graph
@@ -81,41 +91,50 @@ sir
Visualizing the outputs
-```{r}
+```{r sir-figures}
summary(sir)
plot(sir)
+plot_incidence(sir)
```
## SEIR model with a fully connected graph
The SEIR model is similar to the SIR model but includes an exposed state. Here, we simulate a population of 10,000 agents with a 0.01 prevalence, a 0.6 transmission rate, a 0.5 recovery rate, and 7 days-incubation period. The population is fully connected, meaning agents can transmit the disease to any other agent:
-```{r}
+```{r seir-conn}
model_seirconn <- ModelSEIRCONN(
name = "COVID-19",
prevalence = 0.01,
n = 10000,
- contact_rate = 4,
+ contact_rate = 10,
incubation_days = 7,
- transmission_rate = 0.6,
- recovery_rate = 0.5
-) |> add_virus(virus("COVID-19", 0.01, 0.6, 0.5, 7), .5)
+ transmission_rate = 0.1,
+ 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
+ ))
set.seed(132)
run(model_seirconn, ndays = 100)
-model_seirconn
+summary(model_seirconn)
```
Computing some key statistics
-```{r}
+```{r seir-conn-figures}
plot(model_seirconn)
repnum <- get_reproductive_number(model_seirconn)
head(plot(repnum))
-plot_incidence(model_seirconn)
head(plot_generation_time(model_seirconn))
```
@@ -132,7 +151,7 @@ in a small-world network. Each agent is connected to eight other agents. One
percent of the population has the virus, with an 80% chance of transmission.
Infected individuals recover at a 0.3 rate:
-```{r}
+```{r logit-model}
# Simulating a population of 100,000 agents
set.seed(2223)
n <- 100000
@@ -212,7 +231,7 @@ nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE)
`epiworldR` supports running multiple simulations using the `run_multiple` function. The following code simulates 50 SIR models with 1000 agents each. Each agent is connected to ten other agents. One percent of the population has the virus, with a 90% chance of transmission. Infected individuals recover at a 0.1 rate. The results are saved in a `data.frame`:
-```{r}
+```{r multiple-example}
model_sir <- ModelSIRCONN(
name = "COVID-19",
prevalence = 0.01,
@@ -255,9 +274,9 @@ citation("epiworldR")
# Existing Alternatives
-Several alternatives to `epiworldR` exist and provide researchers with a range of options, each with its own unique features and strengths, enabling the exploration and analysis of infectious disease dynamics through agent-based modeling. Below is a manually curated table of existing alternatives including ABM [@ABM], abmR [@abmR], cystiSim [@cystiSim], villager [@villager], and RNetLogo [@RNetLogo].
+Several alternatives to `epiworldR` exist and provide researchers with a range of options, each with its own unique features and strengths, enabling the exploration and analysis of infectious disease dynamics through agent-based modeling. Below is a manually curated table of existing alternatives, including ABM [@ABM], abmR [@abmR], cystiSim [@cystiSim], villager [@villager], and RNetLogo [@RNetLogo].
-| Package | Multiple Viruses | Multiple Tools | Multiple Runs | Global Actions | Built-In Epi Models | Dependencies | Activity |
+| Package | Multiple Viruses | Multiple Tools | Multiple Runs | Global Actions | Built-In Epi Models | Dependencies | Activity |
|:--------|:--------|:--------|:--------|:--------|---------|:--------|:--------|
| [**epiworldR**](https://cran.r-project.org/package=epiworldR) | yes | yes | yes | yes | yes | [](https://CRAN.R-project.org/package=epiworldR) | [](https://github.com/UofUEpiBio/epiworldR) |
| [**ABM**](https://cran.r-project.org/package=ABM) | \- | \- | \- | yes | yes | [](https://CRAN.R-project.org/package=ABM) | [](https://github.com/junlingm/ABM) |
@@ -277,4 +296,4 @@ You may want to check out other R packages for agent-based modeling: [`ABM`](htt
## Code of Conduct
-Please note that 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.
+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 210b70a3..270a5067 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ status](https://www.r-pkg.org/badges/version/epiworldR)](https://CRAN.R-project.
downloads](https://cranlogs.r-pkg.org/badges/grand-total/epiworldR)](https://cran.r-project.org/package=epiworldR)
[](https://github.com/UofUEpiBio/epiworldR/blob/master/LICENSE.md)
-[](https://codecov.io/gh/UofUEpiBio/epiworldR)
+[](https://app.codecov.io/gh/UofUEpiBio/epiworldR)
This R package is a wrapper of the C++ library
@@ -43,6 +43,24 @@ From the package’s description:
> Furthermore, epiworldR is ideal for simulation studies featuring large
> populations.
+Current available models:
+
+1. `ModelDiffNet`
+2. `ModelSEIR`
+3. `ModelSEIRCONN`
+4. `ModelSEIRD`
+5. `ModelSEIRDCONN`
+6. `ModelSEIRMixing`
+7. `ModelSIR`
+8. `ModelSIRCONN`
+9. `ModelSIRD`
+10. `ModelSIRDCONN`
+11. `ModelSIRLogit`
+12. `ModelSIRMixing`
+13. `ModelSIS`
+14. `ModelSISD`
+15. `ModelSURV`
+
## Installation
You can install the development version of epiworldR from
@@ -60,16 +78,15 @@ install.packages("epiworldR")
# Examples
-This R package includes several popular epidemiological models including
-SIS, SIR, and SEIR using either a fully connected graph (similar
-to a compartmental model) or a user-defined network. Here are some
-examples:
+to a compartmental model) or a user-defined network.
## SIR model using a random graph
@@ -97,9 +114,6 @@ sir <- ModelSIR(
#> |Running the model...
#> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> | done.
-```
-
-``` r
sir
#> ________________________________________________________________________________
@@ -123,15 +137,15 @@ summary(sir)
#> Number of entities : 0
#> Days (duration) : 50 (of 50)
#> Number of viruses : 1
-#> Last run elapsed t : 165.00ms
-#> Last run speed : 30.13 million agents x day / second
+#> Last run elapsed t : 147.00ms
+#> Last run speed : 33.83 million agents x day / second
#> Rewiring : off
#>
#> Global events:
#> (none)
#>
#> Virus(es):
-#> - COVID-19 (baseline prevalence: 1.00%)
+#> - COVID-19
#>
#> Tool(s):
#> (none)
@@ -149,13 +163,16 @@ summary(sir)
#> - Susceptible 0.91 0.09 0.00
#> - Infected 0.00 0.70 0.30
#> - Recovered 0.00 0.00 1.00
+plot(sir)
```
+
+
``` r
-plot(sir)
+plot_incidence(sir)
```
-
+
## SEIR model with a fully connected graph
@@ -170,11 +187,20 @@ model_seirconn <- ModelSEIRCONN(
name = "COVID-19",
prevalence = 0.01,
n = 10000,
- contact_rate = 4,
+ contact_rate = 10,
incubation_days = 7,
- transmission_rate = 0.6,
- recovery_rate = 0.5
-) |> add_virus(virus("COVID-19", 0.01, 0.6, 0.5, 7), .5)
+ transmission_rate = 0.1,
+ 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
+ ))
set.seed(132)
run(model_seirconn, ndays = 100)
@@ -182,15 +208,48 @@ run(model_seirconn, ndays = 100)
#> Running the model...
#> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> done.
-```
-
-``` r
-model_seirconn
+summary(model_seirconn)
+#> ________________________________________________________________________________
#> ________________________________________________________________________________
-#> Susceptible-Exposed-Infected-Removed (SEIR) (connected)
-#> It features 10000 agents, 2 virus(es), and 0 tool(s).
-#> The model has 4 states.
-#> The final distribution is: 634 Susceptible, 5 Exposed, 0 Infected, and 9361 Recovered.
+#> 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 : 80.00ms
+#> Last run speed : 12.38 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
+#> - (1) Exposed : 200 -> 0
+#> - (2) Infected : 0 -> 1
+#> - (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
+#> - Recovered 0.00 0.00 0.00 1.00
```
Computing some key statistics
@@ -199,7 +258,7 @@ Computing some key statistics
plot(model_seirconn)
```
-
+
``` r
@@ -208,36 +267,27 @@ repnum <- get_reproductive_number(model_seirconn)
head(plot(repnum))
```
-
+
- #> virus_id virus date avg n sd lb ub
- #> 1 0 COVID-19 0 2.762500 80 2.082135 1.0 7.025
- #> 2 0 COVID-19 2 3.250000 24 2.862805 0.0 9.850
- #> 3 0 COVID-19 3 3.294118 17 2.663755 0.4 9.400
- #> 4 0 COVID-19 4 2.666667 18 2.351470 0.0 7.875
- #> 5 0 COVID-19 5 1.878788 33 1.745666 0.0 5.800
- #> 6 0 COVID-19 6 1.794118 34 1.533058 0.0 4.350
+ #> 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
-``` r
-
-plot_incidence(model_seirconn)
-```
-
-
-
-``` r
-head(plot_generation_time(model_seirconn))
-```
+ head(plot_generation_time(model_seirconn))
-
+
- #> date avg n sd ci_lower ci_upper virus virus_id
- #> 1 2 8.400000 20 6.227274 2.475 24.150 COVID-19 0
- #> 2 3 8.750000 16 6.547264 2.375 21.750 COVID-19 0
- #> 3 4 7.625000 16 5.302515 2.375 19.000 COVID-19 0
- #> 4 5 5.888889 27 3.178453 2.000 12.700 COVID-19 0
- #> 5 6 10.148148 27 5.586410 2.000 21.400 COVID-19 0
- #> 6 7 8.458333 24 6.064717 2.000 20.825 COVID-19 0
+ #> 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
## SIR Logit
@@ -288,13 +338,10 @@ run(model_logit, 50)
#> |Running the model...
#> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> | done.
-```
-
-``` r
plot(model_logit)
```
-
+
``` r
@@ -307,9 +354,6 @@ rn <- get_reproductive_number(model_logit)
) |> prop.table())[,2]
#> 0 1
#> 0.12984 0.14201
-```
-
-``` r
# Looking into the agents
get_agents(model_logit)
@@ -351,9 +395,6 @@ sir <- ModelSIR(
#> |Running the model...
#> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> | done.
-```
-
-``` r
# Transmission network
net <- get_transmissions(sir)
@@ -362,9 +403,6 @@ net <- get_transmissions(sir)
library(epiworldR)
library(netplot)
#> Loading required package: grid
-```
-
-``` r
x <- igraph::graph_from_edgelist(
as.matrix(net[,2:3]) + 1
)
@@ -403,9 +441,6 @@ run_multiple(model_sir, ndays = 100, nsims = 50, saver = saver, nthread = 2)
#> _________________________________________________________________________
#> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> done.
-```
-
-``` r
# Retrieving the results
ans <- run_multiple_get_results(model_sir)
@@ -418,9 +453,6 @@ head(ans$total_hist)
#> 4 1 1 1 Susceptible 977
#> 5 1 1 1 Infected 22
#> 6 1 1 1 Recovered 1
-```
-
-``` r
head(ans$reproductive)
#> sim_num virus_id virus source source_exposure_date rt
#> 1 1 0 COVID-19 976 9 0
@@ -429,14 +461,11 @@ head(ans$reproductive)
#> 4 1 0 COVID-19 314 9 0
#> 5 1 0 COVID-19 41 9 0
#> 6 1 0 COVID-19 32 9 0
-```
-
-``` r
plot(ans$reproductive)
```
-
+
# Tutorials
@@ -461,7 +490,7 @@ citation("epiworldR")
#> And the actual R package:
#>
#> Meyer D, Vega Yon G (2024). _epiworldR: Fast Agent-Based Epi Models_.
-#> R package version 0.2-1, .
+#> R package version 0.3-2, .
#>
#> To see these entries in BibTeX format, use 'print(,
#> bibtex=TRUE)', 'toBibtex(.)', or set
@@ -474,17 +503,17 @@ Several alternatives to `epiworldR` exist and provide researchers with a
range of options, each with its own unique features and strengths,
enabling the exploration and analysis of infectious disease dynamics
through agent-based modeling. Below is a manually curated table of
-existing alternatives including ABM \[@ABM\], abmR \[@abmR\], cystiSim
+existing alternatives, including ABM \[@ABM\], abmR \[@abmR\], cystiSim
\[@cystiSim\], villager \[@villager\], and RNetLogo \[@RNetLogo\].
-| Package | Multiple Viruses | Multiple Tools | Multiple Runs | Global Actions | Built-In Epi Models | Dependencies | Activity |
-|:--------------------------------------------------------------|:-----------------|:---------------|:--------------|:---------------|---------------------|:---------------------------------------------------------------------------------------------------------|:-----------------------------------------------------------------------------------------------------------------------|
-| [**epiworldR**](https://cran.r-project.org/package=epiworldR) | yes | yes | yes | yes | yes | [](https://CRAN.R-project.org/package=epiworldR) | [](https://github.com/UofUEpiBio/epiworldR) |
-| [**ABM**](https://cran.r-project.org/package=ABM) | \- | \- | \- | yes | yes | [](https://CRAN.R-project.org/package=ABM) | [](https://github.com/junlingm/ABM) |
-| [**abmR**](https://cran.r-project.org/package=abmR) | \- | \- | yes | \- | \- | [](https://CRAN.R-project.org/package=abmR) | [](https://github.com/bgoch5/abmR) |
-| [**cystiSim**](https://cran.r-project.org/package=cystiSim) | \- | yes | yes | \- | \- | [](https://CRAN.R-project.org/package=cystiSim) | [](https://github.com/brechtdv/cystiSim) |
-| [**villager**](https://cran.r-project.org/package=villager) | \- | \- | \- | yes | \- | [](https://CRAN.R-project.org/package=villager) | [](https://github.com/zizroc/villager) |
-| [**RNetLogo**](https://cran.r-project.org/package=RNetLogo) | \- | yes | yes | yes | \- | [](https://CRAN.R-project.org/package=RNetLogo) | [](https://github.com/cran/RNetLogo) |
+| Package | Multiple Viruses | Multiple Tools | Multiple Runs | Global Actions | Built-In Epi Models | Dependencies | Activity |
+|:---|:---|:---|:---|:---|----|:---|:---|
+| [**epiworldR**](https://cran.r-project.org/package=epiworldR) | yes | yes | yes | yes | yes | [](https://CRAN.R-project.org/package=epiworldR) | [](https://github.com/UofUEpiBio/epiworldR) |
+| [**ABM**](https://cran.r-project.org/package=ABM) | \- | \- | \- | yes | yes | [](https://CRAN.R-project.org/package=ABM) | [](https://github.com/junlingm/ABM) |
+| [**abmR**](https://cran.r-project.org/package=abmR) | \- | \- | yes | \- | \- | [](https://CRAN.R-project.org/package=abmR) | [](https://github.com/bgoch5/abmR) |
+| [**cystiSim**](https://cran.r-project.org/package=cystiSim) | \- | yes | yes | \- | \- | [](https://CRAN.R-project.org/package=cystiSim) | [](https://github.com/brechtdv/cystiSim) |
+| [**villager**](https://cran.r-project.org/package=villager) | \- | \- | \- | yes | \- | [](https://CRAN.R-project.org/package=villager) | [](https://github.com/zizroc/villager) |
+| [**RNetLogo**](https://cran.r-project.org/package=RNetLogo) | \- | yes | yes | yes | \- | [](https://CRAN.R-project.org/package=RNetLogo) | [](https://github.com/cran/RNetLogo) |
# Other ABM R packages
@@ -502,7 +531,6 @@ target="_blank">RNetLogo.
## Code of Conduct
-Please note that the epiworldR project is released with a [Contributor
-Code of
+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/docker/Makefile b/docker/Makefile
index dffd5d55..b161464e 100644
--- a/docker/Makefile
+++ b/docker/Makefile
@@ -1,2 +1,7 @@
all:
- docker build -t uofuepibio/epiworldr:debug -f Dockerfile .
\ No newline at end of file
+ docker build -t uofuepibio/epiworldr:debug -f Dockerfile .
+
+podman:
+ podman build -t uofuepibio/epiworldr:debug -f Dockerfile .
+
+.PHONY: all podman
\ No newline at end of file
diff --git a/inst/include/epiworld/agent-bones.hpp b/inst/include/epiworld/agent-bones.hpp
index 1e83c069..c4e75dea 100644
--- a/inst/include/epiworld/agent-bones.hpp
+++ b/inst/include/epiworld/agent-bones.hpp
@@ -102,9 +102,9 @@ class Agent {
std::vector< ToolPtr > tools;
epiworld_fast_uint n_tools = 0u;
- std::vector< Agent * > sampled_agents;
+ std::vector< Agent * > sampled_agents = {};
size_t sampled_agents_n = 0u;
- std::vector< size_t > sampled_agents_left;
+ std::vector< size_t > sampled_agents_left = {};
size_t sampled_agents_left_n = 0u;
int date_last_build_sample = -99;
@@ -218,6 +218,7 @@ class Agent {
int get_id() const; ///< Id of the individual
VirusPtr & get_virus();
+ const VirusPtr & get_virus() const;
ToolPtr & get_tool(int i);
Tools get_tools();
@@ -263,6 +264,8 @@ class Agent {
bool has_virus(epiworld_fast_uint t) const;
bool has_virus(std::string name) const;
bool has_virus(const Virus & v) const;
+ bool has_entity(epiworld_fast_uint t) const;
+ bool has_entity(std::string name) const;
void print(Model * model, bool compressed = false) const;
diff --git a/inst/include/epiworld/agent-events-meat.hpp b/inst/include/epiworld/agent-events-meat.hpp
index 8da043fb..cfa9c275 100644
--- a/inst/include/epiworld/agent-events-meat.hpp
+++ b/inst/include/epiworld/agent-events-meat.hpp
@@ -192,7 +192,7 @@ inline void default_add_entity(Event & a, Model *)
if (p->get_n_entities() > e->size()) // Slower search through the agent
{
for (size_t i = 0u; i < e->size(); ++i)
- if(e->operator[](i)->get_id() == p->get_id())
+ if(static_cast(e->operator[](i)) == p->get_id())
throw std::logic_error("An entity cannot be reassigned to an agent.");
}
else // Slower search through the entity
@@ -253,11 +253,15 @@ inline void default_rm_entity(Event & a, Model * m)
// When we move the end entity to the new location, the
// moved entity needs to reflect the change, i.e., where the
// entity will now be located in the agent
- size_t agent_location_in_last_entity = p->entities_locations[p->n_entities];
- Entity * last_entity = &m->get_entities()[p->entities[p->n_entities]]; ///< Last entity of the agent
+ size_t agent_location_in_last_entity =
+ p->entities_locations[p->n_entities];
+
+ Entity * last_entity =
+ &m->get_entity(p->entities[p->n_entities]); ///< Last entity of the agent
// The end entity will be located where the removed was
- last_entity->agents_location[agent_location_in_last_entity] = idx_entity_in_agent;
+ last_entity->agents_location[agent_location_in_last_entity] =
+ idx_entity_in_agent;
// We now make the swap
std::swap(
@@ -274,10 +278,13 @@ inline void default_rm_entity(Event & a, Model * m)
// moved agent needs to reflect the change, i.e., where the
// agent will now be located in the entity
size_t entity_location_in_last_agent = e->agents_location[e->n_agents];
- Agent * last_agent = &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity
+
+ Agent * last_agent =
+ &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity
// The end entity will be located where the removed was
- last_agent->entities_locations[entity_location_in_last_agent] = idx_agent_in_entity;
+ last_agent->entities_locations[entity_location_in_last_agent] =
+ idx_agent_in_entity;
// We now make the swap
std::swap(
diff --git a/inst/include/epiworld/agent-meat.hpp b/inst/include/epiworld/agent-meat.hpp
index a5006b5d..27d4785e 100644
--- a/inst/include/epiworld/agent-meat.hpp
+++ b/inst/include/epiworld/agent-meat.hpp
@@ -242,7 +242,7 @@ inline void Agent::add_entity(
-1, -1
);
- // default_add_entity(a, model); /* passing model makes nothing */
+ default_add_entity(a, model); /* passing model makes nothing */
}
@@ -332,12 +332,19 @@ inline void Agent::rm_entity(
"There is entity to remove here!"
);
- CHECK_COALESCE_(state_new, model->entities[entity_idx].state_post, state);
- CHECK_COALESCE_(queue, model->entities[entity_idx].queue_post, Queue::NoOne);
+ CHECK_COALESCE_(state_new, model->get_entity(entity_idx).state_post, state);
+ CHECK_COALESCE_(queue, model->get_entity(entity_idx).queue_post, Queue::NoOne);
model->events_add(
- this, nullptr, nullptr, model->entities[entity_idx], state_new, queue,
- default_rm_entity, entities_locations[entity_idx], entity_idx
+ this,
+ nullptr,
+ nullptr,
+ &model->get_entity(entity_idx),
+ state_new,
+ queue,
+ default_rm_entity,
+ entities_locations[entity_idx],
+ entity_idx
);
}
@@ -354,22 +361,35 @@ inline void Agent::rm_entity(
int entity_idx = -1;
for (size_t i = 0u; i < n_entities; ++i)
{
- if (entities[i] == entity->get_id())
+ if (static_cast(entities[i]) == entity.get_id())
+ {
entity_idx = i;
+ break;
+ }
}
if (entity_idx == -1)
throw std::logic_error(
- "The agent " + std::to_string(id) + " is not associated with entity \"" +
- entity.get_name() + "\"."
+ std::string("The agent ") +
+ std::to_string(id) +
+ std::string(" is not associated with entity \"") +
+ entity.get_name() +
+ std::string("\".")
);
CHECK_COALESCE_(state_new, entity.state_post, state);
CHECK_COALESCE_(queue, entity.queue_post, Queue::NoOne);
model->events_add(
- this, nullptr, nullptr, entities[entity_idx], state_new, queue,
- default_rm_entity, entities_locations[entity_idx], entity_idx
+ this,
+ nullptr,
+ nullptr,
+ &model->entities[entity.get_id()],
+ state_new,
+ queue,
+ default_rm_entity,
+ entities_locations[entity_idx],
+ entity_idx
);
}
@@ -435,6 +455,11 @@ inline VirusPtr & Agent::get_virus() {
return virus;
}
+template
+inline const VirusPtr & Agent::get_virus() const {
+ return virus;
+}
+
template
inline Tools Agent::get_tools() {
@@ -610,6 +635,10 @@ inline void Agent::reset()
this->tools.clear();
n_tools = 0u;
+ this->entities.clear();
+ this->entities_locations.clear();
+ this->n_entities = 0u;
+
this->state = 0u;
this->state_prev = 0u;
@@ -677,6 +706,30 @@ inline bool Agent::has_virus(const Virus & virus) const
}
+template
+inline bool Agent::has_entity(epiworld_fast_uint t) const
+{
+
+ for (auto & entity : entities)
+ if (entity == t)
+ return true;
+
+ return false;
+
+}
+
+template
+inline bool Agent::has_entity(std::string name) const
+{
+
+ for (auto & entity : entities)
+ if (model->get_entity(entity).get_name() == name)
+ return true;
+
+ return false;
+
+}
+
template
inline void Agent::print(
Model * model,
@@ -788,19 +841,25 @@ inline const Entities_const Agent::get_entities() const
template
inline const Entity & Agent::get_entity(size_t i) const
{
+ if (n_entities == 0)
+ throw std::range_error("Agent id " + std::to_string(id) + " has no entities.");
+
if (i >= n_entities)
throw std::range_error("Trying to get to an agent's entity outside of the range.");
- return model->entities[entities[i]];
+ return model->get_entity(entities[i]);
}
template
inline Entity & Agent::get_entity(size_t i)
{
+ if (n_entities == 0)
+ throw std::range_error("Agent id " + std::to_string(id) + " has no entities.");
+
if (i >= n_entities)
throw std::range_error("Trying to get to an agent's entity outside of the range.");
- return model->entities[entities[i]];
+ return model->get_entity(entities[i]);
}
template
diff --git a/inst/include/epiworld/agentssample-bones.hpp b/inst/include/epiworld/agentssample-bones.hpp
index c8462c58..e0d1285a 100644
--- a/inst/include/epiworld/agentssample-bones.hpp
+++ b/inst/include/epiworld/agentssample-bones.hpp
@@ -30,10 +30,10 @@ class AgentsSample {
size_t sample_size = 0u;
- std::vector< Agent* > * agents = nullptr; ///< Pointer to sample of agents
+ std::vector< Agent* >* agents = nullptr; ///< Pointer to sample of agents
size_t * agents_n = nullptr; ///< Size of sample of agents
- std::vector< size_t > * agents_left = nullptr; ///< Pointer to agents left (iota)
+ std::vector< size_t >* agents_left = nullptr; ///< Pointer to agents left (iota)
size_t * agents_left_n = nullptr; ///< Size of agents left
Model * model = nullptr; ///< Extracts runif() and (if the case) population.
@@ -49,7 +49,7 @@ class AgentsSample {
public:
// Not available (for now)
- AgentsSample() = delete; ///< Default constructor
+ AgentsSample() = delete; ///< Default constructor
AgentsSample(const AgentsSample & a) = delete; ///< Copy constructor
AgentsSample(AgentsSample && a) = delete; ///< Move constructor
@@ -60,13 +60,17 @@ class AgentsSample {
);
AgentsSample(
- Model * model, Entity & entity_, size_t n,
+ Model * model,
+ Entity & entity_,
+ size_t n,
std::vector< size_t > states_ = {},
bool truncate = false
);
AgentsSample(
- Model * model, Agent & agent_, size_t n,
+ Model * model,
+ Agent & agent_,
+ size_t n,
std::vector< size_t > states_ = {},
bool truncate = false
);
@@ -187,9 +191,6 @@ inline AgentsSample::AgentsSample(
agents = &agent_.sampled_agents;
agents_n = &agent_.sampled_agents_n;
- agents_left = &agent_.sampled_agents_left;
- agents_left_n = &agent_.sampled_agents_left_n;
-
// Computing the cumulative sum of counts across entities
size_t agents_in_entities = 0;
Entities entities_a = agent->get_entities();
@@ -228,31 +229,37 @@ inline AgentsSample::AgentsSample(
agents->resize(n);
size_t i_obs = 0u;
- for (size_t i = 0u; i < agents_in_entities; ++i)
+ for (size_t i = 0u; i < sample_size; ++i)
{
+
+ // Sampling a single agent from the set of entities
int jth = std::floor(model->runif() * agents_in_entities);
for (size_t e = 0u; e < cum_agents_count.size(); ++e)
{
+
// Are we in the limit?
if (jth <= cum_agents_count[e])
{
size_t agent_idx = 0u;
if (e == 0) // From the first group
- agent_idx = entities_a[e][jth];
+ agent_idx = entities_a[e][jth]->get_id();
else
- agent_idx = entities_a[e][jth - cum_agents_count[e - 1]];
+ agent_idx = entities_a[e][jth - cum_agents_count[e - 1]]->get_id();
- // Getting the state
- size_t state = model->population[agent_idx].get_state();
// Checking if states was specified
if (states.size())
{
+
+ // Getting the state
+ size_t state = model->population[agent_idx].get_state();
+
if (std::find(states.begin(), states.end(), state) != states.end())
continue;
+
}
- agents->operator[](i_obs++) = agent_idx;
+ agents->operator[](i_obs++) = &(model->population[agent_idx]);
break;
}
diff --git a/inst/include/epiworld/config.hpp b/inst/include/epiworld/config.hpp
index 852243d9..d9251061 100644
--- a/inst/include/epiworld/config.hpp
+++ b/inst/include/epiworld/config.hpp
@@ -75,7 +75,7 @@ using MixerFun = std::function*,VirusPtr,Model
template
using MutFun = std::function*,Virus&,Model*)>;
-template
+template
using PostRecoveryFun = std::function*,Virus&,Model*)>;
template
@@ -90,25 +90,25 @@ using GlobalFun = std::function*)>;
template
struct Event;
-template
-using ActionFun = std::function&,Model*)>;
+template
+using EventFun = std::function&,Model*)>;
/**
* @brief Decides how to distribute viruses at initialization
*/
-template
+template
using VirusToAgentFun = std::function&,Model*)>;
/**
* @brief Decides how to distribute tools at initialization
*/
-template
+template
using ToolToAgentFun = std::function&,Model*)>;
/**
* @brief Decides how to distribute entities at initialization
*/
-template
+template
using EntityToAgentFun = std::function&,Model*)>;
/**
@@ -116,7 +116,7 @@ using EntityToAgentFun = std::function&,Model*)>;
*
* @tparam TSeq
*/
-template
+template
struct Event {
Agent * agent;
VirusPtr virus;
@@ -124,7 +124,7 @@ struct Event {
Entity * entity;
epiworld_fast_int new_state;
epiworld_fast_int queue;
- ActionFun call;
+ EventFun call;
int idx_agent;
int idx_object;
public:
@@ -151,7 +151,7 @@ struct Event {
Entity * entity_,
epiworld_fast_int new_state_,
epiworld_fast_int queue_,
- ActionFun call_,
+ EventFun call_,
int idx_agent_,
int idx_object_
) : agent(agent_), virus(virus_), tool(tool_), entity(entity_),
diff --git a/inst/include/epiworld/entity-bones.hpp b/inst/include/epiworld/entity-bones.hpp
index 6e2ac160..b36190f5 100644
--- a/inst/include/epiworld/entity-bones.hpp
+++ b/inst/include/epiworld/entity-bones.hpp
@@ -24,8 +24,6 @@ class Entity {
friend void default_add_entity(Event & a, Model * m);
friend void default_rm_entity(Event & a, Model * m);
private:
-
- Model * model;
int id = -1;
std::vector< size_t > agents; ///< Vector of agents
@@ -49,7 +47,7 @@ class Entity {
///@}
int max_capacity = -1;
- std::string entity_name = "Unknown entity";
+ std::string entity_name = "Unnamed entity";
std::vector< epiworld_double > location = {0.0}; ///< An arbitrary vector for location
@@ -59,18 +57,30 @@ class Entity {
epiworld_fast_int queue_init = 0; ///< Change of state when added to agent.
epiworld_fast_int queue_post = 0; ///< Change of state when removed from agent.
+ EntityToAgentFun dist_fun = nullptr;
+
public:
- // Entity() = delete;
- // Entity(Entity & e) = delete;
- // Entity(const Entity & e);
- // Entity(Entity && e);
- Entity(std::string name) : entity_name(name) {};
- // Entity & operator=(const Entity & e);
+ /**
+ * @brief Constructs an Entity object.
+ *
+ * This constructor initializes an Entity object with the specified parameters.
+ *
+ * @param name The name of the entity.
+ * @param fun A function pointer to a function that maps the entity to an agent.
+ */
+ Entity(
+ std::string name,
+ EntityToAgentFun fun = nullptr
+ ) :
+ entity_name(name),
+ dist_fun(fun)
+ {};
+
void add_agent(Agent & p, Model * model);
void add_agent(Agent * p, Model * model);
- void rm_agent(size_t idx);
+ void rm_agent(size_t idx, Model * model);
size_t size() const noexcept;
void set_location(std::vector< epiworld_double > loc);
std::vector< epiworld_double > & get_location();
@@ -81,7 +91,7 @@ class Entity {
typename std::vector< Agent * >::const_iterator begin() const;
typename std::vector< Agent * >::const_iterator end() const;
- Agent * operator[](size_t i);
+ size_t operator[](size_t i);
int get_id() const noexcept;
const std::string & get_name() const noexcept;
@@ -96,6 +106,20 @@ class Entity {
bool operator==(const Entity & other) const;
bool operator!=(const Entity & other) const {return !operator==(other);};
+ /**
+ * @name Entity distribution
+ *
+ * @details These functions are used for distributing agents among entities.
+ * The idea is to have a flexible way of distributing agents among entities.
+
+ */
+ void distribute(Model * model);
+
+ std::vector< size_t > & get_agents();
+
+ void print() const;
+ void set_distribution(EntityToAgentFun fun);
+
};
diff --git a/inst/include/epiworld/entity-distribute-meat.hpp b/inst/include/epiworld/entity-distribute-meat.hpp
new file mode 100644
index 00000000..d589a368
--- /dev/null
+++ b/inst/include/epiworld/entity-distribute-meat.hpp
@@ -0,0 +1,157 @@
+#ifndef EPIWORLD_ENTITY_DISTRIBUTE_MEAT_HPP
+#define EPIWORLD_ENTITY_DISTRIBUTE_MEAT_HPP
+
+
+template
+/**
+ * Distributes an entity to unassigned agents in the model.
+ *
+ * @param prevalence The proportion of agents to distribute the entity to.
+ * @param as_proportion Flag indicating whether the prevalence is a proportion
+ * @param to_unassigned Flag indicating whether to distribute the entity only
+ * to unassigned agents.
+ * @return An EntityToAgentFun object that distributes the entity to unassigned
+ * agents.
+ */
+inline EntityToAgentFun distribute_entity_randomly(
+ epiworld_double prevalence,
+ bool as_proportion,
+ bool to_unassigned
+)
+{
+
+ return [prevalence, as_proportion, to_unassigned](
+ Entity & e, Model * m
+ ) -> void {
+
+
+ // Preparing the sampling space
+ std::vector< size_t > idx;
+ if (to_unassigned)
+ {
+ for (const auto & a: m->get_agents())
+ if (a.get_n_entities() == 0)
+ idx.push_back(a.get_id());
+ }
+ else
+ {
+
+ for (const auto & a: m->get_agents())
+ idx.push_back(a.get_id());
+
+ }
+
+ size_t n = idx.size();
+
+ // Figuring out how many to sample
+ int n_to_sample;
+ if (as_proportion)
+ {
+ n_to_sample = static_cast(std::floor(prevalence * n));
+ if (n_to_sample > static_cast(n))
+ --n_to_sample;
+
+ } else
+ {
+ n_to_sample = static_cast(prevalence);
+ if (n_to_sample > static_cast(n))
+ throw std::range_error("There are only " + std::to_string(n) +
+ " individuals in the population. Cannot add the entity to " +
+ std::to_string(n_to_sample));
+ }
+
+ int n_left = n;
+ for (int i = 0; i < n_to_sample; ++i)
+ {
+ int loc = static_cast(
+ floor(m->runif() * n_left--)
+ );
+
+ // Correcting for possible overflow
+ if ((loc > 0) && (loc >= n_left))
+ loc = n_left - 1;
+
+ m->get_agent(idx[loc]).add_entity(e, m);
+
+ std::swap(idx[loc], idx[n_left]);
+
+ }
+
+ };
+
+}
+
+template
+/**
+ * Distributes an entity to a range of agents.
+ *
+ * @param from The starting index of the range.
+ * @param to The ending index of the range.
+ * @param to_unassigned Flag indicating whether to distribute the entity only
+ * to unassigned agents.
+ * @return A lambda function that distributes the entity to the specified range
+ * of agents.
+ */
+inline EntityToAgentFun distribute_entity_to_range(
+ int from,
+ int to,
+ bool to_unassigned = false
+ ) {
+
+ if (to_unassigned)
+ {
+
+ return [from, to](Entity & e, Model * m) -> void {
+
+ auto & agents = m->get_agents();
+ for (size_t i = from; i < to; ++i)
+ {
+ if (agents[i].get_n_entities() == 0)
+ e.add_agent(&agents[i], m);
+ else
+ throw std::logic_error(
+ "Agent " + std::to_string(i) + " already has an entity."
+ );
+ }
+
+ return;
+
+ };
+
+ }
+ else
+ {
+
+ return [from, to](Entity & e, Model * m) -> void {
+
+ auto & agents = m->get_agents();
+ for (size_t i = from; i < to; ++i)
+ {
+ e.add_agent(&agents[i], m);
+ }
+
+ return;
+
+ };
+
+ }
+}
+
+
+template
+inline EntityToAgentFun distribute_entity_to_set(
+ std::vector< size_t > & idx
+ ) {
+
+ return [idx](Entity & e, Model * m) -> void {
+
+ for (const auto & i: idx)
+ {
+ e.add_agent(&m->get_agent(i), m);
+ }
+
+ };
+
+}
+
+#endif
\ No newline at end of file
diff --git a/inst/include/epiworld/entity-meat.hpp b/inst/include/epiworld/entity-meat.hpp
index 4c5a30d4..7c1200f6 100644
--- a/inst/include/epiworld/entity-meat.hpp
+++ b/inst/include/epiworld/entity-meat.hpp
@@ -23,7 +23,7 @@ inline void Entity::add_agent(
}
template
-inline void Entity::rm_agent(size_t idx)
+inline void Entity::rm_agent(size_t idx, Model * model)
{
if (idx >= n_agents)
throw std::out_of_range(
@@ -31,7 +31,7 @@ inline void Entity::rm_agent(size_t idx)
" out of " + std::to_string(n_agents)
);
- model->population[idx].rm_entity(*this);
+ model->get_agents()[agents[idx]].rm_entity(*this, model);
return;
}
@@ -89,12 +89,15 @@ inline typename std::vector< Agent * >::const_iterator Entity::end()
}
template
-inline Agent * Entity::operator[](size_t i)
+size_t Entity::operator[](size_t i)
{
if (n_agents <= i)
- throw std::logic_error("There are not that many agents in this entity.");
+ throw std::logic_error(
+ "There are not that many agents in this entity. " +
+ std::to_string(n_agents) + " <= " + std::to_string(i)
+ );
- return &model->get_agents()[i];
+ return i;
}
template
@@ -164,6 +167,13 @@ inline void Entity::reset()
sampled_agents_n = 0u;
sampled_agents_left.clear();
sampled_agents_left_n = 0u;
+
+ this->agents.clear();
+ this->n_agents = 0u;
+ this->agents_location.clear();
+
+ return;
+
}
template
@@ -216,4 +226,41 @@ inline bool Entity::operator==(const Entity & other) const
}
+template
+inline void Entity::distribute(Model * model)
+{
+
+ if (dist_fun)
+ {
+
+ dist_fun(*this, model);
+
+ }
+
+}
+
+template
+inline std::vector< size_t > & Entity::get_agents()
+{
+ return agents;
+}
+
+template
+inline void Entity::print() const
+{
+
+ printf_epiworld(
+ "Entity '%s' (id %i) with %i agents.\n",
+ this->entity_name.c_str(),
+ static_cast(id),
+ static_cast(n_agents)
+ );
+}
+
+template
+inline void Entity::set_distribution(EntityToAgentFun fun)
+{
+ dist_fun = fun;
+}
+
#endif
\ No newline at end of file
diff --git a/inst/include/epiworld/epiworld-macros.hpp b/inst/include/epiworld/epiworld-macros.hpp
index 0fdf87b3..149a4961 100644
--- a/inst/include/epiworld/epiworld-macros.hpp
+++ b/inst/include/epiworld/epiworld-macros.hpp
@@ -101,4 +101,11 @@
[](epiworld::Model* m) -> void
+#define EPI_NEW_ENTITYTOAGENTFUN(funname,tseq) inline void \
+ (funname)(epiworld::Entity & e, epiworld::Model * m)
+
+#define EPI_NEW_ENTITYTOAGENTFUN_LAMBDA(funname,tseq) \
+ epiworld::EntityToAgentFun funname = \
+ [](epiworld::Entity & e, epiworld::Model * m) -> void
+
#endif
diff --git a/inst/include/epiworld/epiworld.hpp b/inst/include/epiworld/epiworld.hpp
index 5d7d7bdc..7e8f6624 100644
--- a/inst/include/epiworld/epiworld.hpp
+++ b/inst/include/epiworld/epiworld.hpp
@@ -16,11 +16,10 @@
#ifndef EPIWORLD_HPP
#define EPIWORLD_HPP
-
/* Versioning */
#define EPIWORLD_VERSION_MAJOR 0
-#define EPIWORLD_VERSION_MINOR 1
-#define EPIWORLD_VERSION_PATCH 1
+#define EPIWORLD_VERSION_MINOR 3
+#define EPIWORLD_VERSION_PATCH 2
static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR;
static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR;
@@ -61,14 +60,17 @@ namespace epiworld {
#include "viruses-bones.hpp"
#include "virus-bones.hpp"
+ #include "virus-distribute-meat.hpp"
#include "virus-meat.hpp"
#include "tools-bones.hpp"
#include "tool-bones.hpp"
+ #include "tool-distribute-meat.hpp"
#include "tool-meat.hpp"
#include "entity-bones.hpp"
+ #include "entity-distribute-meat.hpp"
#include "entity-meat.hpp"
#include "entities-bones.hpp"
@@ -79,6 +81,9 @@ namespace epiworld {
#include "agentssample-bones.hpp"
+ #include "groupsampler-bones.hpp"
+ #include "groupsampler-meat.hpp"
+
#include "models/models.hpp"
}
diff --git a/inst/include/epiworld/groupsampler-bones.hpp b/inst/include/epiworld/groupsampler-bones.hpp
new file mode 100644
index 00000000..88ec690a
--- /dev/null
+++ b/inst/include/epiworld/groupsampler-bones.hpp
@@ -0,0 +1,59 @@
+#ifndef GROUPSAMPLER_BONES_HPP
+#define GROUPSAMPLER_BONES_HPP
+
+/**
+ * @brief Weighted sampling of groups
+ */
+template
+class GroupSampler {
+
+private:
+
+ std::vector< double > contact_matrix; ///< Contact matrix between groups
+ std::vector< size_t > group_sizes; ///< Sizes of the groups
+ std::vector< double > cumulate; ///< Cumulative sum of the contact matrix (row-major for faster access)
+
+ /**
+ * @brief Get the index of the contact matrix
+ *
+ * The matrix is a vector stored in column-major order.
+ *
+ * @param i Index of the row
+ * @param j Index of the column
+ * @return Index of the contact matrix
+ */
+ inline int idx(const int i, const int j, bool rowmajor = false) const
+ {
+
+ if (rowmajor)
+ return i * group_sizes.size() + j;
+
+ return j * group_sizes.size() + i;
+
+ }
+
+public:
+
+ GroupSampler() {};
+
+ GroupSampler(
+ const std::vector< double > & contact_matrix_,
+ const std::vector< size_t > & group_sizes_,
+ bool normalize = true
+ );
+
+ int sample_1(
+ Model * model,
+ const int origin_group
+ );
+
+ void sample_n(
+ Model * model,
+ std::vector< int > & sample,
+ const int origin_group,
+ const int nsamples
+ );
+
+};
+
+#endif
\ No newline at end of file
diff --git a/inst/include/epiworld/groupsampler-meat.hpp b/inst/include/epiworld/groupsampler-meat.hpp
new file mode 100644
index 00000000..6e960965
--- /dev/null
+++ b/inst/include/epiworld/groupsampler-meat.hpp
@@ -0,0 +1,84 @@
+#ifndef GROUPSAMPLER_MEAT_HPP
+#define GROUPSAMPLER_MEAT_HPP
+
+template
+inline GroupSampler::GroupSampler(
+ const std::vector< double > & contact_matrix_,
+ const std::vector< size_t > & group_sizes_,
+ bool normalize
+ ): contact_matrix(contact_matrix_), group_sizes(group_sizes_) {
+
+
+ this->cumulate.resize(contact_matrix.size());
+ std::fill(cumulate.begin(), cumulate.end(), 0.0);
+
+ // Cumulative sum
+ for (size_t j = 0; j < group_sizes.size(); ++j)
+ {
+ for (size_t i = 0; i < group_sizes.size(); ++i)
+ cumulate[idx(i, j, true)] +=
+ cumulate[idx(i, j - 1, true)] +
+ contact_matrix[idx(i, j)];
+ }
+
+ if (normalize)
+ {
+ for (size_t i = 0; i < group_sizes.size(); ++i)
+ {
+ double sum = 0.0;
+ for (size_t j = 0; j < group_sizes.size(); ++j)
+ sum += contact_matrix[idx(i, j, true)];
+ for (size_t j = 0; j < group_sizes.size(); ++j)
+ contact_matrix[idx(i, j, true)] /= sum;
+ }
+ }
+
+ };
+
+template
+int GroupSampler::sample_1(
+ Model * model,
+ const int origin_group
+ )
+{
+
+ // Random number
+ double r = model->runif();
+
+ // Finding the group
+ size_t j = 0;
+ while (r > cumulate[idx(origin_group, j, true)])
+ ++j;
+
+ // Adjusting the prob
+ r = r - (j == 0 ? 0.0 : cumulate[idx(origin_group, j - 1, true)]);
+
+ int res = static_cast(
+ std::floor(r * group_sizes[j])
+ );
+
+ // Making sure we are not picking outside of the group
+ if (res >= static_cast(group_sizes[j]))
+ res = static_cast(group_sizes[j]) - 1;
+
+ return model->get_entities()[j][res]->get_id();
+
+}
+
+template
+void GroupSampler::sample_n(
+ Model * model,
+ std::vector< int > & sample,
+ const int origin_group,
+ const int nsamples
+)
+{
+
+ for (int i = 0; i < nsamples; ++i)
+ sample[i] = sample_1(model, origin_group);
+
+ return;
+
+}
+
+#endif
\ No newline at end of file
diff --git a/inst/include/epiworld/model-bones.hpp b/inst/include/epiworld/model-bones.hpp
index 3b51e2a8..4029f12a 100644
--- a/inst/include/epiworld/model-bones.hpp
+++ b/inst/include/epiworld/model-bones.hpp
@@ -135,14 +135,7 @@ class Model {
bool directed = false;
std::vector< VirusPtr > viruses = {};
- std::vector< epiworld_double > prevalence_virus = {}; ///< Initial prevalence_virus of each virus
- std::vector< bool > prevalence_virus_as_proportion = {};
- std::vector< VirusToAgentFun > viruses_dist_funs = {};
-
std::vector< ToolPtr > tools = {};
- std::vector< epiworld_double > prevalence_tool = {};
- std::vector< bool > prevalence_tool_as_proportion = {};
- std::vector< ToolToAgentFun > tools_dist_funs = {};
std::vector< Entity > entities = {};
std::vector< Entity > entities_backup = {};
@@ -183,7 +176,7 @@ class Model {
void dist_tools();
void dist_virus();
- // void dist_entities();
+ void dist_entities();
std::chrono::time_point time_start;
std::chrono::time_point time_end;
@@ -227,7 +220,7 @@ class Model {
Entity * entity_,
epiworld_fast_int new_state_,
epiworld_fast_int queue_,
- ActionFun call_,
+ EventFun call_,
int idx_agent_,
int idx_object_
);
@@ -258,6 +251,7 @@ class Model {
std::vector array_double_tmp;
std::vector * > array_virus_tmp;
+ std::vector< int > array_int_tmp;
Model();
Model(const Model & m);
@@ -337,16 +331,12 @@ class Model {
* indicating number of individuals.
*/
///@{
- void add_virus(Virus & v, epiworld_double preval);
- void add_virus_n(Virus & v, epiworld_fast_uint preval);
- void add_virus_fun(Virus & v, VirusToAgentFun fun);
- void add_tool(Tool & t, epiworld_double preval);
- void add_tool_n(Tool & t, epiworld_fast_uint preval);
- void add_tool_fun(Tool & t, ToolToAgentFun fun);
+ void add_virus(Virus & v);
+ void add_tool(Tool & t);
void add_entity(Entity e);
void rm_virus(size_t virus_pos);
void rm_tool(size_t tool_pos);
- void rm_entity(size_t entity_pos);
+ void rm_entity(size_t entity_id);
///@}
/**
@@ -360,6 +350,20 @@ class Model {
* @param skip How many rows to skip.
*/
void load_agents_entities_ties(std::string fn, int skip);
+
+ /**
+ * @brief Associate agents-entities from data
+ */
+ void load_agents_entities_ties(
+ const std::vector & agents_ids,
+ const std::vector & entities_ids
+ );
+
+ void load_agents_entities_ties(
+ const int * agents_id,
+ const int * entities_id,
+ size_t n
+ );
/**
* @name Accessing population of the model
@@ -391,6 +395,8 @@ class Model {
std::vector< Agent > & get_agents(); ///< Returns a reference to the vector of agents.
+ Agent & get_agent(size_t i);
+
std::vector< epiworld_fast_uint > get_agents_states() const; ///< Returns a vector with the states of the agents.
std::vector< Viruses_const > get_agents_viruses() const; ///< Returns a const vector with the viruses of the agents.
@@ -399,6 +405,8 @@ class Model {
std::vector< Entity > & get_entities();
+ Entity & get_entity(size_t entity_id, int * entity_pos = nullptr);
+
Model & agents_smallworld(
epiworld_fast_uint n = 1000,
epiworld_fast_uint k = 5,
@@ -525,8 +533,6 @@ class Model {
virtual void reset();
const Model & print(bool lite = false) const;
- Model && clone() const;
-
/**
* @name Manage state (states) in the model
*
@@ -685,8 +691,6 @@ class Model {
///@}
const std::vector< VirusPtr > & get_viruses() const;
- const std::vector< epiworld_double > & get_prevalence_virus() const;
- const std::vector< bool > & get_prevalence_virus_as_proportion() const;
const std::vector< ToolPtr > & get_tools() const;
Virus & get_virus(size_t id);
Tool & get_tool(size_t id);
diff --git a/inst/include/epiworld/model-meat-print.hpp b/inst/include/epiworld/model-meat-print.hpp
index f17f0c00..1a66a85d 100644
--- a/inst/include/epiworld/model-meat-print.hpp
+++ b/inst/include/epiworld/model-meat-print.hpp
@@ -143,6 +143,8 @@ inline const Model