From c1f9303ae2dc73b0d0fc4838971f19669ecaa1e1 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Thu, 7 Sep 2023 14:14:54 -0600 Subject: [PATCH] Adding new docker image for debugging. Fixed difference in rec rates --- .Rbuildignore | 1 + Makefile | 11 +- R/ModelSEIRD.R | 10 +- docker/Dockerfile | 19 +++ docker/Makefile | 2 + inst/include/epiworld/models/seir.hpp | 6 +- inst/include/epiworld/models/seird.hpp | 95 ++++++------- playground/transition_matrix_checker.Rmd | 171 ++++++++++++----------- 8 files changed, 175 insertions(+), 140 deletions(-) create mode 100644 docker/Dockerfile create mode 100644 docker/Makefile diff --git a/.Rbuildignore b/.Rbuildignore index 90f87b85..1a99162b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -33,3 +33,4 @@ README\.html \.gitattributes paper\..+ +docker \ No newline at end of file diff --git a/Makefile b/Makefile index 8303d860..a9796723 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,15 @@ +# Capture the current value of the version of the package in DESCRIPTION +VERSION := $(shell grep Version DESCRIPTION | sed -e 's/Version: //') + build: docs clean cd .. && rm -f epiworldR_*.tar.gz & R CMD build epiworldR debug: clean - EPI_CONFIG="-DEPI_DEBUG -Wall -pedantic -g" R CMD INSTALL . + docker run --rm -ti -w/mnt -v $(PWD):/mnt uofuepibio/epiworldr:debug make docker-debug + +docker-debug: + EPI_CONFIG="-DEPI_DEBUG -Wall -pedantic -g" R CMD INSTALL \ + --no-docs --build . install: build which R @@ -26,7 +33,7 @@ clean: docs: Rscript --vanilla -e 'devtools::document()' -.PHONY: build update check clean docs +.PHONY: build update check clean docs docker-debug checkv: build R CMD check --as-cran --use-valgrind epiworldR*.tar.gz diff --git a/R/ModelSEIRD.R b/R/ModelSEIRD.R index 87ee93d8..326193b5 100644 --- a/R/ModelSEIRD.R +++ b/R/ModelSEIRD.R @@ -45,8 +45,14 @@ ModelSEIRD <- function( ) { structure( - ModelSEIRD_cpp(name, prevalence, transmission_rate, incubation_days, - recovery_rate, death_rate), + ModelSEIRD_cpp( + name = name, + prevalence = prevalence, + transmission_rate = transmission_rate, + incubation_days = incubation_days, + recovery_rate = recovery_rate, + death_rate = death_rate + ), class = c("epiworld_seird", "epiworld_model") ) diff --git a/docker/Dockerfile b/docker/Dockerfile new file mode 100644 index 00000000..e82145ec --- /dev/null +++ b/docker/Dockerfile @@ -0,0 +1,19 @@ +FROM rocker/r-base:latest + +RUN \ + echo 'options(repos=c(CRAN="https://cloud.r-project.org"))' >> ~/.Rprofile && \ + Rscript --vanilla -e 'getOption("repos")' + + +RUN Rscript --vanilla -e "install.packages('cpp11', dependencies = TRUE, repos = 'https://cloud.r-project.org')" + + + +RUN mkdir ~/.R && \ + echo "CXXFLAGS=-g -O0" > ~/.R/Makevars && \ + echo "CXX14FLAGS=-g -O0" >> ~/.R/Makevars && \ + echo "CXX17FLAGS=-g -O0" >> ~/.R/Makevars && \ + echo "CXX11FLAGS=-g -O0" >> ~/.R/Makevars && \ + echo "SAFE_FLAGS=-g -O0" >> ~/.R/Makevars + +CMD ["bash"] \ No newline at end of file diff --git a/docker/Makefile b/docker/Makefile new file mode 100644 index 00000000..dffd5d55 --- /dev/null +++ b/docker/Makefile @@ -0,0 +1,2 @@ +all: + docker build -t uofuepibio/epiworldr:debug -f Dockerfile . \ No newline at end of file diff --git a/inst/include/epiworld/models/seir.hpp b/inst/include/epiworld/models/seir.hpp index 25f24085..e310a169 100644 --- a/inst/include/epiworld/models/seir.hpp +++ b/inst/include/epiworld/models/seir.hpp @@ -14,14 +14,13 @@ template class ModelSEIR : public epiworld::Model { -private: + +public: static const int SUSCEPTIBLE = 0; static const int EXPOSED = 1; static const int INFECTED = 2; static const int REMOVED = 3; -public: - ModelSEIR() {}; ModelSEIR( @@ -99,6 +98,7 @@ inline ModelSEIR::ModelSEIR( virus.set_prob_infecting(&model("Transmission rate")); virus.set_incubation(&model("Incubation days")); + virus.set_prob_recovery(&model("Recovery rate")); // Adding the tool and the virus model.add_virus(virus, prevalence); diff --git a/inst/include/epiworld/models/seird.hpp b/inst/include/epiworld/models/seird.hpp index adf121e2..a89e038b 100644 --- a/inst/include/epiworld/models/seird.hpp +++ b/inst/include/epiworld/models/seird.hpp @@ -15,15 +15,14 @@ template class ModelSEIRD : public epiworld::Model { -private: + +public: static const int SUSCEPTIBLE = 0; static const int EXPOSED = 1; static const int INFECTED = 2; static const int REMOVED = 3; static const int DECEASED = 4; -public: - ModelSEIRD() {}; ModelSEIRD( @@ -66,65 +65,60 @@ class ModelSEIRD : public epiworld::Model ) -> void { auto state = p->get_state(); - - if (state == ModelSEIRD::INFECTED) + + // 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)); - // 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)); - - } + // 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) - { - printf_epiworld( - "[epi-debug] agent %i has 0 possible events!!\n", - static_cast(p->get_id()) - ); - throw std::logic_error("Zero events in exposed."); - } + 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; + 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 + { + size_t which_v = std::ceil(which / 2); + p->rm_agent_by_virus(which_v, m); - // Running the roulette - int which = roulette(n_events, m); - - if (which < 0) - return; + } else { - // 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); - - } else { - - size_t which_v = std::floor(which / 2); - p->rm_virus(which_v, m); - - } + size_t which_v = std::floor(which / 2); + p->rm_virus(which_v, m); - return ; + } + + return ; - } else - throw std::logic_error("This function can only be applied to exposed or infected individuals. (SEIRD)") ; + return; @@ -165,6 +159,7 @@ inline ModelSEIRD::ModelSEIRD( virus.set_prob_infecting(&model("Transmission rate")); virus.set_incubation(&model("Incubation days")); virus.set_prob_death(&model("Death rate")); + virus.set_prob_recovery(&model("Recovery rate")); // Adding the tool and the virus model.add_virus(virus, prevalence); diff --git a/playground/transition_matrix_checker.Rmd b/playground/transition_matrix_checker.Rmd index 8720b49f..817920e0 100644 --- a/playground/transition_matrix_checker.Rmd +++ b/playground/transition_matrix_checker.Rmd @@ -6,23 +6,25 @@ output: html_document --- # SEIRD Example -```{r} -#' model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, -#' transmission_rate = 0.9, recovery_rate = 0.2, incubation_days = 4, -#' death_rate = 0.1) -#' -#' # Adding a small world population -#' agents_smallworld( -#' model_seird, -#' n = 500000, -#' k = 5, -#' d = FALSE, -#' p = .01 -#' ) -#' -#' # Running and printing -#' run(model_seird, ndays = 100, seed = 1912) -#' model_seird +```{r eval=FALSE} +library(epiworldR) + +model_seird <- ModelSEIRD(name = "COVID-19", prevalence = 0.01, +transmission_rate = 0.9, recovery_rate = 0.2, incubation_days = 4, +death_rate = 0.1) + +# Adding a small world population +agents_smallworld( + model_seird, + n = 500000, + k = 5, + d = FALSE, + p = .01 + ) + +# Running and printing +run(model_seird, ndays = 100, seed = 1912) +model_seird # Transition Probabilities: # - Susceptible 0.96 0.04 0.00 0.00 0.00 @@ -33,36 +35,38 @@ output: html_document pdie <- 0.1 precovery <- 0.2 -pnone <- (1-pdie)*(1-precovery) -patmost_one <- pdie*(1-precovery) + precovery*(1-pdie) +pnone <- (1 - pdie) * (1 - precovery) +patmost_one <- pdie * (1 - precovery) + precovery * (1 - pdie) p_die_given_infected <- pdie * (1-precovery) / (pnone + patmost_one) # 0.08163265 -p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.1836735 +p_recovery_given_infected <- precovery * (1 - pdie) / (pnone + patmost_one) # 0.1836735 +p_neither <- (1 - pdie) * (1 - precovery) / (pnone + patmost_one) # 0.7346939 ``` $P(Dying | Infected) = \frac{p_{die} * (1-p_{recovery})}{(1-p_{die})*(1-p_{recovery})+p_{die}*(1-p_{recovery}) + p_{recovery}*(1-p_{die})} = \frac{0.01 * (1-0.1)}{(1-0.01)*(1-0.1) + 0.01*(1-0.01)+0.1*(1-0.01)} = 0.009 \approx 0.01$ # SEIRDCONN Example -```{r} -#' # An example with COVID-19 -#' model_seirdconn <- ModelSEIRDCONN( -#' name = "COVID-19", -#' prevalence = 0.01, -#' n = 100000, -#' contact_rate = 2, -#' incubation_days = 7, -#' transmission_rate = 0.5, -#' recovery_rate = 0.3, -#' death_rate = 0.01 -#' ) -#' -#' # Running and printing -#' run(model_seirdconn, ndays = 100, seed = 1912) -#' model_seirdconn -#' -#' plot(model_seirdconn) - +```{r eval = FALSE} +# An example with COVID-19 +model_seirdconn <- ModelSEIRDCONN( + name = "COVID-19", + prevalence = 0.01, + n = 100000, + contact_rate = 2, + incubation_days = 7, + transmission_rate = 0.5, + recovery_rate = 0.3, + death_rate = 0.01 +) + +# Running and printing +run(model_seirdconn, ndays = 100, seed = 1912) +model_seirdconn + +plot(model_seirdconn) + +get_transition_probability(model_seirdconn) # Transition Probabilities: # - Susceptible 0.97 0.03 0.00 0.00 0.00 # - Exposed 0.00 0.86 0.14 0.00 0.00 @@ -111,28 +115,28 @@ p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.47 ``` # SIRD Example -```{r} -#' model_sird <- ModelSIRD( -#' name = "COVID-19", -#' prevalence = 0.01, -#' transmission_rate = 0.9, -#' recovery_rate = 0.2, -#' death_rate = 0.01 -#' ) -#' -#' # Adding a small world population -#' agents_smallworld( -#' model_sird, -#' n = 100000, -#' k = 5, -#' d = FALSE, -#' p = .01 -#' ) -#' -#' # Running and printing -#' run(model_sird, ndays = 100, seed = 1912) -#' model_sird - +```{r eval = FALSE} +model_sird <- ModelSIRD( + name = "COVID-19", + prevalence = 0.01, + transmission_rate = 0.9, + recovery_rate = 0.2, + death_rate = 0.01 +) + +# Adding a small world population +agents_smallworld( + model_sird, + n = 500000, + k = 5, + d = FALSE, + p = .01 + ) + +# Running and printing +run(model_sird, ndays = 100, seed = 1912) +model_sird +get_transition_probability(model_sird) # Transition Probabilities: # - Susceptible 0.96 0.04 0.00 0.00 # - Infected 0.00 0.80 0.20 0.01 @@ -149,28 +153,29 @@ p_recovery_given_infected <- precovery * (1-pdie) / (pnone + patmost_one) # 0.19 ``` # SISD Example -```{r} -#' model_sisd <- ModelSISD( -#' name = "COVID-19", -#' prevalence = 0.01, -#' transmission_rate = 0.9, -#' recovery_rate = 0.2, -#' death_rate = 0.01 -#' ) -#' -#' # Adding a small world population -#' agents_smallworld( -#' model_sisd, -#' n = 100000, -#' k = 5, -#' d = FALSE, -#' p = .01 -#' ) -#' -#' # Running and printing -#' run(model_sisd, ndays = 100, seed = 1912) -#' model_sisd - +```{r eval = FALSE} +model_sisd <- ModelSISD( + name = "COVID-19", + prevalence = 0.01, + transmission_rate = 0.9, + recovery_rate = 0.2, + death_rate = 0.01 + ) + +# Adding a small world population +agents_smallworld( + model_sisd, + n = 500000, + k = 5, + d = FALSE, + p = .01 + ) + +# Running and printing +run(model_sisd, ndays = 100, seed = 1912) +model_sisd + +get_transition_probability(model_sisd) # Transition Probabilities: # - Susceptible 0.57 0.43 0.00 # - Infected 0.20 0.79 0.01