Merge pull request #519 from CodeForPhilly/439-doubling-time-fit-mitigation

Fit the back-cast first admission date with mitigation in mind

Author: mdbecker

src/penn_chime/models.py

@@ -64,21 +64,40 @@ def __init__(self, p: Parameters):
             logger.info('Using doubling_time: %s', p.doubling_time)
 
             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
-            raw = self.run_projection(p, [(self.beta, p.n_days)])
-            self.i_day = i_day = int(get_argmin_ds(raw["census_hospitalized"], p.current_hospitalized))
+            if p.mitigation_date is None:
+                self.i_day = 0 # seed to the full length
+                raw = self.run_projection(p, [(self.beta, p.n_days)])
+                self.i_day = i_day = int(get_argmin_ds(raw["census_hospitalized"], p.current_hospitalized))
 
-            self.raw = self.run_projection(p, self.gen_policy(p))
+                self.raw = self.run_projection(p, self.gen_policy(p))
+
+                logger.info('Set i_day = %s', i_day)
+            else:
+                projections = {}
+                best_i_day = -1
+                best_i_day_loss = float('inf')
+                for i_day in range(p.n_days):
+                    self.i_day = i_day
+                    raw = self.run_projection(p, self.gen_policy(p))
+
+                    # Don't fit against results that put the peak before the present day
+                    if raw["census_hospitalized"].argmax() < i_day:
+                        continue
+
+                    loss = get_loss(raw["census_hospitalized"][i_day], p.current_hospitalized)
+                    if loss < best_i_day_loss:
+                        best_i_day_loss = loss
+                        best_i_day = i_day
+                        self.raw = raw
+
+                self.i_day = best_i_day
 
-            logger.info('Set i_day = %s', i_day)
-            p.date_first_hospitalized = p.current_date - timedelta(days=i_day)
             logger.info(
                 'Estimated date_first_hospitalized: %s; current_date: %s; i_day: %s',
-                p.date_first_hospitalized,
+                p.current_date - timedelta(days=self.i_day),
                 p.current_date,
                 self.i_day)