Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

14 higher contact rate error #15

Merged
merged 4 commits into from
Apr 7, 2024
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Updating version and adding test
gvegayon committed Apr 7, 2024
commit 22c7dcfea12c4cde74587738be0acd503c0dc0e7
9 changes: 5 additions & 4 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: epiworldR
Type: Package
Title: Fast Agent-Based Epi Models
Version: 0.0-4
Date: 2024-02-07
Version: 0.0-9999
Date: 2024-04-07
Authors@R: c(
person("Derek", "Meyer", role=c("aut","cre"),
email="[email protected]", comment = c(ORCID = "0009-0005-1350-6988")),
@@ -21,7 +21,7 @@ URL: https://github.com/UofUEpiBio/epiworldR,
https://uofuepibio.github.io/epiworldR-workshop/
BugReports: https://github.com/UofUEpiBio/epiworldR/issues
License: MIT + file LICENSE
RoxygenNote: 7.2.3
RoxygenNote: 7.3.1
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
LinkingTo: cpp11
@@ -30,7 +30,8 @@ Suggests:
rmarkdown,
tinytest,
netplot,
igraph
igraph,
data.table
Imports:
utils
VignetteBuilder: knitr
9 changes: 8 additions & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,11 @@
# epiworldR 0.0-3.9000 (development version)
# epiworldR 0.0-9999 (development version)

* Force model to update agents' states when running a simulation.
This was causing issues when calling `run_multiple()` after a single
call of `run()`. Reported on [14](https://github.com/UofUEpiBio/epiworldR/issues/14).


# epiworldR 0.0-4

* Added missing checks of tool class when adding a model with `add_too_n`.

16 changes: 9 additions & 7 deletions README.Rmd
Original file line number Diff line number Diff line change
@@ -13,12 +13,6 @@ knitr::opts_chunk$set(
)
```

## Versions
The virtual INSNA Sunbelt 2023 session can be found here: https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbelt2023-virtual

The in-person INSNA Sunbelt 2023 session can be found here:
https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbetl2023-inperson

# epiworldR

<!-- badges: start -->
@@ -205,7 +199,7 @@ net <- get_transmissions(sir)

# Plotting
library(epiworldR)
library(epiworldR)
library(netplot)
x <- igraph::graph_from_edgelist(
as.matrix(net[,2:3]) + 1
)
@@ -242,6 +236,14 @@ head(ans$reproductive)
plot(ans$reproductive)
```

# Tutorials

- The virtual INSNA Sunbelt 2023 session can be found here: https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbelt2023-virtual

- The in-person INSNA Sunbelt 2023 session can be found here:
https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbetl2023-inperson


# Existing Alternatives

Several alternatives to `epiworldR` exist and provide researchers with a range of options, each with its own unique features and strengths, enabling the exploration and analysis of infectious disease dynamics through agent-based modeling. Below is a manually curated table of existing alternatives including ABM [@ABM], abmR [@abmR], cystiSim [@cystiSim], villager [@villager], and RNetLogo [@RNetLogo].
44 changes: 22 additions & 22 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,14 +1,6 @@

<!-- README.md is generated from README.Rmd. Please edit that file -->

## Versions

The virtual INSNA Sunbelt 2023 session can be found here:
<https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbelt2023-virtual>

The in-person INSNA Sunbelt 2023 session can be found here:
<https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbetl2023-inperson>

# epiworldR

<!-- badges: start -->
@@ -127,11 +119,11 @@ summary(sir)
#> Number of entities : 0
#> Days (duration) : 50 (of 50)
#> Number of viruses : 1
#> Last run elapsed t : 456.00ms
#> Last run speed : 10.94 million agents x day / second
#> Last run elapsed t : 155.00ms
#> Last run speed : 32.12 million agents x day / second
#> Rewiring : off
#>
#> Global actions:
#> Global events:
#> (none)
#>
#> Virus(es):
@@ -303,16 +295,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), 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
#> Agent: 0, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 1, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 2, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 3, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 4, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 5, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 6, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 7, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 8, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
#> Agent: 9, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
#> ... 99990 more agents ...
```

@@ -346,7 +338,7 @@ net <- get_transmissions(sir)

# Plotting
library(epiworldR)
library(epiworldR)
library(netplot)
#> Loading required package: grid
x <- igraph::graph_from_edgelist(
as.matrix(net[,2:3]) + 1
@@ -412,6 +404,14 @@ plot(ans$reproductive)

<img src="man/figures/README-unnamed-chunk-6-1.png" width="100%" />

# Tutorials

- The virtual INSNA Sunbelt 2023 session can be found here:
<https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbelt2023-virtual>

- The in-person INSNA Sunbelt 2023 session can be found here:
<https://github.com/UofUEpiBio/epiworldR-workshop/tree/sunbetl2023-inperson>

# Existing Alternatives

Several alternatives to `epiworldR` exist and provide researchers with a
7 changes: 7 additions & 0 deletions inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp
Original file line number Diff line number Diff line change
@@ -83,6 +83,13 @@ inline LFMCMCProposalFun<TData> make_proposal_norm_reflective(

}

#ifdef EPI_DEBUG
for (auto & p : params_now)
if (p < lb || p > ub)
throw std::range_error("The parameter is out of bounds.");
#endif


return;

};
24 changes: 11 additions & 13 deletions inst/include/epiworld/model-meat.hpp
Original file line number Diff line number Diff line change
@@ -2039,21 +2039,19 @@ inline void Model<TSeq>::reset() {
}
#endif

}
else
{
for (auto & p : population)
p.reset();
}

#ifdef EPI_DEBUG
for (auto & a: population)
{
if (a.get_state() != 0u)
throw std::logic_error("Model::reset population doesn't match."
"Some agents are not in the baseline state.");
}
#endif
for (auto & p : population)
p.reset();

#ifdef EPI_DEBUG
for (auto & a: population)
{
if (a.get_state() != 0u)
throw std::logic_error("Model::reset population doesn't match."
"Some agents are not in the baseline state.");
}
#endif

if (entities_backup.size() != 0)
{
28 changes: 28 additions & 0 deletions inst/include/epiworld/models/sirlogit.hpp
Original file line number Diff line number Diff line change
@@ -3,6 +3,34 @@
#ifndef EPIWORLD_MODELS_SIRLOGIT_HPP
#define EPIWORLD_MODELS_SIRLOGIT_HPP


/**
* @brief Template for a Susceptible-Infected-Removed (SIR) model
*
* @details
* In this model, infection and recoveru probabilities are computed
* using a logit model. Particularly, the probability of infection
* is computed as:
*
* \f[
* \frac{1}{1 + \exp\left(-\left(\beta_0 E_i + \sum_{i=1}^{n} \beta_i x_i\right)\right)}
* \f]
*
* where \f$\beta_0\f$ is the exposure coefficient and \f$E_i\f$ is the exposure
* number, \f$\beta_i\f$ are the
* coefficients for the features \f$x_i\f$ of the agents, and \f$n\f$ is the
* number of features. The probability of recovery is computed as:
*
* \f[
* \frac{1}{1 + \exp\left(-\left(\sum_{i=1}^{n} \beta_i x_i\right)\right)}
* \f]
*
* where \f$\beta_i\f$ are the coefficients for the features \f$x_i\f$ of the agents,
* and \f$n\f$ is the number of features.
*
* @param TSeq Type of the sequence (e.g. std::vector, std::deque)

*/
template<typename TSeq = EPI_DEFAULT_TSEQ>
class ModelSIRLogit : public epiworld::Model<TSeq>
{
111 changes: 31 additions & 80 deletions inst/tinytest/test-multiple.R
Original file line number Diff line number Diff line change
@@ -1,74 +1,20 @@
# result <- suppressWarnings(suppressMessages(invisible({

library(epiworldR)

# library(data.table)
# library(EpiEstim)
# library(dplyr)
# library(ggplot2)
# library(tidyverse)
# library(EpiNow2)
# })))


model_seircon=ModelSEIRCONN(name="covid",n=50000,prevalence = 0.001,contact_rate = 20,transmission_rate = 0.5,recovery_rate = 1/7,incubation_days = 3)





# run(model_seircon,ndays=50,seed=1912)
#> _________________________________________________________________________
#> |Running the model...
#> |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> | done.
# summary(model_seircon)
#> ________________________________________________________________________________
#> ________________________________________________________________________________
#> SIMULATION STUDY
#>
#> Name of the model : Susceptible-Exposed-Infected-Removed (SEIR) (connected)
#> Population size : 50000
#> Agents' data : (none)
#> Number of entities : 0
#> Days (duration) : 50 (of 50)
#> Number of viruses : 1
#> Last run elapsed t : 473.00ms
#> Last run speed : 5.28 million agents x day / second
#> Rewiring : off
#>
#> Global events:
#> (none)
#>
#> Virus(es):
#> - covid (baseline prevalence: 0.10%)
#>
#> Tool(s):
#> (none)
#>
#> Model parameters:
#> - Avg. Incubation days : 3.0000
#> - Contact rate : 20.0000
#> - Prob. Recovery : 0.1429
#> - Prob. Transmission : 0.5000
#>
#> Distribution of the population at time 50:
#> - (0) Susceptible : 49950 -> 0
#> - (1) Exposed : 50 -> 0
#> - (2) Infected : 0 -> 140
#> - (3) Recovered : 0 -> 49860
#>
#> Transition Probabilities:
#> - Susceptible 0.60 0.40 0.00 0.00
#> - Exposed 0.00 0.66 0.34 0.00
#> - Infected 0.00 0.00 0.86 0.14
#> - Recovered 0.00 0.00 0.00 1.00



# plot(model_seircon)


n <- 5e3
days <- 50
nsims <- 50

model_seircon=ModelSEIRCONN(
name="covid",
n = n,
prevalence = 0.001,
contact_rate = 20,
transmission_rate = 0.5,
recovery_rate = 1/7,
incubation_days = 3
)

run(model_seircon,ndays=days,seed=1912)

saver <- make_saver(
"total_hist",
@@ -77,17 +23,22 @@ saver <- make_saver(
"reproductive",
"generation"
)
run_multiple(model_seircon,ndays=50,nsim=100,seed=1972,saver=saver)
#> Starting multiple runs (100) using 1 thread(s)
#> _________________________________________________________________________
#> _________________________________________________________________________
#> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
#> done.

run_multiple(
model_seircon,
ndays=days,
nsim=nsims,
seed=1972,
saver=saver,
nthreads = 2L
)

res1 <- run_multiple_get_results(model_seircon)
#> Warning in run_multiple_get_results(model_seircon): When retrieving the saved
#> results, for the case of transmission, there were no observations.
#> Error in attributes(.Data) <- c(attributes(.Data), attrib): all attributes must have names [4 does not]
res1$reproductive
#> Error in eval(expr, envir, enclos): object 'res1' not found

res1 <- lapply(res1, data.table::as.data.table)

avg_infected <- res1$total_hist[date==50 & state=="Infected"]$counts |> mean()
avg_recovered <- res1$total_hist[date==50 & state=="Recovered"]$counts |> mean()

expect_true(abs(avg_infected - 13.8) < 50)
expect_true(abs(avg_recovered - 4986) < 50)
10 changes: 5 additions & 5 deletions man/ModelDiffNet.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading
Loading