added SVI as a factor to consider within the population

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

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

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']

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
@@ -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

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
 
 

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':

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

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)
+

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)