```
import pandas as pd # Loading the NYC dataset
import numpy as np # performing fast calculations to build our simulation
import plotly.express as px # Quick visuals for data exploration
from ipywidgets import widgets # sliders, buttons and visual layout of the dashboard
import plotly.graph_objects as go # plotly objects to place on our graph
import math # python mathematical functions
from scipy.integrate import odeint
# We'll set some easy defaults for this function so that we immediately have outputs which make some sense when we execute this method.
def calculate_values_SIR(removal_rate = .1
= 1.1
, infection_rate = 1000
, population = 150
, days = 1
, initial_infection
# we've added this additional parameter which allows us to cut the
# infection rate by 2 at a given day to simulate drastic policy action.
= 50
, intervention_day = .46
, intervention_effectiveness
):
# Define a series of ordinary differential equations (math-fancy way of saying equations that model the
# rate of change of something) notice that all of these values are similar to the ones we were using
# before, we're basically going to pass everything from our current method into this method.
def sir(previous_day
, day
, population
, infection_rate
, removal_rate
, intervention_day
, intervention_effectiveness):
# Initialise our values like we previously did. in this Ordinary Differential Equation Solver
# these values will correspond to the previous day values.
= previous_day[0] + previous_day[3]
susceptible = previous_day[1]
infected = previous_day[2]
recovered
= previous_day[1] + previous_day[3]
noi_infected = previous_day[4]
noi_removed
# Like we did before, check to see if we are past the intervention day, if we are, let's reduce
# the infection rate by the intervention effectiveness
if day >= intervention_day:
= infection_rate * max((1-intervention_effectiveness),0)
inter_infection_rate else:
= infection_rate
inter_infection_rate
# Calculate our daily changes for each of the categories:
# susceptible
= -min(inter_infection_rate * infected * susceptible / population, susceptible)
dsdt # Infected
= (inter_infection_rate * infected * susceptible / population) - (removal_rate * infected)
didt # Removed
= removal_rate * infected
drdt
# If we're after the intervention day, calculate the values as if we had never made an intervention
if day >= intervention_day:
= previous_day[0]
noi_susceptible
= (min(infection_rate * noi_infected * noi_susceptible / population, noi_susceptible))
noi_dsdt =(infection_rate * noi_infected * noi_susceptible / population) - (removal_rate * noi_infected)
noi_didt = removal_rate * noi_infected
noi_drdt
= noi_didt - didt
noi_didt = max(noi_drdt - drdt, 0 )
noi_drdt = dsdt - noi_didt - noi_drdt
dsdt
else:
= 0
noi_dsdt = 0
noi_didt = 0
noi_drdt
# Return our daily change values
return [dsdt, didt, drdt, noi_didt, noi_drdt]
# Set the values we calculate for
= np.rint(np.linspace(0,days,days))
x_days
# initialize our model with the strting conditions
=[population # susceptible start
starting_condition# infection start
, initial_infection 0 # removed start
, 0 # no intevention infected start
, 0] # no intervention removed start
,
# run a function which will simulate our equation system over a number of
# days with our starting positions and the the parameters we give it.
# this basically replaces the for loop from our previous implementation.
= odeint(sir
s
, starting_condition
, x_days=(population
, args
, infection_rate
, removal_rate
, intervention_day
, intervention_effectiveness))
# Round our calculated values to integers and then transpose our output
# from a list of rows to a list of columns.
= np.rint(s).T
dat
# Return the desired values
return {'infected': dat[1]
'removed' : dat[2]
, 'susceptible': dat[0]
, 'infected_noi': dat[3]
, 'removed_noi': dat[4] } ,
```

# Python 105 - Make the SIR model more efficient and add additional parameters

First let’s make the model more efficient. Fortunately, after exploring the functions and calculations that make up this model we can quite easily translate the calculations into a more mathematical system which computes much faster! To do this we will use a function provided by the scientific python library scipy which allows us to solve ordinary differential equations (ODEs). Don’t worry about the maths terminology, we’ve actually already done all the solveing, all we need to do is re-write our calculation method to use the more efficient library calculations.

## Display the new faster calculation method

