Simulation¶
In [1]:
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.cm as cm
import pandas as pd
import matplotlib.ticker as mticker
import matplotlib
from matplotlib.pyplot import figure
from matplotlib.ticker import FormatStrFormatter
import numpy as np
from pathlib import Path
Path("./content/sample_data").mkdir(parents=True, exist_ok=True)
Lockdown only delays the virus spread¶
In [2]:
# Set policy efficiencies
test_cap = 0.001
hygiene = 0.30
distancing = 0.74
lockdown = 0.96
quarantine = 0.96
# ASSUMPTIONS:
# Virus parameters
transmission_days = 16.33/2.40
exposed_days = 5.27
recovery_days = 16.33
r0 = 2.64
transmission_rate = (1/ transmission_days)
incubation_rate = (1/ exposed_days)
recovery_rate = (1/ recovery_days)
# Country parameters
hospital_cap = 6/10000
death_rate_with_med = 0.6 # 60% chance to die with ICU treatment
death_rate_without_med = 1 # 100% chance to die without ICU treatment
general_dr = 0.02 # 2% death rate
icu_case_rate = 0.06 # 6% need ICU treatment
ser_case_rate = 0.15 # 15% of serious cases
# Initialize
susceptible0 = 0.99
infected0 = 0.01
exposed0 = 0
recovered0 = 0
population = 1000000
policy_start = 30
policy_end = 60
days = 300
def virus_model(days,policy):
if policy == "hygiene":
pol = hygiene
elif policy == "distancing":
pol = distancing
elif policy == "lockdown":
pol = lockdown
elif policy == "quarantine":
pol = quarantine
elif policy == "None":
pol = None
# set values for the first date of the simulation
susceptible = susceptible0
infected = infected0
exposed = exposed0
recovered = recovered0
dead = 0
# Scale to population size
susc= [population * susceptible0]
exp = [population * exposed0]
inf = [population * infected0]
rec = [population * recovered0]
re = [r0]
d = [0]
for day in range(days):
if policy_start < day < policy_end and pol != None:
transmission_rate = (1/ transmission_days) * (1-pol)
else:
transmission_rate = (1/ transmission_days)
# Simulate flow between compartments
S2E = susceptible * infected * transmission_rate
E2I = exposed * incubation_rate
I2D = general_dr*infected/recovery_days
# Account different death rates depending on condition
if (icu_case_rate*infected < hospital_cap):
I2D += icu_case_rate*infected*death_rate_with_med/recovery_days
else:
I2D += hospital_cap*death_rate_with_med/recovery_days + (icu_case_rate*infected - hospital_cap) * death_rate_without_med/recovery_days
I2R = (1-general_dr)*infected*recovery_rate
# Update the compartments
exposed += S2E
susceptible -= S2E
infected += E2I
exposed -= E2I
recovered += I2R
infected -= I2R
infected -= I2D
dead += I2D
re.append((transmission_rate/ recovery_rate)*susceptible)
susc.append(population* susceptible)
exp.append(population* exposed)
inf.append(population* infected)
rec.append(population* recovered)
d.append(population* dead)
return susc, exp, inf, rec, d, re
result_h = virus_model(days,'hygiene')
result_d = virus_model(days,'distancing')
result_l = virus_model(days,'lockdown')
result_q = virus_model(days,'quarantine')
result_0 = virus_model(days,'None')
m = [m for m in range(days+1)]
# Plot the graphs
# Here we plot only lockdown
sns.set_style("ticks")
matplotlib.rcParams['figure.figsize'] = 8, 6
sns.despine()
formatter = mticker.ScalarFormatter(useMathText=True)
f = mticker.ScalarFormatter(useOffset=False, useMathText=True)
g = lambda x,pos : "${}$".format(f._formatSciNotation('%10e' % x))
plots = ['Susceptible','Exposed', 'Infected','Recovered','Death cases']
abbrev = ['S','E', 'I','R','D']
for num,policy in enumerate(plots):
plt.plot(m, result_0[num], label="No policy", color ='grey',linewidth=2.5)
plt.plot(m, result_l[num], label="Lockdown", color ='green',linewidth=2.5)
plt.axvline(x=policy_start,linestyle='--', linewidth=1.5, label='Policy start')
plt.axvline(x=policy_end,linestyle='--', linewidth=1.5, color='red', label='Policy end')
plt.tick_params(labelsize=22)
plt.ylabel(policy, fontsize=22)
plt.xlabel('Days', fontsize=22)
plt.legend(fontsize=22)
sns.despine()
plt.gca().yaxis.set_major_formatter(mticker.FuncFormatter(g))
plt.savefig('./content/sample_data/{}.pdf'.format(abbrev[num]),bbox_inches='tight')
plt.show()
Policy Efficacy check¶
Find optimal policy¶
In [3]:
# Set policy efficiencies
test_cap = 0.001
hygiene = 0.30
distancing = 0.74
lockdown = 0.96
quarantine = 0.96
# ASSUMPTIONS:
# Virus parameters
transmission_days = 16.33/2.40
exposed_days = 5.27
recovery_days = 16.33
r0 = 2.64
transmission_rate0 = (1/ transmission_days)
incubation_rate = (1/ exposed_days)
recovery_rate = (1/ recovery_days)
# Country parameters
hospital_cap = 6/10000
death_rate_with_med = 0.6 # 60% chance to die with ICU treatment
death_rate_without_med = 1 # 100% chance to die without ICU treatment
general_dr = 0.02 # 2% death rate
icu_case_rate = 0.06 # 6% need ICU treatment
ser_case_rate = 0.15 # 15% of serious cases
# Initialize
susceptible0 = 0.99
infected0 = 0.01
exposed0 = 0
recovered0 = 0
population = 1000000
# Policy costs:
lockdown_cost = 0.1 # 10% of yearly GDP
quarantine_cost = 0.05 # 5% of yearly GDP
mask_cost = 0.002 # 2$
distancing_cost = 0.05 # 5% of yearly GDP
def evaluate(inf, rec, d, policy, cont_tracing):
''' Calculate penalty points '''
dd_point = -7000 # death
inf_point = -0.3 # infection
ct_point = -6.4 # contact_tracing
result = sum(inf)*inf_point + d[-1]*dd_point + sum(policy) + sum(cont_tracing)*ct_point
return result
def virus_model(hygiene_list, distancing_list,lockdown_list,quarantine_list):
assert len(hygiene_list) == len(distancing_list) == len(lockdown_list) == len(quarantine_list)
duration = len(hygiene_list)
susceptible = susceptible0
infected = infected0
exposed = exposed0
recovered = recovered0
dead = 0
susc= [population * susceptible0]
exp = [population * exposed0]
inf = [population * infected0]
rec = [population * recovered0]
cont_tracing = [] #for calculating contact tracing price
re = [r0]
d = [0]
policy = []
for day in range(duration):
# mitigate transmission rate (Re) by applying policies:
transmission_rate = transmission_rate0 * (1- hygiene_list[day]*hygiene) * (1- distancing_list[day]*distancing) * (1- lockdown_list[day]*lockdown) * (1- quarantine_list[day]*quarantine)
# calculate loss per day for each policy
ld_point = - population * lockdown_cost * 30 /365 # account for 30k$ of yearly GDP
qr_point = - population * quarantine_cost * 30 /365 # account for 30k$ of yearly GDP
hg_point = - population * mask_cost
dst_point = - population * distancing_cost * 30 /365 # account for 30k$ of yearly GDP
policy_point = hg_point * hygiene_list[day] + dst_point * distancing_list[day] + ld_point*lockdown_list[day] + qr_point*quarantine_list[day]
# Simulate the flow between compartments
S2E = susceptible * infected * transmission_rate
E2I = exposed * incubation_rate
I2D = general_dr*infected/recovery_days
# Account different death rates depending on condition
if (icu_case_rate*infected < hospital_cap):
I2D += icu_case_rate*infected*death_rate_with_med/recovery_days
else:
I2D += hospital_cap*death_rate_with_med/recovery_days + (icu_case_rate*infected - hospital_cap) * death_rate_without_med/recovery_days
I2R = (1-general_dr)*infected*recovery_rate
# Update the compartments
exposed += S2E
susceptible -= S2E
infected += E2I
exposed -= E2I
recovered += I2R
infected -= I2R
infected -= I2D
dead += I2D
re.append((transmission_rate/ recovery_rate)*susceptible)
susc.append(population* susceptible)
exp.append(population* exposed)
inf.append(population* infected)
rec.append(population* recovered)
d.append(population* dead)
policy.append(policy_point)
# Append number of exposed people for each date for contact-tracing loss estimation in evaluate()
cont_tracing.append(S2E*quarantine_list[day]*population)
return susc, exp, inf, rec, d, re, policy, cont_tracing
# h1, h2, h3, d1, d2, d3, l1, l2, l3, q1, q2, q3
# numbers indicate month. Ex: l2 = 0.5 -> half lockdown is applied in 2nd month
result = np.zeros((3,3,3,3,3,3,3,3,3,3,3,3), dtype = float)
for index, value in np.ndenumerate(result):
h1, h2, h3, d1, d2, d3, l1, l2, l3, q1, q2, q3 = index
# Filter out unreasonable cases: Lockdown and social distancing in the same month
# OR Distacing and (distabncing + tracing) in the same month
if (d1!=0 and q1!=0) or (d2!=0 and q2!=0) or (d3!=0 and q3!=0) or \
(l1!=0 and (q1!=0 or d1!=0)) or (l2!=0 and (d2!=0 or q2!=0)) or (l3!=0 and (d3!=0 or q3!=0)):
result[h1, h2, h3, d1, d2, d3, l1, l2, l3, q1, q2, q3] = None
continue
# Divide by two to get the policy strength: 0, 50%, 100%
hygiene_list = [h1/2]*30 + [h2/2]*30 + [h3/2]*30
distancing_list = [d1/2]*30 + [d2/2]*30 + [d3/2]*30
lockdown_list = [l1/2]*30 + [l2/2]*30 + [l3/2]*30
quarantine_list = [q1/2]*30 + [q2/2]*30 + [q3/2]*30
# Get the compartments
susc, exp, inf, rec, d, re, policy,cont_tracing = virus_model(hygiene_list, distancing_list,lockdown_list,quarantine_list)
# Evaluate the loss
point = evaluate(inf, rec, d, policy,cont_tracing)
result[h1, h2, h3, d1, d2, d3, l1, l2, l3, q1, q2, q3] = point
#no policy case
print(result[0,0,0,0,0,0,0,0,0,0,0,0])
-197972647.11065745
In [4]:
# Optimal policy loss
loss = np.nanmax(result)
arg = np.unravel_index(np.nanargmax(result), result.shape)
print(loss, arg)
-4525982.639901729 (2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2)
Plot outcomes of different combinations¶
In [5]:
index = []
index.append([0,0,0,0,0,0,0,0,0,0,0,0])
index.append([2,2,2,0,0,0,0,0,0,0,0,0]) # hygiene
index.append([0,0,0,2,2,2,0,0,0,0,0,0]) # distancing
index.append([0,0,0,0,0,0,2,2,2,0,0,0]) # lockdown
index.append([0,0,0,0,0,0,0,0,0,2,2,2]) # quarantine
index.append([2,0,0,0,0,0,0,0,0,2,2,2]) # optimal combination
In [6]:
ind = np.unique(np.asarray(index), axis = 0)
ind.shape
Out[6]:
(6, 12)
In [7]:
# Set policy efficiencies
test_cap = 0.001
hygiene = 0.30
distancing = 0.74
lockdown = 0.96
quarantine = 0.96
# ASSUMPTIONS:
# Virus parameters
transmission_days = 16.33/2.40
exposed_days = 5.27
recovery_days = 16.33
r0 = 2.64
transmission_rate0 = (1/ transmission_days)
incubation_rate = (1/ exposed_days)
recovery_rate = (1/ recovery_days)
# Country parameters
hospital_cap = 6/10000
death_rate_with_med = 0.6 # 60% chance to die with ICU treatment
death_rate_without_med = 1 # 100% chance to die without ICU treatment
general_dr = 0.02 # 2% death rate
icu_case_rate = 0.06 # 6% need ICU treatment
ser_case_rate = 0.15 # 15% of serious cases
# Initialize
susceptible0 = 0.99
infected0 = 0.01
exposed0 = 0
recovered0 = 0
population = 1000000
# Policy costs:
lockdown_cost = 0.1 # 10% of yearly GDP
quarantine_cost = 0.05 # 5% of yearly GDP
mask_cost = 0.002 # 2$
distancing_cost = 0.05 # 5% of yearly GDP
def evaluate(inf, rec, d, policy, cont_tracing):
''' Calculate penalty points '''
dd_point = -7000
inf_point = -0.3
contact_tracing_point = -6.4
result = [-1000*(sum(inf[:i])*inf_point + d[i]*dd_point + sum(policy[:i])+ contact_tracing_point*sum(cont_tracing[:i])) for i in range(len(d))]
return result
def virus_model(hygiene_list, distancing_list,lockdown_list,quarantine_list):
assert len(hygiene_list) == len(distancing_list) == len(lockdown_list) == len(quarantine_list)
duration = len(hygiene_list)
susceptible = susceptible0
infected = infected0
exposed = exposed0
recovered = recovered0
dead = 0
susc= [population * susceptible0]
exp = [population * exposed0]
inf = [population * infected0]
rec = [population * recovered0]
cont_tracing = []
re = [r0]
d = [0]
policy = []
testing = 0
for day in range(duration):
# mitigate transmission rate (Re) by applying policies:
transmission_rate = transmission_rate0 * (1- hygiene_list[day]*hygiene) * (1- distancing_list[day]*distancing) * (1- lockdown_list[day]*lockdown) * (1- quarantine_list[day]*quarantine)
# calculate loss per day for each policy
ld_point = - population * lockdown_cost * 30 /365 # account for 30k$ of yearly GDP
qr_point = - population * quarantine_cost * 30 /365 # account for 30k$ of yearly GDP
hg_point = - population * mask_cost
dst_point = - population * distancing_cost * 30 /365 # account for 30k$ of yearly GDP
policy_point = hg_point * hygiene_list[day] + dst_point * distancing_list[day] + ld_point*lockdown_list[day] + qr_point*quarantine_list[day]
# Simulate the flow between compartments
S2E = susceptible * infected * transmission_rate
E2I = exposed * incubation_rate
I2D = general_dr*infected/recovery_days
# Account different death rates depending on condition
if (icu_case_rate*infected < hospital_cap):
I2D += icu_case_rate*infected*death_rate_with_med/recovery_days
else:
I2D += hospital_cap*death_rate_with_med/recovery_days + (icu_case_rate*infected - hospital_cap) * death_rate_without_med/recovery_days
I2R = (1-general_dr)*infected*recovery_rate
# Update the compartments
exposed += S2E
susceptible -= S2E
infected += E2I
exposed -= E2I
recovered += I2R
infected -= I2R
infected -= I2D
dead += I2D
# Append number of exposed people for each date for contact-tracing loss estimation in evaluate()
re.append((transmission_rate/ recovery_rate)*susceptible)
susc.append(population* susceptible)
exp.append(population* exposed)
inf.append(population* infected)
rec.append(population* recovered)
d.append(population* dead)
policy.append(policy_point)
cont_tracing.append(S2E*quarantine_list[day]*population)
return susc, exp, inf, rec, d, re, policy, cont_tracing
In [8]:
# Plot
figure(figsize=(10, 6), dpi=80)
legends = ['No policy (worst)', 'Hygiene & Mask', 'Distancing', 'Lockdown', 'Distancing + tracing', 'Optimal policy']
al = [0.7]*5+[1]
col = cm.rainbow(np.linspace(0, 1, 6))
plt.yscale('log')
for i,indice in enumerate(index):
h1, h2, h3, d1, d2, d3, l1, l2, l3, q1, q2, q3 = indice
hygiene_list = [h1/2]*30 + [h2/2]*30 + [h3/2]*30
distancing_list = [d1/2]*30 + [d2/2]*30 + [d3/2]*30
lockdown_list = [l1/2]*30 + [l2/2]*30 + [l3/2]*30
quarantine_list = [q1/2]*30 + [q2/2]*30 + [q3/2]*30
susc, exp, inf, rec, d, re, policy,cont_tracing = virus_model(hygiene_list, distancing_list,lockdown_list,quarantine_list)
point = evaluate(inf, rec, d, policy,cont_tracing)
m = [m for m in range(91)]
plt.plot(m, point, label=legends[i], alpha = al[i], color = col[i])
plt.plot()
plt.legend(loc = 'upper left', fontsize = 20)
plt.ylim([10**9, 10**12])
plt.xlim([0, 91])
plt.xlabel('Day(s) from start', fontsize = 20)
plt.ylabel('Loss ($)', fontsize = 20)
sns.despine()
plt.rcParams["figure.figsize"] = (10, 10)
plt.tick_params(labelsize=20)
plt.savefig('./content/sample_data/loss.pdf')