diff --git a/.github/workflows/test_model.yml b/.github/workflows/test_model.yml index 9ea0c79f..395ced2e 100644 --- a/.github/workflows/test_model.yml +++ b/.github/workflows/test_model.yml @@ -29,6 +29,7 @@ jobs: - name: run tests run: | poetry run -C model pytest \ + --mpl --mpl-default-tolerance=10 \ --cov=pyrenew --cov-report term --cov-report xml model - name: Upload coverage to Codecov diff --git a/model/Makefile b/model/Makefile index a6560b69..89dd88b6 100644 --- a/model/Makefile +++ b/model/Makefile @@ -2,9 +2,10 @@ install: poetry install test: - poetry run pytest + poetry run pytest --mpl --mpl-default-tolerance=10 -docs: docs/pyrenew_demo.md docs/getting-started.md +docs: docs/pyrenew_demo.md docs/getting-started.md \ + docs/example-with-datasets.md docs/pyrenew_demo.md: docs/pyrenew_demo.qmd poetry run quarto render docs/pyrenew_demo.qmd @@ -12,9 +13,30 @@ docs/pyrenew_demo.md: docs/pyrenew_demo.qmd docs/getting-started.md: docs/getting-started.qmd poetry run quarto render docs/getting-started.qmd +docs/example-with-datasets.md: docs/example-with-datasets.qmd + poetry run quarto render docs/example-with-datasets.qmd + +docs/py: docs/notebooks + jupyter nbconvert --to python docs/pyrenew_demo.ipynb + jupyter nbconvert --to python docs/getting-started.ipynb + jupyter nbconvert --to python docs/example-with-datasets.ipynb + +docs/notebooks: + quarto convert docs/pyrenew_demo.qmd --output docs/pyrenew_demo.ipynb + quarto convert docs/getting-started.qmd --output docs/getting-started.ipynb + quarto convert docs/example-with-datasets.qmd --output \ + docs/example-with-datasets.ipynb + +test_images: + echo "Generating reference images for tests" + poetry run pytest --mpl-generate-path=src/test/baseline + clean: rm -rf docs/*_files/ - rm -f docs/getting-started.ipynb - rm -f docs/pyrenew_demo.ipynb + for i in .py, .md, .ipynb; do \ + rm -f docs/getting-started.$$i; \ + rm -f docs/pyrenew_demo.$$i; \ + rm -f docs/example-with-datasets.$$i; \ + done -.PHONY: install test docs clean +.PHONY: install test docs clean docs/notebooks docs/py test_images diff --git a/model/docs/example-with-datasets.md b/model/docs/example-with-datasets.md new file mode 100644 index 00000000..9bee20aa --- /dev/null +++ b/model/docs/example-with-datasets.md @@ -0,0 +1,472 @@ +# Fitting a Hospital Admissions-only Model + + +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. + +## Model definition + +In this section, we provide the formal definition of the model. The +hospitalization 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: + +$$ +\begin{align*} +h(t) & \sim \text{NegativeNinomial}\left(\text{concentration} = 1, \text{mean} = H(t)\right) \\ +H(t) & = \omega(t) p_\mathrm{hosp}(t) \sum_{\tau = 0}^{T_d} d(\tau) I(t-\tau) +\end{align*} +$$ + +Were $d(\tau)$ is the infection to hospitalization interval, $I(t)$ is +the number of latent infections at time $t$, $p_\mathrm{hosp}(t)$ is the +infection to hospitalization rate, and $\omega(t)$ is the weekday effect +at time $t$; the last section provides an example building a weekly +effect `RandomVariable`. + +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 +hospitalization rate. The latent infections are modeled as a renewal +process: + +$$ +\begin{align*} +I(t) &= R(t) \times \sum_{\tau < t} I(\tau) g(t - \tau) \\ +I(0) &\sim \text{LogNormal}(\mu = \log(80/0.05), \sigma = 1.5) +\end{align*} +$$ + +The reproductive number $R(t)$ is modeled as a random walk process: + +$$ +\begin{align*} +R(t) & = R(t-1) + \epsilon\\ +\log{\epsilon} & \sim \text{Normal}(\mu=0, \sigma=0.1) \\ +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 +import polars as pl +from pyrenew import datasets + +dat = datasets.load_wastewater() +dat.head(5) +``` + +
+shape: (5, 14) + +| t | lab_wwtp_unique_id | log_conc | date | lod_sewage | below_LOD | daily_hosp_admits | daily_hosp_admits_for_eval | pop | forecast_date | hosp_calibration_time | site | ww_pop | inf_per_capita | +|-----|--------------------|----------|--------------|------------|-----------|-------------------|----------------------------|-----|---------------|-----------------------|------|----------|----------------| +| i64 | i64 | f64 | str | f64 | i64 | i64 | i64 | f64 | str | i64 | i64 | f64 | f64 | +| 1 | 1 | null | "2023-10-30" | null | null | 6 | 6 | 1e6 | "2024-02-05" | 90 | 1 | 400000.0 | 0.000663 | +| 1 | 2 | null | "2023-10-30" | null | null | 6 | 6 | 1e6 | "2024-02-05" | 90 | 1 | 400000.0 | 0.000663 | +| 1 | 3 | null | "2023-10-30" | null | null | 6 | 6 | 1e6 | "2024-02-05" | 90 | 2 | 200000.0 | 0.000663 | +| 1 | 4 | null | "2023-10-30" | null | null | 6 | 6 | 1e6 | "2024-02-05" | 90 | 3 | 100000.0 | 0.000663 | +| 1 | 5 | null | "2023-10-30" | null | null | 6 | 6 | 1e6 | "2024-02-05" | 90 | 4 | 50000.0 | 0.000663 | + +
+ +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 +# 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) +``` + +
+shape: (5, 2) + +| date | daily_hosp_admits | +|--------------|-------------------| +| str | i64 | +| "2023-10-30" | 6 | +| "2023-10-31" | 8 | +| "2023-11-01" | 4 | +| "2023-11-02" | 8 | +| "2023-11-03" | 4 | + +
+ +Let’s take a look at the daily prevalence of hospital admissions. + +``` python +import matplotlib.pyplot as plt + +# Rotating the x-axis labels, and only showing ~10 labels +ax = plt.gca() +ax.xaxis.set_major_locator(plt.MaxNLocator(nbins=10)) +ax.xaxis.set_tick_params(rotation=45) +plt.plot(dat["date"].to_numpy(), dat["daily_hosp_admits"].to_numpy()) +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 hospitalization +interval. + +``` python +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 = gen_int["probability_mass"].to_numpy() +inf_hosp_int = inf_hosp_int["probability_mass"].to_numpy() + +# Taking a pick at the first 5 elements of each +gen_int[:5], inf_hosp_int[: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) +axs[1].set_title("Infection to hospitalization interval") +``` + + Text(0.5, 1.0, 'Infection to hospitalization interval') + + + +With these two in hand, we can start building the model. First, we will +define the latent hospital admissions: + +``` python +from pyrenew import latent, deterministic +import jax.numpy as jnp +import numpyro.distributions as dist + +inf_hosp_int = deterministic.DeterministicPMF(inf_hosp_int) + +hosp_rate = latent.InfectHospRate( + dist=dist.LogNormal(jnp.log(0.05), 0.1) +) + +latent_hosp = latent.HospitalAdmissions( + infection_to_admission_interval=inf_hosp_int, + infect_hosp_rate_dist=hosp_rate, + ) +``` + + /mnt/c/Users/xrd4/Documents/repos/msr/model/.venv/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html + from .autonotebook import tqdm as notebook_tqdm + An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu. + +The `inf_hosp_int` is a `DeterministicPMF` object that takes the +infection to hospitalization interval as input. The `hosp_rate` is an +`InfectHospRate` object that takes the infection to hospitalization rate +as input. The `HospitalAdmissions` class is a `RandomVariable` that +takes two distributions as inputs: the infection to admission interval +and the infection to hospitalization rate. Now, we can define the rest +of the other components: + +``` python +from pyrenew import model, process, observation + +# Infection process +latent_inf = latent.Infections() +I0 = latent.Infections0(I0_dist=dist.LogNormal(loc=jnp.log(80/.05), scale=1.5)) + +# Generation interval and Rt +gen_int = deterministic.DeterministicPMF(gen_int) +rtproc = process.RtRandomWalkProcess( + Rt_rw_dist=dist.Normal(0, 0.1) +) + +# The observation model +obs = observation.NegativeBinomialObservation(concentration_prior=1.0) +``` + +Notice all the components are `RandomVariable` instances. We can now +build the model: + +``` python +hosp_model = model.HospitalizationsModel( + latent_infections=latent_inf, + latent_hospitalizations=latent_hosp, + I0=I0, + gen_int=gen_int, + Rt_process=rtproc, + observation_process=obs, +) +``` + +Let’s simulate to check if the model is working: + +``` python +import numpyro as npro +import numpy as np + +timeframe = 120 + +np.random.seed(223) +with npro.handlers.seed(rng_seed = np.random.randint(1, timeframe)): + sim_data = hosp_model.sample(n_timepoints=timeframe) +``` + +``` python +import matplotlib.pyplot as plt + +fig, axs = plt.subplots(1, 2) + +# Rt plot +axs[0].plot(range(0, timeframe + 1), sim_data.Rt) +axs[0].set_ylabel('Rt') + +# Infections plot +axs[1].plot(range(0, timeframe + 1), sim_data.sampled_admissions) +axs[1].set_ylabel('Infections') +axs[1].set_yscale('log') + +fig.suptitle('Basic renewal model') +fig.supxlabel('Time') +plt.tight_layout() +plt.show() +``` + +![Rt and +Infections](example-with-datasets_files/figure-commonmark/basic-fig-output-1.png) + +## Fitting the model + +We can fit the model to the data. We will use the `run` method of the +model object. The two inputs this model requires are `n_timepoints` and +`observed_hospitalizations` + +``` python +import jax + +hosp_model.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat.shape[0] - 1, + observed_hospitalizations=dat["daily_hosp_admits"].to_numpy(), + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), +) +``` + +We can use the `plot_posterior` method to visualize the results[^1]: + +``` python +out = hosp_model.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat["daily_hosp_admits"].to_numpy(), +) +``` + +![Hospital Admissions posterior +distribution](example-with-datasets_files/figure-commonmark/output-hospital-admissions-output-1.png) + +The first half of the model is not looking good. The reason is that the +infection to hospitalization interval PMF makes it unlikely to observe +admissions from the beginning. The following section shows how to fix +this. + +## Padding the model + +We can use the padding argument to solve the overestimation of +hospitalizations in the first half of the model. By setting +`padding > 0`, the model then assumes that the first `padding` +observations are missing; thus, only observations after `padding` will +count towards the likelihood of the model. In practice, the model will +extend the estimated Rt and latent infections by `padding` days, given +time to adjust to the observed data. The following code will add 21 days +of missing data at the beginning of the model and re-estimate it with +`padding = 21`: + +``` python +days_to_impute = 21 + +dat2 = dat["daily_hosp_admits"].to_numpy() + +# Add 21 Nas to the beginning of dat2 +dat2 = np.hstack((np.repeat(np.nan, days_to_impute), dat2)) + +hosp_model.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat2.shape[0] - 1, + observed_hospitalizations=dat2, + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), + padding=days_to_impute, # Padding the model +) +``` + +And plotting the results: + +``` python +out = hosp_model.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat2, +) +``` + +![Hospital Admissions posterior +distribution](example-with-datasets_files/figure-commonmark/output-admissions-with-padding-output-1.png) + +We can also take a look at the latent infections: + +``` python +out2 = hosp_model.plot_posterior( + var="latent_infections", + ylab="Latent Infections" +) +``` + +![Hospital Admissions posterior +distribution](example-with-datasets_files/figure-commonmark/output-infections-with-padding-output-1.png) + +## Round 2: Incorporating weekday 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 +weekday effect distribution. To do this, we will create a new instance +of `RandomVariable` to model the weekday effect. The weekday effect will +be a truncated normal distribution with a mean of 1.0 and a standard +deviation of 0.5. The distribution will be truncated between 0.1 and +10.0. The weekday effect will be repeated for the number of weeks in the +dataset. + +``` python +from pyrenew import metaclass +import numpyro as npro + +class WeekdayEffect(metaclass.RandomVariable): + """Weekday effect distribution""" + def __init__(self, len: int): + """ Initialize the weekday effect distribution + Parameters + ---------- + len : int + The number of observations + """ + self.nweeks = jnp.ceil(len/7).astype(int) + self.len = len + + @staticmethod + def validate(): + return None + + def sample(self, **kwargs): + ans = npro.sample( + name="weekday_effect", + fn=npro.distributions.TruncatedNormal( + loc=1.0, scale=.5, low=0.1, high=10.0 + ), + sample_shape=(7,) + ) + + return jnp.tile(ans, self.nweeks)[:self.len] + +# Initializing the weekday effect +weekday_effect = WeekdayEffect(dat.shape[0]) +``` + +Notice that the instance’s `nweeks` and `len` members are passed during +construction. Trying to compute the number of weeks and the length of +the dataset in the `validate` method will raise a `jit` error in `jax` +as the shape and size of elements are not known during the validation +step, which happens before the model is run. With the new weekday +effect, we can rebuild the latent hospitalization model: + +``` python +latent_hosp_wday_effect = latent.HospitalAdmissions( + infection_to_admission_interval=inf_hosp_int, + infect_hosp_rate_dist=hosp_rate, + weekday_effect_dist=weekday_effect, + ) + +hosp_model_weekday = model.HospitalizationsModel( + latent_infections=latent_inf, + latent_hospitalizations=latent_hosp_wday_effect, + I0=I0, + gen_int=gen_int, + Rt_process=rtproc, + observation_process=obs, +) +``` + +Running the model (with the same padding as before): + +``` python +hosp_model_weekday.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat2.shape[0] - 1, + observed_hospitalizations=dat2, + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), + padding=days_to_impute, +) +``` + +And plotting the results: + +``` python +out = hosp_model_weekday.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat2, +) +``` + +![Hospital Admissions posterior +distribution](example-with-datasets_files/figure-commonmark/output-admissions-padding-and-weekday-output-1.png) + +[^1]: The output is captured to avoid `quarto` from displaying the + output twice. diff --git a/model/docs/example-with-datasets.qmd b/model/docs/example-with-datasets.qmd new file mode 100644 index 00000000..6c58a08c --- /dev/null +++ b/model/docs/example-with-datasets.qmd @@ -0,0 +1,381 @@ +--- +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. + +## Model definition + +In this section, we provide the formal definition of the model. The hospitalization 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: + +$$ +\begin{align*} +h(t) & \sim \text{NegativeNinomial}\left(\text{concentration} = 1, \text{mean} = H(t)\right) \\ +H(t) & = \omega(t) p_\mathrm{hosp}(t) \sum_{\tau = 0}^{T_d} d(\tau) I(t-\tau) +\end{align*} +$$ + +Were $d(\tau)$ is the infection to hospitalization interval, $I(t)$ is the number of latent infections at time $t$, $p_\mathrm{hosp}(t)$ is the infection to hospitalization rate, and $\omega(t)$ is the weekday effect at time $t$; the last section provides an example building a weekly effect `RandomVariable`. + +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 hospitalization rate. The latent infections are modeled as a renewal process: + +$$ +\begin{align*} +I(t) &= R(t) \times \sum_{\tau < t} I(\tau) g(t - \tau) \\ +I(0) &\sim \text{LogNormal}(\mu = \log(80/0.05), \sigma = 1.5) +\end{align*} +$$ + +The reproductive number $R(t)$ is modeled as a random walk process: + +$$ +\begin{align*} +R(t) & = R(t-1) + \epsilon\\ +\log{\epsilon} & \sim \text{Normal}(\mu=0, \sigma=0.1) \\ +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: plot-hospital admissions +import matplotlib.pyplot as plt + +# Rotating the x-axis labels, and only showing ~10 labels +ax = plt.gca() +ax.xaxis.set_major_locator(plt.MaxNLocator(nbins=10)) +ax.xaxis.set_tick_params(rotation=45) +plt.plot(dat["date"].to_numpy(), dat["daily_hosp_admits"].to_numpy()) +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 hospitalization interval. + +```{python} +#| label: data-extract +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 = gen_int["probability_mass"].to_numpy() +inf_hosp_int = inf_hosp_int["probability_mass"].to_numpy() + +# Taking a pick at the first 5 elements of each +gen_int[:5], inf_hosp_int[: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) +axs[1].set_title("Infection to hospitalization interval") +``` + +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 +import jax.numpy as jnp +import numpyro.distributions as dist + +inf_hosp_int = deterministic.DeterministicPMF(inf_hosp_int) + +hosp_rate = latent.InfectHospRate( + dist=dist.LogNormal(jnp.log(0.05), 0.1) +) + +latent_hosp = latent.HospitalAdmissions( + infection_to_admission_interval=inf_hosp_int, + infect_hosp_rate_dist=hosp_rate, + ) +``` + +The `inf_hosp_int` is a `DeterministicPMF` object that takes the infection to hospitalization interval as input. The `hosp_rate` is an `InfectHospRate` object that takes the infection to hospitalization rate as input. The `HospitalAdmissions` class is a `RandomVariable` that takes two distributions as inputs: the infection to admission interval and the infection to hospitalization rate. Now, we can define the rest of the other components: + +```{python} +#| label: initializing-rest-of-model +from pyrenew import model, process, observation + +# Infection process +latent_inf = latent.Infections() +I0 = latent.Infections0(I0_dist=dist.LogNormal(loc=jnp.log(80/.05), scale=1.5)) + +# Generation interval and Rt +gen_int = deterministic.DeterministicPMF(gen_int) +rtproc = process.RtRandomWalkProcess( + Rt_rw_dist=dist.Normal(0, 0.1) +) + +# The observation model +obs = observation.NegativeBinomialObservation(concentration_prior=1.0) +``` + +Notice all the components are `RandomVariable` instances. We can now build the model: + +```{python} +#| label: init-model +hosp_model = model.HospitalizationsModel( + latent_infections=latent_inf, + latent_hospitalizations=latent_hosp, + I0=I0, + gen_int=gen_int, + Rt_process=rtproc, + observation_process=obs, +) +``` + +Let's simulate to check if the model is working: + +```{python} +#| label: simulation +import numpyro as npro +import numpy as np + +timeframe = 120 + +np.random.seed(223) +with npro.handlers.seed(rng_seed = np.random.randint(1, timeframe)): + sim_data = hosp_model.sample(n_timepoints=timeframe) +``` + +```{python} +#| label: basic-fig +#| fig-cap: Rt and Infections +#| fig-cap-location: bottom +import matplotlib.pyplot as plt + +fig, axs = plt.subplots(1, 2) + +# Rt plot +axs[0].plot(range(0, timeframe + 1), sim_data.Rt) +axs[0].set_ylabel('Rt') + +# Infections plot +axs[1].plot(range(0, timeframe + 1), sim_data.sampled_admissions) +axs[1].set_ylabel('Infections') +axs[1].set_yscale('log') + +fig.suptitle('Basic renewal model') +fig.supxlabel('Time') +plt.tight_layout() +plt.show() +``` + +## Fitting the model + +We can fit the model to the data. We will use the `run` method of the model object. The two inputs this model requires are `n_timepoints` and `observed_hospitalizations` + + +```{python} +#| label: model-fit +import jax + +hosp_model.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat.shape[0] - 1, + observed_hospitalizations=dat["daily_hosp_admits"].to_numpy(), + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), +) +``` + +We can use the `plot_posterior` method to visualize the results[^capture]: + +[^capture]: The output is captured to avoid `quarto` from displaying the output twice. + +```{python} +#| label: output-hospital admissions +#| fig-cap: Hospital Admissions posterior distribution +#| fig-cap-location: bottom +out = hosp_model.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat["daily_hosp_admits"].to_numpy(), +) +``` + +The first half of the model is not looking good. The reason is that the infection to hospitalization interval PMF makes it unlikely to observe admissions from the beginning. The following section shows how to fix this. + +## Padding the model + +We can use the padding argument to solve the overestimation of hospitalizations in the first half of the model. By setting `padding > 0`, the model then assumes that the first `padding` observations are missing; thus, only observations after `padding` will count towards the likelihood of the model. In practice, the model will extend the estimated Rt and latent infections by `padding` days, given time to adjust to the observed data. The following code will add 21 days of missing data at the beginning of the model and re-estimate it with `padding = 21`: + +```{python} +#| label: model-fit-padding +days_to_impute = 21 + +dat2 = dat["daily_hosp_admits"].to_numpy() + +# Add 21 Nas to the beginning of dat2 +dat2 = np.hstack((np.repeat(np.nan, days_to_impute), dat2)) + +hosp_model.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat2.shape[0] - 1, + observed_hospitalizations=dat2, + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), + padding=days_to_impute, # Padding the model +) +``` + +And plotting the results: + +```{python} +#| label: output-admissions-with-padding +#| fig-cap: Hospital Admissions posterior distribution +#| fig-cap-location: bottom +out = hosp_model.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat2, +) +``` + +We can also take a look at the latent infections: + +```{python} +#| label: output-infections-with-padding +#| fig-cap: Hospital Admissions posterior distribution +#| fig-cap-location: bottom +out2 = hosp_model.plot_posterior( + var="latent_infections", + ylab="Latent Infections" +) +``` + +## Round 2: Incorporating weekday 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 weekday effect distribution. To do this, we will create a new instance of `RandomVariable` to model the weekday effect. The weekday effect will be a truncated normal distribution with a mean of 1.0 and a standard deviation of 0.5. The distribution will be truncated between 0.1 and 10.0. The weekday effect will be repeated for the number of weeks in the dataset. + +```{python} +#| label: weekly-effect +from pyrenew import metaclass +import numpyro as npro + +class WeekdayEffect(metaclass.RandomVariable): + """Weekday effect distribution""" + def __init__(self, len: int): + """ Initialize the weekday effect distribution + Parameters + ---------- + len : int + The number of observations + """ + self.nweeks = jnp.ceil(len/7).astype(int) + self.len = len + + @staticmethod + def validate(): + return None + + def sample(self, **kwargs): + ans = npro.sample( + name="weekday_effect", + fn=npro.distributions.TruncatedNormal( + loc=1.0, scale=.5, low=0.1, high=10.0 + ), + sample_shape=(7,) + ) + + return jnp.tile(ans, self.nweeks)[:self.len] + +# Initializing the weekday effect +weekday_effect = WeekdayEffect(dat.shape[0]) +``` + +Notice that the instance's `nweeks` and `len` members are passed during construction. Trying to compute the number of weeks and the length of the dataset in the `validate` method will raise a `jit` error in `jax` as the shape and size of elements are not known during the validation step, which happens before the model is run. With the new weekday effect, we can rebuild the latent hospitalization model: + +```{python} +#| label: latent-hosp-weekday +latent_hosp_wday_effect = latent.HospitalAdmissions( + infection_to_admission_interval=inf_hosp_int, + infect_hosp_rate_dist=hosp_rate, + weekday_effect_dist=weekday_effect, + ) + +hosp_model_weekday = model.HospitalizationsModel( + latent_infections=latent_inf, + latent_hospitalizations=latent_hosp_wday_effect, + I0=I0, + gen_int=gen_int, + Rt_process=rtproc, + observation_process=obs, +) +``` + +Running the model (with the same padding as before): + +```{python} +#| label: model-2-run +hosp_model_weekday.run( + num_samples=2000, + num_warmup=2000, + n_timepoints=dat2.shape[0] - 1, + observed_hospitalizations=dat2, + rng_key=jax.random.PRNGKey(54), + mcmc_args=dict(progress_bar=False), + padding=days_to_impute, +) +``` + +And plotting the results: + +```{python} +#| label: output-admissions-padding-and-weekday +#| fig-cap: Hospital Admissions posterior distribution +#| fig-cap-location: bottom +#| +out = hosp_model_weekday.plot_posterior( + var="predicted_hospitalizations", + ylab="Hospital Admissions", + obs_signal=dat2, +) +``` diff --git a/model/docs/example-with-datasets_files/figure-commonmark/basic-fig-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/basic-fig-output-1.png new file mode 100644 index 00000000..5208cbd1 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/basic-fig-output-1.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/data-extract-output-2.png b/model/docs/example-with-datasets_files/figure-commonmark/data-extract-output-2.png new file mode 100644 index 00000000..ae28a9d7 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/data-extract-output-2.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-padding-and-weekday-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-padding-and-weekday-output-1.png new file mode 100644 index 00000000..f66b7c30 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-padding-and-weekday-output-1.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-with-padding-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-with-padding-output-1.png new file mode 100644 index 00000000..561bde1a Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/output-admissions-with-padding-output-1.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/output-hospital-admissions-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/output-hospital-admissions-output-1.png new file mode 100644 index 00000000..3f52b5e7 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/output-hospital-admissions-output-1.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/output-infections-with-padding-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/output-infections-with-padding-output-1.png new file mode 100644 index 00000000..b94f6e82 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/output-infections-with-padding-output-1.png differ diff --git a/model/docs/example-with-datasets_files/figure-commonmark/plot-hospital-admissions-output-1.png b/model/docs/example-with-datasets_files/figure-commonmark/plot-hospital-admissions-output-1.png new file mode 100644 index 00000000..687f8236 Binary files /dev/null and b/model/docs/example-with-datasets_files/figure-commonmark/plot-hospital-admissions-output-1.png differ diff --git a/model/docs/getting-started.md b/model/docs/getting-started.md index 905219a8..c15e510e 100644 --- a/model/docs/getting-started.md +++ b/model/docs/getting-started.md @@ -75,9 +75,7 @@ initialize the five components: ``` python # (1) The generation interval (deterministic) -gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), -) +gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) # (2) Initial infections (inferred with a prior) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -89,7 +87,7 @@ rt_proc = RtRandomWalkProcess() latent_infections = Infections() # (5) The observed infections process (with mean at the latent infections) -observed_infections = PoissonObservation() +observation_process = PoissonObservation() ``` With these five pieces, we can build the basic renewal model: @@ -100,7 +98,7 @@ model1 = RtInfectionsRenewalModel( I0 = I0, Rt_process = rt_proc, latent_infections = latent_infections, - observed_infections = observed_infections, + observation_process = observation_process, ) ``` @@ -116,7 +114,7 @@ flowchart TB i0["(2) I0\n(Infections0)"] rt["(3) rt_proc\n(RtRandomWalkProcess)"] inf["(4) latent_infections\n(Infections)"] - obs["(5) observed_infections\n(PoissonObservation)"] + obs["(5) observation_process\n(PoissonObservation)"] model1["model1\n(RtInfectionsRenewalModel)"] @@ -143,13 +141,13 @@ sim_data 1.271196 , 1.3189521, 1.3054799, 1.3165426, 1.291952 , 1.3026639, 1.2619467, 1.2852622, 1.3121517, 1.2888998, 1.2641873, 1.2580931, 1.2545817, 1.3092988, 1.2488269, 1.2397509, 1.2071848, 1.2334517, - 1.21868 ], dtype=float32), latent=Array([ 2.3215084, 3.0415602, 4.0327816, 5.180868 , 4.381411 , + 1.21868 ], dtype=float32), latent_infections=Array([ 2.3215084, 3.0415602, 4.0327816, 5.180868 , 4.381411 , 4.978916 , 5.750626 , 6.3024273, 6.66758 , 7.354823 , 7.8755097, 8.416656 , 9.63394 , 10.973988 , 12.043082 , 13.516833 , 14.911659 , 16.75407 , 18.053928 , 20.318869 , 22.975292 , 25.166464 , 27.34265 , 30.13236 , 33.126217 , 37.89362 , 40.11695 , 43.784634 , 46.754696 , 51.974545 , - 55.642136 ], dtype=float32), observed=Array([ 1, 2, 3, 5, 4, 4, 7, 4, 8, 4, 7, 3, 8, 12, 13, 18, 14, + 55.642136 ], dtype=float32), sampled_infections=Array([ 1, 2, 3, 5, 4, 4, 7, 4, 8, 4, 7, 3, 8, 12, 13, 18, 14, 20, 17, 18, 28, 27, 36, 37, 26, 31, 40, 27, 48, 54, 60], dtype=int32)) The `sample()` method of the `RtInfectionsRenewalModel` returns a list @@ -161,11 +159,11 @@ import matplotlib.pyplot as plt fig, axs = plt.subplots(1, 2) # Rt plot -axs[0].plot(range(0, 31), sim_data[0]) +axs[0].plot(range(0, 31), sim_data.Rt) axs[0].set_ylabel('Rt') # Infections plot -axs[1].plot(range(0, 31), sim_data[1]) +axs[1].plot(range(0, 31), sim_data.sampled_infections) axs[1].set_ylabel('Infections') fig.suptitle('Basic renewal model') @@ -187,7 +185,7 @@ import jax model1.run( num_warmup=2000, num_samples=1000, - observed_infections=sim_data.observed, + observed_infections=sim_data.sampled_infections, n_timepoints = len(sim_data[1])-1, rng_key=jax.random.PRNGKey(54), mcmc_args=dict(progress_bar=False), @@ -203,7 +201,7 @@ samps = model1.spread_draws([('Rt', 'time')]) fig, ax = plt.subplots(figsize=[4, 5]) -ax.plot(sim_data[0]) +ax.plot(sim_data.Rt) samp_ids = np.random.randint(size=25, low=0, high=999) for samp_id in samp_ids: sub_samps = samps.filter(pl.col("draw") == samp_id).sort(pl.col('time')) @@ -243,7 +241,7 @@ flowchart LR models((Model\nmetaclass)) subgraph observations[Observations module] - obs["observed_infections\n(PoissonObservation)"] + obs["observation_process\n(PoissonObservation)"] end subgraph latent[Latent module] diff --git a/model/docs/getting-started.qmd b/model/docs/getting-started.qmd index 12c96486..8ca96a11 100644 --- a/model/docs/getting-started.qmd +++ b/model/docs/getting-started.qmd @@ -53,9 +53,7 @@ The basic renewal model defines five components: generation interval, initial in ```{python} #| label: creating-elements # (1) The generation interval (deterministic) -gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), -) +gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) # (2) Initial infections (inferred with a prior) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -67,7 +65,7 @@ rt_proc = RtRandomWalkProcess() latent_infections = Infections() # (5) The observed infections process (with mean at the latent infections) -observed_infections = PoissonObservation() +observation_process = PoissonObservation() ``` With these five pieces, we can build the basic renewal model: @@ -79,7 +77,7 @@ model1 = RtInfectionsRenewalModel( I0 = I0, Rt_process = rt_proc, latent_infections = latent_infections, - observed_infections = observed_infections, + observation_process = observation_process, ) ``` @@ -92,7 +90,7 @@ flowchart TB i0["(2) I0\n(Infections0)"] rt["(3) rt_proc\n(RtRandomWalkProcess)"] inf["(4) latent_infections\n(Infections)"] - obs["(5) observed_infections\n(PoissonObservation)"] + obs["(5) observation_process\n(PoissonObservation)"] model1["model1\n(RtInfectionsRenewalModel)"] @@ -125,11 +123,11 @@ import matplotlib.pyplot as plt fig, axs = plt.subplots(1, 2) # Rt plot -axs[0].plot(range(0, 31), sim_data[0]) +axs[0].plot(range(0, 31), sim_data.Rt) axs[0].set_ylabel('Rt') # Infections plot -axs[1].plot(range(0, 31), sim_data[1]) +axs[1].plot(range(0, 31), sim_data.sampled_infections) axs[1].set_ylabel('Infections') fig.suptitle('Basic renewal model') @@ -147,7 +145,7 @@ import jax model1.run( num_warmup=2000, num_samples=1000, - observed_infections=sim_data.observed, + observed_infections=sim_data.sampled_infections, n_timepoints = len(sim_data[1])-1, rng_key=jax.random.PRNGKey(54), mcmc_args=dict(progress_bar=False), @@ -165,7 +163,7 @@ samps = model1.spread_draws([('Rt', 'time')]) fig, ax = plt.subplots(figsize=[4, 5]) -ax.plot(sim_data[0]) +ax.plot(sim_data.Rt) samp_ids = np.random.randint(size=25, low=0, high=999) for samp_id in samp_ids: sub_samps = samps.filter(pl.col("draw") == samp_id).sort(pl.col('time')) @@ -195,7 +193,7 @@ flowchart LR models((Model\nmetaclass)) subgraph observations[Observations module] - obs["observed_infections\n(PoissonObservation)"] + obs["observation_process\n(PoissonObservation)"] end subgraph latent[Latent module] diff --git a/model/docs/getting-started_files/figure-commonmark/basic-fig-output-1.png b/model/docs/getting-started_files/figure-commonmark/basic-fig-output-1.png index a751b6c4..6aec59a4 100644 Binary files a/model/docs/getting-started_files/figure-commonmark/basic-fig-output-1.png and b/model/docs/getting-started_files/figure-commonmark/basic-fig-output-1.png differ diff --git a/model/docs/getting-started_files/figure-commonmark/output-rt-output-1.png b/model/docs/getting-started_files/figure-commonmark/output-rt-output-1.png index c67d4967..af600fea 100644 Binary files a/model/docs/getting-started_files/figure-commonmark/output-rt-output-1.png and b/model/docs/getting-started_files/figure-commonmark/output-rt-output-1.png differ diff --git a/model/docs/pyrenew_demo.md b/model/docs/pyrenew_demo.md index 3aace074..76eb0444 100644 --- a/model/docs/pyrenew_demo.md +++ b/model/docs/pyrenew_demo.md @@ -8,21 +8,22 @@ Assuming you’ve already installed Python and pip, you’ll need to first install `pyrenew`: ``` python -python3 -m pip install "pyrenew" +pip install pyrenew ``` You’ll also need working installations of `matplotlib`, `numpy`, `jax`, `numpyro`, and `polars`: ``` python -python -m pip install "matplotlib" "numpy" "jax" "numpyro" "polars" +pip install matplotlib numpy jax numpyro polars ``` -Run the following import section to call external modules and functions -necessary to run the `pyrenew` demo. The `import` statement imports the -module and the `as` statement renames the module for use within this -script. The `from` statement imports a specific function from a module -(named after the `.`) within a package (named before the `.`). +To begin, run the following import section to call external modules and +functions necessary to run the `pyrenew` demo. The `import` statement +imports the module and the `as` statement renames the module for use +within this script. The `from` statement imports a specific function +from a module (named after the `.`) within a package (named before the +`.`). ``` python import matplotlib as mpl @@ -37,7 +38,21 @@ import numpyro.distributions as dist ``` python from pyrenew.process import SimpleRandomWalkProcess ``` -To understand the simple random walk process underlying the sampling within the renewal process model, we first examine a single random walk path. Using the `sample` method from an instance of the `SimpleRandomWalkProcess` class, we first create an instance of the `SimpleRandomWalkProcess` class with a normal distribution of mean = 0 and standard deviation = 0.0001 as its input. Next, the `with` statement sets the seed for the random number generator for the duration of the block that follows. Inside the `with` block, the `q_samp = q.sample(duration=100)` generates the sample instance over a duration of 100 time units. Finally, this single random walk process is visualized using `matplot.pyplot` to plot the exponential of the sample instance. + + An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu. + +To understand the simple random walk process underlying the sampling +within the renewal process model, we first examine a single random walk +path. Using the `sample` method from an instance of the +`SimpleRandomWalkProcess` class, we first create an instance of the +`SimpleRandomWalkProcess` class with a normal distribution of mean = 0 +and standard deviation = 0.0001 as its input. Next, the `with` statement +sets the seed for the random number generator for the duration of the +block that follows. Inside the `with` block, the +`q_samp = q.sample(duration=100)` generates the sample instance over a +duration of 100 time units. Finally, this single random walk process is +visualized using `matplot.pyplot` to plot the exponential of the sample +instance. ``` python np.random.seed(3312) @@ -50,45 +65,60 @@ plt.plot(np.exp(q_samp[0])) ![](pyrenew_demo_files/figure-commonmark/fig-randwalk-output-1.png) +Next, import several additional functions from the `latent` module of +the `pyrenew` package to model infections, hospital admissions, initial +infections, and hospitalization rate due to infection. -Next, import several additional functions from the `latent` module of the `pyrenew` package to model infections, hospital admissions, initial infections, and hospitalization rate due to infection. - -```{python} +``` python from pyrenew.latent import ( Infections, HospitalAdmissions, Infections0, InfectHospRate, ) ``` -Additionally, import several classes from Pyrenew, including a Poisson observation process, determininstic PMF and variable classes, the Pyrenew hospitalization model, and a renewal modle (Rt) random walk process: +Additionally, import several classes from Pyrenew, including a Poisson +observation process, determininstic PMF and variable classes, the +Pyrenew hospitalization model, and a renewal modle (Rt) random walk +process: -```{python} +``` python from pyrenew.observation import PoissonObservation from pyrenew.deterministic import DeterministicPMF, DeterministicVariable from pyrenew.model import HospitalizationsModel from pyrenew.process import RtRandomWalkProcess ``` -To initialize a model run, we first define initial conditions, including: +To initialize the model, we first define initial conditions, including: -1) deterministic generation time, defined as an instance of the `DeterministicPMF` class, which gives the probability of each possible outcome for a discrete random variable given as a JAX NumPy array of four possible outcomes +1) deterministic generation time, defined as an instance of the + `DeterministicPMF` class, which gives the probability of each + possible outcome for a discrete random variable given as a JAX NumPy + array of four possible outcomes -2) initial infections at the start of simulation as a log-normal distribution with mean = 0 and standard deviation = 1 +2) initial infections at the start of simulation as a log-normal + distribution with mean = 0 and standard deviation = 1 -3) latent infections as an instance of the `Infections` class with default settings +3) latent infections as an instance of the `Infections` class with + default settings -4) latent hospitalization process, modeled by first defining the time interval from infections to hospitalizations as a `DeterministicPMF` input with 18 possible outcomes and corresponding probabilities given by the values in the array. The `HospitalAdmissions` function then takes in this defined time interval, as well as defining the rate at which infections are admitted to the hospital due to infection, modeled as a log-normal distribution with mean = `jnp.log(0.05)` and standard deviation = 0.05. +4) latent hospitalization process, modeled by first defining the time + interval from infections to hospitalizations as a `DeterministicPMF` + input with 18 possible outcomes and corresponding probabilities + given by the values in the array. The `HospitalAdmissions` function + then takes in this defined time interval, as well as defining the + rate at which infections are admitted to the hospital due to + infection, modeled as a log-normal distribution with mean = + `jnp.log(0.05)` and standard deviation = 0.05. -5) hospitalization observation process, modeled with a Poisson distribution +5) hospitalization observation process, modeled with a Poisson + distribution -6) an Rt random walk process with default settings +6) an Rt random walk process with default settings ``` python # Initializing model components: # 1) A deterministic generation time -gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) +gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) # 2) Initial infections I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -100,7 +130,7 @@ latent_infections = Infections() # First, define a deterministic infection to hosp pmf inf_hosp_int = DeterministicPMF( - (jnp.array([0, 0, 0,0,0,0,0,0,0,0,0,0,0, 0.25, 0.5, 0.1, 0.1, 0.05]),), + jnp.array([0, 0, 0,0,0,0,0,0,0,0,0,0,0, 0.25, 0.5, 0.1, 0.1, 0.05]), ) latent_hospitalizations = HospitalAdmissions( @@ -111,14 +141,15 @@ latent_hospitalizations = HospitalAdmissions( ) # 5) An observation process for the hospitalizations -observed_hospitalizations = PoissonObservation() +obs_process = PoissonObservation() # 6) A random walk process (it could be deterministic using # pyrenew.process.DeterministicProcess()) Rt_process = RtRandomWalkProcess() ``` -The `HospitalizationsModel` is then initialized using the initial conditions just defined: +The `HospitalizationsModel` is then initialized using the initial +conditions just defined: ``` python # Initializing the model @@ -126,12 +157,15 @@ hospmodel = HospitalizationsModel( gen_int=gen_int, I0=I0, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=observed_hospitalizations, + observation_process=obs_process, latent_infections=latent_infections, Rt_process=Rt_process ) ``` +Next, we sample from the `hospmodel` for 30 time steps and view the +output of a single run: + ``` python with seed(rng_seed=np.random.randint(1, 60)): x = hospmodel.sample(n_timepoints=30) @@ -143,92 +177,115 @@ x 1.2466507, 1.2800207, 1.2749145, 1.2619376, 1.2189837, 1.2192641, 1.2290158, 1.2128737, 1.1908046, 1.2174997, 1.1941082, 1.2084603, 1.1965215, 1.2248698, 1.2308019, 1.2426206, 1.2131014, 1.207159 , - 1.1837622], dtype=float32), infections=Array([0.05214045, 0.06867922, 0.08761451, 0.11476436, 0.09757317, + 1.1837622], dtype=float32), latent_infections=Array([0.05214045, 0.06867922, 0.08761451, 0.11476436, 0.09757317, 0.10547114, 0.1167062 , 0.13010225, 0.13824694, 0.14372033, 0.15924728, 0.17601486, 0.19236736, 0.21483542, 0.23664482, 0.25865382, 0.27503362, 0.30029488, 0.3289544 , 0.35262382, 0.37418258, 0.41274938, 0.43839005, 0.47672123, 0.50913286, 0.5625195 , 0.6113282 , 0.67092246, 0.7138808 , 0.77217466, - 0.819254 ], dtype=float32), IHR=Array(0.04929917, dtype=float32), latent=Array([0. , 0. , 0. , 0. , 0. , + 0.819254 ], dtype=float32), IHR=Array(0.04929917, dtype=float32), latent_admissions=Array([0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.00064262, 0.0021317 , 0.00302979, 0.00416974, 0.0049305 , 0.00487205, 0.00530097, 0.00576412, 0.00624666, 0.00665578, 0.00711595, 0.0078055 , 0.00854396, 0.00939666, 0.01042083, 0.0114624 , 0.01246538, - 0.01345188], dtype=float32), sampled=Array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.01345188], dtype=float32), sampled_admissions=Array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], dtype=int32)) +Visualizations of the single model output show (top) infections over the +30 time steps, (middle) hospitalizations over the 30 time steps, and +(bottom) + ``` python fig, ax = plt.subplots(nrows=3, sharex=True) -ax[0].plot(x.infections) +ax[0].plot(x.latent_infections) ax[0].set_ylim([1/5, 5]) -ax[1].plot(x.latent) -ax[2].plot(x.sampled, 'o') +ax[1].plot(x.latent_admissions) +ax[2].plot(x.sampled_admissions, 'o') for axis in ax[:-1]: axis.set_yscale("log") ``` ![](pyrenew_demo_files/figure-commonmark/fig-hosp-output-1.png) +To fit the `hospmodel` to the simulated data, we call `hospmodel.run()`, +an MCMC algorithm, with the arguments generated in `hospmodel` object, +using 1000 warmup stepts and 1000 samples to draw from the posterior +distribution of the model parameters. The model is run for +`len(x.sampled_admissions)-1` time steps with the seed set by +`jax.random.PRNGKey()` + ``` python # from numpyro.infer import MCMC, NUTS hospmodel.run( num_warmup=1000, num_samples=1000, - observed_hospitalizations=x.sampled, - n_timepoints = len(x.sampled)-1, + observed_hospitalizations=x.sampled_admissions, + n_timepoints = len(x.sampled_admissions)-1, rng_key=jax.random.PRNGKey(54), mcmc_args=dict(progress_bar=False), ) ``` +Print a summary of the model: + ``` python hospmodel.print_summary() ``` mean std median 5.0% 95.0% n_eff r_hat - I0 1.26 1.09 0.96 0.09 2.41 1114.64 1.00 - IHR 0.05 0.00 0.05 0.05 0.05 2747.53 1.00 - Rt0 1.23 0.17 1.23 0.95 1.52 1533.77 1.00 - Rt_transformed_rw_diffs[0] 0.00 0.03 -0.00 -0.04 0.04 1574.38 1.00 - Rt_transformed_rw_diffs[1] 0.00 0.03 0.00 -0.04 0.04 2557.65 1.00 - Rt_transformed_rw_diffs[2] 0.00 0.02 0.00 -0.04 0.04 2245.16 1.00 - Rt_transformed_rw_diffs[3] 0.00 0.02 0.00 -0.03 0.04 2423.80 1.00 - Rt_transformed_rw_diffs[4] 0.00 0.02 0.00 -0.03 0.04 2461.65 1.00 - Rt_transformed_rw_diffs[5] 0.00 0.02 0.00 -0.04 0.04 2363.57 1.00 - Rt_transformed_rw_diffs[6] 0.00 0.02 0.00 -0.04 0.04 1720.93 1.00 - Rt_transformed_rw_diffs[7] -0.00 0.02 -0.00 -0.04 0.04 3851.66 1.00 - Rt_transformed_rw_diffs[8] -0.00 0.03 -0.00 -0.04 0.04 1824.74 1.00 - Rt_transformed_rw_diffs[9] 0.00 0.02 0.00 -0.04 0.04 1739.32 1.00 - Rt_transformed_rw_diffs[10] 0.00 0.02 0.00 -0.04 0.04 1944.43 1.00 - Rt_transformed_rw_diffs[11] -0.00 0.03 0.00 -0.04 0.04 1558.14 1.00 - Rt_transformed_rw_diffs[12] 0.00 0.02 0.00 -0.04 0.04 2182.35 1.00 - Rt_transformed_rw_diffs[13] -0.00 0.03 0.00 -0.04 0.04 1175.35 1.00 - Rt_transformed_rw_diffs[14] -0.00 0.03 -0.00 -0.04 0.04 1540.25 1.00 - Rt_transformed_rw_diffs[15] -0.00 0.03 -0.00 -0.04 0.04 2367.82 1.00 - Rt_transformed_rw_diffs[16] -0.00 0.02 -0.00 -0.04 0.04 1636.30 1.00 - Rt_transformed_rw_diffs[17] 0.00 0.03 0.00 -0.04 0.04 1978.96 1.00 - Rt_transformed_rw_diffs[18] 0.00 0.02 -0.00 -0.04 0.04 1589.27 1.00 - Rt_transformed_rw_diffs[19] 0.00 0.03 -0.00 -0.04 0.04 1691.06 1.00 - Rt_transformed_rw_diffs[20] -0.00 0.02 -0.00 -0.04 0.04 2562.99 1.00 - Rt_transformed_rw_diffs[21] 0.00 0.02 -0.00 -0.04 0.04 2352.40 1.00 - Rt_transformed_rw_diffs[22] 0.00 0.03 0.00 -0.04 0.04 1971.40 1.00 - Rt_transformed_rw_diffs[23] 0.00 0.02 0.00 -0.04 0.04 2013.90 1.00 - Rt_transformed_rw_diffs[24] 0.00 0.03 0.00 -0.04 0.04 2022.94 1.00 - Rt_transformed_rw_diffs[25] -0.00 0.02 -0.00 -0.04 0.03 1981.62 1.00 - Rt_transformed_rw_diffs[26] 0.00 0.03 0.00 -0.04 0.05 2696.36 1.00 - Rt_transformed_rw_diffs[27] -0.00 0.03 0.00 -0.04 0.04 2003.38 1.00 - Rt_transformed_rw_diffs[28] -0.00 0.02 -0.00 -0.04 0.04 1843.27 1.00 - Rt_transformed_rw_diffs[29] -0.00 0.03 -0.00 -0.04 0.04 1780.88 1.00 + I0 1.27 1.10 0.97 0.10 2.42 1132.34 1.00 + IHR 0.05 0.00 0.05 0.05 0.05 2306.45 1.00 + Rt0 1.23 0.17 1.23 0.93 1.48 1327.22 1.00 + Rt_transformed_rw_diffs[0] -0.00 0.02 -0.00 -0.04 0.04 1404.95 1.00 + Rt_transformed_rw_diffs[1] 0.00 0.03 0.00 -0.04 0.04 2280.86 1.00 + Rt_transformed_rw_diffs[2] -0.00 0.02 -0.00 -0.04 0.04 2119.83 1.00 + Rt_transformed_rw_diffs[3] 0.00 0.02 -0.00 -0.04 0.04 2196.86 1.00 + Rt_transformed_rw_diffs[4] 0.00 0.02 -0.00 -0.03 0.04 2391.45 1.00 + Rt_transformed_rw_diffs[5] 0.00 0.03 0.00 -0.04 0.04 2043.02 1.00 + Rt_transformed_rw_diffs[6] 0.00 0.02 0.00 -0.04 0.04 1514.40 1.00 + Rt_transformed_rw_diffs[7] -0.00 0.02 -0.00 -0.04 0.04 2619.69 1.00 + Rt_transformed_rw_diffs[8] 0.00 0.03 0.00 -0.04 0.04 1883.84 1.00 + Rt_transformed_rw_diffs[9] 0.00 0.03 0.00 -0.04 0.04 2015.66 1.00 + Rt_transformed_rw_diffs[10] 0.00 0.02 0.00 -0.04 0.04 2045.47 1.00 + Rt_transformed_rw_diffs[11] -0.00 0.03 0.00 -0.04 0.04 1615.10 1.00 + Rt_transformed_rw_diffs[12] 0.00 0.02 0.00 -0.04 0.04 2206.32 1.00 + Rt_transformed_rw_diffs[13] 0.00 0.03 0.00 -0.04 0.04 1175.93 1.00 + Rt_transformed_rw_diffs[14] -0.00 0.03 -0.00 -0.04 0.04 1606.26 1.00 + Rt_transformed_rw_diffs[15] -0.00 0.03 -0.00 -0.04 0.04 2344.62 1.00 + Rt_transformed_rw_diffs[16] -0.00 0.02 0.00 -0.04 0.04 1522.33 1.00 + Rt_transformed_rw_diffs[17] 0.00 0.03 0.00 -0.04 0.04 2157.17 1.00 + Rt_transformed_rw_diffs[18] -0.00 0.02 -0.00 -0.04 0.04 1594.95 1.00 + Rt_transformed_rw_diffs[19] 0.00 0.03 -0.00 -0.04 0.04 1698.70 1.00 + Rt_transformed_rw_diffs[20] 0.00 0.02 0.00 -0.04 0.04 1726.18 1.00 + Rt_transformed_rw_diffs[21] 0.00 0.02 -0.00 -0.04 0.04 2386.35 1.00 + Rt_transformed_rw_diffs[22] 0.00 0.03 0.00 -0.04 0.04 2028.63 1.00 + Rt_transformed_rw_diffs[23] 0.00 0.02 0.00 -0.04 0.03 1669.71 1.00 + Rt_transformed_rw_diffs[24] 0.00 0.02 0.00 -0.04 0.04 2126.33 1.00 + Rt_transformed_rw_diffs[25] -0.00 0.02 -0.00 -0.04 0.04 2119.74 1.00 + Rt_transformed_rw_diffs[26] 0.00 0.03 0.00 -0.04 0.04 2657.91 1.00 + Rt_transformed_rw_diffs[27] -0.00 0.03 0.00 -0.04 0.04 1939.30 1.00 + Rt_transformed_rw_diffs[28] -0.00 0.02 -0.00 -0.04 0.04 1737.84 1.00 + Rt_transformed_rw_diffs[29] -0.00 0.03 -0.00 -0.04 0.04 2105.55 1.00 Number of divergences: 0 +Next, we will use the `spread_draws` function from the +`pyrenew.mcmcutils` module to process the MCMC samples. The +`spread_draws` function reformats the samples drawn from the +`mcmc.get_samples()` from the `hospmodel`. The samples are simulated Rt +values over time. + ``` python from pyrenew.mcmcutils import spread_draws samps = spread_draws(hospmodel.mcmc.get_samples(), [("Rt", "time")]) ``` +We visualize these samples below, with individual possible Rt estimates +over time shown in light blue, and the overall mean estimate Rt shown in +dark blue. + ``` python import numpy as np import polars as pl diff --git a/model/docs/pyrenew_demo.qmd b/model/docs/pyrenew_demo.qmd index 2cc39404..6bea20dd 100644 --- a/model/docs/pyrenew_demo.qmd +++ b/model/docs/pyrenew_demo.qmd @@ -8,14 +8,14 @@ This demo simulates some basic renewal process data and then fits to it using `p Assuming you've already installed Python and pip, you’ll need to first install `pyrenew`: -```{python} +```python pip install pyrenew ``` You’ll also need working installations of `matplotlib`, `numpy`, `jax`, `numpyro`, and `polars`: -```{python} +```python pip install matplotlib numpy jax numpyro polars ``` @@ -86,9 +86,7 @@ To initialize the model, we first define initial conditions, including: # Initializing model components: # 1) A deterministic generation time -gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) +gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) # 2) Initial infections I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -100,7 +98,7 @@ latent_infections = Infections() # First, define a deterministic infection to hosp pmf inf_hosp_int = DeterministicPMF( - (jnp.array([0, 0, 0,0,0,0,0,0,0,0,0,0,0, 0.25, 0.5, 0.1, 0.1, 0.05]),), + jnp.array([0, 0, 0,0,0,0,0,0,0,0,0,0,0, 0.25, 0.5, 0.1, 0.1, 0.05]), ) latent_hospitalizations = HospitalAdmissions( @@ -111,7 +109,7 @@ latent_hospitalizations = HospitalAdmissions( ) # 5) An observation process for the hospitalizations -observed_hospitalizations = PoissonObservation() +obs_process = PoissonObservation() # 6) A random walk process (it could be deterministic using # pyrenew.process.DeterministicProcess()) @@ -126,7 +124,7 @@ hospmodel = HospitalizationsModel( gen_int=gen_int, I0=I0, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=observed_hospitalizations, + observation_process=obs_process, latent_infections=latent_infections, Rt_process=Rt_process ) @@ -146,23 +144,23 @@ Visualizations of the single model output show (top) infections over the 30 time #| label: fig-hosp #| fig-cap: Infections fig, ax = plt.subplots(nrows=3, sharex=True) -ax[0].plot(x.infections) +ax[0].plot(x.latent_infections) ax[0].set_ylim([1/5, 5]) -ax[1].plot(x.latent) -ax[2].plot(x.sampled, 'o') +ax[1].plot(x.latent_admissions) +ax[2].plot(x.sampled_admissions, 'o') for axis in ax[:-1]: axis.set_yscale("log") ``` -To fit the `hospmodel` to the simulated data, we call `hospmodel.run()`, an MCMC algorithm, with the arguments generated in `hospmodel` object, using 1000 warmup stepts and 1000 samples to draw from the posterior distribution of the model parameters. The model is run for `len(x.sampled)-1` time steps with the seed set by `jax.random.PRNGKey()` +To fit the `hospmodel` to the simulated data, we call `hospmodel.run()`, an MCMC algorithm, with the arguments generated in `hospmodel` object, using 1000 warmup stepts and 1000 samples to draw from the posterior distribution of the model parameters. The model is run for `len(x.sampled_admissions)-1` time steps with the seed set by `jax.random.PRNGKey()` ```{python} # from numpyro.infer import MCMC, NUTS hospmodel.run( num_warmup=1000, num_samples=1000, - observed_hospitalizations=x.sampled, - n_timepoints = len(x.sampled)-1, + observed_hospitalizations=x.sampled_admissions, + n_timepoints = len(x.sampled_admissions)-1, rng_key=jax.random.PRNGKey(54), mcmc_args=dict(progress_bar=False), ) diff --git a/model/docs/pyrenew_demo_files/figure-commonmark/fig-randwalk-output-1.png b/model/docs/pyrenew_demo_files/figure-commonmark/fig-randwalk-output-1.png index 9a794124..e3a6cc17 100644 Binary files a/model/docs/pyrenew_demo_files/figure-commonmark/fig-randwalk-output-1.png and b/model/docs/pyrenew_demo_files/figure-commonmark/fig-randwalk-output-1.png differ diff --git a/model/docs/pyrenew_demo_files/figure-commonmark/fig-sampled-rt-output-1.png b/model/docs/pyrenew_demo_files/figure-commonmark/fig-sampled-rt-output-1.png index e31e0699..7f94e989 100644 Binary files a/model/docs/pyrenew_demo_files/figure-commonmark/fig-sampled-rt-output-1.png and b/model/docs/pyrenew_demo_files/figure-commonmark/fig-sampled-rt-output-1.png differ diff --git a/model/poetry.lock b/model/poetry.lock index 395d7445..6ae2435b 100644 --- a/model/poetry.lock +++ b/model/poetry.lock @@ -1487,17 +1487,17 @@ testing = ["pytest", "pytest-benchmark"] [[package]] name = "polars" -version = "0.20.19" +version = "0.20.21" description = "Blazingly fast DataFrame library" optional = false python-versions = ">=3.8" files = [ - {file = "polars-0.20.19-cp38-abi3-macosx_10_12_x86_64.whl", hash = "sha256:1d5a63001c8e47376f35b6e389e7f88e6984c721bd9468b01998807bad3c0f2f"}, - {file = "polars-0.20.19-cp38-abi3-macosx_11_0_arm64.whl", hash = "sha256:6790e9e8ce0bb8dcefa4e494b5e4ee5eca01df0c7d7151dcc45a7b24aa58939b"}, - {file = "polars-0.20.19-cp38-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b9243752dd2f0c9077b2caa285e4b05d3f80558ad4b8ae448dc41a0b9d5f10a5"}, - {file = "polars-0.20.19-cp38-abi3-manylinux_2_24_aarch64.whl", hash = "sha256:a0d57955cdc3610ec644ffd249c94c3102ce32df31661b759ab33a6bb4fae65c"}, - {file = "polars-0.20.19-cp38-abi3-win_amd64.whl", hash = "sha256:20ff4990cbdc05fc5a526902561a16c66ead4d740ff8e7b3bf52d4f3f7def60a"}, - {file = "polars-0.20.19.tar.gz", hash = "sha256:58f9c3766c6fa3d7193364cfcf93c6c640e6c8e6a49ad8c6fd48b326ad8060f2"}, + {file = "polars-0.20.21-cp38-abi3-macosx_10_12_x86_64.whl", hash = "sha256:f1de9b7190111b0d6c3731157da26252697ba27f2eef201d01e121c6873e713e"}, + {file = "polars-0.20.21-cp38-abi3-macosx_11_0_arm64.whl", hash = "sha256:18847d198cbdb133c35847d096f0b72783c9cc81df0878959dba60495a99a1c7"}, + {file = "polars-0.20.21-cp38-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:32d5e7f5a1bf529ca0a3c14416a8802e746534c4220dc8af42b39ea4962e729a"}, + {file = "polars-0.20.21-cp38-abi3-manylinux_2_24_aarch64.whl", hash = "sha256:f8196b42791ca6d1de9448e985f0d32d783166a43928bd294486fc37e1c78310"}, + {file = "polars-0.20.21-cp38-abi3-win_amd64.whl", hash = "sha256:584d3fe8d14516de2d733b9a355ae27a1be943bb2d54e76826887a0e3e122ee3"}, + {file = "polars-0.20.21.tar.gz", hash = "sha256:1ef8c44dc1a0dd79f3c1a6f73677cddac6e5aca3093cc9fb666f2e098f9b3bef"}, ] [package.extras] @@ -1671,6 +1671,28 @@ pytest = ">=4.6" [package.extras] testing = ["fields", "hunter", "process-tests", "pytest-xdist", "virtualenv"] +[[package]] +name = "pytest-mpl" +version = "0.17.0" +description = "pytest plugin to help with testing figures output from Matplotlib" +optional = false +python-versions = ">=3.6" +files = [ + {file = "pytest-mpl-0.17.0.tar.gz", hash = "sha256:fbef05d7e664b4b33452fb698ac188e522791f327de34e7ea37dbdfe9d52cac6"}, + {file = "pytest_mpl-0.17.0-py3-none-any.whl", hash = "sha256:2ee7cf902ef071e45ae14fe190b3d8a9839c2fd8933e0514bf58df310d8ed66b"}, +] + +[package.dependencies] +Jinja2 = "*" +matplotlib = "*" +packaging = "*" +Pillow = "*" +pytest = "*" + +[package.extras] +docs = ["matplotlib (==3.6)", "mpl-sphinx-theme (>=3.6.0.dev0)", "sphinx", "sphinx-design"] +test = ["pytest-cov"] + [[package]] name = "python-dateutil" version = "2.9.0.post0" @@ -2345,4 +2367,4 @@ files = [ [metadata] lock-version = "2.0" python-versions = "^3.10" -content-hash = "efdc5cb1ee45f3e7239a90ec4210574bff65fcd21449b3020b918e8763ea81b3" +content-hash = "3fb2c4d09caa9055098619152c9a6f4727fd6a5ae25750e7d1040afba9d6b31e" diff --git a/model/pyproject.toml b/model/pyproject.toml index e6a749f1..fdc6a6a4 100755 --- a/model/pyproject.toml +++ b/model/pyproject.toml @@ -5,6 +5,9 @@ description = "Pyrenew: a package for Bayesian renewal modeling with JAX and Num authors = ["Dylan H. Morris "] license = "Apache 2.0" readme = "README.md" +include = [{path = "datasets/*.tsv"}] +exclude = [{path = "datasets/*.rds"}] + [tool.poetry.dependencies] python = "^3.10" @@ -25,6 +28,7 @@ numpydoc = "^1.6.0" nbclient = "^0.10.0" nbformat = "^5.10.0" pytest-cov = "^5.0.0" +pytest-mpl = "^0.17.0" [tool.numpydoc_validation] checks = [ diff --git a/model/src/pyrenew/datasets/__init__.py b/model/src/pyrenew/datasets/__init__.py new file mode 100644 index 00000000..bf31218e --- /dev/null +++ b/model/src/pyrenew/datasets/__init__.py @@ -0,0 +1,11 @@ +from pyrenew.datasets.generation_interval import load_generation_interval +from pyrenew.datasets.infection_admission_interval import ( + load_infection_admission_interval, +) +from pyrenew.datasets.wastewater import load_wastewater + +__all__ = [ + "load_wastewater", + "load_infection_admission_interval", + "load_generation_interval", +] diff --git a/model/src/pyrenew/datasets/data-raw/README.md b/model/src/pyrenew/datasets/data-raw/README.md new file mode 100644 index 00000000..d3053f86 --- /dev/null +++ b/model/src/pyrenew/datasets/data-raw/README.md @@ -0,0 +1,3 @@ +# Raw data + +This folder contains the scripts used to download raw data from the wastewater project and save it as tsv for the pyrenew project. diff --git a/model/src/pyrenew/datasets/data-raw/generation_interval.py b/model/src/pyrenew/datasets/data-raw/generation_interval.py new file mode 100644 index 00000000..75b4ea7a --- /dev/null +++ b/model/src/pyrenew/datasets/data-raw/generation_interval.py @@ -0,0 +1,26 @@ +# Dataset from: +# https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/generation_interval.csv +import os + +import polars as pl + +gen_int = pl.read_csv( + "https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/generation_interval.csv", +) + +# Building path to save the file +path = os.path.join( + "src", + "pyrenew", + "datasets", + "gen_int.tsv", +) + +os.makedirs(os.path.dirname(path), exist_ok=True) + +gen_int.write_csv( + file=path, + separator="\t", + include_header=True, + null_value="", +) diff --git a/model/src/pyrenew/datasets/data-raw/infection_admission_interval.py b/model/src/pyrenew/datasets/data-raw/infection_admission_interval.py new file mode 100644 index 00000000..6e95383e --- /dev/null +++ b/model/src/pyrenew/datasets/data-raw/infection_admission_interval.py @@ -0,0 +1,27 @@ +# Dataset from: +# https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/inf_to_hosp.csv +import os + +import polars as pl + +# Read CSV file from URL +url = "https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/inf_to_hosp.csv" +infection_admission_interval = pl.read_csv(url) + +# Building path to save the file +path = os.path.join( + "src", + "pyrenew", + "datasets", + "infection_admission_interval.tsv", +) + +os.makedirs(os.path.dirname(path), exist_ok=True) + +# Save as TSV +infection_admission_interval.write_csv( + file=path, + separator="\t", + include_header=True, + null_value="", +) diff --git a/model/src/pyrenew/datasets/data-raw/wastewater.R b/model/src/pyrenew/datasets/data-raw/wastewater.R new file mode 100644 index 00000000..5ff804f1 --- /dev/null +++ b/model/src/pyrenew/datasets/data-raw/wastewater.R @@ -0,0 +1,15 @@ +load( + "https://github.com/CDCgov/wastewater-informed-covid-forecasting/blob/292526383ece582f10823fc939c7e590ca349c6d/cfaforecastrenewalww/data/example_df.rda" + ) +datasets_dir <- file.path("src", "pyrenew", "datasets") +dir.create(datasets_dir) + +# Saving as TSV +write.table( + x = example_df, + file = file.path(datasets_dir, "wastewater.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE, + na = "" + ) diff --git a/model/src/pyrenew/datasets/generation_interval.py b/model/src/pyrenew/datasets/generation_interval.py new file mode 100644 index 00000000..11c8f9ea --- /dev/null +++ b/model/src/pyrenew/datasets/generation_interval.py @@ -0,0 +1,31 @@ +from importlib.resources import files + +import polars as pl + + +def load_generation_interval() -> pl.DataFrame: + """Load the generation interval dataset + + This dataset contains the generation interval distribution for COVID-19. + + Returns + ------- + pl.DataFrame + The generation interval dataset + + Notes + ----- + + This dataset was downloaded directly from: + https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/generation_interval.csv + + The dataset contains the following columns: + - `timepoint` + - `probability_mass` + """ + + # Load the dataset + return pl.read_csv( + source=files("pyrenew.datasets") / "generation_interval.tsv", + separator="\t", + ) diff --git a/model/src/pyrenew/datasets/generation_interval.tsv b/model/src/pyrenew/datasets/generation_interval.tsv new file mode 100644 index 00000000..f74292b5 --- /dev/null +++ b/model/src/pyrenew/datasets/generation_interval.tsv @@ -0,0 +1,16 @@ +timepoint probability_mass +1 0.16174201341157 +2 0.320625889390807 +3 0.242283351248733 +4 0.134652559431048 +5 0.0689218448080527 +6 0.0345873307236478 +7 0.0175305676468834 +8 0.00900188393721926 +9 0.0048814763363626 +10 0.00258678762113042 +11 0.00144458017530763 +12 0.000825330494803486 +13 0.000478236415587845 +14 0.000270544021118098 +15 0.000167604337729198 diff --git a/model/src/pyrenew/datasets/infection_admission_interval.py b/model/src/pyrenew/datasets/infection_admission_interval.py new file mode 100644 index 00000000..e04e4c37 --- /dev/null +++ b/model/src/pyrenew/datasets/infection_admission_interval.py @@ -0,0 +1,32 @@ +from importlib.resources import files + +import polars as pl + + +def load_infection_admission_interval() -> pl.DataFrame: + """Load the infection to admission interval + + This dataset contains the infection to admission interval distribution for + COVID-19. + + Returns + ------- + pl.DataFrame + The infection to admission interval dataset + + Notes + ----- + + This dataset was downloaded directly from: + https://raw.githubusercontent.com/CDCgov/wastewater-informed-covid-forecasting/0962c5d1652787479ac72caebf076ab55fe4e10c/input/saved_pmfs/inf_to_hosp.csv + + The dataset contains the following columns: + - `timepoint` + - `probability_mass` + """ + + # Load the dataset + return pl.read_csv( + source=files("pyrenew.datasets") / "infection_admission_interval.tsv", + separator="\t", + ) diff --git a/model/src/pyrenew/datasets/infection_admission_interval.tsv b/model/src/pyrenew/datasets/infection_admission_interval.tsv new file mode 100644 index 00000000..2168696f --- /dev/null +++ b/model/src/pyrenew/datasets/infection_admission_interval.tsv @@ -0,0 +1,56 @@ +timepoint probability_mass +0 0.0 +1 0.00469384736487552 +2 0.0145200073436112 +3 0.0278627741704387 +4 0.0423656492135518 +5 0.0558071445014868 +6 0.0665713169684116 +7 0.0737925805176124 +8 0.0772854627892072 +9 0.0773666390616176 +10 0.0746515449009949 +11 0.0698761436052596 +12 0.0637663813017696 +13 0.0569581929821651 +14 0.0499600186601535 +15 0.0431457477049282 +16 0.0367662806214045 +17 0.0309702535668237 +18 0.0258273785539499 +19 0.0213504646948306 +20 0.0175141661880584 +21 0.0142698211023571 +22 0.0115565159519833 +23 0.00930888979824423 +24 0.00746229206759215 +25 0.00595605679409682 +26 0.00473519993107751 +27 0.00375117728281841 +28 0.00296198928038098 +29 0.00233187862772459 +30 0.00183079868293457 +31 0.00143377454057296 +32 0.00107076258525208 +33 0.000773006742366448 +34 0.000539573690886396 +35 0.000364177599116743 +36 0.000237727628685579 +37 0.000150157714457011 +38 0.0000918283319498657 +39 0.0000544079947589854 +40 0.0000312548818921465 +41 0.0000174202619730274 +42 9.42698047424713e-6 +43 4.95614149002087e-6 +44 2.53275674485913e-6 +45 1.25854819834554e-6 +46 6.08116579596933e-7 +47 2.85572858589747e-7 +48 1.30129404249734e-7 +49 5.73280599448306e-8 +50 2.4219376577964e-8 +51 9.6316861194457e-9 +52 3.43804936850951e-9 +53 9.34806280366888e-10 +54 0.0 diff --git a/model/src/pyrenew/datasets/wastewater.py b/model/src/pyrenew/datasets/wastewater.py new file mode 100644 index 00000000..d7e09c6d --- /dev/null +++ b/model/src/pyrenew/datasets/wastewater.py @@ -0,0 +1,47 @@ +"""This module loads the package dataset named 'wastewater' and provides +functions to manipulate the data. It uses the 'polars' library""" + +from importlib.resources import files + +import polars as pl + + +def load_wastewater() -> pl.DataFrame: + """Load the wastewater dataset + + This dataset contains simulated entries of COVID-19 wastewater concentration + data. The dataset is used to demonstrate the use of the wastewater-informed + COVID-19 forecasting model. + + Returns + ------- + pl.DataFrame + The wastewater dataset + + Notes + ----- + + This dataset was downloaded directly from: + https://github.com/CDCgov/wastewater-informed-covid-forecasting/blob/292526383ece582f10823fc939c7e590ca349c6d/cfaforecastrenewalww/data/example_df.rda + + The dataset contains the following columns: + - `lab_wwtp_unique_id` + - `log_conc` + - `date` + - `lod_sewage` + - `below_LOD` + - `daily_hosp_admits` + - `daily_hosp_admits_for_eval` + - `pop` + - `forecast_date` + - `hosp_calibration_time` + - `site` + - `ww_pop` + - `inf_per_capita` + """ + + # Load the dataset + return pl.read_csv( + source=files("pyrenew.datasets") / "wastewater.tsv", + separator="\t", + ) diff --git a/model/src/pyrenew/datasets/wastewater.tsv b/model/src/pyrenew/datasets/wastewater.tsv new file mode 100644 index 00000000..3bc6ffb6 --- /dev/null +++ b/model/src/pyrenew/datasets/wastewater.tsv @@ -0,0 +1,636 @@ +t lab_wwtp_unique_id log_conc date lod_sewage below_LOD daily_hosp_admits daily_hosp_admits_for_eval pop forecast_date hosp_calibration_time site ww_pop inf_per_capita +1 1 2023-10-30 6 6 1e+06 2024-02-05 90 1 4e+05 0.000663426572737469 +1 2 2023-10-30 6 6 1e+06 2024-02-05 90 1 4e+05 0.000663426572737469 +1 3 2023-10-30 6 6 1e+06 2024-02-05 90 2 2e+05 0.000663426572737469 +1 4 2023-10-30 6 6 1e+06 2024-02-05 90 3 1e+05 0.000663426572737469 +1 5 2023-10-30 6 6 1e+06 2024-02-05 90 4 50000 0.000663426572737469 +2 1 2023-10-31 8 8 1e+06 2024-02-05 90 1 4e+05 0.00068300916315395 +2 2 2023-10-31 8 8 1e+06 2024-02-05 90 1 4e+05 0.00068300916315395 +2 3 8.36481035293498 2023-10-31 4.01075017260572 0 8 8 1e+06 2024-02-05 90 2 2e+05 0.00068300916315395 +2 4 2023-10-31 8 8 1e+06 2024-02-05 90 3 1e+05 0.00068300916315395 +2 5 2023-10-31 8 8 1e+06 2024-02-05 90 4 50000 0.00068300916315395 +3 1 2023-11-01 4 4 1e+06 2024-02-05 90 1 4e+05 0.000720865160318388 +3 2 2023-11-01 4 4 1e+06 2024-02-05 90 1 4e+05 0.000720865160318388 +3 3 2023-11-01 4 4 1e+06 2024-02-05 90 2 2e+05 0.000720865160318388 +3 4 2023-11-01 4 4 1e+06 2024-02-05 90 3 1e+05 0.000720865160318388 +3 5 2023-11-01 4 4 1e+06 2024-02-05 90 4 50000 0.000720865160318388 +4 1 8.06994711561617 2023-11-02 3.80549506716936 0 8 8 1e+06 2024-02-05 90 1 4e+05 0.00076104521535863 +4 2 2023-11-02 8 8 1e+06 2024-02-05 90 1 4e+05 0.00076104521535863 +4 3 2023-11-02 8 8 1e+06 2024-02-05 90 2 2e+05 0.00076104521535863 +4 4 2023-11-02 8 8 1e+06 2024-02-05 90 3 1e+05 0.00076104521535863 +4 5 2023-11-02 8 8 1e+06 2024-02-05 90 4 50000 0.00076104521535863 +5 1 2023-11-03 4 4 1e+06 2024-02-05 90 1 4e+05 0.000803295962116898 +5 2 2023-11-03 4 4 1e+06 2024-02-05 90 1 4e+05 0.000803295962116898 +5 3 2023-11-03 4 4 1e+06 2024-02-05 90 2 2e+05 0.000803295962116898 +5 4 2023-11-03 4 4 1e+06 2024-02-05 90 3 1e+05 0.000803295962116898 +5 5 2023-11-03 4 4 1e+06 2024-02-05 90 4 50000 0.000803295962116898 +6 1 2023-11-04 6 6 1e+06 2024-02-05 90 1 4e+05 0.000848573477258153 +6 2 2023-11-04 6 6 1e+06 2024-02-05 90 1 4e+05 0.000848573477258153 +6 3 2023-11-04 6 6 1e+06 2024-02-05 90 2 2e+05 0.000848573477258153 +6 4 2023-11-04 6 6 1e+06 2024-02-05 90 3 1e+05 0.000848573477258153 +6 5 2023-11-04 6 6 1e+06 2024-02-05 90 4 50000 0.000848573477258153 +7 1 2023-11-05 4 4 1e+06 2024-02-05 90 1 4e+05 0.000897005722285224 +7 2 2023-11-05 4 4 1e+06 2024-02-05 90 1 4e+05 0.000897005722285224 +7 3 2023-11-05 4 4 1e+06 2024-02-05 90 2 2e+05 0.000897005722285224 +7 4 2023-11-05 4 4 1e+06 2024-02-05 90 3 1e+05 0.000897005722285224 +7 5 2023-11-05 4 4 1e+06 2024-02-05 90 4 50000 0.000897005722285224 +8 1 2023-11-06 8 8 1e+06 2024-02-05 90 1 4e+05 0.000891349326433801 +8 2 2023-11-06 8 8 1e+06 2024-02-05 90 1 4e+05 0.000891349326433801 +8 3 2023-11-06 8 8 1e+06 2024-02-05 90 2 2e+05 0.000891349326433801 +8 4 2023-11-06 8 8 1e+06 2024-02-05 90 3 1e+05 0.000891349326433801 +8 5 2023-11-06 8 8 1e+06 2024-02-05 90 4 50000 0.000891349326433801 +9 1 2023-11-07 8 8 1e+06 2024-02-05 90 1 4e+05 0.000931363580626765 +9 2 2023-11-07 8 8 1e+06 2024-02-05 90 1 4e+05 0.000931363580626765 +9 3 2023-11-07 8 8 1e+06 2024-02-05 90 2 2e+05 0.000931363580626765 +9 4 2023-11-07 8 8 1e+06 2024-02-05 90 3 1e+05 0.000931363580626765 +9 5 7.30101935886591 2023-11-07 3.86712344199376 0 8 8 1e+06 2024-02-05 90 4 50000 0.000931363580626765 +10 1 2023-11-08 5 5 1e+06 2024-02-05 90 1 4e+05 0.000962241003942569 +10 2 8.62289790917382 2023-11-08 3.46396345552081 0 5 5 1e+06 2024-02-05 90 1 4e+05 0.000962241003942569 +10 3 2023-11-08 5 5 1e+06 2024-02-05 90 2 2e+05 0.000962241003942569 +10 4 2023-11-08 5 5 1e+06 2024-02-05 90 3 1e+05 0.000962241003942569 +10 5 2023-11-08 5 5 1e+06 2024-02-05 90 4 50000 0.000962241003942569 +11 1 2023-11-09 11 11 1e+06 2024-02-05 90 1 4e+05 0.000994116009532259 +11 2 2023-11-09 11 11 1e+06 2024-02-05 90 1 4e+05 0.000994116009532259 +11 3 2023-11-09 11 11 1e+06 2024-02-05 90 2 2e+05 0.000994116009532259 +11 4 2023-11-09 11 11 1e+06 2024-02-05 90 3 1e+05 0.000994116009532259 +11 5 2023-11-09 11 11 1e+06 2024-02-05 90 4 50000 0.000994116009532259 +12 1 8.63601927420382 2023-11-10 3.80549506716936 0 10 10 1e+06 2024-02-05 90 1 4e+05 0.00102753251025335 +12 2 2023-11-10 10 10 1e+06 2024-02-05 90 1 4e+05 0.00102753251025335 +12 3 2023-11-10 10 10 1e+06 2024-02-05 90 2 2e+05 0.00102753251025335 +12 4 2023-11-10 10 10 1e+06 2024-02-05 90 3 1e+05 0.00102753251025335 +12 5 2023-11-10 10 10 1e+06 2024-02-05 90 4 50000 0.00102753251025335 +13 1 8.7769160557112 2023-11-11 3.80549506716936 0 6 6 1e+06 2024-02-05 90 1 4e+05 0.00106205276321297 +13 2 2023-11-11 6 6 1e+06 2024-02-05 90 1 4e+05 0.00106205276321297 +13 3 8.07085011775403 2023-11-11 4.01075017260572 0 6 6 1e+06 2024-02-05 90 2 2e+05 0.00106205276321297 +13 4 2023-11-11 6 6 1e+06 2024-02-05 90 3 1e+05 0.00106205276321297 +13 5 2023-11-11 6 6 1e+06 2024-02-05 90 4 50000 0.00106205276321297 +14 1 2023-11-12 9 9 1e+06 2024-02-05 90 1 4e+05 0.00109775949551849 +14 2 2023-11-12 9 9 1e+06 2024-02-05 90 1 4e+05 0.00109775949551849 +14 3 2023-11-12 9 9 1e+06 2024-02-05 90 2 2e+05 0.00109775949551849 +14 4 2023-11-12 9 9 1e+06 2024-02-05 90 3 1e+05 0.00109775949551849 +14 5 2023-11-12 9 9 1e+06 2024-02-05 90 4 50000 0.00109775949551849 +15 1 2023-11-13 7 7 1e+06 2024-02-05 90 1 4e+05 0.00118548736868214 +15 2 2023-11-13 7 7 1e+06 2024-02-05 90 1 4e+05 0.00118548736868214 +15 3 2023-11-13 7 7 1e+06 2024-02-05 90 2 2e+05 0.00118548736868214 +15 4 2023-11-13 7 7 1e+06 2024-02-05 90 3 1e+05 0.00118548736868214 +15 5 7.56559007519856 2023-11-13 3.86712344199376 0 7 7 1e+06 2024-02-05 90 4 50000 0.00118548736868214 +16 1 2023-11-14 7 7 1e+06 2024-02-05 90 1 4e+05 0.00123489074410856 +16 2 9.14945816486916 2023-11-14 3.46396345552081 0 7 7 1e+06 2024-02-05 90 1 4e+05 0.00123489074410856 +16 3 7.79259216965746 2023-11-14 4.01075017260572 0 7 7 1e+06 2024-02-05 90 2 2e+05 0.00123489074410856 +16 4 2023-11-14 7 7 1e+06 2024-02-05 90 3 1e+05 0.00123489074410856 +16 5 10.9889396162503 2023-11-14 3.86712344199376 0 7 7 1e+06 2024-02-05 90 4 50000 0.00123489074410856 +17 1 2023-11-15 13 13 1e+06 2024-02-05 90 1 4e+05 0.0012971223308044 +17 2 2023-11-15 13 13 1e+06 2024-02-05 90 1 4e+05 0.0012971223308044 +17 3 8.08765357721804 2023-11-15 4.01075017260572 0 13 13 1e+06 2024-02-05 90 2 2e+05 0.0012971223308044 +17 4 2023-11-15 13 13 1e+06 2024-02-05 90 3 1e+05 0.0012971223308044 +17 5 2023-11-15 13 13 1e+06 2024-02-05 90 4 50000 0.0012971223308044 +18 1 8.51517036359834 2023-11-16 3.80549506716936 0 7 7 1e+06 2024-02-05 90 1 4e+05 0.00136244887607943 +18 2 2023-11-16 7 7 1e+06 2024-02-05 90 1 4e+05 0.00136244887607943 +18 3 2023-11-16 7 7 1e+06 2024-02-05 90 2 2e+05 0.00136244887607943 +18 4 2023-11-16 7 7 1e+06 2024-02-05 90 3 1e+05 0.00136244887607943 +18 5 9.93514217966295 2023-11-16 3.86712344199376 0 7 7 1e+06 2024-02-05 90 4 50000 0.00136244887607943 +19 1 2023-11-17 9 9 1e+06 2024-02-05 90 1 4e+05 0.00143067886419521 +19 2 2023-11-17 9 9 1e+06 2024-02-05 90 1 4e+05 0.00143067886419521 +19 3 2023-11-17 9 9 1e+06 2024-02-05 90 2 2e+05 0.00143067886419521 +19 4 2023-11-17 9 9 1e+06 2024-02-05 90 3 1e+05 0.00143067886419521 +19 5 2023-11-17 9 9 1e+06 2024-02-05 90 4 50000 0.00143067886419521 +20 1 2023-11-18 8 8 1e+06 2024-02-05 90 1 4e+05 0.0015024490166559 +20 2 2023-11-18 8 8 1e+06 2024-02-05 90 1 4e+05 0.0015024490166559 +20 3 2023-11-18 8 8 1e+06 2024-02-05 90 2 2e+05 0.0015024490166559 +20 4 2023-11-18 8 8 1e+06 2024-02-05 90 3 1e+05 0.0015024490166559 +20 5 2023-11-18 8 8 1e+06 2024-02-05 90 4 50000 0.0015024490166559 +21 1 8.91358538794599 2023-11-19 3.80549506716936 0 7 7 1e+06 2024-02-05 90 1 4e+05 0.001577887132097 +21 2 2023-11-19 7 7 1e+06 2024-02-05 90 1 4e+05 0.001577887132097 +21 3 2023-11-19 7 7 1e+06 2024-02-05 90 2 2e+05 0.001577887132097 +21 4 2023-11-19 7 7 1e+06 2024-02-05 90 3 1e+05 0.001577887132097 +21 5 2023-11-19 7 7 1e+06 2024-02-05 90 4 50000 0.001577887132097 +22 1 2023-11-20 7 7 1e+06 2024-02-05 90 1 4e+05 0.00157949362228245 +22 2 2023-11-20 7 7 1e+06 2024-02-05 90 1 4e+05 0.00157949362228245 +22 3 2023-11-20 7 7 1e+06 2024-02-05 90 2 2e+05 0.00157949362228245 +22 4 2023-11-20 7 7 1e+06 2024-02-05 90 3 1e+05 0.00157949362228245 +22 5 2023-11-20 7 7 1e+06 2024-02-05 90 4 50000 0.00157949362228245 +23 1 2023-11-21 11 11 1e+06 2024-02-05 90 1 4e+05 0.00164502885114781 +23 2 2023-11-21 11 11 1e+06 2024-02-05 90 1 4e+05 0.00164502885114781 +23 3 8.40842353113087 2023-11-21 4.01075017260572 0 11 11 1e+06 2024-02-05 90 2 2e+05 0.00164502885114781 +23 4 9.0135852111729 2023-11-21 3.57608017908556 0 11 11 1e+06 2024-02-05 90 3 1e+05 0.00164502885114781 +23 5 2023-11-21 11 11 1e+06 2024-02-05 90 4 50000 0.00164502885114781 +24 1 2023-11-22 17 17 1e+06 2024-02-05 90 1 4e+05 0.00169799949177993 +24 2 2023-11-22 17 17 1e+06 2024-02-05 90 1 4e+05 0.00169799949177993 +24 3 2023-11-22 17 17 1e+06 2024-02-05 90 2 2e+05 0.00169799949177993 +24 4 8.17923398023512 2023-11-22 3.57608017908556 0 17 17 1e+06 2024-02-05 90 3 1e+05 0.00169799949177993 +24 5 2023-11-22 17 17 1e+06 2024-02-05 90 4 50000 0.00169799949177993 +25 1 2023-11-23 5 5 1e+06 2024-02-05 90 1 4e+05 0.00175271129382294 +25 2 2023-11-23 5 5 1e+06 2024-02-05 90 1 4e+05 0.00175271129382294 +25 3 2023-11-23 5 5 1e+06 2024-02-05 90 2 2e+05 0.00175271129382294 +25 4 8.80015944143164 2023-11-23 3.57608017908556 0 5 5 1e+06 2024-02-05 90 3 1e+05 0.00175271129382294 +25 5 2023-11-23 5 5 1e+06 2024-02-05 90 4 50000 0.00175271129382294 +26 1 9.11537097218074 2023-11-24 3.80549506716936 0 14 14 1e+06 2024-02-05 90 1 4e+05 0.00180990448808453 +26 2 2023-11-24 14 14 1e+06 2024-02-05 90 1 4e+05 0.00180990448808453 +26 3 2023-11-24 14 14 1e+06 2024-02-05 90 2 2e+05 0.00180990448808453 +26 4 2023-11-24 14 14 1e+06 2024-02-05 90 3 1e+05 0.00180990448808453 +26 5 2023-11-24 14 14 1e+06 2024-02-05 90 4 50000 0.00180990448808453 +27 1 2023-11-25 9 9 1e+06 2024-02-05 90 1 4e+05 0.00186897576775536 +27 2 2023-11-25 9 9 1e+06 2024-02-05 90 1 4e+05 0.00186897576775536 +27 3 8.98637293859266 2023-11-25 4.01075017260572 0 9 9 1e+06 2024-02-05 90 2 2e+05 0.00186897576775536 +27 4 2023-11-25 9 9 1e+06 2024-02-05 90 3 1e+05 0.00186897576775536 +27 5 2023-11-25 9 9 1e+06 2024-02-05 90 4 50000 0.00186897576775536 +28 1 2023-11-26 10 10 1e+06 2024-02-05 90 1 4e+05 0.0019300544630776 +28 2 2023-11-26 10 10 1e+06 2024-02-05 90 1 4e+05 0.0019300544630776 +28 3 8.65889243616303 2023-11-26 4.01075017260572 0 10 10 1e+06 2024-02-05 90 2 2e+05 0.0019300544630776 +28 4 2023-11-26 10 10 1e+06 2024-02-05 90 3 1e+05 0.0019300544630776 +28 5 2023-11-26 10 10 1e+06 2024-02-05 90 4 50000 0.0019300544630776 +29 1 2023-11-27 16 16 1e+06 2024-02-05 90 1 4e+05 0.00205882259115435 +29 2 2023-11-27 16 16 1e+06 2024-02-05 90 1 4e+05 0.00205882259115435 +29 3 2023-11-27 16 16 1e+06 2024-02-05 90 2 2e+05 0.00205882259115435 +29 4 2023-11-27 16 16 1e+06 2024-02-05 90 3 1e+05 0.00205882259115435 +29 5 2023-11-27 16 16 1e+06 2024-02-05 90 4 50000 0.00205882259115435 +30 1 2023-11-28 11 11 1e+06 2024-02-05 90 1 4e+05 0.00213837716378528 +30 2 2023-11-28 11 11 1e+06 2024-02-05 90 1 4e+05 0.00213837716378528 +30 3 2023-11-28 11 11 1e+06 2024-02-05 90 2 2e+05 0.00213837716378528 +30 4 2023-11-28 11 11 1e+06 2024-02-05 90 3 1e+05 0.00213837716378528 +30 5 2023-11-28 11 11 1e+06 2024-02-05 90 4 50000 0.00213837716378528 +31 1 2023-11-29 15 15 1e+06 2024-02-05 90 1 4e+05 0.0022349415558482 +31 2 2023-11-29 15 15 1e+06 2024-02-05 90 1 4e+05 0.0022349415558482 +31 3 2023-11-29 15 15 1e+06 2024-02-05 90 2 2e+05 0.0022349415558482 +31 4 2023-11-29 15 15 1e+06 2024-02-05 90 3 1e+05 0.0022349415558482 +31 5 10.3643033355416 2023-11-29 3.86712344199376 0 15 15 1e+06 2024-02-05 90 4 50000 0.0022349415558482 +32 1 2023-11-30 19 19 1e+06 2024-02-05 90 1 4e+05 0.00233594266690393 +32 2 2023-11-30 19 19 1e+06 2024-02-05 90 1 4e+05 0.00233594266690393 +32 3 2023-11-30 19 19 1e+06 2024-02-05 90 2 2e+05 0.00233594266690393 +32 4 2023-11-30 19 19 1e+06 2024-02-05 90 3 1e+05 0.00233594266690393 +32 5 2023-11-30 19 19 1e+06 2024-02-05 90 4 50000 0.00233594266690393 +33 1 2023-12-01 16 16 1e+06 2024-02-05 90 1 4e+05 0.00244114199543507 +33 2 2023-12-01 16 16 1e+06 2024-02-05 90 1 4e+05 0.00244114199543507 +33 3 2023-12-01 16 16 1e+06 2024-02-05 90 2 2e+05 0.00244114199543507 +33 4 2023-12-01 16 16 1e+06 2024-02-05 90 3 1e+05 0.00244114199543507 +33 5 9.06315348681911 2023-12-01 3.86712344199376 0 16 16 1e+06 2024-02-05 90 4 50000 0.00244114199543507 +34 1 2023-12-02 17 17 1e+06 2024-02-05 90 1 4e+05 0.00255137419422608 +34 2 2023-12-02 17 17 1e+06 2024-02-05 90 1 4e+05 0.00255137419422608 +34 3 2023-12-02 17 17 1e+06 2024-02-05 90 2 2e+05 0.00255137419422608 +34 4 2023-12-02 17 17 1e+06 2024-02-05 90 3 1e+05 0.00255137419422608 +34 5 2023-12-02 17 17 1e+06 2024-02-05 90 4 50000 0.00255137419422608 +35 1 2023-12-03 15 15 1e+06 2024-02-05 90 1 4e+05 0.00266681808967801 +35 2 9.62773574720005 2023-12-03 3.46396345552081 0 15 15 1e+06 2024-02-05 90 1 4e+05 0.00266681808967801 +35 3 8.98340694902573 2023-12-03 4.01075017260572 0 15 15 1e+06 2024-02-05 90 2 2e+05 0.00266681808967801 +35 4 2023-12-03 15 15 1e+06 2024-02-05 90 3 1e+05 0.00266681808967801 +35 5 8.68576411206245 2023-12-03 3.86712344199376 0 15 15 1e+06 2024-02-05 90 4 50000 0.00266681808967801 +36 1 2023-12-04 14 14 1e+06 2024-02-05 90 1 4e+05 0.00229463001235797 +36 2 2023-12-04 14 14 1e+06 2024-02-05 90 1 4e+05 0.00229463001235797 +36 3 2023-12-04 14 14 1e+06 2024-02-05 90 2 2e+05 0.00229463001235797 +36 4 8.82036499084619 2023-12-04 3.57608017908556 0 14 14 1e+06 2024-02-05 90 3 1e+05 0.00229463001235797 +36 5 2023-12-04 14 14 1e+06 2024-02-05 90 4 50000 0.00229463001235797 +37 1 9.74824268632507 2023-12-05 3.80549506716936 0 22 22 1e+06 2024-02-05 90 1 4e+05 0.00232400753101634 +37 2 2023-12-05 22 22 1e+06 2024-02-05 90 1 4e+05 0.00232400753101634 +37 3 2023-12-05 22 22 1e+06 2024-02-05 90 2 2e+05 0.00232400753101634 +37 4 2023-12-05 22 22 1e+06 2024-02-05 90 3 1e+05 0.00232400753101634 +37 5 2023-12-05 22 22 1e+06 2024-02-05 90 4 50000 0.00232400753101634 +38 1 9.38915688760292 2023-12-06 3.80549506716936 0 25 25 1e+06 2024-02-05 90 1 4e+05 0.00227031881582818 +38 2 9.77665706884518 2023-12-06 3.46396345552081 0 25 25 1e+06 2024-02-05 90 1 4e+05 0.00227031881582818 +38 3 2023-12-06 25 25 1e+06 2024-02-05 90 2 2e+05 0.00227031881582818 +38 4 2023-12-06 25 25 1e+06 2024-02-05 90 3 1e+05 0.00227031881582818 +38 5 2023-12-06 25 25 1e+06 2024-02-05 90 4 50000 0.00227031881582818 +39 1 2023-12-07 15 15 1e+06 2024-02-05 90 1 4e+05 0.00221548429304801 +39 2 2023-12-07 15 15 1e+06 2024-02-05 90 1 4e+05 0.00221548429304801 +39 3 2023-12-07 15 15 1e+06 2024-02-05 90 2 2e+05 0.00221548429304801 +39 4 2023-12-07 15 15 1e+06 2024-02-05 90 3 1e+05 0.00221548429304801 +39 5 9.68167643822757 2023-12-07 3.86712344199376 0 15 15 1e+06 2024-02-05 90 4 50000 0.00221548429304801 +40 1 2023-12-08 21 21 1e+06 2024-02-05 90 1 4e+05 0.00216555268800312 +40 2 2023-12-08 21 21 1e+06 2024-02-05 90 1 4e+05 0.00216555268800312 +40 3 2023-12-08 21 21 1e+06 2024-02-05 90 2 2e+05 0.00216555268800312 +40 4 2023-12-08 21 21 1e+06 2024-02-05 90 3 1e+05 0.00216555268800312 +40 5 2023-12-08 21 21 1e+06 2024-02-05 90 4 50000 0.00216555268800312 +41 1 2023-12-09 30 30 1e+06 2024-02-05 90 1 4e+05 0.00211674308322386 +41 2 2023-12-09 30 30 1e+06 2024-02-05 90 1 4e+05 0.00211674308322386 +41 3 2023-12-09 30 30 1e+06 2024-02-05 90 2 2e+05 0.00211674308322386 +41 4 2023-12-09 30 30 1e+06 2024-02-05 90 3 1e+05 0.00211674308322386 +41 5 2023-12-09 30 30 1e+06 2024-02-05 90 4 50000 0.00211674308322386 +42 1 2023-12-10 19 19 1e+06 2024-02-05 90 1 4e+05 0.00206916837552496 +42 2 2023-12-10 19 19 1e+06 2024-02-05 90 1 4e+05 0.00206916837552496 +42 3 2023-12-10 19 19 1e+06 2024-02-05 90 2 2e+05 0.00206916837552496 +42 4 2023-12-10 19 19 1e+06 2024-02-05 90 3 1e+05 0.00206916837552496 +42 5 6.64193078483244 2023-12-10 3.86712344199376 0 19 19 1e+06 2024-02-05 90 4 50000 0.00206916837552496 +43 1 2023-12-11 14 14 1e+06 2024-02-05 90 1 4e+05 0.00198425201386512 +43 2 2023-12-11 14 14 1e+06 2024-02-05 90 1 4e+05 0.00198425201386512 +43 3 2023-12-11 14 14 1e+06 2024-02-05 90 2 2e+05 0.00198425201386512 +43 4 2023-12-11 14 14 1e+06 2024-02-05 90 3 1e+05 0.00198425201386512 +43 5 2023-12-11 14 14 1e+06 2024-02-05 90 4 50000 0.00198425201386512 +44 1 2023-12-12 14 14 1e+06 2024-02-05 90 1 4e+05 0.00193420137020504 +44 2 2023-12-12 14 14 1e+06 2024-02-05 90 1 4e+05 0.00193420137020504 +44 3 2023-12-12 14 14 1e+06 2024-02-05 90 2 2e+05 0.00193420137020504 +44 4 2023-12-12 14 14 1e+06 2024-02-05 90 3 1e+05 0.00193420137020504 +44 5 2023-12-12 14 14 1e+06 2024-02-05 90 4 50000 0.00193420137020504 +45 1 2023-12-13 21 21 1e+06 2024-02-05 90 1 4e+05 0.00187928821226127 +45 2 2023-12-13 21 21 1e+06 2024-02-05 90 1 4e+05 0.00187928821226127 +45 3 2023-12-13 21 21 1e+06 2024-02-05 90 2 2e+05 0.00187928821226127 +45 4 2023-12-13 21 21 1e+06 2024-02-05 90 3 1e+05 0.00187928821226127 +45 5 7.99235103936222 2023-12-13 3.86712344199376 0 21 21 1e+06 2024-02-05 90 4 50000 0.00187928821226127 +46 1 2023-12-14 18 18 1e+06 2024-02-05 90 1 4e+05 0.00182588119935607 +46 2 2023-12-14 18 18 1e+06 2024-02-05 90 1 4e+05 0.00182588119935607 +46 3 2023-12-14 18 18 1e+06 2024-02-05 90 2 2e+05 0.00182588119935607 +46 4 2023-12-14 18 18 1e+06 2024-02-05 90 3 1e+05 0.00182588119935607 +46 5 10.4835008686585 2023-12-14 3.86712344199376 0 18 18 1e+06 2024-02-05 90 4 50000 0.00182588119935607 +47 1 2023-12-15 21 21 1e+06 2024-02-05 90 1 4e+05 0.00177445076407144 +47 2 2023-12-15 21 21 1e+06 2024-02-05 90 1 4e+05 0.00177445076407144 +47 3 2023-12-15 21 21 1e+06 2024-02-05 90 2 2e+05 0.00177445076407144 +47 4 2023-12-15 21 21 1e+06 2024-02-05 90 3 1e+05 0.00177445076407144 +47 5 2023-12-15 21 21 1e+06 2024-02-05 90 4 50000 0.00177445076407144 +48 1 2023-12-16 17 17 1e+06 2024-02-05 90 1 4e+05 0.0017246593197175 +48 2 2023-12-16 17 17 1e+06 2024-02-05 90 1 4e+05 0.0017246593197175 +48 3 2023-12-16 17 17 1e+06 2024-02-05 90 2 2e+05 0.0017246593197175 +48 4 2023-12-16 17 17 1e+06 2024-02-05 90 3 1e+05 0.0017246593197175 +48 5 7.77021794286154 2023-12-16 3.86712344199376 0 17 17 1e+06 2024-02-05 90 4 50000 0.0017246593197175 +49 1 2023-12-17 17 17 1e+06 2024-02-05 90 1 4e+05 0.00167645012593016 +49 2 2023-12-17 17 17 1e+06 2024-02-05 90 1 4e+05 0.00167645012593016 +49 3 2023-12-17 17 17 1e+06 2024-02-05 90 2 2e+05 0.00167645012593016 +49 4 2023-12-17 17 17 1e+06 2024-02-05 90 3 1e+05 0.00167645012593016 +49 5 2023-12-17 17 17 1e+06 2024-02-05 90 4 50000 0.00167645012593016 +50 1 2023-12-18 20 20 1e+06 2024-02-05 90 1 4e+05 0.00158824729708959 +50 2 2023-12-18 20 20 1e+06 2024-02-05 90 1 4e+05 0.00158824729708959 +50 3 2023-12-18 20 20 1e+06 2024-02-05 90 2 2e+05 0.00158824729708959 +50 4 2023-12-18 20 20 1e+06 2024-02-05 90 3 1e+05 0.00158824729708959 +50 5 2023-12-18 20 20 1e+06 2024-02-05 90 4 50000 0.00158824729708959 +51 1 2023-12-19 25 25 1e+06 2024-02-05 90 1 4e+05 0.00153779928056313 +51 2 2023-12-19 25 25 1e+06 2024-02-05 90 1 4e+05 0.00153779928056313 +51 3 2023-12-19 25 25 1e+06 2024-02-05 90 2 2e+05 0.00153779928056313 +51 4 2023-12-19 25 25 1e+06 2024-02-05 90 3 1e+05 0.00153779928056313 +51 5 7.35970606877435 2023-12-19 3.86712344199376 0 25 25 1e+06 2024-02-05 90 4 50000 0.00153779928056313 +52 1 2023-12-20 23 23 1e+06 2024-02-05 90 1 4e+05 0.00148313981255738 +52 2 2023-12-20 23 23 1e+06 2024-02-05 90 1 4e+05 0.00148313981255738 +52 3 2023-12-20 23 23 1e+06 2024-02-05 90 2 2e+05 0.00148313981255738 +52 4 2023-12-20 23 23 1e+06 2024-02-05 90 3 1e+05 0.00148313981255738 +52 5 2023-12-20 23 23 1e+06 2024-02-05 90 4 50000 0.00148313981255738 +53 1 2023-12-21 15 15 1e+06 2024-02-05 90 1 4e+05 0.00143040338345784 +53 2 2023-12-21 15 15 1e+06 2024-02-05 90 1 4e+05 0.00143040338345784 +53 3 2023-12-21 15 15 1e+06 2024-02-05 90 2 2e+05 0.00143040338345784 +53 4 2023-12-21 15 15 1e+06 2024-02-05 90 3 1e+05 0.00143040338345784 +53 5 4.94160988713828 2023-12-21 3.86712344199376 0 15 15 1e+06 2024-02-05 90 4 50000 0.00143040338345784 +54 1 2023-12-22 19 19 1e+06 2024-02-05 90 1 4e+05 0.00138017763674767 +54 2 2023-12-22 19 19 1e+06 2024-02-05 90 1 4e+05 0.00138017763674767 +54 3 2023-12-22 19 19 1e+06 2024-02-05 90 2 2e+05 0.00138017763674767 +54 4 2023-12-22 19 19 1e+06 2024-02-05 90 3 1e+05 0.00138017763674767 +54 5 2023-12-22 19 19 1e+06 2024-02-05 90 4 50000 0.00138017763674767 +55 1 2023-12-23 14 14 1e+06 2024-02-05 90 1 4e+05 0.00133208879651221 +55 2 2023-12-23 14 14 1e+06 2024-02-05 90 1 4e+05 0.00133208879651221 +55 3 2023-12-23 14 14 1e+06 2024-02-05 90 2 2e+05 0.00133208879651221 +55 4 2023-12-23 14 14 1e+06 2024-02-05 90 3 1e+05 0.00133208879651221 +55 5 2023-12-23 14 14 1e+06 2024-02-05 90 4 50000 0.00133208879651221 +56 1 2023-12-24 15 15 1e+06 2024-02-05 90 1 4e+05 0.00128602958318709 +56 2 2023-12-24 15 15 1e+06 2024-02-05 90 1 4e+05 0.00128602958318709 +56 3 8.63762486215468 2023-12-24 4.01075017260572 0 15 15 1e+06 2024-02-05 90 2 2e+05 0.00128602958318709 +56 4 2023-12-24 15 15 1e+06 2024-02-05 90 3 1e+05 0.00128602958318709 +56 5 2023-12-24 15 15 1e+06 2024-02-05 90 4 50000 0.00128602958318709 +57 1 2023-12-25 18 18 1e+06 2024-02-05 90 1 4e+05 0.00127286824470429 +57 2 2023-12-25 18 18 1e+06 2024-02-05 90 1 4e+05 0.00127286824470429 +57 3 2023-12-25 18 18 1e+06 2024-02-05 90 2 2e+05 0.00127286824470429 +57 4 2023-12-25 18 18 1e+06 2024-02-05 90 3 1e+05 0.00127286824470429 +57 5 7.76148032124528 2023-12-25 3.86712344199376 0 18 18 1e+06 2024-02-05 90 4 50000 0.00127286824470429 +58 1 2023-12-26 17 17 1e+06 2024-02-05 90 1 4e+05 0.00123306532431121 +58 2 2023-12-26 17 17 1e+06 2024-02-05 90 1 4e+05 0.00123306532431121 +58 3 9.05510064379279 2023-12-26 4.01075017260572 0 17 17 1e+06 2024-02-05 90 2 2e+05 0.00123306532431121 +58 4 2023-12-26 17 17 1e+06 2024-02-05 90 3 1e+05 0.00123306532431121 +58 5 2023-12-26 17 17 1e+06 2024-02-05 90 4 50000 0.00123306532431121 +59 1 2023-12-27 15 15 1e+06 2024-02-05 90 1 4e+05 0.00120113303290099 +59 2 2023-12-27 15 15 1e+06 2024-02-05 90 1 4e+05 0.00120113303290099 +59 3 2023-12-27 15 15 1e+06 2024-02-05 90 2 2e+05 0.00120113303290099 +59 4 2023-12-27 15 15 1e+06 2024-02-05 90 3 1e+05 0.00120113303290099 +59 5 5.87429657105247 2023-12-27 3.86712344199376 0 15 15 1e+06 2024-02-05 90 4 50000 0.00120113303290099 +60 1 2023-12-28 15 15 1e+06 2024-02-05 90 1 4e+05 0.00117034513242189 +60 2 2023-12-28 15 15 1e+06 2024-02-05 90 1 4e+05 0.00117034513242189 +60 3 8.76052661905076 2023-12-28 4.01075017260572 0 15 15 1e+06 2024-02-05 90 2 2e+05 0.00117034513242189 +60 4 2023-12-28 15 15 1e+06 2024-02-05 90 3 1e+05 0.00117034513242189 +60 5 2023-12-28 15 15 1e+06 2024-02-05 90 4 50000 0.00117034513242189 +61 1 2023-12-29 16 16 1e+06 2024-02-05 90 1 4e+05 0.00114037114655214 +61 2 2023-12-29 16 16 1e+06 2024-02-05 90 1 4e+05 0.00114037114655214 +61 3 2023-12-29 16 16 1e+06 2024-02-05 90 2 2e+05 0.00114037114655214 +61 4 2023-12-29 16 16 1e+06 2024-02-05 90 3 1e+05 0.00114037114655214 +61 5 2023-12-29 16 16 1e+06 2024-02-05 90 4 50000 0.00114037114655214 +62 1 2023-12-30 14 14 1e+06 2024-02-05 90 1 4e+05 0.00111147947392433 +62 2 2023-12-30 14 14 1e+06 2024-02-05 90 1 4e+05 0.00111147947392433 +62 3 8.78382910218133 2023-12-30 4.01075017260572 0 14 14 1e+06 2024-02-05 90 2 2e+05 0.00111147947392433 +62 4 2023-12-30 14 14 1e+06 2024-02-05 90 3 1e+05 0.00111147947392433 +62 5 2023-12-30 14 14 1e+06 2024-02-05 90 4 50000 0.00111147947392433 +63 1 2023-12-31 21 21 1e+06 2024-02-05 90 1 4e+05 0.00108359782801036 +63 2 2023-12-31 21 21 1e+06 2024-02-05 90 1 4e+05 0.00108359782801036 +63 3 2023-12-31 21 21 1e+06 2024-02-05 90 2 2e+05 0.00108359782801036 +63 4 2023-12-31 21 21 1e+06 2024-02-05 90 3 1e+05 0.00108359782801036 +63 5 7.38108580547644 2023-12-31 3.86712344199376 0 21 21 1e+06 2024-02-05 90 4 50000 0.00108359782801036 +64 1 2024-01-01 15 15 1e+06 2024-02-05 90 1 4e+05 0.000958544149071965 +64 2 2024-01-01 15 15 1e+06 2024-02-05 90 1 4e+05 0.000958544149071965 +64 3 2024-01-01 15 15 1e+06 2024-02-05 90 2 2e+05 0.000958544149071965 +64 4 2024-01-01 15 15 1e+06 2024-02-05 90 3 1e+05 0.000958544149071965 +64 5 2024-01-01 15 15 1e+06 2024-02-05 90 4 50000 0.000958544149071965 +65 1 8.82362479595895 2024-01-02 3.80549506716936 0 14 14 1e+06 2024-02-05 90 1 4e+05 0.000920852323403138 +65 2 2024-01-02 14 14 1e+06 2024-02-05 90 1 4e+05 0.000920852323403138 +65 3 2024-01-02 14 14 1e+06 2024-02-05 90 2 2e+05 0.000920852323403138 +65 4 2024-01-02 14 14 1e+06 2024-02-05 90 3 1e+05 0.000920852323403138 +65 5 5.90928660433236 2024-01-02 3.86712344199376 0 14 14 1e+06 2024-02-05 90 4 50000 0.000920852323403138 +66 1 2024-01-03 16 16 1e+06 2024-02-05 90 1 4e+05 0.000869425096381393 +66 2 2024-01-03 16 16 1e+06 2024-02-05 90 1 4e+05 0.000869425096381393 +66 3 2024-01-03 16 16 1e+06 2024-02-05 90 2 2e+05 0.000869425096381393 +66 4 2024-01-03 16 16 1e+06 2024-02-05 90 3 1e+05 0.000869425096381393 +66 5 2024-01-03 16 16 1e+06 2024-02-05 90 4 50000 0.000869425096381393 +67 1 2024-01-04 14 14 1e+06 2024-02-05 90 1 4e+05 0.00082017174434338 +67 2 2024-01-04 14 14 1e+06 2024-02-05 90 1 4e+05 0.00082017174434338 +67 3 2024-01-04 14 14 1e+06 2024-02-05 90 2 2e+05 0.00082017174434338 +67 4 2024-01-04 14 14 1e+06 2024-02-05 90 3 1e+05 0.00082017174434338 +67 5 2024-01-04 14 14 1e+06 2024-02-05 90 4 50000 0.00082017174434338 +68 1 2024-01-05 7 7 1e+06 2024-02-05 90 1 4e+05 0.000774375898920149 +68 2 2024-01-05 7 7 1e+06 2024-02-05 90 1 4e+05 0.000774375898920149 +68 3 2024-01-05 7 7 1e+06 2024-02-05 90 2 2e+05 0.000774375898920149 +68 4 2024-01-05 7 7 1e+06 2024-02-05 90 3 1e+05 0.000774375898920149 +68 5 2024-01-05 7 7 1e+06 2024-02-05 90 4 50000 0.000774375898920149 +69 1 2024-01-06 12 12 1e+06 2024-02-05 90 1 4e+05 0.000731171823371798 +69 2 2024-01-06 12 12 1e+06 2024-02-05 90 1 4e+05 0.000731171823371798 +69 3 2024-01-06 12 12 1e+06 2024-02-05 90 2 2e+05 0.000731171823371798 +69 4 2024-01-06 12 12 1e+06 2024-02-05 90 3 1e+05 0.000731171823371798 +69 5 6.47983174817185 2024-01-06 3.86712344199376 0 12 12 1e+06 2024-02-05 90 4 50000 0.000731171823371798 +70 1 2024-01-07 16 16 1e+06 2024-02-05 90 1 4e+05 0.000690411943226136 +70 2 2024-01-07 16 16 1e+06 2024-02-05 90 1 4e+05 0.000690411943226136 +70 3 2024-01-07 16 16 1e+06 2024-02-05 90 2 2e+05 0.000690411943226136 +70 4 2024-01-07 16 16 1e+06 2024-02-05 90 3 1e+05 0.000690411943226136 +70 5 6.38750048014559 2024-01-07 3.86712344199376 0 16 16 1e+06 2024-02-05 90 4 50000 0.000690411943226136 +71 1 2024-01-08 15 15 1e+06 2024-02-05 90 1 4e+05 0.000823597224512825 +71 2 2024-01-08 15 15 1e+06 2024-02-05 90 1 4e+05 0.000823597224512825 +71 3 2024-01-08 15 15 1e+06 2024-02-05 90 2 2e+05 0.000823597224512825 +71 4 2024-01-08 15 15 1e+06 2024-02-05 90 3 1e+05 0.000823597224512825 +71 5 2024-01-08 15 15 1e+06 2024-02-05 90 4 50000 0.000823597224512825 +72 1 8.46414793720768 2024-01-09 3.80549506716936 0 14 14 1e+06 2024-02-05 90 1 4e+05 0.00080748258927143 +72 2 2024-01-09 14 14 1e+06 2024-02-05 90 1 4e+05 0.00080748258927143 +72 3 2024-01-09 14 14 1e+06 2024-02-05 90 2 2e+05 0.00080748258927143 +72 4 2024-01-09 14 14 1e+06 2024-02-05 90 3 1e+05 0.00080748258927143 +72 5 2024-01-09 14 14 1e+06 2024-02-05 90 4 50000 0.00080748258927143 +73 1 2024-01-10 14 14 1e+06 2024-02-05 90 1 4e+05 0.000826618103107369 +73 2 2024-01-10 14 14 1e+06 2024-02-05 90 1 4e+05 0.000826618103107369 +73 3 2024-01-10 14 14 1e+06 2024-02-05 90 2 2e+05 0.000826618103107369 +73 4 2024-01-10 14 14 1e+06 2024-02-05 90 3 1e+05 0.000826618103107369 +73 5 2024-01-10 14 14 1e+06 2024-02-05 90 4 50000 0.000826618103107369 +74 1 2024-01-11 11 11 1e+06 2024-02-05 90 1 4e+05 0.000846403532211776 +74 2 2024-01-11 11 11 1e+06 2024-02-05 90 1 4e+05 0.000846403532211776 +74 3 8.55794315455975 2024-01-11 4.01075017260572 0 11 11 1e+06 2024-02-05 90 2 2e+05 0.000846403532211776 +74 4 2024-01-11 11 11 1e+06 2024-02-05 90 3 1e+05 0.000846403532211776 +74 5 5.7255505343884 2024-01-11 3.86712344199376 0 11 11 1e+06 2024-02-05 90 4 50000 0.000846403532211776 +75 1 2024-01-12 8 8 1e+06 2024-02-05 90 1 4e+05 0.000865255577508535 +75 2 2024-01-12 8 8 1e+06 2024-02-05 90 1 4e+05 0.000865255577508535 +75 3 8.76014735348211 2024-01-12 4.01075017260572 0 8 8 1e+06 2024-02-05 90 2 2e+05 0.000865255577508535 +75 4 7.61019504798106 2024-01-12 3.57608017908556 0 8 8 1e+06 2024-02-05 90 3 1e+05 0.000865255577508535 +75 5 2024-01-12 8 8 1e+06 2024-02-05 90 4 50000 0.000865255577508535 +76 1 2024-01-13 17 17 1e+06 2024-02-05 90 1 4e+05 0.000884800635366129 +76 2 2024-01-13 17 17 1e+06 2024-02-05 90 1 4e+05 0.000884800635366129 +76 3 7.56863796591067 2024-01-13 4.01075017260572 0 17 17 1e+06 2024-02-05 90 2 2e+05 0.000884800635366129 +76 4 2024-01-13 17 17 1e+06 2024-02-05 90 3 1e+05 0.000884800635366129 +76 5 2024-01-13 17 17 1e+06 2024-02-05 90 4 50000 0.000884800635366129 +77 1 8.73956950752577 2024-01-14 3.80549506716936 0 10 10 1e+06 2024-02-05 90 1 4e+05 0.000904924166613176 +77 2 2024-01-14 10 10 1e+06 2024-02-05 90 1 4e+05 0.000904924166613176 +77 3 2024-01-14 10 10 1e+06 2024-02-05 90 2 2e+05 0.000904924166613176 +77 4 2024-01-14 10 10 1e+06 2024-02-05 90 3 1e+05 0.000904924166613176 +77 5 2024-01-14 10 10 1e+06 2024-02-05 90 4 50000 0.000904924166613176 +78 1 2024-01-15 11 11 1e+06 2024-02-05 90 1 4e+05 0.000861799654799829 +78 2 2024-01-15 11 11 1e+06 2024-02-05 90 1 4e+05 0.000861799654799829 +78 3 8.58159558797247 2024-01-15 4.01075017260572 0 11 11 1e+06 2024-02-05 90 2 2e+05 0.000861799654799829 +78 4 2024-01-15 11 11 1e+06 2024-02-05 90 3 1e+05 0.000861799654799829 +78 5 2024-01-15 11 11 1e+06 2024-02-05 90 4 50000 0.000861799654799829 +79 1 2024-01-16 15 15 1e+06 2024-02-05 90 1 4e+05 0.000871209204542765 +79 2 2024-01-16 15 15 1e+06 2024-02-05 90 1 4e+05 0.000871209204542765 +79 3 2024-01-16 15 15 1e+06 2024-02-05 90 2 2e+05 0.000871209204542765 +79 4 2024-01-16 15 15 1e+06 2024-02-05 90 3 1e+05 0.000871209204542765 +79 5 2024-01-16 15 15 1e+06 2024-02-05 90 4 50000 0.000871209204542765 +80 1 2024-01-17 14 14 1e+06 2024-02-05 90 1 4e+05 0.000869185211053884 +80 2 2024-01-17 14 14 1e+06 2024-02-05 90 1 4e+05 0.000869185211053884 +80 3 2024-01-17 14 14 1e+06 2024-02-05 90 2 2e+05 0.000869185211053884 +80 4 2024-01-17 14 14 1e+06 2024-02-05 90 3 1e+05 0.000869185211053884 +80 5 2024-01-17 14 14 1e+06 2024-02-05 90 4 50000 0.000869185211053884 +81 1 2024-01-18 8 8 1e+06 2024-02-05 90 1 4e+05 0.000867034121663692 +81 2 2024-01-18 8 8 1e+06 2024-02-05 90 1 4e+05 0.000867034121663692 +81 3 2024-01-18 8 8 1e+06 2024-02-05 90 2 2e+05 0.000867034121663692 +81 4 2024-01-18 8 8 1e+06 2024-02-05 90 3 1e+05 0.000867034121663692 +81 5 7.87700698821362 2024-01-18 3.86712344199376 0 8 8 1e+06 2024-02-05 90 4 50000 0.000867034121663692 +82 1 2024-01-19 8 8 1e+06 2024-02-05 90 1 4e+05 0.000865467129392603 +82 2 2024-01-19 8 8 1e+06 2024-02-05 90 1 4e+05 0.000865467129392603 +82 3 2024-01-19 8 8 1e+06 2024-02-05 90 2 2e+05 0.000865467129392603 +82 4 2024-01-19 8 8 1e+06 2024-02-05 90 3 1e+05 0.000865467129392603 +82 5 6.44736463346304 2024-01-19 3.86712344199376 0 8 8 1e+06 2024-02-05 90 4 50000 0.000865467129392603 +83 1 2024-01-20 5 5 1e+06 2024-02-05 90 1 4e+05 0.000863956372546081 +83 2 8.78221834364671 2024-01-20 3.46396345552081 0 5 5 1e+06 2024-02-05 90 1 4e+05 0.000863956372546081 +83 3 2024-01-20 5 5 1e+06 2024-02-05 90 2 2e+05 0.000863956372546081 +83 4 2024-01-20 5 5 1e+06 2024-02-05 90 3 1e+05 0.000863956372546081 +83 5 2024-01-20 5 5 1e+06 2024-02-05 90 4 50000 0.000863956372546081 +84 1 8.73564069126942 2024-01-21 3.80549506716936 0 7 7 1e+06 2024-02-05 90 1 4e+05 0.000862533776479147 +84 2 2024-01-21 7 7 1e+06 2024-02-05 90 1 4e+05 0.000862533776479147 +84 3 2024-01-21 7 7 1e+06 2024-02-05 90 2 2e+05 0.000862533776479147 +84 4 2024-01-21 7 7 1e+06 2024-02-05 90 3 1e+05 0.000862533776479147 +84 5 2024-01-21 7 7 1e+06 2024-02-05 90 4 50000 0.000862533776479147 +85 1 8.1183762055312 2024-01-22 3.80549506716936 0 9 9 1e+06 2024-02-05 90 1 4e+05 0.000828059442079499 +85 2 2024-01-22 9 9 1e+06 2024-02-05 90 1 4e+05 0.000828059442079499 +85 3 8.15121391750941 2024-01-22 4.01075017260572 0 9 9 1e+06 2024-02-05 90 2 2e+05 0.000828059442079499 +85 4 2024-01-22 9 9 1e+06 2024-02-05 90 3 1e+05 0.000828059442079499 +85 5 2024-01-22 9 9 1e+06 2024-02-05 90 4 50000 0.000828059442079499 +86 1 2024-01-23 12 12 1e+06 2024-02-05 90 1 4e+05 0.000821602526519418 +86 2 2024-01-23 12 12 1e+06 2024-02-05 90 1 4e+05 0.000821602526519418 +86 3 2024-01-23 12 12 1e+06 2024-02-05 90 2 2e+05 0.000821602526519418 +86 4 2024-01-23 12 12 1e+06 2024-02-05 90 3 1e+05 0.000821602526519418 +86 5 2024-01-23 12 12 1e+06 2024-02-05 90 4 50000 0.000821602526519418 +87 1 2024-01-24 6 6 1e+06 2024-02-05 90 1 4e+05 0.000809336965969484 +87 2 2024-01-24 6 6 1e+06 2024-02-05 90 1 4e+05 0.000809336965969484 +87 3 2024-01-24 6 6 1e+06 2024-02-05 90 2 2e+05 0.000809336965969484 +87 4 2024-01-24 6 6 1e+06 2024-02-05 90 3 1e+05 0.000809336965969484 +87 5 7.68024834932647 2024-01-24 3.86712344199376 0 6 6 1e+06 2024-02-05 90 4 50000 0.000809336965969484 +88 1 2024-01-25 14 14 1e+06 2024-02-05 90 1 4e+05 0.00079714225263375 +88 2 2024-01-25 14 14 1e+06 2024-02-05 90 1 4e+05 0.00079714225263375 +88 3 2024-01-25 14 14 1e+06 2024-02-05 90 2 2e+05 0.00079714225263375 +88 4 2024-01-25 14 14 1e+06 2024-02-05 90 3 1e+05 0.00079714225263375 +88 5 2024-01-25 14 14 1e+06 2024-02-05 90 4 50000 0.00079714225263375 +89 1 2024-01-26 10 10 1e+06 2024-02-05 90 1 4e+05 0.000785402041327787 +89 2 2024-01-26 10 10 1e+06 2024-02-05 90 1 4e+05 0.000785402041327787 +89 3 2024-01-26 10 10 1e+06 2024-02-05 90 2 2e+05 0.000785402041327787 +89 4 2024-01-26 10 10 1e+06 2024-02-05 90 3 1e+05 0.000785402041327787 +89 5 2024-01-26 10 10 1e+06 2024-02-05 90 4 50000 0.000785402041327787 +90 1 8.9193378512175 2024-01-27 3.80549506716936 0 6 6 1e+06 2024-02-05 90 1 4e+05 0.00077384613079932 +90 2 2024-01-27 6 6 1e+06 2024-02-05 90 1 4e+05 0.00077384613079932 +90 3 8.13157806042153 2024-01-27 4.01075017260572 0 6 6 1e+06 2024-02-05 90 2 2e+05 0.00077384613079932 +90 4 2024-01-27 6 6 1e+06 2024-02-05 90 3 1e+05 0.00077384613079932 +90 5 2024-01-27 6 6 1e+06 2024-02-05 90 4 50000 0.00077384613079932 +91 1 2024-01-28 9 1e+06 2024-02-05 90 1 4e+05 0.000762484020063359 +91 2 2024-01-28 9 1e+06 2024-02-05 90 1 4e+05 0.000762484020063359 +91 3 2024-01-28 9 1e+06 2024-02-05 90 2 2e+05 0.000762484020063359 +91 4 2024-01-28 9 1e+06 2024-02-05 90 3 1e+05 0.000762484020063359 +91 5 6.79403430928 2024-01-28 3.86712344199376 0 9 1e+06 2024-02-05 90 4 50000 0.000762484020063359 +92 1 2024-01-29 6 1e+06 2024-02-05 90 1 4e+05 0.000799766434578809 +92 2 2024-01-29 6 1e+06 2024-02-05 90 1 4e+05 0.000799766434578809 +92 3 2024-01-29 6 1e+06 2024-02-05 90 2 2e+05 0.000799766434578809 +92 4 2024-01-29 6 1e+06 2024-02-05 90 3 1e+05 0.000799766434578809 +92 5 2024-01-29 6 1e+06 2024-02-05 90 4 50000 0.000799766434578809 +93 1 2024-01-30 12 1e+06 2024-02-05 90 1 4e+05 0.000795992564634426 +93 2 2024-01-30 12 1e+06 2024-02-05 90 1 4e+05 0.000795992564634426 +93 3 2024-01-30 12 1e+06 2024-02-05 90 2 2e+05 0.000795992564634426 +93 4 2024-01-30 12 1e+06 2024-02-05 90 3 1e+05 0.000795992564634426 +93 5 2024-01-30 12 1e+06 2024-02-05 90 4 50000 0.000795992564634426 +94 1 2024-01-31 8 1e+06 2024-02-05 90 1 4e+05 0.00080150456741431 +94 2 2024-01-31 8 1e+06 2024-02-05 90 1 4e+05 0.00080150456741431 +94 3 2024-01-31 8 1e+06 2024-02-05 90 2 2e+05 0.00080150456741431 +94 4 2024-01-31 8 1e+06 2024-02-05 90 3 1e+05 0.00080150456741431 +94 5 2024-01-31 8 1e+06 2024-02-05 90 4 50000 0.00080150456741431 +95 1 2024-02-01 9 1e+06 2024-02-05 90 1 4e+05 0.000807169781043701 +95 2 2024-02-01 9 1e+06 2024-02-05 90 1 4e+05 0.000807169781043701 +95 3 2024-02-01 9 1e+06 2024-02-05 90 2 2e+05 0.000807169781043701 +95 4 2024-02-01 9 1e+06 2024-02-05 90 3 1e+05 0.000807169781043701 +95 5 2024-02-01 9 1e+06 2024-02-05 90 4 50000 0.000807169781043701 +96 1 2024-02-02 7 1e+06 2024-02-05 90 1 4e+05 0.000812483294133183 +96 2 2024-02-02 7 1e+06 2024-02-05 90 1 4e+05 0.000812483294133183 +96 3 2024-02-02 7 1e+06 2024-02-05 90 2 2e+05 0.000812483294133183 +96 4 2024-02-02 7 1e+06 2024-02-05 90 3 1e+05 0.000812483294133183 +96 5 2024-02-02 7 1e+06 2024-02-05 90 4 50000 0.000812483294133183 +97 1 2024-02-03 7 1e+06 2024-02-05 90 1 4e+05 0.000817868754252128 +97 2 2024-02-03 7 1e+06 2024-02-05 90 1 4e+05 0.000817868754252128 +97 3 2024-02-03 7 1e+06 2024-02-05 90 2 2e+05 0.000817868754252128 +97 4 2024-02-03 7 1e+06 2024-02-05 90 3 1e+05 0.000817868754252128 +97 5 2024-02-03 7 1e+06 2024-02-05 90 4 50000 0.000817868754252128 +98 1 2024-02-04 9 1e+06 2024-02-05 90 1 4e+05 0.000823297300314913 +98 2 2024-02-04 9 1e+06 2024-02-05 90 1 4e+05 0.000823297300314913 +98 3 2024-02-04 9 1e+06 2024-02-05 90 2 2e+05 0.000823297300314913 +98 4 2024-02-04 9 1e+06 2024-02-05 90 3 1e+05 0.000823297300314913 +98 5 2024-02-04 9 1e+06 2024-02-05 90 4 50000 0.000823297300314913 +99 1 2024-02-05 8 1e+06 2024-02-05 90 1 4e+05 0.000849209201833048 +99 2 2024-02-05 8 1e+06 2024-02-05 90 1 4e+05 0.000849209201833048 +99 3 2024-02-05 8 1e+06 2024-02-05 90 2 2e+05 0.000849209201833048 +99 4 2024-02-05 8 1e+06 2024-02-05 90 3 1e+05 0.000849209201833048 +99 5 2024-02-05 8 1e+06 2024-02-05 90 4 50000 0.000849209201833048 +100 1 2024-02-06 7 1e+06 2024-02-05 90 1 4e+05 0.000858305930407811 +100 2 2024-02-06 7 1e+06 2024-02-05 90 1 4e+05 0.000858305930407811 +100 3 2024-02-06 7 1e+06 2024-02-05 90 2 2e+05 0.000858305930407811 +100 4 2024-02-06 7 1e+06 2024-02-05 90 3 1e+05 0.000858305930407811 +100 5 2024-02-06 7 1e+06 2024-02-05 90 4 50000 0.000858305930407811 +101 1 2024-02-07 4 1e+06 2024-02-05 90 1 4e+05 0.000871552235784833 +101 2 2024-02-07 4 1e+06 2024-02-05 90 1 4e+05 0.000871552235784833 +101 3 2024-02-07 4 1e+06 2024-02-05 90 2 2e+05 0.000871552235784833 +101 4 2024-02-07 4 1e+06 2024-02-05 90 3 1e+05 0.000871552235784833 +101 5 2024-02-07 4 1e+06 2024-02-05 90 4 50000 0.000871552235784833 +102 1 2024-02-08 8 1e+06 2024-02-05 90 1 4e+05 0.000885062046142915 +102 2 2024-02-08 8 1e+06 2024-02-05 90 1 4e+05 0.000885062046142915 +102 3 2024-02-08 8 1e+06 2024-02-05 90 2 2e+05 0.000885062046142915 +102 4 2024-02-08 8 1e+06 2024-02-05 90 3 1e+05 0.000885062046142915 +102 5 2024-02-08 8 1e+06 2024-02-05 90 4 50000 0.000885062046142915 +103 1 2024-02-09 7 1e+06 2024-02-05 90 1 4e+05 0.000898653655708608 +103 2 2024-02-09 7 1e+06 2024-02-05 90 1 4e+05 0.000898653655708608 +103 3 2024-02-09 7 1e+06 2024-02-05 90 2 2e+05 0.000898653655708608 +103 4 2024-02-09 7 1e+06 2024-02-05 90 3 1e+05 0.000898653655708608 +103 5 2024-02-09 7 1e+06 2024-02-05 90 4 50000 0.000898653655708608 +104 1 2024-02-10 7 1e+06 2024-02-05 90 1 4e+05 0.000912511967918031 +104 2 2024-02-10 7 1e+06 2024-02-05 90 1 4e+05 0.000912511967918031 +104 3 2024-02-10 7 1e+06 2024-02-05 90 2 2e+05 0.000912511967918031 +104 4 2024-02-10 7 1e+06 2024-02-05 90 3 1e+05 0.000912511967918031 +104 5 2024-02-10 7 1e+06 2024-02-05 90 4 50000 0.000912511967918031 +105 1 2024-02-11 12 1e+06 2024-02-05 90 1 4e+05 0.000926626990018999 +105 2 2024-02-11 12 1e+06 2024-02-05 90 1 4e+05 0.000926626990018999 +105 3 2024-02-11 12 1e+06 2024-02-05 90 2 2e+05 0.000926626990018999 +105 4 2024-02-11 12 1e+06 2024-02-05 90 3 1e+05 0.000926626990018999 +105 5 2024-02-11 12 1e+06 2024-02-05 90 4 50000 0.000926626990018999 +106 1 2024-02-12 7 1e+06 2024-02-05 90 1 4e+05 0.000913341840074717 +106 2 2024-02-12 7 1e+06 2024-02-05 90 1 4e+05 0.000913341840074717 +106 3 2024-02-12 7 1e+06 2024-02-05 90 2 2e+05 0.000913341840074717 +106 4 2024-02-12 7 1e+06 2024-02-05 90 3 1e+05 0.000913341840074717 +106 5 2024-02-12 7 1e+06 2024-02-05 90 4 50000 0.000913341840074717 +107 1 2024-02-13 7 1e+06 2024-02-05 90 1 4e+05 0.000922896656445218 +107 2 2024-02-13 7 1e+06 2024-02-05 90 1 4e+05 0.000922896656445218 +107 3 2024-02-13 7 1e+06 2024-02-05 90 2 2e+05 0.000922896656445218 +107 4 2024-02-13 7 1e+06 2024-02-05 90 3 1e+05 0.000922896656445218 +107 5 2024-02-13 7 1e+06 2024-02-05 90 4 50000 0.000922896656445218 +108 1 2024-02-14 7 1e+06 2024-02-05 90 1 4e+05 0.000927528395146316 +108 2 2024-02-14 7 1e+06 2024-02-05 90 1 4e+05 0.000927528395146316 +108 3 2024-02-14 7 1e+06 2024-02-05 90 2 2e+05 0.000927528395146316 +108 4 2024-02-14 7 1e+06 2024-02-05 90 3 1e+05 0.000927528395146316 +108 5 2024-02-14 7 1e+06 2024-02-05 90 4 50000 0.000927528395146316 +109 1 2024-02-15 8 1e+06 2024-02-05 90 1 4e+05 0.000932136792469634 +109 2 2024-02-15 8 1e+06 2024-02-05 90 1 4e+05 0.000932136792469634 +109 3 2024-02-15 8 1e+06 2024-02-05 90 2 2e+05 0.000932136792469634 +109 4 2024-02-15 8 1e+06 2024-02-05 90 3 1e+05 0.000932136792469634 +109 5 2024-02-15 8 1e+06 2024-02-05 90 4 50000 0.000932136792469634 +110 1 2024-02-16 4 1e+06 2024-02-05 90 1 4e+05 0.000937028946642744 +110 2 2024-02-16 4 1e+06 2024-02-05 90 1 4e+05 0.000937028946642744 +110 3 2024-02-16 4 1e+06 2024-02-05 90 2 2e+05 0.000937028946642744 +110 4 2024-02-16 4 1e+06 2024-02-05 90 3 1e+05 0.000937028946642744 +110 5 2024-02-16 4 1e+06 2024-02-05 90 4 50000 0.000937028946642744 +111 1 2024-02-17 7 1e+06 2024-02-05 90 1 4e+05 0.00094197897176175 +111 2 2024-02-17 7 1e+06 2024-02-05 90 1 4e+05 0.00094197897176175 +111 3 2024-02-17 7 1e+06 2024-02-05 90 2 2e+05 0.00094197897176175 +111 4 2024-02-17 7 1e+06 2024-02-05 90 3 1e+05 0.00094197897176175 +111 5 2024-02-17 7 1e+06 2024-02-05 90 4 50000 0.00094197897176175 +112 1 2024-02-18 5 1e+06 2024-02-05 90 1 4e+05 0.000947000679861772 +112 2 2024-02-18 5 1e+06 2024-02-05 90 1 4e+05 0.000947000679861772 +112 3 2024-02-18 5 1e+06 2024-02-05 90 2 2e+05 0.000947000679861772 +112 4 2024-02-18 5 1e+06 2024-02-05 90 3 1e+05 0.000947000679861772 +112 5 2024-02-18 5 1e+06 2024-02-05 90 4 50000 0.000947000679861772 +113 1 2024-02-19 10 1e+06 2024-02-05 90 1 4e+05 0.00102938326874019 +113 2 2024-02-19 10 1e+06 2024-02-05 90 1 4e+05 0.00102938326874019 +113 3 2024-02-19 10 1e+06 2024-02-05 90 2 2e+05 0.00102938326874019 +113 4 2024-02-19 10 1e+06 2024-02-05 90 3 1e+05 0.00102938326874019 +113 5 2024-02-19 10 1e+06 2024-02-05 90 4 50000 0.00102938326874019 +114 1 2024-02-20 3 1e+06 2024-02-05 90 1 4e+05 0.00104861781078916 +114 2 2024-02-20 3 1e+06 2024-02-05 90 1 4e+05 0.00104861781078916 +114 3 2024-02-20 3 1e+06 2024-02-05 90 2 2e+05 0.00104861781078916 +114 4 2024-02-20 3 1e+06 2024-02-05 90 3 1e+05 0.00104861781078916 +114 5 2024-02-20 3 1e+06 2024-02-05 90 4 50000 0.00104861781078916 +115 1 2024-02-21 10 1e+06 2024-02-05 90 1 4e+05 0.00108392374299587 +115 2 2024-02-21 10 1e+06 2024-02-05 90 1 4e+05 0.00108392374299587 +115 3 2024-02-21 10 1e+06 2024-02-05 90 2 2e+05 0.00108392374299587 +115 4 2024-02-21 10 1e+06 2024-02-05 90 3 1e+05 0.00108392374299587 +115 5 2024-02-21 10 1e+06 2024-02-05 90 4 50000 0.00108392374299587 +116 1 2024-02-22 6 1e+06 2024-02-05 90 1 4e+05 0.00112041404060539 +116 2 2024-02-22 6 1e+06 2024-02-05 90 1 4e+05 0.00112041404060539 +116 3 2024-02-22 6 1e+06 2024-02-05 90 2 2e+05 0.00112041404060539 +116 4 2024-02-22 6 1e+06 2024-02-05 90 3 1e+05 0.00112041404060539 +116 5 2024-02-22 6 1e+06 2024-02-05 90 4 50000 0.00112041404060539 +117 1 2024-02-23 9 1e+06 2024-02-05 90 1 4e+05 0.00115746811066289 +117 2 2024-02-23 9 1e+06 2024-02-05 90 1 4e+05 0.00115746811066289 +117 3 2024-02-23 9 1e+06 2024-02-05 90 2 2e+05 0.00115746811066289 +117 4 2024-02-23 9 1e+06 2024-02-05 90 3 1e+05 0.00115746811066289 +117 5 2024-02-23 9 1e+06 2024-02-05 90 4 50000 0.00115746811066289 +118 1 2024-02-24 7 1e+06 2024-02-05 90 1 4e+05 0.00119582408925606 +118 2 2024-02-24 7 1e+06 2024-02-05 90 1 4e+05 0.00119582408925606 +118 3 2024-02-24 7 1e+06 2024-02-05 90 2 2e+05 0.00119582408925606 +118 4 2024-02-24 7 1e+06 2024-02-05 90 3 1e+05 0.00119582408925606 +118 5 2024-02-24 7 1e+06 2024-02-05 90 4 50000 0.00119582408925606 +119 1 2024-02-25 7 1e+06 2024-02-05 90 1 4e+05 0.00123545884649584 +119 2 2024-02-25 7 1e+06 2024-02-05 90 1 4e+05 0.00123545884649584 +119 3 2024-02-25 7 1e+06 2024-02-05 90 2 2e+05 0.00123545884649584 +119 4 2024-02-25 7 1e+06 2024-02-05 90 3 1e+05 0.00123545884649584 +119 5 2024-02-25 7 1e+06 2024-02-05 90 4 50000 0.00123545884649584 +120 1 2024-02-26 6 1e+06 2024-02-05 90 1 4e+05 0.00131953405788595 +120 2 2024-02-26 6 1e+06 2024-02-05 90 1 4e+05 0.00131953405788595 +120 3 2024-02-26 6 1e+06 2024-02-05 90 2 2e+05 0.00131953405788595 +120 4 2024-02-26 6 1e+06 2024-02-05 90 3 1e+05 0.00131953405788595 +120 5 2024-02-26 6 1e+06 2024-02-05 90 4 50000 0.00131953405788595 +121 1 2024-02-27 7 1e+06 2024-02-05 90 1 4e+05 0.00137118304332127 +121 2 2024-02-27 7 1e+06 2024-02-05 90 1 4e+05 0.00137118304332127 +121 3 2024-02-27 7 1e+06 2024-02-05 90 2 2e+05 0.00137118304332127 +121 4 2024-02-27 7 1e+06 2024-02-05 90 3 1e+05 0.00137118304332127 +121 5 2024-02-27 7 1e+06 2024-02-05 90 4 50000 0.00137118304332127 +122 1 2024-02-28 14 1e+06 2024-02-05 90 1 4e+05 0.00143382583178741 +122 2 2024-02-28 14 1e+06 2024-02-05 90 1 4e+05 0.00143382583178741 +122 3 2024-02-28 14 1e+06 2024-02-05 90 2 2e+05 0.00143382583178741 +122 4 2024-02-28 14 1e+06 2024-02-05 90 3 1e+05 0.00143382583178741 +122 5 2024-02-28 14 1e+06 2024-02-05 90 4 50000 0.00143382583178741 +123 1 2024-02-29 13 1e+06 2024-02-05 90 1 4e+05 0.00149928295317067 +123 2 2024-02-29 13 1e+06 2024-02-05 90 1 4e+05 0.00149928295317067 +123 3 2024-02-29 13 1e+06 2024-02-05 90 2 2e+05 0.00149928295317067 +123 4 2024-02-29 13 1e+06 2024-02-05 90 3 1e+05 0.00149928295317067 +123 5 2024-02-29 13 1e+06 2024-02-05 90 4 50000 0.00149928295317067 +124 1 2024-03-01 9 1e+06 2024-02-05 90 1 4e+05 0.00156736770858743 +124 2 2024-03-01 9 1e+06 2024-02-05 90 1 4e+05 0.00156736770858743 +124 3 2024-03-01 9 1e+06 2024-02-05 90 2 2e+05 0.00156736770858743 +124 4 2024-03-01 9 1e+06 2024-02-05 90 3 1e+05 0.00156736770858743 +124 5 2024-03-01 9 1e+06 2024-02-05 90 4 50000 0.00156736770858743 +125 1 2024-03-02 15 1e+06 2024-02-05 90 1 4e+05 0.00163860318776472 +125 2 2024-03-02 15 1e+06 2024-02-05 90 1 4e+05 0.00163860318776472 +125 3 2024-03-02 15 1e+06 2024-02-05 90 2 2e+05 0.00163860318776472 +125 4 2024-03-02 15 1e+06 2024-02-05 90 3 1e+05 0.00163860318776472 +125 5 2024-03-02 15 1e+06 2024-02-05 90 4 50000 0.00163860318776472 +126 1 2024-03-03 10 1e+06 2024-02-05 90 1 4e+05 0.00171309115137499 +126 2 2024-03-03 10 1e+06 2024-02-05 90 1 4e+05 0.00171309115137499 +126 3 2024-03-03 10 1e+06 2024-02-05 90 2 2e+05 0.00171309115137499 +126 4 2024-03-03 10 1e+06 2024-02-05 90 3 1e+05 0.00171309115137499 +126 5 2024-03-03 10 1e+06 2024-02-05 90 4 50000 0.00171309115137499 +127 1 2024-03-04 10 1e+06 2024-02-05 90 1 4e+05 0.00163994741751005 +127 2 2024-03-04 10 1e+06 2024-02-05 90 1 4e+05 0.00163994741751005 +127 3 2024-03-04 10 1e+06 2024-02-05 90 2 2e+05 0.00163994741751005 +127 4 2024-03-04 10 1e+06 2024-02-05 90 3 1e+05 0.00163994741751005 +127 5 2024-03-04 10 1e+06 2024-02-05 90 4 50000 0.00163994741751005 diff --git a/model/src/pyrenew/deterministic/__init__.py b/model/src/pyrenew/deterministic/__init__.py index 87a504c0..66bbf5ad 100644 --- a/model/src/pyrenew/deterministic/__init__.py +++ b/model/src/pyrenew/deterministic/__init__.py @@ -2,5 +2,6 @@ from pyrenew.deterministic.deterministic import DeterministicVariable from pyrenew.deterministic.deterministicpmf import DeterministicPMF +from pyrenew.deterministic.process import DeterministicProcess -__all__ = ["DeterministicVariable", "DeterministicPMF"] +__all__ = ["DeterministicVariable", "DeterministicPMF", "DeterministicProcess"] diff --git a/model/src/pyrenew/deterministic/deterministic.py b/model/src/pyrenew/deterministic/deterministic.py index 5c07a7cc..a23cf93e 100644 --- a/model/src/pyrenew/deterministic/deterministic.py +++ b/model/src/pyrenew/deterministic/deterministic.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +from jax.typing import ArrayLike from pyrenew.metaclass import RandomVariable @@ -9,7 +10,7 @@ class DeterministicVariable(RandomVariable): def __init__( self, - vars: tuple, + vars: ArrayLike, label: str = "a_random_variable", ) -> None: """Default constructor @@ -33,14 +34,14 @@ def __init__( return None @staticmethod - def validate(vars: tuple) -> None: + def validate(vars: ArrayLike) -> None: """ Validates inputted to DeterministicPMF Parameters ---------- - vars : tuple - A tuple with arraylike objects. + vars : ArrayLike + An ArrayLike object. Returns ------- @@ -49,10 +50,10 @@ def validate(vars: tuple) -> None: Raises ------ Exception - If the inputted vars object is not a tuple. + If the inputted vars object is not a ArrayLike. """ - if not isinstance(vars, tuple): - raise Exception("vars is not a tuple") + if not isinstance(vars, ArrayLike): + raise Exception("vars is not a ArrayLike") return None @@ -74,4 +75,4 @@ def sample( Containing the stored values during construction. """ - return self.vars + return (self.vars,) diff --git a/model/src/pyrenew/deterministic/deterministicpmf.py b/model/src/pyrenew/deterministic/deterministicpmf.py index 0ed71323..1a0c9450 100644 --- a/model/src/pyrenew/deterministic/deterministicpmf.py +++ b/model/src/pyrenew/deterministic/deterministicpmf.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +from jax.typing import ArrayLike from pyrenew.deterministic.deterministic import DeterministicVariable from pyrenew.distutil import validate_discrete_dist_vector from pyrenew.metaclass import RandomVariable @@ -11,7 +12,7 @@ class DeterministicPMF(RandomVariable): def __init__( self, - vars: tuple, + vars: ArrayLike, label: str = "a_random_variable", tol: float = 1e-20, ) -> None: @@ -33,32 +34,28 @@ def __init__( Passed to pyrenew.distutil.validate_discrete_dist_vector. Defaults to 1e-20. - Returns ------- None """ + vars = validate_discrete_dist_vector( + discrete_dist=vars, + tol=tol, + ) - vars2 = list(vars) - for i in range(0, len(vars2)): - vars2[i] = validate_discrete_dist_vector( - discrete_dist=vars2[i], - tol=tol, - ) - - self.basevar = DeterministicVariable(tuple(vars2), label) + self.basevar = DeterministicVariable(vars, label) return None @staticmethod - def validate(vars: tuple) -> None: + def validate(vars: ArrayLike) -> None: """ Validates inputted to DeterministicPMF Parameters ---------- - vars : tuple - A tuple with arraylike objects. + vars : ArrayLike + An ArrayLike object. Returns ------- diff --git a/model/src/pyrenew/deterministic/process.py b/model/src/pyrenew/deterministic/process.py new file mode 100644 index 00000000..c46206a6 --- /dev/null +++ b/model/src/pyrenew/deterministic/process.py @@ -0,0 +1,37 @@ +import jax.numpy as jnp +from pyrenew.deterministic.deterministic import DeterministicVariable + + +class DeterministicProcess(DeterministicVariable): + """A deterministic process (degenerate) random variable. + Useful to pass fixed quantities over time.""" + + __init__ = DeterministicVariable.__init__ + + def sample( + self, + n_timepoints: int, + **kwargs, + ) -> tuple: + """Retrieve the value of the deterministic Rv + Parameters + ---------- + n_timepoints : int + Number of timepoints to sample. + **kwargs : dict, optional + Ignored. + + Returns + ------- + tuple + Containing the stored values during construction. + """ + + res, *_ = super().sample(**kwargs) + + dif = n_timepoints - res.shape[0] + + if dif > 0: + return (jnp.hstack([res, jnp.repeat(res[-1], dif)]),) + + return (res[:n_timepoints],) diff --git a/model/src/pyrenew/latent/hospitaladmissions.py b/model/src/pyrenew/latent/hospitaladmissions.py index 7ff6af4b..e04a7b25 100644 --- a/model/src/pyrenew/latent/hospitaladmissions.py +++ b/model/src/pyrenew/latent/hospitaladmissions.py @@ -158,7 +158,7 @@ class HospitalAdmissions(RandomVariable): .. math:: - H(t) = \omega(t) p_\mathrm{hosp}(t) \sum_{\\tau = 0}^{T_d} d(\tau) I(t-\tau) + H(t) = \omega(t) p_\mathrm{hosp}(t) \sum_{\tau = 0}^{T_d} d(\tau) I(t-\tau) Where :math:`T_d` is the maximum delay from infection to hospitalization that we consider. @@ -196,9 +196,9 @@ def __init__( """ if weekday_effect_dist is None: - weekday_effect_dist = DeterministicVariable((1,)) + weekday_effect_dist = DeterministicVariable(1) if hosp_report_prob_dist is None: - hosp_report_prob_dist = DeterministicVariable((1,)) + hosp_report_prob_dist = DeterministicVariable(1) HospitalAdmissions.validate( infect_hosp_rate_dist, diff --git a/model/src/pyrenew/latent/infection_functions.py b/model/src/pyrenew/latent/infection_functions.py index ec5dfad8..8d0b62fc 100755 --- a/model/src/pyrenew/latent/infection_functions.py +++ b/model/src/pyrenew/latent/infection_functions.py @@ -44,7 +44,7 @@ def sample_infections_rt( """ incidence_func = new_convolve_scanner(reversed_generation_interval_pmf) - latest, all_infections = jax.lax.scan(incidence_func, I0, Rt) + latest, all_infections = jax.lax.scan(f=incidence_func, init=I0, xs=Rt) return all_infections diff --git a/model/src/pyrenew/latent/infections.py b/model/src/pyrenew/latent/infections.py index a9d0543e..a0effa46 100644 --- a/model/src/pyrenew/latent/infections.py +++ b/model/src/pyrenew/latent/infections.py @@ -24,11 +24,27 @@ class Infections(RandomVariable): - """Latent infections + r"""Latent infections + + This class samples infections given Rt, initial infections, and generation + interval. + + Notes + ----- + The mathematical model is given by: + + .. math:: + + I(t) = R(t) \times \sum_{\tau < t} I(\tau) g(t-\tau) + + where :math:`I(t)` is the number of infections at time :math:`t`, + :math:`R(t)` is the reproduction number at time :math:`t`, and + :math:`g(t-\tau)` is the generation interval. Methods ------- sample(Rt, I0, gen_int, **kwargs) + """ def __init__( @@ -67,9 +83,8 @@ def sample( gen_int: ArrayLike, **kwargs, ) -> tuple: - """Samples infections given Rt - - A random variable representing the pmf of the generation interval. + """Samples infections given Rt, initial infections, and generation + interval. Parameters ---------- diff --git a/model/src/pyrenew/mcmcutils.py b/model/src/pyrenew/mcmcutils.py index 733d7c75..343a9489 100644 --- a/model/src/pyrenew/mcmcutils.py +++ b/model/src/pyrenew/mcmcutils.py @@ -4,8 +4,10 @@ Utilities to deal with MCMC outputs """ +import matplotlib.pyplot as plt import numpy as np import polars as pl +from jax.typing import ArrayLike def spread_draws( @@ -74,3 +76,79 @@ def spread_draws( pass return df + + +def plot_posterior( + var: str, + draws: pl.DataFrame, + obs_signal: ArrayLike = None, + ylab: str = None, + xlab: str = "Time", + samples: int = 50, + figsize: list = [4, 5], + draws_col: str = "darkblue", + obs_col: str = "black", +) -> plt.Figure: + """ + Plot the posterior distribution of a variable + + Parameters + ---------- + var : str + Name of the variable to plot + model : Model + Model object + obs_signal : ArrayLike, optional + Observed signal to plot as reference + ylab : str, optional + Label for the y-axis + xlab : str, optional + Label for the x-axis + samples : int, optional + Number of samples to plot + figsize : list, optional + Size of the figure + draws_col : str, optional + Color of the draws + obs_col : str, optional + + Returns + ------- + plt.Figure + """ + + if ylab is None: + ylab = var + + fig, ax = plt.subplots(figsize=figsize) + + # Reference signal (if any) + if obs_signal is not None: + ax.plot(obs_signal, color=obs_col) + + samp_ids = np.random.randint(size=samples, low=0, high=999) + + for samp_id in samp_ids: + sub_samps = draws.filter(pl.col("draw") == samp_id).sort( + pl.col("time") + ) + ax.plot( + sub_samps.select("time").to_numpy(), + sub_samps.select(var).to_numpy(), + color=draws_col, + alpha=0.1, + ) + + # Some labels + ax.set_xlabel(xlab) + ax.set_ylabel(ylab) + + # Adding a legend + ax.plot([], [], color=draws_col, alpha=0.9, label="Posterior samples") + + if obs_signal is not None: + ax.plot([], [], color=obs_col, label="Observed signal") + + ax.legend() + + return fig diff --git a/model/src/pyrenew/metaclass.py b/model/src/pyrenew/metaclass.py index 5b51584d..ef2626c5 100644 --- a/model/src/pyrenew/metaclass.py +++ b/model/src/pyrenew/metaclass.py @@ -7,9 +7,10 @@ from abc import ABCMeta, abstractmethod import jax +import matplotlib.pyplot as plt import polars as pl from numpyro.infer import MCMC, NUTS -from pyrenew.mcmcutils import spread_draws +from pyrenew.mcmcutils import plot_posterior, spread_draws def _assert_sample_and_rtype( @@ -285,3 +286,27 @@ def spread_draws(self, variables_names: list) -> pl.DataFrame: """ return spread_draws(self.mcmc.get_samples(), variables_names) + + def plot_posterior( + self, + var: list, + obs_signal: jax.typing.ArrayLike = None, + xlab: str = None, + ylab: str = "Signal", + samples: int = 50, + figsize: list = [4, 5], + draws_col: str = "darkblue", + obs_col: str = "black", + ) -> plt.Figure: + """A wrapper of pyrenew.mcmcutils.plot_posterior""" + return plot_posterior( + var=var, + draws=self.spread_draws([(var, "time")]), + xlab=xlab, + ylab=ylab, + samples=samples, + obs_signal=obs_signal, + figsize=figsize, + draws_col=draws_col, + obs_col=obs_col, + ) diff --git a/model/src/pyrenew/model/hospitalizations.py b/model/src/pyrenew/model/hospitalizations.py index 0bddca40..e59535b0 100644 --- a/model/src/pyrenew/model/hospitalizations.py +++ b/model/src/pyrenew/model/hospitalizations.py @@ -2,8 +2,8 @@ from collections import namedtuple +import jax.numpy as jnp from numpy.typing import ArrayLike -from pyrenew.deterministic import DeterministicVariable from pyrenew.metaclass import Model, RandomVariable, _assert_sample_and_rtype from pyrenew.model.rtinfectionsrenewal import RtInfectionsRenewalModel @@ -11,10 +11,10 @@ "HospModelSample", [ "Rt", - "infections", + "latent_infections", "IHR", - "latent", - "sampled", + "latent_admissions", + "sampled_admissions", ], defaults=[None, None, None, None, None], ) @@ -25,13 +25,13 @@ ---------- Rt : float or None The reproduction number over time. Defaults to None. -infections : ArrayLike or None +latent_infections : ArrayLike or None The estimated number of new infections over time. Defaults to None. IHR : float or None The infected hospitalization rate. Defaults to None. -latent : ArrayLike or None +latent_admissions : ArrayLike or None The estimated latent hospitalizations. Defaults to None. -sampled : ArrayLike or None +sampled_admissions : ArrayLike or None The sampled or observed hospital admissions. Defaults to None. Notes @@ -55,7 +55,7 @@ def __init__( gen_int: RandomVariable, I0: RandomVariable, Rt_process: RandomVariable, - observed_hospitalizations: RandomVariable, + observation_process: RandomVariable, ) -> None: """ Default constructor @@ -72,7 +72,7 @@ def __init__( Initial infections (passed to RtInfectionsRenewalModel) Rt_process : RandomVariable Rt process (passed to RtInfectionsRenewalModel). - observed_hospitalizations : RandomVariable, optional + observation_process : RandomVariable, optional Observation process for the hospitalizations. Returns @@ -87,19 +87,19 @@ def __init__( gen_int=gen_int, I0=I0, latent_infections=latent_infections, - observed_infections=DeterministicVariable((0,)), + observation_process=None, Rt_process=Rt_process, ) HospitalizationsModel.validate( - latent_hospitalizations, observed_hospitalizations + latent_hospitalizations, observation_process ) self.latent_hospitalizations = latent_hospitalizations - self.observed_hospitalizations = observed_hospitalizations + self.observation_process = observation_process @staticmethod - def validate(latent_hospitalizations, observed_hospitalizations) -> None: + def validate(latent_hospitalizations, observation_process) -> None: """ Verifies types and status (RV) of latent and observed hospitalizations @@ -107,7 +107,7 @@ def validate(latent_hospitalizations, observed_hospitalizations) -> None: ---------- latent_hospitalizations : ArrayLike The latent process for the hospitalizations. - observed_hospitalizations : ArrayLike + observation_process : ArrayLike The observed hospitalizations. Returns @@ -119,7 +119,7 @@ def validate(latent_hospitalizations, observed_hospitalizations) -> None: _assert_sample_and_rtype : Perform type-checking and verify RV """ _assert_sample_and_rtype(latent_hospitalizations, skip_if_none=False) - _assert_sample_and_rtype(observed_hospitalizations, skip_if_none=False) + _assert_sample_and_rtype(observation_process, skip_if_none=False) return None def sample_hospitalizations_latent( @@ -160,6 +160,7 @@ def sample_hospitalizations_obs( self, predicted: ArrayLike, observed_hospitalizations: ArrayLike, + name: str | None = None, **kwargs, ) -> tuple: """Sample number of hospitalizations @@ -170,6 +171,8 @@ def sample_hospitalizations_obs( The predicted hospitalizations. observed_hospitalizations : ArrayLike The observed hospitalization data (to fit). + name : str, optional + Name of the random variable. Defaults to None. **kwargs : dict, optional Additional keyword arguments passed through to internal sample_hospitalizations_obs calls, should there be any. @@ -178,23 +181,23 @@ def sample_hospitalizations_obs( ------- tuple - See Also - -------- - observed_hospitalizations.sample : For sampling observed hospitalizations - Notes ----- TODO: Include example(s) here. """ - return self.observed_hospitalizations.sample( - predicted=predicted, obs=observed_hospitalizations, **kwargs + return self.observation_process.sample( + predicted=predicted, + obs=observed_hospitalizations, + name=name, + **kwargs, ) def sample( self, n_timepoints: int, observed_hospitalizations: ArrayLike | None = None, + padding: int = 0, **kwargs, ) -> HospModelSample: """ @@ -207,6 +210,9 @@ def sample( observed_hospitalizations : ArrayLike, optional The observed hospitalization data (passed to the basic renewal model). Defaults to None (simulation, rather than fit). + padding : int, optional + Number of padding timepoints to add to the beginning of the + simulation. Defaults to 0. **kwargs : dict, optional Additional keyword arguments passed through to internal sample() calls, should there be any. @@ -227,9 +233,10 @@ def sample( """ # Getting the initial quantities from the basic model - Rt, infections, *_ = self.basic_renewal.sample( + basic_model = self.basic_renewal.sample( n_timepoints=n_timepoints, observed_infections=None, + padding=padding, **kwargs, ) @@ -239,21 +246,38 @@ def sample( latent, *_, ) = self.sample_hospitalizations_latent( - infections=infections, + infections=basic_model.latent_infections, **kwargs, ) # Sampling the hospitalizations - sampled, *_ = self.sample_hospitalizations_obs( - predicted=latent, - observed_hospitalizations=observed_hospitalizations, - **kwargs, - ) + if self.observation_process is not None: + if (observed_hospitalizations is not None) and (padding > 0): + sampled_na = jnp.repeat(jnp.nan, padding) + + sampled_observed, *_ = self.sample_hospitalizations_obs( + predicted=latent[padding:], + observed_hospitalizations=observed_hospitalizations[ + padding: + ], + **kwargs, + ) + + sampled = jnp.hstack([sampled_na, sampled_observed]) + + else: + sampled, *_ = self.sample_hospitalizations_obs( + predicted=latent, + observed_hospitalizations=observed_hospitalizations, + **kwargs, + ) + else: + sampled = None return HospModelSample( - Rt=Rt, - infections=infections, + Rt=basic_model.Rt, + latent_infections=basic_model.latent_infections, IHR=IHR, - latent=latent, - sampled=sampled, + latent_admissions=latent, + sampled_admissions=sampled, ) diff --git a/model/src/pyrenew/model/rtinfectionsrenewal.py b/model/src/pyrenew/model/rtinfectionsrenewal.py index 253948ff..154847b3 100644 --- a/model/src/pyrenew/model/rtinfectionsrenewal.py +++ b/model/src/pyrenew/model/rtinfectionsrenewal.py @@ -3,13 +3,14 @@ from collections import namedtuple from typing import Optional +import jax.numpy as jnp from numpy.typing import ArrayLike from pyrenew.metaclass import Model, RandomVariable, _assert_sample_and_rtype # Output class of the RtInfectionsRenewalModel RtInfectionsRenewalSample = namedtuple( "InfectModelSample", - ["Rt", "latent", "observed"], + ["Rt", "latent_infections", "sampled_infections"], defaults=[None, None, None], ) RtInfectionsRenewalSample.__doc__ = """ @@ -19,10 +20,10 @@ ---------- Rt : float or None The reproduction number over time. Defaults to None. -latent : ArrayLike or None +latent_infections : ArrayLike or None The estimated latent infections. Defaults to None. -observed : ArrayLike or None - The observed infections. Defaults to None. +sampled_infections : ArrayLike or None + The sampled infections. Defaults to None. Notes ----- @@ -42,7 +43,7 @@ def __init__( gen_int: RandomVariable, I0: RandomVariable, Rt_process: RandomVariable, - observed_infections: RandomVariable, + observation_process: RandomVariable, ) -> None: """Default constructor @@ -58,7 +59,7 @@ def __init__( Rt_process : RandomVariable The sample function of the process should return a tuple where the first element is the drawn Rt. - observed_infections : RandomVariable + observation_process : RandomVariable Infections observation process (e.g., pyrenew.observations.Poisson.). @@ -71,14 +72,14 @@ def __init__( gen_int=gen_int, i0=I0, latent_infections=latent_infections, - observed_infections=observed_infections, + observation_process=observation_process, Rt_process=Rt_process, ) self.gen_int = gen_int self.i0 = I0 self.latent_infections = latent_infections - self.observed_infections = observed_infections + self.observation_process = observation_process self.Rt_process = Rt_process @staticmethod @@ -86,7 +87,7 @@ def validate( gen_int, i0, latent_infections, - observed_infections, + observation_process, Rt_process, ) -> None: """ @@ -111,7 +112,7 @@ def validate( _assert_sample_and_rtype(gen_int, skip_if_none=False) _assert_sample_and_rtype(i0, skip_if_none=False) _assert_sample_and_rtype(latent_infections, skip_if_none=False) - _assert_sample_and_rtype(observed_infections, skip_if_none=True) + _assert_sample_and_rtype(observation_process, skip_if_none=True) _assert_sample_and_rtype(Rt_process, skip_if_none=False) return None @@ -211,6 +212,7 @@ def sample_infections_obs( self, predicted: ArrayLike, observed_infections: Optional[ArrayLike] = None, + name: str | None = None, **kwargs, ) -> tuple: """Sample number of hospitalizations @@ -221,6 +223,8 @@ def sample_infections_obs( The predicted infecteds. observed_hospitalizations : ArrayLike, optional The observed values of hospital admissions, if any, for inference. Defaults to None. + name : str, optional + Name of the random variable passed to the RandomVariable. **kwargs : dict, optional Additional keyword arguments passed through to internal sample() calls, should there be any. @@ -229,17 +233,14 @@ def sample_infections_obs( ------- tuple - See Also - -------- - observed_infections.sample : For sampling observed infections - Notes ----- TODO: Include example(s) here. """ - return self.observed_infections.sample( + return self.observation_process.sample( predicted=predicted, obs=observed_infections, + name=name, **kwargs, ) @@ -247,6 +248,7 @@ def sample( self, n_timepoints: int, observed_infections: Optional[ArrayLike] = None, + padding: int = 0, **kwargs, ) -> RtInfectionsRenewalSample: """Sample from the Basic Renewal Model @@ -257,6 +259,9 @@ def sample( Number of timepoints to sample. observed_infections : ArrayLike, optional Observed infections. + padding : int, optional + Number of padding timepoints to add to the beginning of the + simulation. Defaults to 0. **kwargs : dict, optional Additional keyword arguments passed through to internal sample() calls, if any @@ -292,12 +297,29 @@ def sample( ) # Using the predicted infections to sample from the observation process - observed, *_ = self.sample_infections_obs( - predicted=latent, observed_infections=observed_infections, **kwargs - ) + if self.observation_process is not None: + if (observed_infections is not None) and (padding > 0): + sampled_pad = jnp.repeat(jnp.nan, padding) + + sampled_obs, *_ = self.sample_infections_obs( + predicted=latent[padding:], + observed_infections=observed_infections[padding:], + **kwargs, + ) + + sampled = jnp.hstack([sampled_pad, sampled_obs]) + + else: + sampled, *_ = self.sample_infections_obs( + predicted=latent, + observed_infections=observed_infections, + **kwargs, + ) + else: + sampled = None return RtInfectionsRenewalSample( Rt=Rt, - latent=latent, - observed=observed, + latent_infections=latent, + sampled_infections=sampled, ) diff --git a/model/src/pyrenew/observation/negativebinomial.py b/model/src/pyrenew/observation/negativebinomial.py index e4db6ee6..f03538c8 100644 --- a/model/src/pyrenew/observation/negativebinomial.py +++ b/model/src/pyrenew/observation/negativebinomial.py @@ -25,6 +25,7 @@ def __init__( concentration_prior: dist.Distribution | ArrayLike, concentration_suffix: str | None = "_concentration", parameter_name="negbinom_rv", + eps: float = 1e-10, ) -> None: """Default constructor @@ -41,6 +42,9 @@ def __init__( Suffix for the numpy variable. parameter_name : str, optional Name for the numpy variable. + eps : float, optional + Small value to add to the predicted mean to prevent numerical + instability. Returns ------- @@ -59,10 +63,12 @@ def __init__( self.parameter_name = parameter_name self.concentration_suffix = concentration_suffix + self.eps = eps def sample( self, predicted: ArrayLike, + name: str | None = None, obs: ArrayLike | None = None, **kwargs, ) -> tuple: @@ -74,6 +80,9 @@ def sample( Mean parameter of the negative binomial distribution. obs : ArrayLike, optional Observed data, by default None. + name : str, optional + Name of the random variable if other than that defined during + construction, by default None (self.parameter_name). **kwargs : dict, optional Additional keyword arguments passed through to internal sample calls, should there be any. @@ -81,12 +90,17 @@ def sample( ------- tuple """ + concentration = self.sample_prior() + + if name is None: + name = self.parameter_name + return ( numpyro.sample( - self.parameter_name, - dist.NegativeBinomial2( - mean=predicted, - concentration=self.sample_prior(), + name=name, + fn=dist.NegativeBinomial2( + mean=predicted + self.eps, + concentration=concentration, ), obs=obs, ), diff --git a/model/src/pyrenew/observation/poisson.py b/model/src/pyrenew/observation/poisson.py index 0600199c..5feaf0e8 100644 --- a/model/src/pyrenew/observation/poisson.py +++ b/model/src/pyrenew/observation/poisson.py @@ -45,6 +45,7 @@ def sample( self, predicted: ArrayLike, obs: ArrayLike | None = None, + name: str | None = None, **kwargs, ) -> tuple: """Sample from the Poisson process @@ -55,6 +56,8 @@ def sample( Rate parameter of the Poisson distribution. obs : ArrayLike, optional Observed data. Defaults to None. + name : str, optional + Name of the random variable. Defaults to None. **kwargs : dict, optional Additional keyword arguments passed through to internal sample calls, should there be any. @@ -62,10 +65,14 @@ def sample( ------- tuple """ + + if name is None: + name = self.parameter_name + return ( numpyro.sample( - self.parameter_name, - dist.Poisson(rate=predicted + self.eps), + name=name, + fn=dist.Poisson(rate=predicted + self.eps), obs=obs, ), ) diff --git a/model/src/pyrenew/process/rtrandomwalk.py b/model/src/pyrenew/process/rtrandomwalk.py index c6cb32c9..e8f31ea1 100644 --- a/model/src/pyrenew/process/rtrandomwalk.py +++ b/model/src/pyrenew/process/rtrandomwalk.py @@ -8,7 +8,18 @@ class RtRandomWalkProcess(RandomVariable): - """Rt Randomwalk Process""" + r"""Rt Randomwalk Process + + Notes + ----- + + The process is defined as follows: + + .. math:: + + Rt(0) &\sim \text{Rt0_dist} \\ + Rt(t) &\sim \text{Rt_transform}(\text{Rt_transformed_rw}(t)) + """ def __init__( self, @@ -102,7 +113,9 @@ def sample( Rt0_trans = self.Rt_transform(Rt0) Rt_trans_proc = SimpleRandomWalkProcess(self.Rt_rw_dist) Rt_trans_ts, *_ = Rt_trans_proc.sample( - duration=n_timepoints, name="Rt_transformed_rw", init=Rt0_trans + duration=n_timepoints, + name="Rt_transformed_rw", + init=Rt0_trans, ) Rt = npro.deterministic("Rt", self.Rt_transform.inverse(Rt_trans_ts)) diff --git a/model/src/test/baseline/test_model_basicrenewal_plot.png b/model/src/test/baseline/test_model_basicrenewal_plot.png new file mode 100644 index 00000000..80a0d656 Binary files /dev/null and b/model/src/test/baseline/test_model_basicrenewal_plot.png differ diff --git a/model/src/test/test_datasets.py b/model/src/test/test_datasets.py new file mode 100644 index 00000000..a8a3875a --- /dev/null +++ b/model/src/test/test_datasets.py @@ -0,0 +1,36 @@ +import numpy.testing as testing +from pyrenew.datasets import ( + load_generation_interval, + load_infection_admission_interval, + load_wastewater, +) + + +def test_loading_wastewater(): + """Test that the wastewater dataset can be properly loaded""" + df = load_wastewater() + assert len(df) > 0 + assert df.shape == (635, 14) + + testing.assert_approx_equal(df["daily_hosp_admits"].mean(), 12.8888, 3) + testing.assert_approx_equal(df["lod_sewage"].mean(), 3.841025, 3) + + +def test_gen_int(): + """Test that the generation interval dataset can be properly loaded""" + df = load_generation_interval() + assert len(df) > 0 + assert df.shape == (15, 2) + + testing.assert_approx_equal(df["probability_mass"].mean(), 0.0666666, 3) + + +def test_infection_admission_interval(): + """Test that the infection to admission interval dataset can be properly loaded""" + df = load_infection_admission_interval() + assert len(df) > 0 + assert df.shape == (55, 2) + + testing.assert_approx_equal( + df["probability_mass"].mean(), 0.01818181818, 3 + ) diff --git a/model/src/test/test_deterministic.py b/model/src/test/test_deterministic.py new file mode 100644 index 00000000..c52e0f90 --- /dev/null +++ b/model/src/test/test_deterministic.py @@ -0,0 +1,51 @@ +import jax.numpy as jnp +import numpy.testing as testing +from pyrenew.deterministic import ( + DeterministicPMF, + DeterministicProcess, + DeterministicVariable, +) + + +def test_deterministic(): + """Test the DeterministicVariable, DeterministicPMF, and + DeterministicProcess classes in the deterministic module. + """ + + var1 = DeterministicVariable( + jnp.array( + [ + 1, + ] + ) + ) + var2 = DeterministicPMF(jnp.array([0.25, 0.25, 0.2, 0.3])) + var3 = DeterministicProcess(jnp.array([1, 2, 3, 4])) + + testing.assert_array_equal( + var1.sample()[0], + jnp.array( + [ + 1, + ] + ), + ) + testing.assert_array_equal( + var2.sample()[0], + jnp.array([0.25, 0.25, 0.2, 0.3]), + ) + testing.assert_array_equal( + var3.sample(n_timepoints=5)[0], + jnp.array([1, 2, 3, 4, 4]), + ) + + testing.assert_array_equal( + var3.sample(n_timepoints=3)[0], + jnp.array( + [ + 1, + 2, + 3, + ] + ), + ) diff --git a/model/src/test/test_latent_hospitalizations.py b/model/src/test/test_latent_hospitalizations.py index 9a49d9ee..056ce1e9 100644 --- a/model/src/test/test_latent_hospitalizations.py +++ b/model/src/test/test_latent_hospitalizations.py @@ -32,31 +32,30 @@ def test_hospitalizations_sample(): # Testing the hospitalizations inf_hosp = DeterministicPMF( - ( - jnp.array( - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0.25, - 0.5, - 0.1, - 0.1, - 0.05, - ] - ), + jnp.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0.25, + 0.5, + 0.1, + 0.1, + 0.05, + ] ), ) + hosp1 = HospitalAdmissions( infection_to_admission_interval=inf_hosp, infect_hosp_rate_dist=InfectHospRate( diff --git a/model/src/test/test_model_basic_renewal.py b/model/src/test/test_model_basic_renewal.py index 83aa27a5..2afebe2e 100644 --- a/model/src/test/test_model_basic_renewal.py +++ b/model/src/test/test_model_basic_renewal.py @@ -2,10 +2,12 @@ import jax import jax.numpy as jnp +import matplotlib.pyplot as plt import numpy as np import numpyro as npro import numpyro.distributions as dist import polars as pl +import pytest from pyrenew.deterministic import DeterministicPMF, DeterministicVariable from pyrenew.latent import Infections, Infections0 from pyrenew.model import RtInfectionsRenewalModel @@ -19,9 +21,7 @@ def test_model_basicrenewal_no_obs_model(): from the perspective of the infections. It returns expected, not sampled. """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -34,7 +34,7 @@ def test_model_basicrenewal_no_obs_model(): I0=I0, latent_infections=latent_infections, Rt_process=rt, - observed_infections=DeterministicVariable((1,)), + observation_process=DeterministicVariable(1), ) # Sampling and fitting model 0 (with no obs for infections) @@ -46,7 +46,7 @@ def test_model_basicrenewal_no_obs_model(): num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(272), - observed_infections=model0_samp.observed, + observed_infections=model0_samp.sampled_infections, n_timepoints=30, ) @@ -68,10 +68,91 @@ def test_model_basicrenewal_with_obs_model(): from the perspective of the infections. It returns sampled, not expected. """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) + + I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) + + latent_infections = Infections() + + observed_infections = PoissonObservation() + + rt = RtRandomWalkProcess() + + model1 = RtInfectionsRenewalModel( + I0=I0, + gen_int=gen_int, + latent_infections=latent_infections, + observation_process=observed_infections, + Rt_process=rt, + ) + + # Sampling and fitting model 1 (with obs infections) + np.random.seed(2203) + with npro.handlers.seed(rng_seed=np.random.randint(1, 600)): + model1_samp = model1.sample(n_timepoints=30) + + model1.run( + num_warmup=500, + num_samples=500, + rng_key=jax.random.PRNGKey(22), + observed_infections=model1_samp.sampled_infections, + n_timepoints=30, + ) + + inf = model1.spread_draws(["latent_infections"]) + inf_mean = ( + inf.group_by("draw") + .agg(pl.col("latent_infections").mean()) + .sort(pl.col("draw")) + ) + + # For now the assertion is only about the expected number of rows + # It should be about the MCMC inference. + assert inf_mean.to_numpy().shape[0] == 500 + + +@pytest.mark.mpl_image_compare +def test_model_basicrenewal_plot() -> plt.Figure: + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) + + I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) + + latent_infections = Infections() + + observed_infections = PoissonObservation() + + rt = RtRandomWalkProcess() + + model1 = RtInfectionsRenewalModel( + I0=I0, + gen_int=gen_int, + latent_infections=latent_infections, + observation_process=observed_infections, + Rt_process=rt, ) + # Sampling and fitting model 1 (with obs infections) + np.random.seed(2203) + with npro.handlers.seed(rng_seed=np.random.randint(1, 600)): + model1_samp = model1.sample(n_timepoints=30) + + model1.run( + num_warmup=500, + num_samples=500, + rng_key=jax.random.PRNGKey(22), + observed_infections=model1_samp.sampled_infections, + n_timepoints=30, + ) + + return model1.plot_posterior( + var="latent_infections", + obs_signal=model1_samp.sampled_infections, + ) + + +def test_model_basicrenewal_padding() -> None: + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) + I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) latent_infections = Infections() @@ -84,7 +165,7 @@ def test_model_basicrenewal_with_obs_model(): I0=I0, gen_int=gen_int, latent_infections=latent_infections, - observed_infections=observed_infections, + observation_process=observed_infections, Rt_process=rt, ) @@ -93,15 +174,21 @@ def test_model_basicrenewal_with_obs_model(): with npro.handlers.seed(rng_seed=np.random.randint(1, 600)): model1_samp = model1.sample(n_timepoints=30) + new_obs = jnp.hstack( + [jnp.repeat(jnp.nan, 5), model1_samp.sampled_infections[5:]], + ) + model1.run( num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(22), - observed_infections=model1_samp.observed, + observed_infections=new_obs, n_timepoints=30, + padding=5, ) inf = model1.spread_draws(["latent_infections"]) + inf_mean = ( inf.group_by("draw") .agg(pl.col("latent_infections").mean()) diff --git a/model/src/test/test_model_hospitalizations.py b/model/src/test/test_model_hospitalizations.py index 11c7844a..34dba6ab 100644 --- a/model/src/test/test_model_hospitalizations.py +++ b/model/src/test/test_model_hospitalizations.py @@ -41,40 +41,37 @@ def test_model_hosp_no_obs_model(): Hospitalization model runs """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) latent_infections = Infections() Rt_process = RtRandomWalkProcess() inf_hosp = DeterministicPMF( - ( - jnp.array( - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0.25, - 0.5, - 0.1, - 0.1, - 0.05, - ], - ), + jnp.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0.25, + 0.5, + 0.1, + 0.1, + 0.05, + ] ), ) + latent_hospitalizations = HospitalAdmissions( infection_to_admission_interval=inf_hosp, hospitalizations_predicted_varname="observed_hospitalizations", @@ -89,7 +86,7 @@ def test_model_hosp_no_obs_model(): Rt_process=Rt_process, latent_infections=latent_infections, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=DeterministicVariable((0,)), + observation_process=DeterministicVariable(0), ) # Sampling and fitting model 0 (with no obs for infections) @@ -101,7 +98,7 @@ def test_model_hosp_no_obs_model(): num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(272), - obs_mean=model0_samp.sampled, + obs_mean=model0_samp.sampled_admissions, n_timepoints=30, ) @@ -122,9 +119,7 @@ def test_model_hosp_with_obs_model(): Checks that the random Hospitalization model runs """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -133,29 +128,27 @@ def test_model_hosp_with_obs_model(): observed_hospitalizations = PoissonObservation() inf_hosp = DeterministicPMF( - ( - jnp.array( - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0.25, - 0.5, - 0.1, - 0.1, - 0.05, - ], - ), + jnp.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0.25, + 0.5, + 0.1, + 0.1, + 0.05, + ], ), ) @@ -172,7 +165,7 @@ def test_model_hosp_with_obs_model(): Rt_process=Rt_process, latent_infections=latent_infections, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=observed_hospitalizations, + observation_process=observed_hospitalizations, ) # Sampling and fitting model 0 (with no obs for infections) @@ -184,7 +177,7 @@ def test_model_hosp_with_obs_model(): num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(272), - observed_hospitalizations=model1_samp.sampled, + observed_hospitalizations=model1_samp.sampled_admissions, n_timepoints=30, ) @@ -205,9 +198,7 @@ def test_model_hosp_with_obs_model_weekday_phosp_2(): Checks that the random Hospitalization model runs """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -216,29 +207,27 @@ def test_model_hosp_with_obs_model_weekday_phosp_2(): observed_hospitalizations = PoissonObservation() inf_hosp = DeterministicPMF( - ( - jnp.array( - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0.25, - 0.5, - 0.1, - 0.1, - 0.05, - ], - ), + jnp.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0.25, + 0.5, + 0.1, + 0.1, + 0.05, + ], ), ) @@ -266,7 +255,7 @@ def test_model_hosp_with_obs_model_weekday_phosp_2(): Rt_process=Rt_process, latent_infections=latent_infections, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=observed_hospitalizations, + observation_process=observed_hospitalizations, ) # Sampling and fitting model 0 (with no obs for infections) @@ -278,7 +267,7 @@ def test_model_hosp_with_obs_model_weekday_phosp_2(): num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(272), - observed_hospitalizations=model1_samp.sampled, + observed_hospitalizations=model1_samp.sampled_admissions, n_timepoints=30, ) @@ -299,9 +288,7 @@ def test_model_hosp_with_obs_model_weekday_phosp(): Checks that the random Hospitalization model runs """ - gen_int = DeterministicPMF( - (jnp.array([0.25, 0.25, 0.25, 0.25]),), - ) + gen_int = DeterministicPMF(jnp.array([0.25, 0.25, 0.25, 0.25])) I0 = Infections0(I0_dist=dist.LogNormal(0, 1)) @@ -310,29 +297,27 @@ def test_model_hosp_with_obs_model_weekday_phosp(): observed_hospitalizations = PoissonObservation() inf_hosp = DeterministicPMF( - ( - jnp.array( - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0.25, - 0.5, - 0.1, - 0.1, - 0.05, - ], - ), + jnp.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0.25, + 0.5, + 0.1, + 0.1, + 0.05, + ], ), ) @@ -342,7 +327,7 @@ def test_model_hosp_with_obs_model_weekday_phosp(): weekday = weekday / weekday.sum() weekday = weekday[:31] - weekday = DeterministicVariable((weekday,)) + weekday = DeterministicVariable(weekday) hosp_report_prob_dist = jnp.array([0.9, 0.8, 0.7, 0.7, 0.6, 0.4]) hosp_report_prob_dist = jnp.tile(hosp_report_prob_dist, 10) @@ -350,9 +335,7 @@ def test_model_hosp_with_obs_model_weekday_phosp(): hosp_report_prob_dist = hosp_report_prob_dist[:31] - hosp_report_prob_dist = DeterministicVariable( - vars=(hosp_report_prob_dist,) - ) + hosp_report_prob_dist = DeterministicVariable(vars=hosp_report_prob_dist) latent_hospitalizations = HospitalAdmissions( infection_to_admission_interval=inf_hosp, @@ -369,7 +352,7 @@ def test_model_hosp_with_obs_model_weekday_phosp(): Rt_process=Rt_process, latent_infections=latent_infections, latent_hospitalizations=latent_hospitalizations, - observed_hospitalizations=observed_hospitalizations, + observation_process=observed_hospitalizations, ) # Sampling and fitting model 0 (with no obs for infections) @@ -377,12 +360,18 @@ def test_model_hosp_with_obs_model_weekday_phosp(): with npro.handlers.seed(rng_seed=np.random.randint(1, 600)): model1_samp = model1.sample(n_timepoints=30) + obs = jnp.hstack( + [jnp.repeat(jnp.nan, 5), model1_samp.sampled_admissions[5:]] + ) + + # Running with padding model1.run( num_warmup=500, num_samples=500, rng_key=jax.random.PRNGKey(272), - observed_hospitalizations=model1_samp.sampled, + observed_hospitalizations=obs, n_timepoints=30, + padding=5, ) inf = model1.spread_draws(["predicted_hospitalizations"])