Merge pull request #373 from amath-idm/debug-treatment

Fix screening bug
InstituteforDiseaseModeling · Oct 26, 2022 · 2772dcf · 2772dcf
2 parents d89048c + 3153089
commit 2772dcf
Show file tree

Hide file tree

Showing 11 changed files with 255 additions and 108 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -8,9 +8,17 @@ All notable changes to the codebase are documented in this file. Changes that ma
    :local:
    :depth: 1
 
+
+Version 0.3.1 (2022-10-26)
+--------------------------
+- Fixes bug with screening
+- Increases coverage of baseline test
+- *GitHub info*: PR `373 <https://github.com/amath-idm/hpvsim/pull/373>`__
+
+
 Version 0.3.0 (2022-10-26)
 --------------------------
-- Implements multiscale modeling, but currently disabled
+- Implements multiscale modeling
 - Minor release tidying
 - *GitHub info*: PR `365 <https://github.com/amath-idm/hpvsim/pull/365>`__
 

diff --git a/hpvsim/analysis.py b/hpvsim/analysis.py
@@ -12,6 +12,7 @@
 from . import plotting as hppl
 from . import defaults as hpd
 from . import parameters as hppar
+from . import interventions as hpi
 from .settings import options as hpo # For setting global options
 
 
@@ -1034,17 +1035,16 @@ class cancer_detection(Analyzer):
         treat_prob: Probability of receiving treatment for those with symptom-detected cancer
     '''
 
-    def __init__(self, symp_prob=0.01, treat_prob=0.01, **kwargs):
+    def __init__(self, symp_prob=0.01, treat_prob=0.01, product=None, **kwargs):
         super().__init__(**kwargs)
         self.symp_prob = symp_prob
         self.treat_prob = treat_prob
+        self.product = product or hpi.radiation()
 
     def initialize(self, sim):
         super().initialize(sim)
         self.dt = sim['dt']
 
-        # Add entries to results
-
 
     def apply(self, sim):
         '''
@@ -1066,12 +1066,8 @@ def apply(self, sim):
                 sim.people.date_detected_cancer[is_detected_inds] = sim.t
                 treat_probs = np.full(len(is_detected_inds), self.treat_prob)
                 treat_inds = is_detected_inds[hpu.binomial_arr(treat_probs)]
-                if 'cancer_treatment' in sim.pars['treat_pars'].keys():
-                    new_dur_cancer = hpu.sample(**sim.pars['treat_pars']['cancer_treatment']['dur'], size=len(treat_inds))
-                    sim.people.date_dead_cancer[treat_inds] += np.ceil(new_dur_cancer / self.dt)
-                    sim.people.treated[treat_inds] = True
-                    sim.people.date_treated[treat_inds] = sim.t
-                    new_treatments = len(treat_inds)
+                if len(treat_inds)>0:
+                    self.product.administer(sim, treat_inds)
 
         # Update flows
         sim.people.flows['detected_cancers'] = new_detections

diff --git a/hpvsim/interventions.py b/hpvsim/interventions.py
@@ -309,11 +309,12 @@ class RoutineDelivery(Intervention):
     '''
     Base class for any intervention that uses routine delivery; handles interpolation of input years.
     '''
-    def __init__(self, years=None, start_year=None, end_year=None, prob=None):
+    def __init__(self, years=None, start_year=None, end_year=None, prob=None, annual_prob=True):
         self.years      = years
         self.start_year = start_year
         self.end_year   = end_year
         self.prob       = sc.promotetoarray(prob)
+        self.annual_prob = annual_prob # Determines whether the probability is annual or per timestep
         return
 
     def initialize(self, sim):
@@ -353,8 +354,8 @@ def initialize(self, sim):
         else:
             self.prob = sc.smoothinterp(self.yearvec, self.years, self.prob, smoothness=0)
 
-        # Lastly, adjust the annual probability by the sim's timestep
-        self.prob = self.prob*sim['dt']
+        # Lastly, adjust the annual probability by the sim's timestep, if it's an annual probability
+        if self.annual_prob: self.prob = self.prob*sim['dt']
 
         return
 
@@ -363,10 +364,11 @@ class CampaignDelivery(Intervention):
     '''
     Base class for any intervention that uses campaign delivery; handles interpolation of input years.
     '''
