From bdb611241407d598e929181174ca3531ec05fdf5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 9 Sep 2025 19:37:40 +0000 Subject: [PATCH 1/8] Initial plan From d3f00ed36d3ec75799e19a4fded03b3e75e499b7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 9 Sep 2025 19:49:11 +0000 Subject: [PATCH 2/8] Complete MkDocs migration with working build Co-authored-by: damonbayer <8095080+damonbayer@users.noreply.github.com> --- .github/workflows/website.yaml | 4 +- docs/Makefile | 35 +- docs/developer_documentation.md | 109 ++++ docs/index.md | 125 ++++ docs/javascripts/mathjax.js | 16 + docs/reference/additional.md | 25 + docs/reference/datasets.md | 13 + docs/reference/deterministic.md | 3 + docs/reference/distributions.md | 3 + docs/reference/index.md | 15 + docs/reference/latent.md | 21 + docs/reference/model.md | 9 + docs/reference/observation.md | 9 + docs/reference/process.md | 3 + docs/reference/randomvariable.md | 3 + docs/stylesheets/extra.css | 35 ++ docs/tutorials/basic_renewal_model.qmd | 407 +++++++++++++ docs/tutorials/day_of_the_week.qmd | 426 +++++++++++++ docs/tutorials/extending_pyrenew.qmd | 297 ++++++++++ docs/tutorials/getting_started.qmd | 60 ++ docs/tutorials/hospital_admissions_model.qmd | 591 +++++++++++++++++++ docs/tutorials/index.md | 12 + docs/tutorials/periodic_effects.qmd | 105 ++++ hook_scripts/pre-commit-md-placeholder.sh | 6 +- mkdocs.yml | 110 ++++ pyproject.toml | 11 +- uv.lock | 363 ++++++++---- 27 files changed, 2659 insertions(+), 157 deletions(-) create mode 100644 docs/developer_documentation.md create mode 100644 docs/index.md create mode 100644 docs/javascripts/mathjax.js create mode 100644 docs/reference/additional.md create mode 100644 docs/reference/datasets.md create mode 100644 docs/reference/deterministic.md create mode 100644 docs/reference/distributions.md create mode 100644 docs/reference/index.md create mode 100644 docs/reference/latent.md create mode 100644 docs/reference/model.md create mode 100644 docs/reference/observation.md create mode 100644 docs/reference/process.md create mode 100644 docs/reference/randomvariable.md create mode 100644 docs/stylesheets/extra.css create mode 100644 docs/tutorials/basic_renewal_model.qmd create mode 100644 docs/tutorials/day_of_the_week.qmd create mode 100644 docs/tutorials/extending_pyrenew.qmd create mode 100644 docs/tutorials/getting_started.qmd create mode 100644 docs/tutorials/hospital_admissions_model.qmd create mode 100644 docs/tutorials/index.md create mode 100644 docs/tutorials/periodic_effects.qmd create mode 100644 mkdocs.yml diff --git a/.github/workflows/website.yaml b/.github/workflows/website.yaml index 2ed4f2d5..186fbd29 100644 --- a/.github/workflows/website.yaml +++ b/.github/workflows/website.yaml @@ -41,7 +41,7 @@ jobs: - name: Render documentation index from README run: | - cp README.md docs/source/index.md + cp README.md docs/index.md # Building - name: Build HTML using uv @@ -56,7 +56,7 @@ jobs: with: name: github-pages retention-days: 7 - path: docs/build/html/ + path: docs/site/ deploy: # Deploy to the github-pages environment diff --git a/docs/Makefile b/docs/Makefile index 3c246103..5d67e470 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -1,20 +1,26 @@ -# Minimal makefile for Sphinx documentation +# Makefile for MkDocs documentation # # You can set these variables from the command line, and also # from the environment for the first two. -SPHINXOPTS ?= -SPHINXBUILD ?= sphinx-build -SOURCEDIR = source -BUILDDIR = build +MKDOCSBUILD ?= mkdocs +SOURCEDIR = . +BUILDDIR = site # Put it first so that "make" without argument is like "make help". help: - @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + @echo "Please use 'make ' where is one of" + @echo " html to make standalone HTML files" + @echo " serve to serve the documentation locally" + @echo " tutorials to build all tutorial formats" + @echo " tutorials_md to build tutorial markdown files" + @echo " tutorials_ipynb to build tutorial notebook files" + @echo " tutorials_py to build tutorial Python files" + @echo " tutorials_clean to clean tutorial build artifacts" -.PHONY: help Makefile +.PHONY: help html serve Makefile -TUTORIALSDIR = $(SOURCEDIR)/tutorials +TUTORIALSDIR = tutorials .PHONY: tutorials tutorials_md tutorials_ipynb tutorials_py tutorials_clean echo_tutorial_source @@ -23,11 +29,11 @@ $(TUTORIALSDIR)/%.md: $(TUTORIALSDIR)/%.qmd quarto render $< --to gfm # Pattern rule to build a notebook file from a qmd file -$(TUTORIALS_DIR)/%.ipynb: $(TUTORIALS_DIR)/%.qmd +$(TUTORIALSDIR)/%.ipynb: $(TUTORIALSDIR)/%.qmd quarto convert $< --output $@ # Pattern rule to build a Python file from a notebook file -$(TUTORIALS_DIR)/%.py: $(TUTORIALS_DIR)/%.ipynb +$(TUTORIALSDIR)/%.py: $(TUTORIALSDIR)/%.ipynb sed -i.bak 's/"cell_type": "raw"/"cell_type": "markdown"/g' $< jupyter nbconvert --to python $< rm $<.bak @@ -62,9 +68,8 @@ tutorials_clean: .PHONY: html html: $(TUTORIALS_MD) - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + @$(MKDOCSBUILD) build -# Catch-all target: route all unknown targets to Sphinx using the new -# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -.DEFAULT: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) +.PHONY: serve +serve: $(TUTORIALS_MD) + @$(MKDOCSBUILD) serve diff --git a/docs/developer_documentation.md b/docs/developer_documentation.md new file mode 100644 index 00000000..a6b4524a --- /dev/null +++ b/docs/developer_documentation.md @@ -0,0 +1,109 @@ +# Developer Documentation + +**Note: this document is a work in progress. Contrbitions to all sections are welcome.** + +## GitHub Workflow + +- You should request reviews for pull requests from `PyRenew-devs`, unless you have a good reason to request reviews from a smaller group. +- Reviews from all of `PyRenew-devs` are encouraged, but an approving review from a [codeowner](https://github.com/CDCgov/PyRenew/blob/main/.github/CODEOWNERS) ([@dylanhmorris](https://github.com/dylanhmorris) or [@damonbayer](https://github.com/damonbayer) is required before a pull request can be merged to `main`. +- For CDC contributors: if your pull request has not received a review at the time of the next standup, use standup to find a reviewer. +- External contributors should expect to receive a review within a few days of creating a pull request. +- If you create a draft pull request, indicate what, if anything, about the current pull request should be reviewed. +- Only mark a pull request as “ready for review” if you think it is ready to be merged. This indicates that a thorough, all-encompassing review should be given. + +## Installation for Developers + +- `uv sync --extra dev` +- `pre-commit install` + +## Coding Conventions + +A variety of coding conventions are enforced by automated tools in continuous integration ([black](https://github.com/psf/black), [isort](https://github.com/PyCQA/isort), [ruff](https://github.com/astral-sh/ruff), [numpydoc-validation](https://github.com/numpy/numpydoc)) via [pre-commit](https://github.com/pre-commit/pre-commit) hooks. + +## PyRenew Principles + +- Variable naming conventions + + - Use the `data_` prefix for (potentially) observed data. + - Use the `_rv` suffix for random variables. + - Use the `observed_` for the output of sample statements where `obs` is a `data_` prefixed object. + - Thus, code which may reasonably written like `infections = infections.sample(x, obs=infections)` should instead be written `observed_infections = infections_rv.sample(x, obs=data_infections)`. + +- Class conventions + + - Composing models is discouraged. + - Returning anything from `Model.sample` is discouraged. Instead, sample from models using `Predictive` or our `prior_predictive` or `posterior_predictive` functions. + - Using `numpyro.deterministic` within a `RandomVariable` is discouraged. Only use at the `numpyro.deterministic` `Model` level. If something might need to be recorded from a `RandomVariable`, it should be returned from the `RandomVariable` so it can be recorded at the `Model` level. + - Using default site names in a `RandomVariable` is discouraged. Only use default site names at the `Model` level. + - Use `DeterministicVariable`s instead of constants within a model. + +- `scan` conventions + + - Use `jax.lax.scan` for any scan whose iterations are deterministic, i.e. iterations contain no internal calls to `RandomVariable.sample()` or `numpyro.sample()`. + - Use `numpyro.scan` for any scan whose the iterations are stochastic, i.e. the iterations potentially include calls to `RandomVariable.sample()` or `numpyro.sample()`. + +- Multidimensional array conventions + + - In a multidimensional array of timeseries, time is always the first dimension. By default, `jax.lax.scan()` and `numpyro.contrib.control_flow.scan()` build output arrays by augmenting the first dimension, and variables are often scanned over time, making default output of scan over time sensible. + +## Adding Documentation to MkDocs + +The project uses MkDocs to generate documentation. To learn more about MkDocs, visit the [MkDocs documentation](https://www.mkdocs.org/). + +### New module or classes +We use Sphinx's `automodule` functionality to generate module API references. When you add a new module, make sure it gets added to the documentation. + +For example, under the `./docs/source/pyrenew_reference`, the `index.md` file lists the `HospitalAdmissions` module by adding the following entry: + +````markdown +## Hospital Admissions + +```{eval-rst} +.. automodule:: pyrenew.latent.hospitaladmissions + :members: + :undoc-members: + :show-inheritance: +``` +```` + +This entry tells Sphinx to generate documentation for the `HospitalAdmissions` module and its members. + +### New tutorials + +`PyRenew` tutorials are [quarto documents](https://quarto.org) located under [./docs/tutorials](https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials). Tutorials are automatically rendered using GitHub actions (see the Workflow file [here](https://github.com/CDCgov/PyRenew/actions/workflows/website.yaml)). + +To make the new tutorial available in the website, developers should follow these steps: + +1. Create a new `quarto` file in the `./docs/tutorials` directory. For instance, the `example_with_datasets.qmd` file was added to the repository. +2. Add an entry in the `mkdocs.yml` navigation section under the Tutorials section, for example: + +````yaml +nav: + - Tutorials: + - tutorials/index.md + - Getting Started: tutorials/getting_started.md + - Example with Datasets: tutorials/example_with_datasets.md +```` + +3. Add an `md` entry with the same basename as the `quarto` file in the `./docs/tutorials` directory. For instance, the `example_with_datasets.md` file was added to the repository. This last step can be done running the bash script [./hook_scripts/pre-commit-md-placeholder.sh](https://github.com/CDCgov/PyRenew/blob/main/hook_scripts/pre-commit-md-placeholder.sh). Note the script should be executed by `pre-commit`. + +### Adding new pages + +MkDocs also allows adding arbitrary pages. For instance, all the `PyRenew` tutorials are additional documentation. The steps to add a new page are: + +1. Create a `md` file in the appropriate directory. For example, this file about development was added under `./docs/developer_documentation.md`. +2. Make sure the new `md` file is included in the navigation section of `mkdocs.yml`. Here is how it looks: + +````yaml +nav: + - Home: index.md + - API Reference: + - reference/index.md + # ... other reference pages + - Tutorials: + - tutorials/index.md + # ... other tutorial pages + - Developer Documentation: developer_documentation.md +```` + +The last entry is the `developer_documentation` page. diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 00000000..948ad00f --- /dev/null +++ b/docs/index.md @@ -0,0 +1,125 @@ +# PyRenew: A Package for Bayesian Renewal Modeling with JAX and NumPyro. + +⚠️ This is a work in progress ⚠️ + +`pyrenew` is a flexible tool for simulation and statistical inference of epidemiological models, emphasizing renewal models. Built on top of the [numpyro](https://num.pyro.ai/) Python library, `pyrenew` provides core components for model building, including pre-defined models for processing various types of observational processes. To start, visit the tutorials section on the project's website [here](https://cdcgov.github.io/PyRenew/tutorials/index.html). + +The following diagram illustrates the composition of the `HospitalAdmissionsModel` class. Notably, all components are modular and can be replaced with custom implementations. + +```mermaid +flowchart LR + + %% Elements + rt_proc["Random Walk Rt\nProcess (latent)"]; + latent_inf["Latent Infections"] + latent_ihr["Infection to Hosp.\nrate (latent)"] + neg_binom["Observation process\n(hospitalizations)"] + latent_hosp["Latent Hospitalizations"]; + i0["Initial infections\n(latent)"]; + gen_int["Generation\ninterval (fixed)"]; + hosp_int["Hospitalization\ninterval (fixed)"]; + + %% Models + basic_model(("Infections\nModel")); + admin_model(("Hospital Admissions\nModel")); + + %% Latent infections + rt_proc --> latent_inf; + i0 --> latent_inf; + gen_int --> latent_inf; + latent_inf --> basic_model + + %% Hospitalizations + hosp_int --> latent_hosp + + neg_binom --> admin_model; + latent_ihr --> latent_hosp; + basic_model --> admin_model; + latent_hosp --> admin_model; +``` + +## Installation + +Install via pip with + +```bash +pip install git+https://github.com/CDCgov/PyRenew@main +``` + +## Models Implemented With PyRenew + +- [CDCgov/pyrenew-covid-wastewater](https://github.com/CDCgov/pyrenew-covid-wastewater): _Models and infrastructure for forecasting COVID-19 hospitalizations using wastewater data with PyRenew._ +- [CDCgov/pyrenew-flu-light](https://github.com/CDCgov/pyrenew-flu-light/): _An instantiation in PyRenew of an influenza forecasting model used in the 2023-24 respiratory season._ + +## Resources + +* [The MSR Website](https://cdcgov.github.io/PyRenew/tutorials/index.html) provides general documentation and tutorials on using MSR. +* [The Model Equations Sheet](https://github.com/CDCgov/PyRenew/blob/main/equations.md) describe the mathematics of the renewal processes and models MSR supports. +* Additional reading on renewal processes in epidemiology + * [_Semi-mechanistic Bayesian modelling of COVID-19 with renewal processes_](https://academic.oup.com/jrsssa/article-pdf/186/4/601/54770289/qnad030.pdf) + * [_Unifying incidence and prevalence under a time-varying general branching process_](https://link.springer.com/content/pdf/10.1007/s00285-023-01958-w.pdf) + +## General Disclaimer + +This repository was created for use by CDC programs to collaborate on public health related projects in support of the [CDC mission](https://www.cdc.gov/about/cdc/index.html#cdc_about_cio_mission-our-mission). GitHub is not hosted by the CDC, but is a third party website used by CDC and its partners to share information and collaborate on software. CDC use of GitHub does not imply an endorsement of any one particular service, product, or enterprise. + +## Public Domain Standard Notice + +This repository constitutes a work of the United States Government and is not +subject to domestic copyright protection under 17 USC § 105. This repository is in +the public domain within the United States, and copyright and related rights in +the work worldwide are waived through the [CC0 1.0 Universal public domain dedication](https://creativecommons.org/publicdomain/zero/1.0/). +All contributions to this repository will be released under the CC0 dedication. By +submitting a pull request you are agreeing to comply with this waiver of +copyright interest. + +## License Standard Notice + +This repository is licensed under ASL v2 or later. + +This source code in this repository is free: you can redistribute it and/or modify it under +the terms of the Apache Software License version 2, or (at your option) any +later version. + +This source code in this repository is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A +PARTICULAR PURPOSE. See the Apache Software License for more details. + +You should have received a copy of the Apache Software License along with this +program. If not, see [http://www.apache.org/licenses/LICENSE-2.0.html](http://www.apache.org/licenses/LICENSE-2.0.html). + +The source code forked from other open source projects will inherit its license. + +## Privacy Standard Notice + +This repository contains only non-sensitive, publicly available data and +information. All material and community participation is covered by the +[Disclaimer](https://github.com/CDCgov/template/blob/master/DISCLAIMER.md) +and [Code of Conduct](https://github.com/CDCgov/template/blob/master/code-of-conduct.md). +For more information about CDC's privacy policy, please visit [http://www.cdc.gov/other/privacy.html](https://www.cdc.gov/other/privacy.html). + +## Contributing Standard Notice + +Anyone is encouraged to contribute to the repository by [forking](https://help.github.com/articles/fork-a-repo) +and submitting a pull request. (If you are new to GitHub, you might start with a +[basic tutorial](https://help.github.com/articles/set-up-git).) By contributing +to this project, you grant a world-wide, royalty-free, perpetual, irrevocable, +non-exclusive, transferable license to all users under the terms of the +[Apache Software License v2](http://www.apache.org/licenses/LICENSE-2.0.html) or +later. + +All comments, messages, pull requests, and other submissions received through +CDC including this GitHub page may be subject to applicable federal law, including but not limited to the Federal Records Act, and may be archived. Learn more at [http://www.cdc.gov/other/privacy.html](http://www.cdc.gov/other/privacy.html). + +## Records Management Standard Notice + +This repository is not a source of government records but is a copy to increase +collaboration and collaborative potential. All government records will be +published through the [CDC web site](http://www.cdc.gov). + +## Additional Standard Notices + +Please refer to [CDC's Template Repository](https://github.com/CDCgov/template) +for more information about [contributing to this repository](https://github.com/CDCgov/template/blob/master/CONTRIBUTING.md), +[public domain notices and disclaimers](https://github.com/CDCgov/template/blob/master/DISCLAIMER.md), +and [code of conduct](https://github.com/CDCgov/template/blob/master/code-of-conduct.md). diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js new file mode 100644 index 00000000..09fbd880 --- /dev/null +++ b/docs/javascripts/mathjax.js @@ -0,0 +1,16 @@ +window.MathJax = { + tex: { + inlineMath: [["\\(", "\\)"]], + displayMath: [["\\[", "\\]"]], + processEscapes: true, + processEnvironments: true + }, + options: { + ignoreHtmlClass: ".*|", + processHtmlClass: "arithmatex" + } +}; + +document$.subscribe(() => { + MathJax.typesetPromise() +}) \ No newline at end of file diff --git a/docs/reference/additional.md b/docs/reference/additional.md new file mode 100644 index 00000000..2d32a660 --- /dev/null +++ b/docs/reference/additional.md @@ -0,0 +1,25 @@ +# Additional modules + +## Metaclass Module + +::: pyrenew.metaclass + +## Convolution Utility Module + +::: pyrenew.convolve + +## Mathematics Utilities Module + +::: pyrenew.math + +## Transformations Module + +::: pyrenew.transformation + +## Regression Module + +::: pyrenew.regression + +## Distributions Utility Module + +::: pyrenew.distutil diff --git a/docs/reference/datasets.md b/docs/reference/datasets.md new file mode 100644 index 00000000..16ed3879 --- /dev/null +++ b/docs/reference/datasets.md @@ -0,0 +1,13 @@ +# Datasets + +## Generation interval + +::: pyrenew.datasets.generation_interval + +## Infection to admission interval + +::: pyrenew.datasets.infection_admission_interval + +## Simulated Wastewater Data + +::: pyrenew.datasets.wastewater diff --git a/docs/reference/deterministic.md b/docs/reference/deterministic.md new file mode 100644 index 00000000..f0bdc8a1 --- /dev/null +++ b/docs/reference/deterministic.md @@ -0,0 +1,3 @@ +# Deterministic + +::: pyrenew.deterministic diff --git a/docs/reference/distributions.md b/docs/reference/distributions.md new file mode 100644 index 00000000..3aa5621d --- /dev/null +++ b/docs/reference/distributions.md @@ -0,0 +1,3 @@ +# Distributions + +::: pyrenew.distributions diff --git a/docs/reference/index.md b/docs/reference/index.md new file mode 100644 index 00000000..0a1a9e01 --- /dev/null +++ b/docs/reference/index.md @@ -0,0 +1,15 @@ +# API Reference + +This section provides detailed documentation for all classes, functions, and modules in the PyRenew package. + +## Modules + +- [Model](model.md) - Core renewal model implementations +- [Latent](latent.md) - Latent variable components +- [Process](process.md) - Stochastic process implementations +- [Random Variable](randomvariable.md) - Random variable abstractions +- [Observation](observation.md) - Observation model components +- [Distributions](distributions.md) - Probability distributions +- [Deterministic](deterministic.md) - Deterministic transformations +- [Datasets](datasets.md) - Built-in datasets +- [Additional](additional.md) - Additional utilities and helpers diff --git a/docs/reference/latent.md b/docs/reference/latent.md new file mode 100644 index 00000000..7ad29c89 --- /dev/null +++ b/docs/reference/latent.md @@ -0,0 +1,21 @@ +# Latent + +## Hospital Admissions + +::: pyrenew.latent.hospitaladmissions + +## Infections + +::: pyrenew.latent.infections + +## Infection Functions + +::: pyrenew.latent.infection_functions + +## Infection Initialization Process + +::: pyrenew.latent.infection_initialization_process + +## Infection Initialization Method + +::: pyrenew.latent.infection_initialization_method diff --git a/docs/reference/model.md b/docs/reference/model.md new file mode 100644 index 00000000..aa28cd4b --- /dev/null +++ b/docs/reference/model.md @@ -0,0 +1,9 @@ +# Model + +## Hospital Admissions + +::: pyrenew.model.admissionsmodel + +## Reproduction Number Renewal Infections + +::: pyrenew.model.rtinfectionsrenewalmodel diff --git a/docs/reference/observation.md b/docs/reference/observation.md new file mode 100644 index 00000000..9cdc11e0 --- /dev/null +++ b/docs/reference/observation.md @@ -0,0 +1,9 @@ +# Observation + +## Negative Binomial + +::: pyrenew.observation.negativebinomial + +## Poisson + +::: pyrenew.observation.poisson diff --git a/docs/reference/process.md b/docs/reference/process.md new file mode 100644 index 00000000..a0610052 --- /dev/null +++ b/docs/reference/process.md @@ -0,0 +1,3 @@ +# Random Process + +::: pyrenew.process diff --git a/docs/reference/randomvariable.md b/docs/reference/randomvariable.md new file mode 100644 index 00000000..14b1e60f --- /dev/null +++ b/docs/reference/randomvariable.md @@ -0,0 +1,3 @@ +# Random Variables + +::: pyrenew.randomvariable diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css new file mode 100644 index 00000000..4b51ffb6 --- /dev/null +++ b/docs/stylesheets/extra.css @@ -0,0 +1,35 @@ +/* Custom styles for PyRenew documentation */ + +/* Code block styling */ +.highlight { + background-color: var(--md-code-bg-color); + border-radius: 0.2rem; +} + +/* Table styling */ +table { + border-collapse: collapse; + width: 100%; +} + +table th, +table td { + border: 1px solid var(--md-default-fg-color--lighter); + padding: 0.5rem; + text-align: left; +} + +table th { + background-color: var(--md-default-bg-color--light); + font-weight: bold; +} + +/* Admonition styling adjustments */ +.admonition { + margin: 1rem 0; +} + +/* MathJax styling */ +.MathJax { + font-size: 1em !important; +} \ No newline at end of file diff --git a/docs/tutorials/basic_renewal_model.qmd b/docs/tutorials/basic_renewal_model.qmd new file mode 100644 index 00000000..040001b2 --- /dev/null +++ b/docs/tutorials/basic_renewal_model.qmd @@ -0,0 +1,407 @@ +--- +title: "Fitting a basic renewal model" +format: gfm +engine: jupyter +--- + +This document will show the steps to build a simple renewal model featuring a latent infection process, a random walk Rt process, and an observation process for the reported infections. + +We start by loading the needed components to build a basic renewal model: + +```{python} +# | label: loading-pkgs +# | output: false +# | warning: false +import jax.numpy as jnp +import numpy as np +import numpyro +import numpyro.distributions as dist +from pyrenew.process import RandomWalk +from pyrenew.latent import ( + Infections, + InfectionInitializationProcess, + InitializeInfectionsZeroPad, +) +from pyrenew.observation import PoissonObservation +from pyrenew.deterministic import DeterministicPMF +from pyrenew.model import RtInfectionsRenewalModel +from pyrenew.metaclass import RandomVariable +from pyrenew.randomvariable import DistributionalVariable, TransformedVariable +import pyrenew.transformation as t +from numpyro.infer.reparam import LocScaleReparam +``` +By default, XLA (which is used by JAX for compilation) considers all CPU cores as one device. Depending on your system's configuration, we recommend using numpyro's [set_host_device_count()](https://num.pyro.ai/en/stable/utilities.html#set-host-device-count) function to set the number of devices available for parallel computing. Here, we set the device count to be 2. +```{python} +# | label: set-device-count +numpyro.set_host_device_count(2) +``` + +## Architecture of `RtInfectionsRenewalModel` + +`pyrenew` leverages `numpyro`'s flexibility to build models via composition. As a principle, most objects in `pyrenew` can be treated as random variables we can sample. At the top-level `pyrenew` has two metaclasses from which most objects inherit: `RandomVariable` and `Model`. From them, the following four sub-modules arise: + +- The `process` sub-module, +- The `deterministic` sub-module, +- The `observation` sub-module, +- The `latent` sub-module, and +- The `models` sub-module + +The first four are collections of instances of `RandomVariable`, and the last is a collection of instances of `Model`. The following diagram shows a detailed view of how metaclasses, modules, and classes interact to create the `RtInfectionsRenewalModel` instantiated in the next section: + +```mermaid +%%| label: overview-of-RtInfectionsRenewalModel +%%| include: true +flowchart LR + rand((RandomVariable\nmetaclass)) + models((Model\nmetaclass)) + + subgraph observations[Observations module] + obs["infection_obs_process_rv\n(PoissonObservation)"] + end + + subgraph latent[Latent module] + inf["latent_infections_rv\n(Infections)"] + i0["I0_rv\n(DistributionalVariable)"] + end + + subgraph process[Process module] + rt["Rt_process_rv\n(Custom class built using RandomWalk)"] + end + + subgraph deterministic[Deterministic module] + detpmf["gen_int_rv\n(DeterministicPMF)"] + end + + subgraph model[Model module] + model1["model1\n(RtInfectionsRenewalModel)"] + end + + rand-->|Inherited by|observations + rand-->|Inherited by|process + rand-->|Inherited by|latent + rand-->|Inherited by|deterministic + models-->|Inherited by|model + + detpmf-->|Composes|model1 + i0-->|Composes|model1 + rt-->|Composes|model1 + obs-->|Composes|model1 + inf-->|Composes|model1 +``` + + +The pyrenew package models the real-time reproductive number $\mathcal{R}(t)$, the average number of secondary infections caused by an infected individual, as a renewal process model. Our basic renewal process model defines five components: + +(1) generation interval, the times between infections + +(2) initial infections, occurring prior to time $t = 0$ + +(3) $\mathcal{R}(t)$, the time-varying reproductive number, + +(4) latent infections, i.e., those infections which are known to exist but are not observed (or not observable), and + +(5) observed infections, a subset of underlying true infections that are reported, perhaps via hospital admissions, physician's office visits, or routine biosurveillance. + +To initialize these five components within the renewal modeling framework, we estimate each component with: + +(1) In this example, the generation interval is not estimated but passed as a deterministic instance of `RandomVariable` + +(2) an instance of the `InfectionInitializationProcess` class, where the number of latent infections immediately before the renewal process begins follows a log-normal distribution with mean = 0 and standard deviation = 1. By specifying `InitializeInfectionsZeroPad`, the latent infections before this time are assumed to be 0. + +(3) A process to represent $\mathcal{R}(t)$ as a random walk on the log scale, with an inferred initial value and a fixed Normal step-size distribution. For this, we construct a custom `RandomVariable`, `MyRt`. + +(4) an instance of the `Infections` class with default values, and + +(5) an instance of the `PoissonObservation` class with default values + +```{python} +# | label: creating-elements +# (1) The generation interval (deterministic) +pmf_array = jnp.array([0.4, 0.3, 0.2, 0.1]) +gen_int = DeterministicPMF(name="gen_int", value=pmf_array) + +# (2) Initial infections (inferred with a prior) +I0 = InfectionInitializationProcess( + "I0_initialization", + DistributionalVariable(name="I0", distribution=dist.LogNormal(2.5, 1)), + InitializeInfectionsZeroPad(pmf_array.size), +) + + +# (3) The random walk on log Rt, with an inferred s.d. Here, we +# construct a custom RandomVariable. +class MyRt(RandomVariable): + def validate(self): + pass + + def sample(self, n: int, **kwargs) -> tuple: + sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) + + rt_rv = TransformedVariable( + name="log_rt_random_walk", + base_rv=RandomWalk( + name="log_rt", + step_rv=DistributionalVariable( + name="rw_step_rv", distribution=dist.Normal(0, sd_rt) + ), + ), + transforms=t.ExpTransform(), + ) + rt_init_rv = DistributionalVariable( + name="init_log_rt", distribution=dist.Normal(0, 0.2) + ) + init_rt = rt_init_rv.sample() + + return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) + + +rt_proc = MyRt() + +# (4) Latent infection process (which will use 1 and 2) +latent_infections = Infections() + +# (5) The observed infections process (with mean at the latent infections) +observation_process = PoissonObservation("poisson_rv") +``` + +With these five pieces, we can build the basic renewal model as an instance of the `RtInfectionsRenewalModel` class: + +```{python} +# | label: model-creation +model1 = RtInfectionsRenewalModel( + gen_int_rv=gen_int, + I0_rv=I0, + Rt_process_rv=rt_proc, + latent_infections_rv=latent_infections, + infection_obs_process_rv=observation_process, +) +``` + +The following diagram summarizes how the modules interact via composition; notably, `gen_int`, `I0`, `rt_proc`, `latent_infections`, and `observed_infections` are instances of `RandomVariable`, which means these can be easily replaced to generate a different instance of the `RtInfectionsRenewalModel` class: + +```mermaid +%%| label: overview-of-RtInfectionsRenewalModel +%%| include: true +flowchart TB + genint["(1) gen_int\n(DetermnisticPMF)"] + i0["(2) I0\n(InfectionInitializationProcess)"] + rt["(3) rt_proc\n(MyRt, the custom RV defined above)"] + inf["(4) latent_infections\n(Infections)"] + obs["(5) observation_process\n(PoissonObservation)"] + + model1["model1\n(RtInfectionsRenewalModel)"] + + i0-->|Composes|model1 + genint-->|Composes|model1 + rt-->|Composes|model1 + obs-->|Composes|model1 + inf-->|Composes|model1 +``` + +Using `numpyro`, we can simulate data using the `sample()` member function of `RtInfectionsRenewalModel`. The `sample()` method of the `RtInfectionsRenewalModel` class returns a list composed of the `Rt` and `infections` sequences, called `sim_data`: + +```{python} +# | label: simulate +with numpyro.handlers.seed(rng_seed=53): + sim_data = model1.sample(n_datapoints=40) + +sim_data +``` + +To understand what has been accomplished here, visualize an $\mathcal{R}(t)$ sample path (left panel) and infections over time (right panel): + +```{python} +# | label: fig-basic +# | fig-cap: Rt and Infections +import matplotlib.pyplot as plt + +fig, axs = plt.subplots(1, 2) + +# Rt plot +axs[0].plot(sim_data.Rt) +axs[0].set_ylabel("Rt") + +# Infections plot +axs[1].plot(sim_data.observed_infections) +axs[1].set_ylabel("Infections") + +fig.suptitle("Basic renewal model") +fig.supxlabel("Time") +plt.tight_layout() +plt.show() +``` + +To fit the model, we can use the `run()` method of the `RtInfectionsRenewalModel` class (an inherited method from the metaclass `Model`). `model1.run()` will call the `run` method of the `model1` object, which will generate an instance of model MCMC simulation, with 2000 warm-up iterations for the MCMC algorithm, used to tune the parameters of the MCMC algorithm to improve efficiency of the sampling process. From the posterior distribution of the model parameters, 1000 samples will be drawn and used to estimate the posterior distributions and compute summary statistics. Observed data is provided to the model using the `sim_data` object previously generated. `mcmc_args` provides additional arguments for the MCMC algorithm. + +```{python} +# | label: model-fit +import jax + +model1.run( + num_warmup=2000, + num_samples=1000, + data_observed_infections=sim_data.observed_infections, + rng_key=jax.random.key(54), + mcmc_args=dict(progress_bar=False, num_chains=2), +) +``` + +Now, let's investigate the output, particularly the posterior distribution of the $\mathcal{R}(t)$ estimates: + +```{python} +# | label: fig-output-rt +# | fig-cap: Rt posterior distribution +import arviz as az + +# Create arviz inference data object +idata = az.from_numpyro( + posterior=model1.mcmc, +) + +# Extract Rt signal samples across chains +rt = az.extract(idata.posterior["Rt"], num_samples=100)["Rt"].values + + +# Plot Rt signal +fig, ax = plt.subplots(1, 1, figsize=(8, 6)) +ax.plot( + np.arange(rt.shape[0]), + rt, + color="skyblue", + alpha=0.10, +) +ax.plot([], [], color="skyblue", alpha=0.05, label="Rt Posterior Samples") +ax.plot( + np.arange(rt.shape[0]), + rt.mean(axis=1), + color="black", + linewidth=2.0, + linestyle="--", + label="Sample Mean", +) +ax.legend(loc="best") +ax.set_ylabel(r"$\mathscr{R}_t$ Signal", fontsize=20) +ax.set_xlabel("Days", fontsize=20) +plt.show() +``` +We can use the `get_samples` method to extract samples from the model +```{python} +# | label: get-model-fitted-samples +Rt_samp = model1.mcmc.get_samples()["Rt"] +latent_infection_samp = model1.mcmc.get_samples()["all_latent_infections"] +``` + +We can also convert the fitted model to [ArviZ](https://www.arviz.org/) InferenceData object and use ArviZ package to extarct samples, calculate statistics, create model diagnostics and visualizations. + +```{python} +# | label: convert-inference-data +import arviz as az + +idata = az.from_numpyro(model1.mcmc) +``` + +and use the InferenceData to compute the model-fit diagnostics. Here, we show diagnostic summary for the first 10 effective reproduction number $\mathcal{R}(t)$. + +```{python} +# | label: diagnostics +diagnostic_stats_summary = az.summary( + idata.posterior["Rt"][::, ::, 4:], # ignore nan padding + kind="diagnostics", +) + +print(diagnostic_stats_summary) +``` + +Below we use `plot_trace` to inspect the trace of the first 10 inferred $\mathcal{R}(t)$ values. + +```{python} +# | label: fig-trace-Rt +# | fig-cap: Trace plot of Rt posterior distribution +plt.rcParams["figure.constrained_layout.use"] = True + +az.plot_trace( + idata.posterior, + var_names=["Rt"], + coords={"Rt_dim_0": np.arange(4, 14)}, + compact=False, +) +plt.show() +``` + + +We inspect the posterior distribution of $\mathcal{R}(t)$ by plotting the 90% and 50% highest density intervals: + +```{python} +# | label: fig-hdi-Rt +# | fig-cap: High density interval for Effective Reproduction Number +x_data = idata.posterior["Rt_dim_0"][4:] +y_data = idata.posterior["Rt"][::, ::, 4:] + +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.9, + color="C0", + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.5, + color="C0", + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +# Add mean of the posterior to the figure +median_ts = y_data.median(dim=["chain", "draw"]) +axes.plot(x_data, median_ts, color="C0", label="Median") +axes.legend() +axes.set_title("Posterior Effective Reproduction Number", fontsize=10) +axes.set_xlabel("Time", fontsize=10) +axes.set_ylabel("$\\mathcal{R}(t)$", fontsize=10) +plt.show() +``` + +and latent infections: + +```{python} +# | label: fig-hdi-latent-infections +# | fig-cap: High density interval for Latent Infections +x_data = idata.posterior["all_latent_infections_dim_0"] +y_data = idata.posterior["all_latent_infections"] + +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.9, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.5, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +# plot the posterior median +median_ts = y_data.median(dim=["chain", "draw"]) +axes.plot(x_data, median_ts, color="C0", label="Median") + +axes.legend() +axes.set_title("Posterior Latent Infections", fontsize=10) +axes.set_xlabel("Time", fontsize=10) +axes.set_ylabel("Latent Infections", fontsize=10) +plt.show() +``` diff --git a/docs/tutorials/day_of_the_week.qmd b/docs/tutorials/day_of_the_week.qmd new file mode 100644 index 00000000..018ce4e4 --- /dev/null +++ b/docs/tutorials/day_of_the_week.qmd @@ -0,0 +1,426 @@ +--- +title: Implementing a Day of the Week Effect +format: gfm +engine: jupyter +--- + +This document illustrates how to leverage the time-aware arrays to create a a day of the week effect. We use the same model designed in the hospital admissions-only tutorial. + +## Recap: Hospital Admissions Model + +In the ["Fitting a hospital admissions-only model" tutorial](https://cdcgov.github.io/PyRenew/tutorials/hospital_admissions_model.html){target="_blank"}, we built a fairly complicated model that included pre-defined random variables as well as a custom random variable representing the reproductive number. In this tutorial, we will focus on adding a new component: the day-of-the-week effect. We start by reproducing the model without the day-of-the-week effect: + +1. We load the data: + +```{python} +# | label: setup +# | code-fold: true +# Setup +import numpyro +import polars as pl +from pyrenew import datasets + +# Setting the number of devices +numpyro.set_host_device_count(2) + +# Loading and processing the data +dat = ( + datasets.load_wastewater() + .group_by("date") + .first() + .select(["date", "daily_hosp_admits"]) + .sort("date") + .head(90) +) + +daily_hosp_admits = dat["daily_hosp_admits"].to_numpy() +dates = dat["date"].to_numpy() + +# Loading additional datasets +gen_int = datasets.load_generation_interval() +inf_hosp_int = datasets.load_infection_admission_interval() + +# We only need the probability_mass column of each dataset +gen_int_array = gen_int["probability_mass"].to_numpy() +gen_int = gen_int_array +inf_hosp_int_array = inf_hosp_int["probability_mass"].to_numpy() +``` + +2. Next, we defined the model's components: + +```{python} +# | label: latent-hosp +# | code-fold: true +from pyrenew import latent, deterministic, randomvariable +import jax.numpy as jnp +import numpyro.distributions as dist + +inf_hosp_int = deterministic.DeterministicPMF( + name="inf_hosp_int", value=inf_hosp_int_array +) + +hosp_rate = randomvariable.DistributionalVariable( + name="IHR", distribution=dist.LogNormal(jnp.log(0.05), jnp.log(1.1)) +) + +latent_hosp = latent.HospitalAdmissions( + infection_to_admission_interval_rv=inf_hosp_int, + infection_hospitalization_ratio_rv=hosp_rate, +) + +from pyrenew import model, process, observation, metaclass, transformation +from pyrenew.latent import ( + InfectionInitializationProcess, + InitializeInfectionsExponentialGrowth, +) + + +# Infection process +latent_inf = latent.Infections() +n_initialization_points = max(gen_int_array.size, inf_hosp_int_array.size) - 1 + +I0 = InfectionInitializationProcess( + "I0_initialization", + randomvariable.DistributionalVariable( + name="I0", + distribution=dist.LogNormal(loc=jnp.log(100), scale=jnp.log(1.75)), + ), + InitializeInfectionsExponentialGrowth( + n_initialization_points, + deterministic.DeterministicVariable(name="rate", value=0.05), + ), +) + +# Generation interval and Rt +gen_int = deterministic.DeterministicPMF(name="gen_int", value=gen_int) +``` + +including the Rt effect: + +```{python} +# | label: Rt-process +# | code-fold: true +class MyRt(metaclass.RandomVariable): + def __init__(self, sd_rv): + self.sd_rv = sd_rv + + def validate(self): + pass + + def sample(self, n: int, **kwargs) -> tuple: + # Standard deviation of the random walk + sd_rt = self.sd_rv() + + # Random walk step + step_rv = randomvariable.DistributionalVariable( + name="rw_step_rv", distribution=dist.Normal(0, sd_rt) + ) + + rt_init_rv = randomvariable.DistributionalVariable( + name="init_log_rt", distribution=dist.Normal(0, 0.2) + ) + + # Random walk process + base_rv = process.RandomWalk( + name="log_rt", + step_rv=step_rv, + ) + + # Transforming the random walk to the Rt scale + rt_rv = randomvariable.TransformedVariable( + name="Rt_rv", + base_rv=base_rv, + transforms=transformation.ExpTransform(), + ) + init_rt = rt_init_rv.sample() + + return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) + + +rtproc = MyRt( + randomvariable.DistributionalVariable( + name="Rt_random_walk_sd", distribution=dist.HalfNormal(0.025) + ) +) +``` + +3. We defined the observation model: + +```{python} +# | label: obs-model +# | code-fold: true +# we place a log-Normal prior on the concentration +# parameter of the negative binomial. +nb_conc_rv = randomvariable.TransformedVariable( + "concentration", + randomvariable.DistributionalVariable( + name="concentration_raw", + distribution=dist.TruncatedNormal(loc=0, scale=1, low=0.01), + ), + transformation.PowerTransform(-2), +) + +# now we define the observation process +obs = observation.NegativeBinomialObservation( + "negbinom_rv", + concentration_rv=nb_conc_rv, +) +``` + +4. And finally, we build the model: + +```{python} +# | label: init-model +# | code-fold: true +hosp_model = model.HospitalAdmissionsModel( + latent_infections_rv=latent_inf, + latent_hosp_admissions_rv=latent_hosp, + I0_rv=I0, + gen_int_rv=gen_int, + Rt_process_rv=rtproc, + hosp_admission_obs_process_rv=obs, +) +``` + +Here is what the model looks like without the day-of-the-week effect: + +```{python} +# | label: fig-output-admissions-padding-and-weekday +# | fig-cap: Hospital Admissions posterior distribution without weekday effect +# | code-fold: true +import jax +import numpy as np + +# Model without weekday effect +hosp_model.run( + num_samples=2000, + num_warmup=2000, + data_observed_hosp_admissions=daily_hosp_admits, + rng_key=jax.random.key(54), + mcmc_args=dict(progress_bar=False), +) +``` + +```{python} +# | code-fold: true +import arviz as az +import matplotlib.pyplot as plt + + +# Retrieve the posterior samples from the model +ppc_samples = hosp_model.posterior_predictive( + n_datapoints=daily_hosp_admits.size +) + +# Create an InferenceData object from model +idata = az.from_numpyro( + posterior=hosp_model.mcmc, + posterior_predictive=ppc_samples, +) + +# Use a time series plot (plot_ts) from arviz for plotting +axes = az.plot_ts( + idata, + y="negbinom_rv", + y_hat="negbinom_rv", + num_samples=200, + y_kwargs={ + "color": "blue", + "linewidth": 1.0, + "marker": "o", + "linestyle": "solid", + }, + y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, + y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, + backend_kwargs={"figsize": (8, 6)}, + textsize=15.0, +) +ax = axes[0][0] +ax.set_xlabel("Time", fontsize=20) +ax.set_ylabel("Hospital Admissions", fontsize=20) +handles, labels = ax.get_legend_handles_labels() +ax.legend( + handles, ["Observed", "Sample Mean", "Posterior Samples"], loc="best" +) +plt.show() +``` + + + + + +## Round 2: Incorporating day-of-the-week effects + +We will re-use the infection to admission interval and infection to hospitalization rate from the previous model. But we will also add a day-of-the-week effect. To do this, we will add two additional arguments to the latent hospital admissions random variable: `day_of_the_week_rv` (a `RandomVariable`) and `obs_data_first_day_of_the_week` (an `int` mapping days of the week from 0:6, zero being Monday). The `day_of_the_week_rv`'s sample method should return a vector of length seven; those values are then broadcasted to match the length of the dataset. Moreover, since the observed data may start in a weekday other than Monday, the `obs_data_first_day_of_the_week` argument is used to offset the day-of-the-week effect. + +For this example, the effect will be passed as a scaled Dirichlet distribution. It will consist of a `TransformedVariable` that samples an array of length seven from numpyro's `distributions.Dirichlet` and applies a `transformation.AffineTransform` to scale it by seven. [^note-other-examples]: + +[^note-other-examples]: A similar weekday effect is implemented in its own module, with example code here: {doc}`periodic_effects` + +```{python} +# | label: weekly-effect +# Instantiating the day-of-the-week effect +dayofweek_effect = randomvariable.TransformedVariable( + name="dayofweek_effect", + base_rv=randomvariable.DistributionalVariable( + name="dayofweek_effect_raw", + distribution=dist.Dirichlet(jnp.ones(7)), + ), + transforms=transformation.AffineTransform( + loc=0, scale=7, domain=jnp.array([0, 1]) + ), +) +``` + +Now, by re-defining the latent hospital admissions random variable with the day-of-the-week effect, we can build a model that includes this effect. Since the day-of-the-week effect takes into account the first day of the dataset, we need to determine the day of the week of the first observation. We can do this by converting the first date in the dataset to a `datetime` object and extracting the day of the week: + +```{python} +# | label: latent-hosp-weekday +# Figuring out the day of the week of the first observation +import datetime as dt + +first_dow_in_data = dates[0].astype(dt.datetime).weekday() +first_dow_in_data # zero + +# Re-defining the latent hospital admissions RV, now with the +# day-of-the-week effect +latent_hosp_wday_effect = latent.HospitalAdmissions( + infection_to_admission_interval_rv=inf_hosp_int, + infection_hospitalization_ratio_rv=hosp_rate, + day_of_week_effect_rv=dayofweek_effect, + # Concidirently, this is zero + obs_data_first_day_of_the_week=first_dow_in_data, +) + +# New model with day-of-the-week effect +hosp_model_dow = model.HospitalAdmissionsModel( + latent_infections_rv=latent_inf, + latent_hosp_admissions_rv=latent_hosp_wday_effect, + I0_rv=I0, + gen_int_rv=gen_int, + Rt_process_rv=rtproc, + hosp_admission_obs_process_rv=obs, +) +``` + +Running the model: + +```{python} +# | label: model-2-run-weekday +# Model with weekday effect +hosp_model_dow.run( + num_samples=2000, + num_warmup=2000, + data_observed_hosp_admissions=daily_hosp_admits, + rng_key=jax.random.key(54), + mcmc_args=dict(progress_bar=False), +) +``` + +As a result, we can see the posterior distribution of our novel day-of-the-week effect: + +```{python} +# | label: fig-output-day-of-week +# | fig-cap: Day of the week effect +# Create an InferenceData object from hosp_model_dow +dow_idata = az.from_numpyro( + posterior=hosp_model_dow.mcmc, +) + +# Extract day of week effect (DOW) +dow_effect_raw = dow_idata.posterior["dayofweek_effect_raw"].squeeze().T +indices = np.random.choice(dow_effect_raw.shape[1], size=200, replace=False) +dow_plot_samples = dow_effect_raw[:, indices] +fig, ax = plt.subplots(1, 1, figsize=(8, 6)) +ax.plot( + np.arange(dow_effect_raw.shape[0]), + dow_plot_samples, + color="skyblue", + alpha=0.10, +) +ax.plot([], [], color="skyblue", alpha=0.10, label="DOW Posterior Samples") +ax.plot( + np.arange(dow_effect_raw.shape[0]), + dow_plot_samples.mean(dim="draw"), + color="black", + linewidth=2.0, + linestyle="--", + label="Sample Mean", +) +ax.legend(loc="best") +ax.set_ylabel("Effect", fontsize=20) +ax.set_xlabel("Day Of Week", fontsize=20) +plt.show() +``` + +The new model with the day-of-the-week effect can be compared to the previous model without the effect. Finally, let's reproduce the figure without the day-of-the-week effect, and then plot the new model with the effect: + +```{python} +# | label: fig-output-admissions-original +# | fig-cap: Hospital Admissions posterior distribution without weekday effect +# Without weekday effect (from earlier) +axes = az.plot_ts( + idata, + y="negbinom_rv", + y_hat="negbinom_rv", + num_samples=200, + y_kwargs={ + "color": "blue", + "linewidth": 1.0, + "marker": "o", + "linestyle": "solid", + }, + y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, + y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, + backend_kwargs={"figsize": (8, 6)}, + textsize=15.0, +) +ax = axes[0][0] +ax.set_xlabel("Time", fontsize=20) +ax.set_ylabel("Hospital Admissions", fontsize=20) +handles, labels = ax.get_legend_handles_labels() +ax.legend( + handles, + ["Observed", "Posterior Predictive", "Samples wo/ WDE"], + loc="best", +) +plt.show() +``` + +```{python} +# | label: fig-output-admissions-wof +# | fig-cap: Hospital Admissions posterior distribution with weekday effect +# Figure with weekday effect +ppc_samples = hosp_model_dow.posterior_predictive( + n_datapoints=daily_hosp_admits.size +) +idata = az.from_numpyro( + posterior=hosp_model_dow.mcmc, + posterior_predictive=ppc_samples, +) + +axes = az.plot_ts( + idata, + y="negbinom_rv", + y_hat="negbinom_rv", + num_samples=200, + y_kwargs={ + "color": "blue", + "linewidth": 1.0, + "marker": "o", + "linestyle": "solid", + }, + y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, + y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, + backend_kwargs={"figsize": (8, 6)}, + textsize=15.0, +) +ax = axes[0][0] +ax.set_xlabel("Time", fontsize=20) +ax.set_ylabel("Hospital Admissions", fontsize=20) +handles, labels = ax.get_legend_handles_labels() +ax.legend( + handles, ["Observed", "Posterior Predictive", "Samples w/ WDE"], loc="best" +) +plt.show() +``` diff --git a/docs/tutorials/extending_pyrenew.qmd b/docs/tutorials/extending_pyrenew.qmd new file mode 100644 index 00000000..77a83a53 --- /dev/null +++ b/docs/tutorials/extending_pyrenew.qmd @@ -0,0 +1,297 @@ +--- +title: Extending pyrenew +format: gfm +--- + +This tutorial illustrates how to extend `pyrenew` with custom `RandomVariable` classes. We will use the `InfectionsWithFeedback` class as an example. The `InfectionsWithFeedback` class is a `RandomVariable` that models the number of infections at time $t$ as a function of the number of infections at time $t - \tau$ and the reproduction number at time $t$. The reproduction number at time $t$ is a function of the *unadjusted* reproduction number at time $t - \tau$ and the number of infections at time $t - \tau$: + +$$ +\begin{align*} +I(t) & = \mathcal{R}(t)\sum_{\tau=1}^{T_g}I(t - \tau)g(\tau) \\ +\mathcal{R}(t) & = \mathcal{R}^u(t)\exp\left(-\gamma(t)\sum_{\tau=1}^{T_f}I(t - \tau)f(\tau)\right) +\end{align*} +$$ + +Where $\mathcal{R}^u(t)$ is the unadjusted reproduction number, $g(t)$ is the generation interval, $\gamma(t)$ is the infection feedback strength, and $f(t)$ is the infection feedback pmf. + +## The expected outcome + +Before we start, let's simulate the model with the original `InfectionsWithFeedback` class. To keep it simple, we will simulate the model with no observation process, in other words, only with latent infections. The following code-chunk loads the required libraries and defines the model components: + +```{python} +# | label: setup +import jax +import jax.numpy as jnp +import numpy as np +import numpyro +import numpyro.distributions as dist +from pyrenew.deterministic import DeterministicPMF, DeterministicVariable +from pyrenew.latent import InfectionsWithFeedback +from pyrenew.model import RtInfectionsRenewalModel +from pyrenew.process import RandomWalk +from pyrenew.metaclass import RandomVariable +from pyrenew.randomvariable import DistributionalVariable, TransformedVariable +from pyrenew.latent import ( + InfectionInitializationProcess, + InitializeInfectionsExponentialGrowth, +) +import pyrenew.transformation as t +``` + +The following code-chunk defines the model components. Notice that for both the generation interval and the infection feedback, we use a deterministic PMF with equal probabilities: + +```{python} +# | label: model-components +gen_int_array = jnp.array([0.25, 0.5, 0.15, 0.1]) + +gen_int = DeterministicPMF(name="gen_int", value=gen_int_array) +feedback_strength = DeterministicVariable(name="feedback_strength", value=0.01) + + +I0 = InfectionInitializationProcess( + "I0_initialization", + DistributionalVariable(name="I0", distribution=dist.LogNormal(0, 1)), + InitializeInfectionsExponentialGrowth( + gen_int_array.size, + DeterministicVariable(name="rate", value=0.05), + ), +) + +latent_infections = InfectionsWithFeedback( + infection_feedback_strength=feedback_strength, + infection_feedback_pmf=gen_int, +) + + +class MyRt(RandomVariable): + def validate(self): + pass + + def sample(self, n: int, **kwargs) -> tuple: + sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) + + rt_rv = TransformedVariable( + name="log_rt_random_walk", + base_rv=RandomWalk( + name="log_rt", + step_rv=DistributionalVariable( + name="rw_step_rv", distribution=dist.Normal(0, 0.025) + ), + ), + transforms=t.ExpTransform(), + ) + rt_init_rv = DistributionalVariable( + name="init_log_rt", distribution=dist.Normal(0, 0.2) + ) + init_rt = rt_init_rv.sample() + + return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) +``` + +With all the components defined, we can build the model: + +```{python} +# | label: build1 +model0 = RtInfectionsRenewalModel( + gen_int_rv=gen_int, + I0_rv=I0, + latent_infections_rv=latent_infections, + Rt_process_rv=MyRt(), + infection_obs_process_rv=None, +) +``` + +And simulate from it: + +```{python} +# | label: simulate1 +# Sampling and fitting model 0 (with no obs for infections) +with numpyro.handlers.seed(rng_seed=223): + model0_samp = model0.sample(n_datapoints=30) +``` + +```{python} +# | label: fig-simulate1 +# | fig-cap: Simulated infections with no observation process +import matplotlib.pyplot as plt + +fig, ax = plt.subplots() +ax.plot(model0_samp.latent_infections) +ax.set_xlabel("Time") +ax.set_ylabel("Infections") +plt.show() +``` + +## Pyrenew's random variable class + +### Fundamentals + +All instances of PyRenew's `RandomVariable` should have at least three functions: `__init__()`, `validate()`, and `sample()`. The `__init__()` function is the constructor and initializes the class. The `validate()` function checks if the class is correctly initialized. Finally, the `sample()` method contains the core of the class; it should return a tuple or named tuple. The following is a minimal example of a `RandomVariable` class based on `numpyro.distributions.Normal`: + +```{python} +from pyrenew.metaclass import RandomVariable + + +class MyNormal(RandomVariable): + def __init__(self, loc, scale): + self.validate(scale) + self.loc = loc + self.scale = scale + return None + + @staticmethod + def validate(self): + if self.scale <= 0: + raise ValueError("Scale must be positive") + return None + + def sample(self, **kwargs): + return (dist.Normal(loc=self.loc, scale=self.scale),) +``` + +The `@staticmethod` decorator exposes the `validate` function to be used outside the class. Next, we show how to build a more complex `RandomVariable` class; the `InfectionsWithFeedback` class. + +### The `InfectionsWithFeedback` class + +Although returning namedtuples is not strictly required, they are the recommended return type, as they make the code more readable. The following code-chunk shows how to create a named tuple for the `InfectionsWithFeedback` class: + +```{python} +# | label: data-class +from collections import namedtuple + +# Creating a tuple to store the output +InfFeedbackSample = namedtuple( + typename="InfFeedbackSample", + field_names=["post_initialization_infections", "rt"], + defaults=(None, None), +) +``` + +The next step is to create the actual class. The bulk of its implementation lies in the function `pyrenew.latent.compute_infections_from_rt_with_feedback()`. We will also use the `pyrenew.arrayutils.pad_edges_to_match()` function to ensure the passed vectors match their lengths. The following code-chunk shows most of the implementation of the `InfectionsWithFeedback` class: + +```{python} +# | label: new-model-def +# | code-line-numbers: true +# Creating the class +from pyrenew.metaclass import RandomVariable +from pyrenew.latent import compute_infections_from_rt_with_feedback +from pyrenew import arrayutils as au +from jax.typing import ArrayLike +import jax.numpy as jnp + + +class InfFeedback(RandomVariable): + """Latent infections""" + + def __init__( + self, + infection_feedback_strength: RandomVariable, + infection_feedback_pmf: RandomVariable, + ) -> None: + """Constructor""" + + self.infection_feedback_strength = infection_feedback_strength + self.infection_feedback_pmf = infection_feedback_pmf + + return None + + def validate(self): + """ + Generally, this method should be more meaningful, but we will skip it for now + """ + return None + + def sample( + self, + Rt: ArrayLike, + I0: ArrayLike, + gen_int: ArrayLike, + **kwargs, + ) -> tuple: + """Sample infections with feedback""" + + # Generation interval + gen_int_rev = jnp.flip(gen_int) + + # Baseline infections + I0_vec = I0[-gen_int_rev.size :] + + # Sampling inf feedback strength and adjusting the shape + inf_feedback_strength = self.infection_feedback_strength( + **kwargs, + ) + + inf_feedback_strength = jnp.atleast_1d(inf_feedback_strength) + + inf_feedback_strength, _ = au.pad_edges_to_match( + x=inf_feedback_strength, + y=Rt, + ) + + # Sampling inf feedback and adjusting the shape + inf_feedback_pmf = self.infection_feedback_pmf(**kwargs) + inf_fb_pmf_rev = jnp.flip(inf_feedback_pmf) + + # Generating the infections with feedback + all_infections, Rt_adj = compute_infections_from_rt_with_feedback( + I0=I0_vec, + Rt_raw=Rt, + infection_feedback_strength=inf_feedback_strength, + reversed_generation_interval_pmf=gen_int_rev, + reversed_infection_feedback_pmf=inf_fb_pmf_rev, + ) + + # Storing adjusted Rt for future use + numpyro.deterministic("Rt_adjusted", Rt_adj) + + # Preparing theoutput + + return InfFeedbackSample( + post_initialization_infections=all_infections, + rt=Rt_adj, + ) +``` + +The core of the class is implemented in the `sample()` method. Things to highlight from the above code: + +1. **Arguments of `sample`**: The `InfFeedback` class will be used within `RtInfectionsRenewalModel` to generate latent infections. During the sampling process, `InfFeedback()` will receive the reproduction number, the initial number of infections, and the generation interval. `RandomVariable()` calls are expected to include the `**kwargs` argument, even if unused. + +2. **Saving computed quantities**: Since `Rt_adj` is not generated via `numpyro.sample()`, we use `numpyro.deterministic()` to record the quantity to a site; allowing us to access it later. + +3. **Return type of `InfFeedback()`**: As said before, the `sample()` method should return a tuple or named tuple. In our case, we return a named tuple `InfFeedbackSample` with two fields: `infections` and `rt`. + +```{python} +# | label: simulation2 +latent_infections2 = InfFeedback( + infection_feedback_strength=feedback_strength, + infection_feedback_pmf=gen_int, +) + +model1 = RtInfectionsRenewalModel( + gen_int_rv=gen_int, + I0_rv=I0, + latent_infections_rv=latent_infections2, + Rt_process_rv=MyRt(), + infection_obs_process_rv=None, +) + +# Sampling and fitting model 0 (with no obs for infections) +with numpyro.handlers.seed(rng_seed=223): + model1_samp = model1.sample(n_datapoints=30) +``` + +Comparing `model0` with `model1`, these two should match: + +```{python} +# | label: fig-model0-vs-model1 +# | fig-cap: Comparing latent infections from model 0 and model 1 +import matplotlib.pyplot as plt + +fig, ax = plt.subplots(ncols=2) +ax[0].plot(model0_samp.latent_infections) +ax[1].plot(model1_samp.latent_infections) +ax[0].set_xlabel("Time (model 0)") +ax[1].set_xlabel("Time (model 1)") +ax[0].set_ylabel("Infections") +plt.show() +``` diff --git a/docs/tutorials/getting_started.qmd b/docs/tutorials/getting_started.qmd new file mode 100644 index 00000000..24490f91 --- /dev/null +++ b/docs/tutorials/getting_started.qmd @@ -0,0 +1,60 @@ +--- +title: "Getting started with pyrenew" +format: gfm +engine: jupyter +--- + +`pyrenew` is a flexible tool for simulating and making statistical inferences from epidemiologic models, with an emphasis on renewal models. Built on `numpyro`, `pyrenew` provides core components for model building and pre-defined models for processing various observational processes. + +## Prerequisites + +This tutorial assumes some pre-existing knowledge of infectious disease dynamics and Python programming. Before you dive in, we recommend: + + - Installing Python3 (use tools like [pyenv](https://realpython.com/intro-to-pyenv/) or [compile and install](https://ubuntuhandbook.org/index.php/2023/05/install-python-3-12-ubuntu/) from the [release page](https://www.python.org/downloads/)) + - Familiarity with installing and loading modules in python, and with virtual environment management (we recommend [uv](https://docs.astral.sh/uv/)) + - Familiarity with the concept of a [class](https://realpython.com/python-classes/) and [metaclass](https://realpython.com/courses/python-metaclasses/) in python + - Familiarity with Bayesian inference, and a working understanding of MCMC methods used to fit Bayesian models to data (some resources are available [here](https://mc-stan.org/docs/2_18/reference-manual/effective-sample-size-section.html), and [here](https://xcelab.net/rm/)) + +## Installing pyrenew + +You’ll need to install `pyrenew` using either uv or pip. To install `pyrenew` using `uv`, run the following command from within the directory containing the `pyrenew` project: + +```bash +uv sync +``` + +To install `pyrenew` using pip, run the following command: + +```bash +pip install git+https://github.com/CDCgov/PyRenew@main +``` + + +## The fundamentals + +`pyrenew`'s core components are the metaclasses `RandomVariable` and `Model` (in Python, a _metaclass_ is a class whose instances are also classes, where a _class_ is a template for making objects). Within the `pyrenew` package, a `RandomVariable` is a quantity that models can estimate and sample from, **including deterministic quantities**. The benefit of this design is that the definition of the `sample()` function can be arbitrary, allowing the user to either sample from a distribution using `numpyro.sample()`, compute fixed quantities (like a mechanistic equation), or return a fixed value (like a pre-computed PMF.) For instance, when estimating a PMF, the `RandomVariable` sampling function may roughly be defined as: + +```python +# define a new class called MyRandVar that inherits from the RandomVariable class +class MyRandVar(RandomVariable): + #define a method called sample that returns an object of type ArrayLike + def sample(...) -> ArrayLike: + # calls sample function from NumPyro package + return numpyro.sample(...) +``` + +Whereas, in some other cases, we may instead use a fixed quantity for that variable (like a pre-computed PMF), where the `RandomVariable`'s sample function could instead be defined as: + +```python +# instead define MyRandVar to still inherit from the RandVariable class +class MyRandVar(RandomVariable): + #define sample method that still returns an ArrayLike object + def sample(...) -> ArrayLike: + #sampling method is a pre-computed PMF, a JAX NumPy array with explicit elements + return jax.numpy.array([0.2, 0.7, 0.1]) +``` + +Thus, when a `Model` samples from `MyRandVar`, it could be either adding random variables to be estimated (first case) or just retrieving some quantity needed for other calculations (second case.) + +The `Model` metaclass provides basic functionality for estimating and simulation. Like `RandomVariable`, the `Model` metaclass has a `sample()` method that defines the model structure. Ultimately, models can be nested (or inherited), providing a straightforward way to add layers of complexity. +At this stage, the `Model` metaclass consist of two model classes `RtInfectionsRenewalModel` which is basic renewal model consisting of infections and reproduction numbers and `HospitalAdmissionsModel` which includes basic renewal model and hospital admissions. In the subsequent sections, we provide examples of fitting each of these models. diff --git a/docs/tutorials/hospital_admissions_model.qmd b/docs/tutorials/hospital_admissions_model.qmd new file mode 100644 index 00000000..dc3aee8f --- /dev/null +++ b/docs/tutorials/hospital_admissions_model.qmd @@ -0,0 +1,591 @@ +--- +title: Fitting a hospital admissions-only model +format: gfm +engine: jupyter +--- + +This document illustrates how a hospital admissions-only model can be fitted using data from the Pyrenew package, particularly the wastewater dataset. The CFA wastewater team created this dataset, which contains simulated data. + +We begin by loading `numpyro` and configuring the device count to 2 to enable running MCMC chains in parallel. By default, XLA (which is used by JAX for compilation) considers all CPU cores as one device. Depending on your system's configuration, we recommend using numpyro's [set_host_device_count()](https://num.pyro.ai/en/stable/utilities.html#set-host-device-count) function to set the number of devices available for parallel computing. + +```{python} +# | label: numpyro setup +import numpyro + +numpyro.set_host_device_count(2) +``` +## Model definition + +In this section, we provide the formal definition of the model. The hospital admissions model is a semi-mechanistic model that describes the number of observed hospital admissions as a function of a set of latent variables. Mainly, the observed number of hospital admissions is discretely distributed with location at the number of latent hospital admissions: + +$$ +h(t) \sim \text{HospDist}\left(H(t)\right) +$$ + +Where $h(t)$ is the observed number of hospital admissions at time $t$, and $H(t)$ is the number of latent hospital admissions at time $t$. The distribution $\text{HospDist}$ is discrete. For this example, we will use a negative binomial distribution with an inferred concentration. + +$$ +\begin{align*} +h(t) & \sim \mathrm{NegativeBinomial}\left(\mathrm{mean} = H(t), \mathrm{concentration} = k\right) \\ +H(t) & = p_\mathrm{hosp}(t) \sum_{\tau = 0}^{T_d} d(\tau) I(t-\tau) \\ +\log[p_\mathrm{hosp}(t)] & \sim \mathrm{Normal}(\mu=\log(0.05), \sigma=\log(1.1)) \\ +\log(k) & \sim \mathrm{Normal}(\mu=\log(1), \sigma=\log(10)) +\end{align*} +$$ + +Where $d(\tau)$ is the infection to hospital admission interval, $I(t)$ is the number of latent infections at time $t$, and $p_\mathrm{hosp}(t)$ is the infection to admission rate. + +The number of latent hospital admissions at time $t$ is a function of the number of latent infections at time $t$ and the infection to admission rate. The latent infections are modeled as a renewal process: + +$$ +\begin{align*} +I(t) &= \mathcal{R}(t) \times \sum_{\tau < t} I(\tau) g(t - \tau) \\ +\log[I(0)] &\sim \text{Normal}(\mu=\log(100), \sigma=\log(1.75)) +\end{align*} +$$ + +The reproduction number $\mathcal{R}(t)$ is modeled as a random walk in logarithmic space, i.e.: + +$$ +\begin{align*} +\log[\mathcal{R}(t)] & = \log[\mathcal{R}(t-1)] + \epsilon\\ +\epsilon & \sim \text{Normal}(\mu=0, \sigma=0.025) \\ +\mathcal{R}(0) &\sim \text{TruncatedNormal}(\text{loc}=1.2, \text{scale}=0.2, \text{min}=0) +\end{align*} +$$ + + +## Data processing + +We start by loading the data and inspecting the first five rows. + +```{python} +# | label: data-inspect +import polars as pl +from pyrenew import datasets + +dat = datasets.load_wastewater() +dat.head(5) +``` + +The data shows one entry per site, but the way it was simulated, the number of admissions is the same across sites. Thus, we will only keep the first observation per day. + +```{python} +# | label: aggregation +# Keeping the first observation of each date +dat = dat.group_by("date").first().select(["date", "daily_hosp_admits"]) + +# Now, sorting by date +dat = dat.sort("date") + +# Keeping the first 90 days +dat = dat.head(90) + +dat.head(5) +``` + +Let's take a look at the daily prevalence of hospital admissions. + +```{python} +# | label: fig-plot-hospital-admissions +# | fig-cap: Daily hospital admissions from the simulated data +import matplotlib.pyplot as plt +import matplotlib.dates as mdates + +daily_hosp_admits = dat["daily_hosp_admits"].to_numpy() +dates = dat["date"].to_numpy() +ax = plt.gca() +ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d")) +ax.xaxis.set_major_locator(mdates.DayLocator(interval=7)) +ax.set_xlim(dates[0], dates[-1]) +plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") +plt.plot(dates, daily_hosp_admits, "-o") +plt.xlabel("Date") +plt.ylabel("Admissions") +plt.show() +``` + +## Building the model + +First, we will extract two datasets we will use as deterministic quantities: the generation interval and the infection to hospital admission interval. + +```{python} +# | label: fig-data-extract +# | fig-cap: Generation interval and infection to hospital admission interval +gen_int = datasets.load_generation_interval() +inf_hosp_int = datasets.load_infection_admission_interval() + +# We only need the probability_mass column of each dataset +gen_int_array = gen_int["probability_mass"].to_numpy() +gen_int = gen_int_array +inf_hosp_int_array = inf_hosp_int["probability_mass"].to_numpy() + +# Taking a peek at the first 5 elements of each +gen_int[:5], inf_hosp_int_array[:5] + +# Visualizing both quantities side by side +fig, axs = plt.subplots(1, 2) + +axs[0].plot(gen_int) +axs[0].set_title("Generation interval") +axs[1].plot(inf_hosp_int_array) +axs[1].set_title("Infection to hospital admission interval") +plt.show() +``` + +With these two in hand, we can start building the model. First, we will define the latent hospital admissions: + +```{python} +# | label: latent-hosp +from pyrenew import latent, deterministic, metaclass, randomvariable +import jax.numpy as jnp +import numpyro.distributions as dist + +inf_hosp_int = deterministic.DeterministicPMF( + name="inf_hosp_int", value=inf_hosp_int_array +) + +hosp_rate = randomvariable.DistributionalVariable( + name="IHR", distribution=dist.LogNormal(jnp.log(0.05), jnp.log(1.1)) +) + +latent_hosp = latent.HospitalAdmissions( + infection_to_admission_interval_rv=inf_hosp_int, + infection_hospitalization_ratio_rv=hosp_rate, +) +``` + +The `inf_hosp_int` is a `DeterministicPMF` object that takes the infection to hospital admission interval as input. The `hosp_rate` is a `DistributionalVariable` object that takes a numpyro distribution to represent the infection to hospital admission rate. The `HospitalAdmissions` class is a `RandomVariable` that takes two distributions as inputs: the infection to admission interval and the infection to hospital admission rate. Now, we can define the rest of the other components: + +```{python} +# | label: initializing-rest-of-model +from pyrenew import model, process, observation, metaclass, transformation +from pyrenew.latent import ( + InfectionInitializationProcess, + InitializeInfectionsExponentialGrowth, +) + + +# Infection process +latent_inf = latent.Infections() +n_initialization_points = max(gen_int_array.size, inf_hosp_int_array.size) - 1 +I0 = InfectionInitializationProcess( + "I0_initialization", + randomvariable.DistributionalVariable( + name="I0", + distribution=dist.LogNormal(loc=jnp.log(100), scale=jnp.log(1.75)), + ), + InitializeInfectionsExponentialGrowth( + n_initialization_points, + deterministic.DeterministicVariable(name="rate", value=0.05), + ), +) + +# Generation interval and Rt +gen_int = deterministic.DeterministicPMF(name="gen_int", value=gen_int) + + +class MyRt(metaclass.RandomVariable): + def validate(self): + pass + + def sample(self, n: int, **kwargs) -> tuple: + sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) + + rt_rv = randomvariable.TransformedVariable( + name="log_rt_random_walk", + base_rv=process.RandomWalk( + name="log_rt", + step_rv=randomvariable.DistributionalVariable( + name="rw_step_rv", distribution=dist.Normal(0, 0.025) + ), + ), + transforms=transformation.ExpTransform(), + ) + rt_init_rv = randomvariable.DistributionalVariable( + name="init_log_rt", distribution=dist.Normal(0, 0.2) + ) + init_rt = rt_init_rv.sample() + + return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) + + +rtproc = MyRt() + +# The observation model + +# we place a log-Normal prior on the concentration +# parameter of the negative binomial. +nb_conc_rv = randomvariable.TransformedVariable( + "concentration", + randomvariable.DistributionalVariable( + name="concentration_raw", + distribution=dist.TruncatedNormal(loc=0, scale=1, low=0.01), + ), + transformation.PowerTransform(-2), +) + +# now we define the observation process +obs = observation.NegativeBinomialObservation( + "negbinom_rv", + concentration_rv=nb_conc_rv, +) +``` + +Notice all the components are `RandomVariable` instances. We can now build the model: + +```{python} +# | label: init-model +hosp_model = model.HospitalAdmissionsModel( + latent_infections_rv=latent_inf, + latent_hosp_admissions_rv=latent_hosp, + I0_rv=I0, + gen_int_rv=gen_int, + Rt_process_rv=rtproc, + hosp_admission_obs_process_rv=obs, +) +``` + +Let's simulate from the prior predictive distribution to check that the model is working: + +```{python} +# | label: simulation +import numpy as np + +timeframe = 120 + +with numpyro.handlers.seed(rng_seed=223): + simulated_data = hosp_model.sample(n_datapoints=timeframe) +``` + +```{python} +# | label: fig-basic +# | fig-cap: Simulated Rt and Admissions +import matplotlib.pyplot as plt + +fig, axs = plt.subplots(1, 2) + +# Rt plot +axs[0].plot(simulated_data.Rt) +axs[0].set_ylabel("Simulated Rt") + +# Admissions plot +axs[1].plot(simulated_data.observed_hosp_admissions, "-o") +axs[1].set_ylabel("Simulated Admissions") + +fig.suptitle("Basic renewal model") +fig.supxlabel("Time") +plt.tight_layout() +plt.show() +``` + +## Fitting the model + +We now fit the model, not to these simulated data, but rather to the dataset we retrieved above. We use the `run` method of the `Model` object: + +```{python} +# | label: model-fit +import jax + +hosp_model.run( + num_samples=1000, + num_warmup=1000, + data_observed_hosp_admissions=daily_hosp_admits, + rng_key=jax.random.key(54), + mcmc_args=dict(progress_bar=False, num_chains=2), +) +``` + +We can use `arviz` to visualize the model fit. Here, we plot the observed values against the inferred latent values (i.e. the mean of the negative binomial observation process)[^capture]: + +[^capture]: The output is captured to avoid `quarto` from displaying the output twice. + +```{python} +# | label: fig-output-hospital-admissions +# | fig-cap: Latent hospital admissions posterior samples (gray) and observed admissions timeseries (red). +import arviz as az + +ppc_samples = hosp_model.posterior_predictive( + n_datapoints=daily_hosp_admits.size +) +idata = az.from_numpyro( + posterior=hosp_model.mcmc, + posterior_predictive=ppc_samples, +) + +axes = az.plot_ts( + idata, + y="negbinom_rv", + y_hat="negbinom_rv", + num_samples=200, + y_kwargs={ + "color": "blue", + "linewidth": 1.0, + "marker": "o", + "linestyle": "solid", + }, + y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, + y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, + backend_kwargs={"figsize": (8, 6)}, + textsize=15.0, +) +ax = axes[0][0] +ax.set_xlabel("Time", fontsize=20) +ax.set_ylabel("Hospital Admissions", fontsize=20) +handles, labels = ax.get_legend_handles_labels() +ax.legend( + handles, ["Observed", "Sample Mean", "Posterior Samples"], loc="best" +) +plt.show() +``` + +## Results exploration and MCMC diagnostics + +To explore further, We can use [ArviZ](https://www.arviz.org/) to visualize the results. Let's start by loading the module and converting the fitted model to ArviZ `InferenceData` object: + +```{python} +# | label: convert-inferenceData + +idata = az.from_numpyro(hosp_model.mcmc) +``` + +We obtain the summary of model diagnostics and print the diagnostics for `latent_hospital_admissions[1]` + +```{python} +# | label: diagnostics +# | warning: false +diagnostic_stats_summary = az.summary( + idata.posterior, + kind="diagnostics", +) + +print(diagnostic_stats_summary.loc["latent_hospital_admissions[1]"]) +``` + +Below we plot 90% and 50% highest density intervals for latent hospital admissions using [plot_hdi](https://python.arviz.org/en/stable/api/generated/arviz.plot_hdi.html): + +```{python} +# | label: fig-output-admission-distribution +# | fig-cap: Hospital Admissions posterior distribution +x_data = idata.posterior["latent_hospital_admissions_dim_0"] +y_data = idata.posterior["latent_hospital_admissions"] + + +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.9, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.5, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +# Add the posterior median to the figure +median_ts = y_data.median(dim=["chain", "draw"]) + +axes.plot(x_data, median_ts, color="C0", label="Median") +axes.legend() +axes.set_title("Posterior Hospital Admissions", fontsize=10) +axes.set_xlabel("Time", fontsize=10) +axes.set_ylabel("Hospital Admissions", fontsize=10) +plt.show() +``` + +We can also look at credible intervals for the posterior distribution of latent infections: + +```{python} +# | label: fig-output-infections-distribution +# | fig-cap: Posterior Latent Infections +x_data = ( + idata.posterior["all_latent_infections_dim_0"] - n_initialization_points +) +y_data = idata.posterior["all_latent_infections"] + +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.9, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + x_data, + y_data, + hdi_prob=0.5, + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +# Add the posterior median to the figure +median_ts = y_data.median(dim=["chain", "draw"]) +axes.plot(x_data, median_ts, color="C0", label="Median") +axes.legend() +``` + + +## Predictive checks and forecasting +We can use the `Model`'s `posterior_predictive` and `prior_predictive` methods to generate posterior and prior predictive samples for observed admissions. + +```{python} +# | label: demonstrate-use-of-predictive-methods + +idata = az.from_numpyro( + hosp_model.mcmc, + posterior_predictive=hosp_model.posterior_predictive( + n_datapoints=len(daily_hosp_admits) + ), + prior=hosp_model.prior_predictive( + n_datapoints=len(daily_hosp_admits), + numpyro_predictive_args={"num_samples": 1000}, + ), +) +``` + +We will use `plot_lm` method from ArviZ to plot the posterior predictive distribution against the actual observed data below: + +```{python} +# | label: fig-posterior-predictive +# | fig-cap: Hospital Admissions posterior distribution with plot_lm +fig, ax = plt.subplots() +az.plot_lm( + "negbinom_rv", + idata=idata, + kind_pp="hdi", + y_kwargs={"color": "black"}, + y_hat_fill_kwargs={"color": "C0"}, + axes=ax, +) + +ax.set_title("Posterior Predictive Plot") +ax.set_ylabel("Hospital Admissions") +ax.set_xlabel("Days") +plt.show() +``` + + +By increasing `n_datapoints`, we can perform forecasting using the posterior predictive distribution. + +```{python} +# | label: posterior-predictive-distribution +n_forecast_points = 28 +idata = az.from_numpyro( + hosp_model.mcmc, + posterior_predictive=hosp_model.posterior_predictive( + n_datapoints=len(daily_hosp_admits) + n_forecast_points, + ), + prior=hosp_model.prior_predictive( + n_datapoints=len(daily_hosp_admits), + numpyro_predictive_args={"num_samples": 1000}, + ), +) +``` + +Below we plot the prior predictive distributions using equal tailed Bayesian credible intervals: + +```{python} +# | label: fig-output-prior-predictive +# | fig-cap: Prior Predictive Admissions +def compute_eti(dataset, eti_prob): + eti_bdry = dataset.quantile( + ((1 - eti_prob) / 2, 1 / 2 + eti_prob / 2), dim=("chain", "draw") + ) + return eti_bdry.T + + +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + idata.prior_predictive["negbinom_rv_dim_0"], + hdi_data=compute_eti(idata.prior_predictive["negbinom_rv"], 0.9), + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + idata.prior_predictive["negbinom_rv_dim_0"], + hdi_data=compute_eti(idata.prior_predictive["negbinom_rv"], 0.5), + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +plt.scatter( + idata.observed_data["negbinom_rv_dim_0"], + idata.observed_data["negbinom_rv"], + color="black", +) + +axes.set_title("Prior Predictive Admissions", fontsize=10) +axes.set_xlabel("Time", fontsize=10) +axes.set_ylabel("Observed Admissions", fontsize=10) +plt.yscale("log") +plt.show() +``` + +And now we plot the posterior predictive distributions with a `{python} n_forecast_points`-day-ahead forecast: +```{python} +# | label: fig-output-posterior-predictive-forecast +# | fig-cap: Posterior predictive admissions, including a forecast. +x_data = idata.posterior_predictive["negbinom_rv_dim_0"] +y_data = idata.posterior_predictive["negbinom_rv"] +fig, axes = plt.subplots(figsize=(6, 5)) +az.plot_hdi( + x_data, + hdi_data=compute_eti(y_data, 0.9), + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.3}, + ax=axes, +) + +az.plot_hdi( + x_data, + hdi_data=compute_eti(y_data, 0.5), + color="C0", + smooth=False, + fill_kwargs={"alpha": 0.6}, + ax=axes, +) + +# Add median of the posterior to the figure +median_ts = y_data.median(dim=["chain", "draw"]) + +plt.plot( + x_data, + median_ts, + color="C0", + label="Median", +) +plt.scatter( + idata.observed_data["negbinom_rv_dim_0"], + idata.observed_data["negbinom_rv"], + color="black", +) +axes.legend() +axes.set_title( + "Posterior Predictive Admissions, including a forecast", fontsize=10 +) +axes.set_xlabel("Time", fontsize=10) +axes.set_ylabel("Hospital Admissions", fontsize=10) +plt.show() +``` diff --git a/docs/tutorials/index.md b/docs/tutorials/index.md new file mode 100644 index 00000000..b6c1672e --- /dev/null +++ b/docs/tutorials/index.md @@ -0,0 +1,12 @@ +# Tutorials + +This section contains tutorials that demonstrate how to use the `pyrenew` package. The source code for the tutorials can be found in the project repository: . + +## Tutorial List + +- [Getting Started](getting_started.md) - Introduction to PyRenew +- [Basic Renewal Model](basic_renewal_model.md) - Building your first renewal model +- [Hospital Admissions Model](hospital_admissions_model.md) - Using real-world hospital data +- [Extending PyRenew](extending_pyrenew.md) - Adding custom components +- [Periodic Effects](periodic_effects.md) - Modeling periodic patterns +- [Day of the Week Effects](day_of_the_week.md) - Weekly seasonality modeling diff --git a/docs/tutorials/periodic_effects.qmd b/docs/tutorials/periodic_effects.qmd new file mode 100644 index 00000000..7f3cf5b7 --- /dev/null +++ b/docs/tutorials/periodic_effects.qmd @@ -0,0 +1,105 @@ +--- +title: Periodic effects and broadcasting +format: gfm +engine: jupyter +--- + +A common practice in time series forecasting is featuring periodic effects such as day of the week or weekend effects. `pyrenew` has built-in functions and classes to incorporate such effects into renewal models. Notably, the `arrayutils.PeriodicBroadcaster` class. + +This document demonstrates how to use these functions and classes to model periodic effects and broadcasting in renewal models: Repeating elements of a sequence and repeating a sequence as a whole (tiling). + +## Repeated units + +The `RtPeriodicDiff` and `RtWeeklyDiff` classes use `PeriodicBroadcaster` to repeat each vector element to represent Rt values held constant within a period. Let's show it's usage with a basic renewal model. We will use the same model structure shown in the `getting_started` tutorial: + +```{python} +# | label: loading-pkgs +# | output: false +# | warning: false +import jax.numpy as jnp +import numpy as np +import numpyro +from pyrenew import process, deterministic +``` + +```{python} +# The random process for Rt +rt_proc = process.RtWeeklyDiffARProcess( + name="rt_weekly_diff", + offset=0, + log_rt_rv=deterministic.DeterministicVariable( + name="log_rt", value=jnp.array([0.1, 0.2]) + ), + autoreg_rv=deterministic.DeterministicVariable( + name="autoreg", value=jnp.array([0.7]) + ), + periodic_diff_sd_rv=deterministic.DeterministicVariable( + name="periodic_diff_sd", value=jnp.array([0.1]) + ), +) +``` + +```{python} +with numpyro.handlers.seed(rng_seed=20): + sim_data = rt_proc(duration=30) + +# Plotting the Rt values +import matplotlib.pyplot as plt + +plt.step(np.arange(len(sim_data)), sim_data, where="post") +plt.xlabel("Time") +plt.ylabel("Rt") +plt.title("Simulated Rt values") + +# Adding bands to mark weeks +for i in range(0, 30, 7): + plt.axvline(i, color="black", linestyle="--", alpha=0.5) +plt.show() +``` + +The implementation of the `RtWeeklyDiffARProcess` (which is an instance of `RtPeriodicDiffARProcess`), uses `repeat_until_n` to repeating values: `repeat_until_n(..., period_size=7)`. The `RtWeeklyDiff` class is a particular case of `RtPeriodicDiff` with a period size of seven. + +## Repeated sequences (tiling) + +The `PeriodicBroadcaster` class can also be used to repeat a sequence as a whole. For example, we can use this to model a periodic effect such as day-of-the-week effect. Let's use the `DayOfWeekEffect` random variable. We will sample from a scaled Dirchlet distribution such that the sum of the samples is 7: + +```{python} +import numpyro.distributions as dist +from pyrenew import transformation, randomvariable + +# Building the transformed prior: Dirichlet * 7 +mysimplex = dist.TransformedDistribution( + dist.Dirichlet(concentration=jnp.ones(7)), + transformation.AffineTransform(loc=0, scale=7.0), +) + +# Constructing the day of week effect +dayofweek = process.DayOfWeekEffect( + offset=0, + quantity_to_broadcast=randomvariable.DistributionalVariable( + name="simp", distribution=mysimplex + ), +) +``` + +Like before, we can use the `sample` method to generate samples from the day of week effect: + +```{python} +with numpyro.handlers.seed(rng_seed=20): + sim_data = dayofweek(duration=30) + +# Plotting the effect values +import matplotlib.pyplot as plt + +plt.step(np.arange(len(sim_data)), sim_data, where="post") +plt.xlabel("Time") +plt.ylabel("Effect size") +plt.title("Simulated Day of Week Effect values") + +# Adding bands to mark weeks +for i in range(0, 30, 7): + plt.axvline(i, color="black", linestyle="--", alpha=0.5) +plt.show() +``` + +Like the example from the previous section, the `DayOfWeekEffect` class uses `PeriodicBroadcaster` to broadcast the data, this time using `broadcast_type = "tile"`, repeating the sequence as a whole. The `DayOfWeekEffect` class is a particular case of `PeriodicEffect` with a period size of seven. diff --git a/hook_scripts/pre-commit-md-placeholder.sh b/hook_scripts/pre-commit-md-placeholder.sh index 61e4dc04..aa64d7f4 100755 --- a/hook_scripts/pre-commit-md-placeholder.sh +++ b/hook_scripts/pre-commit-md-placeholder.sh @@ -1,9 +1,9 @@ #!/usr/bin/env sh -TUTORIALS_DIR=docs/source/tutorials +TUTORIALS_DIR=docs/tutorials -# For each *qmd file in docs/source/tutorials, create a -# corresponding md file under docs/source/tutorials +# For each *qmd file in docs/tutorials, create a +# corresponding md file under docs/tutorials COUNTER=0 # Removing everything under docs/source/tutorials/*.md with diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 00000000..beb9907f --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,110 @@ +site_name: "Multisignal Renewal Models for Epi Inference" +site_description: "PyRenew, a Python package for multi-signal Bayesian renewal modeling with JAX and NumPyro." +site_url: "https://cdcgov.github.io/PyRenew/" +repo_url: "https://github.com/CDCgov/PyRenew" +repo_name: "CDCgov/PyRenew" +edit_uri: "edit/main/docs/" + +nav: + - Home: index.md + - API Reference: + - reference/index.md + - Model: reference/model.md + - Latent: reference/latent.md + - Process: reference/process.md + - Random Variable: reference/randomvariable.md + - Observation: reference/observation.md + - Distributions: reference/distributions.md + - Deterministic: reference/deterministic.md + - Datasets: reference/datasets.md + - Additional: reference/additional.md + - Tutorials: + - tutorials/index.md + - Getting Started: tutorials/getting_started.md + - Basic Renewal Model: tutorials/basic_renewal_model.md + - Day of the Week Effects: tutorials/day_of_the_week.md + - Extending PyRenew: tutorials/extending_pyrenew.md + - Hospital Admissions Model: tutorials/hospital_admissions_model.md + - Periodic Effects: tutorials/periodic_effects.md + - Developer Documentation: developer_documentation.md + +theme: + name: material + icon: + repo: fontawesome/brands/github + features: + - navigation.tabs + - navigation.sections + - navigation.expand + - navigation.top + - search.highlight + - search.suggest + - content.code.copy + - content.action.edit + palette: + # Palette toggle for light mode + - scheme: default + primary: blue + accent: blue + toggle: + icon: material/brightness-7 + name: Switch to dark mode + # Palette toggle for dark mode + - scheme: slate + primary: blue + accent: blue + toggle: + icon: material/brightness-4 + name: Switch to light mode + +plugins: + - search + - mkdocstrings: + handlers: + python: + options: + docstring_style: numpy + show_source: true + show_root_heading: true + show_root_full_path: false + show_object_full_path: false + show_category_heading: true + show_if_no_docstring: false + annotations_path: brief + line_length: 80 + show_signature_annotations: true + separate_signature: true + +markdown_extensions: + - admonition + - pymdownx.details + - pymdownx.superfences: + custom_fences: + - name: mermaid + class: mermaid + format: !!python/name:pymdownx.superfences.fence_code_format + - pymdownx.highlight: + anchor_linenums: true + line_spans: __span + pygments_lang_class: true + - pymdownx.inlinehilite + - pymdownx.snippets + - pymdownx.tabbed: + alternate_style: true + - attr_list + - md_in_html + - toc: + permalink: true + - pymdownx.arithmatex: + generic: true + +extra_javascript: + - javascripts/mathjax.js + - https://polyfill.io/v3/polyfill.min.js?features=es6 + - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + +extra_css: + - stylesheets/extra.css + +watch: + - pyrenew/ \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 029f2cbf..b93d62b1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,14 +27,11 @@ docs = [ "nbformat>=5.10.0", "numpydoc>=1.7.0", "arviz>=0.19.0", - "sphinx>=7.2.6", - "sphinxcontrib-mermaid>=0.9.2", - "sphinx-autodoc-typehints>=2.1.0", - "sphinx-book-theme>=1.1.2", + "mkdocs>=1.5.0", + "mkdocs-material>=9.0.0", + "mkdocstrings[python]>=0.24.0", "ipywidgets>=8.1.3", - "jupyter>=1.0.0", - "myst-parser>=4.0.0", - "sphinxcontrib-katex>=0.9.10" + "jupyter>=1.0.0" ] test = [ "pytest>=8.3.2", diff --git a/uv.lock b/uv.lock index 7aa3a626..e0fa725b 100644 --- a/uv.lock +++ b/uv.lock @@ -1,23 +1,11 @@ version = 1 -revision = 2 +revision = 3 requires-python = ">=3.13" resolution-markers = [ "python_full_version >= '3.14'", "python_full_version < '3.14'", ] -[[package]] -name = "accessible-pygments" -version = "0.0.5" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "pygments" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/bc/c1/bbac6a50d02774f91572938964c582fff4270eee73ab822a4aeea4d8b11b/accessible_pygments-0.0.5.tar.gz", hash = "sha256:40918d3e6a2b619ad424cb91e556bd3bd8865443d9f22f1dcdf79e33c8046872", size = 1377899, upload-time = "2024-05-10T11:23:10.216Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/8d/3f/95338030883d8c8b91223b4e21744b04d11b161a3ef117295d8241f50ab4/accessible_pygments-0.0.5-py3-none-any.whl", hash = "sha256:88ae3211e68a1d0b011504b2ffc1691feafce124b845bd072ab6f9f66f34d4b7", size = 1395903, upload-time = "2024-05-10T11:23:08.421Z" }, -] - [[package]] name = "alabaster" version = "1.0.0" @@ -190,6 +178,20 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/b7/b8/3fe70c75fe32afc4bb507f75563d39bc5642255d1d94f1f23604725780bf/babel-2.17.0-py3-none-any.whl", hash = "sha256:4d0b53093fdfb4b21c92b5213dba5a1b23885afa8383709427046b21c366e5f2", size = 10182537, upload-time = "2025-02-01T15:17:37.39Z" }, ] +[[package]] +name = "backrefs" +version = "5.9" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/eb/a7/312f673df6a79003279e1f55619abbe7daebbb87c17c976ddc0345c04c7b/backrefs-5.9.tar.gz", hash = "sha256:808548cb708d66b82ee231f962cb36faaf4f2baab032f2fbb783e9c2fdddaa59", size = 5765857, upload-time = "2025-06-22T19:34:13.97Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/19/4d/798dc1f30468134906575156c089c492cf79b5a5fd373f07fe26c4d046bf/backrefs-5.9-py310-none-any.whl", hash = "sha256:db8e8ba0e9de81fcd635f440deab5ae5f2591b54ac1ebe0550a2ca063488cd9f", size = 380267, upload-time = "2025-06-22T19:34:05.252Z" }, + { url = "https://files.pythonhosted.org/packages/55/07/f0b3375bf0d06014e9787797e6b7cc02b38ac9ff9726ccfe834d94e9991e/backrefs-5.9-py311-none-any.whl", hash = "sha256:6907635edebbe9b2dc3de3a2befff44d74f30a4562adbb8b36f21252ea19c5cf", size = 392072, upload-time = "2025-06-22T19:34:06.743Z" }, + { url = "https://files.pythonhosted.org/packages/9d/12/4f345407259dd60a0997107758ba3f221cf89a9b5a0f8ed5b961aef97253/backrefs-5.9-py312-none-any.whl", hash = "sha256:7fdf9771f63e6028d7fee7e0c497c81abda597ea45d6b8f89e8ad76994f5befa", size = 397947, upload-time = "2025-06-22T19:34:08.172Z" }, + { url = "https://files.pythonhosted.org/packages/10/bf/fa31834dc27a7f05e5290eae47c82690edc3a7b37d58f7fb35a1bdbf355b/backrefs-5.9-py313-none-any.whl", hash = "sha256:cc37b19fa219e93ff825ed1fed8879e47b4d89aa7a1884860e2db64ccd7c676b", size = 399843, upload-time = "2025-06-22T19:34:09.68Z" }, + { url = "https://files.pythonhosted.org/packages/fc/24/b29af34b2c9c41645a9f4ff117bae860291780d73880f449e0b5d948c070/backrefs-5.9-py314-none-any.whl", hash = "sha256:df5e169836cc8acb5e440ebae9aad4bf9d15e226d3bad049cf3f6a5c20cc8dc9", size = 411762, upload-time = "2025-06-22T19:34:11.037Z" }, + { url = "https://files.pythonhosted.org/packages/41/ff/392bff89415399a979be4a65357a41d92729ae8580a66073d8ec8d810f98/backrefs-5.9-py39-none-any.whl", hash = "sha256:f48ee18f6252b8f5777a22a00a09a85de0ca931658f1dd96d4406a34f3748c60", size = 380265, upload-time = "2025-06-22T19:34:12.405Z" }, +] + [[package]] name = "beautifulsoup4" version = "4.13.5" @@ -551,6 +553,30 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/cf/58/8acf1b3e91c58313ce5cb67df61001fc9dcd21be4fadb76c1a2d540e09ed/fqdn-1.5.1-py3-none-any.whl", hash = "sha256:3a179af3761e4df6eb2e026ff9e1a3033d3587bf980a0b1b2e1e5d08d7358014", size = 9121, upload-time = "2021-03-11T07:16:28.351Z" }, ] +[[package]] +name = "ghp-import" +version = "2.1.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "python-dateutil" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/d9/29/d40217cbe2f6b1359e00c6c307bb3fc876ba74068cbab3dde77f03ca0dc4/ghp-import-2.1.0.tar.gz", hash = "sha256:9c535c4c61193c2df8871222567d7fd7e5014d835f97dc7b7439069e2413d343", size = 10943, upload-time = "2022-05-02T15:47:16.11Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/f7/ec/67fbef5d497f86283db54c22eec6f6140243aae73265799baaaa19cd17fb/ghp_import-2.1.0-py3-none-any.whl", hash = "sha256:8337dd7b50877f163d4c0289bc1f1c7f127550241988d568c1db512c4324a619", size = 11034, upload-time = "2022-05-02T15:47:14.552Z" }, +] + +[[package]] +name = "griffe" +version = "1.14.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "colorama" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/ec/d7/6c09dd7ce4c7837e4cdb11dce980cb45ae3cd87677298dc3b781b6bce7d3/griffe-1.14.0.tar.gz", hash = "sha256:9d2a15c1eca966d68e00517de5d69dd1bc5c9f2335ef6c1775362ba5b8651a13", size = 424684, upload-time = "2025-09-05T15:02:29.167Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/2a/b1/9ff6578d789a89812ff21e4e0f80ffae20a65d5dd84e7a17873fe3b365be/griffe-1.14.0-py3-none-any.whl", hash = "sha256:0e9d52832cccf0f7188cfe585ba962d2674b241c01916d780925df34873bceb0", size = 144439, upload-time = "2025-09-05T15:02:27.511Z" }, +] + [[package]] name = "h11" version = "0.16.0" @@ -1121,15 +1147,12 @@ wheels = [ ] [[package]] -name = "markdown-it-py" -version = "3.0.0" +name = "markdown" +version = "3.9" source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "mdurl" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/38/71/3b932df36c1a044d397a1f92d1cf91ee0a503d91e470cbd670aa66b07ed0/markdown-it-py-3.0.0.tar.gz", hash = "sha256:e3f60a94fa066dc52ec76661e37c851cb232d92f9886b15cb560aaada2df8feb", size = 74596, upload-time = "2023-06-03T06:41:14.443Z" } +sdist = { url = "https://files.pythonhosted.org/packages/8d/37/02347f6d6d8279247a5837082ebc26fc0d5aaeaf75aa013fcbb433c777ab/markdown-3.9.tar.gz", hash = "sha256:d2900fe1782bd33bdbbd56859defef70c2e78fc46668f8eb9df3128138f2cb6a", size = 364585, upload-time = "2025-09-04T20:25:22.885Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/42/d7/1ec15b46af6af88f19b8e5ffea08fa375d433c998b8a7639e76935c14f1f/markdown_it_py-3.0.0-py3-none-any.whl", hash = "sha256:355216845c60bd96232cd8d8c40e8f9765cc86f46880e43a8fd22dc1a1a8cab1", size = 87528, upload-time = "2023-06-03T06:41:11.019Z" }, + { url = "https://files.pythonhosted.org/packages/70/ae/44c4a6a4cbb496d93c6257954260fe3a6e91b7bed2240e5dad2a717f5111/markdown-3.9-py3-none-any.whl", hash = "sha256:9f4d91ed810864ea88a6f32c07ba8bee1346c0cc1f6b1f9f6c822f2a9667d280", size = 107441, upload-time = "2025-09-04T20:25:21.784Z" }, ] [[package]] @@ -1220,33 +1243,141 @@ wheels = [ ] [[package]] -name = "mdit-py-plugins" -version = "0.5.0" +name = "mergedeep" +version = "1.3.4" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/3a/41/580bb4006e3ed0361b8151a01d324fb03f420815446c7def45d02f74c270/mergedeep-1.3.4.tar.gz", hash = "sha256:0096d52e9dad9939c3d975a774666af186eda617e6ca84df4c94dec30004f2a8", size = 4661, upload-time = "2021-02-05T18:55:30.623Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/2c/19/04f9b178c2d8a15b076c8b5140708fa6ffc5601fb6f1e975537072df5b2a/mergedeep-1.3.4-py3-none-any.whl", hash = "sha256:70775750742b25c0d8f36c55aed03d24c3384d17c951b3175d898bd778ef0307", size = 6354, upload-time = "2021-02-05T18:55:29.583Z" }, +] + +[[package]] +name = "mistune" +version = "3.1.4" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d7/02/a7fb8b21d4d55ac93cdcde9d3638da5dd0ebdd3a4fed76c7725e10b81cbe/mistune-3.1.4.tar.gz", hash = "sha256:b5a7f801d389f724ec702840c11d8fc48f2b33519102fc7ee739e8177b672164", size = 94588, upload-time = "2025-08-29T07:20:43.594Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/7a/f0/8282d9641415e9e33df173516226b404d367a0fc55e1a60424a152913abc/mistune-3.1.4-py3-none-any.whl", hash = "sha256:93691da911e5d9d2e23bc54472892aff676df27a75274962ff9edc210364266d", size = 53481, upload-time = "2025-08-29T07:20:42.218Z" }, +] + +[[package]] +name = "mkdocs" +version = "1.6.1" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "markdown-it-py" }, + { name = "click" }, + { name = "colorama", marker = "sys_platform == 'win32'" }, + { name = "ghp-import" }, + { name = "jinja2" }, + { name = "markdown" }, + { name = "markupsafe" }, + { name = "mergedeep" }, + { name = "mkdocs-get-deps" }, + { name = "packaging" }, + { name = "pathspec" }, + { name = "pyyaml" }, + { name = "pyyaml-env-tag" }, + { name = "watchdog" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/b2/fd/a756d36c0bfba5f6e39a1cdbdbfdd448dc02692467d83816dff4592a1ebc/mdit_py_plugins-0.5.0.tar.gz", hash = "sha256:f4918cb50119f50446560513a8e311d574ff6aaed72606ddae6d35716fe809c6", size = 44655, upload-time = "2025-08-11T07:25:49.083Z" } +sdist = { url = "https://files.pythonhosted.org/packages/bc/c6/bbd4f061bd16b378247f12953ffcb04786a618ce5e904b8c5a01a0309061/mkdocs-1.6.1.tar.gz", hash = "sha256:7b432f01d928c084353ab39c57282f29f92136665bdd6abf7c1ec8d822ef86f2", size = 3889159, upload-time = "2024-08-30T12:24:06.899Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/fb/86/dd6e5db36df29e76c7a7699123569a4a18c1623ce68d826ed96c62643cae/mdit_py_plugins-0.5.0-py3-none-any.whl", hash = "sha256:07a08422fc1936a5d26d146759e9155ea466e842f5ab2f7d2266dd084c8dab1f", size = 57205, upload-time = "2025-08-11T07:25:47.597Z" }, + { url = "https://files.pythonhosted.org/packages/22/5b/dbc6a8cddc9cfa9c4971d59fb12bb8d42e161b7e7f8cc89e49137c5b279c/mkdocs-1.6.1-py3-none-any.whl", hash = "sha256:db91759624d1647f3f34aa0c3f327dd2601beae39a366d6e064c03468d35c20e", size = 3864451, upload-time = "2024-08-30T12:24:05.054Z" }, ] [[package]] -name = "mdurl" -version = "0.1.2" +name = "mkdocs-autorefs" +version = "1.4.3" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/d6/54/cfe61301667036ec958cb99bd3efefba235e65cdeb9c84d24a8293ba1d90/mdurl-0.1.2.tar.gz", hash = "sha256:bb413d29f5eea38f31dd4754dd7377d4465116fb207585f97bf925588687c1ba", size = 8729, upload-time = "2022-08-14T12:40:10.846Z" } +dependencies = [ + { name = "markdown" }, + { name = "markupsafe" }, + { name = "mkdocs" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/51/fa/9124cd63d822e2bcbea1450ae68cdc3faf3655c69b455f3a7ed36ce6c628/mkdocs_autorefs-1.4.3.tar.gz", hash = "sha256:beee715b254455c4aa93b6ef3c67579c399ca092259cc41b7d9342573ff1fc75", size = 55425, upload-time = "2025-08-26T14:23:17.223Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/b3/38/89ba8ad64ae25be8de66a6d463314cf1eb366222074cfda9ee839c56a4b4/mdurl-0.1.2-py3-none-any.whl", hash = "sha256:84008a41e51615a49fc9966191ff91509e3c40b939176e643fd50a5c2196b8f8", size = 9979, upload-time = "2022-08-14T12:40:09.779Z" }, + { url = "https://files.pythonhosted.org/packages/9f/4d/7123b6fa2278000688ebd338e2a06d16870aaf9eceae6ba047ea05f92df1/mkdocs_autorefs-1.4.3-py3-none-any.whl", hash = "sha256:469d85eb3114801d08e9cc55d102b3ba65917a869b893403b8987b601cf55dc9", size = 25034, upload-time = "2025-08-26T14:23:15.906Z" }, ] [[package]] -name = "mistune" -version = "3.1.4" +name = "mkdocs-get-deps" +version = "0.2.0" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/d7/02/a7fb8b21d4d55ac93cdcde9d3638da5dd0ebdd3a4fed76c7725e10b81cbe/mistune-3.1.4.tar.gz", hash = "sha256:b5a7f801d389f724ec702840c11d8fc48f2b33519102fc7ee739e8177b672164", size = 94588, upload-time = "2025-08-29T07:20:43.594Z" } +dependencies = [ + { name = "mergedeep" }, + { name = "platformdirs" }, + { name = "pyyaml" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/98/f5/ed29cd50067784976f25ed0ed6fcd3c2ce9eb90650aa3b2796ddf7b6870b/mkdocs_get_deps-0.2.0.tar.gz", hash = "sha256:162b3d129c7fad9b19abfdcb9c1458a651628e4b1dea628ac68790fb3061c60c", size = 10239, upload-time = "2023-11-20T17:51:09.981Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/7a/f0/8282d9641415e9e33df173516226b404d367a0fc55e1a60424a152913abc/mistune-3.1.4-py3-none-any.whl", hash = "sha256:93691da911e5d9d2e23bc54472892aff676df27a75274962ff9edc210364266d", size = 53481, upload-time = "2025-08-29T07:20:42.218Z" }, + { url = "https://files.pythonhosted.org/packages/9f/d4/029f984e8d3f3b6b726bd33cafc473b75e9e44c0f7e80a5b29abc466bdea/mkdocs_get_deps-0.2.0-py3-none-any.whl", hash = "sha256:2bf11d0b133e77a0dd036abeeb06dec8775e46efa526dc70667d8863eefc6134", size = 9521, upload-time = "2023-11-20T17:51:08.587Z" }, +] + +[[package]] +name = "mkdocs-material" +version = "9.6.19" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "babel" }, + { name = "backrefs" }, + { name = "click" }, + { name = "colorama" }, + { name = "jinja2" }, + { name = "markdown" }, + { name = "mkdocs" }, + { name = "mkdocs-material-extensions" }, + { name = "paginate" }, + { name = "pygments" }, + { name = "pymdown-extensions" }, + { name = "requests" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/44/94/eb0fca39b19c2251b16bc759860a50f232655c4377116fa9c0e7db11b82c/mkdocs_material-9.6.19.tar.gz", hash = "sha256:80e7b3f9acabfee9b1f68bd12c26e59c865b3d5bbfb505fd1344e970db02c4aa", size = 4038202, upload-time = "2025-09-07T17:46:40.468Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/02/23/a2551d1038bedc2771366f65ff3680bb3a89674cd7ca6140850c859f1f71/mkdocs_material-9.6.19-py3-none-any.whl", hash = "sha256:7492d2ac81952a467ca8a10cac915d6ea5c22876932f44b5a0f4f8e7d68ac06f", size = 9240205, upload-time = "2025-09-07T17:46:36.484Z" }, +] + +[[package]] +name = "mkdocs-material-extensions" +version = "1.3.1" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/79/9b/9b4c96d6593b2a541e1cb8b34899a6d021d208bb357042823d4d2cabdbe7/mkdocs_material_extensions-1.3.1.tar.gz", hash = "sha256:10c9511cea88f568257f960358a467d12b970e1f7b2c0e5fb2bb48cab1928443", size = 11847, upload-time = "2023-11-22T19:09:45.208Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/5b/54/662a4743aa81d9582ee9339d4ffa3c8fd40a4965e033d77b9da9774d3960/mkdocs_material_extensions-1.3.1-py3-none-any.whl", hash = "sha256:adff8b62700b25cb77b53358dad940f3ef973dd6db797907c49e3c2ef3ab4e31", size = 8728, upload-time = "2023-11-22T19:09:43.465Z" }, +] + +[[package]] +name = "mkdocstrings" +version = "0.30.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "jinja2" }, + { name = "markdown" }, + { name = "markupsafe" }, + { name = "mkdocs" }, + { name = "mkdocs-autorefs" }, + { name = "pymdown-extensions" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/e2/0a/7e4776217d4802009c8238c75c5345e23014a4706a8414a62c0498858183/mkdocstrings-0.30.0.tar.gz", hash = "sha256:5d8019b9c31ddacd780b6784ffcdd6f21c408f34c0bd1103b5351d609d5b4444", size = 106597, upload-time = "2025-07-22T23:48:45.998Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/de/b4/3c5eac68f31e124a55d255d318c7445840fa1be55e013f507556d6481913/mkdocstrings-0.30.0-py3-none-any.whl", hash = "sha256:ae9e4a0d8c1789697ac776f2e034e2ddd71054ae1cf2c2bb1433ccfd07c226f2", size = 36579, upload-time = "2025-07-22T23:48:44.152Z" }, +] + +[package.optional-dependencies] +python = [ + { name = "mkdocstrings-python" }, +] + +[[package]] +name = "mkdocstrings-python" +version = "1.18.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "griffe" }, + { name = "mkdocs-autorefs" }, + { name = "mkdocstrings" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/95/ae/58ab2bfbee2792e92a98b97e872f7c003deb903071f75d8d83aa55db28fa/mkdocstrings_python-1.18.2.tar.gz", hash = "sha256:4ad536920a07b6336f50d4c6d5603316fafb1172c5c882370cbbc954770ad323", size = 207972, upload-time = "2025-08-28T16:11:19.847Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/d5/8f/ce008599d9adebf33ed144e7736914385e8537f5fc686fdb7cceb8c22431/mkdocstrings_python-1.18.2-py3-none-any.whl", hash = "sha256:944fe6deb8f08f33fa936d538233c4036e9f53e840994f6146e8e94eb71b600d", size = 138215, upload-time = "2025-08-28T16:11:18.176Z" }, ] [[package]] @@ -1285,23 +1416,6 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/51/c0/00c9809d8b9346eb238a6bbd5f83e846a4ce4503da94a4c08cb7284c325b/multipledispatch-1.0.0-py3-none-any.whl", hash = "sha256:0c53cd8b077546da4e48869f49b13164bebafd0c2a5afceb6bb6a316e7fb46e4", size = 12818, upload-time = "2023-06-27T16:45:09.418Z" }, ] -[[package]] -name = "myst-parser" -version = "4.0.1" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "docutils" }, - { name = "jinja2" }, - { name = "markdown-it-py" }, - { name = "mdit-py-plugins" }, - { name = "pyyaml" }, - { name = "sphinx" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/66/a5/9626ba4f73555b3735ad86247a8077d4603aa8628537687c839ab08bfe44/myst_parser-4.0.1.tar.gz", hash = "sha256:5cfea715e4f3574138aecbf7d54132296bfd72bb614d31168f48c477a830a7c4", size = 93985, upload-time = "2025-02-12T10:53:03.833Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/5f/df/76d0321c3797b54b60fef9ec3bd6f4cfd124b9e422182156a1dd418722cf/myst_parser-4.0.1-py3-none-any.whl", hash = "sha256:9134e88959ec3b5780aedf8a99680ea242869d012e8821db3126d427edc9c95d", size = 84579, upload-time = "2025-02-12T10:53:02.078Z" }, -] - [[package]] name = "nbclient" version = "0.10.2" @@ -1492,6 +1606,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/20/12/38679034af332785aac8774540895e234f4d07f7545804097de4b666afd8/packaging-25.0-py3-none-any.whl", hash = "sha256:29572ef2b1f17581046b3a2227d5c611fb25ec70ca1ba8554b24b0e69331a484", size = 66469, upload-time = "2025-04-19T11:48:57.875Z" }, ] +[[package]] +name = "paginate" +version = "0.5.7" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/ec/46/68dde5b6bc00c1296ec6466ab27dddede6aec9af1b99090e1107091b3b84/paginate-0.5.7.tar.gz", hash = "sha256:22bd083ab41e1a8b4f3690544afb2c60c25e5c9a63a30fa2f483f6c60c8e5945", size = 19252, upload-time = "2024-08-25T14:17:24.139Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/90/96/04b8e52da071d28f5e21a805b19cb9390aa17a47462ac87f5e2696b9566d/paginate-0.5.7-py2.py3-none-any.whl", hash = "sha256:b885e2af73abcf01d9559fd5216b57ef722f8c42affbb63942377668e35c7591", size = 13746, upload-time = "2024-08-25T14:17:22.55Z" }, +] + [[package]] name = "pandas" version = "2.3.2" @@ -1537,6 +1660,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/16/32/f8e3c85d1d5250232a5d3477a2a28cc291968ff175caeadaf3cc19ce0e4a/parso-0.8.5-py2.py3-none-any.whl", hash = "sha256:646204b5ee239c396d040b90f9e272e9a8017c630092bf59980beb62fd033887", size = 106668, upload-time = "2025-08-23T15:15:25.663Z" }, ] +[[package]] +name = "pathspec" +version = "0.12.1" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/ca/bc/f35b8446f4531a7cb215605d100cd88b7ac6f44ab3fc94870c120ab3adbf/pathspec-0.12.1.tar.gz", hash = "sha256:a482d51503a1ab33b1c67a6c3813a26953dbdc71c31dacaef9a838c4e29f5712", size = 51043, upload-time = "2023-12-10T22:30:45Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/cc/20/ff623b09d963f88bfde16306a54e12ee5ea43e9b597108672ff3a408aad6/pathspec-0.12.1-py3-none-any.whl", hash = "sha256:a0d503e138a4c123b27490a4f7beda6a01c6f288df0e4a8b79c7eb0dc7b4cc08", size = 31191, upload-time = "2023-12-10T22:30:43.14Z" }, +] + [[package]] name = "pexpect" version = "4.9.0" @@ -1700,31 +1832,25 @@ wheels = [ ] [[package]] -name = "pydata-sphinx-theme" -version = "0.15.4" +name = "pygments" +version = "2.19.2" source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "accessible-pygments" }, - { name = "babel" }, - { name = "beautifulsoup4" }, - { name = "docutils" }, - { name = "packaging" }, - { name = "pygments" }, - { name = "sphinx" }, - { name = "typing-extensions" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/67/ea/3ab478cccacc2e8ef69892c42c44ae547bae089f356c4b47caf61730958d/pydata_sphinx_theme-0.15.4.tar.gz", hash = "sha256:7762ec0ac59df3acecf49fd2f889e1b4565dbce8b88b2e29ee06fdd90645a06d", size = 2400673, upload-time = "2024-06-25T19:28:45.041Z" } +sdist = { url = "https://files.pythonhosted.org/packages/b0/77/a5b8c569bf593b0140bde72ea885a803b82086995367bf2037de0159d924/pygments-2.19.2.tar.gz", hash = "sha256:636cb2477cec7f8952536970bc533bc43743542f70392ae026374600add5b887", size = 4968631, upload-time = "2025-06-21T13:39:12.283Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/e7/d3/c622950d87a2ffd1654208733b5bd1c5645930014abed8f4c0d74863988b/pydata_sphinx_theme-0.15.4-py3-none-any.whl", hash = "sha256:2136ad0e9500d0949f96167e63f3e298620040aea8f9c74621959eda5d4cf8e6", size = 4640157, upload-time = "2024-06-25T19:28:42.383Z" }, + { url = "https://files.pythonhosted.org/packages/c7/21/705964c7812476f378728bdf590ca4b771ec72385c533964653c68e86bdc/pygments-2.19.2-py3-none-any.whl", hash = "sha256:86540386c03d588bb81d44bc3928634ff26449851e99741617ecb9037ee5ec0b", size = 1225217, upload-time = "2025-06-21T13:39:07.939Z" }, ] [[package]] -name = "pygments" -version = "2.19.2" +name = "pymdown-extensions" +version = "10.16.1" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/b0/77/a5b8c569bf593b0140bde72ea885a803b82086995367bf2037de0159d924/pygments-2.19.2.tar.gz", hash = "sha256:636cb2477cec7f8952536970bc533bc43743542f70392ae026374600add5b887", size = 4968631, upload-time = "2025-06-21T13:39:12.283Z" } +dependencies = [ + { name = "markdown" }, + { name = "pyyaml" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/55/b3/6d2b3f149bc5413b0a29761c2c5832d8ce904a1d7f621e86616d96f505cc/pymdown_extensions-10.16.1.tar.gz", hash = "sha256:aace82bcccba3efc03e25d584e6a22d27a8e17caa3f4dd9f207e49b787aa9a91", size = 853277, upload-time = "2025-07-28T16:19:34.167Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/c7/21/705964c7812476f378728bdf590ca4b771ec72385c533964653c68e86bdc/pygments-2.19.2-py3-none-any.whl", hash = "sha256:86540386c03d588bb81d44bc3928634ff26449851e99741617ecb9037ee5ec0b", size = 1225217, upload-time = "2025-06-21T13:39:07.939Z" }, + { url = "https://files.pythonhosted.org/packages/e4/06/43084e6cbd4b3bc0e80f6be743b2e79fbc6eed8de9ad8c629939fa55d972/pymdown_extensions-10.16.1-py3-none-any.whl", hash = "sha256:d6ba157a6c03146a7fb122b2b9a121300056384eafeec9c9f9e584adfdb2a32d", size = 266178, upload-time = "2025-07-28T16:19:31.401Z" }, ] [[package]] @@ -1759,16 +1885,13 @@ docs = [ { name = "ipywidgets" }, { name = "jupyter" }, { name = "matplotlib" }, - { name = "myst-parser" }, + { name = "mkdocs" }, + { name = "mkdocs-material" }, + { name = "mkdocstrings", extra = ["python"] }, { name = "nbclient" }, { name = "nbformat" }, { name = "numpydoc" }, { name = "pyyaml" }, - { name = "sphinx" }, - { name = "sphinx-autodoc-typehints" }, - { name = "sphinx-book-theme" }, - { name = "sphinxcontrib-katex" }, - { name = "sphinxcontrib-mermaid" }, ] test = [ { name = "pytest" }, @@ -1787,7 +1910,9 @@ requires-dist = [ { name = "jax", specifier = ">=0.7.0" }, { name = "jupyter", marker = "extra == 'docs'", specifier = ">=1.0.0" }, { name = "matplotlib", marker = "extra == 'docs'", specifier = ">=3.8.3" }, - { name = "myst-parser", marker = "extra == 'docs'", specifier = ">=4.0.0" }, + { name = "mkdocs", marker = "extra == 'docs'", specifier = ">=1.5.0" }, + { name = "mkdocs-material", marker = "extra == 'docs'", specifier = ">=9.0.0" }, + { name = "mkdocstrings", extras = ["python"], marker = "extra == 'docs'", specifier = ">=0.24.0" }, { name = "nbclient", marker = "extra == 'docs'", specifier = ">=0.10.0" }, { name = "nbconvert", marker = "extra == 'dev'", specifier = ">=7.16.4" }, { name = "nbformat", marker = "extra == 'docs'", specifier = ">=5.10.0" }, @@ -1800,11 +1925,6 @@ requires-dist = [ { name = "pytest-mpl", marker = "extra == 'test'", specifier = ">=0.17.0" }, { name = "pyyaml", marker = "extra == 'docs'", specifier = ">=6.0.0" }, { name = "scipy", marker = "extra == 'test'", specifier = ">=1.14.1" }, - { name = "sphinx", marker = "extra == 'docs'", specifier = ">=7.2.6" }, - { name = "sphinx-autodoc-typehints", marker = "extra == 'docs'", specifier = ">=2.1.0" }, - { name = "sphinx-book-theme", marker = "extra == 'docs'", specifier = ">=1.1.2" }, - { name = "sphinxcontrib-katex", marker = "extra == 'docs'", specifier = ">=0.9.10" }, - { name = "sphinxcontrib-mermaid", marker = "extra == 'docs'", specifier = ">=0.9.2" }, ] provides-extras = ["dev", "docs", "test"] @@ -1924,6 +2044,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/fa/de/02b54f42487e3d3c6efb3f89428677074ca7bf43aae402517bc7cca949f3/PyYAML-6.0.2-cp313-cp313-win_amd64.whl", hash = "sha256:8388ee1976c416731879ac16da0aff3f63b286ffdd57cdeb95f3f2e085687563", size = 156446, upload-time = "2024-08-06T20:33:04.33Z" }, ] +[[package]] +name = "pyyaml-env-tag" +version = "1.1" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "pyyaml" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/eb/2e/79c822141bfd05a853236b504869ebc6b70159afc570e1d5a20641782eaa/pyyaml_env_tag-1.1.tar.gz", hash = "sha256:2eb38b75a2d21ee0475d6d97ec19c63287a7e140231e4214969d0eac923cd7ff", size = 5737, upload-time = "2025-05-13T15:24:01.64Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/04/11/432f32f8097b03e3cd5fe57e88efb685d964e2e5178a48ed61e841f7fdce/pyyaml_env_tag-1.1-py3-none-any.whl", hash = "sha256:17109e1a528561e32f026364712fee1264bc2ea6715120891174ed1b980d2e04", size = 4722, upload-time = "2025-05-13T15:23:59.629Z" }, +] + [[package]] name = "pyzmq" version = "27.0.2" @@ -2244,31 +2376,6 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/31/53/136e9eca6e0b9dc0e1962e2c908fbea2e5ac000c2a2fbd9a35797958c48b/sphinx-8.2.3-py3-none-any.whl", hash = "sha256:4405915165f13521d875a8c29c8970800a0141c14cc5416a38feca4ea5d9b9c3", size = 3589741, upload-time = "2025-03-02T22:31:56.836Z" }, ] -[[package]] -name = "sphinx-autodoc-typehints" -version = "3.2.0" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "sphinx" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/93/68/a388a9b8f066cd865d9daa65af589d097efbfab9a8c302d2cb2daa43b52e/sphinx_autodoc_typehints-3.2.0.tar.gz", hash = "sha256:107ac98bc8b4837202c88c0736d59d6da44076e65a0d7d7d543a78631f662a9b", size = 36724, upload-time = "2025-04-25T16:53:25.872Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/f7/c7/8aab362e86cbf887e58be749a78d20ad743e1eb2c73c2b13d4761f39a104/sphinx_autodoc_typehints-3.2.0-py3-none-any.whl", hash = "sha256:884b39be23b1d884dcc825d4680c9c6357a476936e3b381a67ae80091984eb49", size = 20563, upload-time = "2025-04-25T16:53:24.492Z" }, -] - -[[package]] -name = "sphinx-book-theme" -version = "1.1.4" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "pydata-sphinx-theme" }, - { name = "sphinx" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/45/19/d002ed96bdc7738c15847c730e1e88282d738263deac705d5713b4d8fa94/sphinx_book_theme-1.1.4.tar.gz", hash = "sha256:73efe28af871d0a89bd05856d300e61edce0d5b2fbb7984e84454be0fedfe9ed", size = 439188, upload-time = "2025-02-20T16:32:32.581Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/51/9e/c41d68be04eef5b6202b468e0f90faf0c469f3a03353f2a218fd78279710/sphinx_book_theme-1.1.4-py3-none-any.whl", hash = "sha256:843b3f5c8684640f4a2d01abd298beb66452d1b2394cd9ef5be5ebd5640ea0e1", size = 433952, upload-time = "2025-02-20T16:32:31.009Z" }, -] - [[package]] name = "sphinxcontrib-applehelp" version = "2.0.0" @@ -2305,31 +2412,6 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/c2/42/4c8646762ee83602e3fb3fbe774c2fac12f317deb0b5dbeeedd2d3ba4b77/sphinxcontrib_jsmath-1.0.1-py2.py3-none-any.whl", hash = "sha256:2ec2eaebfb78f3f2078e73666b1415417a116cc848b72e5172e596c871103178", size = 5071, upload-time = "2019-01-21T16:10:14.333Z" }, ] -[[package]] -name = "sphinxcontrib-katex" -version = "0.9.11" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "sphinx" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/9a/cd/71d3c7fc1a4e40a03048ce68787e4e4b6f612dc6b95acc89f5010d4cda7a/sphinxcontrib_katex-0.9.11.tar.gz", hash = "sha256:2d32b210d20bbee451a51d19645f6fef5f5568b9a54e2813affb96efa661238a", size = 100476, upload-time = "2025-06-27T12:04:12.703Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/33/14/c62d4553dbf27156b6549ceed1a9768f801c979237a2cf0c3d074e16272f/sphinxcontrib_katex-0.9.11-py3-none-any.whl", hash = "sha256:b0a318d86b8f43314ccdd481ce0db90283ac51f51afd9479969f36a0c291f44e", size = 97956, upload-time = "2025-06-27T12:04:11.478Z" }, -] - -[[package]] -name = "sphinxcontrib-mermaid" -version = "1.0.0" -source = { registry = "https://pypi.org/simple" } -dependencies = [ - { name = "pyyaml" }, - { name = "sphinx" }, -] -sdist = { url = "https://files.pythonhosted.org/packages/97/69/bf039237ad260073e8c02f820b3e00dc34f3a2de20aff7861e6b19d2f8c5/sphinxcontrib_mermaid-1.0.0.tar.gz", hash = "sha256:2e8ab67d3e1e2816663f9347d026a8dee4a858acdd4ad32dd1c808893db88146", size = 15153, upload-time = "2024-10-12T16:33:03.863Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/cd/c8/784b9ac6ea08aa594c1a4becbd0dbe77186785362e31fd633b8c6ae0197a/sphinxcontrib_mermaid-1.0.0-py3-none-any.whl", hash = "sha256:60b72710ea02087f212028feb09711225fbc2e343a10d34822fe787510e1caa3", size = 9597, upload-time = "2024-10-12T16:33:02.303Z" }, -] - [[package]] name = "sphinxcontrib-qthelp" version = "2.0.0" @@ -2473,6 +2555,27 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl", hash = "sha256:e6b01673c0fa6a13e374b50871808eb3bf7046c4b125b216f6bf1cc604cff0dc", size = 129795, upload-time = "2025-06-18T14:07:40.39Z" }, ] +[[package]] +name = "watchdog" +version = "6.0.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/db/7d/7f3d619e951c88ed75c6037b246ddcf2d322812ee8ea189be89511721d54/watchdog-6.0.0.tar.gz", hash = "sha256:9ddf7c82fda3ae8e24decda1338ede66e1c99883db93711d8fb941eaa2d8c282", size = 131220, upload-time = "2024-11-01T14:07:13.037Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/68/98/b0345cabdce2041a01293ba483333582891a3bd5769b08eceb0d406056ef/watchdog-6.0.0-cp313-cp313-macosx_10_13_universal2.whl", hash = "sha256:490ab2ef84f11129844c23fb14ecf30ef3d8a6abafd3754a6f75ca1e6654136c", size = 96480, upload-time = "2024-11-01T14:06:42.952Z" }, + { url = "https://files.pythonhosted.org/packages/85/83/cdf13902c626b28eedef7ec4f10745c52aad8a8fe7eb04ed7b1f111ca20e/watchdog-6.0.0-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:76aae96b00ae814b181bb25b1b98076d5fc84e8a53cd8885a318b42b6d3a5134", size = 88451, upload-time = "2024-11-01T14:06:45.084Z" }, + { url = "https://files.pythonhosted.org/packages/fe/c4/225c87bae08c8b9ec99030cd48ae9c4eca050a59bf5c2255853e18c87b50/watchdog-6.0.0-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:a175f755fc2279e0b7312c0035d52e27211a5bc39719dd529625b1930917345b", size = 89057, upload-time = "2024-11-01T14:06:47.324Z" }, + { url = "https://files.pythonhosted.org/packages/a9/c7/ca4bf3e518cb57a686b2feb4f55a1892fd9a3dd13f470fca14e00f80ea36/watchdog-6.0.0-py3-none-manylinux2014_aarch64.whl", hash = "sha256:7607498efa04a3542ae3e05e64da8202e58159aa1fa4acddf7678d34a35d4f13", size = 79079, upload-time = "2024-11-01T14:06:59.472Z" }, + { url = "https://files.pythonhosted.org/packages/5c/51/d46dc9332f9a647593c947b4b88e2381c8dfc0942d15b8edc0310fa4abb1/watchdog-6.0.0-py3-none-manylinux2014_armv7l.whl", hash = "sha256:9041567ee8953024c83343288ccc458fd0a2d811d6a0fd68c4c22609e3490379", size = 79078, upload-time = "2024-11-01T14:07:01.431Z" }, + { url = "https://files.pythonhosted.org/packages/d4/57/04edbf5e169cd318d5f07b4766fee38e825d64b6913ca157ca32d1a42267/watchdog-6.0.0-py3-none-manylinux2014_i686.whl", hash = "sha256:82dc3e3143c7e38ec49d61af98d6558288c415eac98486a5c581726e0737c00e", size = 79076, upload-time = "2024-11-01T14:07:02.568Z" }, + { url = "https://files.pythonhosted.org/packages/ab/cc/da8422b300e13cb187d2203f20b9253e91058aaf7db65b74142013478e66/watchdog-6.0.0-py3-none-manylinux2014_ppc64.whl", hash = "sha256:212ac9b8bf1161dc91bd09c048048a95ca3a4c4f5e5d4a7d1b1a7d5752a7f96f", size = 79077, upload-time = "2024-11-01T14:07:03.893Z" }, + { url = "https://files.pythonhosted.org/packages/2c/3b/b8964e04ae1a025c44ba8e4291f86e97fac443bca31de8bd98d3263d2fcf/watchdog-6.0.0-py3-none-manylinux2014_ppc64le.whl", hash = "sha256:e3df4cbb9a450c6d49318f6d14f4bbc80d763fa587ba46ec86f99f9e6876bb26", size = 79078, upload-time = "2024-11-01T14:07:05.189Z" }, + { url = "https://files.pythonhosted.org/packages/62/ae/a696eb424bedff7407801c257d4b1afda455fe40821a2be430e173660e81/watchdog-6.0.0-py3-none-manylinux2014_s390x.whl", hash = "sha256:2cce7cfc2008eb51feb6aab51251fd79b85d9894e98ba847408f662b3395ca3c", size = 79077, upload-time = "2024-11-01T14:07:06.376Z" }, + { url = "https://files.pythonhosted.org/packages/b5/e8/dbf020b4d98251a9860752a094d09a65e1b436ad181faf929983f697048f/watchdog-6.0.0-py3-none-manylinux2014_x86_64.whl", hash = "sha256:20ffe5b202af80ab4266dcd3e91aae72bf2da48c0d33bdb15c66658e685e94e2", size = 79078, upload-time = "2024-11-01T14:07:07.547Z" }, + { url = "https://files.pythonhosted.org/packages/07/f6/d0e5b343768e8bcb4cda79f0f2f55051bf26177ecd5651f84c07567461cf/watchdog-6.0.0-py3-none-win32.whl", hash = "sha256:07df1fdd701c5d4c8e55ef6cf55b8f0120fe1aef7ef39a1c6fc6bc2e606d517a", size = 79065, upload-time = "2024-11-01T14:07:09.525Z" }, + { url = "https://files.pythonhosted.org/packages/db/d9/c495884c6e548fce18a8f40568ff120bc3a4b7b99813081c8ac0c936fa64/watchdog-6.0.0-py3-none-win_amd64.whl", hash = "sha256:cbafb470cf848d93b5d013e2ecb245d4aa1c8fd0504e863ccefa32445359d680", size = 79070, upload-time = "2024-11-01T14:07:10.686Z" }, + { url = "https://files.pythonhosted.org/packages/33/e8/e40370e6d74ddba47f002a32919d91310d6074130fe4e17dabcafc15cbf1/watchdog-6.0.0-py3-none-win_ia64.whl", hash = "sha256:a1914259fa9e1454315171103c6a30961236f508b9b623eae470268bbcc6a22f", size = 79067, upload-time = "2024-11-01T14:07:11.845Z" }, +] + [[package]] name = "wcwidth" version = "0.2.13" From c705a5495900355c6ea907b73463c906ac78607e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 9 Sep 2025 19:52:18 +0000 Subject: [PATCH 3/8] Complete cleanup of Sphinx files and finalize MkDocs migration Co-authored-by: damonbayer <8095080+damonbayer@users.noreply.github.com> --- docs/.gitignore | 1 + docs/README.md | 3 - docs/make.bat | 35 -- docs/source/_static/pyrenew.css | 11 - docs/source/conf.py | 118 ---- docs/source/developer_documentation.md | 110 ---- docs/source/general/ctoc.md | 15 - docs/source/index.md | 5 - docs/source/pyrenew_reference/additional.md | 60 -- docs/source/pyrenew_reference/datasets.md | 28 - .../source/pyrenew_reference/deterministic.md | 8 - .../source/pyrenew_reference/distributions.md | 8 - docs/source/pyrenew_reference/index.md | 23 - docs/source/pyrenew_reference/latent.md | 50 -- docs/source/pyrenew_reference/model.md | 23 - docs/source/pyrenew_reference/observation.md | 23 - docs/source/pyrenew_reference/process.md | 8 - .../pyrenew_reference/randomvariable.md | 8 - docs/source/tutorials/.gitignore | 10 - docs/source/tutorials/basic_renewal_model.md | 5 - docs/source/tutorials/basic_renewal_model.qmd | 407 ------------ docs/source/tutorials/day_of_the_week.qmd | 426 ------------- docs/source/tutorials/extending_pyrenew.qmd | 297 --------- docs/source/tutorials/getting_started.qmd | 60 -- .../tutorials/hospital_admissions_model.md | 5 - .../tutorials/hospital_admissions_model.qmd | 591 ------------------ docs/source/tutorials/index.md | 14 - docs/source/tutorials/periodic_effects.qmd | 105 ---- docs/tutorials/basic_renewal_model.md | 5 + .../{source => }/tutorials/day_of_the_week.md | 2 +- .../tutorials/extending_pyrenew.md | 2 +- .../{source => }/tutorials/getting_started.md | 2 +- docs/tutorials/hospital_admissions_model.md | 5 + .../tutorials/periodic_effects.md | 2 +- 34 files changed, 15 insertions(+), 2460 deletions(-) delete mode 100644 docs/README.md delete mode 100644 docs/make.bat delete mode 100644 docs/source/_static/pyrenew.css delete mode 100644 docs/source/conf.py delete mode 100644 docs/source/developer_documentation.md delete mode 100644 docs/source/general/ctoc.md delete mode 100644 docs/source/index.md delete mode 100644 docs/source/pyrenew_reference/additional.md delete mode 100644 docs/source/pyrenew_reference/datasets.md delete mode 100644 docs/source/pyrenew_reference/deterministic.md delete mode 100644 docs/source/pyrenew_reference/distributions.md delete mode 100644 docs/source/pyrenew_reference/index.md delete mode 100644 docs/source/pyrenew_reference/latent.md delete mode 100644 docs/source/pyrenew_reference/model.md delete mode 100644 docs/source/pyrenew_reference/observation.md delete mode 100644 docs/source/pyrenew_reference/process.md delete mode 100644 docs/source/pyrenew_reference/randomvariable.md delete mode 100644 docs/source/tutorials/.gitignore delete mode 100644 docs/source/tutorials/basic_renewal_model.md delete mode 100644 docs/source/tutorials/basic_renewal_model.qmd delete mode 100644 docs/source/tutorials/day_of_the_week.qmd delete mode 100644 docs/source/tutorials/extending_pyrenew.qmd delete mode 100644 docs/source/tutorials/getting_started.qmd delete mode 100644 docs/source/tutorials/hospital_admissions_model.md delete mode 100644 docs/source/tutorials/hospital_admissions_model.qmd delete mode 100644 docs/source/tutorials/index.md delete mode 100644 docs/source/tutorials/periodic_effects.qmd create mode 100644 docs/tutorials/basic_renewal_model.md rename docs/{source => }/tutorials/day_of_the_week.md (58%) rename docs/{source => }/tutorials/extending_pyrenew.md (57%) rename docs/{source => }/tutorials/getting_started.md (58%) create mode 100644 docs/tutorials/hospital_admissions_model.md rename docs/{source => }/tutorials/periodic_effects.md (57%) diff --git a/docs/.gitignore b/docs/.gitignore index 567609b1..a303fff2 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -1 +1,2 @@ build/ +site/ diff --git a/docs/README.md b/docs/README.md deleted file mode 100644 index 193f2ce0..00000000 --- a/docs/README.md +++ /dev/null @@ -1,3 +0,0 @@ -# Sphinx website source code - -The files under this folder contain the source code for the project's Sphinx website. The website is automatically built and published by GitHub Actions. The website is hosted at . diff --git a/docs/make.bat b/docs/make.bat deleted file mode 100644 index 747ffb7b..00000000 --- a/docs/make.bat +++ /dev/null @@ -1,35 +0,0 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=sphinx-build -) -set SOURCEDIR=source -set BUILDDIR=build - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The 'sphinx-build' command was not found. Make sure you have Sphinx - echo.installed, then set the SPHINXBUILD environment variable to point - echo.to the full path of the 'sphinx-build' executable. Alternatively you - echo.may add the Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.https://www.sphinx-doc.org/ - exit /b 1 -) - -if "%1" == "" goto help - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% - -:end -popd diff --git a/docs/source/_static/pyrenew.css b/docs/source/_static/pyrenew.css deleted file mode 100644 index eb7823d8..00000000 --- a/docs/source/_static/pyrenew.css +++ /dev/null @@ -1,11 +0,0 @@ -/* @import url('https://your-font-hosting-url.com/path/to/iosevka.css'); - -body { - font-family: 'Iosevka', monospace; -} */ - -/* center equations by default - when using KaTeX rendering, - see https://github.com/KaTeX/KaTeX/discussions/3923 -*/ -div.math { justify-content: center } diff --git a/docs/source/conf.py b/docs/source/conf.py deleted file mode 100644 index f5a679ae..00000000 --- a/docs/source/conf.py +++ /dev/null @@ -1,118 +0,0 @@ -# numpydoc ignore=GL08 - -import os -import sys - -# Configuration file for the Sphinx documentation builder. -# -# For the full list of built-in configuration values, see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html - -# -- Project information ----------------------------------------------------- -# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information - -project = "Multisignal Renewal Models for Epi Inference" -copyright = "2024, CDC CFA" -author = "CDC Center for Forecasting Analytics" -release = "0.1.0" - -# -- General configuration --------------------------------------------------- - -sys.path.insert(0, os.path.abspath("../../pyrenew")) - -# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration - -extensions = [ - "sphinx.ext.autodoc", - "sphinx.ext.autosectionlabel", - "sphinx.ext.autosummary", - "sphinx.ext.intersphinx", - "sphinx.ext.doctest", - "sphinx.ext.napoleon", # numpydoc - "sphinx.ext.duration", - "sphinx.ext.githubpages", - "sphinx.ext.inheritance_diagram", - "sphinx.ext.todo", - "sphinx.ext.viewcode", - "sphinxcontrib.katex", - "sphinxcontrib.mermaid", - "sphinx_autodoc_typehints", - "myst_parser", -] - -# Mock heavy / optional dependencies so docs can build in minimal envs (CI, RTD) -autodoc_mock_imports = [ - "jax", - "jaxlib", - "numpyro", - "polars", - "arviz", -] - -# Simplifies printing of type hints -set_type_checking_flag = True -typehints_fully_qualified = False - -# Avoid appending the full module name to the class name -add_module_names = False - -templates_path = ["_templates"] -exclude_patterns = [] - -# Default deph for documentation -toc_deph = 2 - -# We don't want that explicitly -todo_include_todos = False - -# Numpydocs checks -# numpydoc_validation_checks = {"SA01", "EX01"} - -# -- Options for HTML output ------------------------------------------------- -# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output - -html_theme = "sphinx_book_theme" -html_theme_options = { - "repository_url": "https://github.com/CDCgov/PyRenew", - "use_edit_page_button": True, - "use_issues_button": True, - "use_repository_button": True, - "repository_branch": "main", - "path_to_docs": "docs/source", - "use_download_button": True, -} - -html_static_path = ["_static"] -html_css_files = ["pyrenew.css"] - -html_sidebars = { - "**": [ - "navbar-logo.html", - "search-field.html", - "sbt-sidebar-nav.html", - ] -} - -master_doc = "general/ctoc" - -myst_fence_as_directive = ["mermaid"] -myst_enable_extensions = ["amsmath", "dollarmath"] - -intersphinx_mapping = { - "python": ("https://docs.python.org/3/", None), - "matplotlib": ("https://matplotlib.org/stable/", None), - "numpy": ("https://numpy.org/doc/stable/", None), - "numpyro": ("https://num.pyro.ai/en/latest/", None), - "jax": ("https://jax.readthedocs.io/en/latest/", None), - "polars": ("https://docs.pola.rs/api/python/stable/", None), -} - -napoleon_preprocess_types = True -autodoc_typehints = "description" -autodoc_typehints_format = "short" -autodoc_type_aliases = { - "ArrayLike": ":obj:`ArrayLike `", - "RandomVariable": ":class:`RandomVariable `", - "Any": ":obj:`Any `", -} -napoleon_type_aliases = autodoc_type_aliases diff --git a/docs/source/developer_documentation.md b/docs/source/developer_documentation.md deleted file mode 100644 index 403336b1..00000000 --- a/docs/source/developer_documentation.md +++ /dev/null @@ -1,110 +0,0 @@ -# Developer Documentation - -**Note: this document is a work in progress. Contrbitions to all sections are welcome.** - -## GitHub Workflow - -- You should request reviews for pull requests from `PyRenew-devs`, unless you have a good reason to request reviews from a smaller group. -- Reviews from all of `PyRenew-devs` are encouraged, but an approving review from a [codeowner](https://github.com/CDCgov/PyRenew/blob/main/.github/CODEOWNERS) ([@dylanhmorris](https://github.com/dylanhmorris) or [@damonbayer](https://github.com/damonbayer) is required before a pull request can be merged to `main`. -- For CDC contributors: if your pull request has not received a review at the time of the next standup, use standup to find a reviewer. -- External contributors should expect to receive a review within a few days of creating a pull request. -- If you create a draft pull request, indicate what, if anything, about the current pull request should be reviewed. -- Only mark a pull request as “ready for review” if you think it is ready to be merged. This indicates that a thorough, all-encompassing review should be given. - -## Installation for Developers - -- `uv sync --extra dev` -- `pre-commit install` - -## Coding Conventions - -A variety of coding conventions are enforced by automated tools in continuous integration ([black](https://github.com/psf/black), [isort](https://github.com/PyCQA/isort), [ruff](https://github.com/astral-sh/ruff), [numpydoc-validation](https://github.com/numpy/numpydoc)) via [pre-commit](https://github.com/pre-commit/pre-commit) hooks. - -## PyRenew Principles - -- Variable naming conventions - - - Use the `data_` prefix for (potentially) observed data. - - Use the `_rv` suffix for random variables. - - Use the `observed_` for the output of sample statements where `obs` is a `data_` prefixed object. - - Thus, code which may reasonably written like `infections = infections.sample(x, obs=infections)` should instead be written `observed_infections = infections_rv.sample(x, obs=data_infections)`. - -- Class conventions - - - Composing models is discouraged. - - Returning anything from `Model.sample` is discouraged. Instead, sample from models using `Predictive` or our `prior_predictive` or `posterior_predictive` functions. - - Using `numpyro.deterministic` within a `RandomVariable` is discouraged. Only use at the `numpyro.deterministic` `Model` level. If something might need to be recorded from a `RandomVariable`, it should be returned from the `RandomVariable` so it can be recorded at the `Model` level. - - Using default site names in a `RandomVariable` is discouraged. Only use default site names at the `Model` level. - - Use `DeterministicVariable`s instead of constants within a model. - -- `scan` conventions - - - Use `jax.lax.scan` for any scan whose iterations are deterministic, i.e. iterations contain no internal calls to `RandomVariable.sample()` or `numpyro.sample()`. - - Use `numpyro.scan` for any scan whose the iterations are stochastic, i.e. the iterations potentially include calls to `RandomVariable.sample()` or `numpyro.sample()`. - -- Multidimensional array conventions - - - In a multidimensional array of timeseries, time is always the first dimension. By default, `jax.lax.scan()` and `numpyro.contrib.control_flow.scan()` build output arrays by augmenting the first dimension, and variables are often scanned over time, making default output of scan over time sensible. - -## Adding Documentation to Sphinx - -The project uses Sphinx to generate documentation. To learn more about Sphinx, visit the [Sphinx documentation](https://www.sphinx-doc.org/en/master/). - -### New module or classes -We use Sphinx's `automodule` functionality to generate module API references. When you add a new module, make sure it gets added to the documentation. - -For example, under the `./docs/source/pyrenew_reference`, the `index.md` file lists the `HospitalAdmissions` module by adding the following entry: - -````markdown -## Hospital Admissions - -```{eval-rst} -.. automodule:: pyrenew.latent.hospitaladmissions - :members: - :undoc-members: - :show-inheritance: -``` -```` - -This entry tells Sphinx to generate documentation for the `HospitalAdmissions` module and its members. - -### New tutorials - -`PyRenew` tutorials are [quarto documents](https://quarto.org) located under [./docs/source/tutorials](https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials). Tutorials are automatically rendered using GitHub actions (see the Workflow file [here](https://github.com/CDCgov/PyRenew/actions/workflows/website.yaml)). - -To make the new tutorial available in the website, developers should follow these steps: - -1. Create a new `quarto` file in the `./docs/source/tutorials` directory. For instance, the `example_with_datasets.qmd` file was added to the repository. -2. Add an entry in the `./docs/source/tutorials/index.md`, for example: - -````markdown -```{toctree} -:maxdepth: 2 -getting_started -example_with_datasets -``` -```` - -3. Add an `md` entry with the same basename as the `quarto` file in the `./docs/source/tutorials` directory. For instance, the `example_with_datasets.md` file was added to the repository. This last step can be done running the bash script [./hook_scripts/pre-commit-md-placeholder.sh](https://github.com/CDCgov/PyRenew/blob/main/hook_scripts/pre-commit-md-placeholder.sh). Note the script should be executed by `pre-commit`. - -### Adding new pages - -Sphinx also allows adding arbitrary pages. For instance, all the `PyRenew` tutorials are additional documentation. The steps to add a new page are: - -1. Create a `md` file in the appropriate directory. For example, this file about development was added under `./docs/source/developer_documentation.md`. -2. Make sure the new `md` file is included in an indexed file, for instance, `./docs/source/general/ctoc.md`. Here is how it looks: - -````markdown -# Complete Table Of Contents - -```{toctree} -:maxdepth: 2 -../index -../pyrenew_reference/index -../tutorials/index -../genindex -../developer_documentation -``` -```` - -The last entry is the `developer_documentation` page. diff --git a/docs/source/general/ctoc.md b/docs/source/general/ctoc.md deleted file mode 100644 index 659c39ec..00000000 --- a/docs/source/general/ctoc.md +++ /dev/null @@ -1,15 +0,0 @@ -# Complete Table Of Contents - -```{toctree} -:maxdepth: 2 - -../index -../pyrenew_reference/index -../tutorials/index -../genindex -../developer_documentation -``` - -```{eval-rst} -.. todo:: Change order and naming of these toctree entities -``` diff --git a/docs/source/index.md b/docs/source/index.md deleted file mode 100644 index ffed5948..00000000 --- a/docs/source/index.md +++ /dev/null @@ -1,5 +0,0 @@ -## Placeholder file - -This `index.md` file is only a placeholder. Do not edit this file. The final `index.md` document is automatically generated via GitHub Actions using the root level `README.md` file: https://github.com/CDCgov/PyRenew/tree/main/README.md - -Any changes made to this file will be discarted during the next build. The GitHub action that builds the docs is: https://github.com/CDCgov/PyRenew/tree/main/.github/workflows/website.yaml diff --git a/docs/source/pyrenew_reference/additional.md b/docs/source/pyrenew_reference/additional.md deleted file mode 100644 index 7559a639..00000000 --- a/docs/source/pyrenew_reference/additional.md +++ /dev/null @@ -1,60 +0,0 @@ -# Additional modules - -## Metaclass Module - -```{eval-rst} -.. automodule:: pyrenew.metaclass - :members: - :undoc-members: - :show-inheritance: -``` - -## Convolution Utility Module - -```{eval-rst} -.. automodule:: pyrenew.convolve - :members: - :undoc-members: - :show-inheritance: -``` - -## Mathematics Utilities Module - -```{eval-rst} -.. automodule:: pyrenew.math - :members: - :undoc-members: - :show-inheritance: -``` - -## Transformations Module - -```{eval-rst} -.. automodule:: pyrenew.transformation - :members: - :undoc-members: - :show-inheritance: -``` - -## Regression Module - -```{eval-rst} -.. automodule:: pyrenew.regression - :members: - :undoc-members: - :show-inheritance: - -``` - -## Distributions Utility Module - -```{eval-rst} -.. automodule:: pyrenew.distutil - :members: - :undoc-members: - :show-inheritance: -``` - -```{eval-rst} -.. todo:: Determine and naming order of these modules. -``` diff --git a/docs/source/pyrenew_reference/datasets.md b/docs/source/pyrenew_reference/datasets.md deleted file mode 100644 index 7bdb8b28..00000000 --- a/docs/source/pyrenew_reference/datasets.md +++ /dev/null @@ -1,28 +0,0 @@ -# Datasets - -## Generation interval - -```{eval-rst} -.. automodule:: pyrenew.datasets.generation_interval - :members: - :undoc-members: - :show-inheritance: -``` - -## Infection to admission interval - -```{eval-rst} -.. automodule:: pyrenew.datasets.infection_admission_interval - :members: - :undoc-members: - :show-inheritance: -``` - -## Simulated Wastewater Data - -```{eval-rst} -.. automodule:: pyrenew.datasets.wastewater - :members: - :undoc-members: - :show-inheritance: -``` diff --git a/docs/source/pyrenew_reference/deterministic.md b/docs/source/pyrenew_reference/deterministic.md deleted file mode 100644 index 5edd9709..00000000 --- a/docs/source/pyrenew_reference/deterministic.md +++ /dev/null @@ -1,8 +0,0 @@ -# Deterministic - -```{eval-rst} -.. automodule:: pyrenew.deterministic - :members: - :undoc-members: - :show-inheritance: -``` diff --git a/docs/source/pyrenew_reference/distributions.md b/docs/source/pyrenew_reference/distributions.md deleted file mode 100644 index d959337e..00000000 --- a/docs/source/pyrenew_reference/distributions.md +++ /dev/null @@ -1,8 +0,0 @@ -# Distributions - -```{eval-rst} -.. automodule:: pyrenew.distributions - :members: - :undoc-members: - :show-inheritance: -``` diff --git a/docs/source/pyrenew_reference/index.md b/docs/source/pyrenew_reference/index.md deleted file mode 100644 index e24f88c7..00000000 --- a/docs/source/pyrenew_reference/index.md +++ /dev/null @@ -1,23 +0,0 @@ -# Reference - -```{toctree} -:maxdepth: 2 - -model -latent -process -randomvariable -observation -distributions -deterministic -datasets -additional -``` - -```{eval-rst} -.. todolist:: -``` - -```{eval-rst} -.. todo:: Determine order and naming of the above toctree entities. -``` diff --git a/docs/source/pyrenew_reference/latent.md b/docs/source/pyrenew_reference/latent.md deleted file mode 100644 index 632dacf8..00000000 --- a/docs/source/pyrenew_reference/latent.md +++ /dev/null @@ -1,50 +0,0 @@ -# Latent - -## Hospital Admissions - -```{eval-rst} -.. automodule:: pyrenew.latent.hospitaladmissions - :members: - :undoc-members: - :show-inheritance: -``` - -## Infections - -```{eval-rst} -.. automodule:: pyrenew.latent.infections - :members: - :undoc-members: - :show-inheritance: -``` - -## Infection Functions - -```{eval-rst} -.. automodule:: pyrenew.latent.infection_functions - :members: - :undoc-members: - :show-inheritance: -``` - -## Infection Initialization Process - -```{eval-rst} -.. automodule:: pyrenew.latent.infection_initialization_process - :members: - :undoc-members: - :show-inheritance: -``` - -## Infection Initialization Method - -```{eval-rst} -.. automodule:: pyrenew.latent.infection_initialization_method - :members: - :undoc-members: - :show-inheritance: -``` - -```{eval-rst} -.. todo:: Determine and naming order of these modules. -``` diff --git a/docs/source/pyrenew_reference/model.md b/docs/source/pyrenew_reference/model.md deleted file mode 100644 index cc28496f..00000000 --- a/docs/source/pyrenew_reference/model.md +++ /dev/null @@ -1,23 +0,0 @@ -# Model - -## Hospital Admissions - -```{eval-rst} -.. automodule:: pyrenew.model.admissionsmodel - :members: - :undoc-members: - :show-inheritance: -``` - -## Reproduction Number Renewal Infections - -```{eval-rst} -.. automodule:: pyrenew.model.rtinfectionsrenewalmodel - :members: - :undoc-members: - :show-inheritance: -``` - -```{eval-rst} -.. todo:: Determine and naming order of these modules. -``` diff --git a/docs/source/pyrenew_reference/observation.md b/docs/source/pyrenew_reference/observation.md deleted file mode 100644 index 86b3c650..00000000 --- a/docs/source/pyrenew_reference/observation.md +++ /dev/null @@ -1,23 +0,0 @@ -# Observation - -## Negative Binomial - -```{eval-rst} -.. automodule:: pyrenew.observation.negativebinomial - :members: - :undoc-members: - :show-inheritance: -``` - -## Poisson - -```{eval-rst} -.. automodule:: pyrenew.observation.poisson - :members: - :undoc-members: - :show-inheritance: -``` - -```{eval-rst} -.. todo:: Determine order of these modules. -``` diff --git a/docs/source/pyrenew_reference/process.md b/docs/source/pyrenew_reference/process.md deleted file mode 100644 index 25e0f4a6..00000000 --- a/docs/source/pyrenew_reference/process.md +++ /dev/null @@ -1,8 +0,0 @@ -# Random Process - -```{eval-rst} -.. automodule:: pyrenew.process - :members: - :undoc-members: - :show-inheritance: -``` diff --git a/docs/source/pyrenew_reference/randomvariable.md b/docs/source/pyrenew_reference/randomvariable.md deleted file mode 100644 index 1c43a1a1..00000000 --- a/docs/source/pyrenew_reference/randomvariable.md +++ /dev/null @@ -1,8 +0,0 @@ -# Random Variables - -```{eval-rst} -.. automodule:: pyrenew.randomvariable - :members: - :undoc-members: - :show-inheritance: -``` diff --git a/docs/source/tutorials/.gitignore b/docs/source/tutorials/.gitignore deleted file mode 100644 index d1fbb83b..00000000 --- a/docs/source/tutorials/.gitignore +++ /dev/null @@ -1,10 +0,0 @@ -*.png -*_files/ -!README.md - -/.quarto/ -_compiled - -*.ipynb -*.py -!*.rst diff --git a/docs/source/tutorials/basic_renewal_model.md b/docs/source/tutorials/basic_renewal_model.md deleted file mode 100644 index c1d9b7e1..00000000 --- a/docs/source/tutorials/basic_renewal_model.md +++ /dev/null @@ -1,5 +0,0 @@ -## Placeholder file - Please do not edit this file directly. - This file is just a placeholder. - For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/basic_renewal_model.qmd diff --git a/docs/source/tutorials/basic_renewal_model.qmd b/docs/source/tutorials/basic_renewal_model.qmd deleted file mode 100644 index 040001b2..00000000 --- a/docs/source/tutorials/basic_renewal_model.qmd +++ /dev/null @@ -1,407 +0,0 @@ ---- -title: "Fitting a basic renewal model" -format: gfm -engine: jupyter ---- - -This document will show the steps to build a simple renewal model featuring a latent infection process, a random walk Rt process, and an observation process for the reported infections. - -We start by loading the needed components to build a basic renewal model: - -```{python} -# | label: loading-pkgs -# | output: false -# | warning: false -import jax.numpy as jnp -import numpy as np -import numpyro -import numpyro.distributions as dist -from pyrenew.process import RandomWalk -from pyrenew.latent import ( - Infections, - InfectionInitializationProcess, - InitializeInfectionsZeroPad, -) -from pyrenew.observation import PoissonObservation -from pyrenew.deterministic import DeterministicPMF -from pyrenew.model import RtInfectionsRenewalModel -from pyrenew.metaclass import RandomVariable -from pyrenew.randomvariable import DistributionalVariable, TransformedVariable -import pyrenew.transformation as t -from numpyro.infer.reparam import LocScaleReparam -``` -By default, XLA (which is used by JAX for compilation) considers all CPU cores as one device. Depending on your system's configuration, we recommend using numpyro's [set_host_device_count()](https://num.pyro.ai/en/stable/utilities.html#set-host-device-count) function to set the number of devices available for parallel computing. Here, we set the device count to be 2. -```{python} -# | label: set-device-count -numpyro.set_host_device_count(2) -``` - -## Architecture of `RtInfectionsRenewalModel` - -`pyrenew` leverages `numpyro`'s flexibility to build models via composition. As a principle, most objects in `pyrenew` can be treated as random variables we can sample. At the top-level `pyrenew` has two metaclasses from which most objects inherit: `RandomVariable` and `Model`. From them, the following four sub-modules arise: - -- The `process` sub-module, -- The `deterministic` sub-module, -- The `observation` sub-module, -- The `latent` sub-module, and -- The `models` sub-module - -The first four are collections of instances of `RandomVariable`, and the last is a collection of instances of `Model`. The following diagram shows a detailed view of how metaclasses, modules, and classes interact to create the `RtInfectionsRenewalModel` instantiated in the next section: - -```mermaid -%%| label: overview-of-RtInfectionsRenewalModel -%%| include: true -flowchart LR - rand((RandomVariable\nmetaclass)) - models((Model\nmetaclass)) - - subgraph observations[Observations module] - obs["infection_obs_process_rv\n(PoissonObservation)"] - end - - subgraph latent[Latent module] - inf["latent_infections_rv\n(Infections)"] - i0["I0_rv\n(DistributionalVariable)"] - end - - subgraph process[Process module] - rt["Rt_process_rv\n(Custom class built using RandomWalk)"] - end - - subgraph deterministic[Deterministic module] - detpmf["gen_int_rv\n(DeterministicPMF)"] - end - - subgraph model[Model module] - model1["model1\n(RtInfectionsRenewalModel)"] - end - - rand-->|Inherited by|observations - rand-->|Inherited by|process - rand-->|Inherited by|latent - rand-->|Inherited by|deterministic - models-->|Inherited by|model - - detpmf-->|Composes|model1 - i0-->|Composes|model1 - rt-->|Composes|model1 - obs-->|Composes|model1 - inf-->|Composes|model1 -``` - - -The pyrenew package models the real-time reproductive number $\mathcal{R}(t)$, the average number of secondary infections caused by an infected individual, as a renewal process model. Our basic renewal process model defines five components: - -(1) generation interval, the times between infections - -(2) initial infections, occurring prior to time $t = 0$ - -(3) $\mathcal{R}(t)$, the time-varying reproductive number, - -(4) latent infections, i.e., those infections which are known to exist but are not observed (or not observable), and - -(5) observed infections, a subset of underlying true infections that are reported, perhaps via hospital admissions, physician's office visits, or routine biosurveillance. - -To initialize these five components within the renewal modeling framework, we estimate each component with: - -(1) In this example, the generation interval is not estimated but passed as a deterministic instance of `RandomVariable` - -(2) an instance of the `InfectionInitializationProcess` class, where the number of latent infections immediately before the renewal process begins follows a log-normal distribution with mean = 0 and standard deviation = 1. By specifying `InitializeInfectionsZeroPad`, the latent infections before this time are assumed to be 0. - -(3) A process to represent $\mathcal{R}(t)$ as a random walk on the log scale, with an inferred initial value and a fixed Normal step-size distribution. For this, we construct a custom `RandomVariable`, `MyRt`. - -(4) an instance of the `Infections` class with default values, and - -(5) an instance of the `PoissonObservation` class with default values - -```{python} -# | label: creating-elements -# (1) The generation interval (deterministic) -pmf_array = jnp.array([0.4, 0.3, 0.2, 0.1]) -gen_int = DeterministicPMF(name="gen_int", value=pmf_array) - -# (2) Initial infections (inferred with a prior) -I0 = InfectionInitializationProcess( - "I0_initialization", - DistributionalVariable(name="I0", distribution=dist.LogNormal(2.5, 1)), - InitializeInfectionsZeroPad(pmf_array.size), -) - - -# (3) The random walk on log Rt, with an inferred s.d. Here, we -# construct a custom RandomVariable. -class MyRt(RandomVariable): - def validate(self): - pass - - def sample(self, n: int, **kwargs) -> tuple: - sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) - - rt_rv = TransformedVariable( - name="log_rt_random_walk", - base_rv=RandomWalk( - name="log_rt", - step_rv=DistributionalVariable( - name="rw_step_rv", distribution=dist.Normal(0, sd_rt) - ), - ), - transforms=t.ExpTransform(), - ) - rt_init_rv = DistributionalVariable( - name="init_log_rt", distribution=dist.Normal(0, 0.2) - ) - init_rt = rt_init_rv.sample() - - return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) - - -rt_proc = MyRt() - -# (4) Latent infection process (which will use 1 and 2) -latent_infections = Infections() - -# (5) The observed infections process (with mean at the latent infections) -observation_process = PoissonObservation("poisson_rv") -``` - -With these five pieces, we can build the basic renewal model as an instance of the `RtInfectionsRenewalModel` class: - -```{python} -# | label: model-creation -model1 = RtInfectionsRenewalModel( - gen_int_rv=gen_int, - I0_rv=I0, - Rt_process_rv=rt_proc, - latent_infections_rv=latent_infections, - infection_obs_process_rv=observation_process, -) -``` - -The following diagram summarizes how the modules interact via composition; notably, `gen_int`, `I0`, `rt_proc`, `latent_infections`, and `observed_infections` are instances of `RandomVariable`, which means these can be easily replaced to generate a different instance of the `RtInfectionsRenewalModel` class: - -```mermaid -%%| label: overview-of-RtInfectionsRenewalModel -%%| include: true -flowchart TB - genint["(1) gen_int\n(DetermnisticPMF)"] - i0["(2) I0\n(InfectionInitializationProcess)"] - rt["(3) rt_proc\n(MyRt, the custom RV defined above)"] - inf["(4) latent_infections\n(Infections)"] - obs["(5) observation_process\n(PoissonObservation)"] - - model1["model1\n(RtInfectionsRenewalModel)"] - - i0-->|Composes|model1 - genint-->|Composes|model1 - rt-->|Composes|model1 - obs-->|Composes|model1 - inf-->|Composes|model1 -``` - -Using `numpyro`, we can simulate data using the `sample()` member function of `RtInfectionsRenewalModel`. The `sample()` method of the `RtInfectionsRenewalModel` class returns a list composed of the `Rt` and `infections` sequences, called `sim_data`: - -```{python} -# | label: simulate -with numpyro.handlers.seed(rng_seed=53): - sim_data = model1.sample(n_datapoints=40) - -sim_data -``` - -To understand what has been accomplished here, visualize an $\mathcal{R}(t)$ sample path (left panel) and infections over time (right panel): - -```{python} -# | label: fig-basic -# | fig-cap: Rt and Infections -import matplotlib.pyplot as plt - -fig, axs = plt.subplots(1, 2) - -# Rt plot -axs[0].plot(sim_data.Rt) -axs[0].set_ylabel("Rt") - -# Infections plot -axs[1].plot(sim_data.observed_infections) -axs[1].set_ylabel("Infections") - -fig.suptitle("Basic renewal model") -fig.supxlabel("Time") -plt.tight_layout() -plt.show() -``` - -To fit the model, we can use the `run()` method of the `RtInfectionsRenewalModel` class (an inherited method from the metaclass `Model`). `model1.run()` will call the `run` method of the `model1` object, which will generate an instance of model MCMC simulation, with 2000 warm-up iterations for the MCMC algorithm, used to tune the parameters of the MCMC algorithm to improve efficiency of the sampling process. From the posterior distribution of the model parameters, 1000 samples will be drawn and used to estimate the posterior distributions and compute summary statistics. Observed data is provided to the model using the `sim_data` object previously generated. `mcmc_args` provides additional arguments for the MCMC algorithm. - -```{python} -# | label: model-fit -import jax - -model1.run( - num_warmup=2000, - num_samples=1000, - data_observed_infections=sim_data.observed_infections, - rng_key=jax.random.key(54), - mcmc_args=dict(progress_bar=False, num_chains=2), -) -``` - -Now, let's investigate the output, particularly the posterior distribution of the $\mathcal{R}(t)$ estimates: - -```{python} -# | label: fig-output-rt -# | fig-cap: Rt posterior distribution -import arviz as az - -# Create arviz inference data object -idata = az.from_numpyro( - posterior=model1.mcmc, -) - -# Extract Rt signal samples across chains -rt = az.extract(idata.posterior["Rt"], num_samples=100)["Rt"].values - - -# Plot Rt signal -fig, ax = plt.subplots(1, 1, figsize=(8, 6)) -ax.plot( - np.arange(rt.shape[0]), - rt, - color="skyblue", - alpha=0.10, -) -ax.plot([], [], color="skyblue", alpha=0.05, label="Rt Posterior Samples") -ax.plot( - np.arange(rt.shape[0]), - rt.mean(axis=1), - color="black", - linewidth=2.0, - linestyle="--", - label="Sample Mean", -) -ax.legend(loc="best") -ax.set_ylabel(r"$\mathscr{R}_t$ Signal", fontsize=20) -ax.set_xlabel("Days", fontsize=20) -plt.show() -``` -We can use the `get_samples` method to extract samples from the model -```{python} -# | label: get-model-fitted-samples -Rt_samp = model1.mcmc.get_samples()["Rt"] -latent_infection_samp = model1.mcmc.get_samples()["all_latent_infections"] -``` - -We can also convert the fitted model to [ArviZ](https://www.arviz.org/) InferenceData object and use ArviZ package to extarct samples, calculate statistics, create model diagnostics and visualizations. - -```{python} -# | label: convert-inference-data -import arviz as az - -idata = az.from_numpyro(model1.mcmc) -``` - -and use the InferenceData to compute the model-fit diagnostics. Here, we show diagnostic summary for the first 10 effective reproduction number $\mathcal{R}(t)$. - -```{python} -# | label: diagnostics -diagnostic_stats_summary = az.summary( - idata.posterior["Rt"][::, ::, 4:], # ignore nan padding - kind="diagnostics", -) - -print(diagnostic_stats_summary) -``` - -Below we use `plot_trace` to inspect the trace of the first 10 inferred $\mathcal{R}(t)$ values. - -```{python} -# | label: fig-trace-Rt -# | fig-cap: Trace plot of Rt posterior distribution -plt.rcParams["figure.constrained_layout.use"] = True - -az.plot_trace( - idata.posterior, - var_names=["Rt"], - coords={"Rt_dim_0": np.arange(4, 14)}, - compact=False, -) -plt.show() -``` - - -We inspect the posterior distribution of $\mathcal{R}(t)$ by plotting the 90% and 50% highest density intervals: - -```{python} -# | label: fig-hdi-Rt -# | fig-cap: High density interval for Effective Reproduction Number -x_data = idata.posterior["Rt_dim_0"][4:] -y_data = idata.posterior["Rt"][::, ::, 4:] - -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.9, - color="C0", - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.5, - color="C0", - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -# Add mean of the posterior to the figure -median_ts = y_data.median(dim=["chain", "draw"]) -axes.plot(x_data, median_ts, color="C0", label="Median") -axes.legend() -axes.set_title("Posterior Effective Reproduction Number", fontsize=10) -axes.set_xlabel("Time", fontsize=10) -axes.set_ylabel("$\\mathcal{R}(t)$", fontsize=10) -plt.show() -``` - -and latent infections: - -```{python} -# | label: fig-hdi-latent-infections -# | fig-cap: High density interval for Latent Infections -x_data = idata.posterior["all_latent_infections_dim_0"] -y_data = idata.posterior["all_latent_infections"] - -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.9, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.5, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -# plot the posterior median -median_ts = y_data.median(dim=["chain", "draw"]) -axes.plot(x_data, median_ts, color="C0", label="Median") - -axes.legend() -axes.set_title("Posterior Latent Infections", fontsize=10) -axes.set_xlabel("Time", fontsize=10) -axes.set_ylabel("Latent Infections", fontsize=10) -plt.show() -``` diff --git a/docs/source/tutorials/day_of_the_week.qmd b/docs/source/tutorials/day_of_the_week.qmd deleted file mode 100644 index 018ce4e4..00000000 --- a/docs/source/tutorials/day_of_the_week.qmd +++ /dev/null @@ -1,426 +0,0 @@ ---- -title: Implementing a Day of the Week Effect -format: gfm -engine: jupyter ---- - -This document illustrates how to leverage the time-aware arrays to create a a day of the week effect. We use the same model designed in the hospital admissions-only tutorial. - -## Recap: Hospital Admissions Model - -In the ["Fitting a hospital admissions-only model" tutorial](https://cdcgov.github.io/PyRenew/tutorials/hospital_admissions_model.html){target="_blank"}, we built a fairly complicated model that included pre-defined random variables as well as a custom random variable representing the reproductive number. In this tutorial, we will focus on adding a new component: the day-of-the-week effect. We start by reproducing the model without the day-of-the-week effect: - -1. We load the data: - -```{python} -# | label: setup -# | code-fold: true -# Setup -import numpyro -import polars as pl -from pyrenew import datasets - -# Setting the number of devices -numpyro.set_host_device_count(2) - -# Loading and processing the data -dat = ( - datasets.load_wastewater() - .group_by("date") - .first() - .select(["date", "daily_hosp_admits"]) - .sort("date") - .head(90) -) - -daily_hosp_admits = dat["daily_hosp_admits"].to_numpy() -dates = dat["date"].to_numpy() - -# Loading additional datasets -gen_int = datasets.load_generation_interval() -inf_hosp_int = datasets.load_infection_admission_interval() - -# We only need the probability_mass column of each dataset -gen_int_array = gen_int["probability_mass"].to_numpy() -gen_int = gen_int_array -inf_hosp_int_array = inf_hosp_int["probability_mass"].to_numpy() -``` - -2. Next, we defined the model's components: - -```{python} -# | label: latent-hosp -# | code-fold: true -from pyrenew import latent, deterministic, randomvariable -import jax.numpy as jnp -import numpyro.distributions as dist - -inf_hosp_int = deterministic.DeterministicPMF( - name="inf_hosp_int", value=inf_hosp_int_array -) - -hosp_rate = randomvariable.DistributionalVariable( - name="IHR", distribution=dist.LogNormal(jnp.log(0.05), jnp.log(1.1)) -) - -latent_hosp = latent.HospitalAdmissions( - infection_to_admission_interval_rv=inf_hosp_int, - infection_hospitalization_ratio_rv=hosp_rate, -) - -from pyrenew import model, process, observation, metaclass, transformation -from pyrenew.latent import ( - InfectionInitializationProcess, - InitializeInfectionsExponentialGrowth, -) - - -# Infection process -latent_inf = latent.Infections() -n_initialization_points = max(gen_int_array.size, inf_hosp_int_array.size) - 1 - -I0 = InfectionInitializationProcess( - "I0_initialization", - randomvariable.DistributionalVariable( - name="I0", - distribution=dist.LogNormal(loc=jnp.log(100), scale=jnp.log(1.75)), - ), - InitializeInfectionsExponentialGrowth( - n_initialization_points, - deterministic.DeterministicVariable(name="rate", value=0.05), - ), -) - -# Generation interval and Rt -gen_int = deterministic.DeterministicPMF(name="gen_int", value=gen_int) -``` - -including the Rt effect: - -```{python} -# | label: Rt-process -# | code-fold: true -class MyRt(metaclass.RandomVariable): - def __init__(self, sd_rv): - self.sd_rv = sd_rv - - def validate(self): - pass - - def sample(self, n: int, **kwargs) -> tuple: - # Standard deviation of the random walk - sd_rt = self.sd_rv() - - # Random walk step - step_rv = randomvariable.DistributionalVariable( - name="rw_step_rv", distribution=dist.Normal(0, sd_rt) - ) - - rt_init_rv = randomvariable.DistributionalVariable( - name="init_log_rt", distribution=dist.Normal(0, 0.2) - ) - - # Random walk process - base_rv = process.RandomWalk( - name="log_rt", - step_rv=step_rv, - ) - - # Transforming the random walk to the Rt scale - rt_rv = randomvariable.TransformedVariable( - name="Rt_rv", - base_rv=base_rv, - transforms=transformation.ExpTransform(), - ) - init_rt = rt_init_rv.sample() - - return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) - - -rtproc = MyRt( - randomvariable.DistributionalVariable( - name="Rt_random_walk_sd", distribution=dist.HalfNormal(0.025) - ) -) -``` - -3. We defined the observation model: - -```{python} -# | label: obs-model -# | code-fold: true -# we place a log-Normal prior on the concentration -# parameter of the negative binomial. -nb_conc_rv = randomvariable.TransformedVariable( - "concentration", - randomvariable.DistributionalVariable( - name="concentration_raw", - distribution=dist.TruncatedNormal(loc=0, scale=1, low=0.01), - ), - transformation.PowerTransform(-2), -) - -# now we define the observation process -obs = observation.NegativeBinomialObservation( - "negbinom_rv", - concentration_rv=nb_conc_rv, -) -``` - -4. And finally, we build the model: - -```{python} -# | label: init-model -# | code-fold: true -hosp_model = model.HospitalAdmissionsModel( - latent_infections_rv=latent_inf, - latent_hosp_admissions_rv=latent_hosp, - I0_rv=I0, - gen_int_rv=gen_int, - Rt_process_rv=rtproc, - hosp_admission_obs_process_rv=obs, -) -``` - -Here is what the model looks like without the day-of-the-week effect: - -```{python} -# | label: fig-output-admissions-padding-and-weekday -# | fig-cap: Hospital Admissions posterior distribution without weekday effect -# | code-fold: true -import jax -import numpy as np - -# Model without weekday effect -hosp_model.run( - num_samples=2000, - num_warmup=2000, - data_observed_hosp_admissions=daily_hosp_admits, - rng_key=jax.random.key(54), - mcmc_args=dict(progress_bar=False), -) -``` - -```{python} -# | code-fold: true -import arviz as az -import matplotlib.pyplot as plt - - -# Retrieve the posterior samples from the model -ppc_samples = hosp_model.posterior_predictive( - n_datapoints=daily_hosp_admits.size -) - -# Create an InferenceData object from model -idata = az.from_numpyro( - posterior=hosp_model.mcmc, - posterior_predictive=ppc_samples, -) - -# Use a time series plot (plot_ts) from arviz for plotting -axes = az.plot_ts( - idata, - y="negbinom_rv", - y_hat="negbinom_rv", - num_samples=200, - y_kwargs={ - "color": "blue", - "linewidth": 1.0, - "marker": "o", - "linestyle": "solid", - }, - y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, - y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, - backend_kwargs={"figsize": (8, 6)}, - textsize=15.0, -) -ax = axes[0][0] -ax.set_xlabel("Time", fontsize=20) -ax.set_ylabel("Hospital Admissions", fontsize=20) -handles, labels = ax.get_legend_handles_labels() -ax.legend( - handles, ["Observed", "Sample Mean", "Posterior Samples"], loc="best" -) -plt.show() -``` - - - - - -## Round 2: Incorporating day-of-the-week effects - -We will re-use the infection to admission interval and infection to hospitalization rate from the previous model. But we will also add a day-of-the-week effect. To do this, we will add two additional arguments to the latent hospital admissions random variable: `day_of_the_week_rv` (a `RandomVariable`) and `obs_data_first_day_of_the_week` (an `int` mapping days of the week from 0:6, zero being Monday). The `day_of_the_week_rv`'s sample method should return a vector of length seven; those values are then broadcasted to match the length of the dataset. Moreover, since the observed data may start in a weekday other than Monday, the `obs_data_first_day_of_the_week` argument is used to offset the day-of-the-week effect. - -For this example, the effect will be passed as a scaled Dirichlet distribution. It will consist of a `TransformedVariable` that samples an array of length seven from numpyro's `distributions.Dirichlet` and applies a `transformation.AffineTransform` to scale it by seven. [^note-other-examples]: - -[^note-other-examples]: A similar weekday effect is implemented in its own module, with example code here: {doc}`periodic_effects` - -```{python} -# | label: weekly-effect -# Instantiating the day-of-the-week effect -dayofweek_effect = randomvariable.TransformedVariable( - name="dayofweek_effect", - base_rv=randomvariable.DistributionalVariable( - name="dayofweek_effect_raw", - distribution=dist.Dirichlet(jnp.ones(7)), - ), - transforms=transformation.AffineTransform( - loc=0, scale=7, domain=jnp.array([0, 1]) - ), -) -``` - -Now, by re-defining the latent hospital admissions random variable with the day-of-the-week effect, we can build a model that includes this effect. Since the day-of-the-week effect takes into account the first day of the dataset, we need to determine the day of the week of the first observation. We can do this by converting the first date in the dataset to a `datetime` object and extracting the day of the week: - -```{python} -# | label: latent-hosp-weekday -# Figuring out the day of the week of the first observation -import datetime as dt - -first_dow_in_data = dates[0].astype(dt.datetime).weekday() -first_dow_in_data # zero - -# Re-defining the latent hospital admissions RV, now with the -# day-of-the-week effect -latent_hosp_wday_effect = latent.HospitalAdmissions( - infection_to_admission_interval_rv=inf_hosp_int, - infection_hospitalization_ratio_rv=hosp_rate, - day_of_week_effect_rv=dayofweek_effect, - # Concidirently, this is zero - obs_data_first_day_of_the_week=first_dow_in_data, -) - -# New model with day-of-the-week effect -hosp_model_dow = model.HospitalAdmissionsModel( - latent_infections_rv=latent_inf, - latent_hosp_admissions_rv=latent_hosp_wday_effect, - I0_rv=I0, - gen_int_rv=gen_int, - Rt_process_rv=rtproc, - hosp_admission_obs_process_rv=obs, -) -``` - -Running the model: - -```{python} -# | label: model-2-run-weekday -# Model with weekday effect -hosp_model_dow.run( - num_samples=2000, - num_warmup=2000, - data_observed_hosp_admissions=daily_hosp_admits, - rng_key=jax.random.key(54), - mcmc_args=dict(progress_bar=False), -) -``` - -As a result, we can see the posterior distribution of our novel day-of-the-week effect: - -```{python} -# | label: fig-output-day-of-week -# | fig-cap: Day of the week effect -# Create an InferenceData object from hosp_model_dow -dow_idata = az.from_numpyro( - posterior=hosp_model_dow.mcmc, -) - -# Extract day of week effect (DOW) -dow_effect_raw = dow_idata.posterior["dayofweek_effect_raw"].squeeze().T -indices = np.random.choice(dow_effect_raw.shape[1], size=200, replace=False) -dow_plot_samples = dow_effect_raw[:, indices] -fig, ax = plt.subplots(1, 1, figsize=(8, 6)) -ax.plot( - np.arange(dow_effect_raw.shape[0]), - dow_plot_samples, - color="skyblue", - alpha=0.10, -) -ax.plot([], [], color="skyblue", alpha=0.10, label="DOW Posterior Samples") -ax.plot( - np.arange(dow_effect_raw.shape[0]), - dow_plot_samples.mean(dim="draw"), - color="black", - linewidth=2.0, - linestyle="--", - label="Sample Mean", -) -ax.legend(loc="best") -ax.set_ylabel("Effect", fontsize=20) -ax.set_xlabel("Day Of Week", fontsize=20) -plt.show() -``` - -The new model with the day-of-the-week effect can be compared to the previous model without the effect. Finally, let's reproduce the figure without the day-of-the-week effect, and then plot the new model with the effect: - -```{python} -# | label: fig-output-admissions-original -# | fig-cap: Hospital Admissions posterior distribution without weekday effect -# Without weekday effect (from earlier) -axes = az.plot_ts( - idata, - y="negbinom_rv", - y_hat="negbinom_rv", - num_samples=200, - y_kwargs={ - "color": "blue", - "linewidth": 1.0, - "marker": "o", - "linestyle": "solid", - }, - y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, - y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, - backend_kwargs={"figsize": (8, 6)}, - textsize=15.0, -) -ax = axes[0][0] -ax.set_xlabel("Time", fontsize=20) -ax.set_ylabel("Hospital Admissions", fontsize=20) -handles, labels = ax.get_legend_handles_labels() -ax.legend( - handles, - ["Observed", "Posterior Predictive", "Samples wo/ WDE"], - loc="best", -) -plt.show() -``` - -```{python} -# | label: fig-output-admissions-wof -# | fig-cap: Hospital Admissions posterior distribution with weekday effect -# Figure with weekday effect -ppc_samples = hosp_model_dow.posterior_predictive( - n_datapoints=daily_hosp_admits.size -) -idata = az.from_numpyro( - posterior=hosp_model_dow.mcmc, - posterior_predictive=ppc_samples, -) - -axes = az.plot_ts( - idata, - y="negbinom_rv", - y_hat="negbinom_rv", - num_samples=200, - y_kwargs={ - "color": "blue", - "linewidth": 1.0, - "marker": "o", - "linestyle": "solid", - }, - y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, - y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, - backend_kwargs={"figsize": (8, 6)}, - textsize=15.0, -) -ax = axes[0][0] -ax.set_xlabel("Time", fontsize=20) -ax.set_ylabel("Hospital Admissions", fontsize=20) -handles, labels = ax.get_legend_handles_labels() -ax.legend( - handles, ["Observed", "Posterior Predictive", "Samples w/ WDE"], loc="best" -) -plt.show() -``` diff --git a/docs/source/tutorials/extending_pyrenew.qmd b/docs/source/tutorials/extending_pyrenew.qmd deleted file mode 100644 index 77a83a53..00000000 --- a/docs/source/tutorials/extending_pyrenew.qmd +++ /dev/null @@ -1,297 +0,0 @@ ---- -title: Extending pyrenew -format: gfm ---- - -This tutorial illustrates how to extend `pyrenew` with custom `RandomVariable` classes. We will use the `InfectionsWithFeedback` class as an example. The `InfectionsWithFeedback` class is a `RandomVariable` that models the number of infections at time $t$ as a function of the number of infections at time $t - \tau$ and the reproduction number at time $t$. The reproduction number at time $t$ is a function of the *unadjusted* reproduction number at time $t - \tau$ and the number of infections at time $t - \tau$: - -$$ -\begin{align*} -I(t) & = \mathcal{R}(t)\sum_{\tau=1}^{T_g}I(t - \tau)g(\tau) \\ -\mathcal{R}(t) & = \mathcal{R}^u(t)\exp\left(-\gamma(t)\sum_{\tau=1}^{T_f}I(t - \tau)f(\tau)\right) -\end{align*} -$$ - -Where $\mathcal{R}^u(t)$ is the unadjusted reproduction number, $g(t)$ is the generation interval, $\gamma(t)$ is the infection feedback strength, and $f(t)$ is the infection feedback pmf. - -## The expected outcome - -Before we start, let's simulate the model with the original `InfectionsWithFeedback` class. To keep it simple, we will simulate the model with no observation process, in other words, only with latent infections. The following code-chunk loads the required libraries and defines the model components: - -```{python} -# | label: setup -import jax -import jax.numpy as jnp -import numpy as np -import numpyro -import numpyro.distributions as dist -from pyrenew.deterministic import DeterministicPMF, DeterministicVariable -from pyrenew.latent import InfectionsWithFeedback -from pyrenew.model import RtInfectionsRenewalModel -from pyrenew.process import RandomWalk -from pyrenew.metaclass import RandomVariable -from pyrenew.randomvariable import DistributionalVariable, TransformedVariable -from pyrenew.latent import ( - InfectionInitializationProcess, - InitializeInfectionsExponentialGrowth, -) -import pyrenew.transformation as t -``` - -The following code-chunk defines the model components. Notice that for both the generation interval and the infection feedback, we use a deterministic PMF with equal probabilities: - -```{python} -# | label: model-components -gen_int_array = jnp.array([0.25, 0.5, 0.15, 0.1]) - -gen_int = DeterministicPMF(name="gen_int", value=gen_int_array) -feedback_strength = DeterministicVariable(name="feedback_strength", value=0.01) - - -I0 = InfectionInitializationProcess( - "I0_initialization", - DistributionalVariable(name="I0", distribution=dist.LogNormal(0, 1)), - InitializeInfectionsExponentialGrowth( - gen_int_array.size, - DeterministicVariable(name="rate", value=0.05), - ), -) - -latent_infections = InfectionsWithFeedback( - infection_feedback_strength=feedback_strength, - infection_feedback_pmf=gen_int, -) - - -class MyRt(RandomVariable): - def validate(self): - pass - - def sample(self, n: int, **kwargs) -> tuple: - sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) - - rt_rv = TransformedVariable( - name="log_rt_random_walk", - base_rv=RandomWalk( - name="log_rt", - step_rv=DistributionalVariable( - name="rw_step_rv", distribution=dist.Normal(0, 0.025) - ), - ), - transforms=t.ExpTransform(), - ) - rt_init_rv = DistributionalVariable( - name="init_log_rt", distribution=dist.Normal(0, 0.2) - ) - init_rt = rt_init_rv.sample() - - return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) -``` - -With all the components defined, we can build the model: - -```{python} -# | label: build1 -model0 = RtInfectionsRenewalModel( - gen_int_rv=gen_int, - I0_rv=I0, - latent_infections_rv=latent_infections, - Rt_process_rv=MyRt(), - infection_obs_process_rv=None, -) -``` - -And simulate from it: - -```{python} -# | label: simulate1 -# Sampling and fitting model 0 (with no obs for infections) -with numpyro.handlers.seed(rng_seed=223): - model0_samp = model0.sample(n_datapoints=30) -``` - -```{python} -# | label: fig-simulate1 -# | fig-cap: Simulated infections with no observation process -import matplotlib.pyplot as plt - -fig, ax = plt.subplots() -ax.plot(model0_samp.latent_infections) -ax.set_xlabel("Time") -ax.set_ylabel("Infections") -plt.show() -``` - -## Pyrenew's random variable class - -### Fundamentals - -All instances of PyRenew's `RandomVariable` should have at least three functions: `__init__()`, `validate()`, and `sample()`. The `__init__()` function is the constructor and initializes the class. The `validate()` function checks if the class is correctly initialized. Finally, the `sample()` method contains the core of the class; it should return a tuple or named tuple. The following is a minimal example of a `RandomVariable` class based on `numpyro.distributions.Normal`: - -```{python} -from pyrenew.metaclass import RandomVariable - - -class MyNormal(RandomVariable): - def __init__(self, loc, scale): - self.validate(scale) - self.loc = loc - self.scale = scale - return None - - @staticmethod - def validate(self): - if self.scale <= 0: - raise ValueError("Scale must be positive") - return None - - def sample(self, **kwargs): - return (dist.Normal(loc=self.loc, scale=self.scale),) -``` - -The `@staticmethod` decorator exposes the `validate` function to be used outside the class. Next, we show how to build a more complex `RandomVariable` class; the `InfectionsWithFeedback` class. - -### The `InfectionsWithFeedback` class - -Although returning namedtuples is not strictly required, they are the recommended return type, as they make the code more readable. The following code-chunk shows how to create a named tuple for the `InfectionsWithFeedback` class: - -```{python} -# | label: data-class -from collections import namedtuple - -# Creating a tuple to store the output -InfFeedbackSample = namedtuple( - typename="InfFeedbackSample", - field_names=["post_initialization_infections", "rt"], - defaults=(None, None), -) -``` - -The next step is to create the actual class. The bulk of its implementation lies in the function `pyrenew.latent.compute_infections_from_rt_with_feedback()`. We will also use the `pyrenew.arrayutils.pad_edges_to_match()` function to ensure the passed vectors match their lengths. The following code-chunk shows most of the implementation of the `InfectionsWithFeedback` class: - -```{python} -# | label: new-model-def -# | code-line-numbers: true -# Creating the class -from pyrenew.metaclass import RandomVariable -from pyrenew.latent import compute_infections_from_rt_with_feedback -from pyrenew import arrayutils as au -from jax.typing import ArrayLike -import jax.numpy as jnp - - -class InfFeedback(RandomVariable): - """Latent infections""" - - def __init__( - self, - infection_feedback_strength: RandomVariable, - infection_feedback_pmf: RandomVariable, - ) -> None: - """Constructor""" - - self.infection_feedback_strength = infection_feedback_strength - self.infection_feedback_pmf = infection_feedback_pmf - - return None - - def validate(self): - """ - Generally, this method should be more meaningful, but we will skip it for now - """ - return None - - def sample( - self, - Rt: ArrayLike, - I0: ArrayLike, - gen_int: ArrayLike, - **kwargs, - ) -> tuple: - """Sample infections with feedback""" - - # Generation interval - gen_int_rev = jnp.flip(gen_int) - - # Baseline infections - I0_vec = I0[-gen_int_rev.size :] - - # Sampling inf feedback strength and adjusting the shape - inf_feedback_strength = self.infection_feedback_strength( - **kwargs, - ) - - inf_feedback_strength = jnp.atleast_1d(inf_feedback_strength) - - inf_feedback_strength, _ = au.pad_edges_to_match( - x=inf_feedback_strength, - y=Rt, - ) - - # Sampling inf feedback and adjusting the shape - inf_feedback_pmf = self.infection_feedback_pmf(**kwargs) - inf_fb_pmf_rev = jnp.flip(inf_feedback_pmf) - - # Generating the infections with feedback - all_infections, Rt_adj = compute_infections_from_rt_with_feedback( - I0=I0_vec, - Rt_raw=Rt, - infection_feedback_strength=inf_feedback_strength, - reversed_generation_interval_pmf=gen_int_rev, - reversed_infection_feedback_pmf=inf_fb_pmf_rev, - ) - - # Storing adjusted Rt for future use - numpyro.deterministic("Rt_adjusted", Rt_adj) - - # Preparing theoutput - - return InfFeedbackSample( - post_initialization_infections=all_infections, - rt=Rt_adj, - ) -``` - -The core of the class is implemented in the `sample()` method. Things to highlight from the above code: - -1. **Arguments of `sample`**: The `InfFeedback` class will be used within `RtInfectionsRenewalModel` to generate latent infections. During the sampling process, `InfFeedback()` will receive the reproduction number, the initial number of infections, and the generation interval. `RandomVariable()` calls are expected to include the `**kwargs` argument, even if unused. - -2. **Saving computed quantities**: Since `Rt_adj` is not generated via `numpyro.sample()`, we use `numpyro.deterministic()` to record the quantity to a site; allowing us to access it later. - -3. **Return type of `InfFeedback()`**: As said before, the `sample()` method should return a tuple or named tuple. In our case, we return a named tuple `InfFeedbackSample` with two fields: `infections` and `rt`. - -```{python} -# | label: simulation2 -latent_infections2 = InfFeedback( - infection_feedback_strength=feedback_strength, - infection_feedback_pmf=gen_int, -) - -model1 = RtInfectionsRenewalModel( - gen_int_rv=gen_int, - I0_rv=I0, - latent_infections_rv=latent_infections2, - Rt_process_rv=MyRt(), - infection_obs_process_rv=None, -) - -# Sampling and fitting model 0 (with no obs for infections) -with numpyro.handlers.seed(rng_seed=223): - model1_samp = model1.sample(n_datapoints=30) -``` - -Comparing `model0` with `model1`, these two should match: - -```{python} -# | label: fig-model0-vs-model1 -# | fig-cap: Comparing latent infections from model 0 and model 1 -import matplotlib.pyplot as plt - -fig, ax = plt.subplots(ncols=2) -ax[0].plot(model0_samp.latent_infections) -ax[1].plot(model1_samp.latent_infections) -ax[0].set_xlabel("Time (model 0)") -ax[1].set_xlabel("Time (model 1)") -ax[0].set_ylabel("Infections") -plt.show() -``` diff --git a/docs/source/tutorials/getting_started.qmd b/docs/source/tutorials/getting_started.qmd deleted file mode 100644 index 24490f91..00000000 --- a/docs/source/tutorials/getting_started.qmd +++ /dev/null @@ -1,60 +0,0 @@ ---- -title: "Getting started with pyrenew" -format: gfm -engine: jupyter ---- - -`pyrenew` is a flexible tool for simulating and making statistical inferences from epidemiologic models, with an emphasis on renewal models. Built on `numpyro`, `pyrenew` provides core components for model building and pre-defined models for processing various observational processes. - -## Prerequisites - -This tutorial assumes some pre-existing knowledge of infectious disease dynamics and Python programming. Before you dive in, we recommend: - - - Installing Python3 (use tools like [pyenv](https://realpython.com/intro-to-pyenv/) or [compile and install](https://ubuntuhandbook.org/index.php/2023/05/install-python-3-12-ubuntu/) from the [release page](https://www.python.org/downloads/)) - - Familiarity with installing and loading modules in python, and with virtual environment management (we recommend [uv](https://docs.astral.sh/uv/)) - - Familiarity with the concept of a [class](https://realpython.com/python-classes/) and [metaclass](https://realpython.com/courses/python-metaclasses/) in python - - Familiarity with Bayesian inference, and a working understanding of MCMC methods used to fit Bayesian models to data (some resources are available [here](https://mc-stan.org/docs/2_18/reference-manual/effective-sample-size-section.html), and [here](https://xcelab.net/rm/)) - -## Installing pyrenew - -You’ll need to install `pyrenew` using either uv or pip. To install `pyrenew` using `uv`, run the following command from within the directory containing the `pyrenew` project: - -```bash -uv sync -``` - -To install `pyrenew` using pip, run the following command: - -```bash -pip install git+https://github.com/CDCgov/PyRenew@main -``` - - -## The fundamentals - -`pyrenew`'s core components are the metaclasses `RandomVariable` and `Model` (in Python, a _metaclass_ is a class whose instances are also classes, where a _class_ is a template for making objects). Within the `pyrenew` package, a `RandomVariable` is a quantity that models can estimate and sample from, **including deterministic quantities**. The benefit of this design is that the definition of the `sample()` function can be arbitrary, allowing the user to either sample from a distribution using `numpyro.sample()`, compute fixed quantities (like a mechanistic equation), or return a fixed value (like a pre-computed PMF.) For instance, when estimating a PMF, the `RandomVariable` sampling function may roughly be defined as: - -```python -# define a new class called MyRandVar that inherits from the RandomVariable class -class MyRandVar(RandomVariable): - #define a method called sample that returns an object of type ArrayLike - def sample(...) -> ArrayLike: - # calls sample function from NumPyro package - return numpyro.sample(...) -``` - -Whereas, in some other cases, we may instead use a fixed quantity for that variable (like a pre-computed PMF), where the `RandomVariable`'s sample function could instead be defined as: - -```python -# instead define MyRandVar to still inherit from the RandVariable class -class MyRandVar(RandomVariable): - #define sample method that still returns an ArrayLike object - def sample(...) -> ArrayLike: - #sampling method is a pre-computed PMF, a JAX NumPy array with explicit elements - return jax.numpy.array([0.2, 0.7, 0.1]) -``` - -Thus, when a `Model` samples from `MyRandVar`, it could be either adding random variables to be estimated (first case) or just retrieving some quantity needed for other calculations (second case.) - -The `Model` metaclass provides basic functionality for estimating and simulation. Like `RandomVariable`, the `Model` metaclass has a `sample()` method that defines the model structure. Ultimately, models can be nested (or inherited), providing a straightforward way to add layers of complexity. -At this stage, the `Model` metaclass consist of two model classes `RtInfectionsRenewalModel` which is basic renewal model consisting of infections and reproduction numbers and `HospitalAdmissionsModel` which includes basic renewal model and hospital admissions. In the subsequent sections, we provide examples of fitting each of these models. diff --git a/docs/source/tutorials/hospital_admissions_model.md b/docs/source/tutorials/hospital_admissions_model.md deleted file mode 100644 index 3ff4573c..00000000 --- a/docs/source/tutorials/hospital_admissions_model.md +++ /dev/null @@ -1,5 +0,0 @@ -## Placeholder file - Please do not edit this file directly. - This file is just a placeholder. - For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/hospital_admissions_model.qmd diff --git a/docs/source/tutorials/hospital_admissions_model.qmd b/docs/source/tutorials/hospital_admissions_model.qmd deleted file mode 100644 index dc3aee8f..00000000 --- a/docs/source/tutorials/hospital_admissions_model.qmd +++ /dev/null @@ -1,591 +0,0 @@ ---- -title: Fitting a hospital admissions-only model -format: gfm -engine: jupyter ---- - -This document illustrates how a hospital admissions-only model can be fitted using data from the Pyrenew package, particularly the wastewater dataset. The CFA wastewater team created this dataset, which contains simulated data. - -We begin by loading `numpyro` and configuring the device count to 2 to enable running MCMC chains in parallel. By default, XLA (which is used by JAX for compilation) considers all CPU cores as one device. Depending on your system's configuration, we recommend using numpyro's [set_host_device_count()](https://num.pyro.ai/en/stable/utilities.html#set-host-device-count) function to set the number of devices available for parallel computing. - -```{python} -# | label: numpyro setup -import numpyro - -numpyro.set_host_device_count(2) -``` -## Model definition - -In this section, we provide the formal definition of the model. The hospital admissions model is a semi-mechanistic model that describes the number of observed hospital admissions as a function of a set of latent variables. Mainly, the observed number of hospital admissions is discretely distributed with location at the number of latent hospital admissions: - -$$ -h(t) \sim \text{HospDist}\left(H(t)\right) -$$ - -Where $h(t)$ is the observed number of hospital admissions at time $t$, and $H(t)$ is the number of latent hospital admissions at time $t$. The distribution $\text{HospDist}$ is discrete. For this example, we will use a negative binomial distribution with an inferred concentration. - -$$ -\begin{align*} -h(t) & \sim \mathrm{NegativeBinomial}\left(\mathrm{mean} = H(t), \mathrm{concentration} = k\right) \\ -H(t) & = p_\mathrm{hosp}(t) \sum_{\tau = 0}^{T_d} d(\tau) I(t-\tau) \\ -\log[p_\mathrm{hosp}(t)] & \sim \mathrm{Normal}(\mu=\log(0.05), \sigma=\log(1.1)) \\ -\log(k) & \sim \mathrm{Normal}(\mu=\log(1), \sigma=\log(10)) -\end{align*} -$$ - -Where $d(\tau)$ is the infection to hospital admission interval, $I(t)$ is the number of latent infections at time $t$, and $p_\mathrm{hosp}(t)$ is the infection to admission rate. - -The number of latent hospital admissions at time $t$ is a function of the number of latent infections at time $t$ and the infection to admission rate. The latent infections are modeled as a renewal process: - -$$ -\begin{align*} -I(t) &= \mathcal{R}(t) \times \sum_{\tau < t} I(\tau) g(t - \tau) \\ -\log[I(0)] &\sim \text{Normal}(\mu=\log(100), \sigma=\log(1.75)) -\end{align*} -$$ - -The reproduction number $\mathcal{R}(t)$ is modeled as a random walk in logarithmic space, i.e.: - -$$ -\begin{align*} -\log[\mathcal{R}(t)] & = \log[\mathcal{R}(t-1)] + \epsilon\\ -\epsilon & \sim \text{Normal}(\mu=0, \sigma=0.025) \\ -\mathcal{R}(0) &\sim \text{TruncatedNormal}(\text{loc}=1.2, \text{scale}=0.2, \text{min}=0) -\end{align*} -$$ - - -## Data processing - -We start by loading the data and inspecting the first five rows. - -```{python} -# | label: data-inspect -import polars as pl -from pyrenew import datasets - -dat = datasets.load_wastewater() -dat.head(5) -``` - -The data shows one entry per site, but the way it was simulated, the number of admissions is the same across sites. Thus, we will only keep the first observation per day. - -```{python} -# | label: aggregation -# Keeping the first observation of each date -dat = dat.group_by("date").first().select(["date", "daily_hosp_admits"]) - -# Now, sorting by date -dat = dat.sort("date") - -# Keeping the first 90 days -dat = dat.head(90) - -dat.head(5) -``` - -Let's take a look at the daily prevalence of hospital admissions. - -```{python} -# | label: fig-plot-hospital-admissions -# | fig-cap: Daily hospital admissions from the simulated data -import matplotlib.pyplot as plt -import matplotlib.dates as mdates - -daily_hosp_admits = dat["daily_hosp_admits"].to_numpy() -dates = dat["date"].to_numpy() -ax = plt.gca() -ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d")) -ax.xaxis.set_major_locator(mdates.DayLocator(interval=7)) -ax.set_xlim(dates[0], dates[-1]) -plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") -plt.plot(dates, daily_hosp_admits, "-o") -plt.xlabel("Date") -plt.ylabel("Admissions") -plt.show() -``` - -## Building the model - -First, we will extract two datasets we will use as deterministic quantities: the generation interval and the infection to hospital admission interval. - -```{python} -# | label: fig-data-extract -# | fig-cap: Generation interval and infection to hospital admission interval -gen_int = datasets.load_generation_interval() -inf_hosp_int = datasets.load_infection_admission_interval() - -# We only need the probability_mass column of each dataset -gen_int_array = gen_int["probability_mass"].to_numpy() -gen_int = gen_int_array -inf_hosp_int_array = inf_hosp_int["probability_mass"].to_numpy() - -# Taking a peek at the first 5 elements of each -gen_int[:5], inf_hosp_int_array[:5] - -# Visualizing both quantities side by side -fig, axs = plt.subplots(1, 2) - -axs[0].plot(gen_int) -axs[0].set_title("Generation interval") -axs[1].plot(inf_hosp_int_array) -axs[1].set_title("Infection to hospital admission interval") -plt.show() -``` - -With these two in hand, we can start building the model. First, we will define the latent hospital admissions: - -```{python} -# | label: latent-hosp -from pyrenew import latent, deterministic, metaclass, randomvariable -import jax.numpy as jnp -import numpyro.distributions as dist - -inf_hosp_int = deterministic.DeterministicPMF( - name="inf_hosp_int", value=inf_hosp_int_array -) - -hosp_rate = randomvariable.DistributionalVariable( - name="IHR", distribution=dist.LogNormal(jnp.log(0.05), jnp.log(1.1)) -) - -latent_hosp = latent.HospitalAdmissions( - infection_to_admission_interval_rv=inf_hosp_int, - infection_hospitalization_ratio_rv=hosp_rate, -) -``` - -The `inf_hosp_int` is a `DeterministicPMF` object that takes the infection to hospital admission interval as input. The `hosp_rate` is a `DistributionalVariable` object that takes a numpyro distribution to represent the infection to hospital admission rate. The `HospitalAdmissions` class is a `RandomVariable` that takes two distributions as inputs: the infection to admission interval and the infection to hospital admission rate. Now, we can define the rest of the other components: - -```{python} -# | label: initializing-rest-of-model -from pyrenew import model, process, observation, metaclass, transformation -from pyrenew.latent import ( - InfectionInitializationProcess, - InitializeInfectionsExponentialGrowth, -) - - -# Infection process -latent_inf = latent.Infections() -n_initialization_points = max(gen_int_array.size, inf_hosp_int_array.size) - 1 -I0 = InfectionInitializationProcess( - "I0_initialization", - randomvariable.DistributionalVariable( - name="I0", - distribution=dist.LogNormal(loc=jnp.log(100), scale=jnp.log(1.75)), - ), - InitializeInfectionsExponentialGrowth( - n_initialization_points, - deterministic.DeterministicVariable(name="rate", value=0.05), - ), -) - -# Generation interval and Rt -gen_int = deterministic.DeterministicPMF(name="gen_int", value=gen_int) - - -class MyRt(metaclass.RandomVariable): - def validate(self): - pass - - def sample(self, n: int, **kwargs) -> tuple: - sd_rt = numpyro.sample("Rt_random_walk_sd", dist.HalfNormal(0.025)) - - rt_rv = randomvariable.TransformedVariable( - name="log_rt_random_walk", - base_rv=process.RandomWalk( - name="log_rt", - step_rv=randomvariable.DistributionalVariable( - name="rw_step_rv", distribution=dist.Normal(0, 0.025) - ), - ), - transforms=transformation.ExpTransform(), - ) - rt_init_rv = randomvariable.DistributionalVariable( - name="init_log_rt", distribution=dist.Normal(0, 0.2) - ) - init_rt = rt_init_rv.sample() - - return rt_rv.sample(n=n, init_vals=init_rt, **kwargs) - - -rtproc = MyRt() - -# The observation model - -# we place a log-Normal prior on the concentration -# parameter of the negative binomial. -nb_conc_rv = randomvariable.TransformedVariable( - "concentration", - randomvariable.DistributionalVariable( - name="concentration_raw", - distribution=dist.TruncatedNormal(loc=0, scale=1, low=0.01), - ), - transformation.PowerTransform(-2), -) - -# now we define the observation process -obs = observation.NegativeBinomialObservation( - "negbinom_rv", - concentration_rv=nb_conc_rv, -) -``` - -Notice all the components are `RandomVariable` instances. We can now build the model: - -```{python} -# | label: init-model -hosp_model = model.HospitalAdmissionsModel( - latent_infections_rv=latent_inf, - latent_hosp_admissions_rv=latent_hosp, - I0_rv=I0, - gen_int_rv=gen_int, - Rt_process_rv=rtproc, - hosp_admission_obs_process_rv=obs, -) -``` - -Let's simulate from the prior predictive distribution to check that the model is working: - -```{python} -# | label: simulation -import numpy as np - -timeframe = 120 - -with numpyro.handlers.seed(rng_seed=223): - simulated_data = hosp_model.sample(n_datapoints=timeframe) -``` - -```{python} -# | label: fig-basic -# | fig-cap: Simulated Rt and Admissions -import matplotlib.pyplot as plt - -fig, axs = plt.subplots(1, 2) - -# Rt plot -axs[0].plot(simulated_data.Rt) -axs[0].set_ylabel("Simulated Rt") - -# Admissions plot -axs[1].plot(simulated_data.observed_hosp_admissions, "-o") -axs[1].set_ylabel("Simulated Admissions") - -fig.suptitle("Basic renewal model") -fig.supxlabel("Time") -plt.tight_layout() -plt.show() -``` - -## Fitting the model - -We now fit the model, not to these simulated data, but rather to the dataset we retrieved above. We use the `run` method of the `Model` object: - -```{python} -# | label: model-fit -import jax - -hosp_model.run( - num_samples=1000, - num_warmup=1000, - data_observed_hosp_admissions=daily_hosp_admits, - rng_key=jax.random.key(54), - mcmc_args=dict(progress_bar=False, num_chains=2), -) -``` - -We can use `arviz` to visualize the model fit. Here, we plot the observed values against the inferred latent values (i.e. the mean of the negative binomial observation process)[^capture]: - -[^capture]: The output is captured to avoid `quarto` from displaying the output twice. - -```{python} -# | label: fig-output-hospital-admissions -# | fig-cap: Latent hospital admissions posterior samples (gray) and observed admissions timeseries (red). -import arviz as az - -ppc_samples = hosp_model.posterior_predictive( - n_datapoints=daily_hosp_admits.size -) -idata = az.from_numpyro( - posterior=hosp_model.mcmc, - posterior_predictive=ppc_samples, -) - -axes = az.plot_ts( - idata, - y="negbinom_rv", - y_hat="negbinom_rv", - num_samples=200, - y_kwargs={ - "color": "blue", - "linewidth": 1.0, - "marker": "o", - "linestyle": "solid", - }, - y_hat_plot_kwargs={"color": "skyblue", "alpha": 0.05}, - y_mean_plot_kwargs={"color": "black", "linestyle": "--", "linewidth": 2.5}, - backend_kwargs={"figsize": (8, 6)}, - textsize=15.0, -) -ax = axes[0][0] -ax.set_xlabel("Time", fontsize=20) -ax.set_ylabel("Hospital Admissions", fontsize=20) -handles, labels = ax.get_legend_handles_labels() -ax.legend( - handles, ["Observed", "Sample Mean", "Posterior Samples"], loc="best" -) -plt.show() -``` - -## Results exploration and MCMC diagnostics - -To explore further, We can use [ArviZ](https://www.arviz.org/) to visualize the results. Let's start by loading the module and converting the fitted model to ArviZ `InferenceData` object: - -```{python} -# | label: convert-inferenceData - -idata = az.from_numpyro(hosp_model.mcmc) -``` - -We obtain the summary of model diagnostics and print the diagnostics for `latent_hospital_admissions[1]` - -```{python} -# | label: diagnostics -# | warning: false -diagnostic_stats_summary = az.summary( - idata.posterior, - kind="diagnostics", -) - -print(diagnostic_stats_summary.loc["latent_hospital_admissions[1]"]) -``` - -Below we plot 90% and 50% highest density intervals for latent hospital admissions using [plot_hdi](https://python.arviz.org/en/stable/api/generated/arviz.plot_hdi.html): - -```{python} -# | label: fig-output-admission-distribution -# | fig-cap: Hospital Admissions posterior distribution -x_data = idata.posterior["latent_hospital_admissions_dim_0"] -y_data = idata.posterior["latent_hospital_admissions"] - - -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.9, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.5, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -# Add the posterior median to the figure -median_ts = y_data.median(dim=["chain", "draw"]) - -axes.plot(x_data, median_ts, color="C0", label="Median") -axes.legend() -axes.set_title("Posterior Hospital Admissions", fontsize=10) -axes.set_xlabel("Time", fontsize=10) -axes.set_ylabel("Hospital Admissions", fontsize=10) -plt.show() -``` - -We can also look at credible intervals for the posterior distribution of latent infections: - -```{python} -# | label: fig-output-infections-distribution -# | fig-cap: Posterior Latent Infections -x_data = ( - idata.posterior["all_latent_infections_dim_0"] - n_initialization_points -) -y_data = idata.posterior["all_latent_infections"] - -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.9, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - x_data, - y_data, - hdi_prob=0.5, - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -# Add the posterior median to the figure -median_ts = y_data.median(dim=["chain", "draw"]) -axes.plot(x_data, median_ts, color="C0", label="Median") -axes.legend() -``` - - -## Predictive checks and forecasting -We can use the `Model`'s `posterior_predictive` and `prior_predictive` methods to generate posterior and prior predictive samples for observed admissions. - -```{python} -# | label: demonstrate-use-of-predictive-methods - -idata = az.from_numpyro( - hosp_model.mcmc, - posterior_predictive=hosp_model.posterior_predictive( - n_datapoints=len(daily_hosp_admits) - ), - prior=hosp_model.prior_predictive( - n_datapoints=len(daily_hosp_admits), - numpyro_predictive_args={"num_samples": 1000}, - ), -) -``` - -We will use `plot_lm` method from ArviZ to plot the posterior predictive distribution against the actual observed data below: - -```{python} -# | label: fig-posterior-predictive -# | fig-cap: Hospital Admissions posterior distribution with plot_lm -fig, ax = plt.subplots() -az.plot_lm( - "negbinom_rv", - idata=idata, - kind_pp="hdi", - y_kwargs={"color": "black"}, - y_hat_fill_kwargs={"color": "C0"}, - axes=ax, -) - -ax.set_title("Posterior Predictive Plot") -ax.set_ylabel("Hospital Admissions") -ax.set_xlabel("Days") -plt.show() -``` - - -By increasing `n_datapoints`, we can perform forecasting using the posterior predictive distribution. - -```{python} -# | label: posterior-predictive-distribution -n_forecast_points = 28 -idata = az.from_numpyro( - hosp_model.mcmc, - posterior_predictive=hosp_model.posterior_predictive( - n_datapoints=len(daily_hosp_admits) + n_forecast_points, - ), - prior=hosp_model.prior_predictive( - n_datapoints=len(daily_hosp_admits), - numpyro_predictive_args={"num_samples": 1000}, - ), -) -``` - -Below we plot the prior predictive distributions using equal tailed Bayesian credible intervals: - -```{python} -# | label: fig-output-prior-predictive -# | fig-cap: Prior Predictive Admissions -def compute_eti(dataset, eti_prob): - eti_bdry = dataset.quantile( - ((1 - eti_prob) / 2, 1 / 2 + eti_prob / 2), dim=("chain", "draw") - ) - return eti_bdry.T - - -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - idata.prior_predictive["negbinom_rv_dim_0"], - hdi_data=compute_eti(idata.prior_predictive["negbinom_rv"], 0.9), - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - idata.prior_predictive["negbinom_rv_dim_0"], - hdi_data=compute_eti(idata.prior_predictive["negbinom_rv"], 0.5), - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -plt.scatter( - idata.observed_data["negbinom_rv_dim_0"], - idata.observed_data["negbinom_rv"], - color="black", -) - -axes.set_title("Prior Predictive Admissions", fontsize=10) -axes.set_xlabel("Time", fontsize=10) -axes.set_ylabel("Observed Admissions", fontsize=10) -plt.yscale("log") -plt.show() -``` - -And now we plot the posterior predictive distributions with a `{python} n_forecast_points`-day-ahead forecast: -```{python} -# | label: fig-output-posterior-predictive-forecast -# | fig-cap: Posterior predictive admissions, including a forecast. -x_data = idata.posterior_predictive["negbinom_rv_dim_0"] -y_data = idata.posterior_predictive["negbinom_rv"] -fig, axes = plt.subplots(figsize=(6, 5)) -az.plot_hdi( - x_data, - hdi_data=compute_eti(y_data, 0.9), - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.3}, - ax=axes, -) - -az.plot_hdi( - x_data, - hdi_data=compute_eti(y_data, 0.5), - color="C0", - smooth=False, - fill_kwargs={"alpha": 0.6}, - ax=axes, -) - -# Add median of the posterior to the figure -median_ts = y_data.median(dim=["chain", "draw"]) - -plt.plot( - x_data, - median_ts, - color="C0", - label="Median", -) -plt.scatter( - idata.observed_data["negbinom_rv_dim_0"], - idata.observed_data["negbinom_rv"], - color="black", -) -axes.legend() -axes.set_title( - "Posterior Predictive Admissions, including a forecast", fontsize=10 -) -axes.set_xlabel("Time", fontsize=10) -axes.set_ylabel("Hospital Admissions", fontsize=10) -plt.show() -``` diff --git a/docs/source/tutorials/index.md b/docs/source/tutorials/index.md deleted file mode 100644 index d226edc5..00000000 --- a/docs/source/tutorials/index.md +++ /dev/null @@ -1,14 +0,0 @@ -# Tutorials - -This section contains tutorials that demonstrate how to use the `pyrenew` package. The source code for the tutorials can be found in the project repository: . - -```{toctree} -:maxdepth: 2 - -getting_started -basic_renewal_model -hospital_admissions_model -extending_pyrenew -periodic_effects -day_of_the_week -``` diff --git a/docs/source/tutorials/periodic_effects.qmd b/docs/source/tutorials/periodic_effects.qmd deleted file mode 100644 index 7f3cf5b7..00000000 --- a/docs/source/tutorials/periodic_effects.qmd +++ /dev/null @@ -1,105 +0,0 @@ ---- -title: Periodic effects and broadcasting -format: gfm -engine: jupyter ---- - -A common practice in time series forecasting is featuring periodic effects such as day of the week or weekend effects. `pyrenew` has built-in functions and classes to incorporate such effects into renewal models. Notably, the `arrayutils.PeriodicBroadcaster` class. - -This document demonstrates how to use these functions and classes to model periodic effects and broadcasting in renewal models: Repeating elements of a sequence and repeating a sequence as a whole (tiling). - -## Repeated units - -The `RtPeriodicDiff` and `RtWeeklyDiff` classes use `PeriodicBroadcaster` to repeat each vector element to represent Rt values held constant within a period. Let's show it's usage with a basic renewal model. We will use the same model structure shown in the `getting_started` tutorial: - -```{python} -# | label: loading-pkgs -# | output: false -# | warning: false -import jax.numpy as jnp -import numpy as np -import numpyro -from pyrenew import process, deterministic -``` - -```{python} -# The random process for Rt -rt_proc = process.RtWeeklyDiffARProcess( - name="rt_weekly_diff", - offset=0, - log_rt_rv=deterministic.DeterministicVariable( - name="log_rt", value=jnp.array([0.1, 0.2]) - ), - autoreg_rv=deterministic.DeterministicVariable( - name="autoreg", value=jnp.array([0.7]) - ), - periodic_diff_sd_rv=deterministic.DeterministicVariable( - name="periodic_diff_sd", value=jnp.array([0.1]) - ), -) -``` - -```{python} -with numpyro.handlers.seed(rng_seed=20): - sim_data = rt_proc(duration=30) - -# Plotting the Rt values -import matplotlib.pyplot as plt - -plt.step(np.arange(len(sim_data)), sim_data, where="post") -plt.xlabel("Time") -plt.ylabel("Rt") -plt.title("Simulated Rt values") - -# Adding bands to mark weeks -for i in range(0, 30, 7): - plt.axvline(i, color="black", linestyle="--", alpha=0.5) -plt.show() -``` - -The implementation of the `RtWeeklyDiffARProcess` (which is an instance of `RtPeriodicDiffARProcess`), uses `repeat_until_n` to repeating values: `repeat_until_n(..., period_size=7)`. The `RtWeeklyDiff` class is a particular case of `RtPeriodicDiff` with a period size of seven. - -## Repeated sequences (tiling) - -The `PeriodicBroadcaster` class can also be used to repeat a sequence as a whole. For example, we can use this to model a periodic effect such as day-of-the-week effect. Let's use the `DayOfWeekEffect` random variable. We will sample from a scaled Dirchlet distribution such that the sum of the samples is 7: - -```{python} -import numpyro.distributions as dist -from pyrenew import transformation, randomvariable - -# Building the transformed prior: Dirichlet * 7 -mysimplex = dist.TransformedDistribution( - dist.Dirichlet(concentration=jnp.ones(7)), - transformation.AffineTransform(loc=0, scale=7.0), -) - -# Constructing the day of week effect -dayofweek = process.DayOfWeekEffect( - offset=0, - quantity_to_broadcast=randomvariable.DistributionalVariable( - name="simp", distribution=mysimplex - ), -) -``` - -Like before, we can use the `sample` method to generate samples from the day of week effect: - -```{python} -with numpyro.handlers.seed(rng_seed=20): - sim_data = dayofweek(duration=30) - -# Plotting the effect values -import matplotlib.pyplot as plt - -plt.step(np.arange(len(sim_data)), sim_data, where="post") -plt.xlabel("Time") -plt.ylabel("Effect size") -plt.title("Simulated Day of Week Effect values") - -# Adding bands to mark weeks -for i in range(0, 30, 7): - plt.axvline(i, color="black", linestyle="--", alpha=0.5) -plt.show() -``` - -Like the example from the previous section, the `DayOfWeekEffect` class uses `PeriodicBroadcaster` to broadcast the data, this time using `broadcast_type = "tile"`, repeating the sequence as a whole. The `DayOfWeekEffect` class is a particular case of `PeriodicEffect` with a period size of seven. diff --git a/docs/tutorials/basic_renewal_model.md b/docs/tutorials/basic_renewal_model.md new file mode 100644 index 00000000..34e51257 --- /dev/null +++ b/docs/tutorials/basic_renewal_model.md @@ -0,0 +1,5 @@ +## Placeholder file + Please do not edit this file directly. + This file is just a placeholder. + For the source file, see: + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/basic_renewal_model.qmd diff --git a/docs/source/tutorials/day_of_the_week.md b/docs/tutorials/day_of_the_week.md similarity index 58% rename from docs/source/tutorials/day_of_the_week.md rename to docs/tutorials/day_of_the_week.md index f00449a0..a1e77635 100644 --- a/docs/source/tutorials/day_of_the_week.md +++ b/docs/tutorials/day_of_the_week.md @@ -2,4 +2,4 @@ Please do not edit this file directly. This file is just a placeholder. For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/day_of_the_week.qmd + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/day_of_the_week.qmd diff --git a/docs/source/tutorials/extending_pyrenew.md b/docs/tutorials/extending_pyrenew.md similarity index 57% rename from docs/source/tutorials/extending_pyrenew.md rename to docs/tutorials/extending_pyrenew.md index 7da1a07a..8cdf044c 100644 --- a/docs/source/tutorials/extending_pyrenew.md +++ b/docs/tutorials/extending_pyrenew.md @@ -2,4 +2,4 @@ Please do not edit this file directly. This file is just a placeholder. For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/extending_pyrenew.qmd + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/extending_pyrenew.qmd diff --git a/docs/source/tutorials/getting_started.md b/docs/tutorials/getting_started.md similarity index 58% rename from docs/source/tutorials/getting_started.md rename to docs/tutorials/getting_started.md index 8e2d77b7..4f662a22 100644 --- a/docs/source/tutorials/getting_started.md +++ b/docs/tutorials/getting_started.md @@ -2,4 +2,4 @@ Please do not edit this file directly. This file is just a placeholder. For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/getting_started.qmd + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/getting_started.qmd diff --git a/docs/tutorials/hospital_admissions_model.md b/docs/tutorials/hospital_admissions_model.md new file mode 100644 index 00000000..1234d982 --- /dev/null +++ b/docs/tutorials/hospital_admissions_model.md @@ -0,0 +1,5 @@ +## Placeholder file + Please do not edit this file directly. + This file is just a placeholder. + For the source file, see: + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/hospital_admissions_model.qmd diff --git a/docs/source/tutorials/periodic_effects.md b/docs/tutorials/periodic_effects.md similarity index 57% rename from docs/source/tutorials/periodic_effects.md rename to docs/tutorials/periodic_effects.md index a649b468..abfac76b 100644 --- a/docs/source/tutorials/periodic_effects.md +++ b/docs/tutorials/periodic_effects.md @@ -2,4 +2,4 @@ Please do not edit this file directly. This file is just a placeholder. For the source file, see: - https://github.com/CDCgov/PyRenew/tree/main/docs/source/tutorials/periodic_effects.qmd + https://github.com/CDCgov/PyRenew/tree/main/docs/tutorials/periodic_effects.qmd From 65f2d7e24fe913add620b318ae0a801722f7bf8a Mon Sep 17 00:00:00 2001 From: Damon Bayer Date: Tue, 9 Sep 2025 15:05:33 -0500 Subject: [PATCH 4/8] pre-commit --- docs/developer_documentation.md | 2 +- docs/javascripts/mathjax.js | 2 +- docs/reference/index.md | 2 +- docs/stylesheets/extra.css | 2 +- mkdocs.yml | 37 ++++++++++++++++----------------- 5 files changed, 22 insertions(+), 23 deletions(-) diff --git a/docs/developer_documentation.md b/docs/developer_documentation.md index a6b4524a..36da15db 100644 --- a/docs/developer_documentation.md +++ b/docs/developer_documentation.md @@ -102,7 +102,7 @@ nav: # ... other reference pages - Tutorials: - tutorials/index.md - # ... other tutorial pages + # ... other tutorial pages - Developer Documentation: developer_documentation.md ```` diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js index 09fbd880..06dbf38b 100644 --- a/docs/javascripts/mathjax.js +++ b/docs/javascripts/mathjax.js @@ -13,4 +13,4 @@ window.MathJax = { document$.subscribe(() => { MathJax.typesetPromise() -}) \ No newline at end of file +}) diff --git a/docs/reference/index.md b/docs/reference/index.md index 0a1a9e01..edde490d 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -5,7 +5,7 @@ This section provides detailed documentation for all classes, functions, and mod ## Modules - [Model](model.md) - Core renewal model implementations -- [Latent](latent.md) - Latent variable components +- [Latent](latent.md) - Latent variable components - [Process](process.md) - Stochastic process implementations - [Random Variable](randomvariable.md) - Random variable abstractions - [Observation](observation.md) - Observation model components diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css index 4b51ffb6..e5ccfca2 100644 --- a/docs/stylesheets/extra.css +++ b/docs/stylesheets/extra.css @@ -32,4 +32,4 @@ table th { /* MathJax styling */ .MathJax { font-size: 1em !important; -} \ No newline at end of file +} diff --git a/mkdocs.yml b/mkdocs.yml index beb9907f..6f7b5400 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -8,24 +8,24 @@ edit_uri: "edit/main/docs/" nav: - Home: index.md - API Reference: - - reference/index.md - - Model: reference/model.md - - Latent: reference/latent.md - - Process: reference/process.md - - Random Variable: reference/randomvariable.md - - Observation: reference/observation.md - - Distributions: reference/distributions.md - - Deterministic: reference/deterministic.md - - Datasets: reference/datasets.md - - Additional: reference/additional.md + - reference/index.md + - Model: reference/model.md + - Latent: reference/latent.md + - Process: reference/process.md + - Random Variable: reference/randomvariable.md + - Observation: reference/observation.md + - Distributions: reference/distributions.md + - Deterministic: reference/deterministic.md + - Datasets: reference/datasets.md + - Additional: reference/additional.md - Tutorials: - - tutorials/index.md - - Getting Started: tutorials/getting_started.md - - Basic Renewal Model: tutorials/basic_renewal_model.md - - Day of the Week Effects: tutorials/day_of_the_week.md - - Extending PyRenew: tutorials/extending_pyrenew.md - - Hospital Admissions Model: tutorials/hospital_admissions_model.md - - Periodic Effects: tutorials/periodic_effects.md + - tutorials/index.md + - Getting Started: tutorials/getting_started.md + - Basic Renewal Model: tutorials/basic_renewal_model.md + - Day of the Week Effects: tutorials/day_of_the_week.md + - Extending PyRenew: tutorials/extending_pyrenew.md + - Hospital Admissions Model: tutorials/hospital_admissions_model.md + - Periodic Effects: tutorials/periodic_effects.md - Developer Documentation: developer_documentation.md theme: @@ -82,7 +82,6 @@ markdown_extensions: custom_fences: - name: mermaid class: mermaid - format: !!python/name:pymdownx.superfences.fence_code_format - pymdownx.highlight: anchor_linenums: true line_spans: __span @@ -107,4 +106,4 @@ extra_css: - stylesheets/extra.css watch: - - pyrenew/ \ No newline at end of file + - pyrenew/ From ea83f14a68fe35ed03fac134718602561d7dcd81 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 9 Sep 2025 20:32:11 +0000 Subject: [PATCH 5/8] Fix mkdocs.yml path issue in GitHub Actions workflow Co-authored-by: damonbayer <8095080+damonbayer@users.noreply.github.com> --- .github/workflows/website.yaml | 8 ++++---- docs/Makefile | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/website.yaml b/.github/workflows/website.yaml index 186fbd29..a19cece0 100644 --- a/.github/workflows/website.yaml +++ b/.github/workflows/website.yaml @@ -46,9 +46,9 @@ jobs: # Building - name: Build HTML using uv run: | - cd docs - uv run make tutorials_clean - uv run make html + cp README.md docs/index.md + cd docs && uv run make tutorials_clean + uv run mkdocs build # Upload - name: Upload artifacts @@ -56,7 +56,7 @@ jobs: with: name: github-pages retention-days: 7 - path: docs/site/ + path: site/ deploy: # Deploy to the github-pages environment diff --git a/docs/Makefile b/docs/Makefile index 5d67e470..20be0e46 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -68,8 +68,8 @@ tutorials_clean: .PHONY: html html: $(TUTORIALS_MD) - @$(MKDOCSBUILD) build + @cd .. && $(MKDOCSBUILD) build .PHONY: serve serve: $(TUTORIALS_MD) - @$(MKDOCSBUILD) serve + @cd .. && $(MKDOCSBUILD) serve From 9dbbede397c153ac3408d6461f6a516f1ff06a3e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 9 Sep 2025 21:10:06 +0000 Subject: [PATCH 6/8] Fix GitHub Actions workflow to run mkdocs build from root directory Co-authored-by: damonbayer <8095080+damonbayer@users.noreply.github.com> --- .github/workflows/website.yaml | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/.github/workflows/website.yaml b/.github/workflows/website.yaml index a19cece0..1388ef00 100644 --- a/.github/workflows/website.yaml +++ b/.github/workflows/website.yaml @@ -39,16 +39,12 @@ jobs: - name: Install the project run: uv sync --all-extras - - name: Render documentation index from README - run: | - cp README.md docs/index.md - # Building - name: Build HTML using uv run: | cp README.md docs/index.md cd docs && uv run make tutorials_clean - uv run mkdocs build + cd .. && uv run mkdocs build # Upload - name: Upload artifacts From e57f9b4a8cf47997ba22e0f394215cb53181bae0 Mon Sep 17 00:00:00 2001 From: Damon Bayer Date: Wed, 10 Sep 2025 08:10:48 -0500 Subject: [PATCH 7/8] offline --- mkdocs.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/mkdocs.yml b/mkdocs.yml index 6f7b5400..193f1765 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -59,6 +59,7 @@ theme: plugins: - search + - offline - mkdocstrings: handlers: python: From 2085321d941367be63b6e6e18a5f8ed48853064d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 10 Sep 2025 13:35:43 +0000 Subject: [PATCH 8/8] Fix tutorial generation in GitHub Actions workflow Co-authored-by: damonbayer <8095080+damonbayer@users.noreply.github.com> --- .github/workflows/website.yaml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/website.yaml b/.github/workflows/website.yaml index 1388ef00..38e6d21e 100644 --- a/.github/workflows/website.yaml +++ b/.github/workflows/website.yaml @@ -43,8 +43,7 @@ jobs: - name: Build HTML using uv run: | cp README.md docs/index.md - cd docs && uv run make tutorials_clean - cd .. && uv run mkdocs build + cd docs && uv run make html # Upload - name: Upload artifacts