Assignment 4: Modeling Epidemics¶
import time
import EoN as eon
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import scipy
from tqdm.notebook import tqdm
Part 1: Outbreak Modeling [40 Points]¶
def load_flu_network():
# Read the graph
G = nx.read_edgelist(“fludata.txt”, nodetype=int, data=((“weight”, float),))
def simulate_outbreak(G, n_iter, initial_infected=325, tmax=10, beta=0.01, mu=0.5):
G : nx.Graph
n_iter : int
Number of simulations to run.
initial_infected : int
The node that is initially infected
tmax : int
beta : float
mu : float
simulation_runs : list[tuple] –
This is a list like object of tuples of the
form (t, S, I). Each tuple (t, S, I) represents
a simulated run, where t, S, I are vectors representing
the time periods, and the corresponding numbers of susceptible
and infected people.
simulation_runs = []
return simulation_runs
def plot_outbreaks(simulation_runs, save=False):
simulation_runs : list[tuple]
The list like object of tuples returned by `simulate_outbreak`.
save : bool
Whether to save the figure.
plt.savefig(“1_1.png”)
plt.show()
def get_exponent(simulation_run, I_thrsh=100):
simulation_run : tuple
A tuple of (t, S, I) as specified in `simulate_outbreak()`
I_thrsh : int
Threshold of I. We only fit the curve where I<=I_thrsh
tau : float
return tau
def plot_curve_fit(simulation_run, tau, I_thrsh=100, save=False):
simulation_run : tuple
I_thrsh : int
tau : float
save : bool
r2 = 0.0 # R^2 to be computed
plt.savefig("1_2.png")
plt.show()
def calculate_theoretical_taus(G, beta=0.01, mu=0.5):
G : nx.Graph
beta : float
mu : float
tau_rand : float
tau_slide : float
tau_book : float
tau_rand = tau_slide = tau_book = 0.0
return tau_rand, tau_slide, tau_book
def compare_taus(empirical_taus, tau_rand, tau_slide, tau_book, save=False):
This function visually compares the empirical_taus, tau_rand, tau_slide,
and tau_book on a boxplot.
empirical_taus : list[float]
List-like object of floats containing the list of tau's
tau_rand : float
tau_slide : float
tau_book : float
save : bool
plt.savefig("1_3.png")
plt.show()
def calculate_theoretical_endemic_size(G, beta=0.01, mu=0.5):
G : nx.Graph
beta : float
mu : float
theoretical_endemic_size : int
theoretical_endemic_size = 0.0
return theoretical_endemic_size
def compare_endemic_sizes(
empirical_endemic_sizes, theoretical_endemic_size, save=False
empirical_endemic_sizes : list[int]
List-like object of int's representing a distribution of
empirical endemic sizes from multiple simulation runs.
Refer to cell 1.5 to see how it works.
theoretical_endemiz_size : int
save : bool
plt.savefig("1_4.png")
plt.show()
# Do not modify
print(">>>>> Results for Part 1 <<<<<")
G = load_flu_network()
# Generate 10 simulation runs for plotting
simulation_runs = simulate_outbreak(G, 10)
plot_outbreaks(simulation_runs)
tau = get_exponent(simulation_runs[0])
r2 = plot_curve_fit(simulation_runs[0], tau)
print("> Part 1.2 <")
print(f"tau={tau:.2f}, R2={r2:.4f}")
simulation_runs = simulate_outbreak(G, 25)
empirical_taus = []
empirical_endemic_sizes = []
for run in simulation_runs:
empirical_taus.append(get_exponent(run))
empirical_endemic_sizes.append(run[2][-1])
tau_rand, tau_slide, tau_book = calculate_theoretical_taus(G)
print("\n> Part 1.3 <")
print(f"tau_rand={tau_rand:.4f}, tau_slide={tau_slide:.4f}, tau_book={tau_book:.4f}")
print("Empirical distribution of tau")
def describe(x):
"Count={count}\tMean={mean:.3f}\tstd={median:.3f}\n"
"Min={min:.3f}\t25%={p25:.3f}\t50%={p50:.3f}\t75%={p75:<.3f}\tMax={max:.3f}".format(
count=len(x),
mean=np.mean(x),
median=np.std(x),
min=min(x),
p25=np.percentile(x, 25),
p50=np.median(x),
p75=np.percentile(x, 75),
max=max(x),
describe(empirical_taus)
compare_taus(empirical_taus, tau_rand, tau_slide, tau_book)
print("\n> Part 1.4 <")
theoretical_endemic_size = calculate_theoretical_endemic_size(G)
compare_endemic_sizes(empirical_endemic_sizes, theoretical_endemic_size)
print(f"Theoretical endemic size={theoretical_endemic_size:.2f}")
print("Empirical distribution of simulated endemic size")
describe(empirical_endemic_sizes)
1.5 Written Response¶
Part 2: Transmission Rate Variation $\beta$ [25 Points]¶
2.1 Minimum Transmission Rate for Epidemic¶
def simulate_beta_sweep(
beta_min=0.001,
beta_max=0.04,
beta_samples=40,
initial_infected=325,
Generate a list of betas. Run multiple simulations for each beta.
Save the results.
G : nx.Graph
n_sims : int
Number of simulations (or runs) for each beta value
beta_min : float
Minimum beta to simulate
beta_max : float
Maximum beta to simulate
beta_samples : int
The number of betas to simulate. That is, the function will
generate `beta_samples` betas between `beta_min` and `beta_max`.
initial_infect : int
Initial infected node
tmax : int
mu : float
betas : list[float]
The list of betas the function has generated.
beta_runs : list[list[tuple]]
Simulation results corresponding to `betas`. It is a list where
each itemis a list of tuples, representing multiple simulations
for a particular beta. The tuple is in the form of (t, S, I) as
in `simulate_outbreak` in Part 1.
betas = []
beta_runs = []
return betas, beta_runs
def extract_average_tau(beta_runs):
Estimate the average tau value for each beta value.
beta_runs : list[list[tuple]]
See docstring of `simulate_beta_sweep` above.
avg_taus : list[float]
A list like object, where each item is the estimated
average tau values of the corresponding simulations in
`beta_runs` for a particular beta.
avg_taus = []
return avg_taus
def plot_beta_tau_curves(betas, avg_taus, t, save=False):
betas : list[float]
avg_taus : list[float]
t : list[float]
save: bool
plt.savefig("2_1.png")
plt.show()
def extract_average_endemic_size(beta_runs):
Given the simulated data, compute the average size of the endemic state for
each beta value using the simulation data in `beta_runs`.
beta_runs : list[list[tuple]]
See docstring of `simulate_beta_sweep()` above.
avg_ends : list[float]
A list-like object of floats. Each item is the average endemic size
for a particular beta value estimated using the corresponding simulation
results in `beta_runs`.
avg_endemic_sizes = []
return avg_endemic_sizes
def calculate_theoretical_endemic(G, betas, mu=0.5):
Compute the theoretical endemic size for each beta value in `betas`,
assuming random distribution.
Additionally, compute the minimum theoretical beta for
epdemic to occur, assuming random distribution and
arbitrary distribution respectively.
G : nx.Graph
betas : list[float]
List of betas used to compute `theoretical_endemics`.
Each item of the returned `theoretical_endemics` corresponds
to a beta value in the `betas` list at the same index.
mu : float
theoretical_endemics : list[float]
A list-like object of theoretical endemic sizes
corresponding to the betas in the `betas` parameter.
rand_dist_min_beta : float
arb_dist_min_beta : float
theoretical_endemic_sizes = []
rand_dist_min_beta = arb_dist_min_beta = 0.0
return theoretical_endemic_sizes, rand_dist_min_beta, arb_dist_min_beta
def compare_endemic_sizes_vs_beta(
avg_endemic_sizes,
theoretical_endemic_sizes,
rand_dist_min_beta,
arb_dist_min_beta,
save=False,
Plot a figure to visually compare the average endemic sizes from simulation vs
the theoretical endemic sizes for different beta values. Also, it shows the minimum
beta thresholds for an endemic to occur assuming random and arbitrary distributions
respectively.
betas : list[float]
A list-like object of betas corresponding to `avg_endemic_sizes` and `theoretical_endemic_sizes`
avg_endemic_sizes : list[float]
The average endemic sizes as returned by
`extract_average_endemic_size`. Each element
corresponds to a beta value in the `betas` parameter
at the same index.
theoretical_endemic_sizes : list[float]
The theoretical endemic sizes. Each element
corresponds to a beta value in the `betas`
parameter at the same index.
rand_dist_min_beta : float
arb_dist_min_beta : float
save : bool
plt.savefig("2_2.png")
plt.show()
# do not modify
G = load_flu_network()
betas, beta_runs = simulate_beta_sweep(G, 5, beta_samples=20)
avg_taus = extract_average_tau(beta_runs)
times = np.linspace(0, 2.2, 100)
plot_beta_tau_curves(betas, avg_taus, t=times)
avg_endemic_sizes = extract_average_endemic_size(beta_runs)
theoretical_endemics_sizes,
rand_dist_min_beta,
arb_dist_min_beta,
) = calculate_theoretical_endemic(G, betas)
compare_endemic_sizes_vs_beta(
avg_endemic_sizes,
theoretical_endemics_sizes,
rand_dist_min_beta,
arb_dist_min_beta,
# print results
print(">>>>> Results for Part 2 <<<<<")
with np.printoptions(precision=2, suppress=True):
print(f"Avg_taus = {np.array(avg_taus)}\n")
with np.printoptions(precision=0, suppress=True):
print(f"Avg endemic size = {np.array(avg_endemic_sizes)}\n")
print(f"Theo endemic size = {np.array(theoretical_endemics_sizes)}\n")
print(f"Min beta for random distribution = {rand_dist_min_beta:.5f}")
print(f"Min beta for arbitrary distribution = {arb_dist_min_beta:.5f}")
2.3 Written Response¶
def sweep_initial_infected(G, tmax=10, beta=0.01, mu=0.5):
Compute the taus for each node.
G : nx.Graph
tmax : int
beta : float
mu : float
taus : list[float]
List of tau's corresponding to the `nodes` return value below.
nodes : list[int]
nodes = []
return taus, nodes
# Compute centrality metrics
def compute_centrality(G, nodes):
G : nx.Graph
nodes : list[int]
A list of nodes for which we compute centralities.
cent_dict : dict[list[float]]
The keys of the dict are 'deg', 'clo', 'bet' and 'eig'.
The values are lists of floats representing the corresponding
centrality for each node in the `nodes` parameter.
deg_cen = []
clo_cen = []
bet_cen = []
eig_cen = []
cent_dict = {"deg": deg_cen, "clo": clo_cen, "bet": bet_cen, "eig": eig_cen}
return cent_dict
def calculate_person_correlation(taus, cent_dict):
taus : list[float]
cent_dict : dict[list[float]]
See docstring in `compute_centrality`
r_dict : dict[tuple[float, float]]
The keys of the dict are 'deg', 'clo', 'bet' and 'eig'.
The values are tuples of (Pearson coefficient, p-value).
r_deg = r_clo = r_bet = r_eig = (
0.0, # Pearson correlation coefficient
0.0, # p-value
r_dict = {"deg": r_deg, "clo": r_clo, "bet": r_bet, "eig": r_eig}
return r_dict
def plot_centrality_vs_tau(taus, cent_dict, r_dict, save=False):
taus : list[float]
See docstring in `sweep_initial_infected`
cent_dict : dict[list[float]]
See docstring in `compute_centrality`
r_dict : dict[tuple[float, float]]
See docstring in `calculate_person_correlation`
save: bool
# Plot centrality metrics v/s tau - 4 plots
plt.savefig("3_2.png")
plt.show()
# do not modify
start_time = time.time()
G = load_flu_network()
taus, nodes = sweep_initial_infected(G)
cent_dict = compute_centrality(G, nodes)
r_dict = calculate_person_correlation(taus, cent_dict)
plot_centrality_vs_tau(taus, cent_dict, r_dict)
print(f"Number of included nodes = {len(nodes)}")
for k in sorted(r_dict):
coeff, pvalue = r_dict[k]
print(f"{k} = {coeff:.4f}, pvalue = {pvalue:.4f}")
# We don't grade by how long it takes. This is purly informational.
seconds_elapsed = time.time() - start_time
print(f"Elapsed time = {seconds_elapsed/60 : .2f} minutes")
3.3 Written Response¶
Part 4: Knowledge Question [5 Points]¶
(You can do this proof as markdown here or upload an image of the proof on paper. If you upload an image make sure to include the image file with your submission)