<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Chapter 11
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  In this chapter, we develop a model of an epidemic as it spreads in a
susceptible population, and use it to evaluate the effectiveness of
possible interventions.
My presentation of the model in the next few chapters is based on an excellent article by David Smith and Lang Moore, [^1]
Step2:  And then convert the numbers to fractions by dividing by the total
Step3:  For now, let's assume we know the time between contacts and time between
recoveries
Step4:  We can use them to compute the parameters of the model
Step5:  Now we need a System object to store the parameters and initial
conditions. The following function takes the system parameters and returns a new System object
Step6:  The default value for t_end is 14 weeks, about the length of a
semester.
Here's what the System object looks like.
Step7:  The update function
At any point in time, the state of the system is represented by a
State object with three variables, S, I and R. So I'll define an
update function that takes as parameters a State object, the current
time, and a System object
Step8:  The first line uses a feature we have not seen before, multiple
assignment. The value on the right side is a State object that
contains three values. The left side is a sequence of three variable
names. The assignment does just what we want
Step9:  You might notice that this version of update_func does not use one of its parameters, t. I include it anyway because update functions
sometimes depend on time, and it is convenient if they all take the same parameters, whether they need them or not.
Running the simulation
Now we can simulate the model over a sequence of time steps
Step10:  The parameters of run_simulation are the System object and the
update function. The System object contains the parameters, initial
conditions, and values of t0 and t_end.
We can call run_simulation like this
Step11:  The result is the final state of the system
Step12:  This result indicates that after 14 weeks (98 days), about 52% of the
population is still susceptible, which means they were never infected,
less than 1% are actively infected, and 48% have recovered, which means they were infected at some point.
Collecting the results
The previous version of run_simulation only returns the final state,
but we might want to see how the state changes over time. We'll consider two ways to do that
Step13:  First, we create TimeSeries objects to store the results. Notice that
the variables S, I, and R are TimeSeries objects now.
Next we initialize state, t0, and the first elements of S, I and
R.
Inside the loop, we use update_func to compute the state of the system
at the next time step, then use multiple assignment to unpack the
elements of state, assigning each to the corresponding TimeSeries.
At the end of the function, we return the values S, I, and R. This
is the first example we have seen where a function returns more than one
value.
Now we can run the function like this
Step14:  We'll use the following function to plot the results
Step15:  And run it like this
Step16:  Notice that it takes about three weeks (21 days) for the outbreak to get going, and about six weeks (42 days) before it peaks. The fraction of the population that's infected is never very high, but it adds up. In total, almost half the population gets sick.
Now with a TimeFrame
If the number of state variables is small, storing them as separate
TimeSeries objects might not be so bad. But a better alternative is to use a TimeFrame, which is another object defined in the ModSim
library.
A TimeFrame is a kind of a DataFrame, which we used in.
Here's a more concise version of run_simulation using a TimeFrame
Step17:  The first line creates an empty TimeFrame with one column for each
state variable. Then, before the loop starts, we store the initial
conditions in the TimeFrame at t0. Based on the way we've been using
TimeSeries objects, it is tempting to write
Step18:  And plot the results like this
Step19:  As with a DataFrame, we can use the dot operator to select columns
from a TimeFrame.
Summary
Exercises
Exercise  Suppose the time between contacts is 4 days and the recovery time is 5 days.  After 14 weeks, how many students, total, have been infected?
Hint
<ASSISTANT_TASK:>
Python Code:
# install Pint if necessary
try:
    import pint
except ImportError:
    !pip install pint
