diff --git a/.Rbuildignore b/.Rbuildignore index 73b43dc..6380b61 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,6 +7,7 @@ ^pkgdown$ ^README\.Rmd$ ^LICENSE\.md$ +^LICENSE\.note$ ^vignettes/web_only ^vignettes/denim_benchmark.Rmd ^vignettes/denim_vs_deSolve.Rmd @@ -20,4 +21,5 @@ ^extra_benchmark$ ^manuscript_plots$ ^\.vscode$ -^logo\.svg$ \ No newline at end of file +^logo\.svg$ +^CODE_OF_CONDUCT\.md$ diff --git a/.github/CODE_OF_CONDUCT.md b/.github/CODE_OF_CONDUCT.md new file mode 100644 index 0000000..50f91e9 --- /dev/null +++ b/.github/CODE_OF_CONDUCT.md @@ -0,0 +1,126 @@ +# Contributor Covenant Code of Conduct + +## Our Pledge + +We as members, contributors, and leaders pledge to make participation in our +community a harassment-free experience for everyone, regardless of age, body +size, visible or invisible disability, ethnicity, sex characteristics, gender +identity and expression, level of experience, education, socio-economic status, +nationality, personal appearance, race, caste, color, religion, or sexual +identity and orientation. + +We pledge to act and interact in ways that contribute to an open, welcoming, +diverse, inclusive, and healthy community. + +## Our Standards + +Examples of behavior that contributes to a positive environment for our +community include: + +* Demonstrating empathy and kindness toward other people +* Being respectful of differing opinions, viewpoints, and experiences +* Giving and gracefully accepting constructive feedback +* Accepting responsibility and apologizing to those affected by our mistakes, + and learning from the experience +* Focusing on what is best not just for us as individuals, but for the overall + community + +Examples of unacceptable behavior include: + +* The use of sexualized language or imagery, and sexual attention or advances of + any kind +* Trolling, insulting or derogatory comments, and personal or political attacks +* Public or private harassment +* Publishing others' private information, such as a physical or email address, + without their explicit permission +* Other conduct which could reasonably be considered inappropriate in a + professional setting + +## Enforcement Responsibilities + +Community leaders are responsible for clarifying and enforcing our standards of +acceptable behavior and will take appropriate and fair corrective action in +response to any behavior that they deem inappropriate, threatening, offensive, +or harmful. + +Community leaders have the right and responsibility to remove, edit, or reject +comments, commits, code, wiki edits, issues, and other contributions that are +not aligned to this Code of Conduct, and will communicate reasons for moderation +decisions when appropriate. + +## Scope + +This Code of Conduct applies within all community spaces, and also applies when +an individual is officially representing the community in public spaces. +Examples of representing our community include using an official e-mail address, +posting via an official social media account, or acting as an appointed +representative at an online or offline event. + +## Enforcement + +Instances of abusive, harassing, or otherwise unacceptable behavior may be +reported to the community leaders responsible for enforcement at anhptq@oucru.org. +All complaints will be reviewed and investigated promptly and fairly. + +All community leaders are obligated to respect the privacy and security of the +reporter of any incident. + +## Enforcement Guidelines + +Community leaders will follow these Community Impact Guidelines in determining +the consequences for any action they deem in violation of this Code of Conduct: + +### 1. Correction + +**Community Impact**: Use of inappropriate language or other behavior deemed +unprofessional or unwelcome in the community. + +**Consequence**: A private, written warning from community leaders, providing +clarity around the nature of the violation and an explanation of why the +behavior was inappropriate. A public apology may be requested. + +### 2. Warning + +**Community Impact**: A violation through a single incident or series of +actions. + +**Consequence**: A warning with consequences for continued behavior. No +interaction with the people involved, including unsolicited interaction with +those enforcing the Code of Conduct, for a specified period of time. This +includes avoiding interactions in community spaces as well as external channels +like social media. Violating these terms may lead to a temporary or permanent +ban. + +### 3. Temporary Ban + +**Community Impact**: A serious violation of community standards, including +sustained inappropriate behavior. + +**Consequence**: A temporary ban from any sort of interaction or public +communication with the community for a specified period of time. No public or +private interaction with the people involved, including unsolicited interaction +with those enforcing the Code of Conduct, is allowed during this period. +Violating these terms may lead to a permanent ban. + +### 4. Permanent Ban + +**Community Impact**: Demonstrating a pattern of violation of community +standards, including sustained inappropriate behavior, harassment of an +individual, or aggression toward or disparagement of classes of individuals. + +**Consequence**: A permanent ban from any sort of public interaction within the +community. + +## Attribution + +This Code of Conduct is adapted from the [Contributor Covenant][homepage], +version 2.1, available at +. + +Community Impact Guidelines were inspired by +[Mozilla's code of conduct enforcement ladder][https://github.com/mozilla/inclusion]. + +For answers to common questions about this code of conduct, see the FAQ at +. Translations are available at . + +[homepage]: https://www.contributor-covenant.org diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md new file mode 100644 index 0000000..99c342f --- /dev/null +++ b/.github/CONTRIBUTING.md @@ -0,0 +1,60 @@ +--- +editor_options: + markdown: + wrap: sentence +--- + +# Contributing to denim + +This outlines how to propose a change to the package denim. +For a detailed discussion on contributing to this and other tidyverse packages, please see the [development contributing guide](https://rstd.io/tidy-contrib) and our [code review principles](https://code-review.tidyverse.org/). + +## Fixing typos + +You can fix typos, spelling mistakes, or grammatical errors in the documentation directly using the GitHub web interface, as long as the changes are made in the *source* file, meaning: + +- YES: you edit a [roxygen2 comments](https://roxygen2.r-lib.org/articles/roxygen2.html) in a `.R` file below `R/`. + +- NO: you edit an `.Rd` file below `man/`. + +You can find the `.R` file that generates the `.Rd` by reading the comment in the first line. + +## Bigger changes + +If you want to make a bigger change, it's a good idea to first file an issue and make sure someone from the team agrees that it's needed. + +If you've found a bug, please file an issue that illustrates the bug with a reproducible example, preferrably using [reprex](https://www.tidyverse.org/help/#reprex) (this will also help you write a unit test, if needed). +See our guide on [how to create a great issue](https://code-review.tidyverse.org/issues/) for more advice. + +### Pull request process + +- Fork the package and clone onto your computer. + If you haven't done this before, we recommend using `usethis::create_from_github("thinhong/denim", fork = TRUE)`. + +- Install all development dependencies with `devtools::install_dev_deps()`, and then make sure the package passes R CMD check by running `devtools::check()`. + If R CMD check doesn't pass cleanly, it's a good idea to ask for help before continuing. + +- Create a Git branch for your pull request (PR). + We recommend using `usethis::pr_init("brief-description-of-change")`. + +- Make your changes, keep in mind that + + - We use [roxygen2](#0), for documentation. + + - We use [testthat](#0) for unit testings. + Contributions with test cases included are easier to accept. + + - New code should follow the tidyverse [style guide](https://style.tidyverse.org). + You can use the [styler](https://CRAN.R-project.org/package=styler) package to apply these styles, but please don't restyle code that has nothing to do with your PR. + +- Commit to git, and then create a PR by running `usethis::pr_push()`, and following the prompts in your browser. + The title of your PR should briefly describe the change. + The body of your PR should contain `Fixes #issue-number`. + +- For user-facing changes, add a bullet to the top of `NEWS.md` (i.e. just below the first header). + Follow the style described in . + +## Code of Conduct + +Please note that the `denim` project is released with a [Contributor Code of Conduct](CODE_OF_CONDUCT.md). +By contributing to this project you agree to abide by its terms. diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index f99dc87..519199b 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ Version: 1.2.2 -Date: 2025-08-18 05:41:32 UTC -SHA: 0d64dc2b4afd1b166100d94463729a191a8ce96c +Date: 2025-08-21 04:31:08 UTC +SHA: efec9ea304554e3cb0dfd99a39401ae47a5115dd diff --git a/DESCRIPTION b/DESCRIPTION index 91ccf3b..238d328 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,10 +10,14 @@ Authors@R: c( comment = c(ORCID = "0009-0000-2129-435X")), person("Marc", "Choisy", , "mchoisy@oucru.org", role = "aut", comment = c(ORCID = "0000-0002-5187-6390")), - person("Niels", "Lohman", ,role = "ctb"), - person("Bjoern", "Hoehrmann", ,email = "bjoern@hoehrmann.de", role = "ctb"), - person("Florian", "Loitsch", ,role = "ctb"), - person("Ingo", "Berg", ,role = "ctb") + person("Niels", "Lohman", ,role = c("ctb", "cph"), + comment="Author for the bundled JSON C++ library"), + person("Bjoern", "Hoehrmann", ,email = "bjoern@hoehrmann.de", role = c("ctb", "cph"), + comment="Contributor for the bundled JSON C++ library"), + person("Florian", "Loitsch", ,role = c("ctb", "cph"), + comment="Contributor for the bundled JSON C++ library"), + person("Ingo", "Berg", ,role = c("ctb", "cph"), + comment="Author for the bundled MuParser C++ library") ) Description: R package to build and simulate deterministic compartmental models that can be non-Markovian. Length of stay in each compartment can be defined to follow a parametric distribution (d_exponential(), d_gamma(), d_weibull(), d_lognormal()) or a non-parametric distribution (nonparametric()). Other supported types of transition from one compartment to another includes fixed transition (constant()), multinomial (multinomial()), fixed transition probability (transprob()). License: MIT + file LICENSE diff --git a/LICENSE.note b/LICENSE.note new file mode 100644 index 0000000..1ef509f --- /dev/null +++ b/LICENSE.note @@ -0,0 +1,5 @@ +The denim package as a whole is distributed under MIT. The denim package includes other open source software components. The following is a list of these components: + +* muParser: BSD +* JSON: MIT +* diffviewer: written by package authors \ No newline at end of file diff --git a/src/Compartment.cpp b/src/Compartment.cpp index f927666..cb4e512 100644 --- a/src/Compartment.cpp +++ b/src/Compartment.cpp @@ -23,10 +23,9 @@ void Compartment::setOutValues() { outTotals.resize(outTransitions.size(), 0); } -// TODO: set length for nested subCompartment void Compartment::setLengthSubCompartment() { size_t maxLength {0}; - // TODO: if competitive risk -> same subCompartment chain + // if competitive risk -> same subCompartment chain // ignore option to distribute initial values @@ -257,8 +256,6 @@ void Compartment::updateSubCompByDist(size_t iter, size_t outIndex) { for (size_t i = 0; i <= endIndex; ++i) { // if competing risks, compute out population using tmpSubComp instead double out = tmpSubComp[i] * (i >= transProbRef.size() ? 0 : transProbRef[i]); - // TODO: check formulation with 1 - exp(-rate) instead - // double out = tmpSubComp[i] * (i >= transProbRef.size() ? 0 : transProbRef[i]); outTotals[outIndex] += out; subCompartments[0][i] -= out; } @@ -268,8 +265,6 @@ void Compartment::updateSubCompByDist(size_t iter, size_t outIndex) { for (size_t i = 0; i <= endIndex; ++i) { double out = subCompartments[outIndex][i] * (i >= transProbRef.size() ? 0 : transProbRef[i]); - // TODO: check formulation with 1 - exp(-rate) instead - // double out = subCompartments[outIndex][i] * (1 - exp(-(i >= transProbRef.size() ? 0 : transProbRef[i]))); outTotals[outIndex] += out; subCompartments[outIndex][i] -= out; } diff --git a/src/test-denim.cpp b/src/test-denim.cpp index 78fc12f..111f216 100644 --- a/src/test-denim.cpp +++ b/src/test-denim.cpp @@ -44,7 +44,11 @@ context("Lognormal distribution") { test_polymorphism = std::make_shared(1, 0.25); test_that("getTransitionProb") { + // Check the resulting transition prob + expect_true(distr.getTransitionProb(0) == Approx(3.167e-05).margin(1e-06)); + expect_true(distr.getTransitionProb(1) == Approx(0.1098).margin(0.01)); expect_true(distr.getTransitionProb(2) == Approx(0.6106).margin(0.01)); + expect_true(distr.getTransitionProb(3) == Approx(0.8236).margin(0.01)); } test_that("maxDay") { @@ -60,7 +64,11 @@ context("Gamma distribution") { TransitionDiscreteGamma distr(1, 3); test_that("getTransitionProb") { - expect_true(distr.getTransitionProb(4) == Approx(0.4764).margin(0.01)); + expect_true(distr.getTransitionProb(0) == Approx(0.0803).margin(0.001)); + expect_true(distr.getTransitionProb(1) == Approx(0.2642).margin(0.001)); + expect_true(distr.getTransitionProb(2) == Approx(0.3746).margin(0.001)); + expect_true(distr.getTransitionProb(3) == Approx(0.4373).margin(0.001)); + expect_true(distr.getTransitionProb(4) == Approx(0.4764).margin(0.001)); } test_that("getRate") { @@ -76,7 +84,11 @@ context("Exponential distribution") { TransitionDiscreteExponential distr(0.5); test_that("getTransitionProb") { + // expect the same transprob for exponential distribution (i.e. constant hazard rate) + expect_true(distr.getTransitionProb(0) == Approx(0.3934).margin(0.01)); + expect_true(distr.getTransitionProb(1) == Approx(0.3934).margin(0.01)); expect_true(distr.getTransitionProb(2) == Approx(0.3934).margin(0.01)); + expect_true(distr.getTransitionProb(3) == Approx(0.3934).margin(0.01)); } test_that("getRate") { @@ -96,7 +108,10 @@ context("Weibull distribution") { TransitionDiscreteWeibull distr(3, 5); test_that("getTransitionProb") { - expect_true(distr.getTransitionProb(2) == Approx(0.5803).margin(0.01)); + expect_true(distr.getTransitionProb(0) == Approx(0.0041).margin(0.001)); + expect_true(distr.getTransitionProb(1) == Approx(0.1197).margin(0.001)); + expect_true(distr.getTransitionProb(2) == Approx(0.5803).margin(0.001)); + expect_true(distr.getTransitionProb(3) == Approx(0.9598).margin(0.001)); } test_that("getScale") { @@ -210,11 +225,11 @@ context("Testing distributions in model"){ expect_true(comp->getCompTotal()[1] == Approx(2).margin(0.001e-03)); break; case 4: - // check value for compartment V + // check value for compartment D expect_true(comp->getCompTotal()[1] == Approx(3.167e-06).margin(0.001e-06)); break; case 5: - // check value for compartment V + // check value for compartment R expect_true(comp->getCompTotal()[1] == Approx(0.04016).margin(0.001e-03)); break; } diff --git a/supplements/cached_time_scaling.rds b/supplements/cached_time_scaling.rds new file mode 100644 index 0000000..2dc7532 Binary files /dev/null and b/supplements/cached_time_scaling.rds differ diff --git a/tests/testthat/test-denim_dsl.R b/tests/testthat/test-denim_dsl.R index f8f9d3e..8a13a47 100644 --- a/tests/testthat/test-denim_dsl.R +++ b/tests/testthat/test-denim_dsl.R @@ -75,3 +75,158 @@ test_that("Compare result for model definition using DSL vs list", { abs_diff <- round(abs(mod_dsl[, c("S", "I", "R")] - mod_list[, c("S", "I", "R")]), 3) expect_true(all(abs_diff <= tolerance)) }) + + +test_that("Make sure the DSL generated the expected model definition", { + # ----- Set up ---- + initialValues <- c(S = 999, I = 1, R = 0, V = 0, IV = 0, D = 0) + + transitions <- denim_dsl({ + S -> I = beta * S * (I + IV) / N + S -> V = d_exponential(0.01) + 0.2 * I -> D = d_lognormal(2, 0.5) + 0.8 * I -> R = d_gamma(IR_rate, IR_shape) + V -> IV = 0.2 * beta * V * (I + IV) / N + IV -> R = nonparametric(iv_r_dist) + IV -> D = d_weibull(scale = 2, shape = 1.5) + }) + + parameters <- list( + beta = 0.9, + N = 1000, + IR_rate = 1/3, IR_shape = 1.8, + iv_r_dist = c(0, 0.15, 0.15, 0.05, 0.2, 0.2, 0.25) + ) + + simulationDuration <- 10 + timeStep <- 0.01 + + # ----- Try parsing model ---- + mod <- newModel(transitions=transitions, initialValues = initialValues, + parameters = parameters, + simulationDuration=simulationDuration, + errorTolerance = 0.001) + + # make sure the transitions are captured properly + expect_true(all(c( + "S -> I", "S -> V", "0.2 * I -> D", "0.8 * I -> R", + "V -> IV", "IV -> R", "IV -> R") %in% names(mod$transitions))) + + # ---- Make sure RHS is inferred properly ---- + # make sure all the RHS are properly converted to Transition object + expect_true( + all( + sapply( + names(mod$transitions), + \(trans){ + class(mod$transitions[[trans]])[1] == "Transition" + } + ) + ) + ) + + # ---- Check the inferred transition type ---- + expect_true(all( + sapply( + c("S -> I", "V -> IV"), + \(trans){ + mod$transitions[[trans]]$transition == "mathExpression" + }) + )) + expect_equal(mod$transitions$`0.2 * I -> D`$transition, "lognormal") + expect_equal(mod$transitions$`0.8 * I -> R`$transition, "gamma") + expect_equal(mod$transitions$`S -> V`$transition, "exponential") + expect_equal(mod$transitions$`IV -> R`$transition, "nonparametric") + expect_equal(mod$transitions$`IV -> D`$transition, "weibull") + + + # ----- Check whether distributional parameters are inferred ------ + expect_true(is.numeric(mod$transitions$`0.2 * I -> D`$mu) & + is.numeric(mod$transitions$`0.2 * I -> D`$sigma)) + expect_true(is.numeric(mod$transitions$`0.8 * I -> R`$rate) & + is.numeric(mod$transitions$`0.8 * I -> R`$shape)) + # also make sure that the waiting time in nonparametric distribution is captured properly + expect_true(is.numeric(mod$transitions$`IV -> R`$waitingTime)) + + # ----- Make sure distributional parameters are parsed correctly ----- + expect_equal(mod$transitions$`0.8 * I -> R`$rate, + parameters[["IR_rate"]]) + expect_equal(mod$transitions$`0.8 * I -> R`$shape, + parameters[["IR_shape"]]) + expect_setequal(mod$transitions$`IV -> R`$waitingTime, + parameters[["iv_r_dist"]]) +}) + +test_that("Make sure the model definition is parsed to JSON correctly", { + # ----- Set up ---- + transitions <- denim_dsl({ + S -> I = beta* S * (I + IV) / N + S -> V = 2 + 0.2 * I -> D = d_lognormal(mu = 2, sigma = 1/5) + 0.8 * I -> R = d_gamma(rate = IR_rate, shape = IR_shape) + V -> IV = 0.1 * beta * V * (I + IV) / N + IV -> R = d_exponential(2) + }) + initialValues <- c(S = 999, I = 1, R = 0, V = 0, IV = 0, D = 0) + parameters <- c( + beta = 1.2, + N = 1000, + IR_rate = 1/3, IR_shape = 1.8 + ) + simulationDuration <- 10 + timeStep <- 0.01 + + mod <- newModel(transitions=transitions, initialValues = initialValues, + parameters = parameters, + simulationDuration=simulationDuration, + errorTolerance = 0.001) + + # Make sure lower level JSON parsing functions work properly + modelToJson(mod) + strip_ws <- \(x) {gsub("\\s+", "", x)} + + # ----- check newJsonKeyPair function + ivKeyPairs <- newJsonKeyPair(names(initialValues), initialValues) + pmKeyPairs <- newJsonKeyPair(names(parameters), parameters) + json_simdur <- newJsonKeyPair("simulationDuration", simulationDuration) + + expect_setequal(strip_ws(ivKeyPairs), + strip_ws(c("\"S\":999","\"I\":1","\"R\":0","\"V\":0","\"IV\":0","\"D\":0"))) + expect_setequal(strip_ws(pmKeyPairs), + strip_ws(c("\"beta\":1.2","\"N\":1000","\"IR_rate\":0.333333333333333","\"IR_shape\":1.8"))) + expect_equal(json_simdur, '\"simulationDuration\": 10') + + # ---- Check JSON conversion for transition + # math expression + math_trans <- transitionToJson(mod$transitions$`0.2 * I -> D`) + # distribution based transition + dist_trans <- transitionToJson(mod$transitions$`IV -> R`) + expect_equal(math_trans, + '{\"transition\": \"lognormal\", \"mu\": 2, \"sigma\": 0.2, \"dist_init\": 0}') + expect_equal( + dist_trans, + '{\"transition\": \"exponential\", \"rate\": 2, \"dist_init\": 0}' + ) + + + # ----- check the function to make JSON object + ivObj <- newJsonObject(ivKeyPairs, inline = TRUE) + pmObj <- newJsonObject(pmKeyPairs, inline = TRUE) + expect_equal( + ivObj, + '{\"S\": 999, \"I\": 1, \"R\": 0, \"V\": 0, \"IV\": 0, \"D\": 0}') + expect_equal( + pmObj, + '{\"beta\": 1.2, \"N\": 1000, \"IR_rate\": 0.333333333333333, \"IR_shape\": 1.8}') + + # ----- check the function to make nested JSON object + expect_equal( + newJsonNestedObject("initialValues", ivObj), + '\"initialValues\": {\"S\": 999, \"I\": 1, \"R\": 0, \"V\": 0, \"IV\": 0, \"D\": 0}' + ) + expect_equal( + newJsonNestedObject("parameters", pmObj), + '\"parameters\": {\"beta\": 1.2, \"N\": 1000, \"IR_rate\": 0.333333333333333, \"IR_shape\": 1.8}' + ) + +}) diff --git a/vignettes/denim.Rmd b/vignettes/denim.Rmd index ac87f2f..15f6f35 100644 --- a/vignettes/denim.Rmd +++ b/vignettes/denim.Rmd @@ -294,7 +294,13 @@ head(simulation) ![](./images/denim_website-getting-started-extra2.png){width="500"} -```{r echo=FALSE} +```{=html} +
+ + Code for generating the distribution of I->R transition + +``` +```{r} # ------ Code to generate the distribution to demonstrate nonparametric --------- # Helper function compute_dist <- function(dist_func,..., timestep=0.05, error_tolerance=0.0001){ @@ -325,7 +331,7 @@ compute_dist <- function(dist_func,..., timestep=0.05, error_tolerance=0.0001){ } timestep <- 0.05 -# create a multimodal +# create a multimodal distribution first_dist <- compute_dist(pweibull, scale = 5, shape = 5, timestep = timestep) @@ -336,6 +342,9 @@ second_dist <- c(second_dist, rep(0, length(first_dist) - length(second_dist))) ir_dist <- first_dist + second_dist ``` +```{=html} +
+``` Assuming we have the distribution of time for the infected to recover ($P_{IR}$) as followed ```{r} diff --git a/vignettes/denim_benchmark.Rmd b/vignettes/denim_benchmark.Rmd index abccc83..f5ffc3c 100644 --- a/vignettes/denim_benchmark.Rmd +++ b/vignettes/denim_benchmark.Rmd @@ -1,8 +1,8 @@ --- -title: "denim benchmark" +title: "denim performance" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{denim benchmark} + %\VignetteIndexEntry{denim performance} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} python: @@ -10,41 +10,43 @@ python: bibliography: references.bib --- +# Benchmark + To assess denim's performance, we benchmark it against other available packages, listed in the table below. -+------------+--------------------------------------------------------------------+ -| Package | Note | -+============+====================================================================+ -| uSEIR | 2 implementations for uSEIR that will be tested | -| | | -| | - Implementation in Python | -| | | -| | - Implementation in Cython | -+------------+--------------------------------------------------------------------+ -| deSolve | 3 ways to define model for deSolve that will be tested | -| | | -| | - Model in R | -| | | -| | - Model in C | -| | | -| | - Model in Fortran | -+------------+--------------------------------------------------------------------+ -| IONISE | | -+------------+--------------------------------------------------------------------+ -| diffeqr | 2 ways to run model in diffeqr will be tested | -| | | -| | - Define model in R | -| | | -| | - Define model in Julia | -+------------+--------------------------------------------------------------------+ -| odin | A manual implementation of denim's algorithm, using odin is tested | -+------------+--------------------------------------------------------------------+ -| denim | Two methods for defining distributions are tested: | -| | | -| | - Built-in parametric functions | -| | | -| | - Histograms using the `nonparametric()` function | -+------------+--------------------------------------------------------------------+ ++--------------+--------------------------------------------------------------------+ +| Package | Note | ++==============+====================================================================+ +| uSEIR | 2 implementations for uSEIR that will be tested | +| | | +| | - Implementation in Python | +| | | +| | - Implementation in Cython | ++--------------+--------------------------------------------------------------------+ +| deSolve | 3 ways to define model for deSolve that will be tested | +| | | +| | - Model in R | +| | | +| | - Model in C | +| | | +| | - Model in Fortran | ++--------------+--------------------------------------------------------------------+ +| IONISE | | ++--------------+--------------------------------------------------------------------+ +| diffeqr | 2 ways to run model in diffeqr will be tested | +| | | +| | - Define model in R | +| | | +| | - Define model in Julia | ++--------------+--------------------------------------------------------------------+ +| odin | A manual implementation of denim's algorithm, using odin is tested | ++--------------+--------------------------------------------------------------------+ +| denim | Two methods for defining distributions are tested: | +| | | +| | - Built-in parametric functions | +| | | +| | - Histograms using the `nonparametric()` function | ++--------------+--------------------------------------------------------------------+ The benchmarking process was conducted on a Macbook Pro with M2 Pro chip, 16 GBs of RAM and 10 cores (6 performance and 4 efficiency). @@ -52,25 +54,25 @@ The benchmarking process was conducted on a Macbook Pro with M2 Pro chip, 16 GBs Version specifications -+-----------+-----------+-------------------------------------+ -| Name | Version | Note | -+===========+===========+=====================================+ -| R | 4.3.0 | | -+-----------+-----------+-------------------------------------+ -| Python | 3.9.20 | For `uSEIR` Python and Cython | -+-----------+-----------+-------------------------------------+ -| clang | 17.0.0 | C/C++ compiler for: | -| | | | -| | | - `uSEIR` Cython | -| | | | -| | | - `odin` | -| | | | -| | | - `denim` | -| | | | -| | | - `deSolve` (model in C, Fortran) | -+-----------+-----------+-------------------------------------+ -| Julia | 1.9.4 | For `diffeqr` package | -+-----------+-----------+-------------------------------------+ ++--------------+--------------+-------------------------------------+ +| Name | Version | Note | ++==============+==============+=====================================+ +| R | 4.3.0 | | ++--------------+--------------+-------------------------------------+ +| Python | 3.9.20 | For `uSEIR` Python and Cython | ++--------------+--------------+-------------------------------------+ +| clang | 17.0.0 | C/C++ compiler for: | +| | | | +| | | - `uSEIR` Cython | +| | | | +| | | - `odin` | +| | | | +| | | - `denim` | +| | | | +| | | - `deSolve` (model in C, Fortran) | ++--------------+--------------+-------------------------------------+ +| Julia | 1.9.4 | For `diffeqr` package | ++--------------+--------------+-------------------------------------+ R utility packages @@ -127,19 +129,19 @@ Source code: @@ -1237,28 +1239,53 @@ runtime_plot # ggsave("../manuscript_plots/benchmark_all_log_new.png", width = 8, height = 5, dpi = 300) ``` -## 9. Impact of timeStep in denim +It is worth noting that runtime for denim is also dependent on duration of time step (`timeStep` parameter for `sim`). This is discussed in details in the next section. + +------------------------------------------------------------------------ + +# Denim performance details + +## 1. Time - space complexity + +In terms of big-O notation, denim has a time complexity of $O(T S^2 n)$ and space complexity of $O(T S + S^2n)$. Where: + +- $T$ is the total simulated time steps -### 9.1. Run time scaling +- $S$ is the number of states/compartments -It is worth noting that runtime for denim is also dependent on duration of time step (`timeStep` parameter for `sim`). +- $n$ is the maximum length of the sub-compartment chain -The following plot demonstrates how run time changes as value for timeStep changes, using the same model for benchmarking. The values for timeStep being evaluated are `[0.01, 0.02, 0.05, 0.1, 0.25, 0.5, 1]`. +Under the assumption that individuals in a compartment can transition to every other compartment. + +However, under most realistic use cases, individuals can only transfer to a limited number of states in which case the number of transitions denoted $m$ is $m < S^2$. The time complexity would then be $O (T m n)$ and space complexity is $O(T S + mn)$. + +## 2. Impact of timeStep in denim + +### 2.1. Run time scaling + +A change in values for `timeStep` can directly change the value for $T$ and $n$ in the big-O notation, which directly impact the run time of the algorithm. + +The following plot demonstrates how run time changes as value for `timeStep` changes, using the same model for benchmarking. The values for timeStep being evaluated are `[0.01, 0.02, 0.05, 0.1, 0.25, 0.5, 1]`. ```{r echo=FALSE} timeStepDurations <- c(0.01, 0.02, 0.05, 0.1, 0.25, 0.5, 1) -runtime_scaling <- sapply(timeStepDurations, +benchmark_out <- lapply(timeStepDurations, \(x){ - time <- bench::mark( + bench_out <- bench::mark( sim(transitions = denim_model, initialValues = initialValues, parameters = parameters, simulationDuration = sim_duration, timeStep = x) ) - - time$median + bench_out } ) + +runtime_scaling <- sapply(benchmark_out, \(bm){ + bm$median +}) + +runtime_scaling ``` ```{r echo=FALSE, warning=FALSE} @@ -1274,11 +1301,11 @@ text(x = timeStep, y = max(unlist(runtime_scaling)), labels = "timeStep for benc text(x = 1, y = max(unlist(runtime_scaling)), labels = "default timeStep", col = "blueviolet", pos = 2) ``` -### 9.2. Impact on accuracy +### 2.2. Impact on accuracy Aside from run time, duration of time step also impacts the accuracy of `denim`'s output. -The following plots demonstrates how precision is compromised as we increase the duration for time step (with output from `deSolve` used as baseline). The values for timeStep being evaluated are `[0.01, 0.1, 0.5, 1]`. +The following plots demonstrates how precision is compromised as we increase the duration for time step (with output from `deSolve` used as baseline). The values for `timeStep` being evaluated are `[0.01, 0.1, 0.5, 1]`. ```{r echo = FALSE} timeStepDurations <- c(0.01, 0.1, 0.5, 1) @@ -1395,3 +1422,88 @@ reshaped_baseline %>% # glue::glue("../manuscript_plots/denim_dt_pairwise.png"), # width = 8, height = 5, dpi = 300) ``` + +## 3. Run time scaling with more transitions + +To assess the performance of `denim` as the model grows more complex, we measure the run-time as more transitions are added to the model. + +To do that, we consider the following sequential compartmental model where $S$ indicates the number of compartments + +$$ +X_1 \rightarrow X_2 \rightarrow X_3 .... \rightarrow X_S +$$ + +General settings: + +- `simulationDuration=2000` + +- `timeStep=0.01` (i.e., simulate the model for 200000 time steps) + +- The duration for any transition $X_i \rightarrow X_{i+1}$ is set to follow a log-normal distribution with $\mu = 2.2$ and $\sigma = 0.4$ + + - Maximal dwell time is \~13 time units (length of each sub-compartment chain is 1300) + +Under the given setup, the number of transitions $m$ is simply $S-1$ + +We will test the run time for model with `[3, 5, 25, 50, 100]` compartments (i.e., `[2, 4, 24, 49, 99]` transitions under this settings). + +```{r echo=FALSE} +cached_time_scaling <- if (file.exists("../supplements/cached_time_scaling.rds")) { + readRDS("../supplements/cached_time_scaling.rds") +} else { + NULL +} + +compartment_no <- c(3, 5, 25, 50, 100) +parameters <- c(mu = 2.2, sigma = 0.4) +timeStep <- 0.01 +simulationDuration <- 2000 + +runtime_scaling_out <- if(is.null(cached_time_scaling)){ + map_dfr(compartment_no, \(no){ + # build model with given number of compartments + mod <- lapply(1:(no-1), \(curr_comp){ + d_lognormal(mu = "mu", sigma = "sigma") + }) + names(mod) <- paste0("X", 1:(no-1) , "->X", 1:(no-1) + 1) + + # set the initial values + init_vals <- rep(0, no) + names(init_vals) <- paste0("X", 1:no) + init_vals[1] <- 1e6 + + # simulate the model + bench_out <- bench::mark( + sim( + transitions = mod, + initialValues = init_vals, + parameters = parameters, + simulationDuration = simulationDuration, + timeStep = timeStep, + errorTolerance = 0.01 + ) + ) + + bench_out %>% + mutate( + no_comp = no + ) + }) + +}else{ + cached_time_scaling +} +``` + +```{r echo=FALSE} +runtime_scaling_out %>% + ggplot(aes(x = no_comp, y = median)) + + geom_line(color = "cornflowerblue") + + geom_point(color = "cornflowerblue", shape = 1) + + labs( + title = "Run time scaling by compartments", + x = "Number of compartments", + y = "Run time (in seconds)" + ) + + theme_bw() +```