Once we’ve updated the calculation method, we don’t need to make any changes to the visualization code from last week, we can take it as is and plug the new calculation method in to use it.

```
# This slider works with floating point numbers (hence being called Float Slider) and
# allows us to set a variable with this slider. This is going to be the way we set the
# infection rate.
= widgets.FloatSlider(
ir =1.187, # this is the initial value of our slider when it appears
valuemin=0.0, # the minimum value we'll allow
max=5.0, # the maximum value we'll allow
=.001, # by what increments the slider will change when we move it
step='Infection_rate:', # the name of the slider
description=True # Will this slider wait until it stops moving to
continuous_update# call it's update function or does it call the
# update function immediately?
)
= widgets.FloatSlider(
rr =.46,
valuemin=0.1, # this is set to greater than 0 because this is the denominator in the R0 calculation
max=2.0,
=.01,
step='Removal_Rate:',
description=True
continuous_update
)
= widgets.IntSlider(
ii =1,
valuemin=1,
max=50,
=1,
step='Initially Infected:',
description=True
continuous_update
)
= widgets.IntSlider(
ip =1000,
valuemin=500,
max=10_000_000,
=500,
step='Initial Population:',
description=True
continuous_update
)
= widgets.IntSlider(
iday =15,
valuemin=1,
max=150,
=1,
step='Day of intervention (reducing infection rate):',
description=True
continuous_update
)
= widgets.FloatSlider(
ie =.46,
valuemin=0.0,
max=1.0,
=.01,
step='Intervention effectiveness:',
description=True
continuous_update
)
= widgets.HBox(children=[ir, rr, ie])
first_slider_group = widgets.HBox(children=[ii, ip, iday])
second_slider_group
# First, we use the method created above to calculate a model using the initial
# values of the sliders we just created. Given that at this point we haven't
# displayed the sliders yet, their values will be the default values we set above.
= calculate_values_SIR( removal_rate = rr.value
data = ir.value
, infection_rate = ip.value
, population = 150
, days = ii.value
, initial_infection = iday.value
, intervention_day
)
# Next we add all the data traces to the chart
= go.Bar(x = list(range(1,len(data['infected'])))
infected_trace = data['infected']
,y ='Infected'
, name= dict(color='red')
, marker
)
= go.Bar(x = list(range(1, len(data['susceptible'])))
susceptible_trace = data['susceptible']
, y ='Susceptible'
, name= dict(color='rgba(0,0,255,0.5)')
, marker =0.5
, opacity
)
= go.Bar(x =list(range(1, len(data['removed'])))
removed_trace = data['removed']
,y ='Removed'
, name= dict(color='rgba(0,128,0,0.5)')
, marker =0.5)
, opacity
= go.Bar(x = list(range(1,len(data['infected_noi'])))
infected_trace_noi = data['infected_noi']
,y ='Infected No Intervention'
, name= dict(color='rgba(225,0,0,0.5)')
, marker
)
= go.Bar(x =list(range(1, len(data['removed_noi'])))
removed_trace_noi = data['removed_noi']
,y ='Removed No Intervention'
, name= dict(color='rgba(0,225,0,0.5)')
, marker
)
# This trace is interesting as it's basically just drawing a straight line on the
# selected intervention day.
= go.Scatter(x = [iday.value, iday.value]
intervention_day = [0, ip.value]
,y ='Intervention day'
, name= dict(color='darkblue')
, marker = dict(width=5)
, line
)
# We create our figure adding all the traces we created to the data list, and setting some layout values in the layout parameter.
= go.FigureWidget(data=[ infected_trace, infected_trace_noi, removed_trace, susceptible_trace, intervention_day, removed_trace_noi ],
g =go.Layout(
layout={
title'text': f'R0 = {ir.value / rr.value} <br /> Post-Intervention R0: {ir.value/ 2* rr.value} <br />Infection_rate={ir.value} Removal_rate={rr.value}',
'y':.95,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
='stack'
,barmode='x'
,hovermode=900
,height=dict(title='Number of Days')
,xaxis=dict(title='Number of People')
,yaxis
))
# This is to update the x-axis range to show only the days where we have cases
range=[0,np.where(data['infected']==0)[0][0]])
g.update_xaxes(
=1
x
# This method will be called any time one of the sliders is modified. It will re-run our model calculation
# with the new values and update the data for the 4 traces we added to the figure.
def response(change):
=150
num_days
# recalculate the model using the new values defined by the sliders
= calculate_values_SIR(removal_rate = rr.value
pop_values =ir.value
, infection_rate=ii.value
, initial_infection=ip.value
, population=num_days
, days= iday.value
, intervention_day = ie.value
, intervention_effectiveness
)
# Try to find the first day where we have no more infections,
# if that fails and we get an error, use the maximum number of
# days for which we've calculated the model. We use this later
# to update the x-axis range to keep our curve centered.
try:
= max( np.where(pop_values['infected_noi'][iday.value+1:]==0)[0][0] + iday.value +1,
end_infection 'infected']==0)[0][0] )
np.where(pop_values[except IndexError:
= num_days
end_infection
# plotly updates much faster if you send all the updates in one go,
# this construct allows us to open an update session on our chart
# and when we leave it, it will send all the updates to the figure
# at once.
with g.batch_update():
# update the y-axis values from the model
0].y = pop_values['infected']
g.data[1].y = pop_values['infected_noi']
g.data[2].y = pop_values['removed']
g.data[3].y = pop_values['susceptible']
g.data[5].y = pop_values['removed_noi']
g.data[
# Add the intervention day line
4].y = [0 , ip.value]
g.data[4].x = [iday.value, iday.value]
g.data[
# update the title to show the R values as well as the infection rate and the removal rate
= ir.value * max((1-ie.value),0)
ie_rate = {
g.layout.title 'text': f'R0 = {ir.value / rr.value} <br /> Post-Intervention R0: {ie_rate / rr.value} <br /> Infection_rate={ir.value} Removal_rate={rr.value}',
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
# change the x-axis range to keep the infection curve in scope
= dict(range=[0,end_infection])
g.layout.xaxis
# Update each of the widgets and register our update method as the method to
# call when they change.
='value')
ir.observe(response, names='value')
rr.observe(response, names='value')
ii.observe(response, names='value')
ip.observe(response, names='value')
iday.observe(response, names='value')
ie.observe(response, names
# put the widgets, and our chart together into our layout
widgets.VBox([first_slider_group, second_slider_group, g])
```

