%%html
<script>
// AUTORUN ALL CELLS ON NOTEBOOK-LOAD!
require(
['base/js/namespace', 'jquery'],
function(jupyter, $) {
$(jupyter.events).on("kernel_ready.Kernel", function () {
console.log("Auto-running all cells-below...");
jupyter.actions.call('jupyter-notebook:run-all-cells-below');
jupyter.actions.call('jupyter-notebook:save-notebook');
});
}
);
</script>
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
from tqdm import tqdm
import ipywidgets as widgets
from ipywidgets import interact, interactive, fixed, interact_manual, interactive_output, Layout
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import optimize
from IPython.display import display
from matplotlib.gridspec import GridSpec
import matplotlib
matplotlib.rcParams.update({'font.size': 22})
%matplotlib inline
# def get_rate_of_change(SIRvals, t, InitVals):
def get_rate_of_change(y, t, beta, gamma):
S, I, R = y
dydt = [-beta*S*I, beta*S*I - gamma*I, gamma*I]
return dydt
def plot_func(beta,gamma, N, I0, tmax, data, logscale_y=False):
fig=plt.figure(figsize=[12,6])
gs=GridSpec(3,2) # 2 rows, 3 columns
ax1=fig.add_subplot(gs[0,1]) # First row, first column
ax2=fig.add_subplot(gs[1,1]) # First row, second column
ax3=fig.add_subplot(gs[2,1]) # First row, third column
ax4=fig.add_subplot(gs[:,0]) # Second row, span all columns
for ax in [ax1,ax2,ax3]:
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
ax4.set_xlabel('t')
ax4.set_ylabel("Number of Cases")
ax4.grid()
ax1.set_ylabel("Surceptible")
ax2.set_ylabel("Infected")
ax3.set_ylabel("Recovered")
ax4.text(0.5,0.9,"$R_0={0:.1f}$".format(beta/gamma),transform=ax4.transAxes)
fig.subplots_adjust(wspace=0, hspace=0.06)
t = np.linspace(0, tmax, tmax*100)
S0 = N-I0
R0 = 0
D0 = 0
y0 = [S0/N,I0/N,R0/N]
sol = odeint(get_rate_of_change, y0, t, args=(beta,gamma))
ax1.plot(t, sol[:, 0]*N, 'b', label='Surceptible({0:.2f}%)'.format(sol[:, 0][-1]*100))
ax2.plot(t, sol[:, 1]*N, 'r', label='Infected({0:.2f}%)'.format(sol[:, 1][-1]*100))
ax3.plot(t, sol[:, 2]*N, 'g', label='Recovered ({0:.2f}%)'.format(sol[:, 2][-1]*100))
ax4.plot(t, sol[:, 0]*N, 'b', label='Surceptible({0:.2f}%)'.format(sol[:, 0][-1]*100))
ax4.plot(t, sol[:, 1]*N, 'r', label='Infected({0:.2f}%)'.format(sol[:, 1][-1]*100))
ax4.plot(t, sol[:, 2]*N, 'g', label='Recovered ({0:.2f}%)'.format(sol[:, 2][-1]*100))
for ax in [ax1,ax2]:
ax.legend(loc='upper right', frameon=True)
ax3.legend(loc='lower right', frameon=True)
try:
data = pd.read_csv("Data/"+data+".txt")
data["Infected"] = data.TotalInfected.diff()
ydata = data.Infected.values[1:]
xdata = data.Day.values[1:] - data.Day.values[1]
ax2.plot(xdata, ydata, 'o')
except:
pass
if logscale_y:
ax4.set_ylim([1,None])
ax4.set_yscale("log")
dashboard_SIR = interactive(plot_func,
N = widgets.IntText(value=130000, min=10, continuous_update=True,
description ="Population Size",
style={'description_width': 'initial'}),
I0 = widgets.IntText(value=18, min=0, continuous_update=True,
description = "Initial Infected Population",
style={'description_width': 'initial'}),
beta = widgets.FloatSlider(value=1.75,
min=0,
max=2.0,
step=0.01,continuous_update=False, description="Transmission Rate",
style={'description_width': 'initial'}),
gamma = widgets.FloatSlider(value=1.55,
min=0,
max=2.0,
step=0.01,continuous_update=False, description="Recovery Rate",
style={'description_width': 'initial'}),
tmax = widgets.IntSlider(value=60, min = 1 , max= 150, continuous_update=False, description = "Time"),
data = widgets.Dropdown(options=["None","SouthKorea"], value="SouthKorea",description="Select Country",
style={'description_width': 'initial'})
)
init_state = [i.value for i in dashboard_SIR.children[:-1]].copy()
def reset_values(a):
for i, item in enumerate(dashboard_SIR.children[:-1]):
item.value = init_state[i]
reset_button = widgets.Button(description= "Reset")
reset_button.on_click(reset_values)
display(reset_button)
display(dashboard_SIR)
Button(description='Reset', style=ButtonStyle())
interactive(children=(FloatSlider(value=1.75, continuous_update=False, description='Transmission Rate', max=2.…
# def get_rate_of_change(SIRvals, t, InitVals):
def get_rate_of_change_with_deaths(y, t, beta, gamma, omega):
X, Y, Z, D = y
N = X + Y + Z
m = omega/float(1-omega)*gamma
dydt = [-beta*X*Y/N, beta*X*Y/N - gamma*Y - m*Y, gamma*Y, m*Y]
return dydt
def plot_func(tmax,beta,gamma,omega, N, I0, logscale_x,logscale_y):
N = float(N)
X0 = N-I0
Y0 = I0
Z0 = 0
D0 = 0
y0 = [X0,Y0,Z0,D0]
t = np.linspace(0, tmax, tmax*1000)
omega = omega/100.0
m = omega/float(1-omega)*gamma
sol = odeint(get_rate_of_change_with_deaths, y0, t, args=(beta,gamma,omega))
fig=plt.figure(figsize=[12,6])
gs=GridSpec(4,2) # 2 rows, 3 columns
ax1=fig.add_subplot(gs[0,1]) # First row, first column
ax2=fig.add_subplot(gs[1,1]) # First row, second column
ax3=fig.add_subplot(gs[2,1]) # First row, third column
ax4=fig.add_subplot(gs[3,1])
ax5=fig.add_subplot(gs[:,0]) # Second row, span all columns
for ax in [ax1,ax2,ax3,ax4]:
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
ax5.set_xlabel('t')
ax5.set_ylabel("Number of Cases")
ax5.grid()
ax1.set_ylabel("Surceptible")
ax2.set_ylabel("Infected")
ax3.set_ylabel("Recovered")
ax4.set_ylabel("Dead")
ax4.set_xlabel("t")
ax5.text(0.5,0.9,"$R_0={0:.1f}$".format(beta/gamma),transform=ax5.transAxes)
fig.subplots_adjust(wspace=0, hspace=0.06)
ax5.plot(t, sol[:, 0], 'b', label='Surceptible({})'.format(int(sol[:, 0][-1])))
ax5.plot(t, sol[:, 1], 'r', label='Infected({})'.format(int(sol[:, 1][-1])))
ax5.plot(t, sol[:, 2], 'g', label='Recovered ({})'.format(int(sol[:, 2][-1])))
ax5.plot(t, sol[:, 3], 'k', \
label='Dead ({0:.2f}%'.format(sol[:, 3][-1]/N*100)+", {} People".format(int(sol[:, 3][-1])))
ax1.plot(t, sol[:, 0], 'b', label='Surceptible({})'.format(int(sol[:, 0][-1])))
ax2.plot(t, sol[:, 1], 'r', label='Infected({})'.format(int(sol[:, 1][-1])))
ax3.plot(t, sol[:, 2], 'g', label='Recovered ({})'.format(int(sol[:, 2][-1])))
ax4.plot(t, sol[:, 3], 'k', \
label='Dead ({0:.2f}%'.format(sol[:, 3][-1]/N*100)+", {} People".format(int(sol[:, 3][-1])))
for ax in [ax1,ax2,ax3,ax4]:
ax.legend(loc="center right", frameon=True)
ax5.set_ylim([0,None])
if logscale_x:
ax5.set_xscale("log")
if logscale_y:
ax5.set_yscale("log")
ax5.set_ylim([1,None])
ax5.legend(loc='center right', ncol = 1, frameon=True)
ContactNumber = widgets.FloatSlider(value=10, min=0, max=15,step=0.1, continuous_update=False,
description="Contact Number",
style={'description_width': 'auto'})
InfectionProbability = widgets.FloatSlider(value=0.5, min=0, max=1, step=0.01,continuous_update=False,
description="Infection Probabity",
style={'description_width': 'initial'})
def on_value_chage(change):
beta.value = -ContactNumber.value*np.log(1-InfectionProbability.value)
ContactNumber.observe(on_value_chage, names='value')
InfectionProbability.observe(on_value_chage, names='value')
beta = widgets.FloatSlider(value=7,
min=0,
max=10.0,
step=0.01,continuous_update=False, description="Transmission Rate",
style={'description_width': 'initial'})
gamma = widgets.FloatSlider(value=1.3,
min=0,
max=2.0,
step=0.01,continuous_update=False, description="Recovery Rate",
style={'description_width': 'initial'})
omega = widgets.FloatSlider(value=5,
min=0,
max=99.99,
step = 0.1,
continuous_update=False, description="Death Rate (%)",
style={'description_width': 'initial'})
tmax = widgets.IntSlider(value=8, min = 1, max=150, step = 1,
description="Max Time",continuous_update=False,
style={'description_width': 'initial'})
N = widgets.IntText(value=1000, min=10, continuous_update=True,
description ="Population Size",
style={'description_width': 'initial'})
I0 = widgets.IntText(value=1, min=0, continuous_update=True,
description = "Initial Infected Population",
style={'description_width': 'initial'})
logscale_x = widgets.Checkbox(value=False,
description = "Logscale X axis",
style={'description_width': 'initial'})
logscale_y = widgets.Checkbox(value=False,
description = "Logscale Y axis",
style={'description_width': 'initial'})
dashboard_elements = {
'beta' : beta,
'gamma' : gamma,
'omega' :omega,
'tmax' : tmax,
'N' : N,
'I0' : I0,
'logscale_x' : logscale_x,
'logscale_y' : logscale_y
}
dashboard = interactive_output(plot_func, dashboard_elements)
init_values = {key:dashboard_elements[key].value for key in dashboard_elements.keys()}
def reset_values(a):
for key in init_values.keys():
dashboard_elements[key].value = init_values[key]
reset_button = widgets.Button(description= "Reset")
reset_button.on_click(reset_values)
display(widgets.HBox([
widgets.VBox([ContactNumber,InfectionProbability,beta, gamma,omega,tmax]),
widgets.VBox([N,I0,logscale_x,logscale_y])
]), reset_button, dashboard)
HBox(children=(VBox(children=(FloatSlider(value=10.0, continuous_update=False, description='Contact Number', m…
Button(description='Reset', style=ButtonStyle())
Output()
# def calc_instantaneous_mortality_risk(infected_duration, mortality_risk):
# if infected_duration < 6:
# risk = 0
# elif (infected_duration >=6) and (infected_duration<12):
# risk = (mortality_risk/6.0)*(infected_duration-6)/6.0
# elif (infected_duration >=12) and (infected_duration<18):
# risk = (mortality_risk/6.0)*(18-infected_duration)/6.0
# else:
# risk = 0
# return risk
def calc_instantaneous_mortality_risk(infected_duration, mortality_risk,recovery_time):
n = recovery_time/3.0
if infected_duration < n:
risk = 0
elif (infected_duration >=n) and (infected_duration<2*n):
risk = (mortality_risk/n)*(infected_duration-n)/n
elif (infected_duration >=2*n) and (infected_duration<3*n):
risk = (mortality_risk/n)*(3*n-infected_duration)/n
else:
risk = 0
return risk
def prob_toss(p):
return np.random.binomial(1, p)
def generate_mean_with_probability(n, var=2):
return max(0,round(np.random.normal(n,var)))
def run(N,I0,tmax,average_contact_number,transmission_probability,recovery_time,mortality_risk):
recovery_rate = 1/float(recovery_time)
class citizen:
def __init__(self, infected=False, age=30, infected_duration=0):
self.infected = infected
self.recovered = False
self.age = age
self.infected_duration = infected_duration
self.mortality_risk = mortality_risk
self.instantaneous_risk = self.mortality_risk/12
self.dead = False
#Generate population
population = np.array([citizen() for i in range(N)])
#Infect some people
for i in range(I0):
population[np.random.randint(N)].infected = True
fig=plt.figure(figsize=[18,9])
gs=GridSpec(3,2) # 2 rows, 3 columns
ax1=fig.add_subplot(gs[0,1]) # First row, first column
ax2=fig.add_subplot(gs[1,1]) # First row, second column
ax3=fig.add_subplot(gs[2,1]) # First row, third column
ax4=fig.add_subplot(gs[:,0]) # Second row, span all columns
for ax in [ax1,ax2,ax3]:
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
fig.subplots_adjust(wspace=0, hspace=0.06)
D_vals, I_vals, R_vals = [] ,[],[]
tlist = range(tmax)
for t in tqdm(tlist):
for p in population:
if ((not p.infected) or (p.dead)):
pass
else:
contact_number = generate_mean_with_probability(average_contact_number)
contacts = population[[np.random.randint(N) for i in range(contact_number)]]
for contact in contacts:
if (contact.infected==False) and (contact.recovered==False) and (contact.dead==False):
if prob_toss(transmission_probability)==1:
contact.infected = True
p.instantaneous_risk = calc_instantaneous_mortality_risk(p.infected_duration, p.mortality_risk, recovery_time)
if prob_toss(p.instantaneous_risk)==1:
p.dead = True
if not p.dead:
p.infected_duration += 1
if prob_toss(recovery_rate)==1:
p.recovered = True
p.infected = False
I_t = len([i.infected for i in population if i.infected])
R_t = len([i.recovered for i in population if i.recovered])
D_t = len([i.dead for i in population if i.dead])
I_vals.append(I_t)
R_vals.append(R_t)
D_vals.append(D_t)
ax1.plot(t,I_t,'yo')
ax2.plot(t,R_t,'go')
ax3.plot(t,D_t,'ro')
ax4.plot(t,I_t,'yo')
ax4.plot(t,D_t,'ro')
ax4.plot(t,R_t,'go')
ax3.set_xlabel("t")
ax1.set_ylabel("Infected")
ax2.set_ylabel("Recovered")
ax3.set_ylabel("Dead")
ax1.plot(tlist,I_vals,'y-')
ax2.plot(tlist,R_vals,'g-')
ax3.plot(tlist,D_vals,'r-')
ax4.plot(tlist,I_vals,'y-')
ax4.plot(tlist,D_vals,'r-')
ax4.plot(tlist,R_vals,'g-')
N = widgets.IntText(value=1000, description = "Population",
style={'description_width': 'initial'})
I0 = widgets.IntSlider(value=1,min=0,max=100,step=1,description = "Infected",
style={'description_width': 'initial'})
tmax = widgets.IntSlider(value=40,min=0,max=200,step=1,description = "Time",
style={'description_width': 'initial'})
average_contact_number = widgets.FloatSlider(value=2.5, min = 0, max=10,step=0.1,
description = "Averge Number of Contacts",
style={'description_width': 'initial'})
transmission_probability = widgets.FloatSlider(value=0.5, min = 0, max=1,step=0.01,
description = "Transmission Probability",
style={'description_width': 'initial'})
mortality_risk = widgets.FloatSlider(value=0.5, min = 0, max=1,step=0.01,
description = "Mortality Risk",
style={'description_width': 'initial'})
recovery_time = widgets.IntSlider(value=18, min = 1, max=30,step=1, description = "Recovery Time",
style={'description_width': 'initial'})
dashboard = interact_manual(run, N=N, I0=I0,tmax=tmax,
average_contact_number=average_contact_number,
transmission_probability=transmission_probability,
recovery_time=recovery_time, mortality_risk=mortality_risk)
interactive(children=(IntText(value=1000, description='Population', style=DescriptionStyle(description_width='…