# download modsim.py if necessary
from os.path import exists
filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)
# import functions from modsim
from modsim import *
Explanation: Chapter 11
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
init = State(S=89, I=1, R=0)
init
Explanation: In this chapter, we develop a model of an epidemic as it spreads in a
susceptible population, and use it to evaluate the effectiveness of
possible interventions.
My presentation of the model in the next few chapters is based on an excellent article by David Smith and Lang Moore, [^1]: Smith and Moore, "The SIR Model for Spread of Disease," Journal of Online Mathematics and its Applications, December 2001, available at http://modsimpy.com/sir.
The Freshman Plague
Every year at Olin College, about 90 new students come to campus from
around the country and the world. Most of them arrive healthy and happy, but usually at least one brings with them some kind of infectious disease. A few weeks later, predictably, some fraction of the incoming class comes down with what we call "The Freshman Plague".
In this chapter we introduce a well-known model of infectious disease,
the Kermack-McKendrick model, and use it to explain the progression of
the disease over the course of the semester, predict the effect of
possible interventions (like immunization) and design the most effective intervention campaign.
So far we have done our own modeling; that is, we've chosen physical
systems, identified factors that seem important, and made decisions
about how to represent them. In this chapter we start with an existing
model and reverse-engineer it. Along the way, we consider the modeling
decisions that went into it and identify its capabilities and
limitations.
The SIR model
The Kermack-McKendrick model is a simple version of an SIR model,
so-named because it considers three categories of people:
S: People who are "susceptible\", that is, capable of
    contracting the disease if they come into contact with someone who
    is infected.
I: People who are "infectious\", that is, capable of passing
    along the disease if they come into contact with someone
    susceptible.
R: People who are "recovered\". In the basic version of the
    model, people who have recovered are considered to be immune to
    reinfection. That is a reasonable model for some diseases, but not
    for others, so it should be on the list of assumptions to reconsider
    later.
Let's think about how the number of people in each category changes over time. Suppose we know that people with the disease are infectious for a period of 4 days, on average. If 100 people are infectious at a
particular point in time, and we ignore the particular time each one
became infected, we expect about 1 out of 4 to recover on any particular day.
Putting that a different way, if the time between recoveries is 4 days, the recovery rate is about 0.25 recoveries per day, which we'll denote with the Greek letter gamma, $\gamma$, or the variable name gamma.
If the total number of people in the population is $N$, and the fraction currently infectious is $i$, the total number of recoveries we expect per day is $\gamma i N$.
Now let's think about the number of new infections. Suppose we know that each susceptible person comes into contact with 1 person every 3 days, on average, in a way that would cause them to become infected if the other person is infected. We'll denote this contact rate with the Greek letter beta, $\beta$.
It's probably not reasonable to assume that we know $\beta$ ahead of
time, but later we'll see how to estimate it based on data from previous outbreaks.
If $s$ is the fraction of the population that's susceptible, $s N$ is
the number of susceptible people, $\beta s N$ is the number of contacts per day, and $\beta s i N$ is the number of those contacts where the other person is infectious.
In summary:
The number of recoveries we expect per day is $\gamma i N$; dividing by $N$ yields the fraction of the population that recovers in a day, which is $\gamma i$.
The number of new infections we expect per day is $\beta s i N$;
    dividing by $N$ yields the fraction of the population that gets
    infected in a day, which is $\beta s i$.