-    def __init__(self, years, interpolate=True, prob=None):
+    def __init__(self, years, interpolate=True, prob=None, annual_prob=True):
         self.years = sc.promotetoarray(years)
         self.interpolate = interpolate # Whether to space the intervention over the year (if true) or do them all at once (if false)
         self.prob       = sc.promotetoarray(prob)
+        self.annual_prob = annual_prob
         return
 
     def initialize(self, sim):
@@ -381,13 +383,16 @@ def initialize(self, sim):
 
         # Get the probability input into a format compatible with timepoints
         if len(self.prob) == len(self.years) and self.interpolate:
-            self.prob = sc.smoothinterp(np.arange(len(self.timepoints)), np.arange(len(self.years)), self.prob, smoothness=0)*sim['dt']
+            self.prob = sc.smoothinterp(np.arange(len(self.timepoints)), np.arange(len(self.years)), self.prob, smoothness=0)
         elif len(self.prob) == 1:
             self.prob = np.array([self.prob[0]] * len(self.timepoints))
         else:
             errormsg = f'Length of years incompatible with length of probabilities: {len(self.years)} vs {len(self.prob)}'
             raise ValueError(errormsg)
 
+        # Lastly, adjust the annual probability by the sim's timestep, if it's an annual probability
+        if self.annual_prob: self.prob = self.prob*sim['dt']
+
         return
 
 
@@ -751,6 +756,7 @@ def deliver(self, sim):
         '''
         Deliver the diagnostics by finding who's eligible, finding who accepts, and applying the product.
         '''
+        self.outcomes = {k:np.array([], dtype=hpd.default_int) for k in self.product.hierarchy}
         ti = sc.findinds(self.timepoints, sim.t)[0]
         prob            = self.prob[ti] # Get the proportion of people who will be tested this timestep
         eligible_inds   = self.check_eligibility(sim) # Check eligibility
@@ -877,9 +883,9 @@ class routine_triage(BaseTriage, RoutineDelivery):
         triage2 = hpv.routine_triage(product='pos_screen_assessment', eligibility=screen_pos, prob=0.9, start_year=2030)
     '''
     def __init__(self, product=None, prob=None, eligibility=None, age_range=None,
-                         years=None, start_year=None, end_year=None, **kwargs):
+                         years=None, start_year=None, end_year=None, annual_prob=None, **kwargs):
         BaseTriage.__init__(self, product=product, age_range=age_range, eligibility=eligibility, **kwargs)
-        RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years)
+        RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years, annual_prob=annual_prob)
 
     def initialize(self, sim):
         RoutineDelivery.initialize(self, sim) # Initialize this first, as it ensures that prob is interpolated properly
@@ -898,9 +904,9 @@ class campaign_triage(BaseTriage, CampaignDelivery):
         triage1 = hpv.campaign_triage(product='pos_screen_assessment', eligibility=screen_pos, prob=0.9, years=2030)
     '''
     def __init__(self, product=None, age_range=None, sex=None, eligibility=None,
-                 prob=None, years=None, interpolate=None, **kwargs):
+                 prob=None, years=None, interpolate=None, annual_prob=None, **kwargs):
         BaseTriage.__init__(self, product=product, age_range=age_range, sex=sex, eligibility=eligibility, **kwargs)
-        CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate)
+        CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate, annual_prob=annual_prob)
 
     def initialize(self, sim):
         CampaignDelivery.initialize(self, sim) # Initialize this first, as it ensures that prob is interpolated properly
@@ -1117,9 +1123,9 @@ class routine_txvx(BaseTxVx, RoutineDelivery):
     '''
 
     def __init__(self, product=None, prob=None, age_range=None, eligibility=None,
-                 start_year=None, end_year=None, years=None, **kwargs):
+                 start_year=None, end_year=None, years=None, annual_prob=None, **kwargs):
         BaseTxVx.__init__(self, product=product, age_range=age_range, eligibility=eligibility, **kwargs)
-        RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years)
+        RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years, annual_prob=annual_prob)
 
     def initialize(self, sim):
         RoutineDelivery.initialize(self, sim) # Initialize this first, as it ensures that prob is interpolated properly
@@ -1136,9 +1142,9 @@ class campaign_txvx(BaseTxVx, CampaignDelivery):
     '''
 
     def __init__(self, product=None, prob=None, age_range=None, eligibility=None,
-                 years=None, interpolate=True, **kwargs):
+                 years=None, interpolate=True, annual_prob=None, **kwargs):
         BaseScreening.__init__(self, product=product, age_range=age_range, eligibility=eligibility, **kwargs)
