from scipy.stats import beta import matplotlib.pyplot as plt import math import numpy as np import random import seaborn as sns import pandas as pd

sns.set()

master_dict = {} master_dict = {‘test_period’: ‘Test Period’,’test_S’: ‘Susceptible Population’, ‘test_E’: ‘Exposed Population’, ‘test_NE’: ‘Newly Exposed’, ‘test_I’: ‘Infectious Population’, ‘test_NI’: ‘Newly infectious’, ‘test_NI_C’: ‘Total infectious’, ‘test_IC’: ‘Daily confirmed infections’, ‘test_IC_H’: ‘Daily hospitalized’, ‘test_IC_No_H’: ‘Daily Confirmed, Not Hospitalized’, ‘test_IU’: ‘Daily Uncomfirmed Infected’, ‘test_IC_C’: ‘Cumulative confirmed infections’, ‘test_R’: ‘Recovered Population’, ‘test_NR’: ‘Newly Recovered’, ‘test_D’: ‘Dead Population’, ‘test_ND’: ‘Newly Dead’, ‘test_D_C’: ‘Cumulative Dead’}

#conf_perc = 444 / 10_000_000 #print(conf_perc) print(master_dict.get(‘test_I’))

N = 10_000 # population size S = 0 # susceptible E = 0 # exposed I = 0 # infectious H = C_H = 0 # hospitalizations, confirmed and hospitalized C = 250 # confirmed case C_NH = 0 # confirmed not in hospital U = 0 # unconfirmed case R = 0 # recovering D = 0 # dead S_R = 0 # recovered but susceptible again

C = H + C_NH I = C + U

h2i = .0045 # H1N1 ratio for hospitalized to infected h2i_max = .012 h2i_min = .0016

c2i = 1 / 79 # H1N1 ratio for confirmed cases to estimated total infected c2i_min = 1 / 148 c2i_max = 1 / 47

person-to-person spread: primary means of propagation https://www.cdc.gov/coronavirus/2019-ncov/prepare/transmission.html?CDC_AA_refVal=https%3A%2F%2Fwww.cdc.gov%2Fcoronavirus%2F2019-ncov%2Fabout%2Ftransmission.html

rep_contact = 4 rep_contact_std = 2 new_contact = 2 new_contact_std = 1

E_dt_mean = 5.1 # length of incubation period E_dt_max = 11.5 # max CI value of 15.6 E_dt_min = 2.5 # min CI value of 1.9

i_mild = .51 # mild infection rate i_mod = .388 # mod inf rate i_severe = .052 # severe inf rate i_critical = .006 # critical rate i_asym = .044 # asymptomatic rate

spread mostly likely during infectious period, unlikely once recovery begins: https://www.statnews.com/2020/03/09/people-shed-high-levels-of-coronavirus-study-finds-but-most-are-likely-not-infectious-after-recovery-begins/

i_rec_mean = 7 # infectious period, time to begining recovery after exposure to infectious transition i_rec_max = 12 i_rec_min = 4

f_rate_mean = .015 # (.009 default) case fatality rate, https://wwwnc.cdc.gov/eid/article/26/6/20-0320_article f_rate_min = .0025 # max CFR (.0025) f_rate_max = .03 # max CFR (.03 default)

C = 395 quar_num = 500 # number of confirmed cases to trigger quarantine

def dangerPerception(dp_c, dp_q): #print(dp_c/dp_q) if (dp_c/dp_q) <= 0.05: return .05 elif (dp_c/dp_q) <= 0.95: return dp_c / dp_q else: return 0.95

#dang_perc = dangerPerception(C, quar_num) #print(dang_perc)

def behaviorModifier(bm_dp): beta_a = 10 - (10 * bm_dp) beta_b = 10 - beta_a #print(beta_a, beta_b) return np.random.beta(beta_a, beta_b)

#beh_mod = behaviorModifier(dang_perc) #print(beh_mod)

def calc_stdv(cs_mean, cs_min, cs_max): cs_avg = ((cs_mean - cs_min) + (cs_max - cs_mean)) / 2 return cs_avg / 1.65

i_rec_stdv = calc_stdv(i_rec_mean, i_rec_min, i_rec_max) #print(i_rec_stdv)

e_dt_stdv = calc_stdv(E_dt_mean, E_dt_min, E_dt_max) #print(e_dt_stdv)

def offset_normal(m, a, b, r): a_sd = (m - a) / 1.69 b_sd = (b - m) / 1.69 o_n = [] for x in range(r): i = np.random.uniform(a, b) if i < m: av = np.random.normal(m, a_sd) if av > m: av = m - (av - m) if av < 0: av = 0 o_n.append(av) if i >= m: bv = np.random.normal(m, b_sd) if bv < m: bv = m + (m - bv) o_n.append(bv) return o_n

def norm_norm(m, s, r): n_n = [] for x in range(r): i = np.random.normal(m, s) if i < 0: i = 0 n_n.append(i) return n_n

infectious_period_offset = offset_normal(i_rec_mean, i_rec_min, i_rec_max, N)

plt.hist(infectious_period_offset, 50)

plt.show()

