<jupyter_start><jupyter_text>Variying parametersThis notebook contains code to generate some of the data in the paper "Optimizing testing policies for detecting COVID-19 within institutions"by <NAME>, <NAME>, <NAME>, and <NAME>.  Load dependencies<jupyter_code>import networkx
import numpy as np
import pandas as pd
import os
from models import *
from helpers import *<jupyter_output><empty_output><jupyter_text>Parameter setup<jupyter_code>G_normal , G_quarantine = make_graphs(numNodes = 1000, m=9, scale = 100 , plot_degree = True)

base = {
    "G": {"numNodes" : 1000, "m" : 9, "scale":100 , "plot_degree": False}, #Network adjacency matrix (numpy array) or Networkx graph object.
    "beta":0.2,# Rate of transmission (exposure) (global)
    "sigma":1/5.2, # Rate of infection (upon exposure)
    "gamma" : 1/14, # Rate of recovery (upon infection)
    "mu_I": 0.0004, # Rate of infection-related death
    "mu_0": 0,# Rate of baseline death
    "nu":0, # Rate of baseline birth
     "xi":0, # Rate of re-susceptibility (upon recovery)
     "p":0.5, # Probability of interaction outside adjacent nodes
     "Q": None, # Quarantine adjacency matrix (numpy array) or Networkx graph object.
     "beta_D": 0, # Rate of transmission (exposure) for individuals with detected infections (global), changed from 0.155 to set so detected individuals can't transmit
     "sigma_D": 1/5.2,   # Rate of infection (upon exposure) for individuals with detected infections
     "gamma_D":1/14, # Rate of recovery (upon infection) for individuals with detected infections
      "mu_D":0.0004, # Rate of infection-related death for individuals with detected infections
      "theta_E":0, # probability that exposed person is tested (i.e., becomes symptomatic)
      "theta_I":0, # probability that infected person is tested (i.e., becomes symptomatic)
      "phi_E":0, # each detected contact adds phi_E chance for probability of infection for exposed individual
      "phi_I": 0, # each detected contact adds phi_I chance for probability of infection for infected individual
      "psi_E": 1.0, #Probability of positive test results for exposed individuals. 1 =  no false negatives
      "psi_I": 1.0, # Probability of positive test results for exposed individuals. 1 = no false negatives
       "q":0.0, # Probability of quarantined individuals interaction outside adjacent nodes.   changed to 0 to make it strong quarantine
       "initI": 0, # Init number of infectious individuals
       "initE":0, # Init number of exposed  individuals
       "initD_E":0, # Init number of detected exposed individuals
       "initD_I":0, # Init number of detected infectious individuals
       "initR":0, # Init number of recovered individuals
      "initF":0, # Init number of infection -related fatalities
      "p_extern": 1/1000, # Rate of spontanous (external) infection
      "p_periodic": 0, # Fraction of people tested in a periodic subsample tests
      "period": 1, # Period for periodic testing
      "batch" : True, # Test people in batch (True) or test each person with probability p_periodic (False)
      "min_time" : 0, # Minimum time to pass between two tests for the same individual
      "store_Xseries" : False, # Store the full data per node
      "count_non_random" : False, #  True: count tests done to symptomatic or traced people
      "T": 84, # time to run the execution
      "verbose": False, # print log messages during run
       "checkpoints": None, # checkpoints (see SIERPlus documentation - we don't use this)
       "stopping": "1st", # function for stopping policy '1st' means stop at first detection
       "policy": None, # policy function for adaptive policies
       "policyInterval": 1, # period to apply policy function
        "type": "[UNKNOWN]", # type of simulation - useful for aggragating data
        "runTillEnd": True, # True - don't stop even if zero infections (makes sense when external infections > 0),
        "test_recovered": True, # Test people in "recovered" state - makes sense when we stop at first detection, not so much otherwise
        "initT": 0 # Initial time step
}

orig_base = dict(base)
# Testing symptomatic and tracing (unused when we stop at first detection)

symp_ratio = 0.2 # fraction of people that develop symptoms significant enough to be tested and isolated sometime during the time they spend in infected or exposed state
# this number is necessarily guessed, since it depends heavily on the demographics of the workplace
# as well as policies in place for self isolation and testing.
#symp_affin  =symp_ratio* base["sigma"] * base["gamma"] / ( base["gamma"] + base["sigma"] )

symp_affin = 0.2 * base["gamma"]

# Chosen so that symp_affin/sigma + symp_affin/sigma = symp_ratio

# probability of developing infection

symptomatic = {
    "theta_E" : 0 , # probability that exposed person is tested (i.e., becomes symptomatic)
    "theta_I" : symp_affin  # probability that infected person is tested (i.e., becomes symptomatic)
}


# a/sigma = 0.2*(1/sigma / (1/gamma + 1/delta)) = gamma*sigma/