-        CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate)
+        CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate, annual_prob=annual_prob)
 
     def initialize(self, sim):
         CampaignDelivery.initialize(self, sim) # Initialize this first, as it ensures that prob is interpolated properly
@@ -1162,7 +1168,7 @@ def apply(self, sim):
 
 
 #%% Products
-__all__ += ['dx', 'tx', 'vx', 'radiation']
+__all__ += ['dx', 'tx', 'vx', 'radiation', 'default_dx', 'default_tx', 'default_vx']
 
 class Product(hpb.FlexPretty):
     ''' Generic product implementation '''

diff --git a/hpvsim/parameters.py b/hpvsim/parameters.py
@@ -70,17 +70,15 @@ def make_pars(**kwargs):
     pars['n_partner_types'] = 1  # Number of partnership types - reset below
 
     # Basic disease transmission parameters
-    pars['beta_dist']       = dict(dist='neg_binomial', par1=1.0, par2=1.0, step=0.01) # Distribution to draw individual level transmissibility TODO does this get used? if not remove.
     pars['beta']            = 0.25  # Per-act transmission probability; absolute value, calibrated
     pars['transf2m']        = 1.0   # Relative transmissibility of receptive partners in penile-vaginal intercourse; baseline value
     pars['transm2f']        = 3.69  # Relative transmissibility of insertive partners in penile-vaginal intercourse; based on https://doi.org/10.1038/srep10986: "For vaccination types, the risk of male-to-female transmission was higher than that of female-to-male transmission"
 
     # Parameters for disease progression
-    pars['severity_dist'] = dict(dist='lognormal', par1=None, par2=0.1) # Distribution of individual disease severity. Par1 is set to None because the mean is determined as a function of genotype and disease duration
     pars['clinical_cutoffs']    = {'cin1': 0.33, 'cin2':0.67, 'cin3':0.99} # Parameters the control the clinical cliassification of dysplasia
-    pars['cancer_treat_prob']   = 0.1 # probability of receiving cancer treatment given symptom detection
     pars['hpv_control_prob']    = 0.44 # Probability that HPV is controlled latently vs. cleared
     pars['hpv_reactivation']    = 0.025 # Placeholder
+    pars['dur_cancer']          = dict(dist='lognormal', par1=12.0, par2=3.0)  # Duration of untreated invasive cerival cancer before death (years)
 
     # Parameters used to calculate immunity
     pars['imm_init']        = dict(dist='beta', par1=5, par2=3)  # beta distribution for initial level of immunity following infection clearance
@@ -103,20 +101,11 @@ def make_pars(**kwargs):
     pars['vaccine_pars']    = dict()  # Vaccines that are being used; populated during initialization
     pars['vaccine_map']     = dict()  # Reverse mapping from number to vaccine key
 
-    # Screening and treatment parameters
-    pars['screen_pars']     = dict()  # Screening method that is being used; populated during initialization
-    pars['treat_pars'] = dict()  # Treatment method that is being used; populated during initialization
-
-    # Durations
-    pars['dur_cancer']      = dict(dist='lognormal', par1=12.0, par2=3.0)  # Duration of untreated invasive cerival cancer before death (years)
-
     # Parameters determining relative transmissibility at each stage of disease
