diff --git a/src/penn_chime/models.py b/src/penn_chime/models.py index da7311fe..6582b1c3 100644 --- a/src/penn_chime/models.py +++ b/src/penn_chime/models.py @@ -9,7 +9,7 @@ from datetime import date, datetime, timedelta from logging import INFO, basicConfig, getLogger from sys import stdout -from typing import Dict, Generator, Tuple, Sequence,Optional +from typing import Dict, Generator, Tuple, Sequence, Optional import numpy as np import pandas as pd @@ -66,14 +66,13 @@ def __init__(self, p: Parameters): intrinsic_growth_rate = get_growth_rate(p.doubling_time) self.beta = get_beta(intrinsic_growth_rate, gamma, self.susceptible, 0.0) + self.beta_t = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, p.relative_contact_rate) self.i_day = 0 # seed to the full length - self.beta_t = self.beta - self.run_projection(p) + self.run_projection(p, [(self.beta, p.n_days)]) self.i_day = i_day = int(get_argmin_ds(self.census_df, p.current_hospitalized)) - self.beta_t = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, p.relative_contact_rate) - self.run_projection(p) + self.run_projection(p, self.gen_policy(p)) logger.info('Set i_day = %s', i_day) p.date_first_hospitalized = p.current_date - timedelta(days=i_day) @@ -100,7 +99,7 @@ def __init__(self, p: Parameters): self.beta = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, 0.0) self.beta_t = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, p.relative_contact_rate) - self.run_projection(p) + self.run_projection(p, self.gen_policy(p)) loss = self.get_loss() losses[i] = loss @@ -109,7 +108,7 @@ def __init__(self, p: Parameters): intrinsic_growth_rate = get_growth_rate(p.doubling_time) self.beta = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, 0.0) self.beta_t = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, p.relative_contact_rate) - self.run_projection(p) + self.run_projection(p, self.gen_policy(p)) self.population = p.population else: @@ -146,18 +145,35 @@ def __init__(self, p: Parameters): self.daily_growth_rate = get_growth_rate(p.doubling_time) self.daily_growth_rate_t = get_growth_rate(self.doubling_time_t) - def run_projection(self, p): + def gen_policy(self, p: Parameters) -> Sequence[Tuple[float, int]]: + if p.mitigation_date is not None: + mitigation_day = -(p.current_date - p.mitigation_date).days + else: + mitigation_day = 0 + + total_days = self.i_day + p.n_days + + if mitigation_day < -self.i_day: + mitigation_day = -self.i_day + + pre_mitigation_days = self.i_day + mitigation_day + post_mitigation_days = total_days - pre_mitigation_days + + return [ + (self.beta, pre_mitigation_days), + (self.beta_t, post_mitigation_days), + ] + + def run_projection(self, p: Parameters, policy: Sequence[Tuple[float, int]]): self.raw_df = sim_sir_df( self.susceptible, self.infected, p.recovered, self.gamma, -self.i_day, - self.beta, - self.i_day, - self.beta_t, - p.n_days + policy ) + self.dispositions_df = build_dispositions_df(self.raw_df, self.rates, p.market_share, p.current_date) self.admits_df = build_admits_df(self.dispositions_df) self.census_df = build_census_df(self.admits_df, self.days) @@ -221,7 +237,7 @@ def sir( def gen_sir( - s: float, i: float, r: float, gamma: float, i_day: int, *args + s: float, i: float, r: float, gamma: float, i_day: int, policies: Sequence[Tuple[float, int]] ) -> Generator[Tuple[int, float, float, float], None, None]: """Simulate SIR model forward in time yielding tuples. Parameter order has changed to allow multiple (beta, n_days) @@ -230,8 +246,7 @@ def gen_sir( s, i, r = (float(v) for v in (s, i, r)) n = s + i + r d = i_day - while args: - beta, n_days, *args = args + for beta, n_days in policies: for _ in range(n_days): yield d, s, i, r s, i, r = sir(s, i, r, beta, gamma, n) @@ -241,11 +256,11 @@ def gen_sir( def sim_sir_df( s: float, i: float, r: float, - gamma: float, i_day: int, *args + gamma: float, i_day: int, policies: Sequence[Tuple[float, int]] ) -> pd.DataFrame: """Simulate the SIR model forward in time.""" return pd.DataFrame( - data=gen_sir(s, i, r, gamma, i_day, *args), + data=gen_sir(s, i, r, gamma, i_day, policies), columns=("day", "susceptible", "infected", "recovered"), ) diff --git a/src/penn_chime/parameters.py b/src/penn_chime/parameters.py index d9da0472..d6c03a08 100644 --- a/src/penn_chime/parameters.py +++ b/src/penn_chime/parameters.py @@ -55,6 +55,7 @@ def __init__( hospitalized: Disposition, icu: Disposition, relative_contact_rate: float, + mitigation_date: Optional[date] = None, ventilated: Disposition, current_date: date = date.today(), date_first_hospitalized: Optional[date] = None, @@ -68,7 +69,6 @@ def __init__( region: Optional[Regions] = None, ): self.current_hospitalized = Positive(value=current_hospitalized) - self.relative_contact_rate = Rate(value=relative_contact_rate) Rate(value=hospitalized.rate), Rate(value=icu.rate), Rate(value=ventilated.rate) StrictlyPositive(value=hospitalized.days), StrictlyPositive(value=icu.days), @@ -92,6 +92,9 @@ def __init__( self.date_first_hospitalized = OptionalDate(value=date_first_hospitalized) self.doubling_time = OptionalStrictlyPositive(value=doubling_time) + self.relative_contact_rate = Rate(value=relative_contact_rate) + self.mitigation_date = OptionalDate(value=mitigation_date) + self.infectious_days = StrictlyPositive(value=infectious_days) self.market_share = Rate(value=market_share) self.max_y_axis = OptionalStrictlyPositive(value=max_y_axis) diff --git a/src/penn_chime/presentation.py b/src/penn_chime/presentation.py index 6492e1bd..0f50f424 100644 --- a/src/penn_chime/presentation.py +++ b/src/penn_chime/presentation.py @@ -11,6 +11,7 @@ CHANGE_DATE, DATE_FORMAT, DOCS_URL, + EPSILON, FLOAT_INPUT_MIN, FLOAT_INPUT_STEP, ) @@ -207,6 +208,10 @@ def display_sidebar(st, d: Parameters) -> Parameters: st_obj, "Date of first hospitalized case - Enter this date to have chime estimate the initial doubling time", value=d.date_first_hospitalized, ) + mitigation_date_input = DateInput( + st_obj, "Date of social distancing measures effect (may be delayed from implementation)", + value=d.mitigation_date + ) relative_contact_pct_input = PercentInput( st_obj, "Social distancing (% reduction in social contact going forward)", @@ -312,7 +317,15 @@ def display_sidebar(st, d: Parameters) -> Parameters: doubling_time = doubling_time_input() date_first_hospitalized = None - relative_contact_rate = relative_contact_pct_input() + if st.sidebar.checkbox( + "Social distancing measures have been implemented", + value=(d.relative_contact_rate > EPSILON) + ): + mitigation_date = mitigation_date_input() + relative_contact_rate = relative_contact_pct_input() + else: + mitigation_date = None + relative_contact_rate = EPSILON st.sidebar.markdown( "### Severity Parameters [ℹ]({docs_url}/what-is-chime/parameters#severity-parameters)".format( @@ -346,6 +359,7 @@ def display_sidebar(st, d: Parameters) -> Parameters: hospitalized=Disposition(hospitalized_rate, hospitalized_days), icu=Disposition(icu_rate, icu_days), relative_contact_rate=relative_contact_rate, + mitigation_date=mitigation_date, ventilated=Disposition(ventilated_rate, ventilated_days), current_date=current_date, date_first_hospitalized=date_first_hospitalized, diff --git a/src/penn_chime/settings.py b/src/penn_chime/settings.py index 0bd12989..a9ccb347 100644 --- a/src/penn_chime/settings.py +++ b/src/penn_chime/settings.py @@ -16,6 +16,7 @@ def get_defaults(): infectious_days=14, market_share=0.15, n_days=100, + mitigation_date=date.today(), relative_contact_rate=0.3, ventilated=Disposition(0.005, 10), ) diff --git a/tests/conftest.py b/tests/conftest.py index e822d91a..b7cf01f2 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -51,6 +51,7 @@ def DEFAULTS(): doubling_time=4.0, n_days=60, market_share=0.15, + mitigation_date=datetime(year=2020, month=3, day=28), relative_contact_rate=0.3, hospitalized=Disposition(0.025, 7), icu=Disposition(0.0075, 9), @@ -65,6 +66,7 @@ def param(): current_hospitalized=100, doubling_time=6.0, market_share=0.05, + mitigation_date=datetime(year=2020, month=3, day=28), relative_contact_rate=0.15, population=500000, hospitalized=Disposition(0.05, 7), @@ -81,6 +83,7 @@ def halving_param(): current_hospitalized=100, doubling_time=6.0, market_share=0.05, + mitigation_date=datetime(year=2020, month=3, day=28), relative_contact_rate=0.7, population=500000, hospitalized=Disposition(0.05, 7), diff --git a/tests/penn_chime/test_models.py b/tests/penn_chime/test_models.py index a8c4129e..d5e6de20 100644 --- a/tests/penn_chime/test_models.py +++ b/tests/penn_chime/test_models.py @@ -3,11 +3,13 @@ import pytest import pandas as pd import numpy as np +from datetime import timedelta from src.penn_chime.models import ( sir, sim_sir_df, get_growth_rate, + SimSirModel, ) from src.penn_chime.constants import EPSILON @@ -64,7 +66,7 @@ def test_sim_sir(): Rounding to move fast past decimal place issues """ raw_df = sim_sir_df( - 5, 6, 7, 0.1, 0, 0.1, 40, # s # i # r # gamma # i_day # beta1 # n_days1 + 5, 6, 7, 0.1, 0, [(0.1, 40)], # s # i # r # gamma # i_day # beta1 # n_days1 ) first = raw_df.iloc[0, :] @@ -100,6 +102,20 @@ def test_model(model, param): assert model.r_t == 2.307298374881539 assert model.r_naught == 2.7144686763312222 assert model.doubling_time_t == 7.764405988534983 + assert model.i_day == 43 + + +def test_model_first_hosp_fit(param): + param.date_first_hospitalized = param.current_date - timedelta(days=43) + param.doubling_time = None + + my_model = SimSirModel(param) + + assert my_model.intrinsic_growth_rate == 0.12246204830937302 + assert abs(my_model.beta - 4.21501347256401e-07) < EPSILON + assert my_model.r_t == 2.307298374881539 + assert my_model.r_naught == 2.7144686763312222 + assert my_model.doubling_time_t == 7.764405988534983 def test_model_raw_start(model, param):