# Creating the SEIRDS model

In order to prepare to visualize our SEIRDS model we will need to calculate the number of people in the susceptible, exposed, infectious, deceased, and recovered categories for every single day. Given that we were able to create the previous SIR modle with a series of lists, we can do the same here. In fact, the method we just implemented above makes that very easy for us.

```
days ---> [ 1, 2, 3, 4, 5, 6, 7, 8, 9, ..... 196 ]
susceptible ---> [1000, 999, 990, 949, 906, 810, 690, 550, 400, ..... 0 ]
exposed. ---> [ 0, 1, 10, 35, 80, 140, 190, 250, 320, ..... 0 ]
infectious ---> [ 0, 0, 0, 1, 4, 40, 90, 150, 200, ..... 0 ]
recovered ---> [ 0, 0, 0, 5, 10, 20, 30, 50, 80, ..... 900 ]
deceased ---> [ 0, 0, 0, 0, 0, 1, 2, 3, 4, ..... 100 ]
```

All we need to do is set some initial conditions, and identify how to calculate the difference from one day to the next. If we pass that information to the Ordinary Differential Equation solver then it will create the lists for us. Well, it will create a list of daily figures with the numbers for each day (so the transpose of the representation above). We can use that easily by transposing it.

```
import pandas as pd # converting the calculated data into a dataframe for easy display
import numpy as np # performing fast calculations to build our simulation
import plotly.express as px # Quick visuals for data exploration
from ipywidgets import widgets # sliders, buttons and visual layout of the dashboard
import plotly.graph_objects as go # plotly objects to place on our graph
import math # python mathematical functions
from scipy.integrate import odeint
# We'll set some easy defaults for this function so that we immediately have outputs which make some sense when we execute this method.
def calculate_values_SEIRDS(days_to_recover = 16
= 1.1
, infection_rate = 1000
, population = 150
, days = 1
, initial_infection
# These parameters will be unused for now, but when we re-add this functionality later
# we can use them to set up our intervention again
= 50
, intervention_day = .46
, intervention_effectiveness
# Adding in the incubation period to determine the time it takes for exposure
= 14
, incubation_period
# Adding the fatality rate
= .01
, fatality_rate
# Adding the immunity period
= 180
, relapse_period
):
# Define a series of ordinary differential equations (math-fancy way of saying equations that model the
# rate of change of something) notice that all of these values are similar to the ones we were using
# before, we're basically going to pass everything from our current method into this method.
def sir(previous_day
, day
, population
, infection_rate
, days_to_recover
, intervention_day
, intervention_effectiveness
, incubation_period
, fatality_rate
, relapse_period):
= previous_day[0]
susceptible = previous_day[1]
exposed = previous_day[2]
infected = previous_day[3]
recovered = previous_day[4]
fatalities
= 1/incubation_period
incubation = 1/relapse_period
relapse = 1/days_to_recover
recovery_rate
= infection_rate
inter_infection_rate
# Daily change in susceptible
= -min(inter_infection_rate * infected * susceptible / population, susceptible) + (relapse * recovered)
dsdt
# Daily change in exposed
= (inter_infection_rate * infected * susceptible / population) - (incubation * exposed)
dedt
# Daily change in infected
= incubation * exposed - (recovery_rate * infected) - (fatality_rate * infected)
didt
# Daily change in recovered
= recovery_rate * infected - (relapse * recovered)
drdt
# Daily change in fatalities
= fatality_rate * infected
dfdt
return [dsdt, dedt, didt, drdt, dfdt]
= np.rint(np.linspace(0,days,days))
x_days
=[population # susceptible start
starting_condition# exposure start
, initial_infection 0 # infection start
, 0 # recovered
, 0] # fatalities
,
# run a function which will simulate our equation system over a number of
# days with our starting positions and the the parameters we give it.
# this basically replaces the for loop from our previous implementation.
= odeint(sir
s
, starting_condition
, x_days=(population
, args
, infection_rate
, days_to_recover
, intervention_day
, intervention_effectiveness
, incubation_period
, fatality_rate
, relapse_period))
# Round our calculated values to integers and then transpose our output
# from a list of rows to a list of columns.
= np.rint(s).T
dat
# Return the desired values
return { 'infected' : dat[2]
'exposed' : dat[1]
, 'recovered' : dat[3]
, 'susceptible': dat[0]
, 'fatalities' : dat[4] } ,
```

