Skip to content

Commit 5c156ec

Browse files
committed
Fit first admission date against given doubling time - still needs to prevent post-peak results
1 parent 4ba3bed commit 5c156ec

File tree

1 file changed

+23
-9
lines changed

1 file changed

+23
-9
lines changed

src/penn_chime/models.py

Lines changed: 23 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -64,21 +64,35 @@ def __init__(self, p: Parameters):
6464
logger.info('Using doubling_time: %s', p.doubling_time)
6565

6666
intrinsic_growth_rate = get_growth_rate(p.doubling_time)
67-
6867
self.beta = get_beta(intrinsic_growth_rate, gamma, self.susceptible, 0.0)
6968
self.beta_t = get_beta(intrinsic_growth_rate, self.gamma, self.susceptible, p.relative_contact_rate)
7069

71-
self.i_day = 0 # seed to the full length
72-
raw = self.run_projection(p, [(self.beta, p.n_days)])
73-
self.i_day = i_day = int(get_argmin_ds(raw["census_hospitalized"], p.current_hospitalized))
74-
75-
self.raw = self.run_projection(p, self.gen_policy(p))
70+
if p.mitigation_date is None:
71+
self.i_day = 0 # seed to the full length
72+
raw = self.run_projection(p, [(self.beta, p.n_days)])
73+
self.i_day = i_day = int(get_argmin_ds(raw["census_hospitalized"], p.current_hospitalized))
74+
75+
self.raw = self.run_projection(p, self.gen_policy(p))
76+
77+
logger.info('Set i_day = %s', i_day)
78+
else:
79+
projections = {}
80+
best_i_day = -1
81+
best_i_day_loss = float('inf')
82+
for i_day in range(p.n_days):
83+
self.i_day = i_day
84+
raw = self.run_projection(p, self.gen_policy(p))
85+
loss = get_loss(raw["census_hospitalized"][i_day], p.current_hospitalized)
86+
if loss < best_i_day_loss:
87+
best_i_day_loss = loss
88+
best_i_day = i_day
89+
self.raw = raw
90+
91+
self.i_day = best_i_day
7692

77-
logger.info('Set i_day = %s', i_day)
78-
p.date_first_hospitalized = p.current_date - timedelta(days=i_day)
7993
logger.info(
8094
'Estimated date_first_hospitalized: %s; current_date: %s; i_day: %s',
81-
p.date_first_hospitalized,
95+
p.current_date - timedelta(days=self.i_day),
8296
p.current_date,
8397
self.i_day)
8498

0 commit comments

Comments
 (0)