-    pars['rel_trans'] = {}
-    pars['rel_trans']['precin']     = 1 # Baseline value
-    pars['rel_trans']['cin1']       = 1 # Baseline assumption, can be adjusted during calibration
-    pars['rel_trans']['cin2']       = 1 # Baseline assumption, can be adjusted during calibration
-    pars['rel_trans']['cin3']       = 1 # Baseline assumption, can be adjusted during calibration
-    pars['rel_trans']['cancerous']  = 0.5 # Baseline assumption, can be adjusted during calibration
+    pars['rel_trans_cin1']       = 1 # Transmissibility of people with CIN1 compared to those without dysplasia
+    pars['rel_trans_cin2']       = 1 # Transmissibility of people with CIN2 compared to those without dysplasia
+    pars['rel_trans_cin3']       = 1 # Transmissibility of people with CIN3 compared to those without dysplasia
+    pars['rel_trans_cancerous']  = 0.5 # Transmissibility of people with cancer compared to those without dysplasia
 
     # Efficacy of protection
     pars['eff_condoms']     = 0.7  # The efficacy of condoms; https://www.nejm.org/doi/10.1056/NEJMoa053284?url_ver=Z39.88-2003&rfr_id=ori:rid:crossref.org&rfr_dat=cr_pub%20%200www.ncbi.nlm.nih.gov

diff --git a/hpvsim/sim.py b/hpvsim/sim.py
@@ -402,7 +402,7 @@ def init_res(*args, **kwargs):
         for lkey,llab,cstride,g in zip(['total_',''], ['Total ',''], [0.95,np.linspace(0.2,0.8,ng)], [0,ng]):  # key, label, and color stride by level (total vs genotype-specific)
             for flow in hpd.flows:
                 if (flow.by_genotype and lkey=='') or lkey=='total_':
-                    results[f'{lkey + flow.name}'] = init_res(f'{llab} {flow.label}', color=flow.cmap(cstride), n_rows=g)
+                    results[f'{lkey + flow.name}'] = init_res(f'{(llab + flow.label).capitalize()}', color=flow.cmap(cstride), n_rows=g)
 
         # Create stocks
         for llabel,cstride,g in zip(['Total number','Number'], [0.95,np.linspace(0.2,0.8,ng)], [0,ng]):
@@ -684,12 +684,11 @@ def step(self):
         hiv_rel_sus[hiv_inds] *= mod
 
         # Calculate relative transmissibility by stage of infection
-        rel_trans_pars = self['rel_trans']
         rel_trans = people.infectious[:].astype(hpd.default_float)
-        rel_trans[people.cin1] *= rel_trans_pars['cin1']
-        rel_trans[people.cin2] *= rel_trans_pars['cin2']
-        rel_trans[people.cin3] *= rel_trans_pars['cin3']
-        rel_trans[people.cancerous] *= rel_trans_pars['cancerous']
+        rel_trans[people.cin1] *= self['rel_trans_cin1']
+        rel_trans[people.cin2] *= self['rel_trans_cin2']
+        rel_trans[people.cin3] *= self['rel_trans_cin3']
+        rel_trans[people.cancerous] *= self['rel_trans_cancerous']
 
         inf = people.infectious.copy() # calculate transmission based on infectiousness at start of timestep i.e. someone infected in one layer cannot transmit the infection via a different layer in the same timestep
 

diff --git a/hpvsim/version.py b/hpvsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '0.3.0'
+__version__ = '0.3.1'
 __versiondate__ = '2022-10-26'
 __license__ = f'HPVsim {__version__} ({__versiondate__}) — © 2022 by IDM'