This model assumes that the population is closed; that is, no one
arrives or departs, so the size of the population, $N$, is constant.
The SIR equations
If we treat time as a continuous quantity, we can write differential
equations that describe the rates of change for $s$, $i$, and $r$ (where $r$ is the fraction of the population that has recovered):
$$\begin{aligned}
\frac{ds}{dt} &= -\beta s i \
\frac{di}{dt} &= \beta s i - \gamma i\
\frac{dr}{dt} &= \gamma i\end{aligned}$$ 
To avoid cluttering the equations, I leave it implied that $s$ is a function of time, $s(t)$, and likewise for $i$ and $r$.
SIR models are examples of compartment models, so-called because
they divide the world into discrete categories, or compartments, and
describe transitions from one compartment to another. Compartments are
also called stocks and transitions between them are called
flows.
In this example, there are three stocks---susceptible, infectious, and
recovered---and two flows---new infections and recoveries. Compartment
models are often represented visually using stock and flow diagrams (see http://modsimpy.com/stock).
The following figure shows the stock and flow diagram for an SIR
model.
{width="4in"}
Stocks are represented by rectangles, flows by arrows. The widget in the middle of the arrows represents a valve that controls the rate of flow; the diagram shows the parameters that control the valves.
Implementation
For a given physical system, there are many possible models, and for a
given model, there are many ways to represent it. For example, we can
represent an SIR model as a stock-and-flow diagram, as a set of
differential equations, or as a Python program. The process of
representing a model in these forms is called implementation. In
this section, we implement the SIR model in Python.
I'll represent the initial state of the system using a State object
with state variables S, I, and R; they represent the fraction of
the population in each compartment.
We can initialize the State object with the number of people in each compartment, assuming there is one infected student in a class of 90:
End of explanation
from numpy import sum
init /= sum(init)
init
Explanation: And then convert the numbers to fractions by dividing by the total:
End of explanation
tc = 3             # time between contacts in days 
tr = 4             # recovery time in days
Explanation: For now, let's assume we know the time between contacts and time between
recoveries:
End of explanation
beta = 1 / tc      # contact rate in per day
gamma = 1 / tr     # recovery rate in per day
Explanation: We can use them to compute the parameters of the model:
End of explanation
def make_system(beta, gamma):
    init = State(S=89, I=1, R=0)
    init /= sum(init)
    t0 = 0
    t_end = 7 * 14
    return System(init=init, t0=t0, t_end=t_end,
                  beta=beta, gamma=gamma)
Explanation: Now we need a System object to store the parameters and initial
conditions. The following function takes the system parameters and returns a new System object:
End of explanation
system = make_system(beta, gamma)
system
Explanation: The default value for t_end is 14 weeks, about the length of a
semester.
Here's what the System object looks like.
End of explanation
def update_func(state, t, system):
    s, i, r = state
    infected = system.beta * i * s    
    recovered = system.gamma * i
    
    s -= infected
    i += infected - recovered
    r += recovered
    
    return State(S=s, I=i, R=r)
Explanation: The update function
At any point in time, the state of the system is represented by a
State object with three variables, S, I and R. So I'll define an
update function that takes as parameters a State object, the current
time, and a System object:
End of explanation
state = update_func(init, 0, system)
state
Explanation: The first line uses a feature we have not seen before, multiple
assignment. The value on the right side is a State object that
contains three values. The left side is a sequence of three variable
names. The assignment does just what we want: it assigns the three
values from the State object to the three variables, in order.
The variables s, i and r, are lowercase to distinguish them
from the state variables, S, I and R.
The update function computes infected and recovered as a fraction of the population, then updates s, i and r. The return value is a State that contains the updated values.
When we call update_func like this:
End of explanation
from numpy import arange
def run_simulation(system, update_func):
    state = system.init
    for t in arange(system.t0, system.t_end):
        state = update_func(state, t, system)
    return state
Explanation: You might notice that this version of update_func does not use one of its parameters, t. I include it anyway because update functions
sometimes depend on time, and it is convenient if they all take the same parameters, whether they need them or not.
Running the simulation
Now we can simulate the model over a sequence of time steps:
End of explanation
system = make_system(beta, gamma)
final_state = run_simulation(system, update_func)
Explanation: The parameters of run_simulation are the System object and the
update function. The System object contains the parameters, initial
conditions, and values of t0 and t_end.
We can call run_simulation like this:
End of explanation
final_state
Explanation: The result is the final state of the system:
End of explanation
def run_simulation(system, update_func):
    S = TimeSeries()
    I = TimeSeries()
    R = TimeSeries()
    state = system.init
    t0 = system.t0
    S[t0], I[t0], R[t0] = state
    
    for t in arange(system.t0, system.t_end):
        state = update_func(state, t, system)
        S[t+1], I[t+1], R[t+1] = state
    
    return S, I, R
Explanation: This result indicates that after 14 weeks (98 days), about 52% of the
population is still susceptible, which means they were never infected,
less than 1% are actively infected, and 48% have recovered, which means they were infected at some point.
Collecting the results
The previous version of run_simulation only returns the final state,
but we might want to see how the state changes over time. We'll consider two ways to do that: first, using three TimeSeries objects, then using a new object called a TimeFrame.
Here's the first version:
End of explanation
system = make_system(beta, gamma)
S, I, R = run_simulation(system, update_func)
Explanation: First, we create TimeSeries objects to store the results. Notice that
the variables S, I, and R are TimeSeries objects now.
Next we initialize state, t0, and the first elements of S, I and
R.
Inside the loop, we use update_func to compute the state of the system
at the next time step, then use multiple assignment to unpack the
elements of state, assigning each to the corresponding TimeSeries.
At the end of the function, we return the values S, I, and R. This
is the first example we have seen where a function returns more than one
value.
Now we can run the function like this:
End of explanation
def plot_results(S, I, R):
    S.plot(style='--', label='Susceptible')
    I.plot(style='-', label='Infected')
    R.plot(style=':', label='Resistant')
    decorate(xlabel='Time (days)',
             ylabel='Fraction of population')
Explanation: We'll use the following function to plot the results:
End of explanation
plot_results(S, I, R)
Explanation: And run it like this:
End of explanation
def run_simulation(system, update_func):
    frame = TimeFrame(columns=system.init.index)
    frame.loc[system.t0] = system.init
    
    for t in arange(system.t0, system.t_end):
        frame.loc[t+1] = update_func(frame.loc[t], t, system)
    
    return frame
Explanation: Notice that it takes about three weeks (21 days) for the outbreak to get going, and about six weeks (42 days) before it peaks. The fraction of the population that's infected is never very high, but it adds up. In total, almost half the population gets sick.
Now with a TimeFrame
If the number of state variables is small, storing them as separate
TimeSeries objects might not be so bad. But a better alternative is to use a TimeFrame, which is another object defined in the ModSim
library.
A TimeFrame is a kind of a DataFrame, which we used in.
Here's a more concise version of run_simulation using a TimeFrame:
End of explanation
results = run_simulation(system, update_func)
Explanation: The first line creates an empty TimeFrame with one column for each
state variable. Then, before the loop starts, we store the initial
conditions in the TimeFrame at t0. Based on the way we've been using
TimeSeries objects, it is tempting to write:
frame[system.t0] = system.init
But when you use the bracket operator with a TimeFrame or DataFrame, it selects a column, not a row. 
To select a row, we have to use loc, like this:
frame.loc[system.t0] = system.init
Since the value on the right side is a State, the assignment matches
up the index of the State with the columns of the TimeFrame; that
is, it assigns the S value from system.init to the S column of
frame, and likewise with I and R.
We use the same feature to write the loop more concisely, assigning the State we get from update_func directly to the next row of
frame.
Finally, we return frame. We can call this version of run_simulation like this:
End of explanation
plot_results(results.S, results.I, results.R)
Explanation: And plot the results like this:
End of explanation
# Solution
tc = 4      # time between contacts in days 
tr = 5      # recovery time in days
beta = 1 / tc      # contact rate in per day
gamma = 1 / tr     # recovery rate in per day
system = make_system(beta, gamma)
s_0 = system.init.S
final = run_simulation(system, update_func)
s_end = final.S[system.t_end]
s_0 - s_end
Explanation: As with a DataFrame, we can use the dot operator to select columns
from a TimeFrame.
Summary
Exercises
Exercise  Suppose the time between contacts is 4 days and the recovery time is 5 days.  After 14 weeks, how many students, total, have been infected?
Hint: what is the change in S between the beginning and the end of the simulation?
End of explanation
<END_TASK>