That was a fair amount of calculation, but really no different from what we were doing last week and the week before. So before we jump into updating our dashboard to handle these new values, let’s have a quick look at what they look like to check we’re not way off the mark.

```
= calculate_values_SEIRDS( days_to_recover = 16
df = 1.34
, infection_rate = 7_000_000
, population = 600
, days = 1
, initial_infection # , intervention_day = iday.value
= 0.03
, fatality_rate = 180
, relapse_period
)
# create a dataframe out of the data we calculated. Note that we have to call reset_index to
# get a column with the days for the visualization.
= pd.DataFrame(df)
pdf pdf
```

infected | exposed | recovered | susceptible | fatalities | |
---|---|---|---|---|---|

0 | 0.0 | 1.0 | 0.0 | 7000000.0 | 0.0 |

1 | 0.0 | 1.0 | 0.0 | 7000000.0 | 0.0 |

2 | 0.0 | 1.0 | 0.0 | 7000000.0 | 0.0 |

3 | 0.0 | 1.0 | 0.0 | 7000000.0 | 0.0 |

4 | 0.0 | 1.0 | 0.0 | 6999999.0 | 0.0 |

... | ... | ... | ... | ... | ... |

595 | 109454.0 | 139284.0 | 1750564.0 | 464178.0 | 4536521.0 |

596 | 109278.0 | 139061.0 | 1747682.0 | 464177.0 | 4539802.0 |

597 | 109103.0 | 138838.0 | 1744805.0 | 464176.0 | 4543078.0 |

598 | 108928.0 | 138616.0 | 1741933.0 | 464175.0 | 4546348.0 |

599 | 108754.0 | 138394.0 | 1739066.0 | 464173.0 | 4549614.0 |