incubation_period_offset = offset_normal(E_dt_mean, E_dt_min, E_dt_max, N)

plt.hist(incubation_period_offset, 50)

plt.show()

h2i_offset = offset_normal(h2i, h2i_min, h2i_max, N)

plt.hist(h2i_offset, 50)

plt.show()

r_contact_norm = norm_norm(rep_contact, rep_contact_std, N)

plt.hist(r_contact_norm, 50)

plt.show()

n_contact_norm = norm_norm(new_contact, new_contact_std, N)

plt.hist(n_contact_norm, 50)

plt.show()

c2i_offset = offset_normal(c2i, c2i_min, c2i_max, N)

plt.hist(c2i_offset, 50)

plt.show()

f_rate_offset = offset_normal(f_rate_mean, f_rate_min, f_rate_max, N)

Running the main ensemble test here:

quar_num_test = 10 big_S = []

for num in range(1_000):

test_period = 180 # days in test period
test_S = [N] # total susceptible
test_E = [1] # total exposed
test_NE = [0] # new exposed
test_I = [0] # total infectious
test_NI = [0] # new infectious
test_NI_C = [0] # total infectious
test_IC = [0] # confirmed infections
test_IC_H, test_IC_No_H, test_IU = [0], [0], [0]
test_IC_C = [0] # cumulative confirmed infections
test_R = [0] # total recovered
test_NR = [0] # new recovered
test_D = [0] # total dead
test_ND = [0] # new dead
test_D_C = [0] # cumulative dead

inf_per = random.sample(infectious_period_offset, 1)[0]
if inf_per < 1:
    inf_per = 1
inc_per = random.sample(incubation_period_offset, 1)[0]
if inc_per < 1:
    inc_per = 1
h2i_test = random.sample(h2i_offset, 1)[0]
c2i_test = random.sample(c2i_offset, 1)[0]
r_cont_test = random.sample(r_contact_norm, 1)[0]
n_cont_test = random.sample(n_contact_norm, 1)[0]
f_rate_test = random.sample(f_rate_offset, 1)[0]

rec_rate = 1 / inf_per # recovery rate, conversion rate from infectious to recovering
lat_rate = 1 / inc_per # latent rate, conversion rate from exposed to infectious
f_rate_rate = f_rate_test / inf_per # case fatality rate *per day* of infectiousness

for day in range(test_period):
    
    today = day + 1
    yesterday = day

    cur_sus = test_S[yesterday] # today's number of currently susceptible - the population that has not been exposed yet
    confirmed_cases = test_IC_C[yesterday] # confirmed cases based on the the H1N1 study proportions of infectious to confirmed

    dang_per_test = dangerPerception(confirmed_cases, quar_num_test) # testing for danger perception given the number of confirmed cases and max risk or quarantine number
    beh_mod_test = behaviorModifier(dang_per_test) # behavioral modification multiplier given the danger/risk perception

    new_NE = (beh_mod_test * r_cont_test * test_NI[yesterday]) + (beh_mod_test * n_cont_test * test_NI[yesterday]) # personal/close contacts and exposures for the newly infectious
    
    not_new_I = test_I[yesterday] - test_NI[yesterday] - test_IC_H[yesterday]
    if not_new_I < 0:
        not_new_I = 0
    not_new_NE = beh_mod_test * n_cont_test * not_new_I # continued personal/close contacts with non-repeat persons for the infectious at date of infectious > 1
    
    next_NE = new_NE + not_new_NE # new exposed population today
    if next_NE >= cur_sus:
        next_NE = test_S[yesterday]

    next_NI = test_E[yesterday] * lat_rate # new infectious population today
    next_NI_C = test_NI_C[yesterday] + next_NI # cumulative infectous population

    next_E = test_E[yesterday] + next_NE - next_NI # total exposed -> adds new exposures, loses new infectious
    if next_E < 0:
        next_E = 0
    
    next_IC = c2i_test * test_I[yesterday] # current confirmed infectious based on H1N1 ratios and estimated total infectious
    new_IC = next_IC - test_IC[yesterday]
    if new_IC < 0:
        new_IC = 0

    next_IC_C = test_IC_C[yesterday] + new_IC # cumulative confirmed infectious
    
    next_IC_H = h2i_test * test_I[yesterday] # current hospitalized infectious based on H1N1 ratios and estimated total infectious

    next_NR = test_I[yesterday] * rec_rate # new recoveries today
    next_R = test_R[yesterday] + next_NR # cumulative recovery rate

    next_D = test_IC[yesterday] * f_rate_rate
    next_ND = next_D - test_D[yesterday]
    if next_ND < 0:
        next_ND = 0

    next_D_C = test_D_C[yesterday] + next_D

    next_I = test_I[yesterday] + next_NI - next_NR - next_ND
    if next_I < 0:
        next_I = 0
    
    next_S = N - next_E - next_I - next_R - next_D
    if next_S < 0:
        next_S = 0

    test_S.append(next_S)
    test_E.append(next_E)
    test_NI.append(next_NI)
    test_NI_C.append(next_NI_C)
    test_I.append(next_I)
    test_NR.append(next_NR)
    test_R.append(next_R)
    test_ND.append(next_ND)
    test_D.append(next_D)
    test_D_C.append(next_D_C)
    test_IC.append(next_IC)
    test_IC_H.append(next_IC_H)
    test_IC_C.append(next_IC_C)

