Increase range of options for yearly rescaling of HR capabilities (#1297)

* Increase range of options for yearly rescaling of HR capabilities

* Style fix

* Fix cond statement

* refactor to put two resource files for scaling into same folder

* elaborate comment

* move checks to point of reading-in data

* update comments concerning the first scheduled occurennce of DynamicRescalinHRCapabilities

* guarantee correct read-in of types

* tidy formatting of xlsx

* engross comment

* let the default be called something more informative ('no_scaling') and provide 'scaling_by_population_growth' to give easy access to that configuration

* check that there is an entry that is early enough for the first call in the simulation

* add test for using sequence of factors for dynamic scaling

* [OPTIONAL - DISCUSS WITH MM] possible refactor for readability and speed (avoid access DataFrame using .loc and .iloc repeatedely)

* isort linting

* ruff linting

---------

Co-authored-by: Tim Hallett <[email protected]>

Author: marghe-molaro

resources/healthsystem/ResourceFile_HealthSystem_parameters.csv

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:e36cbd225191f893f2de5f0f34adc251046af5ec58daaf7c86b09a6c83c1e31d
-size 379
+oid sha256:0efd27feb7dbb09331cb29ebeb10d8428981c7246c4fdee5ebacf510b9514795
+size 372

resources/healthsystem/human_resources/const_HR_scaling/ResourceFile_const_HR_scaling.xlsx resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_const_HR_scaling.xlsx

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:75ac6ae6b873d2b6e95090da64da98c5646e21a6c485ee6a10d5b8bc0e352950
+size 10780

resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_dynamic_HR_scaling.xlsx

@@ -554,14 +554,24 @@ class HealthSystem(Module):
                           "and custom (user can freely set these factors as parameters in the analysis).",
         ),
 