# Contact tracing - a detected E/I neighbor increases the chance to be tested by 50%
tracing =  {
    "phi_E": 0.8,
    "phi_I": 0.8
}

# Uncomment lines below to include tracing & testing of symptomatic people

#base.update(tracing)
#base.update(symptomatic)<jupyter_output><empty_output><jupyter_text>We set a false negative rate of 67% for exposed individuals and 20% for infected individuals.<jupyter_code># False negatives

base.update({
    "psi_E": 0.33, #Probability of positive test results for exposed individuals. 1 =  no false negatives
    "psi_I": 0.8, # Probability of positive test results for exposed individuals. 1 = no false negatives
})<jupyter_output><empty_output><jupyter_text>Experiments<jupyter_code>def test(days,batches):
    """Modify parameters to test 1/batches of the population every days/batches days"""
    return { "p_periodic" : 1/batches, "period": days/batches }

def product(A,B):
    return [(a,b) for a in A for b in B]

to_do =product([28],[1,2,3,28]) + \
       product([14],[1,2,3,14])


    

base["T"] = 108
to_do = [(28,1),(28,4)]

settings = [{**base, **test(d,k), "type":f"{d}/{k}" } for (d,k) in to_do ]



settings += [{**base , "type": "No testing (56 days)", "T": 56 } ]


base<jupyter_output><empty_output><jupyter_text>Get variants for p_extern, beta<jupyter_code>from fractions import Fraction

std_graph =   {"numNodes" : 1000, "m" : 9, "scale":100 , "plot_degree": False}

options = [(beta,p,std_graph) for beta in np.linspace(0.025,0.3,15) for p in np.linspace(0.00005,0.001,15)]


def name(G):
    if G["m"] <=9:
        return "" 
    else: 
        return " (more edges)"

to_run = [] 
variants = {}
for beta,p,G in options:
    #print(beta,p,G)
    q = Fraction(p).limit_denominator(100000) if (p and (not isinstance(p,str))) else 0
    for d in settings:
        e = dict(d)
        
        v = f"beta={beta},pextern={q}"+name(G)
        if p==0:
            e["initI"] = 1
            v += " (1 initial infection)"
        if isinstance(p,str) and p[0]=='r':
            a = float(p[1:])
            e["initT"]=(0,a)
            e["initI"]=1
            e["p_extern"] = 0
            v += f" (random start 0-{a})"
        else:
            e["p_extern"] = p
        if e["type"][:15] != "Business closed":
            e["G"] = G
            e["beta"] = beta
            
        #e["verbose"] = True
        variants[v] = { "beta":beta, "p_extern": p, "G":G }
        e["variant"] =  v
        to_run += [e]
        
print(len(variants.keys()))
#variants<jupyter_output>225<jupyter_text>This is the time consuming part<jupyter_code># commented out - this is just when debugging
#to_run = [{**base, "p_periodic":1/20, "period":1, "T":200, "p_extern":1/30000 , "type":"Trial", "verbose":False, "min_time":0}]
realizations = 400
#lview = None 
len(to_run)*realizations
if not os.path.exists("data/"):
    os.makedirs("data")

timestamp = datetime.now().strftime('%Y%m%d_%H_%M_%S')
datadir = f"data/{timestamp}_{realizations}"
os.makedirs(datadir)
print(datadir)
def pline(L,n=None):
    T = []
    for a in L:
        if isinstance(a,float):
            if a > 0.001:
                T += [f"{a:.3f}"]
            else:
                q = Fraction(p).limit_denominator(100000) 
                T += [f"{q}"]
        else:
            T += [str(a) ]
    if not n:
        n = max([len(a)+2 for a in T])
    
    return "|".join([a.ljust(n) for a in T])+"\n"
                
        
    
readme = fr"""
# Description

The above is results of simulations done at {timestamp}

We use {len(settings)} parameter settings, each with {len(list(variants.keys()))} variants. We describe the main settings below:

We use {realizations} realizations, for a total of {realizations * len(to_run)} executions.

## Variants:

"""
for v in variants.keys():
    readme += f"   *  {v}\n"

readme += "## Settings \n"

keys = ["type", "p_periodic", "period"]
readme += pline(keys,25)
readme += "|".join(["-"*25]*len(keys))+"\n"
for s in settings:
    readme += pline([s["type"], s["p_periodic"],s["period"] ],25)

readme += f"""

## Baseline parameters
"""    

readme += pline(["parameter","value"],20)
readme += "|".join(["-"*20]*2)+"\n"
for k in base:
    readme += pline([k,base[k]],20)

with open(f"{datadir}/readme.md", 'a') as out:
    out.write(readme + '\n')

from IPython.display import display, Markdown