# if len(big_S) > 0:
#     big_S = np.vstack((big_S, test_IC))
# else:
#     big_S = test_IC

test_num = num
test_name = "Test #" + str(num)
test_1 = test_I
test_2 = test_D
try: df
except NameError: df = None
if df is None:
    df = pd.DataFrame({test_name: test_1})
else:
    df_new = pd.DataFrame({test_name: test_1})
    df = pd.concat([df, df_new], axis=1)

try: df2
except NameError: df2 = None
if df2 is None:
    df2 = pd.DataFrame({test_name: test_2})
else:
    df_new2 = pd.DataFrame({test_name: test_2})
    df2 = pd.concat([df2, df_new2], axis=1)

big_array = np.stack([test_S, test_E, test_I, test_IC, test_IC_H, test_R, test_D]).T

plt.plot(big_array)

plt.legend([‘Sus’, ‘Exp’, ‘Inf’, ‘Conf’, ‘Hosp’, ‘Rec’, ‘Dead’])

plt.xlabel(‘Days’)

plt.show()

lombard_conf_new = [0, 0, 15, 40, 57, 61,67,65,98,128,84,369,270,266,300,431,361,808,769,1280,322,1489,1445,1095,1865,1587,1377,1571,1493,2171, 2380, 3251, 1691, 1555, 1942, 1643, 2543, 2409] hubei_conf = [0, 0, 0, 0, 0, 444,444,549,761,1058,1423,3554,3554,4903,5806,7153,11177,13522,16678,19665,22112,24953,27100,29631,31728,33366,33366,48206,54406,56249,58182,59989,61682,62031,62442,62662,64084,64084,64287,64786,65187,65596,65914,66337,66907,67103,67217,67332,67466,67592,67666,67707,67743,67760,67773,67781,67786,67790,67794,67798,67799,67800,67800] ny_conf = [1, 1,2,2,4,5,12,14,20,37,52,96,155,269,330,464,645,1339,2468,4408,6211,9045,12305,14905,20011,23112,25399]

hubei_fat = [0, 0, 0, 0, 0, 0, 0, 0, 0, 17,17,24,40,52,76,125,125,162,204,249,350,414,479,549,618,699,780,871,974,1068,1068,1310,1457,1596,1696,1789,1921,2029,2144,2144,2346,2346,2495,2563,2615,2641,2682,2727,2761,2803,2835,2871,2902,2931,2959,2986,3008,3024,3046,3056,3062,3075,3085,3099,3111,3122,3130] lombard_fat = [0, 0, 0, 1, 1, 4, 2, 1, 5, 3, 6, 1, 14, 17, 18, 25, 37, 19, 113, 66, 135, 149, 127, 146, 76, 252, 202, 220, 319, 209, 381, 546, 361, 320, 402, 296, 387, 541] ny_deaths = [0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,1,5,7,10,20,22,43,60,99,131,192,280,365,450]

lombard_conf = np.cumsum(lombard_conf_new) # Lombardi region only (pop 10m) - population density probably doesn’t work for calc entire countries lombard_cum = np.cumsum(lombard_fat) italy_norm = [i/1000 for i in lombard_conf] # Lombardi region only hubei_norm = [i/1100 for i in hubei_conf] ny_norm = [i/1880 for i in ny_conf]

hubei_fat_norm = [i/1100 for i in hubei_fat] lombard_fat_norm = [i/1000 for i in lombard_cum] ny_deaths_norm = [i/1880 for i in ny_deaths]

print(hubei_fat_norm)

df.plot(alpha=0.3, legend=False).set(title = f’COVID-19 per 10k Population - 1k Model Runs’, xlabel=’Day’, ylabel=’Per 10k Population’) plt.show()

df_mean = df.mean(axis = 1) df_10 = df.quantile(.1, axis=1) df_90 = df.quantile(.9, axis=1)

df2_mean = df2.mean(axis = 1)

df_mean.plot() df_10.plot() df_90.plot() plt.show()

sns.lineplot(data=df_mean, label=’Model 1’).set(title = ‘COVID-19 Models per 10,000 Simulated’, xlabel = ‘Day’, ylabel = ‘Per 10,000’) df_90.plot(label=’90 percentile’) df_10.plot(label=’10 percentile’)

plt.plot(hubei_norm, label=’Hubei’)

plt.plot(lombard_fat_norm, label=’Lombardy (ITA)’)

plt.plot(ny_deaths_norm, label=’NYC Confirmed Cases’) #sns.lineplot(data=df2_mean, label=’Model 2’)

plt.plot(ny_norm, label=’New York’)

plt.plot(italy_norm, label=’Lombardi (ITA)’)

legend = plt.legend()

legend.get_texts()[0].set_text(‘Hubei’)

plt.legend() plt.show()

plt.plot(big_S_stack.T)

plt.show()

df_mean.to_csv(‘test_df_1.csv’)