<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Modeling and Simulation in Python
Chapter 14
Copyright 2017 Allen Downey
License
Step7:  Code from previous chapters
Step8:  Contact number
Here's the SweepFrame from the previous chapter, with one row for each value of beta and one column for each value of gamma.
Step9:  The following loop shows how we can loop through the columns and rows of the SweepFrame.  With 11 rows and 4 columns, there are 44 elements.
Step11:  Now we can wrap that loop in a function and plot the results.  For each element of the SweepFrame, we have beta, gamma, and frac_infected, and we plot beta/gamma on the x-axis and frac_infected on the y-axis.
Step12:  Here's what it looks like
Step13:  It turns out that the ratio beta/gamma, called the "contact number" is sufficient to predict the total number of infections; we don't have to know beta and gamma separately.
We can see that in the previous plot
Step14:  total_infected is the change in $s$ from the beginning to the end.
Step15:  Now we can plot the analytic results and compare them to the simulations.
Step16:  The agreement is generally good, except for values of c less than 1.
Exercises
Exercise
Step17:  Exercise
<ASSISTANT_TASK:>
Python Code:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *
Explanation: Modeling and Simulation in Python
Chapter 14
Copyright 2017 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
def make_system(beta, gamma):
    Make a system object for the SIR model.
    
    beta: contact rate in days
    gamma: recovery rate in days
    
    returns: System object
    
    init = State(S=89, I=1, R=0)
    init /= np.sum(init)
    t0 = 0
    t_end = 7 * 14
    return System(init=init, t0=t0, t_end=t_end,
                  beta=beta, gamma=gamma)
def update_func(state, t, system):
    Update the SIR model.
    
    state: State (s, i, r)
    t: time
    system: System object
    
    returns: State (sir)
    
    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)
def run_simulation(system, update_func):
    Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    
    init, t0, t_end = system.init, system.t0, system.t_end
    
    frame = TimeFrame(columns=init.index)
    frame.row[t0] = init
    
    for t in linrange(t0, t_end):
        frame.row[t+1] = update_func(frame.row[t], t, system)
    
    return frame
def calc_total_infected(results):
    Fraction of population infected during the simulation.
    
    results: DataFrame with columns S, I, R
    
    returns: fraction of population
    
    return get_first_value(results.S) - get_last_value(results.S)
def sweep_beta(beta_array, gamma):
    Sweep a range of values for beta.
    
    beta_array: array of beta values
    gamma: recovery rate
    
    returns: SweepSeries that maps from beta to total infected
    
    sweep = SweepSeries()
    for beta in beta_array:
        system = make_system(beta, gamma)
        results = run_simulation(system, update_func)
        sweep[system.beta] = calc_total_infected(results)
    return sweep
def sweep_parameters(beta_array, gamma_array):
    Sweep a range of values for beta and gamma.
    
    beta_array: array of infection rates
    gamma_array: array of recovery rates
    
    returns: SweepFrame with one row for each beta
             and one column for each gamma
    
    frame = SweepFrame(columns=gamma_array)
    for gamma in gamma_array:
        frame[gamma] = sweep_beta(beta_array, gamma)
    return frame
Explanation: Code from previous chapters
End of explanation
beta_array = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 , 1.1]
gamma_array = [0.2, 0.4, 0.6, 0.8]
frame = sweep_parameters(beta_array, gamma_array)
frame.head()
frame.shape
Explanation: Contact number
Here's the SweepFrame from the previous chapter, with one row for each value of beta and one column for each value of gamma.
End of explanation
for gamma in frame.columns:
    column = frame[gamma]
    for beta in column.index:
        frac_infected = column[beta]
        print(beta, gamma, frac_infected)
Explanation: The following loop shows how we can loop through the columns and rows of the SweepFrame.  With 11 rows and 4 columns, there are 44 elements.
End of explanation
def plot_sweep_frame(frame):
    Plot the values from a SweepFrame.
    
    For each (beta, gamma), compute the contact number,
    beta/gamma
    
    frame: SweepFrame with one row per beta, one column per gamma
    
    for gamma in frame.columns:
        column = frame[gamma]
        for beta in column.index:
            frac_infected = column[beta]
            plot(beta/gamma, frac_infected, 'ro')
Explanation: Now we can wrap that loop in a function and plot the results.  For each element of the SweepFrame, we have beta, gamma, and frac_infected, and we plot beta/gamma on the x-axis and frac_infected on the y-axis.
End of explanation
plot_sweep_frame(frame)
decorate(xlabel='Contact number (beta/gamma)',
         ylabel='Fraction infected')
savefig('figs/chap14-fig01.pdf')
Explanation: Here's what it looks like:
End of explanation
s_inf_array = linspace(0.0001, 0.9999, 101);
c_array = log(s_inf_array) / (s_inf_array - 1);
Explanation: It turns out that the ratio beta/gamma, called the "contact number" is sufficient to predict the total number of infections; we don't have to know beta and gamma separately.
We can see that in the previous plot: when we plot the fraction infected versus the contact number, the results fall close to a curve.
Analysis
In the book we figured out the relationship between $c$ and $s_{\infty}$ analytically.  Now we can compute it for a range of values:
End of explanation
frac_infected = 1 - s_inf_array
frac_infected_series = Series(frac_infected, index=c_array);
Explanation: total_infected is the change in $s$ from the beginning to the end.
End of explanation
plot_sweep_frame(frame)
plot(frac_infected_series, label='Analysis')
decorate(xlabel='Contact number (c)',
         ylabel='Fraction infected')
savefig('figs/chap14-fig02.pdf')
Explanation: Now we can plot the analytic results and compare them to the simulations.
End of explanation
# Solution
def plot_sweep_frame_difference(frame):
    for gamma in frame.columns:
        column = frame[gamma]
        for beta in column.index:
            frac_infected = column[beta]
            plot(beta - gamma, frac_infected, 'ro')
# Solution
plot_sweep_frame_difference(frame)
decorate(xlabel='Excess infection rate (infections-recoveries per day)',
         ylabel='Fraction infected',
         legend=False)
# Solution
# The results don't fall on a line, which means that if we know the difference between
# `beta` and `gamma`, but not their ratio, that's not enough to predict the fraction infected.
Explanation: The agreement is generally good, except for values of c less than 1.
Exercises
Exercise:  If we didn't know about contact numbers, we might have explored other possibilities, like the difference between beta and gamma, rather than their ratio.
Write a version of plot_sweep_frame, called plot_sweep_frame_difference, that plots the fraction infected versus the difference beta-gamma.
What do the results look like, and what does that imply?
End of explanation
# Solution
frac_infected_series
# Solution
# It looks like the fraction infected is 0.26 when the contact number is about 1.16
Explanation: Exercise: Suppose you run a survey at the end of the semester and find that 26% of students had the Freshman Plague at some point.
What is your best estimate of c?
Hint: if you print frac_infected_series, you can read off the answer.
End of explanation
<END_TASK>