diff --git a/Makefile b/Makefile index a9796723..beeac088 100644 --- a/Makefile +++ b/Makefile @@ -11,9 +11,18 @@ docker-debug: EPI_CONFIG="-DEPI_DEBUG -Wall -pedantic -g" R CMD INSTALL \ --no-docs --build . -install: build - which R - cd ../ && R CMD INSTALL epiworldR_*.tar.gz +install-dev: clean + sed -i -E 's/@useDynLib\s+[a-zA-Z]+/@useDynLib epiworldRdev/g' R/epiworldR-package.R + sed -i -E 's/useDynLib\(+[a-zA-Z]+/useDynLib(epiworldRdev/g' NAMESPACE + sed -i -E 's/^Package:.+/Package: epiworldRdev/g' DESCRIPTION + sed -i -E 's/^library\([a-zA-Z]+\)/library(epiworldRdev)/g' README.* + Rscript --vanilla -e 'roxygen2::roxygenize()' + EPI_DEV=yes R CMD INSTALL .& $(MAKE) clean + +install: + $(MAKE) clean + R CMD INSTALL . + README.md: README.Rmd Rscript --vanilla -e 'rmarkdown::render("README.Rmd")' @@ -29,6 +38,11 @@ check: build clean: Rscript --vanilla -e 'devtools::clean_dll()' + sed -i -E 's/@useDynLib\s+[a-zA-Z]+/@useDynLib epiworldR/g' R/epiworldR-package.R + sed -i -E 's/useDynLib\(+[a-zA-Z]+/useDynLib(epiworldR/g' NAMESPACE + sed -i -E 's/^Package:.+/Package: epiworldR/g' DESCRIPTION + # sed -i -E 's/^\\(name|alias|title)\{[a-zA-Z]+/\\\1{epiworldR/g' man/epiworldR-package.Rd + sed -i -E 's/^library\([a-zA-Z]+\)/library(epiworldR)/g' README.* docs: Rscript --vanilla -e 'devtools::document()' diff --git a/NAMESPACE b/NAMESPACE index 52bf74a0..8acbfbe6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -106,7 +106,6 @@ export(get_agents) export(get_agents_data_ncols) export(get_agents_states) export(get_agents_tools) -export(get_agents_viruses) export(get_generation_time) export(get_hist_tool) export(get_hist_total) @@ -135,6 +134,7 @@ export(globalaction_tool) export(globalaction_tool_logit) export(has_tool) export(has_virus) +export(initial_states) export(make_saver) export(plot_generation_time) export(plot_incidence) diff --git a/R/ModelSEIR.R b/R/ModelSEIR.R index bc1726d5..a081109e 100644 --- a/R/ModelSEIR.R +++ b/R/ModelSEIR.R @@ -12,6 +12,11 @@ #' @export #' @family Models #' @aliases epiworld_seir +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. The user must provide a vector of proportions indicating the following +#' values: (1) Proportion of non-infected agents who are removed, and (2) +#' Proportion of exposed agents to be set as infected. #' @returns #' - The `ModelSEIR`function returns a model of class [epiworld_model]. #' @examples diff --git a/R/ModelSEIRD.R b/R/ModelSEIRD.R index 326193b5..fc1394e4 100644 --- a/R/ModelSEIRD.R +++ b/R/ModelSEIRD.R @@ -11,6 +11,12 @@ #' @param x Object of class SEIRD. #' @param ... Currently ignore. #' @export +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. The user must provide a vector of proportions indicating the following +#' values: (1) Proportion of exposed agents who are infected, (2) +#' proportion of non-infected agents already removed, and (3) proportion of +#' non-ifected agents already deceased. #' @family Models #' @aliases epiworld_seird #' @returns diff --git a/R/ModelSEIRDCONN.R b/R/ModelSEIRDCONN.R index 2432c16a..44aed7ba 100644 --- a/R/ModelSEIRDCONN.R +++ b/R/ModelSEIRDCONN.R @@ -17,6 +17,12 @@ #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. The user must provide a vector of proportions indicating the following +#' values: (1) Proportion of exposed agents who are infected, (2) +#' proportion of non-infected agents already removed, and (3) proportion of +#' non-ifected agents already deceased. #' @family Models #' @aliases epiworld_seirdconn #' @returns diff --git a/R/ModelSIR.R b/R/ModelSIR.R index c2eaec5b..76a38297 100644 --- a/R/ModelSIR.R +++ b/R/ModelSIR.R @@ -13,6 +13,10 @@ #' @param ... Additional arguments passed to [graphics::plot]. #' @export #' @family Models +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. In particular, the user can specify how many of the non-infected +#' agents have been removed at the beginning of the simulation. #' @aliases epiworld_sir #' @returns #' - The `ModelSIR` function returns a model of class [epiworld_model]. diff --git a/R/ModelSIRCONN.R b/R/ModelSIRCONN.R index 69bab0c7..6a92b8a2 100644 --- a/R/ModelSIRCONN.R +++ b/R/ModelSIRCONN.R @@ -10,6 +10,10 @@ #' @param n Number of individuals in the population. #' @export #' @family Models +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. In particular, the user can specify how many of the non-infected +#' agents have been removed at the beginning of the simulation. #' @returns #' - The `ModelSIRCONN`function returns a model of class [epiworld_model]. #' @aliases epiworld_sirconn diff --git a/R/ModelSIRD.R b/R/ModelSIRD.R index e7aba713..08272f17 100644 --- a/R/ModelSIRD.R +++ b/R/ModelSIRD.R @@ -13,6 +13,11 @@ #' @param x Object of class SIR. #' @param ... Additional arguments passed to [graphics::plot]. #' @export +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. The user must provide a vector of proportions indicating the following +#' values: (1) proportion of non-infected agents already removed, and (2) proportion of +#' non-ifected agents already deceased. #' @family Models #' @aliases epiworld_sird #' @returns diff --git a/R/ModelSIRDCONN.R b/R/ModelSIRDCONN.R index 0e745b88..b4b3fb13 100644 --- a/R/ModelSIRDCONN.R +++ b/R/ModelSIRDCONN.R @@ -10,6 +10,11 @@ #' @param ... Currently ignore. #' @param n Number of individuals in the population. #' @export +#' @details +#' The [initial_states] function allows the user to set the initial state of the +#' model. The user must provide a vector of proportions indicating the following +#' values: (1) proportion of non-infected agents already removed, and (2) proportion of +#' non-ifected agents already deceased. #' @family Models #' @returns #' - The `ModelSIRDCONN`function returns a model of class [epiworld_model]. diff --git a/R/cpp11.R b/R/cpp11.R index fd169573..ee345e3c 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -256,6 +256,10 @@ get_network_cpp <- function(model) { .Call(`_epiworldR_get_network_cpp`, model) } +initial_states_cpp <- function(model, proportions) { + .Call(`_epiworldR_initial_states_cpp`, model, proportions) +} + 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) } @@ -427,11 +431,3 @@ get_name_virus_cpp <- function(virus) { set_name_virus_cpp <- function(virus, name) { .Call(`_epiworldR_set_name_virus_cpp`, virus, name) } - -get_agents_viruses_cpp <- function(model) { - .Call(`_epiworldR_get_agents_viruses_cpp`, model) -} - -print_agent_viruses_cpp <- function(viruses) { - .Call(`_epiworldR_print_agent_viruses_cpp`, viruses) -} diff --git a/R/epiworldR-package.R.in b/R/epiworldR-package.R.in new file mode 100644 index 00000000..21204cbb --- /dev/null +++ b/R/epiworldR-package.R.in @@ -0,0 +1,5 @@ +#' @EPIWORLD_NAME@ +#' @useDynLib @EPIWORLD_NAME@, .registration = TRUE +#' @importFrom graphics boxplot plot +"_PACKAGE" + diff --git a/R/model-methods.R b/R/model-methods.R index e6f6f6fe..a4f27bbd 100644 --- a/R/model-methods.R +++ b/R/model-methods.R @@ -351,3 +351,17 @@ get_tool <- function(model, tool_pos) { class = "epiworld_tool" ) } + +#' @export +#' @param proportions Numeric vector. Proportions in which agents will be +#' distributed (see details). +#' @return +#' - `inital_states` returns the model with an updated initial state. +#' @rdname epiworld-methods +initial_states <- function(model, proportions) { + + stopifnot_model(model) + invisible(initial_states_cpp(model, proportions)) + +} + diff --git a/R/virus.R b/R/virus.R index 7770e115..fe49564b 100644 --- a/R/virus.R +++ b/R/virus.R @@ -536,25 +536,6 @@ set_incubation_fun <- function(virus, model, vfun) { } -#' @export -#' @rdname agents_smallworld -#' @returns -#' - `get_agents_viruses` returns a list of class `epiworld_agents_viruses` -#' with `epiworld_viruses` (list of lists). -get_agents_viruses <- function(model) { - - stopifnot_model(model) - - res <- lapply( - get_agents_viruses_cpp(model), - `class<-`, - "epiworld_viruses" - ) - - structure(res, class = c("epiworld_agents_viruses", class(res))) - -} - #' @export #' @rdname virus #' @param max_print Numeric scalar. Maximum number of viruses to print. diff --git a/README.Rmd b/README.Rmd index fd801bc0..fffefcb1 100644 --- a/README.Rmd +++ b/README.Rmd @@ -71,10 +71,10 @@ library(epiworldR) # Creating a SIR model sir <- ModelSIR( - name = "COVID-19", - prevalence = .01, + name = "COVID-19", + prevalence = .01, transmission_rate = .7, - recovery = .3 + recovery = .3 ) |> # Adding a Small world population agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |> @@ -204,9 +204,11 @@ sir <- ModelSIR( net <- get_transmissions(sir) # Plotting +library(epiworldR) library(netplot) -library(igraph) -x <- graph_from_edgelist(as.matrix(net[,2:3]) + 1) +x <- igraph::graph_from_edgelist( + as.matrix(net[,2:3]) + 1 + ) nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE) ``` diff --git a/README.md b/README.md index 39213bd1..05625989 100644 --- a/README.md +++ b/README.md @@ -91,10 +91,10 @@ library(epiworldR) # Creating a SIR model sir <- ModelSIR( - name = "COVID-19", - prevalence = .01, + name = "COVID-19", + prevalence = .01, transmission_rate = .7, - recovery = .3 + recovery = .3 ) |> # Adding a Small world population agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |> @@ -127,8 +127,8 @@ summary(sir) #> Number of entities : 0 #> Days (duration) : 50 (of 50) #> Number of viruses : 1 -#> Last run elapsed t : 196.00ms -#> Last run speed : 25.43 million agents x day / second +#> Last run elapsed t : 166.00ms +#> Last run speed : 30.01 million agents x day / second #> Rewiring : off #> #> Global actions: @@ -303,16 +303,16 @@ rn <- get_reproductive_number(model_logit) # Looking into the agents get_agents(model_logit) #> Agents from the model "Susceptible-Infected-Removed (SIR) (logit)": -#> Agent: 0, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 1, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 2, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 3, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 4, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 5, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 6, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 7, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 8, state: Susceptible (0), Nvirus: 0, NTools: 0, NNeigh: 8 -#> Agent: 9, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8 +#> Agent: 0, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 1, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 2, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 3, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 4, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 5, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 6, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 7, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 8, state: Susceptible (0), Has virus: no, NTools: 0, NNeigh: 8 +#> Agent: 9, state: Recovered (2), Has virus: no, NTools: 0, NNeigh: 8 #> ... 99990 more agents ... ``` @@ -345,21 +345,12 @@ sir <- ModelSIR( net <- get_transmissions(sir) # Plotting +library(epiworldR) library(netplot) #> Loading required package: grid -library(igraph) -#> -#> Attaching package: 'igraph' -#> The following object is masked from 'package:netplot': -#> -#> ego -#> The following objects are masked from 'package:stats': -#> -#> decompose, spectrum -#> The following object is masked from 'package:base': -#> -#> union -x <- graph_from_edgelist(as.matrix(net[,2:3]) + 1) +x <- igraph::graph_from_edgelist( + as.matrix(net[,2:3]) + 1 + ) nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE) ``` @@ -411,10 +402,10 @@ head(ans$reproductive) #> sim_num virus_id virus source source_exposure_date rt #> 1 1 0 COVID-19 767 11 0 #> 2 1 0 COVID-19 835 10 0 -#> 3 1 0 COVID-19 466 9 0 +#> 3 1 0 COVID-19 793 9 0 #> 4 1 0 COVID-19 612 9 0 -#> 5 1 0 COVID-19 793 9 0 -#> 6 1 0 COVID-19 20 8 0 +#> 5 1 0 COVID-19 466 9 0 +#> 6 1 0 COVID-19 920 8 0 plot(ans$reproductive) ``` diff --git a/configure b/configure index 4b715d07..65317b69 100755 --- a/configure +++ b/configure @@ -617,6 +617,7 @@ PACKAGE_URL='' ac_subst_vars='LTLIBOBJS LIBOBJS +EPIWORLD_NAME OPENMP_FLAG CXXCPP OBJEXT @@ -3365,11 +3366,24 @@ printf "%s\n" "found and suitable" >&6; } fi fi +# Checking if the development version is being built +{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking whether to build the dev version" >&5 +printf %s "checking whether to build the dev version... " >&6; } +epiworldname="epiworldR" +if test x"${EPI_DEV}" != x"" ; then + epiworldname="epiworldRdev" +fi ## now use all these OPENMP_FLAG="${openmp_flag}" -ac_config_files="$ac_config_files src/Makevars" +EPIWORLD_NAME="${epiworldname}" + +ac_config_files="$ac_config_files src/Makevars R/epiworldR-package.R" + + +# Use sed to replace the line useDynLib\([a-zA-Z]+ in the NAMESPACE +# file with useDynLib(${epiworldname} cat >confcache <<\_ACEOF # This file is a shell script that caches the results of configure @@ -4072,6 +4086,7 @@ for ac_config_target in $ac_config_targets do case $ac_config_target in "src/Makevars") CONFIG_FILES="$CONFIG_FILES src/Makevars" ;; + "R/epiworldR-package.R") CONFIG_FILES="$CONFIG_FILES R/epiworldR-package.R" ;; *) as_fn_error $? "invalid argument: \`$ac_config_target'" "$LINENO" 5;; esac diff --git a/configure.ac b/configure.ac index 1a541949..97048765 100644 --- a/configure.ac +++ b/configure.ac @@ -169,8 +169,19 @@ if test x"${can_use_openmp}" = x"yes"; then fi fi +# Checking if the development version is being built +AC_MSG_CHECKING([whether to build the dev version]) +epiworldname="epiworldR" +if test x"${EPI_DEV}" != x"" ; then + epiworldname="epiworldRdev" +fi ## now use all these AC_SUBST([OPENMP_FLAG], ["${openmp_flag}"]) -AC_CONFIG_FILES([src/Makevars]) +AC_SUBST([EPIWORLD_NAME], ["${epiworldname}"]) +AC_CONFIG_FILES([src/Makevars R/epiworldR-package.R]) + +# Use sed to replace the line useDynLib\([a-zA-Z]+ in the NAMESPACE +# file with useDynLib(${epiworldname} + AC_OUTPUT \ No newline at end of file diff --git a/inst/include/epiworld/agent-actions-meat.hpp b/inst/include/epiworld/agent-actions-meat.hpp index 20cc829c..ef078cbb 100644 --- a/inst/include/epiworld/agent-actions-meat.hpp +++ b/inst/include/epiworld/agent-actions-meat.hpp @@ -8,9 +8,6 @@ inline void default_add_virus(Action & a, Model * m) Agent * p = a.agent; VirusPtr v = a.virus; - CHECK_COALESCE_(a.new_state, v->state_init, p->get_state()) - CHECK_COALESCE_(a.queue, v->queue_init, 1) - // Has a agent? If so, we need to register the transmission if (v->get_agent()) { @@ -26,30 +23,22 @@ inline void default_add_virus(Action & a, Model * m) } - // Update virus accounting - p->n_viruses++; - size_t n_viruses = p->n_viruses; - - if (n_viruses <= p->viruses.size()) - p->viruses[n_viruses - 1] = std::make_shared< Virus >(*v); - else - p->viruses.push_back(std::make_shared< Virus >(*v)); - - n_viruses--; + p->virus = std::make_shared< Virus >(*v); + p->virus->set_date(m->today()); + p->virus->set_agent(p); - // Notice that both agent and date can be changed in this case - // as only the sequence is a shared_ptr itself. - #ifdef EPI_DEBUG - if (n_viruses >= p->viruses.size()) + // Change of state needs to be recorded and updated on the + // tools. + if (p->state_prev != p->state) { - throw std::logic_error( - "[epi-debug]::default_add_virus Index for new virus out of range." - ); + auto & db = m->get_db(); + db.update_state(p->state_prev, p->state); + + for (size_t i = 0u; i < p->n_tools; ++i) + db.update_tool(p->tools[i]->get_id(), p->state_prev, p->state); } - #endif - p->viruses[n_viruses]->set_agent(p, n_viruses); - p->viruses[n_viruses]->set_date(m->today()); + // Lastly, we increase the daily count of the virus #ifdef EPI_DEBUG m->get_db().today_virus.at(v->get_id()).at(p->state)++; #else @@ -64,9 +53,6 @@ inline void default_add_tool(Action & a, Model * m) Agent * p = a.agent; ToolPtr t = a.tool; - - CHECK_COALESCE_(a.new_state, t->state_init, p->get_state()) - CHECK_COALESCE_(a.queue, t->queue_init, Queue::NoOne) // Update tool accounting p->n_tools++; @@ -82,47 +68,64 @@ inline void default_add_tool(Action & a, Model * m) p->tools[n_tools]->set_date(m->today()); p->tools[n_tools]->set_agent(p, n_tools); + // Change of state needs to be recorded and updated on the + // tools. + if (p->state_prev != p->state) + { + auto & db = m->get_db(); + db.update_state(p->state_prev, p->state); + + if (p->virus) + db.update_virus(p->virus->get_id(), p->state_prev, p->state); + } + m->get_db().today_tool[t->get_id()][p->state]++; + } template inline void default_rm_virus(Action & a, Model * model) { - Agent * p = a.agent; - VirusPtr & v = a.agent->viruses[a.virus->pos_in_agent]; - - CHECK_COALESCE_(a.new_state, v->state_post, p->get_state()) - CHECK_COALESCE_(a.queue, v->queue_post, -Queue::Everyone) + Agent * p = a.agent; + VirusPtr & v = a.virus; - if (--p->n_viruses > 0) - { - // The new virus will change positions - p->viruses[p->n_viruses]->pos_in_agent = v->pos_in_agent; - std::swap( - p->viruses[p->n_viruses], // Moving to the end - p->viruses[v->pos_in_agent] // Moving to the beginning - ); - } - // Calling the virus action over the removed virus v->post_recovery(model); + p->virus = nullptr; + + // Change of state needs to be recorded and updated on the + // tools. + if (p->state_prev != p->state) + { + auto & db = model->get_db(); + db.update_state(p->state_prev, p->state); + + for (size_t i = 0u; i < p->n_tools; ++i) + db.update_tool(p->tools[i]->get_id(), p->state_prev, p->state); + } + + // The counters of the virus only needs to decrease + #ifdef EPI_DEBUG + model->get_db().today_virus.at(v->get_id()).at(p->state_prev)--; + #else + model->get_db().today_virus[v->get_id()][p->state_prev]--; + #endif + + return; } template -inline void default_rm_tool(Action & a, Model * /*m*/) +inline void default_rm_tool(Action & a, Model * m) { Agent * p = a.agent; ToolPtr & t = a.agent->tools[a.tool->pos_in_agent]; - CHECK_COALESCE_(a.new_state, t->state_post, p->get_state()) - CHECK_COALESCE_(a.queue, t->queue_post, Queue::NoOne) - if (--p->n_tools > 0) { p->tools[p->n_tools]->pos_in_agent = t->pos_in_agent; @@ -132,10 +135,49 @@ inline void default_rm_tool(Action & a, Model * /*m*/) ); } + // Change of state needs to be recorded and updated on the + // tools. + if (p->state_prev != p->state) + { + auto & db = m->get_db(); + db.update_state(p->state_prev, p->state); + + if (p->virus) + db.update_virus(p->virus->get_id(), p->state_prev, p->state); + } + + // Lastly, we increase the daily count of the tool + #ifdef EPI_DEBUG + m->get_db().today_tool.at(t->get_id()).at(p->state_prev)--; + #else + m->get_db().today_tool[t->get_id()][p->state_prev]--; + #endif + return; } +template +inline void default_change_state(Action & a, Model * m) +{ + + Agent * p = a.agent; + + if (p->state_prev != p->state) + { + auto & db = m->get_db(); + db.update_state(p->state_prev, p->state); + + if (p->virus) + db.update_virus(p->virus->get_id(), p->state_prev, p->state); + + for (size_t i = 0u; i < p->n_tools; ++i) + db.update_tool(p->tools[i]->get_id(), p->state_prev, p->state); + + } + +} + template inline void default_add_entity(Action & a, Model *) { @@ -143,9 +185,6 @@ inline void default_add_entity(Action & a, Model *) Agent * p = a.agent; Entity * e = a.entity; - CHECK_COALESCE_(a.new_state, e->state_post, p->get_state()) - CHECK_COALESCE_(a.queue, e->queue_post, Queue::NoOne) - // Checking the agent and the entity are not linked if ((p->get_n_entities() > 0) && (e->size() > 0)) { @@ -208,9 +247,6 @@ inline void default_rm_entity(Action & a, Model * m) size_t idx_agent_in_entity = a.idx_agent; size_t idx_entity_in_agent = a.idx_object; - CHECK_COALESCE_(a.new_state, e->state_post, p->get_state()) - CHECK_COALESCE_(a.queue, e->queue_post, Queue::NoOne) - if (--p->n_entities > 0) { diff --git a/inst/include/epiworld/agent-bones.hpp b/inst/include/epiworld/agent-bones.hpp index 422f071f..06bc77ce 100644 --- a/inst/include/epiworld/agent-bones.hpp +++ b/inst/include/epiworld/agent-bones.hpp @@ -52,6 +52,9 @@ inline void default_rm_tool(Action & a, Model * m); template inline void default_rm_entity(Action & a, Model * m); +template +inline void default_change_state(Action & a, Model * m); + /** @@ -63,8 +66,6 @@ template class Agent { friend class Model; friend class Virus; - friend class Viruses; - friend class Viruses_const; friend class Tool; friend class Tools; friend class Tools_const; @@ -77,6 +78,7 @@ class Agent { friend void default_rm_virus(Action & a, Model * m); friend void default_rm_tool(Action & a, Model * m); friend void default_rm_entity(Action & a, Model * m); + friend void default_change_state(Action & a, Model * m); private: Model * model; @@ -95,22 +97,11 @@ class Agent { int state_last_changed = -1; ///< Last time the agent was updated. int id = -1; - std::vector< VirusPtr > viruses; - epiworld_fast_uint n_viruses = 0u; + VirusPtr virus = nullptr; std::vector< ToolPtr > tools; epiworld_fast_uint n_tools = 0u; - ActionFun add_virus_ = default_add_virus; - ActionFun add_tool_ = default_add_tool; - ActionFun add_entity_ = default_add_entity; - - ActionFun rm_virus_ = default_rm_virus; - ActionFun rm_tool_ = default_rm_tool; - ActionFun rm_entity_ = default_rm_entity; - - epiworld_fast_uint action_counter = 0u; - std::vector< Agent * > sampled_agents; size_t sampled_agents_n = 0u; std::vector< size_t > sampled_agents_left; @@ -149,14 +140,14 @@ class Agent { epiworld_fast_int queue = -99 ); - void add_virus( + void set_virus( VirusPtr virus, Model * model, epiworld_fast_int state_new = -99, epiworld_fast_int queue = -99 ); - void add_virus( + void set_virus( Virus virus, Model * model, epiworld_fast_int state_new = -99, @@ -185,14 +176,6 @@ class Agent { ); void rm_virus( - epiworld_fast_uint virus_idx, - Model * model, - epiworld_fast_int state_new = -99, - epiworld_fast_int queue = -99 - ); - - void rm_virus( - VirusPtr & virus, Model * model, epiworld_fast_int state_new = -99, epiworld_fast_int queue = -99 @@ -213,14 +196,6 @@ class Agent { ); void rm_agent_by_virus( - epiworld_fast_uint virus_idx, - Model * model, - epiworld_fast_int state_new = -99, - epiworld_fast_int queue = -99 - ); ///< Agent removed by virus - - void rm_agent_by_virus( - VirusPtr & virus, Model * model, epiworld_fast_int state_new = -99, epiworld_fast_int queue = -99 @@ -242,10 +217,7 @@ class Agent { int get_id() const; ///< Id of the individual - VirusPtr & get_virus(int i); - Viruses get_viruses(); - const Viruses_const get_viruses() const; - size_t get_n_viruses() const noexcept; + VirusPtr & get_virus(); ToolPtr & get_tool(int i); Tools get_tools(); diff --git a/inst/include/epiworld/agent-meat-state.hpp b/inst/include/epiworld/agent-meat-state.hpp index e45c5646..91133f4a 100644 --- a/inst/include/epiworld/agent-meat-state.hpp +++ b/inst/include/epiworld/agent-meat-state.hpp @@ -33,7 +33,7 @@ inline void default_update_susceptible( if (virus == nullptr) return; - p->add_virus(*virus, m); + p->set_virus(*virus, m); return; @@ -42,59 +42,37 @@ inline void default_update_susceptible( template inline void default_update_exposed(Agent * p, Model * m) { - if (p->get_n_viruses() == 0u) + if (p->get_virus() == nullptr) throw std::logic_error( std::string("Using the -default_update_exposed- on agents WITHOUT viruses makes no sense! ") + std::string("Agent id ") + std::to_string(p->get_id()) + std::string(" has no virus registered.") ); - // Odd: Die, Even: Recover - epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { - - // Die - m->array_double_tmp[n_events++] = - v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); - - // Recover - m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); - - } + // Die + auto & virus = p->get_virus(); + m->array_double_tmp[0u] = + virus->get_prob_death(m) * (1.0 - p->get_death_reduction(virus, m)); - #ifdef EPI_DEBUG - if (n_events == 0u) - { - printf_epiworld( - "[epi-debug] agent %i has 0 possible events!!\n", - static_cast(p->get_id()) - ); - throw std::logic_error("Zero events in exposed."); - } - #else - if (n_events == 0u) - return; - #endif + // Recover + m->array_double_tmp[1u] = + 1.0 - (1.0 - virus->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(virus, m)); // Running the roulette - int which = roulette(n_events, m); + int which = roulette(2u, m); if (which < 0) return; // Which roulette happen? - if ((which % 2) == 0) // If odd + if (which == 0u) // If odd { - size_t which_v = std::ceil(which / 2); - p->rm_agent_by_virus(which_v, m); + p->rm_agent_by_virus(m); } else { - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + p->rm_virus(m); } diff --git a/inst/include/epiworld/agent-meat-virus-sampling.hpp b/inst/include/epiworld/agent-meat-virus-sampling.hpp index 82096206..d7289937 100644 --- a/inst/include/epiworld/agent-meat-virus-sampling.hpp +++ b/inst/include/epiworld/agent-meat-virus-sampling.hpp @@ -34,37 +34,31 @@ inline std::function*,Model*)> make_update_susceptible( [](Agent * p, Model * m) -> void { - if (p->get_n_viruses() > 0u) + if (p->get_virus() != nullptr) throw std::logic_error( std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense! ") + std::string("Agent id ") + std::to_string(p->get_id()) + - std::string(" has ") + std::to_string(p->get_n_viruses()) + - std::string(" viruses.") + std::string(" has a virus.") ); // This computes the prob of getting any neighbor variant size_t nviruses_tmp = 0u; for (auto & neighbor: p->get_neighbors()) { - - for (const VirusPtr & v : neighbor->get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; - m->array_virus_tmp[nviruses_tmp++] = &(*v); + auto & v = neighbor->get_virus(); + if (v == nullptr) + continue; + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); - } } // No virus to compute @@ -77,7 +71,7 @@ inline std::function*,Model*)> make_update_susceptible( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; }; @@ -118,12 +112,11 @@ inline std::function*,Model*)> make_update_susceptible( } - if (p->get_n_viruses() > 0u) + if (p->get_virus() != nullptr) throw std::logic_error( std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense! ") + std::string("Agent id ") + std::to_string(p->get_id()) + - std::string(" has ") + std::to_string(p->get_n_viruses()) + - std::string(" viruses.") + std::string(" has a virus.") ); // This computes the prob of getting any neighbor variant @@ -134,26 +127,21 @@ inline std::function*,Model*)> make_update_susceptible( // If the state is in the list, exclude it if (exclude_agent_bool->operator[](neighbor->get_state())) continue; - - for (const VirusPtr & v : neighbor->get_viruses()) - { - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - - #endif + auto & v = neighbor->get_virus(); + if (v == nullptr) + continue; - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); - m->array_virus_tmp[nviruses_tmp++] = &(*v); - - } } // No virus to compute @@ -166,7 +154,7 @@ inline std::function*,Model*)> make_update_susceptible( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; @@ -202,37 +190,37 @@ inline std::function*(Agent*,Model*)> make_sample_virus_ std::function*(Agent*,Model*)> res = [](Agent * p, Model * m) -> Virus* { - if (p->get_n_viruses() > 0u) + if (p->get_virus() != nullptr) throw std::logic_error( std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense! ") + std::string("Agent id ") + std::to_string(p->get_id()) + - std::string(" has ") + std::to_string(p->get_n_viruses()) + - std::string(" viruses.") + std::string(" has a virus.") ); // This computes the prob of getting any neighbor variant size_t nviruses_tmp = 0u; for (auto & neighbor: p->get_neighbors()) { - - for (const VirusPtr & v : neighbor->get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; - m->array_virus_tmp[nviruses_tmp++] = &(*v); + if (neighbor->get_virus() == nullptr) + continue; + + auto & v = neighbor->get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + } // No virus to compute @@ -286,12 +274,11 @@ inline std::function*(Agent*,Model*)> make_sample_virus_ } - if (p->get_n_viruses() > 0u) + if (p->get_virus() != nullptr) throw std::logic_error( std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense! ") + std::string("Agent id ") + std::to_string(p->get_id()) + - std::string(" has ") + std::to_string(p->get_n_viruses()) + - std::string(" viruses.") + std::string(" has a virus.") ); // This computes the prob of getting any neighbor variant @@ -302,25 +289,26 @@ inline std::function*(Agent*,Model*)> make_sample_virus_ // If the state is in the list, exclude it if (exclude_agent_bool->operator[](neighbor->get_state())) continue; - - for (const VirusPtr & v : neighbor->get_viruses()) - { - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif + if (neighbor->get_virus() == nullptr) + continue; + + auto & v = neighbor->get_virus(); - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + } // No virus to compute @@ -363,12 +351,11 @@ template inline Virus * sample_virus_single(Agent * p, Model * m) { - if (p->get_n_viruses() > 0u) + if (p->get_virus() != nullptr) throw std::logic_error( - std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense! ") + + std::string("Using the -default_update_susceptible- on agents WITH viruses makes no sense!") + std::string("Agent id ") + std::to_string(p->get_id()) + - std::string(" has ") + std::to_string(p->get_n_viruses()) + - std::string(" viruses.") + std::string(" has a virus.") ); // This computes the prob of getting any neighbor variant @@ -377,40 +364,42 @@ inline Virus * sample_virus_single(Agent * p, Model * m) { #ifdef EPI_DEBUG int _vcount_neigh = 0; - #endif - for (const VirusPtr & v : neighbor->get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= m->array_virus_tmp.size()) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); - - #ifdef EPI_DEBUG - if ( - (m->array_double_tmp[nviruses_tmp - 1] < 0.0) | - (m->array_double_tmp[nviruses_tmp - 1] > 1.0) - ) - { - printf_epiworld( - "[epi-debug] Agent %i's virus %i has transmission prob outside of [0, 1]: %.4f!\n", - static_cast(neighbor->get_id()), - static_cast(_vcount_neigh++), - m->array_double_tmp[nviruses_tmp - 1] - ); - } - #endif + #endif + + if (neighbor->get_virus() == nullptr) + continue; + + auto & v = neighbor->get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= m->array_virus_tmp.size()) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + + #ifdef EPI_DEBUG + if ( + (m->array_double_tmp[nviruses_tmp - 1] < 0.0) | + (m->array_double_tmp[nviruses_tmp - 1] > 1.0) + ) + { + printf_epiworld( + "[epi-debug] Agent %i's virus %i has transmission prob outside of [0, 1]: %.4f!\n", + static_cast(neighbor->get_id()), + static_cast(_vcount_neigh++), + m->array_double_tmp[nviruses_tmp - 1] + ); + } + #endif - } } diff --git a/inst/include/epiworld/agent-meat.hpp b/inst/include/epiworld/agent-meat.hpp index 47f69626..6bed922b 100644 --- a/inst/include/epiworld/agent-meat.hpp +++ b/inst/include/epiworld/agent-meat.hpp @@ -1,9 +1,6 @@ #ifndef EPIWORLD_PERSON_MEAT_HPP #define EPIWORLD_PERSON_MEAT_HPP -// template -// inline bool IN(Ta & a, std::vector< Ta > & b); - #define CHECK_COALESCE_(proposed_, virus_tool_, alt_) \ if (static_cast(proposed_) == -99) {\ if (static_cast(virus_tool_) == -99) \ @@ -29,36 +26,21 @@ inline Agent::Agent(Agent && p) : state_prev(p.state_prev), state_last_changed(p.state_last_changed), id(p.id), - viruses(std::move(p.viruses)), /// Needs to be adjusted - n_viruses(p.n_viruses), tools(std::move(p.tools)), /// Needs to be adjusted - n_tools(p.n_tools), - add_virus_(std::move(p.add_virus_)), - add_tool_(std::move(p.add_tool_)), - add_entity_(std::move(p.add_entity_)), - rm_virus_(std::move(p.rm_virus_)), - rm_tool_(std::move(p.rm_tool_)), - rm_entity_(std::move(p.rm_entity_)), - action_counter(p.action_counter) + n_tools(p.n_tools) { state = p.state; id = p.id; // Dealing with the virus - - int loc = 0; - for (auto & v : viruses) + if (p.virus != nullptr) { - - // Will create a copy of the virus, with the exeption of - // the virus code - v->agent = this; - v->pos_in_agent = loc++; - + virus = std::move(p.virus); + virus->set_agent(this); } - loc = 0; + int loc = 0; for (auto & t : tools) { @@ -90,17 +72,12 @@ inline Agent::Agent(const Agent & p) : id = p.id; // Dealing with the virus - viruses.reserve(p.get_n_viruses()); - n_viruses = viruses.size(); - for (size_t i = 0u; i < n_viruses; ++i) + if (p.virus != nullptr) { - - // Will create a copy of the virus, with the exeption of - // the virus code - viruses.emplace_back(std::make_shared>(*p.viruses[i])); - viruses.back()->set_agent(this, i); - + virus = std::make_shared>(*p.virus); + virus->set_agent(this); } + tools.reserve(p.get_n_tools()); n_tools = tools.size(); @@ -113,11 +90,6 @@ inline Agent::Agent(const Agent & p) : tools.back()->set_agent(this, i); } - - add_virus_ = p.add_virus_; - add_tool_ = p.add_tool_; - rm_virus_ = p.rm_virus_; - rm_tool_ = p.rm_tool_; } @@ -151,14 +123,12 @@ inline Agent & Agent::operator=( state_last_changed = other_agent.state_last_changed; id = other_agent.id; - // viruses = other_agent.viruses; - n_viruses = other_agent.n_viruses; - viruses.resize(n_viruses, nullptr); - for (size_t i = 0u; i < n_viruses; ++i) + if (other_agent.virus != nullptr) { - viruses[i] = std::make_shared>(*other_agent.viruses[i]); - viruses[i]->set_agent(this, i); - } + virus = std::make_shared>(*other_agent.virus); + virus->set_agent(this); + } else + virus = nullptr; // tools = other_agent.tools; n_tools = other_agent.n_tools; @@ -167,14 +137,6 @@ inline Agent & Agent::operator=( tools[i] = std::make_shared>(*other_agent.tools[i]); tools[i]->set_agent(this, i); } - - add_virus_ = other_agent.add_virus_; - add_tool_ = other_agent.add_tool_; - add_entity_ = other_agent.add_entity_; - rm_virus_ = other_agent.rm_virus_; - rm_tool_ = other_agent.rm_tool_; - rm_entity_ = other_agent.rm_entity_; - action_counter = other_agent.action_counter; return *this; @@ -193,10 +155,12 @@ inline void Agent::add_tool( throw std::range_error("The tool with id: " + std::to_string(tool->get_id()) + " has not been registered. There are only " + std::to_string(model->get_n_tools()) + " included in the model."); - + + CHECK_COALESCE_(state_new, tool->state_init, state); + CHECK_COALESCE_(queue, tool->queue_init, Queue::NoOne); model->actions_add( - this, nullptr, tool, nullptr, state_new, queue, add_tool_, -1, -1 + this, nullptr, tool, nullptr, state_new, queue, default_add_tool, -1, -1 ); } @@ -214,7 +178,7 @@ inline void Agent::add_tool( } template -inline void Agent::add_virus( +inline void Agent::set_virus( VirusPtr virus, Model * model, epiworld_fast_int state_new, @@ -228,14 +192,17 @@ inline void Agent::add_virus( " has not been registered. There are only " + std::to_string(model->get_n_viruses()) + " included in the model."); + CHECK_COALESCE_(state_new, virus->state_init, state); + CHECK_COALESCE_(queue, virus->queue_init, Queue::NoOne); + model->actions_add( - this, virus, nullptr, nullptr, state_new, queue, add_virus_, -1, -1 + this, virus, nullptr, nullptr, state_new, queue, default_add_virus, -1, -1 ); } template -inline void Agent::add_virus( +inline void Agent::set_virus( Virus virus, Model * model, epiworld_fast_int state_new, @@ -243,7 +210,7 @@ inline void Agent::add_virus( ) { VirusPtr virus_ptr = std::make_shared< Virus >(virus); - add_virus(virus_ptr, model, state_new, queue); + set_virus(virus_ptr, model, state_new, queue); } template @@ -255,11 +222,14 @@ inline void Agent::add_entity( ) { + CHECK_COALESCE_(state_new, entity.state_init, state); + CHECK_COALESCE_(queue, entity.queue_init, Queue::NoOne); + if (model != nullptr) { model->actions_add( - this, nullptr, nullptr, &entity, state_new, queue, add_entity_, -1, -1 + this, nullptr, nullptr, &entity, state_new, queue, default_add_entity, -1, -1 ); } @@ -268,11 +238,11 @@ inline void Agent::add_entity( { Action a( - this, nullptr, nullptr, &entity, state_new, queue, add_entity_, + this, nullptr, nullptr, &entity, state_new, queue, default_add_entity, -1, -1 ); - default_add_entity(a, model); /* passing model makes nothing */ + // default_add_entity(a, model); /* passing model makes nothing */ } @@ -287,6 +257,9 @@ inline void Agent::rm_tool( ) { + CHECK_COALESCE_(state_new, tools[tool_idx]->state_post, state); + CHECK_COALESCE_(queue, tools[tool_idx]->queue_post, Queue::NoOne); + if (tool_idx >= n_tools) throw std::range_error( "The Tool you want to remove is out of range. This Agent only has " + @@ -294,7 +267,7 @@ inline void Agent::rm_tool( ); model->actions_add( - this, nullptr, tools[tool_idx], nullptr, state_new, queue, rm_tool_, -1, -1 + this, nullptr, tools[tool_idx], nullptr, state_new, queue, default_rm_tool, -1, -1 ); } @@ -312,63 +285,32 @@ inline void Agent::rm_tool( throw std::logic_error("Cannot remove a virus from another agent!"); model->actions_add( - this, nullptr, tool, nullptr, state_new, queue, rm_tool_, -1, -1 + this, nullptr, tool, nullptr, state_new, queue, default_rm_tool, -1, -1 ); } template inline void Agent::rm_virus( - epiworld_fast_uint virus_idx, Model * model, epiworld_fast_int state_new, epiworld_fast_int queue ) { - if (virus_idx >= n_viruses) - throw std::range_error( - "The Virus you want to remove is out of range. This Agent only has " + - std::to_string(n_viruses) + " viruses." - ); - else if (n_viruses == 0u) - throw std::logic_error( - "There is no virus to remove here!" - ); - #ifdef EPI_DEBUG - if (viruses[virus_idx]->pos_in_agent >= static_cast(n_viruses)) - { + if (virus == nullptr) throw std::logic_error( - "[epi-debug]::rm_virus the position in the agent is wrong." - ); - } - #endif - - model->actions_add( - this, viruses[virus_idx], nullptr, nullptr, state_new, queue, - default_rm_virus, -1, -1 + "There is no virus to remove here!" ); - -} -template -inline void Agent::rm_virus( - VirusPtr & virus, - Model * model, - epiworld_fast_int state_new, - epiworld_fast_int queue -) -{ - - if (virus->agent != this) - throw std::logic_error("Cannot remove a virus from another agent!"); + CHECK_COALESCE_(state_new, virus->state_post, state); + CHECK_COALESCE_(queue, virus->queue_post, Queue::Everyone); model->actions_add( this, virus, nullptr, nullptr, state_new, queue, default_rm_virus, -1, -1 ); - - + } template @@ -390,9 +332,12 @@ 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); + model->actions_add( this, nullptr, nullptr, model->entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + default_rm_entity, entities_locations[entity_idx], entity_idx ); } @@ -419,82 +364,30 @@ inline void Agent::rm_entity( entity.get_name() + "\"." ); + CHECK_COALESCE_(state_new, entity.state_post, state); + CHECK_COALESCE_(queue, entity.queue_post, Queue::NoOne); model->actions_add( this, nullptr, nullptr, entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + default_rm_entity, entities_locations[entity_idx], entity_idx ); } template inline void Agent::rm_agent_by_virus( - epiworld_fast_uint virus_idx, Model * model, epiworld_fast_int state_new, epiworld_fast_int queue ) { - if (state_new == -99) - state_new = state; - - if (virus_idx >= n_viruses) - throw std::range_error( - std::string("The virus trying to remove the agent is out of range. ") + - std::string("This agent has only ") + std::to_string(n_viruses) + - std::string(" and you are trying to remove virus # ") + - std::to_string(virus_idx) + std::string(".") - ); + CHECK_COALESCE_(state_new, virus->state_removed, state); + CHECK_COALESCE_(queue, virus->queue_removed, Queue::Everyone); - // Removing viruses - for (size_t i = 0u; i < n_viruses; ++i) - { - if (i != virus_idx) - rm_virus(i, model); - } - - // Changing state to new_state - VirusPtr & v = viruses[virus_idx]; - epiworld_fast_int dead_state, dead_queue; - v->get_state(nullptr, nullptr, &dead_state); - v->get_queue(nullptr, nullptr, &dead_queue); - - if (queue != -99) - dead_queue = queue; - - change_state( - model, - // Either preserve the current state or apply a new one - (dead_state < 0) ? state : static_cast(dead_state), - - // By default, it will be removed from the queue... unless the user - // says the contrary! - (dead_queue == -99) ? Queue::NoOne : dead_queue - ); - -} - -template -inline void Agent::rm_agent_by_virus( - VirusPtr & virus, - Model * model, - epiworld_fast_int state_new, - epiworld_fast_int queue -) -{ - - if (virus->get_agent() == nullptr) - throw std::logic_error("The virus trying to remove the agent has no host."); - - if (virus->get_agent()->id != id) - throw std::logic_error("Viruses can only remove their hosts'."); - - rm_agent_by_virus( - virus->pos_in_agent, - model, - state_new, - queue - ); + model->actions_add( + this, virus, nullptr, nullptr, state_new, queue, + default_rm_virus, -1, -1 + ); } @@ -538,29 +431,10 @@ inline int Agent::get_id() const } template -inline Viruses Agent::get_viruses() { - - return Viruses(*this); - -} - -template -inline const Viruses_const Agent::get_viruses() const { - - return Viruses_const(*this); - +inline VirusPtr & Agent::get_virus() { + return virus; } -template -inline VirusPtr & Agent::get_virus(int i) { - return viruses.at(i); -} - -template -inline size_t Agent::get_n_viruses() const noexcept -{ - return n_viruses; -} template inline Tools Agent::get_tools() { @@ -588,8 +462,7 @@ template inline void Agent::mutate_virus() { - for (auto & v : viruses) - v->mutate(); + virus->mutate(); } @@ -715,7 +588,8 @@ inline void Agent::change_state( { model->actions_add( - this, nullptr, nullptr, nullptr, new_state, queue, nullptr, -1, -1 + this, nullptr, nullptr, nullptr, new_state, queue, + default_change_state, -1, -1 ); return; @@ -731,8 +605,7 @@ template inline void Agent::reset() { - this->viruses.clear(); - n_viruses = 0u; + this->virus = nullptr; this->tools.clear(); n_tools = 0u; @@ -779,9 +652,8 @@ inline bool Agent::has_tool(const Tool & tool) const template inline bool Agent::has_virus(epiworld_fast_uint t) const { - for (auto & v : viruses) - if (v->get_id() == static_cast(t)) - return true; + if (virus->get_id() == static_cast(t)) + return true; return false; } @@ -790,9 +662,8 @@ template inline bool Agent::has_virus(std::string name) const { - for (auto & v : viruses) - if (v->get_name() == name) - return true; + if (virus->get_name() == name) + return true; return false; @@ -816,15 +687,18 @@ inline void Agent::print( if (compressed) { printf_epiworld( - "Agent: %i, state: %s (%lu), Nvirus: %lu, NTools: %lu, NNeigh: %lu\n", - id, model->states_labels[state].c_str(), state, n_viruses, n_tools, neighbors.size() + "Agent: %i, state: %s (%lu), Has virus: %s, NTools: %lu, NNeigh: %lu\n", + id, model->states_labels[state].c_str(), state, + virus == nullptr ? std::string("no").c_str() : std::string("yes").c_str(), + n_tools, neighbors.size() ); } else { printf_epiworld("Information about agent id %i\n", this->id); printf_epiworld(" State : %s (%lu)\n", model->states_labels[state].c_str(), state); - printf_epiworld(" Virus count : %lu\n", n_viruses); + printf_epiworld(" Has virus : %s\n", virus == nullptr ? + std::string("no").c_str() : std::string("yes").c_str()); printf_epiworld(" Tool count : %lu\n", n_tools); printf_epiworld(" Neigh. count : %lu\n", neighbors.size()); @@ -975,24 +849,21 @@ inline bool Agent::operator==(const Agent & other) const // state_last_changed != other.state_last_changed, // "Agent:: state_last_changed don't match" // ) ///< Last time the agent was updated. - - + EPI_DEBUG_FAIL_AT_TRUE( - n_viruses != other.n_viruses, - "Agent:: n_viruses don't match" - ) - + ((virus == nullptr) && (other.virus != nullptr)) || + ((virus != nullptr) && (other.virus == nullptr)), + "Agent:: virus don't match" + ) - for (size_t i = 0u; i < n_viruses; ++i) + if ((virus != nullptr) && (other.virus != nullptr)) { - EPI_DEBUG_FAIL_AT_TRUE( - *viruses[i] != *other.viruses[i], - "Agent:: viruses[i] don't match" + *virus != *other.virus, + "Agent:: virus doesn't match" ) - } - + EPI_DEBUG_FAIL_AT_TRUE(n_tools != other.n_tools, "Agent:: n_tools don't match") for (size_t i = 0u; i < n_tools; ++i) diff --git a/inst/include/epiworld/agentssample-bones.hpp b/inst/include/epiworld/agentssample-bones.hpp index 165c4985..c8462c58 100644 --- a/inst/include/epiworld/agentssample-bones.hpp +++ b/inst/include/epiworld/agentssample-bones.hpp @@ -41,6 +41,7 @@ class AgentsSample { Agent * agent = nullptr; int sample_type = SAMPLETYPE::AGENT; + std::vector< size_t > states = {}; void sample_n(size_t n); ///< Backbone function for sampling @@ -52,9 +53,23 @@ class AgentsSample { AgentsSample(const AgentsSample & a) = delete; ///< Copy constructor AgentsSample(AgentsSample && a) = delete; ///< Move constructor - AgentsSample(Model & model_, size_t n, bool truncate = false); - AgentsSample(Model * model, Entity & entity_, size_t n, bool truncate = false); - AgentsSample(Model * model, Agent & agent_, size_t n, bool truncate = false); + AgentsSample( + Model & model_, size_t n, + std::vector< size_t > states_ = {}, + bool truncate = false + ); + + AgentsSample( + Model * model, Entity & entity_, size_t n, + std::vector< size_t > states_ = {}, + bool truncate = false + ); + + AgentsSample( + Model * model, Agent & agent_, size_t n, + std::vector< size_t > states_ = {}, + bool truncate = false + ); ~AgentsSample(); @@ -71,9 +86,12 @@ template inline AgentsSample::AgentsSample( Model & model_, size_t n, + std::vector< size_t > states_, bool truncate ) { + states = states_; + if (truncate) { @@ -85,12 +103,12 @@ inline AgentsSample::AgentsSample( "There are only " + std::to_string(model_.size()) + " agents. You cannot " + "sample " + std::to_string(n)); - sample_size = n; - model = &model_; - sample_type = SAMPLETYPE::MODEL; + sample_size = n; + model = &model_; + sample_type = SAMPLETYPE::MODEL; - agents = &model_.sampled_population; - agents_n = &model_.sampled_population_n; + agents = &model_.sampled_population; + agents_n = &model_.sampled_population_n; agents_left = &model_.population_left; agents_left_n = &model_.population_left_n; @@ -105,9 +123,13 @@ template inline AgentsSample::AgentsSample( Model * model, Entity & entity_, - size_t n, bool truncate + size_t n, + std::vector< size_t > states_, + bool truncate ) { + states = states_; + if (truncate) { @@ -119,12 +141,12 @@ inline AgentsSample::AgentsSample( "There are only " + std::to_string(entity_.size()) + " agents. You cannot " + "sample " + std::to_string(n)); - sample_size = n; - model = &entity_.model; - sample_type = SAMPLETYPE::ENTITY; + sample_size = n; + model = &entity_.model; + sample_type = SAMPLETYPE::ENTITY; - agents = &entity_.sampled_agents; - agents_n = &entity_.sampled_agents_n; + agents = &entity_.sampled_agents; + agents_n = &entity_.sampled_agents_n; agents_left = &entity_.sampled_agents_left; agents_left_n = &entity_.sampled_agents_left_n; @@ -152,16 +174,18 @@ inline AgentsSample::AgentsSample( Model * model, Agent & agent_, size_t n, + std::vector< size_t > states_, bool truncate ) { + states = states_; sample_type = SAMPLETYPE::AGENT; - agent = &agent_; + agent = &agent_; - agents = &agent_.sampled_agents; - agents_n = &agent_.sampled_agents_n; + agents = &agent_.sampled_agents; + agents_n = &agent_.sampled_agents_n; agents_left = &agent_.sampled_agents_left; agents_left_n = &agent_.sampled_agents_left_n; @@ -170,7 +194,7 @@ inline AgentsSample::AgentsSample( size_t agents_in_entities = 0; Entities entities_a = agent->get_entities(); - std::vector< int > cum_agents_count(entities_a.size(), 0); + std::vector< size_t > cum_agents_count(entities_a.size(), 0); int idx = -1; for (auto & e : entities_a) { @@ -193,15 +217,18 @@ inline AgentsSample::AgentsSample( } else if (n > agents_in_entities) throw std::logic_error( - "There are only " + std::to_string(agents_in_entities) + " agents. You cannot " + - "sample " + std::to_string(n)); + "There are only " + std::to_string(agents_in_entities) + + " agents. You cannot " + + "sample " + std::to_string(n) + ); sample_size = n; if (agents->size() < n) agents->resize(n); - for (size_t i = 0u; i < n; ++i) + size_t i_obs = 0u; + for (size_t i = 0u; i < agents_in_entities; ++i) { int jth = std::floor(model->runif() * agents_in_entities); for (size_t e = 0u; e < cum_agents_count.size(); ++e) @@ -209,11 +236,24 @@ inline AgentsSample::AgentsSample( // Are we in the limit? if (jth <= cum_agents_count[e]) { + size_t agent_idx = 0u; if (e == 0) // From the first group - agents->operator[](i) = entities_a[e][jth]; + agent_idx = entities_a[e][jth]; else - agents->operator[](i) = entities_a[e][jth - cum_agents_count[e - 1]]; + agent_idx = entities_a[e][jth - cum_agents_count[e - 1]]; + + // Getting the state + size_t state = model->population[agent_idx].get_state(); + + // Checking if states was specified + if (states.size()) + { + if (std::find(states.begin(), states.end(), state) != states.end()) + continue; + } + agents->operator[](i_obs++) = agent_idx; + break; } @@ -275,59 +315,143 @@ template inline void AgentsSample::sample_n(size_t n) { - // Checking if the size of the entity has changed (or hasn't been initialized) - if (sample_type == SAMPLETYPE::MODEL) + // Reducing size + if (states.size()) { + + // Getting the number of agents left + agents_left->clear(); - if (model->size() != agents_left->size()) + if (sample_type == SAMPLETYPE::MODEL) { - agents_left->resize(model->size(), 0u); - std::iota(agents_left->begin(), agents_left->end(), 0u); - } + // Making some room + agents_left->reserve(model->size()); + + // Iterating through the agents in the population + for (size_t a_i = 0u; a_i < model->population.size(); ++a_i) + { - } else if (sample_type == SAMPLETYPE::ENTITY) { + // If the agent is within the selected set of states, + // then we add it to the list of agents left + size_t s = model->population[a_i].get_state(); + if (std::find(states.begin(), states.end(), s) != states.end()) + agents_left->push_back(a_i); - if (entity->size() != agents_left->size()) + } + + } + else if (sample_type == SAMPLETYPE::ENTITY) { - agents_left->resize(entity->size(), 0u); - std::iota(agents_left->begin(), agents_left->end(), 0u); + // Making room + agents_left->reserve(entity->size()); + + // Iterating through the agents in the entity + for (size_t a_i = 0u; a_i < entity->size(); ++a_i) + { + size_t s = model->population[entity->agents[a_i]].get_state(); + if (std::find(states.begin(), states.end(), s) != states.end()) + agents_left->push_back(a_i); + + } } - } + } else { + + // Checking if the size of the entity has changed (or hasn't been initialized) + if (sample_type == SAMPLETYPE::MODEL) + { + + if (model->size() != agents_left->size()) + { + agents_left->resize(model->size(), 0u); + std::iota(agents_left->begin(), agents_left->end(), 0u); + } + + } else if (sample_type == SAMPLETYPE::ENTITY) { + + if (entity->size() != agents_left->size()) + { + + agents_left->resize(entity->size(), 0u); + std::iota(agents_left->begin(), agents_left->end(), 0u); + + } + + } + + } // Restart the counter of agents left *agents_left_n = agents_left->size(); + // Making sure we have enough room for the sample of agents if (agents->size() < sample_size) agents->resize(sample_size, nullptr); if (sample_type == SAMPLETYPE::MODEL) { + #ifdef EPI_DEBUG + std::vector< bool > __sampled(model->size(), true); + for (auto & a_i: *agents_left) + __sampled[a_i] = false; + #endif + for (size_t i = 0u; i < n; ++i) { - size_t ith = agents_left->operator[](model->runif() * ((*agents_left_n)--)); + // Sampling from 0 to (agents_left_n - 1) + size_t ith_ = static_cast(model->runif() * ((*agents_left_n)--)); + + // Getting the id of the agent and adding it to the list of agents + size_t ith = agents_left->operator[](ith_); agents->operator[](i) = &model->population[ith]; + #ifdef EPI_DEBUG + if (__sampled[ith]) + throw std::logic_error("The same agent was sampled twice."); + else + __sampled[ith] = true; + #endif + // Updating list - std::swap(agents_left->operator[](ith), agents_left->operator[](*agents_left_n)); + std::swap( + agents_left->operator[](ith_), + agents_left->operator[](*agents_left_n) + ); } - } else if (sample_type == SAMPLETYPE::ENTITY) { + + } + else if (sample_type == SAMPLETYPE::ENTITY) + { + + #ifdef EPI_DEBUG + std::vector< bool > __sampled(entity->size(), true); + for (auto & a_i: *agents_left) + __sampled[a_i] = false; + #endif for (size_t i = 0u; i < n; ++i) { - size_t ith = agents_left->operator[](model->runif() * (--(*agents_left_n))); - agents->operator[](i) = entity->agents[ith]; + size_t ith_ = static_cast(model->runif() * ((*agents_left_n)--)); + size_t ith = agents_left->operator[](ith_); + agents->operator[](i) = &model->population[entity->agents[ith]]; + + #ifdef EPI_DEBUG + if (__sampled[ith]) + throw std::logic_error("The same agent was sampled twice."); + else + __sampled[ith] = true; + #endif // Updating list - std::swap(agents_left->operator[](ith), agents_left->operator[](*agents_left_n)); + std::swap(agents_left->operator[](ith_), agents_left->operator[](*agents_left_n)); } diff --git a/inst/include/epiworld/config.hpp b/inst/include/epiworld/config.hpp index 8c65e067..633dd0fd 100644 --- a/inst/include/epiworld/config.hpp +++ b/inst/include/epiworld/config.hpp @@ -9,7 +9,7 @@ #define EPIWORLD_MAXNEIGHBORS 1048576 #endif -#ifdef _OPENMP +#if defined(_OPENMP) || defined(__OPENMP) #include // #else // #define omp_get_thread_num() 0 @@ -259,7 +259,9 @@ struct Action { if (a) \ {\ throw EPI_DEBUG_ERROR(std::logic_error, b); \ - } + } + + #define epiexception(a) std::logic_error #else #define EPI_DEBUG_PRINTF(fmt, ...) #define EPI_DEBUG_ERROR(fmt, ...) @@ -271,9 +273,10 @@ struct Action { #define EPI_DEBUG_FAIL_AT_TRUE(a, b) \ if (a) \ return false; + #define epiexception(a) a #endif -#ifdef EPI_DEBUG_NO_THREAD_ID +#if defined(EPI_DEBUG_NO_THREAD_ID) || (!defined(__OPENMP) && !defined(_OPENMP)) #define EPI_GET_THREAD_ID() 0 #else #define EPI_GET_THREAD_ID() omp_get_thread_num() diff --git a/inst/include/epiworld/database-bones.hpp b/inst/include/epiworld/database-bones.hpp index be62c671..a57f37a1 100644 --- a/inst/include/epiworld/database-bones.hpp +++ b/inst/include/epiworld/database-bones.hpp @@ -22,6 +22,9 @@ inline void default_rm_virus(Action & a, Model * m); template inline void default_rm_tool(Action & a, Model * m); +template +inline void default_change_state(Action & a, Model * m); + /** * @brief Statistical data about the process * @@ -34,6 +37,7 @@ class DataBase { friend void default_add_tool(Action & a, Model * m); friend void default_rm_virus(Action & a, Model * m); friend void default_rm_tool(Action & a, Model * m); + friend void default_change_state(Action & a, Model * m); private: Model * model; diff --git a/inst/include/epiworld/entity-meat.hpp b/inst/include/epiworld/entity-meat.hpp index 3220848e..6912ec52 100644 --- a/inst/include/epiworld/entity-meat.hpp +++ b/inst/include/epiworld/entity-meat.hpp @@ -1,28 +1,6 @@ #ifndef EPIWORLD_ENTITY_MEAT_HPP #define EPIWORLD_ENTITY_MEAT_HPP -// template -// inline Entity::Entity(const Entity & e) : -// model(e.model), -// id(e.id), -// agents(0u), -// agents_location(0u), -// n_agents(0), -// sampled_agents(0u), -// sampled_agents_n(0u), -// sampled_agents_left(0u), -// sampled_agents_left_n(0u), -// max_capacity(e.max_capacity), -// entity_name(e.entity_name), -// location(e.location), -// state_init(e.state_init), -// state_post(e.state_post), -// queue_init(e.queue_init), -// queue_post(e.queue_post) -// { - -// } - template inline void Entity::add_agent( Agent & p, diff --git a/inst/include/epiworld/epiworld.hpp b/inst/include/epiworld/epiworld.hpp index 3644e5df..583878ff 100644 --- a/inst/include/epiworld/epiworld.hpp +++ b/inst/include/epiworld/epiworld.hpp @@ -71,8 +71,6 @@ namespace epiworld { #include "models/models.hpp" - - } #endif \ No newline at end of file diff --git a/inst/include/epiworld/model-bones.hpp b/inst/include/epiworld/model-bones.hpp index b913a383..01efeea3 100644 --- a/inst/include/epiworld/model-bones.hpp +++ b/inst/include/epiworld/model-bones.hpp @@ -104,7 +104,6 @@ class Model { bool using_backup = true; std::vector< Agent > population_backup = {}; - /** * @name Auxiliary variables for AgentsSample iterators * @@ -170,8 +169,13 @@ class Model { epiworld_fast_uint ndays = 0; Progress pb; - std::vector< UpdateFun > state_fun = {}; - std::vector< std::string > states_labels = {}; + std::vector< UpdateFun > state_fun = {}; ///< Functions to update states + std::vector< std::string > states_labels = {}; ///< Labels of the states + + /** Function to distribute states. Goes along with the function */ + std::function*)> initial_states_fun = [](Model * /**/) + -> void {}; + epiworld_fast_uint nstates = 0u; bool verbose = true; @@ -221,20 +225,13 @@ class Model { VirusPtr virus_, ToolPtr tool_, Entity * entity_, - epiworld_fast_uint new_state_, + epiworld_fast_int new_state_, epiworld_fast_int queue_, ActionFun call_, int idx_agent_, int idx_object_ ); - /** - * @brief Executes the stored action - * - * @param model_ Model over which it will be executed. - */ - void actions_run(); - /** * @name Tool Mixers * @@ -258,12 +255,13 @@ class Model { public: + std::vector array_double_tmp; std::vector * > array_virus_tmp; Model(); Model(const Model & m); - Model(Model & m) = delete; + Model(Model & m); Model(Model && m); Model & operator=(const Model & m); @@ -401,7 +399,7 @@ class Model { std::vector< Entity > & get_entities(); - void agents_smallworld( + Model & agents_smallworld( epiworld_fast_uint n = 1000, epiworld_fast_uint k = 5, bool d = false, @@ -423,7 +421,7 @@ class Model { void update_state(); void mutate_virus(); void next(); - virtual void run( + virtual Model & run( epiworld_fast_uint ndays, int seed = -1 ); ///< Runs the simulation (after initialization) @@ -438,14 +436,14 @@ class Model { ); ///@} - size_t get_n_viruses() const; - size_t get_n_tools() const; + size_t get_n_viruses() const; ///< Number of viruses in the model + size_t get_n_tools() const; ///< Number of tools in the model epiworld_fast_uint get_ndays() const; epiworld_fast_uint get_n_replicates() const; void set_ndays(epiworld_fast_uint ndays); bool get_verbose() const; - void verbose_off(); - void verbose_on(); + Model & verbose_off(); + Model & verbose_on(); int today() const; ///< The current time of the model /** @@ -525,7 +523,7 @@ class Model { * */ virtual void reset(); - void print(bool lite = false) const; + const Model & print(bool lite = false) const; Model && clone() const; @@ -550,6 +548,19 @@ class Model { void print_state_codes() const; ///@} + /** + * @name Initial states + * + * @details These functions are called before the simulation starts. + * + * @param proportions_ Vector of proportions for each state. + * @param queue_ Vector of queue for each state. + */ + virtual Model & initial_states( + std::vector< double > /*proportions_*/, + std::vector< int > /*queue_*/ + ) {return *this;}; + /** * @name Setting and accessing parameters from the model * @@ -655,7 +666,7 @@ class Model { */ ////@{ void queuing_on(); ///< Activates the queuing system (default.) - void queuing_off(); ///< Deactivates the queuing system. + Model & queuing_off(); ///< Deactivates the queuing system. bool is_queuing_on() const; ///< Query if the queuing system is on. Queue & get_queue(); ///< Retrieve the `Queue` object. ///@} @@ -674,6 +685,8 @@ 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); @@ -704,6 +717,14 @@ class Model { bool operator==(const Model & other) const; bool operator!=(const Model & other) const {return !operator==(other);}; + /** + * @brief Executes the stored action + * + * @param model_ Model over which it will be executed. + */ + void actions_run(); + + }; #endif diff --git a/inst/include/epiworld/model-meat-print.hpp b/inst/include/epiworld/model-meat-print.hpp index 673f5189..a5a419aa 100644 --- a/inst/include/epiworld/model-meat-print.hpp +++ b/inst/include/epiworld/model-meat-print.hpp @@ -2,7 +2,7 @@ #define EPIWORLD_MODEL_MEAT_PRINT_HPP template -inline void Model::print(bool lite) const +inline const Model & Model::print(bool lite) const { // Horizontal line @@ -58,7 +58,7 @@ inline void Model::print(bool lite) const printf_epiworld(" The model hasn't been run yet.\n"); } - return; + return *this; } printf_epiworld("%s\n%s\n\n",line.c_str(), "SIMULATION STUDY"); @@ -324,7 +324,7 @@ inline void Model::print(bool lite) const if (today() != 0) (void) db.transition_probability(true); - return; + return *this; } diff --git a/inst/include/epiworld/model-meat.hpp b/inst/include/epiworld/model-meat.hpp index 431a58c2..ce7e59b8 100644 --- a/inst/include/epiworld/model-meat.hpp +++ b/inst/include/epiworld/model-meat.hpp @@ -1,8 +1,6 @@ #ifndef EPIWORLD_MODEL_MEAT_HPP #define EPIWORLD_MODEL_MEAT_HPP - - /** * @brief Function factory for saving model runs * @@ -150,38 +148,31 @@ inline std::function*)> make_save_run( return saver; } + template inline void Model::actions_add( Agent * agent_, VirusPtr virus_, ToolPtr tool_, Entity * entity_, - epiworld_fast_uint new_state_, + epiworld_fast_int new_state_, epiworld_fast_int queue_, ActionFun call_, int idx_agent_, int idx_object_ ) { - + ++nactions; #ifdef EPI_DEBUG if (nactions == 0) throw std::logic_error("Actions cannot be zero!!"); - - if ((virus_ != nullptr) && idx_agent_ >= 0) - { - if (idx_agent_ >= static_cast(virus_->get_agent()->get_n_viruses())) - throw std::logic_error( - "The virus to add is out of range in the host agent." - ); - } #endif if (nactions > actions.size()) { - actions.push_back( + actions.emplace_back( Action( agent_, virus_, tool_, entity_, new_state_, queue_, call_, idx_agent_, idx_object_ @@ -197,7 +188,7 @@ inline void Model::actions_add( A.virus = virus_; A.tool = tool_; A.entity = entity_; - A.new_state = new_state_; + A.new_state = new_state_; A.queue = queue_; A.call = call_; A.idx_agent = idx_agent_; @@ -213,76 +204,46 @@ template inline void Model::actions_run() { // Making the call - while (nactions != 0u) + size_t nactions_tmp = 0; + while (nactions_tmp < nactions) { - Action a = actions[--nactions]; + Action & a = actions[nactions_tmp++]; Agent * p = a.agent; - // Applying function - if (a.call) - { - a.call(a, this); - } - - // Updating state - if (static_cast(p->state) != a.new_state) - { - - if (a.new_state >= static_cast(nstates)) - throw std::range_error( - "The proposed state " + std::to_string(a.new_state) + " is out of range. " + - "The model currently has " + std::to_string(nstates - 1) + " states."); - - // Figuring out if we need to undo a change - // If the agent has made a change in the state recently, then we - // need to undo the accounting, e.g., if A->B was made, we need to - // undo it and set B->A so that the daily accounting is right. - if (p->state_last_changed == today()) - { - - // Updating accounting - db.update_state(p->state_prev, p->state, true); // Undoing - db.update_state(p->state_prev, a.new_state); - - for (size_t v = 0u; v < p->n_viruses; ++v) - { - db.update_virus(p->viruses[v]->id, p->state, p->state_prev); // Undoing - db.update_virus(p->viruses[v]->id, p->state_prev, a.new_state); - } - - for (size_t t = 0u; t < p->n_tools; ++t) - { - db.update_tool(p->tools[t]->id, p->state, p->state_prev); // Undoing - db.update_tool(p->tools[t]->id, p->state_prev, a.new_state); - } - - // Changing to the new state, we won't update the - // previous state in case we need to undo the change - p->state = a.new_state; - - } else { - - // Updating accounting - db.update_state(p->state, a.new_state); - - for (size_t v = 0u; v < p->n_viruses; ++v) - db.update_virus(p->viruses[v]->id, p->state, a.new_state); + #ifdef EPI_DEBUG + if (a.new_state >= static_cast(nstates)) + throw std::range_error( + "The proposed state " + std::to_string(a.new_state) + " is out of range. " + + "The model currently has " + std::to_string(nstates - 1) + " states."); - for (size_t t = 0u; t < p->n_tools; ++t) - db.update_tool(p->tools[t]->id, p->state, a.new_state); + if (a.new_state < 0) + throw std::range_error( + "The proposed state " + std::to_string(a.new_state) + " is out of range. " + + "The state cannot be negative."); + #endif - // Saving the last state and setting the new one - p->state_prev = p->state; - p->state = a.new_state; + // Undoing the change in the transition matrix + if ((p->state_last_changed == today()) && (static_cast(p->state) != a.new_state)) + { + // Undoing state change in the transition matrix + // The previous state is already recorded + db.update_state(p->state_prev, p->state, true); - // It used to be a day before, but we still - p->state_last_changed = today(); + } else + p->state_prev = p->state; // Recording the previous state - } - + // Applying function after the fact. This way, if there were + // updates, they can be recorded properly, before losing the information + p->state = a.new_state; + if (a.call) + { + a.call(a, this); } + // Registering that the last change was today + p->state_last_changed = today(); + #ifdef EPI_DEBUG if (static_cast(p->state) >= static_cast(nstates)) throw std::range_error( @@ -311,6 +272,9 @@ inline void Model::actions_run() } + // Go back to square 1 + nactions = 0u; + return; } @@ -435,6 +399,7 @@ inline Model::Model(const Model & model) : pb(model.pb), state_fun(model.state_fun), states_labels(model.states_labels), + initial_states_fun(model.initial_states_fun), nstates(model.nstates), verbose(model.verbose), current_date(model.current_date), @@ -475,6 +440,10 @@ inline Model::Model(const Model & model) : } +template +inline Model::Model(Model & model) : + Model(dynamic_cast< const Model & >(model)) {} + template inline Model::Model(Model && model) : name(std::move(model.name)), @@ -515,6 +484,7 @@ inline Model::Model(Model && model) : pb(std::move(model.pb)), state_fun(std::move(model.state_fun)), states_labels(std::move(model.states_labels)), + initial_states_fun(std::move(model.initial_states_fun)), nstates(model.nstates), verbose(model.verbose), current_date(std::move(model.current_date)), @@ -586,6 +556,7 @@ inline Model & Model::operator=(const Model & m) state_fun = m.state_fun; states_labels = m.states_labels; + initial_states_fun = m.initial_states_fun; nstates = m.nstates; verbose = m.verbose; @@ -648,7 +619,7 @@ inline std::vector< Viruses_const > Model::get_agents_viruses() cons std::vector< Viruses_const > viruses(population.size()); for (size_t i = 0u; i < population.size(); ++i) - viruses[i] = population[i].get_viruses(); + viruses[i] = population[i].get_virus(); return viruses; @@ -661,7 +632,7 @@ inline std::vector< Viruses > Model::get_agents_viruses() std::vector< Viruses > viruses(population.size()); for (size_t i = 0u; i < population.size(); ++i) - viruses[i] = population[i].get_viruses(); + viruses[i] = population[i].get_virus(); return viruses; @@ -674,7 +645,7 @@ inline std::vector> & Model::get_entities() } template -inline void Model::agents_smallworld( +inline Model & Model::agents_smallworld( epiworld_fast_uint n, epiworld_fast_uint k, bool d, @@ -684,6 +655,8 @@ inline void Model::agents_smallworld( agents_from_adjlist( rgraph_smallworld(n, k, p, d, *this) ); + + return *this; } template @@ -770,10 +743,9 @@ inline void Model::dist_virus() // Starting first infection int n = size(); - std::vector< size_t > idx(n); - - int n_left = n; + std::vector< size_t > idx(n, 0u); std::iota(idx.begin(), idx.end(), 0); + int n_left = idx.size(); for (size_t v = 0u; v < viruses.size(); ++v) { @@ -811,7 +783,7 @@ inline void Model::dist_virus() Agent & agent = population[idx[loc]]; // Adding action - agent.add_virus( + agent.set_virus( virus, const_cast * >(this), virus->state_init, @@ -1488,7 +1460,7 @@ inline void Model::next() { } template -inline void Model::run( +inline Model & Model::run( epiworld_fast_uint ndays, int seed ) @@ -1594,6 +1566,8 @@ inline void Model::run( chrono_end(); + return *this; + } template @@ -1824,6 +1798,15 @@ inline void Model::update_state() { template inline void Model::mutate_virus() { + // Checking if any virus has mutation + size_t nmutates = 0u; + for (const auto & v: viruses) + if (v->mutation_fun) + nmutates++; + + if (nmutates == 0u) + return; + if (use_queuing) { @@ -1834,9 +1817,8 @@ inline void Model::mutate_virus() { if (queue[++i] == 0) continue; - if (p.n_viruses > 0u) - for (auto & v : p.get_viruses()) - v->mutate(this); + if (p.virus != nullptr) + p.virus->mutate(this); } @@ -1847,9 +1829,8 @@ inline void Model::mutate_virus() { for (auto & p: population) { - if (p.n_viruses > 0u) - for (auto & v : p.get_viruses()) - v->mutate(this); + if (p.virus != nullptr) + p.virus->mutate(this); } @@ -1890,13 +1871,15 @@ inline bool Model::get_verbose() const { } template -inline void Model::verbose_on() { +inline Model & Model::verbose_on() { verbose = true; + return *this; } template -inline void Model::verbose_off() { +inline Model & Model::verbose_off() { verbose = false; + return *this; } template @@ -2105,6 +2088,9 @@ inline void Model::reset() { dist_virus(); dist_tools(); + // Distributing initial state, if specified + initial_states_fun(this); + // Recording the original state (at time 0) and advancing // to time 1 next(); @@ -2531,9 +2517,10 @@ inline void Model::queuing_on() } template -inline void Model::queuing_off() +inline Model & Model::queuing_off() { use_queuing = false; + return *this; } template @@ -2554,6 +2541,18 @@ inline const std::vector< VirusPtr > & Model::get_viruses() const return viruses; } +template +inline const std::vector< epiworld_double > & Model::get_prevalence_virus() const +{ + return prevalence_virus; +} + +template +inline const std::vector< bool > & Model::get_prevalence_virus_as_proportion() const +{ + return prevalence_virus_as_proportion; +} + template const std::vector< ToolPtr > & Model::get_tools() const { diff --git a/inst/include/epiworld/models/diffnet.hpp b/inst/include/epiworld/models/diffnet.hpp index ee9a9967..b7ab0724 100644 --- a/inst/include/epiworld/models/diffnet.hpp +++ b/inst/include/epiworld/models/diffnet.hpp @@ -91,24 +91,25 @@ inline ModelDiffNet::ModelDiffNet( if (neighbor->get_state() == ModelDiffNet::ADOPTER) { - for (const VirusPtr & v : neighbor->get_viruses()) - { - - /* And it is a function of susceptibility_reduction as well */ - double p_i = - (1.0 - agent.get_susceptibility_reduction(v, m)) * - (1.0 - agent.get_transmission_reduction(v, m)) - ; + auto & v = neighbor->get_virus(); - size_t vid = v->get_id(); - if (!stored[vid]) - { - stored[vid] = true; - innovations[vid] = &(*v); - } - exposure[vid] += p_i; - - } + if (v == nullptr) + continue; + + /* And it is a function of susceptibility_reduction as well */ + double p_i = + (1.0 - agent.get_susceptibility_reduction(v, m)) * + (1.0 - agent.get_transmission_reduction(v, m)) + ; + + size_t vid = v->get_id(); + if (!stored[vid]) + { + stored[vid] = true; + innovations[vid] = &(*v); + } + exposure[vid] += p_i; + } @@ -141,7 +142,7 @@ inline ModelDiffNet::ModelDiffNet( return; // Otherwise, it is adopted from any of the neighbors - agent.add_virus( + agent.set_virus( *innovations.at(which), m, ModelDiffNet::ADOPTER diff --git a/inst/include/epiworld/models/globalactions.hpp b/inst/include/epiworld/models/globalactions.hpp index f0e20fe1..432928d4 100644 --- a/inst/include/epiworld/models/globalactions.hpp +++ b/inst/include/epiworld/models/globalactions.hpp @@ -81,7 +81,9 @@ inline std::function*)> globalaction_tool_logit( // Computing the probability using a logit. Uses OpenMP reduction // to sum the coefficients. double p = 0.0; + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp parallel for reduction(+:p) + #endif for (size_t i = 0u; i < coefs.size(); ++i) p += coefs.at(i) * agent(vars[i]); diff --git a/inst/include/epiworld/models/init-functions.hpp b/inst/include/epiworld/models/init-functions.hpp new file mode 100644 index 00000000..660d5db8 --- /dev/null +++ b/inst/include/epiworld/models/init-functions.hpp @@ -0,0 +1,341 @@ +#ifndef EPIWORLD_MODELS_INIT_FUNCTIONS_HPP +#define EPIWORLD_MODELS_INIT_FUNCTIONS_HPP + +/** + * @brief Creates an initial function for the SIR-like models + * The function is used for the initial states of the model. +*/ +template +inline std::function*)> create_init_function_sir( + std::vector< double > proportions_ +) { + + // Checking widths + if (proportions_.size() != 1u) + throw std::invalid_argument( + "The vector of proportions must have a single element." + ); + + // Proportion should be within [0, 1] + if ((proportions_[0] < 0.0) || (proportions_[0] > 1.0)) + throw std::invalid_argument( + "The proportion must be within (0, 1)." + ); + + double prop = proportions_[0u]; + + std::function*)> fun = + [prop] (epiworld::Model * model) -> void { + + // Figuring out information about the viruses + double tot = 0.0; + const auto & vpreval = model->get_prevalence_virus(); + const auto & vprop = model->get_prevalence_virus_as_proportion(); + double n = static_cast(model->size()); + for (size_t i = 0u; i < model->get_n_viruses(); ++i) + { + if (vprop[i]) + tot += vpreval[i]; + else + tot += vpreval[i] / n; + } + + // Putting the total into context + double tot_left = 1.0 - tot; + + // Since susceptible and infected are "fixed," + // we only need to change recovered + size_t nrecovered = prop * tot_left * n; + + epiworld::AgentsSample sample( + *model, + nrecovered, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample) + agent->change_state(model, 2, Queue::NoOne); + + // Running the actions + model->actions_run(); + + return; + + }; + + return fun; + +} + +/** + * @brief Creates an initial function for the SIR-like models + * The function is used for the initial states of the model. +*/ +template +inline std::function*)> create_init_function_sird( + std::vector< double > prop +) { + + // Check length of prop equals two + if (prop.size() != 2u) + throw std::invalid_argument( + "The vector of proportions must have two elements." + ); + + // Check elements in prop are within [0, 1] and sum up to 1 + double tot = 0.0; + for (auto & v : prop) + { + if ((v < 0.0) || (v > 1.0)) + throw std::invalid_argument( + "The proportion must be within (0, 1)." + ); + tot += v; + } + + if (tot >= 1.0) + throw std::invalid_argument( + "The proportions must sum up to 1." + ); + + std::function*)> fun = + [prop] (epiworld::Model * model) -> void { + + // Figuring out information about the viruses + double tot = 0.0; + const auto & vpreval = model->get_prevalence_virus(); + const auto & vprop = model->get_prevalence_virus_as_proportion(); + double n = static_cast(model->size()); + for (size_t i = 0u; i < model->get_n_viruses(); ++i) + { + if (vprop[i]) + tot += vpreval[i]; + else + tot += vpreval[i] / n; + } + + // Putting the total into context + double tot_left = 1.0 - tot; + + // Since susceptible and infected are "fixed," + // we only need to change recovered + size_t nrecovered = prop[0u] * tot_left * n; + size_t ndeceased = prop[01] * tot_left * n; + + epiworld::AgentsSample sample_recover( + *model, + nrecovered, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_recover) + agent->change_state(model, 2, Queue::NoOne); + + epiworld::AgentsSample sample_deceased( + *model, + ndeceased, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_deceased) + agent->change_state(model, 3, Queue::NoOne); + + // Running the actions + model->actions_run(); + + return; + + }; + + return fun; + +} + + +/** + * @brief Creates an initial function for the SEIR-like models + * The function is used for the initial states of the model. +*/ +template +inline std::function*)> create_init_function_seir( + std::vector< double > proportions_ +) { + + // Checking widths + if (proportions_.size() != 2u) { + throw std::invalid_argument("-proportions_- must have two entries."); + } + + // proportions_ are values between 0 and 1, otherwise error + for (auto & v : proportions_) + if ((v < 0.0) || (v > 1.0)) + throw std::invalid_argument( + "-proportions_- must have values between 0 and 1." + ); + + + std::function*)> fun = + [proportions_] (epiworld::Model * model) -> void { + + // Figuring out information about the viruses + double tot = 0.0; + const auto & vpreval = model->get_prevalence_virus(); + const auto & vprop = model->get_prevalence_virus_as_proportion(); + double n = static_cast(model->size()); + for (size_t i = 0u; i < model->get_n_viruses(); ++i) + { + if (vprop[i]) + tot += vpreval[i]; + else + tot += vpreval[i] / n; + } + + // Putting the total into context + double tot_left = 1.0 - tot; + + // Since susceptible and infected are "fixed," + // we only need to change recovered + size_t nexposed = proportions_[0u] * tot * n; + size_t nrecovered = proportions_[1u] * tot_left * n; + + epiworld::AgentsSample sample_suscept( + *model, + nrecovered, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_suscept) + agent->change_state(model, 3, Queue::NoOne); + + epiworld::AgentsSample sample_exposed( + *model, + nexposed, + {1u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_exposed) + agent->change_state(model, 2, Queue::NoOne); + + // Running the actions + model->actions_run(); + + return; + + }; + + return fun; + +} + +/** + * @brief Creates an initial function for the SEIR-like models + * The function is used for the initial states of the model. +*/ +template +inline std::function*)> create_init_function_seird( + std::vector< double > proportions_ +) { + + // Checking widths + if (proportions_.size() != 3u) { + throw std::invalid_argument("-proportions_- must have three entries."); + } + + // proportions_ are values between 0 and 1, otherwise error + for (auto & v : proportions_) + if ((v < 0.0) || (v > 1.0)) + throw std::invalid_argument( + "-proportions_- must have values between 0 and 1." + ); + + // Last first two terms shouldn't add up to more than 1 + if ((proportions_[1u] + proportions_[2u]) > 1.0) + throw std::invalid_argument( + "The last two terms of -proportions_- must add up to less than 1." + ); + + std::function*)> fun = + [proportions_] (epiworld::Model * model) -> void { + + // Figuring out information about the viruses + double tot = 0.0; + const auto & vpreval = model->get_prevalence_virus(); + const auto & vprop = model->get_prevalence_virus_as_proportion(); + double n = static_cast(model->size()); + for (size_t i = 0u; i < model->get_n_viruses(); ++i) + { + if (vprop[i]) + tot += vpreval[i]; + else + tot += vpreval[i] / n; + } + + // Putting the total into context + double tot_left = 1.0 - tot; + + // Since susceptible and infected are "fixed," + // we only need to change recovered + size_t nexposed = proportions_[0u] * tot * n; + size_t nrecovered = proportions_[1u] * tot_left * n; + size_t ndeceased = proportions_[2u] * tot_left * n; + + epiworld::AgentsSample sample_suscept( + *model, + nrecovered, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_suscept) + agent->change_state(model, 3, Queue::NoOne); + + epiworld::AgentsSample sample_exposed( + *model, + nexposed, + {1u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_exposed) + agent->change_state(model, 2, Queue::NoOne); + + // Running the actions + model->actions_run(); + + // Setting the initial states for the deceased + epiworld::AgentsSample sample_deceased( + *model, + ndeceased, + {0u}, + true + ); + + // Setting up the initial states + for (auto & agent : sample_deceased) + agent->change_state(model, 4, Queue::NoOne); + + // Running the actions + model->actions_run(); + + return; + + }; + + return fun; + +} + + +#endif \ No newline at end of file diff --git a/inst/include/epiworld/models/models.hpp b/inst/include/epiworld/models/models.hpp index 2e27bfef..6c01abe4 100644 --- a/inst/include/epiworld/models/models.hpp +++ b/inst/include/epiworld/models/models.hpp @@ -2,7 +2,9 @@ #define EPIWORLD_MODELS_HPP namespace epimodels { - + + #include "init-functions.hpp" + #include "globalactions.hpp" #include "sis.hpp" #include "sir.hpp" diff --git a/inst/include/epiworld/models/seir.hpp b/inst/include/epiworld/models/seir.hpp index e310a169..301796ee 100644 --- a/inst/include/epiworld/models/seir.hpp +++ b/inst/include/epiworld/models/seir.hpp @@ -46,7 +46,7 @@ class ModelSEIR : public epiworld::Model ) -> void { // Getting the virus - auto v = p->get_virus(0); + auto v = p->get_virus(); // Does the agent become infected? if (m->runif() < 1.0/(v->get_incubation(m))) @@ -62,11 +62,22 @@ class ModelSEIR : public epiworld::Model ) -> void { // Does the agent recover? if (m->runif() < (m->par("Recovery rate"))) - p->rm_virus(0, m); + p->rm_virus(m); return; }; + /** + * @brief Set up the initial states of the model. + * @param proportions_ Double vector with the following values: + * - 0: Proportion of non-infected agents who are removed. + * - 1: Proportion of exposed agents to be set as infected. + */ + ModelSEIR & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + }; @@ -132,6 +143,18 @@ inline ModelSEIR::ModelSEIR( } +template +inline ModelSEIR & ModelSEIR::initial_states( + std::vector< double > proportions_, + std::vector< int > /**/ +) { + + Model::initial_states_fun = + create_init_function_seir(proportions_) + ; + + return *this; +} #endif diff --git a/inst/include/epiworld/models/seirconnected.hpp b/inst/include/epiworld/models/seirconnected.hpp index 4a8ed8e4..93f0eeea 100644 --- a/inst/include/epiworld/models/seirconnected.hpp +++ b/inst/include/epiworld/models/seirconnected.hpp @@ -35,7 +35,7 @@ class ModelSEIRCONN : public epiworld::Model epiworld_double recovery_rate ); - void run( + ModelSEIRCONN & run( epiworld_fast_uint ndays, int seed = -1 ); @@ -44,10 +44,20 @@ class ModelSEIRCONN : public epiworld::Model Model * clone_ptr(); + /** + * @brief Set the initial states of the model + * @param proportions_ Double vector with a single element: + * - The proportion of non-infected individuals who have recovered. + */ + ModelSEIRCONN & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + }; template -inline void ModelSEIRCONN::run( +inline ModelSEIRCONN & ModelSEIRCONN::run( epiworld_fast_uint ndays, int seed ) @@ -55,6 +65,8 @@ inline void ModelSEIRCONN::run( Model::run(ndays, seed); + return *this; + } template @@ -149,24 +161,21 @@ inline ModelSEIRCONN::ModelSEIRCONN( if (neighbor.get_state() == ModelSEIRCONN::INFECTED) { - for (const VirusPtr & v : neighbor.get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); + auto & v = neighbor.get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); } } @@ -181,7 +190,7 @@ inline ModelSEIRCONN::ModelSEIRCONN( if (which < 0) return; - p->add_virus( + p->set_virus( *m->array_virus_tmp[which], m, ModelSEIRCONN::EXPOSED @@ -201,7 +210,7 @@ inline ModelSEIRCONN::ModelSEIRCONN( { // Getting the virus - auto & v = p->get_virus(0u); + auto & v = p->get_virus(); // Does the agent become infected? if (m->runif() < 1.0/(v->get_incubation(m))) @@ -219,14 +228,11 @@ inline ModelSEIRCONN::ModelSEIRCONN( // Odd: Die, Even: Recover epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { + const auto & v = p->get_virus(); - // Recover - m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); - - } + // Recover + m->array_double_tmp[n_events++] = + 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); #ifdef EPI_DEBUG if (n_events == 0u) @@ -250,8 +256,7 @@ inline ModelSEIRCONN::ModelSEIRCONN( return; // Which roulette happen? - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + p->rm_virus(m); return ; @@ -327,4 +332,19 @@ inline ModelSEIRCONN::ModelSEIRCONN( } +template +inline ModelSEIRCONN & ModelSEIRCONN::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_seir(proportions_) + ; + + return *this; + +} + #endif diff --git a/inst/include/epiworld/models/seirconnected_logit.hpp b/inst/include/epiworld/models/seirconnected_logit.hpp index 6718f678..397aa95b 100644 --- a/inst/include/epiworld/models/seirconnected_logit.hpp +++ b/inst/include/epiworld/models/seirconnected_logit.hpp @@ -113,7 +113,7 @@ inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( for (auto & p: *_tracked_agents_infected) { - if (p->get_n_viruses() == 0) + if (p->get_virus() == nullptr) throw std::logic_error("Cannot be infected and have no viruses."); } @@ -154,7 +154,7 @@ inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( // Infecting the individual #ifdef EPI_DEBUG - if (_tracked_agents_infected->operator[](which)->get_n_viruses() == 0) + if (_tracked_agents_infected->operator[](which)->get_virus() == nullptr) { printf_epiworld("[epi-debug] date: %i\n", m->today()); @@ -167,8 +167,8 @@ inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( ); } #endif - p->add_virus( - _tracked_agents_infected->operator[](which)->get_virus(0u), + p->set_virus( + _tracked_agents_infected->operator[](which)->get_virus(), ModelSEIRCONNLogit::EXPOSED ); @@ -216,7 +216,7 @@ inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( { *_tracked_ninfected_next -= 1; - p->rm_virus(0, m); + p->rm_virus(m); return; } diff --git a/inst/include/epiworld/models/seird.hpp b/inst/include/epiworld/models/seird.hpp index a89e038b..fe608934 100644 --- a/inst/include/epiworld/models/seird.hpp +++ b/inst/include/epiworld/models/seird.hpp @@ -3,15 +3,7 @@ /** * @brief Template for a Susceptible-Exposed-Infected-Removed-Deceased (SEIRD) model - * - * @param model A Model object where to set up the SIR. - * @param vname std::string Name of the virus - * @param prevalence epiworld_double Initial prevalence the immune system - * @param transmission_rate epiworld_double Transmission rate of the virus - * @param avg_incubation_days epiworld_double Average incubation days of the virus - * @param recovery_rate epiworld_double Recovery rate of the virus. - * @param death_rate epiworld_double Death rate of the virus. - */ +*/ template class ModelSEIRD : public epiworld::Model { @@ -25,6 +17,18 @@ class ModelSEIRD : public epiworld::Model ModelSEIRD() {}; + /** + * @brief Constructor for the SEIRD model. + * + * @tparam TSeq Type of the sequence used in the model. + * @param model Reference to the SEIRD model. + * @param vname Name of the model. + * @param prevalence Prevalence of the disease. + * @param transmission_rate Transmission rate of the disease. + * @param avg_incubation_days Average incubation period of the disease. + * @param recovery_rate Recovery rate of the disease. + * @param death_rate Death rate of the disease. + */ ModelSEIRD( ModelSEIRD & model, std::string vname, @@ -35,6 +39,16 @@ class ModelSEIRD : public epiworld::Model epiworld_double death_rate ); + /** + * @brief Constructor for the SEIRD model. + * + * @param vname Name of the model. + * @param prevalence Initial prevalence of the disease. + * @param transmission_rate Transmission rate of the disease. + * @param avg_incubation_days Average incubation period of the disease. + * @param recovery_rate Recovery rate of the disease. + * @param death_rate Death rate of the disease. + */ ModelSEIRD( std::string vname, epiworld_double prevalence, @@ -50,7 +64,7 @@ class ModelSEIRD : public epiworld::Model ) -> void { // Getting the virus - auto v = p->get_virus(0); + auto v = p->get_virus(); // Does the agent become infected? if (m->runif() < 1.0/(v->get_incubation(m))) @@ -63,23 +77,20 @@ class ModelSEIRD : public epiworld::Model epiworld::UpdateFun update_infected = []( epiworld::Agent * p, epiworld::Model * m ) -> void { - - auto state = p->get_state(); - + // Odd: Die, Even: Recover epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { - - // Die - m->array_double_tmp[n_events++] = - v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); - - // Recover - m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); + + const auto & v = p->get_virus(); - } + // Die + m->array_double_tmp[n_events++] = + v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); + + // Recover + m->array_double_tmp[n_events++] = + 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); + #ifdef EPI_DEBUG if (n_events == 0u) @@ -106,13 +117,11 @@ class ModelSEIRD : public epiworld::Model if ((which % 2) == 0) // If odd { - size_t which_v = std::ceil(which / 2); - p->rm_agent_by_virus(which_v, m); + p->rm_agent_by_virus(m); } else { - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + p->rm_virus(m); } @@ -123,6 +132,12 @@ class ModelSEIRD : public epiworld::Model return; }; + + ModelSEIRD & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + }; @@ -195,6 +210,19 @@ inline ModelSEIRD::ModelSEIRD( } +template +inline ModelSEIRD & ModelSEIRD::initial_states( + std::vector< double > proportions_, + std::vector< int > /**/ +) { + + Model::initial_states_fun = + create_init_function_seird(proportions_) + ; + + return *this; + +} #endif \ No newline at end of file diff --git a/inst/include/epiworld/models/seirdconnected.hpp b/inst/include/epiworld/models/seirdconnected.hpp index bffdd0e3..148801b8 100644 --- a/inst/include/epiworld/models/seirdconnected.hpp +++ b/inst/include/epiworld/models/seirdconnected.hpp @@ -38,7 +38,7 @@ class ModelSEIRDCONN : public epiworld::Model epiworld_double death_rate ); - void run( + ModelSEIRDCONN & run( epiworld_fast_uint ndays, int seed = -1 ); @@ -47,10 +47,21 @@ class ModelSEIRDCONN : public epiworld::Model Model * clone_ptr(); + /** + * @brief Set up the initial states of the model. + * @param proportions_ Double vector with the following values: + * - 0: Proportion of non-infected agents who are removed. + * - 1: Proportion of exposed agents to be set as infected. + */ + ModelSEIRDCONN & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + }; template -inline void ModelSEIRDCONN::run( +inline ModelSEIRDCONN & ModelSEIRDCONN::run( epiworld_fast_uint ndays, int seed ) @@ -58,6 +69,8 @@ inline void ModelSEIRDCONN::run( Model::run(ndays, seed); + return *this; + } template @@ -154,25 +167,23 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( if (neighbor.get_state() == ModelSEIRDCONN::INFECTED) { - for (const VirusPtr & v : neighbor.get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); - - } + const auto & v = neighbor.get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + } } @@ -186,7 +197,7 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( if (which < 0) return; - p->add_virus( + p->set_virus( *m->array_virus_tmp[which], m, ModelSEIRDCONN::EXPOSED @@ -206,7 +217,7 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( { // Getting the virus - auto & v = p->get_virus(0u); + auto & v = p->get_virus(); // Does the agent become infected? if (m->runif() < 1.0/(v->get_incubation(m))) @@ -221,56 +232,50 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( } else if (state == ModelSEIRDCONN::INFECTED) { - - // Odd: Die, Even: Recover - epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { + // Odd: Die, Even: Recover + epiworld_fast_uint n_events = 0u; + const auto & v = p->get_virus(); // Die m->array_double_tmp[n_events++] = - v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); + v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); // Recover m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); - - } - - #ifdef EPI_DEBUG - if (n_events == 0u) - { + 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); + + #ifdef EPI_DEBUG + if (n_events == 0u) + { printf_epiworld( - "[epi-debug] agent %i has 0 possible events!!\n", - static_cast(p->get_id()) + "[epi-debug] agent %i has 0 possible events!!\n", + static_cast(p->get_id()) ); throw std::logic_error("Zero events in exposed."); - } - #else - if (n_events == 0u) + } + #else + if (n_events == 0u) return; - #endif - - - // Running the roulette - int which = roulette(n_events, m); - - if (which < 0) + #endif + + + // Running the roulette + int which = roulette(n_events, m); + + if (which < 0) return; - - // Which roulette happen? - if ((which % 2) == 0) // If odd - { - size_t which_v = std::ceil(which / 2); - p->rm_agent_by_virus(which_v, m); + // Which roulette happen? + if ((which % 2) == 0) // If odd + { + + p->rm_agent_by_virus(m); - } else { + } else { - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + p->rm_virus(m); - } + } return ; @@ -350,4 +355,18 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( } +template +inline ModelSEIRDCONN & ModelSEIRDCONN::initial_states( + std::vector< double > proportions_, + std::vector< int > /**/ +) { + + Model::initial_states_fun = + create_init_function_seird(proportions_) + ; + + return *this; + +} + #endif diff --git a/inst/include/epiworld/models/sir.hpp b/inst/include/epiworld/models/sir.hpp index 25bce99a..ae9ca2e3 100644 --- a/inst/include/epiworld/models/sir.hpp +++ b/inst/include/epiworld/models/sir.hpp @@ -31,6 +31,16 @@ class ModelSIR : public epiworld::Model epiworld_double transmission_rate, epiworld_double recovery_rate ); + + /** + * @brief Set the initial states of the model + * @param proportions_ Double vector with a single element: + * - The proportion of non-infected individuals who have recovered. + */ + ModelSIR & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); }; @@ -89,4 +99,18 @@ inline ModelSIR::ModelSIR( } +template +inline ModelSIR & ModelSIR::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) { + + Model::initial_states_fun = + create_init_function_sir(proportions_) + ; + + return *this; + +} + #endif diff --git a/inst/include/epiworld/models/sirconnected.hpp b/inst/include/epiworld/models/sirconnected.hpp index 84cf844d..dd7ec49e 100644 --- a/inst/include/epiworld/models/sirconnected.hpp +++ b/inst/include/epiworld/models/sirconnected.hpp @@ -11,12 +11,7 @@ class ModelSIRCONN : public epiworld::Model public: - ModelSIRCONN() { - - // tracked_agents_infected.reserve(1e4); - // tracked_agents_infected_next.reserve(1e4); - - }; + ModelSIRCONN() {}; ModelSIRCONN( ModelSIRCONN & model, @@ -37,17 +32,7 @@ class ModelSIRCONN : public epiworld::Model epiworld_double recovery_rate ); - // Tracking who is infected and who is not - // std::vector< epiworld::Agent* > tracked_agents_infected = {}; - // std::vector< epiworld::Agent* > tracked_agents_infected_next = {}; - // std::vector< epiworld_double > tracked_agents_weight = {}; - // std::vector< epiworld_double > tracked_agents_weight_next = {}; - - // int tracked_ninfected = 0; - // int tracked_ninfected_next = 0; - // epiworld_double tracked_current_infect_prob = 0.0; - - void run( + ModelSIRCONN & run( epiworld_fast_uint ndays, int seed = -1 ); @@ -56,24 +41,28 @@ class ModelSIRCONN : public epiworld::Model Model * clone_ptr(); + /** + * @brief Set the initial states of the model + * @param proportions_ Double vector with a single element: + * - The proportion of non-infected individuals who have recovered. + */ + ModelSIRCONN & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + }; template -inline void ModelSIRCONN::run( +inline ModelSIRCONN & ModelSIRCONN::run( epiworld_fast_uint ndays, int seed ) { - // tracked_agents_infected.clear(); - // tracked_agents_infected_next.clear(); - - // tracked_ninfected = 0; - // tracked_ninfected_next = 0; - // tracked_current_infect_prob = 0.0; - Model::run(ndays, seed); + return *this; } @@ -82,14 +71,6 @@ inline void ModelSIRCONN::reset() { Model::reset(); - - // Model::set_rand_binom( - // Model::size(), - // static_cast( - // Model::par("Contact rate"))/ - // static_cast(Model::size()) - // ); - return; } @@ -129,8 +110,6 @@ inline ModelSIRCONN::ModelSIRCONN( ) { - - epiworld::UpdateFun update_susceptible = []( epiworld::Agent * p, epiworld::Model * m ) -> void @@ -177,24 +156,25 @@ inline ModelSIRCONN::ModelSIRCONN( if (neighbor.get_state() == ModelSIRCONN::INFECTED) { - for (const VirusPtr & v : neighbor.get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); + if (neighbor.get_virus() == nullptr) + continue; + + auto & v = neighbor.get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + } } @@ -209,7 +189,7 @@ inline ModelSIRCONN::ModelSIRCONN( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; @@ -228,14 +208,10 @@ inline ModelSIRCONN::ModelSIRCONN( // Odd: Die, Even: Recover epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { - - // Recover - m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); - - } + // Recover + m->array_double_tmp[n_events++] = + 1.0 - (1.0 - p->get_virus()->get_prob_recovery(m)) * + (1.0 - p->get_recovery_enhancer(p->get_virus(), m)); #ifdef EPI_DEBUG if (n_events == 0u) @@ -259,13 +235,14 @@ inline ModelSIRCONN::ModelSIRCONN( return; // Which roulette happen? - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + p->rm_virus(m); return ; } else - throw std::logic_error("This function can only be applied to infected individuals. (SIR)") ; + throw std::logic_error( + "This function can only be applied to infected individuals. (SIR)" + ) ; return; @@ -325,5 +302,19 @@ inline ModelSIRCONN::ModelSIRCONN( } +template +inline ModelSIRCONN & ModelSIRCONN::initial_states( + std::vector< double > proportions_, + std::vector< int > /**/ +) { + + Model::initial_states_fun = + create_init_function_sir(proportions_) + ; + + return *this; + +} + #endif diff --git a/inst/include/epiworld/models/sird.hpp b/inst/include/epiworld/models/sird.hpp index 6109e94a..f9f51846 100644 --- a/inst/include/epiworld/models/sird.hpp +++ b/inst/include/epiworld/models/sird.hpp @@ -3,13 +3,6 @@ /** * @brief Template for a Susceptible-Infected-Removed-Deceased (SIRD) model - * - * @param model A Model object where to set up the SIRD. - * @param vname std::string Name of the virus - * @param initial_prevalence epiworld_double Initial prevalence - * @param initial_efficacy epiworld_double Initial susceptibility_reduction of the immune system - * @param initial_recovery epiworld_double Initial recovery_rate rate of the immune system - * @param initial_death epiworld_double Initial death_rate of the immune system */ template class ModelSIRD : public epiworld::Model @@ -18,6 +11,18 @@ class ModelSIRD : public epiworld::Model ModelSIRD() {}; + + /** + * @brief Constructs a new SIRD model with the given parameters. + * + * @param model The SIRD model to copy from. + * @param vname The name of the vertex associated with this model. + * @param prevalence The initial prevalence of the disease in the population. + * @param transmission_rate The rate at which the disease spreads from infected to susceptible individuals. + * @param recovery_rate The rate at which infected individuals recover and become immune. + * @param death_rate The rate at which infected individuals die. + */ + ///@{ ModelSIRD( ModelSIRD & model, std::string vname, @@ -34,6 +39,18 @@ class ModelSIRD : public epiworld::Model epiworld_double recovery_rate, epiworld_double death_rate ); + ///@} + + /** + * @brief Set the initial states of the model + * @param proportions_ Double vector with two elements: + * - The proportion of non-infected individuals who have recovered. + * - The proportion of non-infected individuals who have died. + */ + ModelSIRD & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); }; @@ -98,4 +115,18 @@ inline ModelSIRD::ModelSIRD( } +template +inline ModelSIRD & ModelSIRD::initial_states( + std::vector< double > proportions_, + std::vector< int > /**/ +) { + + Model::initial_states_fun = + create_init_function_sird(proportions_) + ; + + return *this; + +} + #endif diff --git a/inst/include/epiworld/models/sirdconnected.hpp b/inst/include/epiworld/models/sirdconnected.hpp index 4caa1cbf..68f2db0e 100644 --- a/inst/include/epiworld/models/sirdconnected.hpp +++ b/inst/include/epiworld/models/sirdconnected.hpp @@ -48,7 +48,7 @@ class ModelSIRDCONN : public epiworld::Model // int tracked_ninfected_next = 0; // epiworld_double tracked_current_infect_prob = 0.0; - void run( + ModelSIRDCONN & run( epiworld_fast_uint ndays, int seed = -1 ); @@ -61,21 +61,16 @@ class ModelSIRDCONN : public epiworld::Model }; template -inline void ModelSIRDCONN::run( +inline ModelSIRDCONN & ModelSIRDCONN::run( epiworld_fast_uint ndays, int seed ) { - // tracked_agents_infected.clear(); - // tracked_agents_infected_next.clear(); - - // tracked_ninfected = 0; - // tracked_ninfected_next = 0; - // tracked_current_infect_prob = 0.0; - Model::run(ndays, seed); + return *this; + } template @@ -180,24 +175,21 @@ inline ModelSIRDCONN::ModelSIRDCONN( if (neighbor.get_state() == ModelSIRDCONN::INFECTED) { - for (const VirusPtr & v : neighbor.get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; + const auto & v = neighbor.get_virus(); - m->array_virus_tmp[nviruses_tmp++] = &(*v); + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); } } @@ -212,7 +204,7 @@ inline ModelSIRDCONN::ModelSIRDCONN( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; @@ -229,55 +221,50 @@ inline ModelSIRDCONN::ModelSIRDCONN( { - // Odd: Die, Even: Recover - epiworld_fast_uint n_events = 0u; - for (const auto & v : p->get_viruses()) - { - + // Odd: Die, Even: Recover + epiworld_fast_uint n_events = 0u; + const auto & v = p->get_virus(); + // Die m->array_double_tmp[n_events++] = - v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); + v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); // Recover m->array_double_tmp[n_events++] = - 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); + 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); - } - -#ifdef EPI_DEBUG - if (n_events == 0u) - { - printf_epiworld( - "[epi-debug] agent %i has 0 possible events!!\n", - static_cast(p->get_id()) - ); - throw std::logic_error("Zero events in exposed."); - } -#else - if (n_events == 0u) - return; -#endif - - - // Running the roulette - int which = roulette(n_events, m); - - if (which < 0) - return; - - // Which roulette happen? - if ((which % 2) == 0) // If odd - { + #ifdef EPI_DEBUG + if (n_events == 0u) + { + printf_epiworld( + "[epi-debug] agent %i has 0 possible events!!\n", + static_cast(p->get_id()) + ); + throw std::logic_error("Zero events in exposed."); + } + #else + if (n_events == 0u) + return; + #endif - size_t which_v = std::ceil(which / 2); - p->rm_agent_by_virus(which_v, m); - } else { + // Running the roulette + int which = roulette(n_events, m); - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); + if (which < 0) + return; - } + // Which roulette happen? + if ((which % 2) == 0) // If odd + { + + p->rm_agent_by_virus(m); + + } else { + + p->rm_virus(m); + + } return ; diff --git a/inst/include/epiworld/models/sirlogit.hpp b/inst/include/epiworld/models/sirlogit.hpp index 00eff4e7..81ed8616 100644 --- a/inst/include/epiworld/models/sirlogit.hpp +++ b/inst/include/epiworld/models/sirlogit.hpp @@ -51,7 +51,7 @@ class ModelSIRLogit : public epiworld::Model epiworld_double prevalence ); - void run( + ModelSIRLogit & run( epiworld_fast_uint ndays, int seed = -1 ); @@ -70,13 +70,14 @@ class ModelSIRLogit : public epiworld::Model template -inline void ModelSIRLogit::run( +inline ModelSIRLogit & ModelSIRLogit::run( epiworld_fast_uint ndays, int seed ) { Model::run(ndays, seed); + return *this; } @@ -189,30 +190,30 @@ inline ModelSIRLogit::ModelSIRLogit( for (auto & neighbor: p->get_neighbors()) { - for (const VirusPtr & v : neighbor->get_viruses()) - { - - #ifdef EPI_DEBUG - if (nviruses_tmp >= m->array_virus_tmp.size()) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - baseline + - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) * - coef_exposure - ; - - // Applying the plogis function - m->array_double_tmp[nviruses_tmp] = 1.0/ - (1.0 + std::exp(-m->array_double_tmp[nviruses_tmp])); - - m->array_virus_tmp[nviruses_tmp++] = &(*v); - - } + if (neighbor->get_virus() == nullptr) + continue; + + auto & v = neighbor->get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= m->array_virus_tmp.size()) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + baseline + + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) * + coef_exposure + ; + + // Applying the plogis function + m->array_double_tmp[nviruses_tmp] = 1.0/ + (1.0 + std::exp(-m->array_double_tmp[nviruses_tmp])); + + m->array_virus_tmp[nviruses_tmp++] = &(*v); } @@ -226,7 +227,7 @@ inline ModelSIRLogit::ModelSIRLogit( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; @@ -242,7 +243,9 @@ inline ModelSIRLogit::ModelSIRLogit( // Computing recovery probability once double prob = 0.0; + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:prob) + #endif for (size_t i = 0u; i < _m->coefs_recover.size(); ++i) prob += p->operator[](i) * _m->coefs_recover[i]; @@ -250,7 +253,7 @@ inline ModelSIRLogit::ModelSIRLogit( prob = 1.0/(1.0 + std::exp(-prob)); if (prob > m->runif()) - p->rm_virus(0, m); + p->rm_virus(m); return; diff --git a/inst/include/epiworld/models/surveillance.hpp b/inst/include/epiworld/models/surveillance.hpp index 695f0a7f..81fd83b9 100644 --- a/inst/include/epiworld/models/surveillance.hpp +++ b/inst/include/epiworld/models/surveillance.hpp @@ -122,22 +122,20 @@ inline ModelSURV::ModelSURV( for (auto & neighbor: p->get_neighbors()) { - for (size_t i = 0u; i < neighbor->get_n_viruses(); ++i) - { + auto & v = neighbor->get_virus(); - auto & v = neighbor->get_virus(i); - - /* And it is a function of susceptibility_reduction as well */ - epiworld_double tmp_transmission = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor->get_transmission_reduction(v, m)) - ; - - m->array_double_tmp[nviruses_tmp] = tmp_transmission; - m->array_virus_tmp[nviruses_tmp++] = &(*v); + if (v == nullptr) + continue; - } + /* And it is a function of susceptibility_reduction as well */ + epiworld_double tmp_transmission = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor->get_transmission_reduction(v, m)) + ; + + m->array_double_tmp[nviruses_tmp] = tmp_transmission; + m->array_virus_tmp[nviruses_tmp++] = &(*v); } // No virus to compute on @@ -150,7 +148,7 @@ inline ModelSURV::ModelSURV( if (which < 0) return; - p->add_virus(*m->array_virus_tmp[which], m); + p->set_virus(*m->array_virus_tmp[which], m); return; }; @@ -160,7 +158,7 @@ inline ModelSURV::ModelSURV( [](epiworld::Agent * p, epiworld::Model * m) -> void { - epiworld::VirusPtr & v = p->get_virus(0u); + epiworld::VirusPtr & v = p->get_virus(); epiworld_double p_die = v->get_prob_death(m) * (1.0 - p->get_death_reduction(v, m)); epiworld_fast_uint days_since_exposed = m->today() - v->get_date(); @@ -184,7 +182,7 @@ inline ModelSURV::ModelSURV( // If past days infected + latent, then bye. if (days_since_exposed >= v->get_data()[1u]) { - p->rm_virus(0, m); + p->rm_virus(m); return; } diff --git a/inst/include/epiworld/tool-meat.hpp b/inst/include/epiworld/tool-meat.hpp index e70d2f60..a4925c55 100644 --- a/inst/include/epiworld/tool-meat.hpp +++ b/inst/include/epiworld/tool-meat.hpp @@ -66,7 +66,9 @@ inline ToolFun tool_fun_logit( size_t K = coefs_f.size(); epiworld_double res = 0.0; + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:res) + #endif for (size_t i = 0u; i < K; ++i) res += agent->operator[](vars.at(i)) * coefs_f.at(i); diff --git a/inst/include/epiworld/virus-bones.hpp b/inst/include/epiworld/virus-bones.hpp index 1a34f8e0..8397e453 100644 --- a/inst/include/epiworld/virus-bones.hpp +++ b/inst/include/epiworld/virus-bones.hpp @@ -30,8 +30,6 @@ class Virus { private: Agent * agent = nullptr; - int pos_in_agent = -99; ///< Location in the agent - int agent_exposure_number = -99; std::shared_ptr baseline_sequence = nullptr; std::shared_ptr virus_name = nullptr; @@ -46,7 +44,6 @@ class Virus { VirusFun incubation_fun = nullptr; // Setup parameters - std::vector< epiworld_double * > params = {}; std::vector< epiworld_double > data = {}; epiworld_fast_int state_init = -99; ///< Change of state when added to agent. @@ -67,7 +64,7 @@ class Virus { void set_sequence(TSeq sequence); Agent * get_agent(); - void set_agent(Agent * p, epiworld_fast_uint idx); + void set_agent(Agent * p); void set_date(int d); int get_date() const; diff --git a/inst/include/epiworld/virus-meat.hpp b/inst/include/epiworld/virus-meat.hpp index d5a10ae4..c0fae7a1 100644 --- a/inst/include/epiworld/virus-meat.hpp +++ b/inst/include/epiworld/virus-meat.hpp @@ -65,7 +65,9 @@ inline VirusFun virus_fun_logit( size_t K = coefs_f.size(); epiworld_double res = 0.0; + #if defined(__OPENMP) || defined(_OPENMP) #pragma omp simd reduction(+:res) + #endif for (size_t i = 0u; i < K; ++i) res += agent->operator[](vars.at(i)) * coefs_f.at(i); @@ -128,21 +130,9 @@ inline Agent * Virus::get_agent() } template -inline void Virus::set_agent(Agent * p, epiworld_fast_uint idx) +inline void Virus::set_agent(Agent * p) { - - #ifdef EPI_DEBUG - if (idx >= p->viruses.size()) - { - printf_epiworld( - "[epi-debug]Virus::set_agent id to set up is outside of range." - ); - } - #endif - - agent = p; - pos_in_agent = static_cast(idx); - + agent = p; } template diff --git a/man/ModelSEIR.Rd b/man/ModelSEIR.Rd index 3f43afab..103a8661 100644 --- a/man/ModelSEIR.Rd +++ b/man/ModelSEIR.Rd @@ -40,6 +40,12 @@ The \code{plot} function returns a plot of the SEIR model of class \description{ Susceptible Exposed Infected Recovered model (SEIR) } +\details{ +The \link{initial_states} function allows the user to set the initial state of the +model. The user must provide a vector of proportions indicating the following +values: (1) Proportion of non-infected agents who are removed, and (2) +Proportion of exposed agents to be set as infected. +} \examples{ model_seir <- ModelSEIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4) diff --git a/man/ModelSEIRD.Rd b/man/ModelSEIRD.Rd index 273f1292..eab5a77e 100644 --- a/man/ModelSEIRD.Rd +++ b/man/ModelSEIRD.Rd @@ -49,6 +49,13 @@ The \code{plot} function returns a plot of the SEIRD model of class \description{ Susceptible-Exposed-Infected-Recovered-Deceased model (SEIRD) } +\details{ +The \link{initial_states} function allows the user to set the initial state of the +model. The user must provide a vector of proportions indicating the following +values: (1) Proportion of exposed agents who are infected, (2) +proportion of non-infected agents already removed, and (3) proportion of +non-ifected agents already deceased. +} \examples{ model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, recovery_rate = 0.1, incubation_days = 4, diff --git a/man/ModelSEIRDCONN.Rd b/man/ModelSEIRDCONN.Rd index ad3f0edd..ca06e481 100644 --- a/man/ModelSEIRDCONN.Rd +++ b/man/ModelSEIRDCONN.Rd @@ -56,6 +56,13 @@ The \code{plot} function returns a plot of the SEIRDCONN model of class The SEIRD connected model implements a model where all agents are connected. This is equivalent to a compartmental model (\href{https://en.wikipedia.org/w/index.php?title=Compartmental_models_in_epidemiology&oldid=1155757336#The_SEIR_model}{wiki}). } +\details{ +The \link{initial_states} function allows the user to set the initial state of the +model. The user must provide a vector of proportions indicating the following +values: (1) Proportion of exposed agents who are infected, (2) +proportion of non-infected agents already removed, and (3) proportion of +non-ifected agents already deceased. +} \examples{ # An example with COVID-19 model_seirdconn <- ModelSEIRDCONN( diff --git a/man/ModelSIR.Rd b/man/ModelSIR.Rd index 4e285c38..4a75ad0d 100644 --- a/man/ModelSIR.Rd +++ b/man/ModelSIR.Rd @@ -39,6 +39,11 @@ infection.} \description{ SIR model } +\details{ +The \link{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. +} \examples{ model_sir <- ModelSIR(name = "COVID-19", prevalence = 0.01, transmission_rate = 0.9, recovery_rate = 0.1) diff --git a/man/ModelSIRCONN.Rd b/man/ModelSIRCONN.Rd index 39ecc915..eff21cdb 100644 --- a/man/ModelSIRCONN.Rd +++ b/man/ModelSIRCONN.Rd @@ -48,6 +48,11 @@ The \code{plot} function returns a plot of the SIRCONN model of class \description{ Susceptible Infected Removed model (SIR connected) } +\details{ +The \link{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. +} \examples{ model_sirconn <- ModelSIRCONN( name = "COVID-19", diff --git a/man/ModelSIRD.Rd b/man/ModelSIRD.Rd index 0fcc2757..a3fc448f 100644 --- a/man/ModelSIRD.Rd +++ b/man/ModelSIRD.Rd @@ -41,6 +41,12 @@ infection.} \description{ SIRD model } +\details{ +The \link{initial_states} function allows the user to set the initial state of the +model. The user must provide a vector of proportions indicating the following +values: (1) proportion of non-infected agents already removed, and (2) proportion of +non-ifected agents already deceased. +} \examples{ model_sird <- ModelSIRD( name = "COVID-19", diff --git a/man/ModelSIRDCONN.Rd b/man/ModelSIRDCONN.Rd index 81bae87f..59ac3e6d 100644 --- a/man/ModelSIRDCONN.Rd +++ b/man/ModelSIRDCONN.Rd @@ -51,6 +51,12 @@ The \code{plot} function returns a plot of the SIRDCONN model of class \description{ Susceptible Infected Removed Deceased model (SIRD connected) } +\details{ +The \link{initial_states} function allows the user to set the initial state of the +model. The user must provide a vector of proportions indicating the following +values: (1) proportion of non-infected agents already removed, and (2) proportion of +non-ifected agents already deceased. +} \examples{ model_sirdconn <- ModelSIRDCONN( name = "COVID-19", diff --git a/man/agents_smallworld.Rd b/man/agents_smallworld.Rd index 2bb35a10..afa19178 100644 --- a/man/agents_smallworld.Rd +++ b/man/agents_smallworld.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/agents.R, R/tool.R, R/virus.R +% Please edit documentation in R/agents.R, R/tool.R \name{agents_smallworld} \alias{agents_smallworld} \alias{agents} @@ -13,7 +13,6 @@ \alias{has_tool} \alias{change_state} \alias{get_agents_tools} -\alias{get_agents_viruses} \title{Load agents to a model} \usage{ agents_smallworld(model, n, k, d, p) @@ -35,8 +34,6 @@ has_tool(agent, tool) change_state(agent, model, state_new, queue = -99) get_agents_tools(model) - -get_agents_viruses(model) } \arguments{ \item{model}{Model object of class \link{epiworld_model}.} @@ -102,11 +99,6 @@ indicating whether the agent has the virus/tool or not. \item \code{get_agents_tools} returns a list of class \code{epiworld_agents_tools} with \code{epiworld_tools} (list of lists). } - -\itemize{ -\item \code{get_agents_viruses} returns a list of class \code{epiworld_agents_viruses} -with \code{epiworld_viruses} (list of lists). -} } \description{ These functions provide access to the network of the model. The network is diff --git a/man/epiworld-methods.Rd b/man/epiworld-methods.Rd index 0ea686a8..f0be4980 100644 --- a/man/epiworld-methods.Rd +++ b/man/epiworld-methods.Rd @@ -23,6 +23,7 @@ \alias{get_agents_data_ncols} \alias{get_virus} \alias{get_tool} +\alias{initial_states} \title{Methods for epiworldR objects} \usage{ queuing_on(x) @@ -64,6 +65,8 @@ get_agents_data_ncols(model) get_virus(model, virus_pos) get_tool(model, tool_pos) + +initial_states(model, proportions) } \arguments{ \item{x}{An object of class \code{epiworld_model}.} @@ -91,6 +94,9 @@ in the model} \item{tool_pos}{Integer. Relative location (starting from 0) of the tool in the model} + +\item{proportions}{Numeric vector. Proportions in which agents will be +distributed (see details).} } \value{ \itemize{ @@ -164,6 +170,10 @@ of \code{epiworld_model}. \itemize{ \item \code{get_tool} returns a \link{tool}. } + +\itemize{ +\item \code{inital_states} returns the model with an updated initial state. +} } \description{ The functions described in this section are methods for objects of class diff --git a/man/figures/README-transmission-net-1.png b/man/figures/README-transmission-net-1.png index e4bc3b0a..c9b20e2d 100644 Binary files a/man/figures/README-transmission-net-1.png and b/man/figures/README-transmission-net-1.png differ diff --git a/man/tool.Rd b/man/tool.Rd index 3a2e4172..cef1eecb 100644 --- a/man/tool.Rd +++ b/man/tool.Rd @@ -42,8 +42,6 @@ add_tool_n(model, tool, n) rm_tool(model, tool_pos) -rm_tool(model, tool_pos) - tool_fun_logit(vars, coefs, model) set_susceptibility_reduction(tool, prob) diff --git a/playground/benchmark-seir.R b/playground/benchmark-seir.R new file mode 100644 index 00000000..0c860d18 --- /dev/null +++ b/playground/benchmark-seir.R @@ -0,0 +1,50 @@ +# library(epiworldRdev) +# library(epiworldR) + +library(microbenchmark) + +ns <- c(1e3, 5e3, 1e4, 5e4, 1e5, 5e5) +ans <- vector("list", length(ns)) +names(ans) <- as.character(ns) +for (n in ns) { + + sir <- epiworldR::ModelSEIR( + name = "COVID-19", + prevalence = 0.01, + incubation_days = 7, + transmission_rate = 0.6, + recovery_rate = 0.5 + ) |> + epiworldR::agents_smallworld(n = n, k = 20, p = 0.0, d = FALSE) |> + epiworldR::add_virus( + epiworldR::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 + ) |> + epiworldR::verbose_off() + + + sirfast <- epiworldRdev::ModelSEIR( + name = "COVID-19", + prevalence = 0.01, + incubation_days = 7, + transmission_rate = 0.6, + recovery_rate = 0.5 + ) |> + epiworldRdev::agents_smallworld(n = n, k = 20, p = 0, d = FALSE) |> + epiworldRdev::add_virus( + epiworldRdev::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 + ) |> + epiworldRdev::verbose_off() + + + ans[[as.character(n)]] <- microbenchmark( + old = epiworldR::run(sir, ndays = 100, seed = 1912), + new = epiworldRdev::run(sirfast, ndays = 100, seed = 1912), + times = 10 + ) + + message("Simulation with ", n, " individuals finished.") + +} + +saveRDS(ans, "playground/benchmark-seir.rds") + diff --git a/playground/benchmark-seirconn.R b/playground/benchmark-seirconn.R new file mode 100644 index 00000000..6e6a947b --- /dev/null +++ b/playground/benchmark-seirconn.R @@ -0,0 +1,52 @@ +library(epiworldR) +library(epiworldRfaster) + +library(microbenchmark) + +ns <- c(1e3, 5e3, 1e4, 5e4, 1e5, 5e5) +ans <- vector("list", length(ns)) +names(ans) <- as.character(ns) +for (n in ns) { + + sir <- epiworldR::ModelSEIRCONN( + name = "COVID-19", + prevalence = 0.01, + n = n, + contact_rate = 4, + incubation_days = 7, + transmission_rate = 0.6, + recovery_rate = 0.5 + ) |> + epiworldR::add_virus( + epiworldR::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 + ) |> + epiworldR::verbose_off() + + + sirfast <- epiworldRfaster::ModelSEIRCONN( + name = "COVID-19", + prevalence = 0.01, + n = n, + contact_rate = 4, + incubation_days = 7, + transmission_rate = 0.6, + recovery_rate = 0.5 + ) |> + epiworldRfaster::add_virus( + epiworldRfaster::virus("COVID-19-beta", 0.01, 0.6, 0.5, 7), .2 + ) |> + epiworldRfaster::verbose_off() + + + ans[[as.character(n)]] <- microbenchmark( + old = epiworldR::run(sir, ndays = 100, seed = 1912), + new = epiworldRfaster::run(sirfast, ndays = 100, seed = 1912), + times = 10 + ) + + message("Simulation with ", n, " individuals finished.") + +} + +saveRDS(ans, "playground/benchmark-seirconn.rds") + diff --git a/src/agents.cpp b/src/agents.cpp index e7034c8d..04c6b8d1 100644 --- a/src/agents.cpp +++ b/src/agents.cpp @@ -92,7 +92,7 @@ SEXP add_virus_agent_cpp(SEXP agent, SEXP model, SEXP virus, int state_new, int cpp11::external_pointer> ptr_model(model); cpp11::external_pointer> ptr_virus(virus); - ptr_agent->add_virus(*ptr_virus, &(*ptr_model)); + ptr_agent->set_virus(*ptr_virus, &(*ptr_model)); return agent; diff --git a/src/cpp11.cpp b/src/cpp11.cpp index da431a64..2fe6ced3 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -453,6 +453,13 @@ extern "C" SEXP _epiworldR_get_network_cpp(SEXP model) { return cpp11::as_sexp(get_network_cpp(cpp11::as_cpp>(model))); END_CPP11 } +// model.cpp +SEXP initial_states_cpp(SEXP model, cpp11::doubles proportions); +extern "C" SEXP _epiworldR_initial_states_cpp(SEXP model, SEXP proportions) { + BEGIN_CPP11 + return cpp11::as_sexp(initial_states_cpp(cpp11::as_cpp>(model), cpp11::as_cpp>(proportions))); + END_CPP11 +} // tool.cpp SEXP tool_cpp(std::string name, double susceptibility_reduction, double transmission_reduction, double recovery_enhancer, double death_reduction); extern "C" SEXP _epiworldR_tool_cpp(SEXP name, SEXP susceptibility_reduction, SEXP transmission_reduction, SEXP recovery_enhancer, SEXP death_reduction) { @@ -754,20 +761,6 @@ extern "C" SEXP _epiworldR_set_name_virus_cpp(SEXP virus, SEXP name) { return cpp11::as_sexp(set_name_virus_cpp(cpp11::as_cpp>(virus), cpp11::as_cpp>(name))); END_CPP11 } -// virus.cpp -cpp11::writable::list get_agents_viruses_cpp(SEXP model); -extern "C" SEXP _epiworldR_get_agents_viruses_cpp(SEXP model) { - BEGIN_CPP11 - return cpp11::as_sexp(get_agents_viruses_cpp(cpp11::as_cpp>(model))); - END_CPP11 -} -// virus.cpp -SEXP print_agent_viruses_cpp(SEXP viruses); -extern "C" SEXP _epiworldR_print_agent_viruses_cpp(SEXP viruses) { - BEGIN_CPP11 - return cpp11::as_sexp(print_agent_viruses_cpp(cpp11::as_cpp>(viruses))); - END_CPP11 -} extern "C" { static const R_CallMethodDef CallEntries[] = { @@ -799,7 +792,6 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_get_agents_data_ncols_cpp", (DL_FUNC) &_epiworldR_get_agents_data_ncols_cpp, 1}, {"_epiworldR_get_agents_states_cpp", (DL_FUNC) &_epiworldR_get_agents_states_cpp, 1}, {"_epiworldR_get_agents_tools_cpp", (DL_FUNC) &_epiworldR_get_agents_tools_cpp, 1}, - {"_epiworldR_get_agents_viruses_cpp", (DL_FUNC) &_epiworldR_get_agents_viruses_cpp, 1}, {"_epiworldR_get_generation_time_cpp", (DL_FUNC) &_epiworldR_get_generation_time_cpp, 1}, {"_epiworldR_get_hist_tool_cpp", (DL_FUNC) &_epiworldR_get_hist_tool_cpp, 1}, {"_epiworldR_get_hist_total_cpp", (DL_FUNC) &_epiworldR_get_hist_total_cpp, 1}, @@ -828,10 +820,10 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_globalaction_tool_logit_cpp", (DL_FUNC) &_epiworldR_globalaction_tool_logit_cpp, 5}, {"_epiworldR_has_tool_cpp", (DL_FUNC) &_epiworldR_has_tool_cpp, 2}, {"_epiworldR_has_virus_cpp", (DL_FUNC) &_epiworldR_has_virus_cpp, 2}, + {"_epiworldR_initial_states_cpp", (DL_FUNC) &_epiworldR_initial_states_cpp, 2}, {"_epiworldR_make_saver_cpp", (DL_FUNC) &_epiworldR_make_saver_cpp, 10}, {"_epiworldR_print_agent_cpp", (DL_FUNC) &_epiworldR_print_agent_cpp, 3}, {"_epiworldR_print_agent_tools_cpp", (DL_FUNC) &_epiworldR_print_agent_tools_cpp, 1}, - {"_epiworldR_print_agent_viruses_cpp", (DL_FUNC) &_epiworldR_print_agent_viruses_cpp, 1}, {"_epiworldR_print_cpp", (DL_FUNC) &_epiworldR_print_cpp, 2}, {"_epiworldR_print_global_action_cpp", (DL_FUNC) &_epiworldR_print_global_action_cpp, 1}, {"_epiworldR_print_tool_cpp", (DL_FUNC) &_epiworldR_print_tool_cpp, 1}, diff --git a/src/model.cpp b/src/model.cpp index 79ce4ed3..80e80fe1 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -289,17 +289,29 @@ SEXP get_tool_model_cpp(SEXP model, int tool_pos) { [[cpp11::register]] cpp11::data_frame get_network_cpp(SEXP model) { - external_pointer> modelptr(model); - - std::vector from; - std::vector to; + external_pointer> modelptr(model); + + std::vector from; + std::vector to; - modelptr->write_edgelist(from, to); + modelptr->write_edgelist(from, to); - return cpp11::writable::data_frame({ - "from"_nm = from, - "to"_nm = to - }); + return cpp11::writable::data_frame({ + "from"_nm = from, + "to"_nm = to + }); } +[[cpp11::register]] +SEXP initial_states_cpp(SEXP model, cpp11::doubles proportions) { + + external_pointer> modelptr(model); + + std::vector< double > states_vec(proportions.begin(), proportions.end()); + + modelptr->initial_states(states_vec, std::vector< int >({})); + + return model; + +} diff --git a/src/virus.cpp b/src/virus.cpp index ebc0b857..d1fbb97b 100644 --- a/src/virus.cpp +++ b/src/virus.cpp @@ -276,31 +276,5 @@ SEXP set_name_virus_cpp(SEXP virus, std::string name) { return virus; } -// Function to get agent's viruses using get_agents_viruses() -[[cpp11::register]] -cpp11::writable::list get_agents_viruses_cpp(SEXP model) { - - cpp11::external_pointer> ptr(model); - - cpp11::writable::list viruses; - - for (auto & agent : ptr->get_agents()) - viruses.push_back( - cpp11::external_pointer< Viruses<> >( - new Viruses<>(agent.get_viruses()) - ) - ); - - return viruses; - -} - - -[[cpp11::register]] -SEXP print_agent_viruses_cpp(SEXP viruses) { - external_pointer> vptr(viruses); - vptr->print(); - return viruses; -} #undef WrapVirus \ No newline at end of file