-        'dynamic_HR_scaling_factor': Parameter(
-            Types.REAL, "Factor by which HR capabilities are scaled at regular intervals of 1 year (in addition to the"
-                        " [optional] scaling with population size (controlled by `scale_HR_by_popsize`"
+        'yearly_HR_scaling': Parameter(
+            Types.DICT, "Factors by which HR capabilities are scaled. "
+                        "Each sheet specifies a 'mode' for dynamic HR scaling. The mode to use is determined by the "
+                        "parameter `yearly_HR_scaling_mode`. Each sheet must have the same format, including the same "
+                        "column headers. On each sheet, the first row (for `2010`, when the simulation starts) "
+                        "specifies the initial configuration: `dynamic_HR_scaling_factor` (float) is the factor by "
+                        "which all human resoucres capabilities and multiplied; `scale_HR_by_popsize` (bool) specifies "
+                        "whether the capabilities should (also) grow by the factor by which the population has grown in"
+                        " the last year. Each subsequent row specifies a year where there should be a CHANGE in the "
+                        "configuration. If there are no further rows, then there is no change. But, for example, an"
+                        " additional row of the form ```2015, 1.05, TRUE``` would mean that on 1st January of 2015, "
+                        "2016, 2017, ....(and the rest of the simulation), the capabilities would increase by the "
+                        "product of 1.05 and by the ratio of the population size to that in the year previous."
         ),
 
-        'scale_HR_by_popsize': Parameter(
-            Types.BOOL, "Decide whether to scale HR capabilities by population size every year. Can be used as well as"
-                        " the dynamic_HR_scaling_factor"
+        'yearly_HR_scaling_mode': Parameter(
+            Types.STRING, "Specifies which of the policies in yearly_HR_scaling should be adopted. This corresponds to"
+                          "a worksheet of the file `ResourceFile_dynamic_HR_scaling.xlsx`."
         ),
 
         'tclose_overwrite': Parameter(
@@ -685,6 +695,7 @@ def __init__(
                                        'LCOA_EHP']
         self.arg_policy_name = policy_name
 
+
         self.tclose_overwrite = None
         self.tclose_days_offset_overwrite = None
 
@@ -820,19 +831,36 @@ def read_parameters(self, data_folder):
         self.parameters['const_HR_scaling_table']: Dict = pd.read_excel(
             path_to_resourcefiles_for_healthsystem /
             "human_resources" /
-            "const_HR_scaling" /
+            "scaling_capabilities" /
             "ResourceFile_const_HR_scaling.xlsx",
             sheet_name=None  # all sheets read in
         )
+        # Ensure the mode of HR scaling to be considered in included in the tables loaded
+        assert self.parameters['const_HR_scaling_mode'] in self.parameters['const_HR_scaling_table'], \
+            f"Value of `const_HR_scaling_mode` not recognised: {self.parameters['const_HR_scaling_mode']}"
+
+        self.parameters['yearly_HR_scaling']: Dict = pd.read_excel(
+            path_to_resourcefiles_for_healthsystem /
+            "human_resources" /
+            "scaling_capabilities" /
+            "ResourceFile_dynamic_HR_scaling.xlsx",
+            sheet_name=None,  # all sheets read in
+            dtype={
+                'year': int,
+                'dynamic_HR_scaling_factor': float,
+                'scale_HR_by_popsize': bool
+            }  # Ensure that these column are read as the right type
+        )
+        # Ensure the mode of yearly HR scaling to be considered in included in the tables loaded
+        assert self.parameters['yearly_HR_scaling_mode'] in self.parameters['yearly_HR_scaling'], \
+            f"Value of `yearly_HR_scaling` not recognised: {self.parameters['yearly_HR_scaling_mode']}"
+        # Ensure that a value for the year at the start of the simulation is provided.
+        assert all(2010 in sheet['year'].values for sheet in self.parameters['yearly_HR_scaling'].values())
 
 
     def pre_initialise_population(self):
         """Generate the accessory classes used by the HealthSystem and pass to them the data that has been read."""
 
-        # Ensure the mode of HR scaling to be considered in included in the tables loaded
-        assert self.parameters['const_HR_scaling_mode'] in self.parameters['const_HR_scaling_table'], \
-            f"Value of `const_HR_scaling_mode` not recognised: {self.parameters['const_HR_scaling_mode']}"
-
         # Create dedicated RNGs for separate functions done by the HealthSystem module
         self.rng_for_hsi_queue = np.random.RandomState(self.rng.randint(2 ** 31 - 1))
         self.rng_for_dx = np.random.RandomState(self.rng.randint(2 ** 31 - 1))
@@ -914,10 +942,10 @@ def initialise_simulation(self, sim):
         sim.schedule_event(HealthSystemChangeMode(self),
                            Date(self.parameters["year_mode_switch"], 1, 1))
 
-        # Schedule recurring event which will rescale daily capabilities at regular intervals.
-        # The first event scheduled will only be used to update self.last_year_pop_size parameter,
-        # actual scaling will only take effect from 2011 onwards
-        sim.schedule_event(DynamicRescalingHRCapabilities(self), Date(sim.date) + pd.DateOffset(years=1))
+        # Schedule recurring event which will rescale daily capabilities (at yearly intervals).
+        # The first event scheduled for the start of the simulation is only used to update self.last_year_pop_size,
+        # whilst the actual scaling will only take effect from 2011 onwards.
+        sim.schedule_event(DynamicRescalingHRCapabilities(self), Date(sim.date))
 
     def on_birth(self, mother_id, child_id):
         self.bed_days.on_birth(self.sim.population.props, mother_id, child_id)
@@ -2861,27 +2889,42 @@ class DynamicRescalingHRCapabilities(RegularEvent, PopulationScopeEventMixin):
     """ This event exists to scale the daily capabilities assumed at fixed time intervals"""
     def __init__(self, module):
         super().__init__(module, frequency=DateOffset(years=1))
-        self.last_year_pop_size = self.current_pop_size  # store population size at initiation (when this class is
-        #                                                  created)
+        self.last_year_pop_size = self.current_pop_size  # will store population size at initiation (when this class is
+        #                                                  created, at the start of the simulation)
+
+        # Store the sequence of updates as a dict of the form
+        #   {<year_of_change>: {`dynamic_HR_scaling_factor`: float, `scale_HR_by_popsize`: bool}}
+        self.scaling_values = self.module.parameters['yearly_HR_scaling'][
+            self.module.parameters['yearly_HR_scaling_mode']].set_index("year").to_dict("index")
 
     @property
     def current_pop_size(self) -> float:
         """Returns current population size"""
         df = self.sim.population.props
         return df.is_alive.sum()
 
-    def apply(self, population):
+    def _get_most_recent_year_specified_for_a_change_in_configuration(self) -> int:
+        """Get the most recent year (in the past), for which there is an entry in `parameters['yearly_HR_scaling']`."""
+        years = np.array(list(self.scaling_values.keys()))
+        return years[years <= self.sim.date.year].max()
 
+    def apply(self, population):
+        """Do the scaling on the capabilities based on instruction that is in force at this time."""
+        # Get current population size
         this_year_pop_size = self.current_pop_size
 
-        # Rescale by fixed amount
-        self.module._daily_capabilities *= self.module.parameters['dynamic_HR_scaling_factor']
+        # Get the configuration to apply now (the latest entry in the `parameters['yearly_HR_scaling']`)
+        config = self.scaling_values.get(self._get_most_recent_year_specified_for_a_change_in_configuration())
+
+        # ... Do the rescaling specified for this year by the specified factor
+        self.module._daily_capabilities *= config['dynamic_HR_scaling_factor']
 
-        # Rescale daily capabilities by population size, if this option is included
-        if self.module.parameters['scale_HR_by_popsize']:
-            self.module._daily_capabilities *= this_year_pop_size/self.last_year_pop_size
+        # ... If requested, also do the scaling for the population growth that has occurred since the last year
+        if config['scale_HR_by_popsize']:
+            self.module._daily_capabilities *= this_year_pop_size / self.last_year_pop_size
 
-        self.last_year_pop_size = this_year_pop_size   # Save for next year
+        # Save current population size as that for 'last year'.
+        self.last_year_pop_size = this_year_pop_size
 
 
 class HealthSystemChangeMode(RegularEvent, PopulationScopeEventMixin):

src/tlo/methods/healthsystem.py

@@ -1,7 +1,7 @@
 import heapq as hp
 import os
 from pathlib import Path
-from typing import Set, Tuple
+from typing import Dict, Set, Tuple
 
 import numpy as np
 import pandas as pd
@@ -2244,8 +2244,8 @@ def get_hsi_log(service_availability, randomise_hsi_queue) -> pd.DataFrame:
 
 
 def test_const_HR_scaling_assumption(seed, tmpdir):
-    """Check that we can use the parameter `const_HR_scaling_mode` to manipulate the minutes of time available for healthcare
-    workers."""
+    """Check that we can use the parameter `const_HR_scaling_mode` to manipulate the minutes of time available for
+    healthcare workers."""
 
     def get_capabilities_today(const_HR_scaling_mode: str) -> pd.Series:
         sim = Simulation(start_date=start_date, seed=seed)
@@ -2292,25 +2292,28 @@ def get_initial_capabilities() -> pd.Series:
 
         return sim.modules['HealthSystem'].capabilities_today
 
-    def get_capabilities_after_two_years(dynamic_HR_scaling_factor: float, scale_HR_by_pop_size: bool) -> tuple:
+    def get_capabilities_after_two_updates(dynamic_HR_scaling_factor: float, scale_HR_by_pop_size: bool) -> tuple:
         sim = Simulation(start_date=start_date, seed=seed)
         sim.register(
             demography.Demography(resourcefilepath=resourcefilepath),
             healthsystem.HealthSystem(resourcefilepath=resourcefilepath),
             simplified_births.SimplifiedBirths(resourcefilepath=resourcefilepath),
 
         )
-        sim.modules['HealthSystem'].parameters['dynamic_HR_scaling_factor'] = dynamic_HR_scaling_factor
-        sim.modules['HealthSystem'].parameters['scale_HR_by_popsize'] = scale_HR_by_pop_size
-        sim.make_initial_population(n=100)
-
-        # Ensure simulation lasts long enough so that current capabilities reflect that used in the third year of
-        # simulation (i.e. after two annual updates)
-        sim.simulate(end_date=start_date + pd.DateOffset(years=2, days=1))
+        params = sim.modules['HealthSystem'].parameters
+        df = params['yearly_HR_scaling'][params['yearly_HR_scaling_mode']]
+        df.loc[df['year'] == 2010, 'dynamic_HR_scaling_factor'] = dynamic_HR_scaling_factor
+        df.loc[df['year'] == 2010, 'scale_HR_by_popsize'] = scale_HR_by_pop_size
+        popsize = 100
+        sim.make_initial_population(n=popsize)
 
-        popsize = sim.modules['Demography'].popsize_by_year
+        # Ensure simulation lasts long enough so that current capabilities reflect that used after two updates
+        # (updates occur on 1st Jan, starting in 2010, so simulation should stop on 2nd Jan 2011).
+        sim.simulate(end_date=Date(2011, 1, 2))
 
-        final_popsize_increase = popsize[2012]/popsize[2010]
+        popsize_start = popsize
+        popsize_curr = sim.population.props['is_alive'].sum()
+        final_popsize_increase = popsize_curr / popsize_start
 
         return sim.modules['HealthSystem'].capabilities_today, final_popsize_increase
 
@@ -2321,17 +2324,17 @@ def get_capabilities_after_two_years(dynamic_HR_scaling_factor: float, scale_HR_
     initial_caps = initial_caps[initial_caps != 0]
 
     # Check that dynamic expansion over two years leads to expansion = dynamic_HR_scaling_factor^2
-    caps, final_popsize_increase = get_capabilities_after_two_years(
+    caps, final_popsize_increase = get_capabilities_after_two_updates(
         dynamic_HR_scaling_factor=dynamic_HR_scaling_factor,
         scale_HR_by_pop_size=False
     )
     caps = caps[caps != 0]
     ratio_in_sim = caps/initial_caps
-    expected_value = dynamic_HR_scaling_factor*dynamic_HR_scaling_factor
+    expected_value = dynamic_HR_scaling_factor * dynamic_HR_scaling_factor
     assert np.allclose(ratio_in_sim, expected_value)
 
     # Check that expansion over two years with scaling prop to pop expansion works as expected
-    caps, final_popsize_increase = get_capabilities_after_two_years(
+    caps, final_popsize_increase = get_capabilities_after_two_updates(
         dynamic_HR_scaling_factor=1.0,
         scale_HR_by_pop_size=True
     )
@@ -2341,11 +2344,77 @@ def get_capabilities_after_two_years(dynamic_HR_scaling_factor: float, scale_HR_
     assert np.allclose(ratio_in_sim, expected_value)
 
     # Check that expansion over two years with both fixed scaling and pop expansion scaling works as expected
-    caps, final_popsize_increase = get_capabilities_after_two_years(
+    caps, final_popsize_increase = get_capabilities_after_two_updates(
         dynamic_HR_scaling_factor=dynamic_HR_scaling_factor,
         scale_HR_by_pop_size=True
     )
     caps = caps[caps != 0]
     ratio_in_sim = caps/initial_caps
     expected_value = final_popsize_increase*dynamic_HR_scaling_factor*dynamic_HR_scaling_factor
     assert np.allclose(ratio_in_sim, expected_value)
+
+
+def test_dynamic_HR_scaling_multiple_changes(seed, tmpdir):
+    """Check that we can scale the minutes of time available for healthcare workers with a sequence of factors that
+    apply in different years."""
+
+    def get_initial_capabilities() -> pd.Series:
+        sim = Simulation(start_date=start_date, seed=seed)
+        sim.register(
+            demography.Demography(resourcefilepath=resourcefilepath),
+            healthsystem.HealthSystem(resourcefilepath=resourcefilepath)
+        )
+        sim.make_initial_population(n=100)
+        sim.simulate(end_date=start_date + pd.DateOffset(days=0))
+
+        return sim.modules['HealthSystem'].capabilities_today
+
+    def run_sim(dynamic_HR_scaling_factor: Dict[int, float]) -> tuple:
+        """Run simulation for 10 years, with a sequence of factors that apply, specified in a dict of the form
+        {year: factor_to_apply_this_year_and_subsequent_years_until_next_instruction} (i.e. how the ResourceFile should
+        be structured.)
+        Returns capabilities at the end of the 10-year simulation"""
+
+        sim = Simulation(start_date=start_date, seed=seed)
+        sim.register(
+            demography.Demography(resourcefilepath=resourcefilepath),
+            healthsystem.HealthSystem(resourcefilepath=resourcefilepath),
+            simplified_births.SimplifiedBirths(resourcefilepath=resourcefilepath),
+
+        )
+        params = sim.modules['HealthSystem'].parameters
+        params['yearly_HR_scaling'][params['yearly_HR_scaling_mode']] = pd.DataFrame({
+            'year': dynamic_HR_scaling_factor.keys(),
+            'dynamic_HR_scaling_factor': dynamic_HR_scaling_factor.values(),
+            'scale_HR_by_popsize': False,
+        })
+
+        popsize = 100
+        sim.make_initial_population(n=popsize)
+
+        # Ensure simulation lasts long enough so that current capabilities reflect that used after two updates
+        # (updates occur on 1st Jan, starting in 2010, so simulation should stop on 2nd Jan 2011).
+        sim.simulate(end_date=sim.date + pd.DateOffset(years=10, days=1))
+
+        return sim.modules['HealthSystem'].capabilities_today
+
+    dynamic_HR_scaling_factor = {
+        2010: 1.0,
+        2011: 2.0,
+        # (2012 and 2013) are skipped: implies that the value for 2010 should apply in these years
+        2014: 0.2,
+        2015: 1.0,
+        # (2016, ..., 2020 are skipped: implies that the value for 2015 should apply in the years)
+    }
+    expected_overall_scaling = 2.0 * 2.0 * 2.0 * 0.2
+
+    # Get initial capabilities and remove all officers with no minutes available
+    initial_caps = get_initial_capabilities()
+    initial_caps = initial_caps[initial_caps != 0]
+
+    # Check that dynamic expansion over two years leads to expansion = dynamic_HR_scaling_factor^2
+    caps = run_sim(dynamic_HR_scaling_factor=dynamic_HR_scaling_factor)
+    caps = caps[caps != 0]
+    ratio_in_sim = caps / initial_caps
+
+    assert np.allclose(ratio_in_sim, expected_overall_scaling)