display(Markdown(readme))
import pickle
with open(f'{datadir}/to_run.pickle', 'wb') as handle:
    pickle.dump(to_run, handle, protocol=pickle.HIGHEST_PROTOCOL)
import os
os.chdir(datadir)
!pwd
%%time
!python ../../par_run.py $realizations
os.chdir("/home/jupyter/covidtesting/")<jupyter_output><empty_output><jupyter_text>**End of time consuming part**<jupyter_code>datadir = 'data/20200623_12_17_42_400' #'data/20200622_21_02_55_400'
datas = []
for i in range(3):
    datas.append(pd.read_pickle(f'{datadir}/data_{i}.zip'))
data = pd.concat(datas)
data

data["R"] = data["beta"]*14
data.loc[:,"risk"] = data["meanUndetected1st"]
data["Baseline"] = data["p_extern"]*14*1000

infrequent = data[data.type == "28/1"].groupby(["p_extern","R"]).risk.mean().unstack()
infrequent
frequent = data[data.type == "28/4"].groupby(["p_extern","R"]).risk.mean().unstack()
notesting = data[data.type == "No testing (56 days)"].groupby(["p_extern","R"]).risk.mean().unstack()
notesting

baseline=  data[data.type == "28/4"].groupby(["p_extern","R"]).Baseline.mean().unstack()
panel1 = 100*frequent/baseline
panel1
panel2 = 100*frequent/infrequent
panel3 = 100*frequent/notesting
panel3
import seaborn as sns
L1 = sns.color_palette("RdBu_r", 100).as_hex()
greens = list(reversed(sns.color_palette("BuGn_d",50).as_hex()))
L = greens+L1[:50]
for i in range(50,100):
    L+= [L1[i]]*4
from matplotlib.colors import ListedColormap
cmap = ListedColormap(L)
cmap
from matplotlib.ticker import StrMethodFormatter 
size = 25
titlesize = 30

fig , axes = plt.subplots(1,3, sharex=True, sharey= True, figsize= (20,10))
cbar_ax = fig.add_axes([1.01, 0.2, .03, .7])
axes[0].set_title("a) 28/4 vs. baseline",loc="left",fontsize=int(titlesize))
#axes[0].title.set_text('28/4 vs baseline')
#cmap = sns.diverging_palette(150, 275, sep=1,s=80, l=55, n=25)



color = "white"
sns.heatmap(panel1,ax=axes[0],cbar=False, vmin=0, annot=False, vmax=300, cmap=cmap)# annot=True,fmt='.0f' ,cmap=cmap) # annot_kws={'color':color}

axes[1].set_title('b) 28/4 vs 28/1', loc='left',fontsize=int(titlesize))
sns.heatmap(panel2,ax=axes[1],cbar=False, vmin=0, vmax=300, cmap=cmap)# annot=True,fmt='.0f', cmap=cmap)
axes[2].set_title('c) 28/4 vs No Testing', loc='left',fontsize=int(titlesize))
sns.heatmap(panel3,ax=axes[2],cbar=True, vmin=0, vmax=300,  cbar_ax = cbar_ax ,cmap=cmap) # annot=True,fmt='.0f', cmap=cmap)

cbar_ax.tick_params(labelsize=20)
def f(x, pos):
    if pos % 2: return ""
    p= list(np.linspace(0.00005,0.001,15))[pos]
    return f"{p:.5f}"


def g(x,pos):
    if pos % 2: return ""
    beta = list(np.linspace(0.025,0.3,15))[pos]
    return f"{beta*14:.2f}"
    
import matplotlib 

f_ = matplotlib.ticker.FuncFormatter(f)
g_ = matplotlib.ticker.FuncFormatter(g)

for i in range(3):
    axes[i].set_xlabel('' if i!=1 else 'Internal repdroductive number / R', fontsize=size)
    axes[i].set_ylabel('' if i else 'Probability of infection', fontsize=size)
    axes[i].tick_params(axis='x', labelsize=size ) 
    axes[i].tick_params(axis='y', labelsize=size ) 
    axes[i].yaxis.set_major_formatter(f_)
    axes[i].xaxis.set_major_formatter(g_)

cbar = axes[2].collections[0].colorbar
cbar_ax.set_yticklabels([f"{int(i)}%" for i in cbar.get_ticks()]) # set ticks of your format

    
fig.tight_layout() #rect=[0, 0, .9, 1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1473519027069439, 'fraction_numerical': 0.0335033346410357, 'mean_word_length': 3.8947772657450077, 'pattern_counts': {'":': 62, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '457845', 'license': 'MIT', 'n_tokens_mistral': 4754, 'n_tokens_neox': 4349, 'n_words': 1237, 'path': 'testing_policies/generate_continuous_params.ipynb', 'repo_name': 'boazbk/seirsplus-old', 'size': '12773'}