600 rows × 5 columns

```
= pdf.reset_index()
pdf
# We melt the dataframe (this is basically an unpivot) to make sure that we have one x value and
# one y value column, we can differentiate between the different groups with the variable column
= pdf.melt(id_vars=['index'])
pdf ='index', y='value', color='variable') px.line(pdf, x
```

The data looks good now we can put all of these plots into our dashboard.

```
# This slider works with floating point numbers (hence being called Float Slider) and
# allows us to set a variable with this slider. This is going to be the way we set the
# infection rate.
= widgets.FloatSlider(
ir =1.187, # this is the initial value of our slider when it appears
valuemin=0.0, # the minimum value we'll allow
max=5.0, # the maximum value we'll allow
=.001, # by what increments the slider will change when we move it
step='Infection_rate:', # the name of the slider
description=True # Will this slider wait until it stops moving to
continuous_update# call it's update function or does it call the
# update function immediately?
)
= widgets.IntSlider(
rr =16,
valuemin=1, # this is set to greater than 0 because this is the denominator in the R0 calculation
max=30,
=1,
step='Days to Recover:',
description=True
continuous_update
)
= widgets.IntSlider(
ii =1,
valuemin=1,
max=50,
=1,
step='Initially Infected:',
description=True
continuous_update
)
= widgets.IntSlider(
ip =7_000_000,
valuemin=500,
max=10_000_000,
=500,
step='Initial Population:',
description=True
continuous_update
)
= widgets.IntSlider(
iday =15,
valuemin=1,
max=150,
=1,
step='Day of intervention (reducing infection rate):',
description=True
continuous_update
)
= widgets.FloatSlider(
ie =.46,
valuemin=0.0,
max=1.0,
=.01,
step='Intervention effectiveness:',
description=True
continuous_update
)
= widgets.FloatSlider(
fr =.01,
valuemin=0.0,
max=.5,
=.001,
step='Fatality Rate:',
description=True
continuous_update
)
= widgets.IntSlider(
relr =180,
valuemin=0.0,
max=600,
=1,
step='Immunity Period:',
description=True
continuous_update
)
= widgets.IntSlider(
incu =14,
valuemin=0.0,
max=30,
=1,
step='Incubation Period:',
description=True
continuous_update
)
= widgets.HBox(children=[ir, fr, ip])
first_slider_group = widgets.HBox(children=[rr, relr, incu])
second_slider_group = widgets.HBox(children=[ie, iday])
third_slider_group
# First, we use the method created above to calculate a model using the initial
# values of the sliders we just created. Given that at this point we haven't
# displayed the sliders yet, their values will be the default values we set above.
= calculate_values_SEIRDS( days_to_recover = rr.value
data = ir.value
, infection_rate = ip.value
, population = 200
, days = ii.value
, initial_infection = iday.value
, intervention_day = fr.value
, fatality_rate=relr.value
, relapse_period=incu.value
, incubation_period
)
# Next we add all the data traces to the chart
= go.Bar(x = list(range(1,len(data['infected'])))
infected_trace = data['infected']
,y ='Infected'
, name= dict(color='red')
, marker
)
= go.Bar(x = list(range(1, len(data['susceptible'])))
susceptible_trace = data['susceptible']
, y ='Susceptible'
, name= dict(color='rgba(0,0,255,0.5)')
, marker =0.5
, opacity
)
= go.Bar(x =list(range(1, len(data['recovered'])))
removed_trace = data['recovered']
,y ='Recovered'
, name= dict(color='rgba(0,128,0,0.5)')
, marker =0.5)
, opacity
######
# ADDITION STARTS HERE
######
= go.Bar(x = list(range(1,len(data['exposed'])))
exposed_trace = data['exposed']
,y ='Exposed'
, name= dict(color='rgba(125,125,0,0.5)')
, marker
)
= go.Bar(x =list(range(1, len(data['fatalities'])))
fatality_trace = data['fatalities']
,y ='Fatalities'
, name= dict(color='rgba(115,115,115,0.5)')
, marker
)
######
# ADDITION END HERE
######
# This trace is interesting as it's basically just drawing a straight line on the
# selected intervention day.
# intervention_day = go.Scatter(x = [iday.value, iday.value]
# ,y = [0, ip.value]
# , name='Intervention day'
# , marker = dict(color='darkblue')
# , line = dict(width=5)
# )
= ir.value * max((1-ie.value),0)
ie_rate = fr.value + (1/np.sqrt(rr.value)) + (1/incu.value)
rr_rate # We create our figure adding all the traces we created to the data list, and setting some layout values in the layout parameter.
= go.FigureWidget(data=[fatality_trace, infected_trace, exposed_trace, removed_trace, susceptible_trace],
g =go.Layout(
layout={
title'text': f' Infection_rate={ir.value} Days to recover={rr.value}',
'y':.95,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
='stack'
,barmode='x'
,hovermode=900
,height=dict(title='Number of Days')
,xaxis=dict(title='Number of People')
,yaxis
))
# This is to update the x-axis range to show only the days where we have cases
try:
= [0,np.where(data['exposed'][:1]==0)[0][0]]
rng except:
= [0,len(data['exposed'])]
rng
range=rng)
g.update_xaxes(
=1
x
# This method will be called any time one of the sliders is modified. It will re-run our model calculation
# with the new values and update the data for the 4 traces we added to the figure.
def response(change):
=200
num_days
# recalculate the model using the new values defined by the sliders
= calculate_values_SEIRDS(days_to_recover = rr.value
pop_values =ir.value
, infection_rate=ii.value
, initial_infection=ip.value
, population=num_days
, days= iday.value
, intervention_day = ie.value
, intervention_effectiveness= fr.value
, fatality_rate=relr.value
, relapse_period=incu.value
, incubation_period
)
# Try to find the first day where we have no more infections,
# if that fails and we get an error, use the maximum number of
# days for which we've calculated the model. We use this later
# to update the x-axis range to keep our curve centered.
try:
= np.where(pop_values['exposed'][1:]==0)[0][0]
end_infection except IndexError:
= num_days
end_infection
# plotly updates much faster if you send all the updates in one go,
# this construct allows us to open an update session on our chart
# and when we leave it, it will send all the updates to the figure
# at once.
with g.batch_update():
# update the y-axis values from the model, we changed the order above so we need to adapt it here
# [fatality_trace, infected_trace, exposed, removed_trace, susceptible_trace, intervention_day ]
# ------ 0 ------- --- 1 -------- -- 2 ---- ------- 3 ---- -------- 4 ------- ------- 5 -------
0].y = pop_values['fatalities']
g.data[1].y = pop_values['infected']
g.data[2].y = pop_values['exposed']
g.data[3].y = pop_values['recovered']
g.data[4].y = pop_values['susceptible']
g.data[
# Add the intervention day line
# g.data[5].y = [0 , ip.value]
# g.data[5].x = [iday.value, iday.value]
# update the title to show the R values as well as the infection rate and the removal rate
= ir.value * max((1-ie.value),0)
ie_rate = fr.value + (1/np.sqrt(rr.value)) + (1/incu.value)
rr_rate = {
g.layout.title 'text': f' Infection_rate={ir.value} Days to recover={rr.value}',
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
# change the x-axis range to keep the infection curve in scope
= dict(range=[0,end_infection])
g.layout.xaxis
# Update each of the widgets and register our update method as the method to
# call when they change.
='value')
ir.observe(response, names='value')
rr.observe(response, names='value')
ii.observe(response, names='value')
ip.observe(response, names='value')
iday.observe(response, names='value')
ie.observe(response, names
='value')
fr.observe(response, names='value')
relr.observe(response, names='value')
incu.observe(response, names
# put the widgets, and our chart together into our layout
widgets.VBox([first_slider_group, second_slider_group, g])
```

