Commit 3a3b583a authored by Chris Jewell's avatar Chris Jewell
Browse files

Removed epsilon background term. It was a bad idea!

parent 15938bd6
......@@ -86,9 +86,7 @@ class CovidUKODE: # TODO: add background case importation rate to the UK, e.g.
self.m_select = np.int64((self.times >= holidays[0]) &
(self.times < holidays[1]))
self.max_t = self.m_select.shape[0] - 1
self.bg_select = tf.constant(self.times < bg_max_t, dtype=dtype)
self.bg_max_t = tf.convert_to_tensor(bg_max_t, dtype=dtype)
self.solver = tode.DormandPrince()
def make_h(self, param):
......@@ -101,12 +99,11 @@ class CovidUKODE: # TODO: add background case importation rate to the UK, e.g.
# holiday status as their nearest neighbors in the desired range.
t_idx = tf.clip_by_value(tf.cast(t, tf.int64), 0, self.max_t)
m_switch = tf.gather(self.m_select, t_idx)
epsilon = param['epsilon'] * tf.gather(self.bg_select, t_idx)
infec_rate = param['beta1'] * (
tf.gather(self.M.matvec(I), m_switch) +
param['beta2'] * self.Kbar * self.C.matvec(I / self.N_sum))
infec_rate = S / self.N * (infec_rate + epsilon)
infec_rate = S / self.N * infec_rate
dS = -infec_rate
dE = infec_rate - param['nu'] * E
......
......@@ -68,7 +68,21 @@ def phe_death_hosp_to_death(filename, date_range=['2020-02-02', '2020-03-21']):
return data.dropna(axis=0)
def phe_linelist_timeseries(filename, date_range):
linelist = pd.read_excel(filename)
linelist['age_group'] = (linelist['Age_in_Years'] // 5) # id of 5-year age group
cols = ['patienttable_Specimen_Date', 'UTLA Code', 'age_group']
missing = linelist[cols].apply(lambda x: np.sum(x.isna()), axis=0)
missingness = {k: v for k, v in zip(cols, missing)}
grouped = linelist.groupby(cols)
case_counts = grouped.size()
print("Missing:", missingness)
return case_counts
def spatial_report(cases, utla_geom):
pass
if __name__=='__main__':
......
......@@ -103,15 +103,13 @@ if __name__ == '__main__':
def logp(par):
p = param
p['epsilon'] = par[0]
p['beta1'] = par[1]
p['gamma'] = par[2]
epsilon_logp = tfd.Gamma(concentration=tf.constant(8., tf.float64), rate=tf.constant(1000., tf.float64)).log_prob(p['epsilon'])
p['beta1'] = par[0]
p['gamma'] = par[1]
beta_logp = tfd.Gamma(concentration=tf.constant(1., tf.float64), rate=tf.constant(1., tf.float64)).log_prob(p['beta1'])
gamma_logp = tfd.Gamma(concentration=tf.constant(100., tf.float64), rate=tf.constant(400., tf.float64)).log_prob(p['gamma'])
t, sim, solve = simulator.simulate(p, state_init)
y_logp = covid19uk_logp(y_incr, sim, 0.1)
logp = epsilon_logp + beta_logp + gamma_logp + tf.reduce_sum(y_logp)
logp = beta_logp + gamma_logp + tf.reduce_sum(y_logp)
return logp
def trace_fn(_, pkr):
......@@ -122,7 +120,7 @@ if __name__ == '__main__':
unconstraining_bijector = [tfb.Exp()]
initial_mcmc_state = np.array([0.008, 0.05, 0.25], dtype=np.float64)
initial_mcmc_state = np.array([0.05, 0.25], dtype=np.float64)
print("Initial log likelihood:", logp(initial_mcmc_state))
@tf.function(autograph=False, experimental_compile=True)
......@@ -141,7 +139,7 @@ if __name__ == '__main__':
joint_posterior = tf.zeros([0] + list(initial_mcmc_state.shape), dtype=DTYPE)
scale = np.diag([0.00001, 0.00001, 0.00001])
scale = np.diag([0.1, 0.1])
overall_start = time.perf_counter()
num_covariance_estimation_iterations = 50
......@@ -177,7 +175,6 @@ if __name__ == '__main__':
fig, ax = plt.subplots(1, 3)
ax[0].plot(joint_posterior[:, 0])
ax[1].plot(joint_posterior[:, 1])
ax[2].plot(joint_posterior[:, 2])
plt.show()
print(f"Posterior mean: {np.mean(joint_posterior, axis=0)}")
......
......@@ -15,6 +15,7 @@ from covid.model import CovidUKODE
from covid.rdata import load_age_mixing, load_mobility_matrix, load_population
from covid.util import sanitise_settings, sanitise_parameter, seed_areas, doubling_time
DTYPE = np.float64
def save_sims(sims, la_names, age_groups, filename):
f = h5py.File(filename, 'w')
......@@ -44,12 +45,18 @@ if __name__ == '__main__':
N, n_names = load_population(config['data']['population_size'])
K_tt = K_tt.astype(DTYPE)
K_hh = K_hh.astype(DTYPE)
T = T.astype(DTYPE)
N = N.astype(DTYPE)
param = sanitise_parameter(config['parameter'])
param['epsilon'] = 0.0
settings = sanitise_settings(config['settings'])
case_reports = pd.read_csv(config['data']['reported_cases'])
case_reports['DateVal'] = pd.to_datetime(case_reports['DateVal'])
case_reports = case_reports[case_reports['DateVal'] >= '2020-02-19']
date_range = [case_reports['DateVal'].min(), case_reports['DateVal'].max()]
y = case_reports['CumCases'].to_numpy()
y_incr = y[1:] - y[-1:]
......@@ -67,12 +74,11 @@ if __name__ == '__main__':
state_init = simulator.create_initial_state(init_matrix=seeding)
@tf.function
def prediction(epsilon, beta, gamma):
def prediction(beta, gamma):
sims = tf.TensorArray(tf.float64, size=beta.shape[0])
R0 = tf.TensorArray(tf.float64, size=beta.shape[0])
for i in tf.range(beta.shape[0]):
p = param
p['epsilon'] = epsilon[i]
p['beta1'] = beta[i]
p['gamma'] = gamma[i]
t, sim, solver_results = simulator.simulate(p, state_init)
......@@ -83,7 +89,7 @@ if __name__ == '__main__':
draws = pi_beta.numpy()[np.arange(5000, pi_beta.shape[0], 10), :]
with tf.device('/CPU:0'):
sims, R0 = prediction(draws[:, 0], draws[:, 1], draws[:, 2])
sims, R0 = prediction(draws[:, 0], draws[:, 1])
sims = tf.stack(sims) # shape=[n_sims, n_times, n_states, n_metapops]
save_sims(sims, la_names, age_groups, 'pred_2020-03-15.h5')
......@@ -140,5 +146,5 @@ if __name__ == '__main__':
print("Doubling time:", dub_ci)
# Infectious period
ip = tfs.percentile(1./pi_beta[3000:, 2], q=[2.5, 50, 97.5])
ip = tfs.percentile(1./pi_beta[3000:, 1], q=[2.5, 50, 97.5])
print("Infectious period:", ip)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment