Allow the user to prescribe mitigation date

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"),
     )
 

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)

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,

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),
     )

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),

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):