From 506699b789b10cfcd15a471ce6d9c6b32b1543aa Mon Sep 17 00:00:00 2001 From: Lucas Boettcher Date: Mon, 6 Jul 2020 08:29:52 -0700 Subject: [PATCH] final paper examples --- examples/paper/seeds_final_examples.txt | 20 + .../paper/simulate_idealized_NYC_epidemic.py | 47 +-- ...e_idealized_NYC_epidemic_nointervention.py | 8 +- figs/NYC_plot.py | 351 +++++++++++++----- figs/NYC_plot_nointervention.py | 192 ---------- 5 files changed, 285 insertions(+), 333 deletions(-) create mode 100644 examples/paper/seeds_final_examples.txt delete mode 100644 figs/NYC_plot_nointervention.py diff --git a/examples/paper/seeds_final_examples.txt b/examples/paper/seeds_final_examples.txt new file mode 100644 index 00000000..6fa89020 --- /dev/null +++ b/examples/paper/seeds_final_examples.txt @@ -0,0 +1,20 @@ +112312231 +982525 +12369433 +9478123 +2314675 +123123957 +71278 +2394194 +123412 +98211034 +84912937 +1231 +746184904 +1231414 +716488130 +8165143 +19577134 +93483 +2342315 +914235813 diff --git a/examples/paper/simulate_idealized_NYC_epidemic.py b/examples/paper/simulate_idealized_NYC_epidemic.py index 7e6f3c63..7d0b5247 100644 --- a/examples/paper/simulate_idealized_NYC_epidemic.py +++ b/examples/paper/simulate_idealized_NYC_epidemic.py @@ -67,7 +67,7 @@ set_num_threads(1) # Set random seeds for reproducibility -seed = 942395 +seed = 1231414#1285854#1237102#123746 seed_three_random_states(seed) # contact network ############################################################## @@ -119,7 +119,7 @@ statuses = random_epidemic( network.get_node_count(), network.get_nodes(), - fraction_infected=0.0025) + fraction_infected=0.0015) epidemic_simulator.set_statuses(statuses) @@ -129,7 +129,7 @@ # set the new contact rates on the network # run the kinetic model [kinetic produces the current statuses used as data] network = epidemic_simulator.run( - stop_time = epidemic_simulator.time + 18, + stop_time = epidemic_simulator.time + 19, current_network = network) kinetic_model = epidemic_simulator.kinetic_model @@ -141,40 +141,7 @@ λ_min = 2 # minimum contact rate λ_max = 8 # maximum contact rate -health_service = epidemic_simulator.health_service - -statuses = epidemic_simulator.kinetic_model.current_statuses - -epidemic_simulator = EpidemicSimulator( - network, - community_transmission_rate = community_transmission_rate['mean'], - hospital_transmission_reduction = hospital_transmission_reduction, - static_contact_interval = static_contact_interval, - mean_contact_lifetime = mean_contact_lifetime, - day_inception_rate = λ_max, - night_inception_rate = λ_min, - health_service = health_service, - start_time = epidemic_simulator.time) - -epidemic_simulator.set_statuses(statuses) - -# run the kinetic model -network = epidemic_simulator.run( - stop_time = epidemic_simulator.time + 83, - current_network = network) - -for key in kinetic_model.statuses.keys(): - kinetic_model.statuses[key].extend(epidemic_simulator.kinetic_model.statuses[key]) - - -kinetic_model.times.extend(kinetic_model.times[-1]+epidemic_simulator.kinetic_model.times) - -################################################################################ -# loosening SD intervention #################################################### -################################################################################ - -λ_min = 2 # minimum contact rate -λ_max = 22 # maximum contact rate +network.set_lambdas(λ_min, λ_max) health_service = epidemic_simulator.health_service @@ -195,20 +162,20 @@ # run the kinetic model network = epidemic_simulator.run( - stop_time = epidemic_simulator.time + 90, + stop_time = epidemic_simulator.time + 112, current_network = network) for key in kinetic_model.statuses.keys(): kinetic_model.statuses[key].extend(epidemic_simulator.kinetic_model.statuses[key]) -kinetic_model.times.extend(kinetic_model.times[-1]+epidemic_simulator.kinetic_model.times) +kinetic_model.times.extend(epidemic_simulator.kinetic_model.times) ################################################################################ # save ######################################################################### ################################################################################ if SAVE_FLAG: np.savetxt( - os.path.join(SIMULATION_PATH, 'NYC_interventions_1e5.txt'), + os.path.join(SIMULATION_PATH, 'NYC_interventions_1e5_13.txt'), np.c_[ kinetic_model.times, kinetic_model.statuses['S'], diff --git a/examples/paper/simulate_idealized_NYC_epidemic_nointervention.py b/examples/paper/simulate_idealized_NYC_epidemic_nointervention.py index 8785bdbe..03ead053 100644 --- a/examples/paper/simulate_idealized_NYC_epidemic_nointervention.py +++ b/examples/paper/simulate_idealized_NYC_epidemic_nointervention.py @@ -67,7 +67,7 @@ set_num_threads(1) # Set random seeds for reproducibility -seed = 942395 +seed = 914235813 seed_three_random_states(seed) # contact network ############################################################## @@ -119,7 +119,7 @@ statuses = random_epidemic( network.get_node_count(), network.get_nodes(), - fraction_infected=0.0025) + fraction_infected=0.0015) epidemic_simulator.set_statuses(statuses) @@ -129,7 +129,7 @@ # set the new contact rates on the network # run the kinetic model [kinetic produces the current statuses used as data] network = epidemic_simulator.run( - stop_time = epidemic_simulator.time + 185, + stop_time = epidemic_simulator.time + 132, current_network = network) kinetic_model = epidemic_simulator.kinetic_model @@ -139,7 +139,7 @@ ################################################################################ if SAVE_FLAG: np.savetxt( - os.path.join(SIMULATION_PATH, 'NYC_nointerventions_1e5.txt'), + os.path.join(SIMULATION_PATH, 'NYC_nointerventions_1e5_19.txt'), np.c_[ kinetic_model.times, kinetic_model.statuses['S'], diff --git a/figs/NYC_plot.py b/figs/NYC_plot.py index a2d52a9d..f424875a 100644 --- a/figs/NYC_plot.py +++ b/figs/NYC_plot.py @@ -1,15 +1,12 @@ import os, sys; sys.path.append(os.path.join("..")) -from timeit import default_timer as timer - -import networkx as nx import numpy as np import pandas as pd -import random import datetime as dt import matplotlib.dates as mdates import matplotlib.ticker as ticker import matplotlib.pyplot as plt +from epiforecast.contact_network import ContactNetwork import seaborn as sns @@ -43,8 +40,6 @@ fig_size = [1.5*fig_width, 2*fig_height] rcParams.update({'figure.figsize': fig_size}) -from epiforecast.scenarios import load_edges - def simulation_average(model_data, times, sampling_time = 1): """ Returns daily averages of simulation data. @@ -73,24 +68,14 @@ def simulation_average(model_data, times, sampling_time = 1): return daily_average #%% load network data to get the population size -edges = load_edges(os.path.join('..', 'data', 'networks', 'edge_list_SBM_1e4_nobeds.txt')) +NETWORKS_PATH = os.path.join('..', 'data', 'networks') +edges_filename = os.path.join(NETWORKS_PATH, 'edge_list_SBM_1e5_nobeds.txt') +groups_filename = os.path.join(NETWORKS_PATH, 'node_groups_SBM_1e5_nobeds.json') -contact_network = nx.Graph() -contact_network.add_edges_from(edges) -contact_network = nx.convert_node_labels_to_integers(contact_network) -population = len(contact_network) +network = ContactNetwork.from_files(edges_filename, groups_filename) +population = len(network.graph) #%% -#%% load simulation data (with and without interventions) -simulation_data = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'simulation_data_NYC_1e4_6.txt')) -simulation_data_nointervention = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'simulation_data_nointervention2.txt')) - -times = simulation_data[:,0] -times_nointervention = simulation_data_nointervention[:,0] - -kinetic_model = {'S': simulation_data[:,1], 'E': simulation_data[:,2], 'I': simulation_data[:,3], 'R': simulation_data[:,4], 'H': simulation_data[:,5], 'D': simulation_data[:,6]} -kinetic_model_nointervention = {'S': simulation_data_nointervention[:,1], 'E': simulation_data_nointervention[:,2], 'I': simulation_data_nointervention[:,3], 'R': simulation_data_nointervention[:,4], 'H': simulation_data_nointervention[:,5], 'D': simulation_data_nointervention[:,6]} - #%% load NYC data NYC_data = pd.read_csv(os.path.join('..', 'data', 'NYC_COVID_CASES', 'data_new_york.csv')) NYC_cases = np.asarray([float(x) for x in NYC_data['Cases'].tolist()]) @@ -116,44 +101,11 @@ def simulation_average(model_data, times, sampling_time = 1): NYC_cases_weekly = np.mean(np.append(reported_cases_NYC, (7-len(reported_cases_NYC)%7)*[reported_cases_NYC[-1]]).reshape(-1, 7), axis=1) -#%% determine averages of simulation data -# daily averages of simulation data -# sampling_time = 1 means that we average over 1-day intervals -sampling_time = 7 -daily_average = simulation_average(kinetic_model, times, sampling_time) -cumulative_cases_simulation = (1-np.asarray(daily_average['S'])/population)*1e5 -daily_average = simulation_average(kinetic_model, times, sampling_time) -cumulative_deaths_simulation = np.asarray(daily_average['D'])/population*1e5 - -daily_average_nointervention = simulation_average(kinetic_model_nointervention, times_nointervention, sampling_time) -cumulative_cases_simulation_nointervention = (1-np.asarray(daily_average_nointervention['S'])/population)*1e5 -cumulative_deaths_simulation_nointervention = np.asarray(daily_average_nointervention['D'])/population*1e5 - -cases_simulation = np.ediff1d(cumulative_cases_simulation)/sampling_time -cases_simulation_nointervention = np.ediff1d(cumulative_cases_simulation_nointervention)/sampling_time - -simulation_dates = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time) for i in range(len(cumulative_cases_simulation))]) -simulation_dates_nointervention = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time) for i in range(len(cumulative_deaths_simulation_nointervention))]) - -sampling_time2 = 7 -daily_average2 = simulation_average(kinetic_model, times, sampling_time2) -cumulative_cases_simulation2 = (1-np.asarray(daily_average2['S'])/population)*1e5 -daily_average2 = simulation_average(kinetic_model, times, sampling_time2) -cumulative_deaths_simulation2 = np.asarray(daily_average2['D'])/population*1e5 -daily_average_nointervention2 = simulation_average(kinetic_model_nointervention, times_nointervention, sampling_time2) -cumulative_cases_simulation_nointervention2 = (1-np.asarray(daily_average_nointervention2['S'])/population)*1e5 -cumulative_deaths_simulation_nointervention2 = np.asarray(daily_average_nointervention2['D'])/population*1e5 - -deaths_simulation = np.ediff1d(cumulative_deaths_simulation2)/sampling_time2 -deaths_simulation_nointervention = np.ediff1d(cumulative_deaths_simulation_nointervention2)/sampling_time2 - -simulation_dates2 = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time2) for i in range(len(deaths_simulation))]) -simulation_dates_nointervention2 = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time2) for i in range(len(deaths_simulation_nointervention))]) - -#%% plot results - +#%% plot definitions fig, axs = plt.subplots(nrows = 2, ncols = 2) +# initial simulation data shift +dd = 6 # days # cumulative death panel ax00 = axs[0][0] @@ -162,20 +114,18 @@ def simulation_average(model_data, times, sampling_time = 1): ax00.set_title(r'Deaths per 100,000') -ax00.text(dt.date(2020, 3, 12), 0.9*1000, r'(a)') - -ax00.plot(simulation_dates+dt.timedelta(days = 13), cumulative_deaths_simulation, 'cornflowerblue', alpha = 1) -ax00.plot(simulation_dates_nointervention+dt.timedelta(days = 13), cumulative_deaths_simulation_nointervention, 'Grey') -ax00.bar(NYC_data_date_of_interest_deaths, cumulative_deaths_NYC, facecolor='Grey', edgecolor='Grey', alpha = 0.6, width = 0.00001, align = 'center') +ax00.text(dt.date(2020, 2, 26), 0.9*800, r'(a)') +ax00.bar(NYC_data_date_of_interest_deaths, cumulative_deaths_NYC, facecolor='red', edgecolor='red', alpha = 0.4, width = 0.0001) + ax00.text(dt.date(2020, 4, 29), 75, r'data') -ax00.text(dt.date(2020, 6, 14), 300, r'model') -ax00.text(dt.date(2020, 4, 20), 900, r'no SD') +ax00.text(dt.date(2020, 6, 14), 350, r'model') +ax00.text(dt.date(2020, 4, 28), 650, r'no SD') ax00.set_ylabel("cumulative") -ax00.set_ylim(0,1000) -ax00.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax00.set_ylim(0,800) +ax00.set_xlim([dt.date(2020, 2, 23), dt.date(2020, 7, 12)]) ax00.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) ax00.xaxis.set_major_locator(ticker.MultipleLocator(14)) ax00.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) @@ -192,18 +142,16 @@ def simulation_average(model_data, times, sampling_time = 1): ax01_2 = axs[0][1].twinx() -ax01.text(dt.date(2020, 3, 12), 0.9*60000, r'(b)') +ax01.text(dt.date(2020, 2, 26), 0.9*60000, r'(b)') -ax01.plot(simulation_dates+dt.timedelta(days = 13), cumulative_cases_simulation, 'cornflowerblue') -ax01.plot(simulation_dates_nointervention+dt.timedelta(days = 13), cumulative_cases_simulation_nointervention, 'Grey') ax01_2.bar(NYC_date_of_interest_cases, cumulative_reported_cases_NYC, facecolor='red', edgecolor='red', alpha = 0.4, width = 0.00001, align = 'center') - + ax01.text(dt.date(2020, 4, 29), 6000, r'data') -ax01.text(dt.date(2020, 6, 11), 18000, r'model') -ax01.text(dt.date(2020, 4, 5), 55000, r'no SD') +ax01.text(dt.date(2020, 6, 11), 31000, r'model') +ax01.text(dt.date(2020, 4, 10), 52000, r'no SD') ax01.set_ylim(0,60000) -ax01.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax01.set_xlim([dt.date(2020, 2, 23), dt.date(2020, 7, 12)]) ax01.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) ax01.xaxis.set_major_locator(ticker.MultipleLocator(14)) ax01.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) @@ -211,7 +159,8 @@ def simulation_average(model_data, times, sampling_time = 1): ticker.FuncFormatter(lambda x, p: format(int(x), ','))) -ax01_2.set_ylim(0,9000) +ax01_2.set_ylim(0,6000) +ax01_2.set_yticks([0,2000,4000,6000]) ax01_2.tick_params(axis='y', colors='indianred') ax01_2.get_yaxis().set_major_formatter( ticker.FuncFormatter(lambda x, p: format(int(x), ','))) @@ -224,24 +173,22 @@ def simulation_average(model_data, times, sampling_time = 1): ax10_2 = axs[1][0].twinx() -ax10.text(dt.date(2020, 3, 12), 0.9*30, r'(c)') +ax10.text(dt.date(2020, 2, 26), 0.9*25, r'(c)') -ax10.plot(simulation_dates2+dt.timedelta(days = 13), deaths_simulation, 'cornflowerblue') -ax10.plot(simulation_dates_nointervention2+dt.timedelta(days = 13), deaths_simulation_nointervention, 'Grey') -ax10.fill_between(NYC_data_date_of_interest_deaths[::7]+dt.timedelta(days = 3.5), NYC_death_data_weekly, edgecolor = 'k', facecolor = 'Grey', alpha = 0.4, linewidth = 1.) -ax10.plot(NYC_data_date_of_interest_deaths[::7]+dt.timedelta(days = 3.5), NYC_death_data_weekly, color = 'k', linewidth = 1.) +ax10.fill_between(NYC_data_date_of_interest_deaths[::7]+dt.timedelta(days = 3.5), NYC_death_data_weekly, edgecolor = 'red', facecolor = 'red', alpha = 0.15, linewidth = 1.) +ax10.plot(NYC_data_date_of_interest_deaths[::7]+dt.timedelta(days = 3.5), NYC_death_data_weekly, color = 'red', linewidth = 1.) -ax10.bar(NYC_data_date_of_interest_deaths, reported_deaths_NYC, facecolor='Grey', edgecolor='Grey', alpha = 0.5, width = 0.0001) +ax10.bar(NYC_data_date_of_interest_deaths, reported_deaths_NYC, facecolor='red', edgecolor='red', alpha = 0.4, width = 0.0001) -ax10.text(dt.date(2020, 4, 3), 10.5, r'new deaths', fontsize = 7) +ax10.text(dt.date(2020, 4, 7), 9.5, r'new deaths', fontsize = 7) ax10.text(dt.date(2020, 5, 18), 8, r'7-day average', fontsize = 7) ax10.plot([dt.date(2020, 5, 16), dt.date(2020, 4, 28)], [8.0, 3.2], color = 'k', linewidth = 0.5) ax10.set_ylabel("daily") -ax10.set_ylim(0,30) -ax10.set_yticks([0,10,20,30]) -ax10.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax10.set_ylim(0,25) +#ax10.set_yticks([0,10,20]) +ax10.set_xlim([dt.date(2020, 2, 23), dt.date(2020, 7, 12)]) ax10.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) ax10.xaxis.set_major_locator(ticker.MultipleLocator(14)) ax10.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) @@ -257,37 +204,247 @@ def simulation_average(model_data, times, sampling_time = 1): ax11_2 = axs[1][1].twinx() -ax11.text(dt.date(2020, 3, 12), 0.9*2000, r'(d)') +ax11.text(dt.date(2020, 2, 26), 0.9*2000, r'(d)') -ax11.plot(simulation_dates2+dt.timedelta(days = 13), cases_simulation, 'cornflowerblue') -ax11.plot(simulation_dates_nointervention2+dt.timedelta(days = 13), cases_simulation_nointervention, 'Grey') -ax11_2.fill_between(NYC_date_of_interest_cases[::7]+dt.timedelta(days = 3.5), NYC_cases_weekly, edgecolor = 'red', facecolor = 'red', alpha = 0.15, linewidth = 1.) +ax11_2.fill_between(NYC_date_of_interest_cases[::7]+dt.timedelta(days = 3.5), NYC_cases_weekly, edgecolor = 'red', facecolor = 'red', alpha = 0.15, linewidth = 1.) ax11_2.plot(NYC_date_of_interest_cases[::7]+dt.timedelta(days = 3.5), NYC_cases_weekly, color = 'red', linewidth = 1.) ax11_2.bar(NYC_date_of_interest_cases, reported_cases_NYC, facecolor='red', edgecolor='red', alpha = 0.4, width = 0.0001) -ax11.text(dt.date(2020, 3, 26), 550, r'new cases', fontsize = 7) +ax11.text(dt.date(2020, 3, 29), 800, r'new cases', fontsize = 7) ax11.text(dt.date(2020, 5, 17), 400, r'7-day average', fontsize = 7) -ax11.plot([dt.date(2020, 5, 16), dt.date(2020, 4, 23)], [420, 240], color = 'k', linewidth = 0.5) +ax11.plot([dt.date(2020, 5, 16), dt.date(2020, 4, 28)], [420, 240], color = 'k', linewidth = 0.5) ax11.set_ylim(0,2050) -ax11.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax11.set_xlim([dt.date(2020, 2, 23), dt.date(2020, 7, 12)]) ax11.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) ax11.xaxis.set_major_locator(ticker.MultipleLocator(14)) ax11.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) ax11.get_yaxis().set_major_formatter( - ticker.FuncFormatter(lambda x, p: format(int(x), ','))) - -ax11_2.set_ylim(0,300) -ax11_2.set_yticks([0,100,200,300]) +ticker.FuncFormatter(lambda x, p: format(int(x), ','))) + +ax11_2.set_ylim(0,2000*6000/60000) +ax11_2.set_yticks([0,100,200]) ax11_2.tick_params(axis='y', colors='indianred') ax11_2.get_yaxis().set_major_formatter( - ticker.FuncFormatter(lambda x, p: format(int(x), ','))) - +ticker.FuncFormatter(lambda x, p: format(int(x), ','))) + ax11.yaxis.grid() +cumulative_deaths_simulation_mean_arr = [] +cumulative_deaths_simulation_nointervention_mean_arr = [] +cumulative_cases_simulation_mean_arr = [] +cumulative_cases_simulation_nointervention_mean_arr = [] +deaths_simulation_mean_arr = [] +deaths_simulation_nointervention_mean_arr = [] +cases_simulation_mean_arr = [] +cases_simulation_nointervention_mean_arr = [] + +for i in range(20): + #%% load simulation data (with and without interventions) + simulation_data = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'NYC_interventions_1e5_%d.txt'%i)) + simulation_data_nointervention = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'NYC_nointerventions_1e5_%d.txt'%i)) + + times = simulation_data[:,0] + times_nointervention = simulation_data_nointervention[:,0] + + kinetic_model = {'S': simulation_data[:,1], 'E': simulation_data[:,2], 'I': simulation_data[:,3], 'H': simulation_data[:,4], 'R': simulation_data[:,5], 'D': simulation_data[:,6]} + kinetic_model_nointervention = {'S': simulation_data_nointervention[:,1], 'E': simulation_data_nointervention[:,2], 'I': simulation_data_nointervention[:,3], 'H': simulation_data_nointervention[:,4], 'R': simulation_data_nointervention[:,5], 'D': simulation_data_nointervention[:,6]} + + #%% determine averages of simulation data + # daily averages of simulation data + # sampling_time = 1 means that we average over 1-day intervals + sampling_time = 7 + daily_average = simulation_average(kinetic_model, times, sampling_time) + cumulative_cases_simulation = (1-np.asarray(daily_average['S'])/population)*1e5 + daily_average = simulation_average(kinetic_model, times, sampling_time) + cumulative_deaths_simulation = np.asarray(daily_average['D'])/population*1e5 + + daily_average_nointervention = simulation_average(kinetic_model_nointervention, times_nointervention, sampling_time) + cumulative_cases_simulation_nointervention = (1-np.asarray(daily_average_nointervention['S'])/population)*1e5 + cumulative_deaths_simulation_nointervention = np.asarray(daily_average_nointervention['D'])/population*1e5 + + cases_simulation = np.ediff1d(cumulative_cases_simulation)/sampling_time + cases_simulation_nointervention = np.ediff1d(cumulative_cases_simulation_nointervention)/sampling_time + + simulation_dates = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time) for i in range(len(cumulative_cases_simulation))]) + simulation_dates_nointervention = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time) for i in range(len(cumulative_deaths_simulation_nointervention))]) + + sampling_time2 = 7 + daily_average2 = simulation_average(kinetic_model, times, sampling_time2) + cumulative_cases_simulation2 = (1-np.asarray(daily_average2['S'])/population)*1e5 + daily_average2 = simulation_average(kinetic_model, times, sampling_time2) + cumulative_deaths_simulation2 = np.asarray(daily_average2['D'])/population*1e5 + daily_average_nointervention2 = simulation_average(kinetic_model_nointervention, times_nointervention, sampling_time2) + cumulative_cases_simulation_nointervention2 = (1-np.asarray(daily_average_nointervention2['S'])/population)*1e5 + cumulative_deaths_simulation_nointervention2 = np.asarray(daily_average_nointervention2['D'])/population*1e5 + + deaths_simulation = np.ediff1d(cumulative_deaths_simulation2)/sampling_time2 + deaths_simulation_nointervention = np.ediff1d(cumulative_deaths_simulation_nointervention2)/sampling_time2 + + simulation_dates2 = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time2) for i in range(len(deaths_simulation))]) + simulation_dates_nointervention2 = np.asarray([NYC_date_of_interest_cases[0] + i*dt.timedelta(days=sampling_time2) for i in range(len(deaths_simulation_nointervention))]) + + #%% plot results + + ax00.plot(simulation_dates+dt.timedelta(days = dd), cumulative_deaths_simulation, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax00.plot(simulation_dates_nointervention+dt.timedelta(days = dd), cumulative_deaths_simulation_nointervention, 'Grey', linewidth = 0.8, alpha = 0.2) + + ax01.plot(simulation_dates+dt.timedelta(days = dd), cumulative_cases_simulation, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax01.plot(simulation_dates_nointervention+dt.timedelta(days = dd), cumulative_cases_simulation_nointervention, 'Grey', linewidth = 0.8, alpha = 0.2) + + ax10.plot(simulation_dates2+dt.timedelta(days = dd), deaths_simulation, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax10.plot(simulation_dates_nointervention2+dt.timedelta(days = dd), deaths_simulation_nointervention, 'Grey', linewidth = 0.8, alpha = 0.2) + + ax11.plot(simulation_dates2+dt.timedelta(days = dd), cases_simulation, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax11.plot(simulation_dates_nointervention2+dt.timedelta(days = dd), cases_simulation_nointervention, 'Grey', linewidth = 0.8, alpha = 0.2) + + cumulative_deaths_simulation_mean_arr.append(cumulative_deaths_simulation) + cumulative_deaths_simulation_nointervention_mean_arr.append(cumulative_deaths_simulation_nointervention) + cumulative_cases_simulation_mean_arr.append(cumulative_cases_simulation) + cumulative_cases_simulation_nointervention_mean_arr.append(cumulative_cases_simulation_nointervention) + deaths_simulation_mean_arr.append(deaths_simulation) + deaths_simulation_nointervention_mean_arr.append(deaths_simulation_nointervention) + cases_simulation_mean_arr.append(cases_simulation) + cases_simulation_nointervention_mean_arr.append(cases_simulation_nointervention) + +ax00.plot(simulation_dates+dt.timedelta(days = dd), np.mean(cumulative_deaths_simulation_mean_arr, axis = 0), 'cornflowerblue', linewidth = 1.5) +ax00.plot(simulation_dates_nointervention+dt.timedelta(days = dd), np.mean(cumulative_deaths_simulation_nointervention_mean_arr, axis = 0), 'Grey', linewidth = 1.5) + +ax01.plot(simulation_dates+dt.timedelta(days = dd), np.mean(cumulative_cases_simulation_mean_arr, axis = 0), 'cornflowerblue', linewidth = 1.5) +ax01.plot(simulation_dates_nointervention+dt.timedelta(days = dd), np.mean(cumulative_cases_simulation_nointervention_mean_arr, axis = 0), 'Grey', linewidth = 1.5) + +ax10.plot(simulation_dates2+dt.timedelta(days = dd), np.mean(deaths_simulation_mean_arr, axis = 0), 'cornflowerblue', linewidth = 1.5) +ax10.plot(simulation_dates_nointervention2+dt.timedelta(days = dd), np.mean(deaths_simulation_nointervention_mean_arr, axis = 0), 'Grey', linewidth = 1.5) + +ax11.plot(simulation_dates2+dt.timedelta(days = dd), np.mean(cases_simulation_mean_arr, axis = 0), 'cornflowerblue', linewidth = 1.5) +ax11.plot(simulation_dates_nointervention2+dt.timedelta(days = dd), np.mean(cases_simulation_nointervention_mean_arr, axis = 0), 'Grey', linewidth = 1.5) + plt.tight_layout() plt.margins(0,0) sns.despine(top=True, right=True, left=True) plt.savefig('new_york_cases.pdf', dpi=300, bbox_inches = 'tight', pad_inches = 0.05) + +# set nice figure sizes +fig_width_pt = 368 # Get this from LaTeX using \showthe\columnwidth +golden_mean = (np.sqrt(5.) - 1.) / 2. # Aesthetic ratio +ratio = golden_mean +inches_per_pt = 1. / 72.27 # Convert pt to inches +fig_width = fig_width_pt * inches_per_pt # width in inches +fig_height = 0.75*fig_width*ratio # height in inches +fig_size = [1.5*fig_width, fig_height] +rcParams.update({'figure.figsize': fig_size}) + +# hospitalization plot +#%% + +fig, axs = plt.subplots(ncols = 2) + +# cumulative death panel +ax00 = axs[0] +ax00_2 = axs[0].twinx() + +ax00.set_zorder(ax00_2.get_zorder()+1) +ax00.patch.set_visible(False) + +ax00.set_ylabel("cumulative") + +ax00.text(dt.date(2020, 3, 12), 0.9*1500, r'(a)') + +ax00.text(dt.date(2020, 5, 10), 1350, r'no SD') + +ax00.set_ylim(0,1500) +ax00.set_yticks([0,500,1000,1500]) +ax00.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax00.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) +ax00.xaxis.set_major_locator(ticker.MultipleLocator(14)) +ax00.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) +ax00.get_yaxis().set_major_formatter( + ticker.FuncFormatter(lambda x, p: format(int(x), ','))) + +ax00_2.set_yticks([]) + +ax00.yaxis.grid() + +# cumulative death panel +ax01 = axs[1] + +ax01_2 = axs[1].twinx() + +ax01.set_zorder(ax01_2.get_zorder()+1) +ax01.patch.set_visible(False) + +ax01.set_ylabel("daily") + +ax01.text(dt.date(2020, 3, 12), 0.9*40, r'(b)') + +ax01.text(dt.date(2020, 5, 2), 35, r'no SD') + +ax01.set_ylim(0,40) +ax01.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 9, 7)]) +ax01.set_xticklabels(NYC_date_of_interest_cases[::7], rotation = 45) +ax01.xaxis.set_major_locator(ticker.MultipleLocator(14)) +ax01.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) +ax01.get_yaxis().set_major_formatter( + ticker.FuncFormatter(lambda x, p: format(int(x), ','))) + +ax01_2.set_yticks([]) + +ax01.yaxis.grid() + +fig.suptitle(r'Hospitalization per 100,000', y = 1.05) + +hospitalizations_simulation_mean_arr = [] +hospitalizations_simulation_nointervention_mean_arr = [] +hospitalizations_simulation_daily_mean_arr = [] +hospitalizations_simulation_nointervention_daily_mean_arr = [] + +for i in range(20): + #%% load simulation data (with and without interventions) + simulation_data = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'NYC_interventions_1e5_%d.txt'%i)) + simulation_data_nointervention = np.loadtxt(os.path.join('..', 'data', 'simulation_data', 'NYC_nointerventions_1e5_%d.txt'%i)) + + times = simulation_data[:,0] + times_nointervention = simulation_data_nointervention[:,0] + + kinetic_model = {'S': simulation_data[:,1], 'E': simulation_data[:,2], 'I': simulation_data[:,3], 'H': simulation_data[:,4], 'R': simulation_data[:,5], 'D': simulation_data[:,6]} + kinetic_model_nointervention = {'S': simulation_data_nointervention[:,1], 'E': simulation_data_nointervention[:,2], 'I': simulation_data_nointervention[:,3], 'H': simulation_data_nointervention[:,4], 'R': simulation_data_nointervention[:,5], 'D': simulation_data_nointervention[:,6]} + + #%% determine averages of simulation data + # daily averages of simulation data + # sampling_time = 1 means that we average over 1-day intervals + + sampling_time = 7 + daily_average = simulation_average(kinetic_model, times, sampling_time) + hospitalizations_simulation = np.cumsum(np.asarray(daily_average['H']))/population*1e5 + + hospitalizations_simulation_daily = np.ediff1d(hospitalizations_simulation)/sampling_time + + daily_average_nointervention = simulation_average(kinetic_model_nointervention, times_nointervention, sampling_time) + hospitalizations_simulation_nointervention = np.cumsum(np.asarray(daily_average_nointervention['H']))/population*1e5 + + hospitalizations_simulation_nointervention_daily = np.ediff1d(hospitalizations_simulation_nointervention)/sampling_time + + ax00.plot(simulation_dates+dt.timedelta(days = dd), hospitalizations_simulation, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax00.plot(simulation_dates_nointervention+dt.timedelta(days = dd), hospitalizations_simulation_nointervention, 'Grey', linewidth = 0.8, alpha = 0.2) + + ax01.plot(simulation_dates[:-1]+dt.timedelta(days = dd), hospitalizations_simulation_daily, 'cornflowerblue', linewidth = 0.8, alpha = 0.2) + ax01.plot(simulation_dates_nointervention[:-1]+dt.timedelta(days = dd), hospitalizations_simulation_nointervention_daily, 'Grey', linewidth = 0.8, alpha = 0.2) + + hospitalizations_simulation_mean_arr.append(hospitalizations_simulation) + hospitalizations_simulation_nointervention_mean_arr.append(hospitalizations_simulation_nointervention) + hospitalizations_simulation_daily_mean_arr.append(hospitalizations_simulation_daily) + hospitalizations_simulation_nointervention_daily_mean_arr.append(hospitalizations_simulation_nointervention_daily) + +ax00.plot(simulation_dates+dt.timedelta(days = dd), np.mean(hospitalizations_simulation_mean_arr, axis = 0), 'cornflowerblue') +ax00.plot(simulation_dates_nointervention+dt.timedelta(days = dd), np.mean(hospitalizations_simulation_nointervention_mean_arr, axis = 0), 'Grey') + +ax01.plot(simulation_dates[:-1]+dt.timedelta(days = dd), np.mean(hospitalizations_simulation_daily_mean_arr, axis = 0), 'cornflowerblue') +ax01.plot(simulation_dates_nointervention[:-1]+dt.timedelta(days = dd), np.mean(hospitalizations_simulation_nointervention_daily_mean_arr, axis = 0), 'Grey') + +fig.tight_layout() +plt.margins(0,0) +sns.despine(top=True, right=True, left=True) +plt.savefig('new_york_hospitalizations.pdf', dpi=300, bbox_inches = 'tight', + pad_inches = 0.05) \ No newline at end of file diff --git a/figs/NYC_plot_nointervention.py b/figs/NYC_plot_nointervention.py deleted file mode 100644 index 5191baa1..00000000 --- a/figs/NYC_plot_nointervention.py +++ /dev/null @@ -1,192 +0,0 @@ -import os, sys; sys.path.append(os.path.join("..")) - -from timeit import default_timer as timer - -import networkx as nx -import numpy as np -import pandas as pd -import random -import datetime as dt -import matplotlib.dates as mdates -import matplotlib.ticker as ticker -import matplotlib.pyplot as plt - -from matplotlib import rcParams - -# customized settings -params = { # 'backend': 'ps', - 'font.family': 'serif', - 'font.serif': 'Latin Modern Roman', - 'font.size': 10, - 'axes.labelsize': 'medium', - 'axes.titlesize': 'medium', - 'legend.fontsize': 'medium', - 'xtick.labelsize': 'small', - 'ytick.labelsize': 'small', - 'savefig.dpi': 150, - 'text.usetex': True} -# tell matplotlib about your params -rcParams.update(params) - -# set nice figure sizes -fig_width_pt = 368 # Get this from LaTeX using \showthe\columnwidth -golden_mean = (np.sqrt(5.) - 1.) / 2. # Aesthetic ratio -ratio = golden_mean -inches_per_pt = 1. / 72.27 # Convert pt to inches -fig_width = fig_width_pt * inches_per_pt # width in inches -fig_height = 0.75*fig_width*ratio # height in inches -fig_size = [fig_width, fig_height] -rcParams.update({'figure.figsize': fig_size}) - -from epiforecast.scenarios import load_edges - -def simulation_average(model_data, times, sampling_time = 1): - """ - Returns daily averages of simulation data. - """ - - simulation_data_average = {} - daily_average = {} - - for key in model_data.keys(): - simulation_data_average[key] = [] - daily_average[key] = [] - - tav = 0 - - for i in range(len(times)): - for key in model_data.keys(): - simulation_data_average[key].append(model_data[key][i]) - - if times[i] >= tav: - for key in model_data.keys(): - daily_average[key].append(np.mean(simulation_data_average[key])) - simulation_data_average[key] = [] - - tav += sampling_time - - return daily_average - -edges = load_edges(os.path.join('..', 'data', 'networks', 'edge_list_SBM_1e4_nobeds.txt')) - -contact_network = nx.Graph() -contact_network.add_edges_from(edges) -contact_network = nx.convert_node_labels_to_integers(contact_network) -population = len(contact_network) - -# plot simulation data and empirical data -simulation_data = np.loadtxt('simulation_data_nointervention.txt') - -times = simulation_data[:,0] - -kinetic_model = {'S': simulation_data[:,1], 'E': simulation_data[:,2], 'I': simulation_data[:,3], 'R': simulation_data[:,4], 'H': simulation_data[:,5], 'D': simulation_data[:,6]} - -NYC_data = pd.read_csv(os.path.join('..', 'data', 'NYC_COVID_CASES', 'data_new_york.csv')) -NYC_cases = np.asarray([float(x) for x in NYC_data['Cases'].tolist()]) -NYC_deaths = np.asarray([float(x) for x in NYC_data['Deaths'].tolist()]) -NYC_date_of_interest = np.asarray([dt.datetime.strptime(x, "%m/%d/%Y") for x in NYC_data['DATE_OF_INTEREST'].tolist()]) - -# population of NYC -NYC_population = 8.399e6 - -# fraction reported cases -fraction_reported = 0.13 - -# cumulative cases NYC -cumulative_reported_cases_NYC = 1/fraction_reported*np.cumsum(NYC_cases)/NYC_population -cumulative_deaths_NYC = np.cumsum(NYC_deaths)/NYC_population*1e5 - -reported_cases_NYC = 1/fraction_reported*NYC_cases/NYC_population -reported_deaths_NYC = NYC_deaths/NYC_population*1e5 - -# daily averages of simulation data -# sampling_time = 1 means that we average over 1-day intervals -sampling_time = 2 -daily_average = simulation_average(kinetic_model, times, sampling_time) -cumulative_cases_simulation = 1-np.asarray(daily_average['S'])/population -cumulative_deaths_simulation = np.asarray(daily_average['D'])/population*1e5 - -reported_deaths_simulation = np.asarray(daily_average['D'])/population*1e5 -reported_deaths_simulation = np.ediff1d(reported_deaths_simulation) - -reported_cases_simulation = np.ediff1d(cumulative_cases_simulation) - -simulation_dates = np.asarray([NYC_date_of_interest[0] + i*dt.timedelta(days=sampling_time) for i in range(len(cumulative_cases_simulation))]) - -fig, ax = plt.subplots() - -ax2 = ax.twinx() - -ax.plot(NYC_date_of_interest[::3], cumulative_reported_cases_NYC[::3], marker = 'o', markersize = 3, color = 'k', ls = 'None', label = r'total cases (NYC)') -ax.plot(simulation_dates+dt.timedelta(days = 17), cumulative_cases_simulation, 'k', label = 'total cases (simulation)') - -ax2.plot(NYC_date_of_interest[::3], cumulative_deaths_NYC[::3], marker = 's', markersize = 4, color = 'darkred', markeredgecolor = 'Grey', ls = 'None', label = r'total deaths (NYC)') -ax2.plot(simulation_dates+dt.timedelta(days = 17), cumulative_deaths_simulation, 'darkred', label = 'total deaths (simulation)') - -#ax.text(dt.date(2020, 3, 10), 0.13, 'no SD') -#ax.text(dt.date(2020, 4, 19), 0.03, 'SD intervention') -#ax.text(dt.date(2020, 6, 21), 0.13, 'loosening SD') -# -# -#ax.fill_between([dt.date(2020, 3, 1), dt.date(2020, 3, 26)], 1, 0, color = 'Salmon', alpha = 0.2) -#ax.fill_between([dt.date(2020, 3, 26), dt.date(2020, 6, 15)], 1, 0, color = 'orange', alpha = 0.2) -#ax.fill_between([dt.date(2020, 6, 15), dt.date(2020, 7, 26)], 1, 0, color = 'Salmon', alpha = 0.2) - -ax.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 7, 26)]) -ax.set_xticklabels(NYC_date_of_interest[::7], rotation = 45) -ax.xaxis.set_major_locator(ticker.MultipleLocator(14)) -ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) - -ax.set_ylim(0,1) -ax.set_ylabel(r'proportion of total cases', labelpad = 3) - -ax2.set_ylim(0,1200) -ax2.set_ylabel(r'total deaths/100,000', color = 'darkred') -ax2.tick_params(axis='y', labelcolor = 'darkred') - -ax.legend(frameon = True, loc = 1, fontsize = 7) -ax2.legend(frameon = True, loc = 4, fontsize = 7) - -plt.tight_layout() -plt.margins(0,0) -plt.savefig('new_york_cases_nointervention.png', dpi=300, bbox_inches = 'tight', - pad_inches = 0.05) - -fig, ax = plt.subplots() - -ax2 = ax.twinx() - -ax.plot(NYC_date_of_interest[::3], reported_cases_NYC[::3], marker = 'o', markersize = 3, color = 'k', ls = 'None', label = r'cases (NYC)') -ax.plot(simulation_dates[:-1]+dt.timedelta(days = 17), reported_cases_simulation, 'k', label = 'cases (simulation)') - -ax2.plot(NYC_date_of_interest[::3], reported_deaths_NYC[::3], marker = 's', markersize = 4, color = 'darkred', markeredgecolor = 'Grey', ls = 'None', label = r'deaths (NYC)') -ax2.plot(simulation_dates[:-1]+dt.timedelta(days = 17), reported_deaths_simulation, 'darkred', label = 'deaths (simulation)') - -#ax.text(dt.date(2020, 3, 10), 0.13, 'no SD') -#ax.text(dt.date(2020, 4, 19), 0.03, 'SD intervention') -#ax.text(dt.date(2020, 6, 21), 0.13, 'loosening SD') -# -# -#ax.fill_between([dt.date(2020, 3, 1), dt.date(2020, 3, 26)], 1, 0, color = 'Salmon', alpha = 0.2) -#ax.fill_between([dt.date(2020, 3, 26), dt.date(2020, 6, 15)], 1, 0, color = 'orange', alpha = 0.2) -#ax.fill_between([dt.date(2020, 6, 15), dt.date(2020, 7, 26)], 1, 0, color = 'Salmon', alpha = 0.2) - -ax.set_xlim([dt.date(2020, 3, 8), dt.date(2020, 7, 26)]) -ax.set_xticklabels(NYC_date_of_interest[::7], rotation = 45) -ax.xaxis.set_major_locator(ticker.MultipleLocator(14)) -ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d')) - -ax.set_ylim(0,0.05) -ax.set_ylabel(r'proportion of cases', labelpad = 3) - -ax2.set_ylim(0,100) -ax2.set_ylabel(r'deaths/100,000', color = 'darkred') -ax2.tick_params(axis='y', labelcolor = 'darkred') - -ax.legend(frameon = True, loc = 1, fontsize = 7) -ax2.legend(frameon = True, loc = 4, fontsize = 7) - -plt.tight_layout() -plt.margins(0,0) -plt.savefig('new_york_cases2_nointervention.png', dpi=300, bbox_inches = 'tight', - pad_inches = 0.05)