diff --git a/tests/baseline.json b/tests/baseline.json
@@ -1,61 +1,68 @@
 {
   "summary": {
-    "total_infections": 2131.0,
-    "total_cin1s": 446.0,
-    "total_cin2s": 435.0,
-    "total_cin3s": 217.0,
-    "total_cins": 1098.0,
-    "total_cancers": 33.0,
+    "total_infections": 1291.0,
+    "total_cin1s": 298.0,
+    "total_cin2s": 271.0,
+    "total_cin3s": 162.0,
+    "total_cins": 731.0,
+    "total_cancers": 26.0,
     "total_detected_cancers": 0.0,
-    "total_cancer_deaths": 22.0,
+    "total_cancer_deaths": 10.0,
     "total_detected_cancer_deaths": 0.0,
-    "total_reinfections": 727.0,
-    "total_reactivations": 423.0,
+    "total_reinfections": 431.0,
+    "total_reactivations": 242.0,
     "total_hiv_infections": 0.0,
-    "n_total_susceptible": 68602.0,
-    "n_total_infectious": 4588.0,
-    "n_total_inactive": 12371.0,
-    "n_total_no_dysp": 55375.0,
-    "n_total_cin1": 591.0,
-    "n_total_cin2": 615.0,
-    "n_total_cin3": 631.0,
-    "n_total_cancerous": 314.0,
-    "n_total_infected": 15479.0,
-    "n_total_cin": 1790.0,
-    "n_total_precin": 2638.0,
-    "n_total_latent": 12371.0,
+    "n_total_susceptible": 38774.0,
+    "n_total_infectious": 2884.0,
+    "n_total_inactive": 7234.0,
+    "n_total_no_dysp": 31229.0,
+    "n_total_cin1": 397.0,
+    "n_total_cin2": 393.0,
+    "n_total_cin3": 427.0,
+    "n_total_cancerous": 195.0,
+    "n_total_infected": 9262.0,
+    "n_total_cin": 1173.0,
+    "n_total_precin": 1612.0,
+    "n_total_latent": 7234.0,
     "n_detected_cancer": 0.0,
     "n_screened": 0.0,
     "n_treated": 0.0,
     "n_vaccinated": 0.0,
     "n_tx_vaccinated": 0.0,
     "n_hiv": 0.0,
     "n_art_adherence": 0.0,
-    "total_hpv_incidence": 0.031063234308037667,
-    "total_cin1_incidence": 1989.9165662784992,
-    "total_cin2_incidence": 1940.8379065720787,
-    "total_cin3_incidence": 968.1881051175657,
-    "total_cin_incidence": 4898.942577968143,
-    "total_cancer_incidence": 147.23597911926115,
-    "births": 1235.0,
-    "other_deaths": 393.0,
+    "total_hpv_incidence": 0.03329550729870532,
+    "total_cin1_incidence": 2316.002176109427,
+    "total_cin2_incidence": 2106.1630527706534,
+    "total_cin3_incidence": 1259.0347400326416,
+    "total_cin_incidence": 5681.199968912722,
+    "total_cancer_incidence": 202.06730395585606,
+    "births": 708.0,
+    "other_deaths": 219.0,
     "migration": 0.0,
-    "asr_cancer": 153.95065201479377,
+    "asr_cancer": 189.67339045604984,
     "new_total_vaccinated": 0.0,
     "cum_total_vaccinated": 0.0,
-    "new_doses": 0.0,
-    "cum_doses": 0.0,
-    "new_txvx_doses": 0.0,
-    "new_tx_vaccinated": 0.0,
-    "cum_txvx_doses": 0.0,
-    "cum_tx_vaccinated": 0.0,
+    "new_doses": 6.0,
+    "cum_doses": 107.0,
+    "new_txvx_doses": 1279.0,
+    "new_tx_vaccinated": 1279.0,
+    "cum_txvx_doses": 1279.0,
+    "cum_tx_vaccinated": 1279.0,
     "detected_cancer_incidence": 0.0,
-    "cancer_mortality": 96.78412740310588,
-    "n_alive": 44863.0,
-    "cdr": 0.008760002674809976,
-    "cbr": 0.027528252680382496,
+    "cancer_mortality": 76.53451706719731,
+    "n_alive": 25693.0,
+    "cdr": 0.008523722414665473,
+    "cbr": 0.027556143696726734,
     "hiv_incidence": 0.0,
     "hiv_prevalence": 0.0,
-    "total_hpv_prevalence": 0.10226690145554243
+    "total_hpv_prevalence": 0.11224847234655354,
+    "resources_screen": 63.0,
+    "resources_triage": 1.0,
+    "resources_assign_tx": 0.0,
+    "resources_treat_num": 0.0,
+    "resources_treat_delay": 0.0,
+    "resources_None": 6.0,
+    "resources_routine_txvx": 1279.0
   }
 }
diff --git a/tests/benchmark.json b/tests/benchmark.json
@@ -1,12 +1,14 @@
 {
   "time": {
     "initialize": 0.014,
-    "run": 2.052
+    "run": 2.203
   },
   "parameters": {
-    "n_agents": 20000,
+    "n_agents": 10000,
     "n_genotypes": 2,
-    "n_years": 40
+    "n_years": 40,
+    "n_interventions": 7,
+    "n_analyzers": 0
   },
-  "cpu_performance": 0.9683884406244703
-}
+  "cpu_performance": 0.32539441404329295
+}