From c988d1a935bdb499115c4dcd550b7af7ac17522d Mon Sep 17 00:00:00 2001 From: georgia witchel Date: Thu, 29 Dec 2022 20:24:42 -0700 Subject: [PATCH 1/2] added SVI as a factor to consider within the population --- covasim/defaults.py | 2 ++ covasim/parameters.py | 6 ++++-- covasim/people.py | 3 +++ covasim/population.py | 8 +++++--- examples/run_sim.py | 2 +- 5 files changed, 15 insertions(+), 6 deletions(-) diff --git a/covasim/defaults.py b/covasim/defaults.py index 291fe9920..707c26cd4 100644 --- a/covasim/defaults.py +++ b/covasim/defaults.py @@ -44,6 +44,7 @@ def __init__(self): 'uid', # Int 'age', # Float 'sex', # Float + 'SVI', # Int 'symp_prob', # Float 'severe_prob', # Float 'crit_prob', # Float @@ -219,6 +220,7 @@ def __init__(self): 'severe', ] +default_SVI_data = np.array([]) # Default age data, based on Seattle 2018 census data -- used in population.py default_age_data = np.array([ [ 0, 4, 0.0605], diff --git a/covasim/parameters.py b/covasim/parameters.py index 616c8c9f8..bf6463034 100644 --- a/covasim/parameters.py +++ b/covasim/parameters.py @@ -227,7 +227,7 @@ def reset_layer_pars(pars, layer_keys=None, force=False): return -def get_prognoses(by_age=True, version=None): +def get_prognoses(by_age=True,by_SVI=True, version=None): ''' Return the default parameter values for prognoses @@ -251,7 +251,7 @@ def get_prognoses(by_age=True, version=None): crit_probs = np.array([0.04]), death_probs = np.array([0.01]), ) - else: + else : prognoses = dict( age_cutoffs = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90,]), # Age cutoffs (lower limits) sus_ORs = np.array([0.34, 0.67, 1.00, 1.00, 1.00, 1.00, 1.24, 1.47, 1.47, 1.47]), # Odds ratios for relative susceptibility -- from Zhang et al., https://science.sciencemag.org/content/early/2020/05/04/science.abb8001; 10-20 and 60-70 bins are the average across the ORs @@ -262,6 +262,8 @@ def get_prognoses(by_age=True, version=None): crit_probs = np.array([0.00003, 0.00008, 0.00036, 0.00104, 0.00216, 0.00933, 0.03639, 0.08923, 0.17420, 0.17420]), # Overall probability of developing critical symptoms (derived from Table 1 of https://www.imperial.ac.uk/media/imperial-college/medicine/mrc-gida/2020-03-16-COVID19-Report-9.pdf) death_probs = np.array([0.00002, 0.00002, 0.00010, 0.00032, 0.00098, 0.00265, 0.00766, 0.02439, 0.08292, 0.16190]), # Overall probability of dying -- from O'Driscoll et al., https://www.nature.com/articles/s41586-020-2918-0; last data point from Brazeau et al., https://www.imperial.ac.uk/mrc-global-infectious-disease-analysis/covid-19/report-34-ifr/ ) + + prognoses = relative_prognoses(prognoses) # Convert to conditional probabilities # If version is specified, load old parameters diff --git a/covasim/people.py b/covasim/people.py index 2c45b1faa..e7c5bc0bf 100644 --- a/covasim/people.py +++ b/covasim/people.py @@ -158,6 +158,9 @@ def set_prognoses(self): self.rel_sus[:] = progs['sus_ORs'][inds] # Default susceptibilities self.rel_trans[:] = progs['trans_ORs'][inds] * cvu.sample(**self.pars['beta_dist'], size=len(inds)) # Default transmissibilities, with viral load drawn from a distribution + # Consider SVI + self.symp_prob = self.symp_prob + self.symp_prob*0.65*self.SVI # source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7318979/ + return diff --git a/covasim/population.py b/covasim/population.py index effd89700..65eb534bd 100644 --- a/covasim/population.py +++ b/covasim/population.py @@ -93,7 +93,7 @@ def make_people(sim, popdict=None, die=True, reset=False, recreate=False, verbos people = popdict people.set_pars(sim.pars) else: - people = cvppl.People(sim.pars, uid=popdict['uid'], age=popdict['age'], sex=popdict['sex'], contacts=popdict['contacts']) # List for storing the people + people = cvppl.People(sim.pars, uid=popdict['uid'], age=popdict['age'], sex=popdict['sex'], SVI=popdict['SVI'], contacts=popdict['contacts']) # List for storing the people sc.printv(f'Created {pop_size} people, average age {people.age.mean():0.2f} years', 2, verbose) @@ -140,7 +140,7 @@ def validate_popdict(popdict, pars, verbose=True): return -def make_randpop(pars, use_age_data=True, use_household_data=True, sex_ratio=0.5, microstructure='random', **kwargs): +def make_randpop(pars, use_age_data=True, use_household_data=True, sex_ratio=0.5, SVI_ratio=0.1, microstructure='random', **kwargs): ''' Make a random population, with contacts. @@ -192,7 +192,7 @@ def make_randpop(pars, use_age_data=True, use_household_data=True, sex_ratio=0.5 warnmsg = f'Could not load household size data for requested location "{location}" ({str(E)}), using default' cvm.warn(warnmsg) - # Handle sexes and ages + # Handle sexes, ages, and SVI's uids = np.arange(pop_size, dtype=cvd.default_int) sexes = np.random.binomial(1, sex_ratio, pop_size) age_data_min = age_data[:,0] @@ -202,12 +202,14 @@ def make_randpop(pars, use_age_data=True, use_household_data=True, sex_ratio=0.5 age_data_prob /= age_data_prob.sum() # Ensure it sums to 1 age_bins = cvu.n_multinomial(age_data_prob, pop_size) # Choose age bins ages = age_data_min[age_bins] + age_data_range[age_bins]*np.random.random(pop_size) # Uniformly distribute within this age bin + SVIs = np.random.binomial(1, SVI_ratio, pop_size) # SVI's # Store output popdict = {} popdict['uid'] = uids popdict['age'] = ages popdict['sex'] = sexes + popdict['SVI'] = SVIs # Actually create the contacts if microstructure == 'random': diff --git a/examples/run_sim.py b/examples/run_sim.py index 52c95d90f..058b9a865 100644 --- a/examples/run_sim.py +++ b/examples/run_sim.py @@ -19,7 +19,7 @@ pop_type = 'hybrid', # Population to use -- "hybrid" is random with household, school,and work structure ) -# Optionally add an intervention +#Optionally add an intervention if interv: pars.interventions = cv.change_beta(days=45, changes=0.5) From b95f956ba6568019f75172f2f9ea56e7d285ebb6 Mon Sep 17 00:00:00 2001 From: georgia witchel Date: Mon, 27 Feb 2023 13:22:19 -0800 Subject: [PATCH 2/2] greenspace intervention --- beta_changes.csv | 9 +++ covasim/defaults.py | 1 - covasim/interventions.py | 84 +++++++++++++++++++++++++ covasim/parameters.py | 2 +- examples/example_data.csv | 33 ---------- examples/introduce_greenspace.py | 104 +++++++++++++++++++++++++++++++ examples/run_sim_change_beta.py | 81 ++++++++++++++++++++++++ examples/run_sim_change_trans.py | 0 8 files changed, 279 insertions(+), 35 deletions(-) create mode 100644 beta_changes.csv create mode 100644 examples/introduce_greenspace.py create mode 100644 examples/run_sim_change_beta.py create mode 100644 examples/run_sim_change_trans.py diff --git a/beta_changes.csv b/beta_changes.csv new file mode 100644 index 000000000..533d1a533 --- /dev/null +++ b/beta_changes.csv @@ -0,0 +1,9 @@ +-,0.0 spd,0.01 spd,0.02 spd,0.03 spd,0.04 spd,0.05 spd,0.06 spd,0.07 spd,0.08 spd,0.09 spd +cum_infections,50526.0,51171.0,50952.0,52532.0,50717.0,46211.0,53092.0,48452.0,48836.0,49295.0 +cum_reinfections,4104.0,4442.0,4003.0,4524.0,3823.0,3411.0,4391.0,3708.0,3652.0,3477.0 +cum_infectious,41400.0,41950.0,41839.0,43339.0,41547.0,37344.0,43778.0,39324.0,39724.0,40360.0 +cum_symptomatic,27452.0,27155.0,27458.0,27812.0,27243.0,22255.0,25376.0,22217.0,22994.0,23977.0 +cum_severe,1638.0,1451.0,1351.0,1532.0,1405.0,640.0,749.0,590.0,933.0,837.0 +cum_critical,475.0,445.0,361.0,493.0,418.0,213.0,158.0,145.0,187.0,230.0 +cum_recoveries,23830.0,24619.0,24498.0,25618.0,24192.0,21840.0,26903.0,23505.0,23195.0,24283.0 +cum_deaths,78.0,54.0,61.0,60.0,48.0,36.0,19.0,36.0,22.0,11.0 diff --git a/covasim/defaults.py b/covasim/defaults.py index 707c26cd4..c5b1e51f3 100644 --- a/covasim/defaults.py +++ b/covasim/defaults.py @@ -220,7 +220,6 @@ def __init__(self): 'severe', ] -default_SVI_data = np.array([]) # Default age data, based on Seattle 2018 census data -- used in population.py default_age_data = np.array([ [ 0, 4, 0.0605], diff --git a/covasim/interventions.py b/covasim/interventions.py index 33cf50651..24c35b2ec 100644 --- a/covasim/interventions.py +++ b/covasim/interventions.py @@ -7,6 +7,7 @@ import pandas as pd import pylab as pl import sciris as sc +import random import inspect import datetime as dt from . import misc as cvm @@ -522,6 +523,89 @@ def apply(self, sim): return self.interventions[inds[0]].apply(sim) +#%% Socio-economic interventions + +__all__ += ['introduce_greenspace'] + +class introduce_greenspace(Intervention): + + def __init__(self, pop_size, num_greenspaces=None, custom_spaces=None, specifications=None, symp_prob_dec=0.1 , crit_prob_dec=0.1, severe_prob_dec=0.1, **kwargs): + super().__init__(**kwargs) # Initialize the Intervention object + self.specifications = sc.dcp(specifications) + self.pop_size = sc.dcp(pop_size) + + self.symp_prob_dec = sc.dcp(symp_prob_dec) + self.crit_prob_dec = sc.dcp(crit_prob_dec) + self.severe_prob_dec = sc.dcp(severe_prob_dec) + + #populate greenspaces + basic_greenspace = { + "daily_access": 1000, + "prob_daily_access": .9, + "max_daily_cap": 1000, + "daily_visits": [[]] + } + + if custom_spaces: self.greenspaces = custom_spaces + else: self.greenspaces = [basic_greenspace] + + if num_greenspaces and len(self.greenspaces) < num_greenspaces: + for i in range(num_greenspaces-len(self.greenspaces)): self.greenspaces.append(basic_greenspace) + + #assign IDs + for i in range(len(self.greenspaces)): self.greenspaces[i]['id'] = i + + + self.fitness_improvement = [0 * pop_size] + self.greenspace_visits = [{} for i in range(pop_size)] + self.greenspace_access = [[] for i in range(pop_size)] + + for greenspace in self.greenspaces: + loc = np.random.randint(self.pop_size-greenspace['daily_access']) + for i in range(loc,loc + greenspace['daily_access']): + self.greenspace_access[i] = self.greenspace_access[i] + [greenspace['id']] + + return + + def initialize(self, sim): + ''' Fix days and store beta ''' + + super().initialize() + + return + + def visit_greenspaces(self,sim): + ''' simulate people visiting greenspaces''' + + for person in range(self.pop_size): + for greenspace_id in self.greenspace_access[person]: + + gpda = 0 + for greenspace in self.greenspaces: + if greenspace['id'] == greenspace_id: gpda = greenspace['prob_daily_access'] + + person_visits_greenspace = random.choices([1,0], weights=[gpda,1-gpda], cum_weights=None, k=1)[0] + + if person_visits_greenspace and len(greenspace['daily_visits'][-1]) < greenspace['max_daily_cap']: + self.greenspaces[greenspace_id]['daily_visits'][-1].append(person) + break + + + + def apply(self, sim): + + self.visit_greenspaces(sim) + + for greenspace in self.greenspaces: + for visitor_id in greenspace['daily_visits'][-1]: + sim.people.symp_prob[visitor_id] = sim.people.symp_prob[visitor_id] - sim.people.symp_prob[visitor_id]*self.symp_prob_dec + sim.people.crit_prob[visitor_id] = sim.people.crit_prob[visitor_id] - sim.people.crit_prob[visitor_id]*self.crit_prob_dec + sim.people.severe_prob[visitor_id] = sim.people.severe_prob[visitor_id] - sim.people.severe_prob[visitor_id]*self.severe_prob_dec + greenspace['daily_visits'].append([]) + + return + + #%% Beta interventions __all__+= ['change_beta', 'clip_edges'] diff --git a/covasim/parameters.py b/covasim/parameters.py index bf6463034..d018df8d0 100644 --- a/covasim/parameters.py +++ b/covasim/parameters.py @@ -54,7 +54,7 @@ def make_pars(set_prognoses=False, prog_by_age=True, version=None, **kwargs): # Network parameters, generally initialized after the population has been constructed pars['contacts'] = None # The number of contacts per layer; set by reset_layer_pars() below pars['dynam_layer'] = None # Which layers are dynamic; set by reset_layer_pars() below - pars['beta_layer'] = None # Transmissibility per layer; set by reset_layer_pars() below + pars['beta_l ayer'] = None # Transmissibility per layer; set by reset_layer_pars() below # Basic disease transmission parameters pars['beta_dist'] = dict(dist='neg_binomial', par1=1.0, par2=0.45, step=0.01) # Distribution to draw individual level transmissibility; dispersion from https://www.researchsquare.com/article/rs-29548/v1 diff --git a/examples/example_data.csv b/examples/example_data.csv index dfcfae3f2..e69de29bb 100644 --- a/examples/example_data.csv +++ b/examples/example_data.csv @@ -1,33 +0,0 @@ -date,new_diagnoses,new_tests,new_deaths -2020-03-01,1,24,0 -2020-03-02,3,22,0 -2020-03-03,2,15,0 -2020-03-04,8,40,0 -2020-03-05,20,38,0 -2020-03-06,9,61,0 -2020-03-07,6,43,0 -2020-03-08,13,98,0 -2020-03-09,6,93,0 -2020-03-10,25,170,0 -2020-03-11,28,368,0 -2020-03-12,27,437,0 -2020-03-13,22,291,2 -2020-03-14,43,328,0 -2020-03-15,76,1147,0 -2020-03-16,65,1438,1 -2020-03-17,88,1209,0 -2020-03-18,86,1269,0 -2020-03-19,115,1195,1 -2020-03-20,51,529,0 -2020-03-21,55,482,3 -2020-03-22,95,1106,2 -2020-03-23,74,471,1 -2020-03-24,63,438,4 -2020-03-25,178,1111,2 -2020-03-26,83,621,1 -2020-03-27,140,1059,2 -2020-03-28,137,951,1 -2020-03-29,150,964,0 -2020-03-30,144,1058,1 -2020-04-03,145,1058,1 -2020-04-11,0,0,1 diff --git a/examples/introduce_greenspace.py b/examples/introduce_greenspace.py new file mode 100644 index 000000000..d56df0f5b --- /dev/null +++ b/examples/introduce_greenspace.py @@ -0,0 +1,104 @@ +''' +Simple script for running Covasim scenarios +''' + +import covasim as cv +import csv + +def write_result_dict_to_csv(result_dict,cols): + out = [] + for row in result_dict.keys(): + out.append([row] + list(result_dict[row].values())) + return out + + +def create_greenspaces(num_greenspaces=1, daily_access=[100],prob_daily_access=[0.5],max_daily_cap=[100]): + basic_greenspace = { + "daily_access": 1000, + "prob_daily_access": 0.1, + "max_daily_cap": 1000, + "daily_visits": [[]] + } + + greenspaces = [basic_greenspace for i in range(num_greenspaces)] + + for i in range(len(daily_access)) : greenspaces[i]['daily_access'] = daily_access[i] + for i in range(len(prob_daily_access)) : greenspaces[i]['prob_daily_access'] = prob_daily_access[i] + for i in range(len(max_daily_cap)) : greenspaces[i]['max_daily_cap'] = max_daily_cap[i] + + return greenspaces + + +# Run options +do_plot = 1 +do_show = 1 +verbose = 1 + +# Sim options +basepars = dict( + pop_size = 2000, + verbose = verbose, +) + +# Scenario metaparameters +metapars = dict( + n_runs = 3, # Number of parallel runs; change to 3 for quick, 11 for real + noise = 0.1, # Use noise, optionally + noisepar = 'beta', + rand_seed = 1, + quantiles = {'low':0.1, 'high':0.9}, +) + + +report = ['cum_infections','cum_reinfections','cum_infectious','cum_symptomatic','cum_severe','cum_critical','cum_recoveries','cum_deaths'] +runs = [] + + +# Define the actual scenarios +scenarios = {'baseline': { + 'name':'Baseline', + 'pars': { + 'interventions': None, + } + },} + +for i in range(0,10): + greenspaces = create_greenspaces(num_greenspaces = 4,prob_daily_access=[i/10 for j in range(4)]) + runs.append(str(i/100)+" spd") + scenarios[str(i/100)+" spd"] = { + 'name' : str(i/100) + "spd", + 'pars': { + 'interventions': cv.introduce_greenspace(num_greenspaces=i,custom_spaces=greenspaces,pop_size=2000,symp_prob_dec=1/100) + } + } + + + +# Run the scenarios -- this block is required for parallel processing on Windows +if __name__ == "__main__": + + scens = cv.Scenarios(basepars=basepars, metapars=metapars, scenarios=scenarios) + scens.run(verbose=verbose) + + results = scens.results + final_results = {} + for r in report: + final_report = {} + single_results = results[r] + for run in runs: + single_report_results = single_results[run] + final_report[run] = sum(single_report_results['best']) + print("SINGLE",sum(single_report_results['best'])) + final_results[r] = final_report + + + #summarize the beta changes in a csv file + with open('beta_changes.csv', 'w', encoding='UTF8', newline='') as f: + # create the csv writer + writer = csv.writer(f) + writer.writerow(['-'] + runs) + for row in write_result_dict_to_csv(final_results,report): + writer.writerow(row) + if do_plot: + fig1 = scens.plot(do_show=do_show) + diff --git a/examples/run_sim_change_beta.py b/examples/run_sim_change_beta.py new file mode 100644 index 000000000..59eb1dc10 --- /dev/null +++ b/examples/run_sim_change_beta.py @@ -0,0 +1,81 @@ +''' +Simple script for running Covasim scenarios +''' + +import covasim as cv +import csv + + +# Run options +do_plot = 1 +do_show = 1 +verbose = 1 + +# Sim options +basepars = dict( + pop_size = 2000, + verbose = verbose, +) + +# Scenario metaparameters +metapars = dict( + n_runs = 3, # Number of parallel runs; change to 3 for quick, 11 for real + noise = 0.1, # Use noise, optionally + noisepar = 'beta', + rand_seed = 1, + quantiles = {'low':0.1, 'high':0.9}, +) + +report = ['cum_infections','cum_reinfections','cum_infectious','cum_symptomatic','cum_severe','cum_critical','cum_recoveries','cum_deaths','cum_tests','cum_diagnoses','cum_known_deaths','cum_quarantined','cum_isolated'] +runs = [] +# Define the actual scenarios +start_day = '2020-04-04' +scenarios = {} +for i in range(1,11): + runs.append(str(i)+"beta") + scenarios[str(i)+"beta"] = { + 'name' : str(i) + "beta", + 'pars': { + 'interventions': cv.change_beta(days=start_day, changes=i/100) + } + } + +def write_result_dict_to_csv(result_dict,cols): + out = [] + for row in result_dict.keys(): + out.append([row] + list(result_dict[row].values())) + return out + +# Run the scenarios -- this block is required for parallel processing on Windows +if __name__ == "__main__": + + scens = cv.Scenarios(basepars=basepars, metapars=metapars, scenarios=scenarios) + #scens.run(verbose=verbose) + scens.run() + results = scens.results + final_results = {} + for r in report: + final_report = {} + single_results = results[r] + for run in runs: + single_report_results = single_results[run] + final_report[run] = sum(single_report_results['best']) + print("SINGLE",sum(single_report_results['best'])) + final_results[r] = final_report + + print(final_results) + + #summarize the beta changes in a csv file + with open('beta_changes.csv', 'w', encoding='UTF8', newline='') as f: + # create the csv writer + writer = csv.writer(f) + writer.writerow(['-'] + runs) + for row in write_result_dict_to_csv(final_results,report): + writer.writerow(row) + + + print("RESULTS", scens.results) + if do_plot: + fig1 = scens.plot(do_show=do_show,style='simple') + + diff --git a/examples/run_sim_change_trans.py b/examples/run_sim_change_trans.py new file mode 100644 index 000000000..e69de29bb