# Challenge (Hard): Add notion of intervention day to the SEIRDS model

This is a hard challenge, can you add the concept of an intervention day to the SEIRDS model? some of the peices are already there, you’ll need to modify the model to take into account the intervention day in the same way that we did for the SIR model.

```
import pandas as pd # Loading the NYC dataset
import numpy as np # performing fast calculations to build our simulation
import plotly.express as px # Quick visuals for data exploration
from ipywidgets import widgets # sliders, buttons and visual layout of the dashboard
import plotly.graph_objects as go # plotly objects to place on our graph
import math # python mathematical functions
from scipy.integrate import odeint
# We'll set some easy defaults for this function so that we immediately have outputs which make some sense when we execute this method.
def calculate_values_SEIRS(days_to_recover = 16
= 1.1
, infection_rate = 1000
, population = 150
, days = 1
, initial_infection
# These parameters will be unused for now, but when we re-add this functionality later
# we can use them to set up our intervention again
= 50
, intervention_day = .46
, intervention_effectiveness
# Adding in the incubation period to determine the time it takes for exposure
= 14
, incubation_period
# Adding the fatality rate
= .01
, fatality_rate
# Adding the immunity period
= 180
, relapse_period
):
# Define a series of ordinary differential equations (math-fancy way of saying equations that model the
# rate of change of something) notice that all of these values are similar to the ones we were using
# before, we're basically going to pass everything from our current method into this method.
def sir(previous_day
, day
, population
, infection_rate
, days_to_recover
, intervention_day
, intervention_effectiveness
, incubation_period
, fatality_rate
, relapse_period):
= previous_day[0]
susceptible = previous_day[1]
exposed = previous_day[2]
infected = previous_day[3]
recovered = previous_day[4]
fatalities
= 1/incubation_period
incubation = 1/relapse_period
relapse = 1/days_to_recover
recovery_rate
if day >= intervention_day:
= infection_rate * max((1-intervention_effectiveness),0)
inter_infection_rate else:
= infection_rate
inter_infection_rate
# Daily change in susceptible
= -min(inter_infection_rate * infected * susceptible / population, susceptible) + (relapse * recovered)
dsdt
# Daily change in exposed
= (inter_infection_rate * infected * susceptible / population) - (incubation * exposed)
dedt
# Daily change in infected
= incubation * exposed - (recovery_rate * infected) - (fatality_rate * infected)
didt
# Daily change in recovered
= recovery_rate * infected - (relapse * recovered)
drdt
# Daily change in fatalities
= fatality_rate * infected
dfdt
# if day >= intervention_day:
# noi_susceptible = previous_day[0]
# noi_dsdt = (min(infection_rate * noi_infected * noi_susceptible / population, noi_susceptible))
# noi_didt =(infection_rate * noi_infected * noi_susceptible / population) - (removal_rate * noi_infected)
# noi_drdt = removal_rate * noi_infected
# noi_didt = noi_didt - didt
# noi_drdt = max(noi_drdt - drdt, 0 )
# dsdt = dsdt - noi_didt - noi_drdt
# else:
# noi_dsdt = 0
# noi_didt = 0
# noi_drdt = 0
return [dsdt, dedt, didt, drdt, dfdt]
= np.rint(np.linspace(0,days,days))
x_days
=[population # susceptible start
starting_condition# exposure start
, initial_infection 0 # infection start
, 0 # recovered
, 0] # fatalities
,
# run a function which will simulate our equation system over a number of
# days with our starting positions and the the parameters we give it.
# this basically replaces the for loop from our previous implementation.
= odeint(sir
s
, starting_condition
, x_days=(population
, args
, infection_rate
, days_to_recover
, intervention_day
, intervention_effectiveness
, incubation_period
, fatality_rate
, relapse_period))
# Round our calculated values to integers and then transpose our output
# from a list of rows to a list of columns.
= np.rint(s).T
dat
# Return the desired values
return { 'infected' : dat[2]
'exposed' : dat[1]
, 'recovered' : dat[3]
, 'susceptible': dat[0]
, 'fatalities' : dat[4] }
,
= calculate_values_SEIRS( days_to_recover = 16
df = ir.value
, infection_rate = 7_000_000
, population = 600
, days = ii.value
, initial_infection = iday.value
, intervention_day = 0.03
, fatality_rate = 180
, relapse_period
)
= pd.DataFrame(df).reset_index()
pdf = pdf.melt(id_vars=['index'])
pdf ='index', y='value', color='variable') px.line(pdf, x
```

# Challenge (Very Hard): Draw the situation without intervention on the graph

Can you add in the same hypothetical scenarios as we had in the SIR model to the graph?