From baeafb342d5caf4bab175cd70dd90819b4570605 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Tue, 12 Oct 2021 03:38:06 +1100
Subject: [PATCH 01/33] Set up stochattic probability-based transitions to
 allow for measuring probability of e.g. elimination of a disease

---
 atomica/model.py   | 96 +++++++++++++++++++++++++++++++---------------
 atomica/project.py | 11 +++++-
 2 files changed, 76 insertions(+), 31 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index cc99fdc5..e45d927c 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -228,11 +228,12 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """ A class to wrap up data for one compartment within a cascade network. """
 
-    def __init__(self, pop, name):
+    def __init__(self, pop, name, stochastic:bool = False):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
         self.inlinks = []
+        self.stochastic = stochastic
         self._cached_outflow = None
 
     def unlink(self):
@@ -351,12 +352,24 @@ def resolve_outflows(self, ti: int) -> None:
             rescale = 1 / outflow
         else:
             rescale = 1
-
-        n = rescale * self.vals[ti]
+        
         self._cached_outflow = 0
-        for link in self.outlinks:
-            link.vals[ti] = link._cache * n
-            self._cached_outflow += link.vals[ti]
+        
+        if self.stochastic:
+            #for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
+            stays = max(0., 1. - outflow)
+            rescaled_frac_flows = [rescale * link._cache for link in self.outlinks] + [stays]
+            number_outflows = np.random.multinomial(self[ti], rescaled_frac_flows, size=1)[0]
+            for ln, link in enumerate(self.outlinks):
+                link.vals[ti] = number_outflows[ln]
+                self._cached_outflow += link.vals[ti]
+        else:
+            n = rescale * self.vals[ti]
+            for link in self.outlinks:
+                link.vals[ti] = link._cache * n
+                self._cached_outflow += link.vals[ti]
+        
+        
 
     def update(self, ti: int) -> None:
         """
@@ -398,7 +411,7 @@ def connect(self, dest, par) -> None:
 
 
 class JunctionCompartment(Compartment):
-    def __init__(self, pop, name: str, duration_group: str = None):
+    def __init__(self, pop, name: str, duration_group: str = None, stochastic:bool = False):
         """
 
         A TimedCompartment has a duration group by virtue of having a `.parameter` attribute and a flush link.
@@ -417,7 +430,7 @@ def __init__(self, pop, name: str, duration_group: str = None):
         :param duration_group: Optionally specify a duration group
 
         """
-        super().__init__(pop, name)
+        super().__init__(pop, name, stochastic=stochastic)
         self.duration_group = duration_group  #: Store the name of the duration group, if the junction belongs to one
 
     def connect(self, dest, par) -> None:
@@ -515,13 +528,19 @@ def balance(self, ti: int) -> None:
         # Next, get the total outflow. Note that the parameters are guaranteed to be in proportion units here
         outflow_fractions = [link.parameter.vals[ti] for link in self.outlinks]
         total_outflow = sum(outflow_fractions)
-
-        # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
+        
+        outflow_fractions /= total_outflow #proportionately accounting for the total outflow downscaling
+        
+        # assign outflow stochastically to compartments based on the probabilities
+        if self.stochastic and net_inflow > 0.:
+            outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
+
+        # Finally, assign the inflow to the outflow proportionately
         for frac, link in zip(outflow_fractions, self.outlinks):
             if self.duration_group:
-                link._vals[:, ti] = net_inflow * frac / total_outflow
+                link._vals[:, ti] = net_inflow * frac
             else:
-                link.vals[ti] = net_inflow * frac / total_outflow
+                link.vals[ti] = net_inflow * frac
 
     def initial_flush(self) -> None:
         """
@@ -541,6 +560,10 @@ def initial_flush(self) -> None:
             # Work out the outflow fractions
             outflow_fractions = np.array([link.parameter.vals[0] for link in self.outlinks])
             outflow_fractions /= np.sum(outflow_fractions)
+            
+            # assign outflow stochastically to compartments based on the probabilities
+            if self.stochastic:
+                outflow_fractions = np.random.multinomial(self.vals[0], outflow_fractions, size=1)[0]/self.vals[0]
 
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
@@ -593,6 +616,10 @@ def balance(self, ti: int) -> None:
             outflow_fractions /= total_outflow
             has_residual = False
 
+        # assign outflow stochastically to compartments based on the probabilities
+        if self.stochastic and net_inflow > 0.:
+            outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/(net_inflow * sum(outflow_fractions))
+            
         # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
         for frac, link in zip(outflow_fractions, self.outlinks):
             if link.parameter is None:
@@ -631,6 +658,10 @@ def initial_flush(self) -> None:
                 outflow_fractions /= total_outflow
                 has_residual = False
 
+            # assign outflow stochastically to compartments based on the probabilities
+            if self.stochastic and total_outflow > 0.:
+                outflow_fractions = np.random.multinomial(self.vals[0], outflow_fractions, size=1)[0]/(self.vals[0] * sum(outflow_fractions))
+                
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
             # to the downstream compartment type
@@ -661,8 +692,8 @@ class SourceCompartment(Compartment):
 
     """
 
-    def __init__(self, pop, name):
-        super().__init__(pop, name)
+    def __init__(self, pop, name, stochastic:bool = False):
+        super().__init__(pop, name, stochastic=stochastic)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         super().preallocate(tvec, dt)
@@ -677,8 +708,8 @@ def update(self, ti: int) -> None:
 
 
 class SinkCompartment(Compartment):
-    def __init__(self, pop, name):
-        super().__init__(pop, name)
+    def __init__(self, pop, name, stochastic:bool = False):
+        super().__init__(pop, name, stochastic=stochastic)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         """
@@ -734,7 +765,7 @@ def update(self, ti: int) -> None:
 
 
 class TimedCompartment(Compartment):
-    def __init__(self, pop, name: str, parameter: ParsetParameter):
+    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic:bool = False):
         """
         Instantiate the TimedCompartment
 
@@ -747,7 +778,7 @@ def __init__(self, pop, name: str, parameter: ParsetParameter):
 
         """
 
-        Compartment.__init__(self, pop=pop, name=name)
+        Compartment.__init__(self, pop=pop, name=name, stochastic=stochastic)
         self._vals = None  #: Primary storage, a matrix of size (duration x timesteps). The first row is the one that people leave from, the last row is the one that people arrive in
         self.parameter = parameter  #: The parameter to read the duration from - this needs to be done after parameters are precomputed though, in case the duration is coming from a (constant) function
         self.flush_link = None  #: Reference to the timed outflow link that flushes the compartment. Note that this needs to be set externally because Links are instantiated after Compartments
@@ -1533,7 +1564,7 @@ class Population:
     Each model population must contain a set of compartments with equivalent names.
     """
 
-    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str):
+    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool=False):
         """
         Construct a Population
 
@@ -1552,6 +1583,7 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         self.name = name  #: The code name of the population
         self.label = label  #: The full name/label of the population
         self.type = pop_type  #: The population's type
+        self.stochastic = stochastic #: Whether the population should be handled stochastically as discrete integer people instead of fractions
 
         self.comps = list()  #: List of Compartment objects
         self.characs = list()  #: List of Characteristic objects
@@ -1762,17 +1794,17 @@ def build(self, framework, progset):
         for comp_name in list(comps.index):
             if comps.at[comp_name, "population type"] == self.type:
                 if comp_name in residual_junctions:
-                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "is junction"] == "y":
-                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "duration group"]:
-                    self.comps.append(TimedCompartment(pop=self, name=comp_name, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
+                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic = self.stochastic, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
                 elif comps.at[comp_name, "is source"] == "y":
-                    self.comps.append(SourceCompartment(pop=self, name=comp_name))
+                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic = self.stochastic))
                 elif comps.at[comp_name, "is sink"] == "y":
-                    self.comps.append(SinkCompartment(pop=self, name=comp_name))
+                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic = self.stochastic))
                 else:
-                    self.comps.append(Compartment(pop=self, name=comp_name))
+                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic = self.stochastic))
 
         self.comp_lookup = {comp.name: comp for comp in self.comps}
 
@@ -1848,7 +1880,6 @@ def initialize_compartments(self, parset: ParameterSet, framework, t_init: float
         :param t_init: The year to use for initialization. This should generally be set to the sim start year
 
         """
-
         # Given a set of characteristics and their initial values, compute the initial
         # values for the compartments by solving the set of characteristics simultaneously
 
@@ -1955,7 +1986,11 @@ def report_characteristic(charac, n_indent=0):
 
         # Otherwise, insert the values
         for i, c in enumerate(comps):
-            c[0] = max(0.0, x[i])
+            if self.stochastic:
+                rem = np.random.binomial(1, x[i]-np.floor(x[i]), size=1)[0] if x[i] > np.floor(x[i]) else 0.
+                c[0] = np.floor((max(0.0, x[i]))) + rem #integer values (still represented as floats)
+            else:
+                c[0] = max(0.0, x[i])
 
 
 class Model:
@@ -1981,6 +2016,7 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.program_instructions = sc.dcp(program_instructions)  # program instructions
         self.t = settings.tvec  #: Simulation time vector (this is a brand new instance from the `settings.tvec` property method)
         self.dt = settings.sim_dt  #: Simulation time step
+        self.stochastic = settings.stochastic #: Run with discrete numbers of people per compartment instead of fractions.
 
         self._t_index = 0  # Keeps track of array index for current timepoint data within all compartments.
         self._vars_by_pop = None  # Cache to look up lists of variables by name across populations
@@ -2128,7 +2164,7 @@ def build(self, parset):
 
         # First construct populations
         for k, (pop_name, pop_label, pop_type) in enumerate(zip(parset.pop_names, parset.pop_labels, parset.pop_types)):
-            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type))
+            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, stochastic=self.stochastic))
             self._pop_ids[pop_name] = k
 
         # Expand interactions into matrix form
@@ -2402,7 +2438,7 @@ def update_links(self) -> None:
         """
 
         ti = self._t_index
-
+        
         # First, populate all of the link values without any outflow constraints
         for par in self._exec_order["transition_pars"]:
 
@@ -2442,7 +2478,7 @@ def update_links(self) -> None:
             elif par.units == FS.QUANTITY_TYPE_NUMBER:
                 # Disaggregate proportionally across all source compartment sizes related to all links.
                 converted_amt = transition * (self.dt / par.timescale)  # Number flow in this timestep, so it includes a timescale factor
-
+                
                 if isinstance(par.links[0].source, SourceCompartment):
                     # For a source compartment, the link value should be explicitly set directly
                     # Also, there is guaranteed to only be one link per parameter for outflows from source compartments
diff --git a/atomica/project.py b/atomica/project.py
index 54759ca7..e62898aa 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -44,10 +44,11 @@
 
 
 class ProjectSettings:
-    def __init__(self, sim_start=2000, sim_end=2035, sim_dt=0.25):
+    def __init__(self, sim_start=2000, sim_end=2035, sim_dt=0.25, stochastic=False):
         self._sim_start = sim_start
         self._sim_dt = sim_dt
         self._sim_end = 0.0
+        self._stochastic = stochastic
         self.update_time_vector(end=sim_end)
 
     def __repr__(self):
@@ -66,6 +67,10 @@ def sim_end(self):
     @property
     def sim_dt(self):
         return self._sim_dt
+    
+    @property
+    def stochastic(self):
+        return self._stochastic
 
     @sim_start.setter
     def sim_start(self, sim_start):
@@ -82,6 +87,10 @@ def sim_dt(self, sim_dt):
         assert sim_dt > 0, "Simulation timestep must be greater than 0"
         self._sim_dt = sim_dt
         self.sim_end = self.sim_end  # Call the setter function to change sim_end if it is no longer valid
+        
+    @stochastic.setter
+    def stochastic(self, stochastic:bool=False):
+        self._stochastic = stochastic
 
     @property
     def tvec(self) -> np.ndarray:

From 4037277ca2f9c3de25b380805955f4849b96ae91 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Tue, 12 Oct 2021 11:56:25 +1100
Subject: [PATCH 02/33] add stochastic_rounding

---
 atomica/model.py | 20 +++++++++++++++-----
 atomica/utils.py | 12 ++++++++++++
 2 files changed, 27 insertions(+), 5 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index e45d927c..b29e3d0f 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -23,6 +23,7 @@
 from .parameters import Parameter as ParsetParameter
 from .parameters import ParameterSet as ParameterSet
 import math
+from .utils import stochastic_rounding
 
 model_settings = dict()
 model_settings["tolerance"] = 1e-6
@@ -533,7 +534,10 @@ def balance(self, ti: int) -> None:
         
         # assign outflow stochastically to compartments based on the probabilities
         if self.stochastic and net_inflow > 0.:
-            outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
+            if self.duration_group:
+                raise Exception('Need to implement code here to make sure outflows from EACH part of the duration group are integers')
+            else:
+                outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
 
         # Finally, assign the inflow to the outflow proportionately
         for frac, link in zip(outflow_fractions, self.outlinks):
@@ -697,11 +701,15 @@ def __init__(self, pop, name, stochastic:bool = False):
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         super().preallocate(tvec, dt)
-        self.vals.fill(0.0)  #: Source compartments have an unlimited number of people in them. TODO: If this complicates validation, set to 0.0
+        self.vals.fill(0.0)  #: Source compartments have an unlimited number of people in them. Set to 0.0 to simplify validation.
 
     def resolve_outflows(self, ti: int) -> None:
+        #Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
         for link in self.outlinks:
-            link.vals[ti] = link._cache
+            if self.stochastic:
+                link.vals[ti] = stochastic_rounding(link._cache)
+            else:
+                link.vals[ti] = link._cache
 
     def update(self, ti: int) -> None:
         pass
@@ -881,6 +889,9 @@ def resolve_outflows(self, ti: int) -> None:
         :param ti: Time index at which to update link outflows
 
         """
+        
+        #TODO adjust for stochastic variables
+        assert self.stochastic == False, 'resolve_outflows needs to be updated to work with stochastic transitions'
 
         # First, work out the scale factors as usual
         self.flush_link._cache = 0.0  # At this stage, no outflow goes via the flush link
@@ -1987,8 +1998,7 @@ def report_characteristic(charac, n_indent=0):
         # Otherwise, insert the values
         for i, c in enumerate(comps):
             if self.stochastic:
-                rem = np.random.binomial(1, x[i]-np.floor(x[i]), size=1)[0] if x[i] > np.floor(x[i]) else 0.
-                c[0] = np.floor((max(0.0, x[i]))) + rem #integer values (still represented as floats)
+                c[0] = stochastic_rounding(max(0.0, x[i]))#integer values (still represented as floats)
             else:
                 c[0] = max(0.0, x[i])
 
diff --git a/atomica/utils.py b/atomica/utils.py
index b440bcdd..2617fef0 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -1037,3 +1037,15 @@ def stop_logging() -> None:
             logger.removeHandler(handler)
             # Don't terminate the loop, if by some change there is more than one handler
             # (not supposed to happen though) then we would want to close them all
+
+def stochastic_rounding(x):
+    """
+    Stochastically round a float up or down to an integer-equivalent value (note: returns still as a float)
+    :param x: value to be rounded
+    """
+    floor = np.floor(x)
+    remainder = x - floor
+    if remainder:
+        x = floor + int(np.random.random()<remainder)
+    return x
+    
\ No newline at end of file

From 425223daa015164fa57ce12afba4b36f2d222a45 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Sat, 6 Nov 2021 03:28:56 +1100
Subject: [PATCH 03/33] Random seed for parset sampling, random.randn ->
 default_rng sampling

---
 atomica/parameters.py | 21 +++++++++++++++------
 atomica/project.py    | 23 +++++++++++++++--------
 atomica/utils.py      | 10 +++++++---
 3 files changed, 37 insertions(+), 17 deletions(-)

diff --git a/atomica/parameters.py b/atomica/parameters.py
index de35e9ed..f9c137d2 100644
--- a/atomica/parameters.py
+++ b/atomica/parameters.py
@@ -90,7 +90,7 @@ def interpolate(self, tvec, pop_name: str) -> np.array:
 
         return self.ts[pop_name].interpolate(tvec, method=self._interpolation_method)
 
-    def sample(self, constant: bool) -> None:
+    def sample(self, constant: bool, rand_seed: int = None) -> None:
         """
         Perturb parameter based on uncertainties
 
@@ -99,11 +99,15 @@ def sample(self, constant: bool) -> None:
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
-
+        :param rand_seed: To generate a consistent "random" sample
         """
+        if rand_seed is not None:
+            rng = np.random.default_rng(seed = rand_seed)
 
         for k, ts in self.ts.items():
-            self.ts[k] = ts.sample(constant)
+            if rand_seed is not None:
+                rand_seed = rng.integers(1e15)
+            self.ts[k] = ts.sample(constant, rand_seed)
 
     def smooth(self, tvec, method="smoothinterp", pop_names=None, **kwargs):
         """
@@ -280,7 +284,7 @@ def get_par(self, name: str, pop: str = None) -> Parameter:
         else:
             raise KeyError(f'Parameter "{name}" not found')
 
-    def sample(self, constant=True):
+    def sample(self, constant:bool = True, rand_seed:int = None):
         """
         Return a sampled copy of the ParameterSet
 
@@ -291,8 +295,13 @@ def sample(self, constant=True):
         """
 
         new = sc.dcp(self)
-        for par in new.all_pars():
-            par.sample(constant)
+        if rand_seed is not None:
+            rng = np.random.default_rng(seed = rand_seed)
+
+        for i, par in enumerate(new.all_pars()):
+            if rand_seed is not None:
+                rand_seed = rng.integers(1e15) 
+            par.sample(constant = constant, rand_seed = rand_seed)
         return new
 
     ### SAVE AND LOAD CALIBRATION (Y-FACTORS)
diff --git a/atomica/project.py b/atomica/project.py
index e62898aa..f8f50735 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -516,7 +516,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
 
         return result
 
-    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, parallel=False, max_attempts=None, num_workers=None) -> list:
+    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, rand_seed: int = None, parallel=False, max_attempts=None, num_workers=None) -> list:
         """
         Run sampled simulations
 
@@ -540,6 +540,7 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         :param parallel: If True, run simulations in parallel (on Windows, must have ``if __name__ == '__main__'`` gating the calling code)
         :param max_attempts: Number of retry attempts for bad initializations
         :param num_workers: If ``parallel`` is True, this determines the number of parallel workers to use (default is usually number of CPUs)
+        :param rand_seed: random seed to provide consistent results
         :return: A list of Results that can be passed to `Ensemble.update()`. If multiple instructions are provided, the return value of this
                  function will be a list of lists, where the inner list iterates over different instructions for the same parset/progset samples.
                  It is expected in that case that the Ensemble's mapping function would take in a list of results
@@ -562,18 +563,24 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
             assert (len(result_names) == 1 and not progset) or (len(progset_instructions) == len(result_names)), "Number of result names must match number of instructions"
 
         show_progress = n_samples > 1 and logger.getEffectiveLevel() <= logging.INFO
+        
+        if rand_seed is not None:
+            rng = np.random.default_rng(seed = rand_seed)
+            seed_samples = rng.integers(1e15, size=n_samples)
+        else:
+            seed_samples = [None]*n_samples
 
         if parallel:
             fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
-            results = parallel_progress(fcn, n_samples, show_progress=show_progress, num_workers=num_workers)
+            results = parallel_progress(fcn, seed_samples, show_progress=show_progress, num_workers=num_workers)
         elif show_progress:
             # Print the progress bar if the logging level was INFO or lower
             # This means that the user can still set the logging level higher e.g. WARNING to suppress output from Atomica in general
             # (including any progress bars)
             with Quiet():
-                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts) for _ in tqdm.trange(n_samples)]
+                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rand_seed=seed) for seed in tqdm.tqdm(seed_samples)]
         else:
-            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts) for _ in range(n_samples)]
+            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rand_seed=seed) for seed in seed_samples]
 
         return results
 
@@ -714,8 +721,7 @@ def __setstate__(self, d):
         P = migrate(self)
         self.__dict__ = P.__dict__
 
-
-def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None):
+def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rand_seed: int =  None):
     """
     Internal function to run simulation with sampling
 
@@ -736,6 +742,7 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     :param progset_instructions: A list of instructions to run against a single sample
     :param result_names: A list of result names (strings)
     :param max_attempts: Maximum number of sampling attempts before raising an error
+    :param rand_seed: A random seed used to get a consistent sampled parset
     :return: A list of results that either contains 1 result, or the same number of results as instructions
 
     """
@@ -748,14 +755,14 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     attempts = 0
     while attempts < max_attempts:
         try:
+            sampled_parset = parset.sample(rand_seed)
             if progset:
-                sampled_parset = parset.sample()
                 sampled_progset = progset.sample()
                 results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y) for x, y in zip(progset_instructions, result_names)]
             else:
-                sampled_parset = parset.sample()
                 results = [proj.run_sim(parset=sampled_parset, result_name=y) for y in result_names]
             return results
         except BadInitialization:
             attempts += 1
+            rand_seed += 1 #this will avoid some bad initializations due to parset sampling while remaining consistent
     raise Exception("Failed simulation after %d attempts - something might have gone wrong" % (max_attempts))
diff --git a/atomica/utils.py b/atomica/utils.py
index 2617fef0..41001a9b 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -588,7 +588,7 @@ def interpolate(self, t2: np.array, method="linear", **kwargs) -> np.array:
         interpolator = method(t1, v1, **kwargs)
         return interpolator(t2)
 
-    def sample(self, constant=True):
+    def sample(self, constant:bool = True, rand_seed:int = None):
         """
         Return a sampled copy of the TimeSeries
 
@@ -597,16 +597,20 @@ def sample(self, constant=True):
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
+        :param rand_seed: To generate a consistent "random" sample
+                         
         :return: A copied ``TimeSeries`` with perturbed values
 
         """
 
         if self._sampled:
             raise Exception("Sampling has already been performed - can only sample once")
+            
+        rng = np.random.default_rng(seed = rand_seed)
 
         new = self.copy()
         if self.sigma is not None:
-            delta = self.sigma * np.random.randn(1)[0]
+            delta = self.sigma * rng.standard_normal(1)[0]
             if self.assumption is not None:
                 new.assumption += delta
 
@@ -615,7 +619,7 @@ def sample(self, constant=True):
                 new.vals = [v + delta for v in new.vals]
             else:
                 # Sample again for each data point
-                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * np.random.randn(len(new.vals)))):
+                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * rng.standard_normal(len(new.vals)))):
                     new.vals[i] = v + delta
 
         # Sampling flag only needs to be set if the TimeSeries had data to change

From e3c628520423d8787443f6927af140e787f0780a Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Sat, 6 Nov 2021 03:58:48 +1100
Subject: [PATCH 04/33] default_rng for programs/covouts also

---
 atomica/programs.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/atomica/programs.py b/atomica/programs.py
index 6325d759..5bb61cf3 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1367,25 +1367,29 @@ def n_progs(self) -> int:
 
         return len(self.progs)
 
-    def sample(self) -> None:
+    def sample(self, rand_seed:int = None) -> None:
         """
         Perturb the values entered in the databook
 
         The :class:`Covout` instance is modified in-place. Note that the program outcomes are scalars
         that do not vary over time - therefore, :meth:`Covout.sample()` does not have a ``constant``
         argument.
+        
+        :param rand_seed: used to generate consistent random results
 
         """
 
         if self.sigma is None:
             return
+        
+        rng = np.random.default_rng(seed=rand_seed)
 
         for k, v in self.progs.items():
-            self.progs[k] = v + self.sigma * np.random.randn(1)[0]
+            self.progs[k] = v + self.sigma * rng.standard_normal(1)[0]
         # Perturb the interactions
         if self._interactions:
             for k, v in self.interactions.items():
-                self.interactions[k] = v + self.sigma * np.random.randn(1)[0]
+                self.interactions[k] = v + self.sigma * rng.standard_normal(1)[0]
             tokens = ["%s=%.4f" % ("+".join(k), v) for k, v in self.interactions.items()]
             self.imp_interaction = ",".join(tokens)
 

From 6cbcb6f1bcc7bb04059b656d9db7dafe6d933b49 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Sat, 6 Nov 2021 04:41:14 +1100
Subject: [PATCH 05/33] Update project.py

---
 atomica/project.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/atomica/project.py b/atomica/project.py
index f8f50735..1ab2ac66 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -756,6 +756,7 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     while attempts < max_attempts:
         try:
             sampled_parset = parset.sample(rand_seed)
+            sampled_parset = parset.sample(rand_seed = rand_seed)
             if progset:
                 sampled_progset = progset.sample()
                 results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y) for x, y in zip(progset_instructions, result_names)]

From 752154a5da3c64258acddae440ce40a74bc78ff6 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Sat, 6 Nov 2021 04:41:17 +1100
Subject: [PATCH 06/33] Update project.py

---
 atomica/project.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/atomica/project.py b/atomica/project.py
index 1ab2ac66..6e7a6a56 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -755,7 +755,6 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     attempts = 0
     while attempts < max_attempts:
         try:
-            sampled_parset = parset.sample(rand_seed)
             sampled_parset = parset.sample(rand_seed = rand_seed)
             if progset:
                 sampled_progset = progset.sample()

From 396527603fdb7f28ab89061dd1fb2cd7af0f7960 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 01:33:07 +1100
Subject: [PATCH 07/33] Pass the same random number generator within a given
 model run rather than creating local ones

---
 atomica/model.py      | 67 ++++++++++++++++++++++++++-----------------
 atomica/parameters.py | 21 +++++---------
 atomica/programs.py   | 18 ++++++------
 atomica/project.py    | 29 ++++++++++++-------
 atomica/utils.py      | 19 ++++++------
 5 files changed, 85 insertions(+), 69 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index b29e3d0f..81aeae9d 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -229,12 +229,15 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """ A class to wrap up data for one compartment within a cascade network. """
 
-    def __init__(self, pop, name, stochastic:bool = False):
+    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
         self.inlinks = []
         self.stochastic = stochastic
+        self.rng_sampler = rng_sampler
+        if self.stochastic and rng_sampler is None:
+            self.rng_sampler = np.random.default_rng()
         self._cached_outflow = None
 
     def unlink(self):
@@ -360,7 +363,7 @@ def resolve_outflows(self, ti: int) -> None:
             #for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
             stays = max(0., 1. - outflow)
             rescaled_frac_flows = [rescale * link._cache for link in self.outlinks] + [stays]
-            number_outflows = np.random.multinomial(self[ti], rescaled_frac_flows, size=1)[0]
+            number_outflows = self.rng_sampler.multinomial(self[ti], rescaled_frac_flows, size=1)[0]
             for ln, link in enumerate(self.outlinks):
                 link.vals[ti] = number_outflows[ln]
                 self._cached_outflow += link.vals[ti]
@@ -412,7 +415,7 @@ def connect(self, dest, par) -> None:
 
 
 class JunctionCompartment(Compartment):
-    def __init__(self, pop, name: str, duration_group: str = None, stochastic:bool = False):
+    def __init__(self, pop, name: str, duration_group: str = None, stochastic:bool = False, rng_sampler = None):
         """
 
         A TimedCompartment has a duration group by virtue of having a `.parameter` attribute and a flush link.
@@ -431,7 +434,7 @@ def __init__(self, pop, name: str, duration_group: str = None, stochastic:bool =
         :param duration_group: Optionally specify a duration group
 
         """
-        super().__init__(pop, name, stochastic=stochastic)
+        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
         self.duration_group = duration_group  #: Store the name of the duration group, if the junction belongs to one
 
     def connect(self, dest, par) -> None:
@@ -537,7 +540,7 @@ def balance(self, ti: int) -> None:
             if self.duration_group:
                 raise Exception('Need to implement code here to make sure outflows from EACH part of the duration group are integers')
             else:
-                outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
+                outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
 
         # Finally, assign the inflow to the outflow proportionately
         for frac, link in zip(outflow_fractions, self.outlinks):
@@ -567,7 +570,7 @@ def initial_flush(self) -> None:
             
             # assign outflow stochastically to compartments based on the probabilities
             if self.stochastic:
-                outflow_fractions = np.random.multinomial(self.vals[0], outflow_fractions, size=1)[0]/self.vals[0]
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0]/self.vals[0]
 
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
@@ -622,7 +625,7 @@ def balance(self, ti: int) -> None:
 
         # assign outflow stochastically to compartments based on the probabilities
         if self.stochastic and net_inflow > 0.:
-            outflow_fractions = np.random.multinomial(net_inflow, outflow_fractions, size=1)[0]/(net_inflow * sum(outflow_fractions))
+            outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0]/(net_inflow * sum(outflow_fractions))
             
         # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
         for frac, link in zip(outflow_fractions, self.outlinks):
@@ -664,7 +667,7 @@ def initial_flush(self) -> None:
 
             # assign outflow stochastically to compartments based on the probabilities
             if self.stochastic and total_outflow > 0.:
-                outflow_fractions = np.random.multinomial(self.vals[0], outflow_fractions, size=1)[0]/(self.vals[0] * sum(outflow_fractions))
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0]/(self.vals[0] * sum(outflow_fractions))
                 
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
@@ -696,8 +699,8 @@ class SourceCompartment(Compartment):
 
     """
 
-    def __init__(self, pop, name, stochastic:bool = False):
-        super().__init__(pop, name, stochastic=stochastic)
+    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
+        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         super().preallocate(tvec, dt)
@@ -707,7 +710,7 @@ def resolve_outflows(self, ti: int) -> None:
         #Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
         for link in self.outlinks:
             if self.stochastic:
-                link.vals[ti] = stochastic_rounding(link._cache)
+                link.vals[ti] = stochastic_rounding(link._cache, rng_sampler = self.rng_sampler)
             else:
                 link.vals[ti] = link._cache
 
@@ -716,8 +719,8 @@ def update(self, ti: int) -> None:
 
 
 class SinkCompartment(Compartment):
-    def __init__(self, pop, name, stochastic:bool = False):
-        super().__init__(pop, name, stochastic=stochastic)
+    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
+        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         """
@@ -773,7 +776,7 @@ def update(self, ti: int) -> None:
 
 
 class TimedCompartment(Compartment):
-    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic:bool = False):
+    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic:bool = False, rng_sampler= None):
         """
         Instantiate the TimedCompartment
 
@@ -786,7 +789,7 @@ def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic:bool =
 
         """
 
-        Compartment.__init__(self, pop=pop, name=name, stochastic=stochastic)
+        Compartment.__init__(self, pop=pop, name=name, stochastic=stochastic, rng_sampler=rng_sampler)
         self._vals = None  #: Primary storage, a matrix of size (duration x timesteps). The first row is the one that people leave from, the last row is the one that people arrive in
         self.parameter = parameter  #: The parameter to read the duration from - this needs to be done after parameters are precomputed though, in case the duration is coming from a (constant) function
         self.flush_link = None  #: Reference to the timed outflow link that flushes the compartment. Note that this needs to be set externally because Links are instantiated after Compartments
@@ -1575,7 +1578,7 @@ class Population:
     Each model population must contain a set of compartments with equivalent names.
     """
 
-    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool=False):
+    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool=False, rng_sampler=None):
         """
         Construct a Population
 
@@ -1588,6 +1591,9 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         :param name: The code name of the population
         :param label: The full name of the population
         :param progset: A ProgramSet instance
+        
+        :param stochastic: whether to use stochasticity in determining population movement
+        :param rng_sampler: optional Generator for random numbers to give consistent results between model runs
 
         """
 
@@ -1595,6 +1601,7 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         self.label = label  #: The full name/label of the population
         self.type = pop_type  #: The population's type
         self.stochastic = stochastic #: Whether the population should be handled stochastically as discrete integer people instead of fractions
+        self.rng_sampler = rng_sampler
 
         self.comps = list()  #: List of Compartment objects
         self.characs = list()  #: List of Characteristic objects
@@ -1805,17 +1812,17 @@ def build(self, framework, progset):
         for comp_name in list(comps.index):
             if comps.at[comp_name, "population type"] == self.type:
                 if comp_name in residual_junctions:
-                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "is junction"] == "y":
-                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "duration group"]:
-                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic = self.stochastic, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
+                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
                 elif comps.at[comp_name, "is source"] == "y":
-                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic = self.stochastic))
+                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
                 elif comps.at[comp_name, "is sink"] == "y":
-                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic = self.stochastic))
+                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
                 else:
-                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic = self.stochastic))
+                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
 
         self.comp_lookup = {comp.name: comp for comp in self.comps}
 
@@ -1998,7 +2005,7 @@ def report_characteristic(charac, n_indent=0):
         # Otherwise, insert the values
         for i, c in enumerate(comps):
             if self.stochastic:
-                c[0] = stochastic_rounding(max(0.0, x[i]))#integer values (still represented as floats)
+                c[0] = stochastic_rounding(max(0.0, x[i]), rng_sampler = self.rng_sampler)#integer values (still represented as floats)
             else:
                 c[0] = max(0.0, x[i])
 
@@ -2006,7 +2013,7 @@ def report_characteristic(charac, n_indent=0):
 class Model:
     """ A class to wrap up multiple populations within model and handle cross-population transitions. """
 
-    def __init__(self, settings, framework, parset, progset=None, program_instructions=None):
+    def __init__(self, settings, framework, parset, progset=None, program_instructions=None, rng_sampler=None):
 
         # Note that if a progset is provided and program instructions are not, then programs will not be
         # turned on. However, the progset is still available so that the coverage can still be probed
@@ -2027,6 +2034,11 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.t = settings.tvec  #: Simulation time vector (this is a brand new instance from the `settings.tvec` property method)
         self.dt = settings.sim_dt  #: Simulation time step
         self.stochastic = settings.stochastic #: Run with discrete numbers of people per compartment instead of fractions.
+        
+
+        if rng_sampler is None:
+            rng_sampler = np.random.default_rng()
+        self.rng_sampler = rng_sampler
 
         self._t_index = 0  # Keeps track of array index for current timepoint data within all compartments.
         self._vars_by_pop = None  # Cache to look up lists of variables by name across populations
@@ -2174,7 +2186,7 @@ def build(self, parset):
 
         # First construct populations
         for k, (pop_name, pop_label, pop_type) in enumerate(zip(parset.pop_names, parset.pop_labels, parset.pop_types)):
-            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, stochastic=self.stochastic))
+            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
             self._pop_ids[pop_name] = k
 
         # Expand interactions into matrix form
@@ -2677,7 +2689,7 @@ def update_pars(self) -> None:
                     par.constrain(ti)
 
 
-def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None):
+def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler = None):
     """
     Build and process model
 
@@ -2701,10 +2713,11 @@ def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = N
     :param progset: Optionally provide a :class:`ProgramSet` instance to use programs
     :param program_instructions: Optional :class:`ProgramInstructions` instance. If ``progset`` is specified, then instructions must be provided
     :param name: Optionally specify the name to assign to the output result
+    :param rng_sampler: Optionally specify a random number generator that may have been seeded to generate consistent results
     :return: A :class:`Result` object containing the processed model
 
     """
 
-    m = Model(settings, framework, parset, progset, program_instructions)
+    m = Model(settings, framework, parset, progset, program_instructions, rng_sampler=rng_sampler)
     m.process()
     return Result(model=m, parset=parset, name=name)
diff --git a/atomica/parameters.py b/atomica/parameters.py
index f9c137d2..8c643ba3 100644
--- a/atomica/parameters.py
+++ b/atomica/parameters.py
@@ -90,7 +90,7 @@ def interpolate(self, tvec, pop_name: str) -> np.array:
 
         return self.ts[pop_name].interpolate(tvec, method=self._interpolation_method)
 
-    def sample(self, constant: bool, rand_seed: int = None) -> None:
+    def sample(self, constant: bool, rng_sampler = None) -> None:
         """
         Perturb parameter based on uncertainties
 
@@ -99,15 +99,11 @@ def sample(self, constant: bool, rand_seed: int = None) -> None:
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
-        :param rand_seed: To generate a consistent "random" sample
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
-        if rand_seed is not None:
-            rng = np.random.default_rng(seed = rand_seed)
 
         for k, ts in self.ts.items():
-            if rand_seed is not None:
-                rand_seed = rng.integers(1e15)
-            self.ts[k] = ts.sample(constant, rand_seed)
+            self.ts[k] = ts.sample(constant, rng_sampler)
 
     def smooth(self, tvec, method="smoothinterp", pop_names=None, **kwargs):
         """
@@ -284,24 +280,21 @@ def get_par(self, name: str, pop: str = None) -> Parameter:
         else:
             raise KeyError(f'Parameter "{name}" not found')
 
-    def sample(self, constant:bool = True, rand_seed:int = None):
+    def sample(self, constant:bool = True, rng_sampler = None):
         """
         Return a sampled copy of the ParameterSet
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+        
         :return: A new :class:`ParameterSet` with perturbed values
 
         """
-
         new = sc.dcp(self)
-        if rand_seed is not None:
-            rng = np.random.default_rng(seed = rand_seed)
 
         for i, par in enumerate(new.all_pars()):
-            if rand_seed is not None:
-                rand_seed = rng.integers(1e15) 
-            par.sample(constant = constant, rand_seed = rand_seed)
+            par.sample(constant = constant, rng_sampler = rng_sampler)
         return new
 
     ### SAVE AND LOAD CALIBRATION (Y-FACTORS)
diff --git a/atomica/programs.py b/atomica/programs.py
index 5bb61cf3..3e3d709b 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1191,7 +1191,7 @@ def is_one_off(self) -> bool:
         """
         return "/year" not in self.unit_cost.units
 
-    def sample(self, constant: bool) -> None:
+    def sample(self, constant: bool, rng_sampler = None) -> None:
         """
         Perturb program values based on uncertainties
 
@@ -1201,11 +1201,12 @@ def sample(self, constant: bool) -> None:
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
 
         self.spend_data = self.spend_data.sample(constant)
         self.unit_cost = self.unit_cost.sample(constant)
-        self.capacity_constraint = self.capacity_constraint.sample(constant)
+        self.capacity_constraint = self.capacity_constraint.sample(constant, rng_sampler)
         self.saturation = self.saturation.sample(constant)
         self.coverage = self.coverage.sample(constant)
 
@@ -1367,29 +1368,28 @@ def n_progs(self) -> int:
 
         return len(self.progs)
 
-    def sample(self, rand_seed:int = None) -> None:
+    def sample(self, rng_sampler = None) -> None:
         """
         Perturb the values entered in the databook
 
         The :class:`Covout` instance is modified in-place. Note that the program outcomes are scalars
         that do not vary over time - therefore, :meth:`Covout.sample()` does not have a ``constant``
         argument.
-        
-        :param rand_seed: used to generate consistent random results
-
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         """
 
         if self.sigma is None:
             return
         
-        rng = np.random.default_rng(seed=rand_seed)
+        if rng_sampler is None:
+            rng_sampler = np.random.default_rng()
 
         for k, v in self.progs.items():
-            self.progs[k] = v + self.sigma * rng.standard_normal(1)[0]
+            self.progs[k] = v + self.sigma * rng_sampler.standard_normal(1)[0]
         # Perturb the interactions
         if self._interactions:
             for k, v in self.interactions.items():
-                self.interactions[k] = v + self.sigma * rng.standard_normal(1)[0]
+                self.interactions[k] = v + self.sigma * rng_sampler.standard_normal(1)[0]
             tokens = ["%s=%.4f" % ("+".join(k), v) for k, v in self.interactions.items()]
             self.imp_interaction = ",".join(tokens)
 
diff --git a/atomica/project.py b/atomica/project.py
index 6e7a6a56..897d6779 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -475,7 +475,7 @@ def update_settings(self, sim_start=None, sim_end=None, sim_dt=None):
         """ Modify the project settings, e.g. the simulation time vector. """
         self.settings.update_time_vector(start=sim_start, end=sim_end, dt=sim_dt)
 
-    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None):
+    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng_sampler = None):
         """
         Run a single simulation
 
@@ -488,6 +488,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
         :param progset_instructions: A :class:`ProgramInstructions` instance. Programs will only be used if a instructions are provided
         :param store_results: If True, then the result will automatically be stored in ``self.results``
         :param result_name: Optionally assign a specific name to the result (otherwise, a unique default name will automatically be selected)
+        :param rng_sampler: A random number generator that may have been seeded to generate consistent results
         :return: A :class:`Result` instance
 
         """
@@ -509,7 +510,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
                 k += 1
 
         tm = sc.tic()
-        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name)
+        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name, rng_sampler=rng_sampler)
         logger.info('Elapsed time for running "%s": %ss', self.name, sc.sigfig(sc.toc(tm, output=True), 3))
         if store_results:
             self.results.append(result)
@@ -540,7 +541,7 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         :param parallel: If True, run simulations in parallel (on Windows, must have ``if __name__ == '__main__'`` gating the calling code)
         :param max_attempts: Number of retry attempts for bad initializations
         :param num_workers: If ``parallel`` is True, this determines the number of parallel workers to use (default is usually number of CPUs)
-        :param rand_seed: random seed to provide consistent results
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
         :return: A list of Results that can be passed to `Ensemble.update()`. If multiple instructions are provided, the return value of this
                  function will be a list of lists, where the inner list iterates over different instructions for the same parset/progset samples.
                  It is expected in that case that the Ensemble's mapping function would take in a list of results
@@ -569,18 +570,22 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
             seed_samples = rng.integers(1e15, size=n_samples)
         else:
             seed_samples = [None]*n_samples
+            
+        model_rngs = [np.random.default_rng(seed = seed) for seed in seed_samples] #generate a RNG for each model
+
+        print (seed_samples)
 
         if parallel:
             fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
-            results = parallel_progress(fcn, seed_samples, show_progress=show_progress, num_workers=num_workers)
+            results = parallel_progress(fcn, model_rngs, show_progress=show_progress, num_workers=num_workers)
         elif show_progress:
             # Print the progress bar if the logging level was INFO or lower
             # This means that the user can still set the logging level higher e.g. WARNING to suppress output from Atomica in general
             # (including any progress bars)
             with Quiet():
-                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rand_seed=seed) for seed in tqdm.tqdm(seed_samples)]
+                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng) for rng in tqdm.tqdm(model_rngs)]
         else:
-            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rand_seed=seed) for seed in seed_samples]
+            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng) for rng in model_rngs]
 
         return results
 
@@ -721,7 +726,7 @@ def __setstate__(self, d):
         P = migrate(self)
         self.__dict__ = P.__dict__
 
-def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rand_seed: int =  None):
+def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler =  None):
     """
     Internal function to run simulation with sampling
 
@@ -742,7 +747,7 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     :param progset_instructions: A list of instructions to run against a single sample
     :param result_names: A list of result names (strings)
     :param max_attempts: Maximum number of sampling attempts before raising an error
-    :param rand_seed: A random seed used to get a consistent sampled parset
+    :param rng_sampler: A Generator used to 
     :return: A list of results that either contains 1 result, or the same number of results as instructions
 
     """
@@ -751,18 +756,20 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
 
     if max_attempts is None:
         max_attempts = 50
+    
+    if rng_sampler is None:
+        rng_sampler = np.random.default_rng()
 
     attempts = 0
     while attempts < max_attempts:
         try:
-            sampled_parset = parset.sample(rand_seed = rand_seed)
+            sampled_parset = parset.sample(rng_sampler=rng_sampler)
             if progset:
-                sampled_progset = progset.sample()
+                sampled_progset = progset.sample(rng_sampler=rng_sampler)
                 results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y) for x, y in zip(progset_instructions, result_names)]
             else:
                 results = [proj.run_sim(parset=sampled_parset, result_name=y) for y in result_names]
             return results
         except BadInitialization:
             attempts += 1
-            rand_seed += 1 #this will avoid some bad initializations due to parset sampling while remaining consistent
     raise Exception("Failed simulation after %d attempts - something might have gone wrong" % (max_attempts))
diff --git a/atomica/utils.py b/atomica/utils.py
index 41001a9b..023ffb2e 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -588,7 +588,7 @@ def interpolate(self, t2: np.array, method="linear", **kwargs) -> np.array:
         interpolator = method(t1, v1, **kwargs)
         return interpolator(t2)
 
-    def sample(self, constant:bool = True, rand_seed:int = None):
+    def sample(self, constant:bool = True, rng_sampler = None):
         """
         Return a sampled copy of the TimeSeries
 
@@ -597,7 +597,7 @@ def sample(self, constant:bool = True, rand_seed:int = None):
 
         :param constant: If True, time series will be perturbed by a single constant offset. If False,
                          an different perturbation will be applied to each time specific value independently.
-        :param rand_seed: To generate a consistent "random" sample
+        :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
                          
         :return: A copied ``TimeSeries`` with perturbed values
 
@@ -605,12 +605,13 @@ def sample(self, constant:bool = True, rand_seed:int = None):
 
         if self._sampled:
             raise Exception("Sampling has already been performed - can only sample once")
-            
-        rng = np.random.default_rng(seed = rand_seed)
+        
+        if rng_sampler is None:
+            rng_sampler = np.random.default_rng()
 
         new = self.copy()
         if self.sigma is not None:
-            delta = self.sigma * rng.standard_normal(1)[0]
+            delta = self.sigma * rng_sampler.standard_normal(1)[0]
             if self.assumption is not None:
                 new.assumption += delta
 
@@ -619,7 +620,7 @@ def sample(self, constant:bool = True, rand_seed:int = None):
                 new.vals = [v + delta for v in new.vals]
             else:
                 # Sample again for each data point
-                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * rng.standard_normal(len(new.vals)))):
+                for i, (v, delta) in enumerate(zip(new.vals, self.sigma * rng_sampler.standard_normal(len(new.vals)))):
                     new.vals[i] = v + delta
 
         # Sampling flag only needs to be set if the TimeSeries had data to change
@@ -1042,14 +1043,16 @@ def stop_logging() -> None:
             # Don't terminate the loop, if by some change there is more than one handler
             # (not supposed to happen though) then we would want to close them all
 
-def stochastic_rounding(x):
+def stochastic_rounding(x, rng_sampler = None):
     """
     Stochastically round a float up or down to an integer-equivalent value (note: returns still as a float)
     :param x: value to be rounded
     """
+    sample = np.random.random() if rng_sampler is None else rng_sampler.random() 
+    
     floor = np.floor(x)
     remainder = x - floor
     if remainder:
-        x = floor + int(np.random.random()<remainder)
+        x = floor + int(sample<remainder)
     return x
     
\ No newline at end of file

From ac4e45851d96f08b545adefae9013b5199e94ee7 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 01:54:24 +1100
Subject: [PATCH 08/33] Consistent seeding

---
 atomica/project.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/atomica/project.py b/atomica/project.py
index 897d6779..0beb6d0a 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -475,7 +475,7 @@ def update_settings(self, sim_start=None, sim_end=None, sim_dt=None):
         """ Modify the project settings, e.g. the simulation time vector. """
         self.settings.update_time_vector(start=sim_start, end=sim_end, dt=sim_dt)
 
-    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng_sampler = None):
+    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng = None):
         """
         Run a single simulation
 
@@ -488,10 +488,11 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
         :param progset_instructions: A :class:`ProgramInstructions` instance. Programs will only be used if a instructions are provided
         :param store_results: If True, then the result will automatically be stored in ``self.results``
         :param result_name: Optionally assign a specific name to the result (otherwise, a unique default name will automatically be selected)
-        :param rng_sampler: A random number generator that may have been seeded to generate consistent results
+        :param rng: Optionally a random number generator that may have been seeded to generate consistent results, or a random_seed used to generate a Generator
         :return: A :class:`Result` instance
 
         """
+        rng_sampler = rng if isinstance(rng, np.random._generator.Generator) else np.random.default_rng(rng)
 
         parset = self.parset(parset)
         if progset is not None:
@@ -573,8 +574,6 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
             
         model_rngs = [np.random.default_rng(seed = seed) for seed in seed_samples] #generate a RNG for each model
 
-        print (seed_samples)
-
         if parallel:
             fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
             results = parallel_progress(fcn, model_rngs, show_progress=show_progress, num_workers=num_workers)
@@ -747,7 +746,7 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     :param progset_instructions: A list of instructions to run against a single sample
     :param result_names: A list of result names (strings)
     :param max_attempts: Maximum number of sampling attempts before raising an error
-    :param rng_sampler: A Generator used to 
+    :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
     :return: A list of results that either contains 1 result, or the same number of results as instructions
 
     """
@@ -766,9 +765,9 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
             sampled_parset = parset.sample(rng_sampler=rng_sampler)
             if progset:
                 sampled_progset = progset.sample(rng_sampler=rng_sampler)
-                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y) for x, y in zip(progset_instructions, result_names)]
+                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_sampler) for x, y in zip(progset_instructions, result_names)]
             else:
-                results = [proj.run_sim(parset=sampled_parset, result_name=y) for y in result_names]
+                results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_sampler) for y in result_names]
             return results
         except BadInitialization:
             attempts += 1

From 1e7c94fbbdd9132b46f2f5c3393ebbd1e6c3e6ba Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 03:56:09 +1100
Subject: [PATCH 09/33] Pass generators within kwargs, and handle that option
 in parallel_progress

---
 atomica/project.py | 4 +++-
 atomica/utils.py   | 5 ++++-
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/atomica/project.py b/atomica/project.py
index 0beb6d0a..f783a8d6 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -576,7 +576,9 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
 
         if parallel:
             fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
-            results = parallel_progress(fcn, model_rngs, show_progress=show_progress, num_workers=num_workers)
+            #as multiprocessing does not handle partial functions as compiled functions, need to send the rngs as kwargs in a dictionary, not as args to the partial function
+            model_rng_kwargs = [{'rng_sampler': rng} for rng in model_rngs]
+            results = parallel_progress(fcn, model_rng_kwargs, show_progress=show_progress, num_workers=num_workers)
         elif show_progress:
             # Print the progress bar if the logging level was INFO or lower
             # This means that the user can still set the logging level higher e.g. WARNING to suppress output from Atomica in general
diff --git a/atomica/utils.py b/atomica/utils.py
index 023ffb2e..1fc84410 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -904,7 +904,10 @@ def callback(result, idx):
             pool.apply_async(fcn, callback=partial(callback, idx=i))
     else:
         for i, x in enumerate(inputs):
-            pool.apply_async(fcn, args=(x,), callback=partial(callback, idx=i))
+            if isinstance(x, dict):
+                pool.apply_async(fcn, args=(), kwds=x, callback=partial(callback, idx=i)) #if the list of inputs is given as a dictionary then pass as kwargs (kwds) instead
+            else:
+                pool.apply_async(fcn, args=(x,), callback=partial(callback, idx=i))
 
     pool.close()
     pool.join()

From 930c10e65674d7fe258953997a98e889beb87378 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 05:03:09 +1100
Subject: [PATCH 10/33] Optional acceptance criteria for sampled_sims - to
 ensure that sampled runs are reasonably matched historically

---
 atomica/project.py | 31 ++++++++++++++++++++++++++-----
 1 file changed, 26 insertions(+), 5 deletions(-)

diff --git a/atomica/project.py b/atomica/project.py
index f783a8d6..d1fbcdfc 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -518,7 +518,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
 
         return result
 
-    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, rand_seed: int = None, parallel=False, max_attempts=None, num_workers=None) -> list:
+    def run_sampled_sims(self, parset, progset=None, progset_instructions=None, result_names=None, n_samples: int = 1, rand_seed: int = None, parallel=False, max_attempts=None, num_workers=None, acceptance_criteria=[]) -> list:
         """
         Run sampled simulations
 
@@ -543,6 +543,10 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         :param max_attempts: Number of retry attempts for bad initializations
         :param num_workers: If ``parallel`` is True, this determines the number of parallel workers to use (default is usually number of CPUs)
         :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+        :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters in given years
+            tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't': 2018}]
+
+        
         :return: A list of Results that can be passed to `Ensemble.update()`. If multiple instructions are provided, the return value of this
                  function will be a list of lists, where the inner list iterates over different instructions for the same parset/progset samples.
                  It is expected in that case that the Ensemble's mapping function would take in a list of results
@@ -575,7 +579,7 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         model_rngs = [np.random.default_rng(seed = seed) for seed in seed_samples] #generate a RNG for each model
 
         if parallel:
-            fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts)
+            fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts, acceptance_criteria=acceptance_criteria)
             #as multiprocessing does not handle partial functions as compiled functions, need to send the rngs as kwargs in a dictionary, not as args to the partial function
             model_rng_kwargs = [{'rng_sampler': rng} for rng in model_rngs]
             results = parallel_progress(fcn, model_rng_kwargs, show_progress=show_progress, num_workers=num_workers)
@@ -584,9 +588,9 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
             # This means that the user can still set the logging level higher e.g. WARNING to suppress output from Atomica in general
             # (including any progress bars)
             with Quiet():
-                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng) for rng in tqdm.tqdm(model_rngs)]
+                results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng, acceptance_criteria=acceptance_criteria) for rng in tqdm.tqdm(model_rngs)]
         else:
-            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng) for rng in model_rngs]
+            results = [_run_sampled_sim(self, parset, progset, progset_instructions, result_names, max_attempts=max_attempts, rng_sampler=rng, acceptance_criteria=acceptance_criteria) for rng in model_rngs]
 
         return results
 
@@ -727,7 +731,7 @@ def __setstate__(self, d):
         P = migrate(self)
         self.__dict__ = P.__dict__
 
-def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler =  None):
+def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler =  None, acceptance_criteria = []):
     """
     Internal function to run simulation with sampling
 
@@ -749,6 +753,8 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     :param result_names: A list of result names (strings)
     :param max_attempts: Maximum number of sampling attempts before raising an error
     :param rng_sampler: Optional random number generator that may have been seeded to generate consistent results
+    :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters over a given period of time
+        tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't_range': (2018, 2018.99)}]
     :return: A list of results that either contains 1 result, or the same number of results as instructions
 
     """
@@ -770,6 +776,21 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
                 results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_sampler) for x, y in zip(progset_instructions, result_names)]
             else:
                 results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_sampler) for y in result_names]
+                
+            for criteria in acceptance_criteria:
+                par = criteria['parameter']
+                pop = criteria['population']
+                val = criteria['value']
+                t_range = criteria['t_range']
+                for res in results:
+                    res_par = res.get_variable(par, pop)[0]
+                    t_inds = np.where(np.logical_and(res_par.t>=t_range[0], res_par.t<=t_range[1]))
+                    res_t_val = np.mean(res_par.vals[t_inds]) #annualized so always average
+                    if res_t_val <= val[0] or res_t_val >= val[1]:
+                        print (f'Rejecting run as sampled sim value {res_t_val} outside of acceptable bound {val} for parameter {par}, population {pop} at time {t_range}')
+                                  
+                        raise BadInitialization(f'Rejecting run as sampled sim value {res_t_val} outside of acceptable bound {val} for parameter {par}, population {pop} at time {t_range}')
+                
             return results
         except BadInitialization:
             attempts += 1

From 0b4e49394fae8a650aae216018b8e48ed861e8ac Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 05:03:40 +1100
Subject: [PATCH 11/33] Update project.py

---
 atomica/project.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/atomica/project.py b/atomica/project.py
index d1fbcdfc..9cd5d568 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -787,8 +787,6 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
                     t_inds = np.where(np.logical_and(res_par.t>=t_range[0], res_par.t<=t_range[1]))
                     res_t_val = np.mean(res_par.vals[t_inds]) #annualized so always average
                     if res_t_val <= val[0] or res_t_val >= val[1]:
-                        print (f'Rejecting run as sampled sim value {res_t_val} outside of acceptable bound {val} for parameter {par}, population {pop} at time {t_range}')
-                                  
                         raise BadInitialization(f'Rejecting run as sampled sim value {res_t_val} outside of acceptable bound {val} for parameter {par}, population {pop} at time {t_range}')
                 
             return results

From 0e3bb4590831aaf0055a091c578d4c0e25cb249d Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 8 Nov 2021 23:31:08 +1100
Subject: [PATCH 12/33] Acceptance criteria evaluated during runtime or
 postprocessing as appropriate

Random seeds for runs ensures consistency (some needless rerunning of seeded runs that will fail between scenarios though)
---
 atomica/model.py   | 71 +++++++++++++++++++++++++++++++++++++++++++---
 atomica/project.py | 33 +++++++++------------
 2 files changed, 81 insertions(+), 23 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index 81aeae9d..3ded6a95 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -1212,6 +1212,7 @@ def set_dynamic(self, progset=None) -> None:
             # Pop aggregations and derivatives are handled in `update_pars()` so we set the dynamic flag
             # because the values are modified during integration
             self._is_dynamic = True
+            
 
         if self.deps:  # If there are no dependencies, then we know that this is precompute-only
             for deps in self.deps.values():  # deps is {'dep_name':[dep_objects]}
@@ -2013,7 +2014,7 @@ def report_characteristic(charac, n_indent=0):
 class Model:
     """ A class to wrap up multiple populations within model and handle cross-population transitions. """
 
-    def __init__(self, settings, framework, parset, progset=None, program_instructions=None, rng_sampler=None):
+    def __init__(self, settings, framework, parset, progset=None, program_instructions=None, rng_sampler=None, acceptance_criteria=[]):
 
         # Note that if a progset is provided and program instructions are not, then programs will not be
         # turned on. However, the progset is still available so that the coverage can still be probed
@@ -2035,6 +2036,7 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.dt = settings.sim_dt  #: Simulation time step
         self.stochastic = settings.stochastic #: Run with discrete numbers of people per compartment instead of fractions.
         
+        self.acceptance_criteria = sc.dcp(acceptance_criteria) #: If the model fails to adhere to these criteria, raise a BadInitialization exception
 
         if rng_sampler is None:
             rng_sampler = np.random.default_rng()
@@ -2258,7 +2260,7 @@ def build(self, parset):
                 # Make the parameter dynamic, because it needs to be computed during integration
                 for par in pars:
                     par.set_dynamic()
-
+                    
                 # Also make the variable being aggregated dynamic
                 for var in self._vars_by_pop[pars[0].pop_aggregation[1]]:
                     var.set_dynamic(progset=self.progset)
@@ -2267,6 +2269,8 @@ def build(self, parset):
                 if len(pars[0].pop_aggregation) > 3:
                     for var in self._vars_by_pop[pars[0].pop_aggregation[3]]:
                         var.set_dynamic(progset=self.progset)
+                        
+                
 
         # Set execution order - needs to be done _after_ pop aggregations have been flagged as dynamic
         self._set_exec_order()
@@ -2309,6 +2313,23 @@ def build(self, parset):
             for obj in pop.comps + pop.characs + pop.links:
                 obj.preallocate(self.t, self.dt)
             pop.initialize_compartments(parset, self.framework, self.t[0])
+            
+        #More finally, set up acceptance criteria for those which should be checked at runtime (variables that have ._is_dynamic = True), and those which can only be checked after conclusion
+        if self.acceptance_criteria:
+            for criteria in self.acceptance_criteria:
+                pars = self._vars_by_pop[criteria['parameter']]
+                for par in pars:
+                    if par.pop.name == criteria['population']:
+                        if par._is_dynamic:
+                            criteria['evaluate_at'] = self.t[np.where(self.t > criteria['t_range'][1])][0]
+                        else:
+                            criteria['evaluate_at'] = np.inf #can only evaluate postcompute
+                            
+            #sort by evaluation time - now we only have to look at the first element every timestep
+            self.acceptance_criteria = sorted(self.acceptance_criteria, key=lambda item: item['evaluate_at'])
+            self.acceptance_index = 0 #which index to evaluate
+            
+
 
     def _set_exec_order(self) -> None:
         """
@@ -2433,6 +2454,8 @@ def process(self) -> None:
             self.update_comps()
             self.update_pars()
             self.update_links()
+            if self.acceptance_criteria:
+                self.update_acceptance()
 
         # Update postcompute parameters - note that it needs to be done in execution order
         for par_name in self._exec_order["all_pars"]:
@@ -2440,6 +2463,9 @@ def process(self) -> None:
                 if par.fcn_str and not (par._is_dynamic or par._precompute):
                     par.update()
                     par.constrain()
+                    
+        if self.acceptance_criteria: #call update_acceptance once more for any variables that could only be assessed postcompute
+                self.update_acceptance(evaluate_all = True)
 
         # Clear characteristic internal storage and switch to dynamic computation to save space
         for pop in self.pops:
@@ -2447,6 +2473,7 @@ def process(self) -> None:
                 charac._vals = None
 
         self._program_cache = None  # Drop the program cache afterwards to save space
+        
 
     def update_links(self) -> None:
         """
@@ -2555,6 +2582,41 @@ def update_comps(self) -> None:
         for pop in self.pops:
             for comp in pop.comps:
                 comp.update(ti)
+                
+    def update_acceptance(self, evaluate_all = False) -> None:
+        """
+        Update any acceptance criteria for the model run, raise BadInitialization error if the model run is failing
+        
+        :param evaluate_all: Final evaluation of all acceptance criteria.
+        
+        """
+        
+        time = self.t[self._t_index] #actual time
+
+        #We know that these are sorted by evaluation time so we only need to look at the first index
+        while (self.acceptance_index < len(self.acceptance_criteria) ) and (evaluate_all or time >= self.acceptance_criteria[self.acceptance_index]['evaluate_at']):
+            #time ranges for assessment are in order, so only need to check the first one
+            assessable = self.acceptance_criteria[self.acceptance_index]
+            #we have passed the time window, assess if the criteria was met - either return a fail condition or remove this acceptance criteria
+            a_par = assessable['parameter']
+            a_pop = assessable['population']
+            a_times = assessable['t_range']
+            accept_vals = assessable['value']
+            a_t_inds = np.where(np.logical_and(self.t>=a_times[0], self.t<=a_times[1])) #self.t is the full tvec for the model
+            
+            pars = self._vars_by_pop[a_par]
+            for par in pars: #TODO should be more efficient way to directly access the par values for the pop rather than looping
+                if par.pop.name == a_pop: 
+                    a_t_val = np.mean(par[a_t_inds]) #mean across the time range
+                    if a_t_val <= accept_vals[0] or a_t_val >= accept_vals[1]:
+                        # print(f'{time} Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+                        raise BadInitialization(f'Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+                    # else:
+                        # print(f'{time} ACCEPTING run as sampled sim value {a_t_val} inside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+            assessable['assessed_value'] = a_t_val
+            self.acceptance_index += 1
+                  
+        
 
     def flush_junctions(self) -> None:
         """
@@ -2689,7 +2751,7 @@ def update_pars(self) -> None:
                     par.constrain(ti)
 
 
-def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler = None):
+def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler = None, acceptance_criteria = []):
     """
     Build and process model
 
@@ -2714,10 +2776,11 @@ def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = N
     :param program_instructions: Optional :class:`ProgramInstructions` instance. If ``progset`` is specified, then instructions must be provided
     :param name: Optionally specify the name to assign to the output result
     :param rng_sampler: Optionally specify a random number generator that may have been seeded to generate consistent results
+    
     :return: A :class:`Result` object containing the processed model
 
     """
 
-    m = Model(settings, framework, parset, progset, program_instructions, rng_sampler=rng_sampler)
+    m = Model(settings, framework, parset, progset, program_instructions, rng_sampler=rng_sampler, acceptance_criteria=acceptance_criteria)
     m.process()
     return Result(model=m, parset=parset, name=name)
diff --git a/atomica/project.py b/atomica/project.py
index 9cd5d568..2de48dc0 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -475,7 +475,7 @@ def update_settings(self, sim_start=None, sim_end=None, sim_dt=None):
         """ Modify the project settings, e.g. the simulation time vector. """
         self.settings.update_time_vector(start=sim_start, end=sim_end, dt=sim_dt)
 
-    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng = None):
+    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng = None, acceptance_criteria=[]):
         """
         Run a single simulation
 
@@ -489,6 +489,9 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
         :param store_results: If True, then the result will automatically be stored in ``self.results``
         :param result_name: Optionally assign a specific name to the result (otherwise, a unique default name will automatically be selected)
         :param rng: Optionally a random number generator that may have been seeded to generate consistent results, or a random_seed used to generate a Generator
+        :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters over a given period of time
+        tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't_range': (2018, 2018.99)}]
+        
         :return: A :class:`Result` instance
 
         """
@@ -511,7 +514,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
                 k += 1
 
         tm = sc.tic()
-        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name, rng_sampler=rng_sampler)
+        result = run_model(settings=self.settings, framework=self.framework, parset=parset, progset=progset, program_instructions=progset_instructions, name=result_name, rng_sampler=rng_sampler, acceptance_criteria=acceptance_criteria)
         logger.info('Elapsed time for running "%s": %ss', self.name, sc.sigfig(sc.toc(tm, output=True), 3))
         if store_results:
             self.results.append(result)
@@ -766,29 +769,21 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
     
     if rng_sampler is None:
         rng_sampler = np.random.default_rng()
+    
+    #Set up separate seeds for each attempt and seed a rng for each one (in case scenarios diverge in later timepoints, it's important to START from the same seeds)
+    rand_seeds = rng_sampler.integers(1e15, size=max_attempts)
 
     attempts = 0
     while attempts < max_attempts:
         try:
-            sampled_parset = parset.sample(rng_sampler=rng_sampler)
+            rng_attempt = np.random.default_rng(seed=rand_seeds[attempts])
+            sampled_parset = parset.sample(rng_sampler=rng_attempt)
             if progset:
-                sampled_progset = progset.sample(rng_sampler=rng_sampler)
-                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_sampler) for x, y in zip(progset_instructions, result_names)]
+                sampled_progset = progset.sample(rng_sampler=rng_attempt)
+                results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for x, y in zip(progset_instructions, result_names)]
             else:
-                results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_sampler) for y in result_names]
-                
-            for criteria in acceptance_criteria:
-                par = criteria['parameter']
-                pop = criteria['population']
-                val = criteria['value']
-                t_range = criteria['t_range']
-                for res in results:
-                    res_par = res.get_variable(par, pop)[0]
-                    t_inds = np.where(np.logical_and(res_par.t>=t_range[0], res_par.t<=t_range[1]))
-                    res_t_val = np.mean(res_par.vals[t_inds]) #annualized so always average
-                    if res_t_val <= val[0] or res_t_val >= val[1]:
-                        raise BadInitialization(f'Rejecting run as sampled sim value {res_t_val} outside of acceptable bound {val} for parameter {par}, population {pop} at time {t_range}')
-                
+                results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for y in result_names]
+            
             return results
         except BadInitialization:
             attempts += 1

From 211a7fd1e50ea2848fd6241a5dbca9674f24e2a5 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Wed, 1 Mar 2023 16:21:39 +1100
Subject: [PATCH 13/33] Fix ResidualJunctionCompartment not applying
 stochasticity correctly with a ">" residual link.  - `outflow` wasn't updated
 with the stochastic outflow_fractions, and  - The probabilities need to sum
 to 1 in a multinomial so we give the residual probability to the residual
 link  - Also checks that net_inflow is an integer and throws an error if it
 is not

---
 atomica/model.py | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index 458af3a1..8ccf60b6 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -608,26 +608,37 @@ def balance(self, ti: int) -> None:
                 net_inflow += link.vals[ti]  # If not part of a duration group, get scalar flow from Link.vals
 
         outflow_fractions = np.zeros(len(self.outlinks))
+        residual_link_index = None
         for i, link in enumerate(self.outlinks):
             if link.parameter is not None:
                 outflow_fractions[i] = link.parameter.vals[ti]
             else:
+                residual_link_index = i
                 outflow_fractions[i] = 0
 
         total_outflow = outflow_fractions.sum()
         if total_outflow > 1:
             outflow_fractions /= total_outflow
 
+        # assign outflow stochastically to compartments based on the probabilities
+        if self.stochastic and net_inflow > 0.:
+            if np.abs(np.round(net_inflow)-net_inflow).sum() > 1e-3:
+                raise Exception(f'Net inflow should be an integer for a stochastic run, was: {net_inflow}')
+
+            # Give residual probability to the residual link if there is one, otherwise normalize probabilities
+            if residual_link_index is not None:
+                outflow_fractions[residual_link_index] = max(0, 1 - outflow_fractions.sum())
+            else: # Probabilities have to sum to one as we have to flush the junction, and multinomials probabilities should sum to 1
+                outflow_fractions = outflow_fractions / outflow_fractions.sum()
+
+            outflow_fractions = self.rng_sampler.multinomial(np.round(net_inflow).astype(int), outflow_fractions, size=1)[0]/(net_inflow)
+
         # Calculate outflows
         outflow = net_inflow.reshape(-1, 1) * outflow_fractions.reshape(1, -1)
 
-        # assign outflow stochastically to compartments based on the probabilities
-        if self.stochastic and net_inflow > 0.:
-            outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0]/(net_inflow * sum(outflow_fractions))
-            
         # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
         for frac, link in zip(outflow_fractions, self.outlinks):
-            if link.parameter is None and total_outflow < 1:
+            if link.parameter is None and total_outflow < 1 and not self.stochastic:  # If we are stochastic then we have already taken the residual link into account
                 flow = net_inflow - np.sum(outflow, axis=1)  # Sum after multiplying by outflow fractions to reduce numerical precision errors and enforce conserved quantities more accurately
             else:
                 flow = net_inflow * frac

From 98b7274a7f6ba7a9578051f31389f78e18cd5d04 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Fri, 3 Mar 2023 21:49:49 +1100
Subject: [PATCH 14/33] Whitespace

---
 atomica/model.py      | 2 +-
 atomica/parameters.py | 2 +-
 atomica/programs.py   | 2 +-
 atomica/project.py    | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index 8ccf60b6..fbf220d0 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -228,7 +228,7 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """A class to wrap up data for one compartment within a cascade network."""
 
-    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
+    def __init__(self, pop, name, stochastic: bool = False, rng_sampler = None):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
diff --git a/atomica/parameters.py b/atomica/parameters.py
index dc8cf465..b71da570 100644
--- a/atomica/parameters.py
+++ b/atomica/parameters.py
@@ -90,7 +90,7 @@ def interpolate(self, tvec, pop_name: str) -> np.array:
 
         return self.ts[pop_name].interpolate(tvec, method=self._interpolation_method)
 
-    def sample(self, constant: bool, rng_sampler = None) -> None:
+    def sample(self, constant: bool, rng_sampler=None) -> None:
         """
         Perturb parameter based on uncertainties
 
diff --git a/atomica/programs.py b/atomica/programs.py
index e72804c5..d52c6957 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1198,7 +1198,7 @@ def is_one_off(self) -> bool:
         """
         return "/year" not in self.unit_cost.units
 
-    def sample(self, constant: bool, rng_sampler = None) -> None:
+    def sample(self, constant: bool, rng_sampler=None) -> None:
         """
         Perturb program values based on uncertainties
 
diff --git a/atomica/project.py b/atomica/project.py
index 05e90590..4208003e 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -67,7 +67,7 @@ def sim_end(self):
     @property
     def sim_dt(self):
         return self._sim_dt
-    
+
     @property
     def stochastic(self):
         return self._stochastic

From a54571589bc8c916f85af4d74156339d88e14eab Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Wed, 1 Nov 2023 16:52:25 +1100
Subject: [PATCH 15/33] Add rng_sampler=None to Programset.sample()

---
 atomica/programs.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/atomica/programs.py b/atomica/programs.py
index d52c6957..9cdbc447 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1118,7 +1118,7 @@ def get_outcomes(self, prop_coverage: dict) -> dict:
 
         return {(covout.par, covout.pop): covout.get_outcome(prop_coverage) for covout in self.covouts.values()}
 
-    def sample(self, constant: bool = True):
+    def sample(self, constant: bool = True, rng_sampler = None):
         """
         Perturb programs based on uncertainties
 
@@ -1135,9 +1135,9 @@ def sample(self, constant: bool = True):
 
         new = sc.dcp(self)
         for prog in new.programs.values():
-            prog.sample(constant)
+            prog.sample(constant, rng_sampler=rng_sampler)
         for covout in new.covouts.values():
-            covout.sample()
+            covout.sample(rng_sampler=rng_sampler)
         return new
 
 

From 72f231dd69c33e203a5b644c189424ca299b679f Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 3 Nov 2023 08:42:18 +1100
Subject: [PATCH 16/33] Keep "Output" column in worksheet and don't group by
 output, just print the whole dataframe as it is

---
 atomica/results.py | 27 ++++++++++-----------------
 1 file changed, 10 insertions(+), 17 deletions(-)

diff --git a/atomica/results.py b/atomica/results.py
index 3de0251f..af19b977 100644
--- a/atomica/results.py
+++ b/atomica/results.py
@@ -837,30 +837,23 @@ def _write_df(writer, formats, sheet_name, df, level_ordering):
     for level in level_ordering:
         order[level] = list(dict.fromkeys(df.index.get_level_values(level)))
 
-    required_width = [0] * (len(level_ordering) - 1)
-
-    row = 0
+    df = df.reorder_levels(level_ordering)
+    for i in range(len(level_ordering)):
+        df = df.reindex(order[level_ordering[i]], level=i)
 
     worksheet = writer.book.add_worksheet(sheet_name)
     writer.sheets[sheet_name] = worksheet  # Need to add it to the ExcelWriter for it to behave properly
 
-    for title, table in df.groupby(level=level_ordering[0], sort=False):
-        worksheet.write_string(row, 0, title, formats["center_bold"])
-        row += 1
+    df.index = df.index.set_names([level_substitutions[x] if x in level_substitutions else x.title() for x in df.index.names])
 
-        table.reset_index(level=level_ordering[0], drop=True, inplace=True)  # Drop the title column
-        table = table.reorder_levels(level_ordering[1:])
-        for i in range(1, len(level_ordering)):
-            table = table.reindex(order[level_ordering[i]], level=i - 1)
-        table.index = table.index.set_names([level_substitutions[x] if x in level_substitutions else x.title() for x in table.index.names])
-        table.to_excel(writer, sheet_name, startcol=0, startrow=row, merge_cells=False)
-        row += table.shape[0] + 2
+    row = 0
+    df.to_excel(writer, sheet_name, startcol=0, startrow=row, merge_cells=False)
 
-        required_width[0] = max(required_width[0], len(title))
-        for i in range(0, len(required_width)):
-            required_width[i] = max(required_width[i], max(table.index.get_level_values(i).str.len()))
+    required_width = [0] * len(level_ordering)
+    for i in range(len(required_width)):
+        required_width[i] = max(len(val) for val in order[level_ordering[i]])
 
-    for i in range(0, len(required_width)):
+    for i in range(len(required_width)):
         if required_width[i] > 0:
             worksheet.set_column(i, i, required_width[i] * 1.1 + 1)
 

From 76a1dba3caaa1dbc531e96ec6d26d240d3ab5ff4 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 3 Nov 2023 09:26:01 +1100
Subject: [PATCH 17/33] Also make sure the required_width of a column fits the
 column name

---
 atomica/results.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/atomica/results.py b/atomica/results.py
index af19b977..1d784075 100644
--- a/atomica/results.py
+++ b/atomica/results.py
@@ -849,7 +849,7 @@ def _write_df(writer, formats, sheet_name, df, level_ordering):
     row = 0
     df.to_excel(writer, sheet_name, startcol=0, startrow=row, merge_cells=False)
 
-    required_width = [0] * len(level_ordering)
+    required_width = [len(name) for name in df.index.names]
     for i in range(len(required_width)):
         required_width[i] = max(len(val) for val in order[level_ordering[i]])
 

From 1ee101c4b89c7b9ba633a2977d22fa7ff9e2fd1a Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 6 Nov 2023 11:12:11 +1100
Subject: [PATCH 18/33] Add "constant" to reserved_keywords

---
 atomica/system.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/atomica/system.py b/atomica/system.py
index c72ab22a..46587e6a 100644
--- a/atomica/system.py
+++ b/atomica/system.py
@@ -85,7 +85,7 @@ class FrameworkSettings:
     DEFAULT_SYMBOL_INAPPLICABLE = "N.A."
     DEFAULT_POP_TYPE = "default"
 
-    RESERVED_KEYWORDS = ["t", "flow", "all", "dt", "total"]  # A code_name in the framework cannot be equal to one of these values
+    RESERVED_KEYWORDS = ["t", "flow", "all", "dt", "total", "constant"]  # A code_name in the framework cannot be equal to one of these values
     RESERVED_KEYWORDS += supported_functions.keys()
 
     RESERVED_SYMBOLS = set(":,;/+-*'\" @")  # A code_name in the framework (for characs, comps, pars) cannot contain any of these characters

From 70f34cd28f3b486a78f357db3dd3e2d6996b2d82 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 24 Nov 2023 15:12:05 +1100
Subject: [PATCH 19/33] Use ax = fig.gca() instead of ax = plt.gca() so that
 you can have multiple figures plotting in parallel

---
 atomica/results.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/atomica/results.py b/atomica/results.py
index 1d784075..a1d2d717 100644
--- a/atomica/results.py
+++ b/atomica/results.py
@@ -1221,7 +1221,7 @@ def plot_series(self, fig=None, style="quartile", results=None, outputs=None, po
 
         if fig is None:
             fig = plt.figure()
-        ax = plt.gca()
+        ax = fig.gca()
 
         series_lookup = self._get_series()
 
@@ -1256,7 +1256,7 @@ def plot_series(self, fig=None, style="quartile", results=None, outputs=None, po
 
             proposed_label = "%s (%s)" % (output, these_series[0].unit_string)
             if ax.yaxis.get_label().get_text():
-                assert proposed_label == ax.yaxis.get_label().get_text(), "The outputs being superimposed have different units"
+                assert proposed_label == ax.yaxis.get_label().get_text(), f"The outputs being superimposed have different units: {proposed_label} vs {ax.yaxis.get_label().get_text()} {self} {results} {outputs}"
             else:
                 ax.set_ylabel(proposed_label)
 

From 02c4cb99a16e21f12f31f587bef60b70430607c9 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 24 Nov 2023 15:16:41 +1100
Subject: [PATCH 20/33] Improve error message when different y-labels

---
 atomica/results.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/atomica/results.py b/atomica/results.py
index a1d2d717..0194d309 100644
--- a/atomica/results.py
+++ b/atomica/results.py
@@ -1256,7 +1256,7 @@ def plot_series(self, fig=None, style="quartile", results=None, outputs=None, po
 
             proposed_label = "%s (%s)" % (output, these_series[0].unit_string)
             if ax.yaxis.get_label().get_text():
-                assert proposed_label == ax.yaxis.get_label().get_text(), f"The outputs being superimposed have different units: {proposed_label} vs {ax.yaxis.get_label().get_text()} {self} {results} {outputs}"
+                assert proposed_label == ax.yaxis.get_label().get_text(), f"The outputs being superimposed have different units: {proposed_label} vs {ax.yaxis.get_label().get_text()}"
             else:
                 ax.set_ylabel(proposed_label)
 

From 269f71244bba8bd324e6f0e7d847a2fc5669a9fe Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Thu, 11 Apr 2024 13:32:17 +1000
Subject: [PATCH 21/33] Add test_tox_markovchain.py

---
 tests/test_tox_markovchain.py | 175 ++++++++++++++++++++++++++++++++++
 1 file changed, 175 insertions(+)
 create mode 100644 tests/test_tox_markovchain.py

diff --git a/tests/test_tox_markovchain.py b/tests/test_tox_markovchain.py
new file mode 100644
index 00000000..656e3ebe
--- /dev/null
+++ b/tests/test_tox_markovchain.py
@@ -0,0 +1,175 @@
+import atomica as at
+import numpy as np
+import os
+import matplotlib.pyplot as plt
+import sciris as sc
+import pytest
+
+
+
+dirname = at.parent_dir()
+# np.seterr(all='raise')
+
+
+def test_functional_all_framework_databook():
+    framework_files = sc.getfilelist(folder=dirname, pattern='*framework*.xlsx')
+    databook_files = sc.getfilelist(folder=dirname, pattern='*databook*.xlsx')
+
+    for framework_file in framework_files:
+        for databook_file in databook_files:
+
+            try:
+                P = at.Project(framework=framework_file, databook=databook_file, do_run=False)
+                P.run_sim(result_name='Original')
+            except:
+                # print(f'Failed when running normally so skipping')
+                continue  # Running without MC didn't work so skip
+
+            print(f'\nTesting framework, databook which runs normally: {framework_file}, {databook_file}')
+
+            try:
+                P.settings.stochastic = True
+                P.run_sim(result_name='DTMC')
+            except Exception as e:
+                raise Exception(f'Ran normally but failed to run with P.settings.stochastic = True. Framework, databook = {framework_file}, {databook_file}') \
+                    from e
+
+            print('Success', framework_file, databook_file, '\n')
+
+    print('\nPASS!')
+
+
+def test_compare_deterministic_markovchain():
+
+    seed = 0
+    N = 1000
+
+    y_cutoff = 10  # (deterministic) Compartment values below this will not be considered for the comparisons
+    index_cutoff = 1500
+    y_diff_cutoff = 0
+
+    ratio_cutoff = 0.05
+
+
+
+    def filter_points(t, det, mc, mc_mean):
+        # t_fil, det_fil, mc_fil = sc.dcp(t), sc.dcp(det), sc.dcp(mc)
+
+        indices = np.ones(t.shape)
+        indices = np.logical_and(indices, det >= y_cutoff)
+        indices = np.logical_and(indices, np.arange(len(indices)) <= index_cutoff)
+        indices = np.logical_and(indices, np.abs(det - mc_mean) > y_diff_cutoff)
+
+        return t[indices], det[indices], mc[:, indices], mc_mean[indices]
+
+    def get_array_outputs(results, output, pop):
+        extracted = [None] * len(results)
+        for i, res in enumerate(results):
+            extracted[i] = res.model.get_pop(pop).get_variable(output)[0].vals
+
+        array = np.stack(extracted, axis=0)
+        t = results[0].model.get_pop(pop).get_variable(output)[0].t
+        return t, array
+
+    def get_average_output(results, output, pop):
+        t, array = get_array_outputs(results, output, pop)
+        average = np.mean( array , axis=0)
+
+        return t, average, array
+
+    def do_hypothesis_test(det_array, mc_results):
+        import scipy
+        res = scipy.stats.ttest_1samp(mc_results, det_array, keepdims=True)
+        return res
+
+    def run_simple_percent_test(P, baseline_results, mc_results):
+        baseline_results = sc.promotetolist(baseline_results)
+        mc_results = sc.promotetolist(mc_results)
+
+        compartment_names = [name for name in baseline_results[0].model.framework.comps.index.values]
+
+        max_max_diff = 0
+
+        for compartment_name in compartment_names:
+            for pop in baseline_results[0].pop_names:
+                try: baseline_results[0].model.get_pop(pop).get_variable(compartment_name)[0].vals
+                except: continue  # This population doesn't have this compartment
+
+                t, mc_mean, mc_array = get_average_output(mc_results, compartment_name, pop)
+                _, det_mean, _ = get_average_output(baseline_results, compartment_name, pop)
+
+                t, det_mean, mc_array, mc_mean = filter_points(t, det_mean, mc_array, mc_mean)
+
+                if len(t) == 0:  # No filtered outputs so skip this output, pop combo
+                    continue
+
+                indices = det_mean != 0
+                ratio = np.ones(mc_mean.shape)
+                ratio[indices] = mc_mean[indices] / det_mean[indices]
+                if min(ratio) < 1-ratio_cutoff or max(ratio) > 1+ratio_cutoff:
+                    print('RATIO', compartment_name, pop, min(ratio), max(ratio), list(zip(ratio, det_mean, mc_mean)))
+                    indices = np.logical_or(ratio < 1-ratio_cutoff, ration > 1+ratio_cutoff)
+
+                    raise Exception(f'The Markov Chain runs (N={len(mc_results)} differed too much from the baseline runs (N={len(baseline_results)} at starting at timestep={indices[0]}, t={t[indices[0]]}, with compartment={compartment_name}, population={pop}. Baseline={det_mean[indices[0]]}, MC={mc_mean[indices[0]]}, difference={100*(mc_mean[indices[0]]/det_mean[indices[0]] - 1):.2f}')
+
+                max_diff = np.max(np.abs(ratio - 1))
+                max_max_diff = max(max_diff, max_max_diff)
+
+        print(f'Success: max difference {max_max_diff*100:.2f}%')
+
+
+
+    seed = np.random.default_rng(seed).integers(2**32-1)
+
+    framework_files = sc.getfilelist(folder=dirname, pattern='*framework*.xlsx')
+    databook_files  = sc.getfilelist(folder=dirname, pattern='*databook*.xlsx')
+
+    projects = []
+
+    for framework_file in framework_files:
+        for databook_file in databook_files:
+
+            try:
+                P = at.Project(framework=framework_file, databook=databook_file, do_run=False)
+                res = P.run_sim(result_name='Original')
+
+                baseline_results = [res]
+
+                if 'stochastic' in framework_file or 'stochastic' in databook_file:
+                    for i in range(N):
+                        baseline_results.append(P.run_sim(result_name=f'Original {i}'))
+
+            except:
+                # print(f'Failed when running normally so skipping')
+                continue  # Running without MC didn't work so skip
+
+
+            print(f'\nTrying framework, databook: {framework_file}, {databook_file}')
+
+            try:
+                at.logger.setLevel("WARN")
+                P.settings.stochastic = True
+                mc_results = [None] * N
+                for i in range(N):
+                    this_seed = ( i + seed + hash(framework_file) + hash(databook_file) )% (2**32)
+                    mc_results[i] = P.run_sim(result_name=f'DTMC {i}', rng=this_seed )
+
+                at.logger.setLevel("INFO")
+            except:
+                print('Ran normally but failed to run with P.settings.stochastic = True')
+                continue # No success so continue without printing success
+
+            # plot_comparison(P, baseline_results=baseline_results, mc_results=mc_results)
+            run_simple_percent_test(P, baseline_results=baseline_results, mc_results=mc_results)
+
+
+            print('Success', framework_file, databook_file)
+            projects.append(P)
+
+    print('\nPASS!')
+
+
+
+if __name__ == '__main__':
+    test_functional_all_framework_databook()
+    test_compare_deterministic_markovchain()

From 9dcb653df16727ecd60c8a29261c02bb9ea9f552 Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Fri, 7 Jun 2024 15:25:24 +1000
Subject: [PATCH 22/33] Fix tests

---
 atomica/demos.py              | 54 ++++++++++----------
 tests/test_tox_library.py     |  2 +-
 tests/test_tox_markovchain.py | 95 ++++++++---------------------------
 3 files changed, 49 insertions(+), 102 deletions(-)

diff --git a/atomica/demos.py b/atomica/demos.py
index b56ed912..baee7f61 100644
--- a/atomica/demos.py
+++ b/atomica/demos.py
@@ -11,10 +11,29 @@
 from .scenarios import BudgetScenario
 from .utils import TimeSeries
 
-__all__ = ["demo", "make_demo_scenarios"]
-
-
-def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Project:
+__all__ = ["demo", "demo_projects", "make_demo_scenarios"]
+
+
+demo_projects = [
+    "udt",
+    "udt_dyn",
+    "usdt",
+    "cervicalcancer",
+    "sir",
+    "diabetes",
+    "combined",
+    'service',
+    "hypertension",
+    "hypertension_dyn",
+    "hiv",
+    "hiv_dyn",
+    "tb_simple",
+    "tb_simple_dyn",
+    "tb",
+    "dt",
+]
+
+def demo(which: str = 'sir', do_run: bool = True, addprogs: bool = True) -> Project:
     """
     Return a demo project
 
@@ -25,30 +44,11 @@ def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Proje
 
     """
 
-    options = [
-        "udt",
-        "udt_dyn",
-        "usdt",
-        "cervicalcancer",
-        "sir",
-        "diabetes",
-        "combined",
-        # 'service',
-        "hypertension",
-        "hypertension_dyn",
-        "hiv",
-        "hiv_dyn",
-        "tb_simple",
-        "tb_simple_dyn",
-        "environment",
-        "tb",
-    ]
-
-    dtdict = sc.odict.fromkeys(options, 1.0)
+    dtdict = sc.odict.fromkeys(demo_projects, 1.0)
     dtdict["tb"] = 0.5
 
-    if which is None or which not in options:
-        raise Exception("Supported project types are:\n%s" % ("\n".join(options)))
+    if which not in demo_projects:
+        raise Exception("Supported demonstration projects are:\n%s" % ("\n".join(demos)))
 
     framework = LIBRARY_PATH / f"{which}_framework.xlsx"
     databook = LIBRARY_PATH / f"{which}_databook.xlsx"
@@ -61,7 +61,7 @@ def demo(which: str = None, do_run: bool = True, addprogs: bool = True) -> Proje
     if do_run:
         P.run_sim(store_results=True)
 
-    if addprogs:
+    if addprogs and progbook.exists():
         logger.debug("Loading progbook")
         P.load_progbook(progbook)
 
diff --git a/tests/test_tox_library.py b/tests/test_tox_library.py
index 704742dc..16e21d58 100644
--- a/tests/test_tox_library.py
+++ b/tests/test_tox_library.py
@@ -173,7 +173,7 @@ def run_regenerated_framework(proj):
 # each function call as a separate test, which can be run in parallel
 # simply by adding `pytest -n 4` for instance. Or via tox
 # `tox -- -o -n 4` (because the default config is to use n=2 for TravisCI)
-@pytest.mark.parametrize("model", models)
+@pytest.mark.parametrize("model", at.demo_projects)
 def test_model(model):
 
     framework_file = at.LIBRARY_PATH / f"{model}_framework.xlsx"
diff --git a/tests/test_tox_markovchain.py b/tests/test_tox_markovchain.py
index 656e3ebe..a06ae9a7 100644
--- a/tests/test_tox_markovchain.py
+++ b/tests/test_tox_markovchain.py
@@ -10,36 +10,15 @@
 dirname = at.parent_dir()
 # np.seterr(all='raise')
 
+@pytest.mark.parametrize("demo", at.demo_projects)
+def test_functional_all_framework_databook(demo):
+    # Test that all demo frameworks can be run with stochastic mode turned on
+    P = at.demo(demo, do_run=False)
+    P.settings.stochastic = True
+    P.run_sim(result_name='DTMC')
 
-def test_functional_all_framework_databook():
-    framework_files = sc.getfilelist(folder=dirname, pattern='*framework*.xlsx')
-    databook_files = sc.getfilelist(folder=dirname, pattern='*databook*.xlsx')
-
-    for framework_file in framework_files:
-        for databook_file in databook_files:
-
-            try:
-                P = at.Project(framework=framework_file, databook=databook_file, do_run=False)
-                P.run_sim(result_name='Original')
-            except:
-                # print(f'Failed when running normally so skipping')
-                continue  # Running without MC didn't work so skip
-
-            print(f'\nTesting framework, databook which runs normally: {framework_file}, {databook_file}')
-
-            try:
-                P.settings.stochastic = True
-                P.run_sim(result_name='DTMC')
-            except Exception as e:
-                raise Exception(f'Ran normally but failed to run with P.settings.stochastic = True. Framework, databook = {framework_file}, {databook_file}') \
-                    from e
-
-            print('Success', framework_file, databook_file, '\n')
-
-    print('\nPASS!')
-
-
-def test_compare_deterministic_markovchain():
+@pytest.mark.parametrize("demo", ['sir','udt','tb_simple']) # This test is expensive so only run it for a few lightweight models
+def test_compare_deterministic_markovchain(demo):
 
     seed = 0
     N = 1000
@@ -121,55 +100,23 @@ def run_simple_percent_test(P, baseline_results, mc_results):
 
     seed = np.random.default_rng(seed).integers(2**32-1)
 
-    framework_files = sc.getfilelist(folder=dirname, pattern='*framework*.xlsx')
-    databook_files  = sc.getfilelist(folder=dirname, pattern='*databook*.xlsx')
-
-    projects = []
-
-    for framework_file in framework_files:
-        for databook_file in databook_files:
-
-            try:
-                P = at.Project(framework=framework_file, databook=databook_file, do_run=False)
-                res = P.run_sim(result_name='Original')
-
-                baseline_results = [res]
-
-                if 'stochastic' in framework_file or 'stochastic' in databook_file:
-                    for i in range(N):
-                        baseline_results.append(P.run_sim(result_name=f'Original {i}'))
-
-            except:
-                # print(f'Failed when running normally so skipping')
-                continue  # Running without MC didn't work so skip
-
-
-            print(f'\nTrying framework, databook: {framework_file}, {databook_file}')
-
-            try:
-                at.logger.setLevel("WARN")
-                P.settings.stochastic = True
-                mc_results = [None] * N
-                for i in range(N):
-                    this_seed = ( i + seed + hash(framework_file) + hash(databook_file) )% (2**32)
-                    mc_results[i] = P.run_sim(result_name=f'DTMC {i}', rng=this_seed )
-
-                at.logger.setLevel("INFO")
-            except:
-                print('Ran normally but failed to run with P.settings.stochastic = True')
-                continue # No success so continue without printing success
-
-            # plot_comparison(P, baseline_results=baseline_results, mc_results=mc_results)
-            run_simple_percent_test(P, baseline_results=baseline_results, mc_results=mc_results)
 
+    P = at.demo(demo, do_run=False)
+    res = P.run_sim(result_name='Original')
 
-            print('Success', framework_file, databook_file)
-            projects.append(P)
+    baseline_results = [res]
 
-    print('\nPASS!')
+    P.settings.stochastic = True
+    mc_results = [None] * N
+    with at.Quiet():  # Using at.Quiet should automatically reset the logging level if an exception occurs
+        for i in range(N):
+            this_seed = ( i + seed + hash(demo) )% (2**32)
+            mc_results[i] = P.run_sim(result_name=f'DTMC {i}', rng=this_seed )
 
+    run_simple_percent_test(P, baseline_results=baseline_results, mc_results=mc_results)
 
 
 if __name__ == '__main__':
-    test_functional_all_framework_databook()
-    test_compare_deterministic_markovchain()
+    for demo in at.demo_projects:
+        test_functional_all_framework_databook()
+    test_compare_deterministic_markovchain('sir')

From e39b7806fbabcd299abaeee2474f5163edc4b51d Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Fri, 7 Jun 2024 15:41:34 +1000
Subject: [PATCH 23/33] Fix seed selection

---
 docs/examples/Uncertainty.ipynb | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/examples/Uncertainty.ipynb b/docs/examples/Uncertainty.ipynb
index daef6d9d..fa652784 100644
--- a/docs/examples/Uncertainty.ipynb
+++ b/docs/examples/Uncertainty.ipynb
@@ -95,7 +95,7 @@
     "import matplotlib.pyplot as plt\n",
     "from scipy import stats\n",
     "import sciris as sc\n",
-    "np.random.seed(3)\n",
+    "np.random.seed(5)\n",
     "testdir = '../../tests/'\n",
     "\n",
     "P = at.Project(framework=testdir + 'test_uncertainty_framework.xlsx', databook=testdir + 'test_uncertainty_databook.xlsx')\n",
@@ -146,7 +146,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If there is uncertainty, we can obtain a sampled epidemiological prediction by sampling from all of the quantities containing uncertinty. For example, replacing the value of `214` above with a sample from the distribution $\\ \\mathcal{N}(214,20)$ - along with all of the other quantities. The input to the model is effectively a parameter set derived from the original set of parameters, but with values replaced by sampled values where appropriate. You can obtain a sampled parameter set simply by calling the `ParameterSet.sample()` method:"
+    "If there is uncertainty, we can obtain a sampled epidemiological prediction by sampling from all of the quantities containing uncertainty. For example, replacing the value of `214` above with a sample from the distribution $\\ \\mathcal{N}(214,20)$ - along with all of the other quantities. The input to the model is effectively a parameter set derived from the original set of parameters, but with values replaced by sampled values where appropriate. You can obtain a sampled parameter set simply by calling the `ParameterSet.sample()` method:"
    ]
   },
   {
@@ -1035,9 +1035,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:atomica37]",
+   "display_name": "Python [conda env:atomica312]",
    "language": "python",
-   "name": "conda-env-atomica37-py"
+   "name": "conda-env-atomica312-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1049,7 +1049,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.12.2"
   },
   "toc": {
    "base_numbering": 1,

From 8abb36ac4f77d00a3cd2e56618b2e5ae7de94b9c Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Fri, 7 Jun 2024 15:42:13 +1000
Subject: [PATCH 24/33] Implement very simple example of test reordering for
 further use later

---
 tests/conftest.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/conftest.py b/tests/conftest.py
index 9934aed3..464d0a2c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -6,3 +6,11 @@
 def close_figures():
     yield
     plt.close("all")
+
+# conftest.py
+def pytest_collection_modifyitems(items):
+    for i, item in enumerate(items):
+        if item.name == 'test_model[tb]':
+            items.pop(i)
+            items.insert(0, item)
+            break

From 2a8bc1ef84cad49dcc082bbc78dcfc9df4f86b4e Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Fri, 7 Jun 2024 15:42:36 +1000
Subject: [PATCH 25/33] Formatting pass

---
 atomica/demos.py                       |   5 +-
 atomica/migration.py                   |   3 +-
 atomica/model.py                       | 132 ++++++++++++-------------
 atomica/parameters.py                  |  12 +--
 atomica/programs.py                    |   6 +-
 atomica/project.py                     |  35 +++----
 atomica/utils.py                       |   9 +-
 tests/conftest.py                      |   3 +-
 tests/test_tox_markovchain.py          |  43 ++++----
 tests/test_uncertainty_framework.xlsx  | Bin 27823 -> 28222 bytes
 tests/text_tox_timed_initialization.py |  29 +++---
 11 files changed, 135 insertions(+), 142 deletions(-)

diff --git a/atomica/demos.py b/atomica/demos.py
index baee7f61..c9e30de0 100644
--- a/atomica/demos.py
+++ b/atomica/demos.py
@@ -22,7 +22,7 @@
     "sir",
     "diabetes",
     "combined",
-    'service',
+    "service",
     "hypertension",
     "hypertension_dyn",
     "hiv",
@@ -33,7 +33,8 @@
     "dt",
 ]
 
-def demo(which: str = 'sir', do_run: bool = True, addprogs: bool = True) -> Project:
+
+def demo(which: str = "sir", do_run: bool = True, addprogs: bool = True) -> Project:
     """
     Return a demo project
 
diff --git a/atomica/migration.py b/atomica/migration.py
index 7c893a48..ea536eff 100644
--- a/atomica/migration.py
+++ b/atomica/migration.py
@@ -807,9 +807,8 @@ def _convert_framework_columns(framework):
     framework._validate()
     return framework
 
+
 @migration("ParameterSet", "1.28.3", "1.28.4", "Add initialization atttribute")
 def _parset_add_initialization(parset):
     parset.initialization = None
     return parset
-
-
diff --git a/atomica/model.py b/atomica/model.py
index a6b5ea6d..f664664d 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -230,7 +230,7 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """A class to wrap up data for one compartment within a cascade network."""
 
-    def __init__(self, pop, name, stochastic: bool = False, rng_sampler = None):
+    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
@@ -361,8 +361,8 @@ def resolve_outflows(self, ti: int) -> None:
         self._cached_outflow = 0
 
         if self.stochastic:
-            #for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
-            stays = max(0., 1. - outflow)
+            # for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
+            stays = max(0.0, 1.0 - outflow)
             rescaled_frac_flows = [rescale * link._cache for link in self.outlinks] + [stays]
             number_outflows = self.rng_sampler.multinomial(self[ti], rescaled_frac_flows, size=1)[0]
             for ln, link in enumerate(self.outlinks):
@@ -374,8 +374,6 @@ def resolve_outflows(self, ti: int) -> None:
                 link.vals[ti] = link._cache * n
                 self._cached_outflow += link.vals[ti]
 
-
-
     def update(self, ti: int) -> None:
         """
         Update compartment value
@@ -416,7 +414,7 @@ def connect(self, dest, par) -> None:
 
 
 class JunctionCompartment(Compartment):
-    def __init__(self, pop, name: str, duration_group: str = None, stochastic:bool = False, rng_sampler = None):
+    def __init__(self, pop, name: str, duration_group: str = None, stochastic: bool = False, rng_sampler=None):
         """
 
         A TimedCompartment has a duration group by virtue of having a `.parameter` attribute and a flush link.
@@ -534,14 +532,14 @@ def balance(self, ti: int) -> None:
         outflow_fractions = [link.parameter.vals[ti] for link in self.outlinks]
         total_outflow = sum(outflow_fractions)
 
-        outflow_fractions /= total_outflow #proportionately accounting for the total outflow downscaling
+        outflow_fractions /= total_outflow  # proportionately accounting for the total outflow downscaling
 
         # assign outflow stochastically to compartments based on the probabilities
-        if self.stochastic and net_inflow > 0.:
+        if self.stochastic and net_inflow > 0.0:
             if self.duration_group:
-                raise Exception('Need to implement code here to make sure outflows from EACH part of the duration group are integers')
+                raise Exception("Need to implement code here to make sure outflows from EACH part of the duration group are integers")
             else:
-                outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0]/net_inflow
+                outflow_fractions = self.rng_sampler.multinomial(net_inflow, outflow_fractions, size=1)[0] / net_inflow
 
         # Finally, assign the inflow to the outflow proportionately
         for frac, link in zip(outflow_fractions, self.outlinks):
@@ -571,7 +569,7 @@ def initial_flush(self) -> None:
 
             # assign outflow stochastically to compartments based on the probabilities
             if self.stochastic:
-                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0]/self.vals[0]
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / self.vals[0]
 
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
@@ -623,17 +621,17 @@ def balance(self, ti: int) -> None:
             outflow_fractions /= total_outflow
 
         # assign outflow stochastically to compartments based on the probabilities
-        if self.stochastic and net_inflow > 0.:
-            if np.abs(np.round(net_inflow)-net_inflow).sum() > 1e-3:
-                raise Exception(f'Net inflow should be an integer for a stochastic run, was: {net_inflow}')
+        if self.stochastic and net_inflow > 0.0:
+            if np.abs(np.round(net_inflow) - net_inflow).sum() > 1e-3:
+                raise Exception(f"Net inflow should be an integer for a stochastic run, was: {net_inflow}")
 
             # Give residual probability to the residual link if there is one, otherwise normalize probabilities
             if residual_link_index is not None:
                 outflow_fractions[residual_link_index] = max(0, 1 - outflow_fractions.sum())
-            else: # Probabilities have to sum to one as we have to flush the junction, and multinomials probabilities should sum to 1
+            else:  # Probabilities have to sum to one as we have to flush the junction, and multinomials probabilities should sum to 1
                 outflow_fractions = outflow_fractions / outflow_fractions.sum()
 
-            outflow_fractions = self.rng_sampler.multinomial(np.round(net_inflow).astype(int), outflow_fractions, size=1)[0]/(net_inflow)
+            outflow_fractions = self.rng_sampler.multinomial(np.round(net_inflow).astype(int), outflow_fractions, size=1)[0] / (net_inflow)
 
         # Calculate outflows
         outflow = net_inflow.reshape(-1, 1) * outflow_fractions.reshape(1, -1)
@@ -674,8 +672,8 @@ def initial_flush(self) -> None:
                 has_residual = False
 
             # assign outflow stochastically to compartments based on the probabilities
-            if self.stochastic and total_outflow > 0.:
-                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0]/(self.vals[0] * sum(outflow_fractions))
+            if self.stochastic and total_outflow > 0.0:
+                outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / (self.vals[0] * sum(outflow_fractions))
 
             # Assign the inflow directly to the outflow compartments
             # This is done using the [] indexing on the downstream compartment so it is agnostic
@@ -707,7 +705,7 @@ class SourceCompartment(Compartment):
 
     """
 
-    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
+    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
         super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
@@ -715,10 +713,10 @@ def preallocate(self, tvec: np.array, dt: float) -> None:
         self.vals.fill(0.0)  #: Source compartments have an unlimited number of people in them. Set to 0.0 to simplify validation.
 
     def resolve_outflows(self, ti: int) -> None:
-        #Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
+        # Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
         for link in self.outlinks:
             if self.stochastic:
-                link.vals[ti] = stochastic_rounding(link._cache, rng_sampler = self.rng_sampler)
+                link.vals[ti] = stochastic_rounding(link._cache, rng_sampler=self.rng_sampler)
             else:
                 link.vals[ti] = link._cache
 
@@ -727,7 +725,7 @@ def update(self, ti: int) -> None:
 
 
 class SinkCompartment(Compartment):
-    def __init__(self, pop, name, stochastic:bool = False, rng_sampler = None):
+    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
         super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
@@ -784,7 +782,7 @@ def update(self, ti: int) -> None:
 
 
 class TimedCompartment(Compartment):
-    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic:bool = False, rng_sampler= None):
+    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic: bool = False, rng_sampler=None):
         """
         Instantiate the TimedCompartment
 
@@ -901,8 +899,8 @@ def resolve_outflows(self, ti: int) -> None:
 
         """
 
-        #TODO adjust for stochastic variables
-        assert self.stochastic == False, 'resolve_outflows needs to be updated to work with stochastic transitions'
+        # TODO adjust for stochastic variables
+        assert self.stochastic == False, "resolve_outflows needs to be updated to work with stochastic transitions"
 
         # First, work out the scale factors as usual
         self.flush_link._cache = 0.0  # At this stage, no outflow goes via the flush link
@@ -1221,7 +1219,6 @@ def set_dynamic(self, progset=None) -> None:
             # because the values are modified during integration
             self._is_dynamic = True
 
-
         if self.deps:  # If there are no dependencies, then we know that this is precompute-only
             for deps in self.deps.values():  # deps is {'dep_name':[dep_objects]}
                 for dep in deps:
@@ -1590,7 +1587,7 @@ class Population:
     Each model population must contain a set of compartments with equivalent names.
     """
 
-    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool=False, rng_sampler=None):
+    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool = False, rng_sampler=None):
         """
         Construct a Population
 
@@ -1612,7 +1609,7 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         self.name = name  #: The code name of the population
         self.label = label  #: The full name/label of the population
         self.type = pop_type  #: The population's type
-        self.stochastic = stochastic #: Whether the population should be handled stochastically as discrete integer people instead of fractions
+        self.stochastic = stochastic  #: Whether the population should be handled stochastically as discrete integer people instead of fractions
         self.rng_sampler = rng_sampler
 
         self.comps = list()  #: List of Compartment objects
@@ -1824,17 +1821,17 @@ def build(self, framework, progset):
         for comp_name in list(comps.index):
             if comps.at[comp_name, "population type"] == self.type:
                 if comp_name in residual_junctions:
-                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "is junction"] == "y":
-                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "duration group"]:
-                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
+                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
                 elif comps.at[comp_name, "is source"] == "y":
-                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
                 elif comps.at[comp_name, "is sink"] == "y":
-                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
                 else:
-                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic = self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
 
         self.comp_lookup = {comp.name: comp for comp in self.comps}
 
@@ -2027,7 +2024,7 @@ def report_characteristic(charac, n_indent=0):
         # Otherwise, insert the values
         for i, c in enumerate(comps):
             if self.stochastic:
-                c[0] = stochastic_rounding(max(0.0, x[i]), rng_sampler = self.rng_sampler)#integer values (still represented as floats)
+                c[0] = stochastic_rounding(max(0.0, x[i]), rng_sampler=self.rng_sampler)  # integer values (still represented as floats)
             else:
                 c[0] = max(0.0, x[i])
 
@@ -2055,9 +2052,9 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.program_instructions = sc.dcp(program_instructions)  # program instructions
         self.t = settings.tvec  #: Simulation time vector (this is a brand new instance from the `settings.tvec` property method)
         self.dt = settings.sim_dt  #: Simulation time step
-        self.stochastic = settings.stochastic #: Run with discrete numbers of people per compartment instead of fractions.
+        self.stochastic = settings.stochastic  #: Run with discrete numbers of people per compartment instead of fractions.
 
-        self.acceptance_criteria = sc.dcp(acceptance_criteria) #: If the model fails to adhere to these criteria, raise a BadInitialization exception
+        self.acceptance_criteria = sc.dcp(acceptance_criteria)  #: If the model fails to adhere to these criteria, raise a BadInitialization exception
 
         if rng_sampler is None:
             rng_sampler = np.random.default_rng()
@@ -2288,8 +2285,6 @@ def build(self, parset):
                     for var in self._vars_by_pop[pars[0].pop_aggregation[3]]:
                         var.set_dynamic(progset=self.progset)
 
-
-
         # Set execution order - needs to be done _after_ pop aggregations have been flagged as dynamic
         self._set_exec_order()
 
@@ -2332,22 +2327,20 @@ def build(self, parset):
                 obj.preallocate(self.t, self.dt)
             pop.initialize_compartments(parset, self.framework, self.t[0])
 
-        #More finally, set up acceptance criteria for those which should be checked at runtime (variables that have ._is_dynamic = True), and those which can only be checked after conclusion
+        # More finally, set up acceptance criteria for those which should be checked at runtime (variables that have ._is_dynamic = True), and those which can only be checked after conclusion
         if self.acceptance_criteria:
             for criteria in self.acceptance_criteria:
-                pars = self._vars_by_pop[criteria['parameter']]
+                pars = self._vars_by_pop[criteria["parameter"]]
                 for par in pars:
-                    if par.pop.name == criteria['population']:
+                    if par.pop.name == criteria["population"]:
                         if par._is_dynamic:
-                            criteria['evaluate_at'] = self.t[np.where(self.t > criteria['t_range'][1])][0]
+                            criteria["evaluate_at"] = self.t[np.where(self.t > criteria["t_range"][1])][0]
                         else:
-                            criteria['evaluate_at'] = np.inf #can only evaluate postcompute
-
-            #sort by evaluation time - now we only have to look at the first element every timestep
-            self.acceptance_criteria = sorted(self.acceptance_criteria, key=lambda item: item['evaluate_at'])
-            self.acceptance_index = 0 #which index to evaluate
-
+                            criteria["evaluate_at"] = np.inf  # can only evaluate postcompute
 
+            # sort by evaluation time - now we only have to look at the first element every timestep
+            self.acceptance_criteria = sorted(self.acceptance_criteria, key=lambda item: item["evaluate_at"])
+            self.acceptance_index = 0  # which index to evaluate
 
     def _set_exec_order(self) -> None:
         """
@@ -2482,8 +2475,8 @@ def process(self) -> None:
                     par.update()
                     par.constrain()
 
-        if self.acceptance_criteria: #call update_acceptance once more for any variables that could only be assessed postcompute
-                self.update_acceptance(evaluate_all = True)
+        if self.acceptance_criteria:  # call update_acceptance once more for any variables that could only be assessed postcompute
+            self.update_acceptance(evaluate_all=True)
 
         # Clear characteristic internal storage and switch to dynamic computation to save space
         for pop in self.pops:
@@ -2492,7 +2485,6 @@ def process(self) -> None:
 
         self._program_cache = None  # Drop the program cache afterwards to save space
 
-
     def update_links(self) -> None:
         """
         Update link values
@@ -2613,7 +2605,7 @@ def update_comps(self) -> None:
             for comp in pop.comps:
                 comp.update(ti)
 
-    def update_acceptance(self, evaluate_all = False) -> None:
+    def update_acceptance(self, evaluate_all=False) -> None:
         """
         Update any acceptance criteria for the model run, raise BadInitialization error if the model run is failing
 
@@ -2621,33 +2613,31 @@ def update_acceptance(self, evaluate_all = False) -> None:
 
         """
 
-        time = self.t[self._t_index] #actual time
+        time = self.t[self._t_index]  # actual time
 
-        #We know that these are sorted by evaluation time so we only need to look at the first index
-        while (self.acceptance_index < len(self.acceptance_criteria) ) and (evaluate_all or time >= self.acceptance_criteria[self.acceptance_index]['evaluate_at']):
-            #time ranges for assessment are in order, so only need to check the first one
+        # We know that these are sorted by evaluation time so we only need to look at the first index
+        while (self.acceptance_index < len(self.acceptance_criteria)) and (evaluate_all or time >= self.acceptance_criteria[self.acceptance_index]["evaluate_at"]):
+            # time ranges for assessment are in order, so only need to check the first one
             assessable = self.acceptance_criteria[self.acceptance_index]
-            #we have passed the time window, assess if the criteria was met - either return a fail condition or remove this acceptance criteria
-            a_par = assessable['parameter']
-            a_pop = assessable['population']
-            a_times = assessable['t_range']
-            accept_vals = assessable['value']
-            a_t_inds = np.where(np.logical_and(self.t>=a_times[0], self.t<=a_times[1])) #self.t is the full tvec for the model
+            # we have passed the time window, assess if the criteria was met - either return a fail condition or remove this acceptance criteria
+            a_par = assessable["parameter"]
+            a_pop = assessable["population"]
+            a_times = assessable["t_range"]
+            accept_vals = assessable["value"]
+            a_t_inds = np.where(np.logical_and(self.t >= a_times[0], self.t <= a_times[1]))  # self.t is the full tvec for the model
 
             pars = self._vars_by_pop[a_par]
-            for par in pars: #TODO should be more efficient way to directly access the par values for the pop rather than looping
+            for par in pars:  # TODO should be more efficient way to directly access the par values for the pop rather than looping
                 if par.pop.name == a_pop:
-                    a_t_val = np.mean(par[a_t_inds]) #mean across the time range
+                    a_t_val = np.mean(par[a_t_inds])  # mean across the time range
                     if a_t_val <= accept_vals[0] or a_t_val >= accept_vals[1]:
                         # print(f'{time} Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
-                        raise BadInitialization(f'Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+                        raise BadInitialization(f"Rejecting run as sampled sim value {a_t_val} outside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}")
                     # else:
-                        # print(f'{time} ACCEPTING run as sampled sim value {a_t_val} inside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
-            assessable['assessed_value'] = a_t_val
+                    # print(f'{time} ACCEPTING run as sampled sim value {a_t_val} inside of acceptable bound {accept_vals} for parameter {a_par}, population {a_pop} at time {a_times}')
+            assessable["assessed_value"] = a_t_val
             self.acceptance_index += 1
 
-
-
     def flush_junctions(self) -> None:
         """
         Flush initialization values from junctions
@@ -2784,7 +2774,7 @@ def update_pars(self) -> None:
                     par.constrain(ti)
 
 
-def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler = None, acceptance_criteria = []):
+def run_model(settings, framework, parset: ParameterSet, progset: ProgramSet = None, program_instructions: ProgramInstructions = None, name: str = None, rng_sampler=None, acceptance_criteria=[]):
     """
     Build and process model
 
diff --git a/atomica/parameters.py b/atomica/parameters.py
index e6347f62..d3221497 100644
--- a/atomica/parameters.py
+++ b/atomica/parameters.py
@@ -341,16 +341,16 @@ def from_excel(cls, excelfile: pd.ExcelFile):
         """
         # excelfile = spreadsheet.pandas()
 
-        metadata, value_df = atomica.excel.read_dataframes(excelfile.book['Initialization'])
+        metadata, value_df = atomica.excel.read_dataframes(excelfile.book["Initialization"])
 
         values = {}
-        for k,s in value_df.T.reset_index().T.set_index([0,1]).iterrows():
+        for k, s in value_df.T.reset_index().T.set_index([0, 1]).iterrows():
             v = s.dropna().values
             values[k] = v[0] if len(v) == 1 else v
 
         self = cls(values)
 
-        for k,v in metadata.T.reset_index().T.set_index(0)[1].to_dict().items():
+        for k, v in metadata.T.reset_index().T.set_index(0)[1].to_dict().items():
             setattr(self, k, v)
 
         return self
@@ -509,7 +509,7 @@ def get_par(self, name: str, pop: str = None) -> Parameter:
         else:
             raise KeyError(f'Parameter "{name}" not found')
 
-    def sample(self, constant:bool = True, rng_sampler = None):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Return a sampled copy of the ParameterSet
 
@@ -522,7 +522,7 @@ def sample(self, constant:bool = True, rng_sampler = None):
         new = sc.dcp(self)
 
         for i, par in enumerate(new.all_pars()):
-            par.sample(constant = constant, rng_sampler = rng_sampler)
+            par.sample(constant=constant, rng_sampler=rng_sampler)
         return new
 
     def make_constant(self, year: float):
@@ -535,7 +535,7 @@ def make_constant(self, year: float):
         :param year: Year to use for interpolation
         :return: A copy of the ParameterSet with constant parameters
         """
-        ps = self.copy(f'{self.name} (constant)')
+        ps = self.copy(f"{self.name} (constant)")
         for par in ps.all_pars():
             for ts in par.ts.values():
                 ts.insert(None, ts.interpolate(year))
diff --git a/atomica/programs.py b/atomica/programs.py
index 83d26d4b..8e3edad8 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1152,7 +1152,7 @@ def get_outcomes(self, prop_coverage: dict) -> dict:
 
         return {(covout.par, covout.pop): covout.get_outcome(prop_coverage) for covout in self.covouts.values()}
 
-    def sample(self, constant: bool = True, rng_sampler = None):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Perturb programs based on uncertainties
 
@@ -1409,7 +1409,7 @@ def n_progs(self) -> int:
 
         return len(self.progs)
 
-    def sample(self, rng_sampler = None) -> None:
+    def sample(self, rng_sampler=None) -> None:
         """
         Perturb the values entered in the databook
 
@@ -1421,7 +1421,7 @@ def sample(self, rng_sampler = None) -> None:
 
         if self.sigma is None:
             return
-        
+
         if rng_sampler is None:
             rng_sampler = np.random.default_rng()
 
diff --git a/atomica/project.py b/atomica/project.py
index 7fb27cf4..29c6f6ed 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -87,9 +87,9 @@ def sim_dt(self, sim_dt):
         assert sim_dt > 0, "Simulation timestep must be greater than 0"
         self._sim_dt = sim_dt
         self.sim_end = self.sim_end  # Call the setter function to change sim_end if it is no longer valid
-        
+
     @stochastic.setter
-    def stochastic(self, stochastic:bool=False):
+    def stochastic(self, stochastic: bool = False):
         self._stochastic = stochastic
 
     @property
@@ -484,7 +484,7 @@ def update_settings(self, sim_start=None, sim_end=None, sim_dt=None):
         """Modify the project settings, e.g. the simulation time vector."""
         self.settings.update_time_vector(start=sim_start, end=sim_end, dt=sim_dt)
 
-    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng = None, acceptance_criteria=[]):
+    def run_sim(self, parset=None, progset=None, progset_instructions=None, store_results=False, result_name: str = None, rng=None, acceptance_criteria=[]):
         """
         Run a single simulation
 
@@ -500,7 +500,7 @@ def run_sim(self, parset=None, progset=None, progset_instructions=None, store_re
         :param rng: Optionally a random number generator that may have been seeded to generate consistent results, or a random_seed used to generate a Generator
         :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters over a given period of time
         tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't_range': (2018, 2018.99)}]
-        
+
         :return: A :class:`Result` instance
 
         """
@@ -558,7 +558,7 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
         :param acceptance_criteria: Optional criteria to assert that outputs are within a given tolerance of a given value for given parameters in given years
             tolerance e.g.  = [{'parameter': 'incidence', 'population': 'adults', 'value': (0.02, 0.025), 't': 2018}]
 
-        
+
         :return: A list of Results that can be passed to `Ensemble.update()`. If multiple instructions are provided, the return value of this
                  function will be a list of lists, where the inner list iterates over different instructions for the same parset/progset samples.
                  It is expected in that case that the Ensemble's mapping function would take in a list of results
@@ -581,19 +581,19 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu
             assert (len(result_names) == 1 and not progset) or (len(progset_instructions) == len(result_names)), "Number of result names must match number of instructions"
 
         show_progress = n_samples > 1 and logger.getEffectiveLevel() <= logging.INFO
-        
+
         if rand_seed is not None:
-            rng = np.random.default_rng(seed = rand_seed)
+            rng = np.random.default_rng(seed=rand_seed)
             seed_samples = rng.integers(1e15, size=n_samples)
         else:
-            seed_samples = [None]*n_samples
-            
-        model_rngs = [np.random.default_rng(seed = seed) for seed in seed_samples] #generate a RNG for each model
+            seed_samples = [None] * n_samples
+
+        model_rngs = [np.random.default_rng(seed=seed) for seed in seed_samples]  # generate a RNG for each model
 
         if parallel:
             fcn = functools.partial(_run_sampled_sim, proj=self, parset=parset, progset=progset, progset_instructions=progset_instructions, result_names=result_names, max_attempts=max_attempts, acceptance_criteria=acceptance_criteria)
-            #as multiprocessing does not handle partial functions as compiled functions, need to send the rngs as kwargs in a dictionary, not as args to the partial function
-            model_rng_kwargs = [{'rng_sampler': rng} for rng in model_rngs]
+            # as multiprocessing does not handle partial functions as compiled functions, need to send the rngs as kwargs in a dictionary, not as args to the partial function
+            model_rng_kwargs = [{"rng_sampler": rng} for rng in model_rngs]
             results = parallel_progress(fcn, model_rng_kwargs, show_progress=show_progress, num_workers=num_workers)
         elif show_progress:
             # Print the progress bar if the logging level was INFO or lower
@@ -743,7 +743,8 @@ def __setstate__(self, d):
         P = migrate(self)
         self.__dict__ = P.__dict__
 
-def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler =  None, acceptance_criteria = []):
+
+def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_names: list, max_attempts: int = None, rng_sampler=None, acceptance_criteria=[]):
     """
     Internal function to run simulation with sampling
 
@@ -775,11 +776,11 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
 
     if max_attempts is None:
         max_attempts = 50
-    
+
     if rng_sampler is None:
         rng_sampler = np.random.default_rng()
-    
-    #Set up separate seeds for each attempt and seed a rng for each one (in case scenarios diverge in later timepoints, it's important to START from the same seeds)
+
+    # Set up separate seeds for each attempt and seed a rng for each one (in case scenarios diverge in later timepoints, it's important to START from the same seeds)
     rand_seeds = rng_sampler.integers(1e15, size=max_attempts)
 
     attempts = 0
@@ -792,7 +793,7 @@ def _run_sampled_sim(proj, parset, progset, progset_instructions: list, result_n
                 results = [proj.run_sim(parset=sampled_parset, progset=sampled_progset, progset_instructions=x, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for x, y in zip(progset_instructions, result_names)]
             else:
                 results = [proj.run_sim(parset=sampled_parset, result_name=y, rng=rng_attempt, acceptance_criteria=acceptance_criteria) for y in result_names]
-            
+
             return results
         except BadInitialization:
             attempts += 1
diff --git a/atomica/utils.py b/atomica/utils.py
index 2c7ce036..87b711b1 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -587,7 +587,7 @@ def interpolate(self, t2: np.array, method="linear", **kwargs) -> np.array:
         interpolator = method(t1, v1, **kwargs)
         return interpolator(t2)
 
-    def sample(self, constant:bool = True, rng_sampler = None):
+    def sample(self, constant: bool = True, rng_sampler=None):
         """
         Return a sampled copy of the TimeSeries
 
@@ -847,7 +847,7 @@ def callback(result, idx):
     else:
         for i, x in enumerate(inputs):
             if isinstance(x, dict):
-                pool.apply_async(fcn, args=(), kwds=x, callback=partial(callback, idx=i)) #if the list of inputs is given as a dictionary then pass as kwargs (kwds) instead
+                pool.apply_async(fcn, args=(), kwds=x, callback=partial(callback, idx=i))  # if the list of inputs is given as a dictionary then pass as kwargs (kwds) instead
             else:
                 pool.apply_async(fcn, args=(x,), callback=partial(callback, idx=i))
 
@@ -991,7 +991,8 @@ def stop_logging() -> None:
             # Don't terminate the loop, if by some change there is more than one handler
             # (not supposed to happen though) then we would want to close them all
 
-def stochastic_rounding(x, rng_sampler = None):
+
+def stochastic_rounding(x, rng_sampler=None):
     """
     Stochastically round a float up or down to an integer-equivalent value (note: returns still as a float)
     :param x: value to be rounded
@@ -1001,5 +1002,5 @@ def stochastic_rounding(x, rng_sampler = None):
     floor = np.floor(x)
     remainder = x - floor
     if remainder:
-        x = floor + int(sample<remainder)
+        x = floor + int(sample < remainder)
     return x
diff --git a/tests/conftest.py b/tests/conftest.py
index 464d0a2c..3439bfe3 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -7,10 +7,11 @@ def close_figures():
     yield
     plt.close("all")
 
+
 # conftest.py
 def pytest_collection_modifyitems(items):
     for i, item in enumerate(items):
-        if item.name == 'test_model[tb]':
+        if item.name == "test_model[tb]":
             items.pop(i)
             items.insert(0, item)
             break
diff --git a/tests/test_tox_markovchain.py b/tests/test_tox_markovchain.py
index a06ae9a7..e5d723bd 100644
--- a/tests/test_tox_markovchain.py
+++ b/tests/test_tox_markovchain.py
@@ -6,18 +6,19 @@
 import pytest
 
 
-
 dirname = at.parent_dir()
 # np.seterr(all='raise')
 
+
 @pytest.mark.parametrize("demo", at.demo_projects)
 def test_functional_all_framework_databook(demo):
     # Test that all demo frameworks can be run with stochastic mode turned on
     P = at.demo(demo, do_run=False)
     P.settings.stochastic = True
-    P.run_sim(result_name='DTMC')
+    P.run_sim(result_name="DTMC")
+
 
-@pytest.mark.parametrize("demo", ['sir','udt','tb_simple']) # This test is expensive so only run it for a few lightweight models
+@pytest.mark.parametrize("demo", ["sir", "udt", "tb_simple"])  # This test is expensive so only run it for a few lightweight models
 def test_compare_deterministic_markovchain(demo):
 
     seed = 0
@@ -29,8 +30,6 @@ def test_compare_deterministic_markovchain(demo):
 
     ratio_cutoff = 0.05
 
-
-
     def filter_points(t, det, mc, mc_mean):
         # t_fil, det_fil, mc_fil = sc.dcp(t), sc.dcp(det), sc.dcp(mc)
 
@@ -52,12 +51,13 @@ def get_array_outputs(results, output, pop):
 
     def get_average_output(results, output, pop):
         t, array = get_array_outputs(results, output, pop)
-        average = np.mean( array , axis=0)
+        average = np.mean(array, axis=0)
 
         return t, average, array
 
     def do_hypothesis_test(det_array, mc_results):
         import scipy
+
         res = scipy.stats.ttest_1samp(mc_results, det_array, keepdims=True)
         return res
 
@@ -71,8 +71,10 @@ def run_simple_percent_test(P, baseline_results, mc_results):
 
         for compartment_name in compartment_names:
             for pop in baseline_results[0].pop_names:
-                try: baseline_results[0].model.get_pop(pop).get_variable(compartment_name)[0].vals
-                except: continue  # This population doesn't have this compartment
+                try:
+                    baseline_results[0].model.get_pop(pop).get_variable(compartment_name)[0].vals
+                except:
+                    continue  # This population doesn't have this compartment
 
                 t, mc_mean, mc_array = get_average_output(mc_results, compartment_name, pop)
                 _, det_mean, _ = get_average_output(baseline_results, compartment_name, pop)
@@ -85,24 +87,21 @@ def run_simple_percent_test(P, baseline_results, mc_results):
                 indices = det_mean != 0
                 ratio = np.ones(mc_mean.shape)
                 ratio[indices] = mc_mean[indices] / det_mean[indices]
-                if min(ratio) < 1-ratio_cutoff or max(ratio) > 1+ratio_cutoff:
-                    print('RATIO', compartment_name, pop, min(ratio), max(ratio), list(zip(ratio, det_mean, mc_mean)))
-                    indices = np.logical_or(ratio < 1-ratio_cutoff, ration > 1+ratio_cutoff)
+                if min(ratio) < 1 - ratio_cutoff or max(ratio) > 1 + ratio_cutoff:
+                    print("RATIO", compartment_name, pop, min(ratio), max(ratio), list(zip(ratio, det_mean, mc_mean)))
+                    indices = np.logical_or(ratio < 1 - ratio_cutoff, ration > 1 + ratio_cutoff)
 
-                    raise Exception(f'The Markov Chain runs (N={len(mc_results)} differed too much from the baseline runs (N={len(baseline_results)} at starting at timestep={indices[0]}, t={t[indices[0]]}, with compartment={compartment_name}, population={pop}. Baseline={det_mean[indices[0]]}, MC={mc_mean[indices[0]]}, difference={100*(mc_mean[indices[0]]/det_mean[indices[0]] - 1):.2f}')
+                    raise Exception(f"The Markov Chain runs (N={len(mc_results)} differed too much from the baseline runs (N={len(baseline_results)} at starting at timestep={indices[0]}, t={t[indices[0]]}, with compartment={compartment_name}, population={pop}. Baseline={det_mean[indices[0]]}, MC={mc_mean[indices[0]]}, difference={100*(mc_mean[indices[0]]/det_mean[indices[0]] - 1):.2f}")
 
                 max_diff = np.max(np.abs(ratio - 1))
                 max_max_diff = max(max_diff, max_max_diff)
 
-        print(f'Success: max difference {max_max_diff*100:.2f}%')
-
-
-
-    seed = np.random.default_rng(seed).integers(2**32-1)
+        print(f"Success: max difference {max_max_diff*100:.2f}%")
 
+    seed = np.random.default_rng(seed).integers(2**32 - 1)
 
     P = at.demo(demo, do_run=False)
-    res = P.run_sim(result_name='Original')
+    res = P.run_sim(result_name="Original")
 
     baseline_results = [res]
 
@@ -110,13 +109,13 @@ def run_simple_percent_test(P, baseline_results, mc_results):
     mc_results = [None] * N
     with at.Quiet():  # Using at.Quiet should automatically reset the logging level if an exception occurs
         for i in range(N):
-            this_seed = ( i + seed + hash(demo) )% (2**32)
-            mc_results[i] = P.run_sim(result_name=f'DTMC {i}', rng=this_seed )
+            this_seed = (i + seed + hash(demo)) % (2**32)
+            mc_results[i] = P.run_sim(result_name=f"DTMC {i}", rng=this_seed)
 
     run_simple_percent_test(P, baseline_results=baseline_results, mc_results=mc_results)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     for demo in at.demo_projects:
         test_functional_all_framework_databook()
-    test_compare_deterministic_markovchain('sir')
+    test_compare_deterministic_markovchain("sir")
diff --git a/tests/test_uncertainty_framework.xlsx b/tests/test_uncertainty_framework.xlsx
index 120b173526f99fec154d7b0a722db3f2aa266018..37073ea5173f9baace0cff63e594fdd697ef1536 100644
GIT binary patch
delta 13627
zcmbVzWl$Ykw=ELf-Q8V-ySoN=cX!>m`^GIe1PkuLVM7S+?g4@aC&=aG-1_o;=iYkt
z-rIk=tJkR6y}IUDW31V8e_jUfUIVWqg#knhH^t^cf`RG5fq|ibfq{A3vv@hV+L<~z
z*)e-N*q3NJIPLJDz68#^MAB}~4*EbA&0wKUY1nF2jJ0^7t630=`z!9`o?=$czuXbQ
zq(JlJEipszw^$KRj&X?B3ZMV{=3OerZiuMpzMl!r1gorFkC(lP%~9oqueMFzhX-Ie
zoboT8dGfzHN<3!5VPTjsLIk7Tn^H)lr&vv=Eagqt>MDWCt!^Kw4=G(*YQV&#VK4o#
zpLNQoWPOl{^_`duJ;V6oCL0akO>MtP#Sc}Kqq^5Wt3i#+nz#KOGiaxjo~i!|om7+2
z>_-js`}kam)|?Q8uCKkPMDTHWK~4Zz;l181F7CeRWer}65czIREii}kc%ds2UM?|d
zvZXxG7CYr8_I273VT)$;MW|*tljc2gMDC`PxJXSl1KPXuP-{`U9J>wu&pJ{hcXRKZ
zbprd^Cq7u><E7hXvq}OBmqV#GXNr+j({^JloPXSkwK}bEGvyuh%N*isW$6OIk_o&k
z5Es4HRpcSqq_T>NdSjHl3xgtIdIB@t-rrL0_d4$ezv!`v7S;tBZ85}UN_bNpXqQS+
zH5Bqjr*BfW^<WB=vm2y`yoyJs$XD}BL9<|aGi4Cv?_3q#LzS5pt+)2Qeu@8eFsmK(
ziQIOHN2@>RVD%oE{U=tBSYQwUYP8Vbi41(ik2xx8>gp~y@iEjpgb_r5#0VhtUZ9oA
zpHJKkL{&DBnZJhNq3Dlt*8g;!<clQ888C9SJT18Ko@OMPKD{gFa-G+oQSc+;_%h0A
zoug}FoZ~*f?Dfo&U_JS^o->*I<JP=hqvrrIW9!WkhWrP3P6PuM>)RJVi9Il6bmTm|
zgF38Aksw=|ddLrGc~O7uY@V`S|LLsUSWl2Z(g3m;>oexH!@BJ;X--7cHP#N7=t_C5
za~xJJ!2_XFL3@v#Xrx{9ntDyNW0u$v^V<9+AK%8O39&7)h7Jv{+Dm@=m*w|#*qI{D
zY@cwCHr{R9*B~g(e=kS_sH+^?ityXL{~>oKWfDTwH=6v~TcO#lrQQ?$z7H>0s<51|
z3v#awl1bF-hf<|)1inSmtjiP$|3+`~6JEBC5I)a>H{qiu>SR|z<jyE(rj6#Y@U(R~
z+9JtGd_d<ynYP9hxBFPdVluH+<hrNZmrb6F{R%h<lzPwof@eSoztJ`oSuqa$BrywS
zu7A7Hnh~uZw}F80I5QUUdgAW%{?Fpb$mz}4BcpMwTLAQ2M6NiEJ-ZY4xg@tR5RD3}
zu?yZ20^WeLnknIp#K~|#Rw_II+h5Gh1<EI`lA?mn(;kUX?{IeoxaxPRJ9H%1S1VUR
zc)Kofrf~>+rInjsLi0d!M`XI@kP%s%PcXw*pwQY;<ZAV@QQAaPxq{>}m2tje;qK4S
z>T<a_ie$F4B2tN4hr2sh<qsQ6@l{{~qGG%3x;~p_alAUg*v#xqoPaiU{Hkz}a;If&
zGil6(>>^|IF5Ilx`zD%_TjB^hd8+N@k(!9%J5+?`i;7K+N)M~-^t-L5x*P}oB&pV@
z5sjLs`cSnZi^m4XI4qxdACh|+8P~dd%6sHVla8Al8C_YsT)KcT^98$ZMZmD7!U)0<
z*J~%M{~+VZhqhh@dcZ>_rBHWW3WNQ<=+Yhp0?1;vhMxVAMyMG)LBRd8ng9H`G}cFp
z=Cdc;9l25Vini&ZkCfI@y0!GO=>#)Zc)u&Gv<PD=-gEyfjSwPr-jFhs-<_ip>I3w9
z8I+1Bb@UisZc7YABWF6~Oa9Ox%$X`8#==^T$FNUrn+B-4h5*iYpvL^WX-$KN9$hsN
zx|U{;MBt}>ExbIb*-LNEETOm~pR1f@0913WfNlSLZPU!oJ))9nO>reHUS!hDTqe@3
zMAgNtX-$myN9jo=Z%(;k6!nB#3CBF{N6fU3Pkzc>yAv3ofOB5IsC(M{cRhr;MxP4y
z0f-CH_b^~!uYlLLb^P@;HFX)epA=ORP`V9iUmmL=o@4EQ!HrM(BUHY`5pZsd4~>em
zf&bn4oV8uOwngoBtoB?tP-r8?Z6Jv1#~cx6;J9o0TH(o{{j;_l)n1{KHrH2%1=qsm
z_6ULqTx`sI!_~O*#q@GhQtb10xa5-ZyRqag=h#?RMF7HiEcrpn8_=`d8P_8wRu6Y+
z8h*r`f_nSn9KUK<4!v1Kp@|o+g*KMnFTmg}e`TE^QE|CR>zIqS2YPW#U6$vm;%cp8
zzu_z{OIZBk+*N~}h>m*-w9bdYWG+vU?Et1;F+0^X*YEhlr<gO8vVDYzc|Rh9f?E}Y
z)0Z}k3IK6P5lT1Wngvn#bTijt!pyp{m4U`;k7k14f4!IAFNZTu6uw^tQQ|YYD;Psm
zn~fi&mZ0e1^=rXPn7hQqb&HtSRRs5JB<2@6Y&!%&kKta>Ksa5^oZ!IjW`33<^1WSo
zrR(C_;UH(HS)}7h>!$1_GCbV)Xm>}9LOy0062ujkvp4xr?{UE|Z=v+cm`fO-s1znV
zU3|KN=o?SP{^lwBZ=R;gC+4vJ%~QtBP&R@VIC!0A0?E$S${)KYnPl4aY|VKYm4SB&
zqh{*+A=tqn){C6iC5~_T$WlNep3nArYWXFpA`A|tD;}}{f1?EI)O*quIqOw1*~G1*
z*{9ve<&BPbCol+6v19p_uqmuKpXDGlp?abwz~QGF4)U!UJ&sLYRE^i<X-4^3Y0EI%
zqiB^jGDN*ax&HpIHNjlDq{vOM6`S0ZS2!z$_=8osOJ=n=g52mhG?Qr{kk?~e%7m|8
zq=fE((#%<7Zo?742cMTxm150uH6C@%#Fe2`gXi{0&D8v|!?6Z~A)eoL%ZzNaZ#!!n
zaN|s&y(W#%tTV!J6i5YQa-Ye&^hsOSsp2G{YUGN6!7o&>3***1d}!ZwNK=FA8JM4l
zrb`^@p8At#jSz6;uSr+^B@@d<CF7cdrQ2C)*?%G9_Km9m4w`y(%>~DzU98=bsEY#8
z;cP01apZi_<CJcN`sKhom0a>9gXmVk6+*!B-GY`heZ`Pwq*WE*6Fuh=NC403p&V|~
z)%P;6@fVKDF?9jK$6E{a4<IP-=)AB7bbmZExhgiWgu#k?`Cea)$c?kf3)vJ!*L-2A
zdw75C%8?QY$AG%;<lEo)RGlY%;a51Q#$eR7S&m44C*FTkQxrLb6@V)FdFC1{T~mdu
zG&&-dvpJf9Gqqr21DSS$7NNr~kmQ^LE)*H({MYd9!r_^(mBBcb)wbm-s&d-S37Ok*
z<yGTl-$To1pB7ap?iFyqci#S#hGCh=Qt7&gQV4faO0Q90EhPqL+Id@QoS@sj^IEwz
zrU;KCMCm$ilyO2Ljstc$$lyl=F83EEAzAB-g0y2X8F~j}AWYXwwZW-1F;S0&hTwk>
zpwMXsMPoq?>X!XIQ1wl%i9)f9r*Nl9Fs}tQpb4D{sO|p>#qwv04RkUJ$~TIPa+BW}
zG)H6=8R|9Ir#Q6KXum2<2Gz^|MCdJ&ER%zmeW0RxXlFH)vNe^nMT@ekE=fxNXauo)
zxVpSoQbt3q)qzQsZuPyC0aM4KY>6Jrrwo@KeE)jvDz+h(Lrq=&1Pmx5g$1uVZwyZR
z4h#$e>9+uK{7+&*W<rUnj8v}B{@AD#noMBpgCCyrG^}}C=f{;%+ZW@O!=;&c`Yw#(
z_u!xYp3l`G%+^(*3yfleRU^tWW7o-uyAhMZ%tn1m_q_lJ$MNFo$8+bt?9rD`Jgu;>
z@uMyl4ml6qCRiDB-j(}EfRj`?dE4!(jfb-NT<i(ig+SbQ(JY^=os~zfF(c^Z`|fq?
zL-*oO$+8#cLITyJ7vh@S*rK(peM9u61-;N0pC{6!i(mrOSZm_wYe1>jem03Z%uhqY
zK5|@Cs<JmPQxjloVy=3s2!1X{rj$Cj9d@Z@tN47b)Myrm)!Rh<fP*6?)_^!p19l6o
z88a#x6PZ1`K+JXPb(QVaN|_x7hOWN+4Oscwx>|ic{wQ)#Pw4i0Ncj$FRs>_Yb=$RH
zp;}{Ny(b<zWaDFm0McdZz8Bjg>R+8(G;Xa@J*g>1xTm_3DesEUw(QeZya4{;b>OH-
zpLvg^a8+5QU;2;v03ZSgtMl^M`e0njs2T4w@`pv0s(lB+i1jv(JK5rpjBAxu9F14h
z@EX@#UmHScS-!-}wg;U&A=mmU#*t8?bxw;Ri*ck-Na&wktY;}^mC!;0$(ZYk;AS$D
z6ZyVe!Bu8s@#|5@;b{E0RePwJh*2AOe22Mx2@`{{{oeWL09dXxh#EJ{I!LIq%4qna
zKAe*+g3UBb@^>YC;Knc`{Rp#}gQEa*aExuzb#ugUM4H8ZFM|h!f2l+2oStJSD{A<5
z4Jk_ja%)pGq8i(+5US)?b=@IaMNcFC%0eL`n|^|j#7Ka@>KUh^oA5XBD&~)y%cW@T
zBWr*~XU~_THk@Uq&k>jRF6T7dN(b#|of$>gHTa<cDEF7KuUk*6c6zCDhI&*WqZzak
zxoB#cyp(*>U#YkuU1fvldk$nXD8^Z0^%n5p+(qcieBqD^+ro5StWY6)tWf#5`Aae-
zMrA$Y`|00<X1_NCS%&m}pL*A4`Fo@eu*|Y@)hTFf9Z|D}_O7LY8kXy7@M2IZWc=n0
zUx<fnD2>YFolUy*4PyX<0FLWd&dGxZGH&*w-Hywp$Gc#P&u}5*0&D{H$&BN!U*-XS
zSgoJ<#c_ZYG0aGEYi&+gRIfLwxe1Tme7ttJf$;;HK}>|%JD<A|fP=x@zVlvy{+W5P
z4oS0@4L_o4t~0E@IJjJ`QAqj5+?IhkYW;u<$zw1IwUh$}OZRpHtPPa#Xw&<y^!YQq
z!SKqy_wV(uJ}#FSOr$vdbWIm_lgHC7x!Uq5Bo$B~Z8TzO^boiqk)ZmOdQvdjF*laO
zHUE{wNF=Fo-Ie3kh~*k?ke_|m$0ct%zHr$=v?;19Q$g0LT1e5UL>i~Mz}8XXS6R^Y
zm{dvP65hZM<lowSh-5m5{-({ou)miO&c8kMYQmc~N1QX<5HW7TK*DjxifPU}7M~f+
z51cidMHmmk*z%Wa2HX(=sR`(~r*(}92)Y5oR~((2c9JbyT}@?WIw@)pQcF_V8j65t
z{kXE!`(DL?wD&5rklwz+ulFV`jgF`MK_W@U{)(Nbk1=d4d{b*QEij<~{u+;L7=L$$
z?|v{{ctSHV`!>4l<E4etV*8FJOD$;-#5GIStlO|I1>U{)0W^MqA{!S3CAZd%YYG#O
z#Hi8_?1a#?RJtFO8O(@`=)d`WRBR*li_}WEyQk!+Hlnkk%e`qqf=gLu`+j_VzC^UU
zqtG$&O0B?Ug51K6eSwq)=wT_T=4t%Ny>y>f_Y*M|v|kx-Z~MsoyeY+FSuZbr#N&x`
zN|HB;AB~lug9hWEC~P%fRAMhg@N~s!a=`L^FR}SCa;yG>HKSaSV8%l2ZmcqG4*7H_
z1X{cp53-&&SgGw}sYK)}F{g73#7VU?uka@(XmNmBNpQN;2k*Kc07<UGsw2)U?4gp}
zK>9<QtIqWRsP-QaFS>1fnWeVKiCa(|IK^^5A*hW0n7WeGP2><#%t)gl?ZxER4(mZn
z)APQO1;v%1>W%JDHl%Wl#k|?gNRk6W(jU6%*)HRB$u}HHFqlOm>prDhHa+5f3<_;a
zA2+K%6alKQ!5z3&V`$N+{W!O=gHhz<y#a+j0iAB2oRT!FJX#W9yH#eUq&))B)>J~m
zhula6TH4TOO0!V%w0-Z7`(c!rs$n*0NW85~RJHQBBk;Y#6R<J)IJSl1k=l+>ugeQZ
zJ=bUqcg*_Tgl6*~I~L(U9w}6O)B@!kyzpRmTmaS@Pl-IsQh9p??+HuW7Tn%1NU0Zh
zy}8{C{QY$`Y$iM%8z3HzZQkvAA4e8&=<5=fOM?}*h#LOkww^GdKaBmq(Bc10hOz(X
z@L|RQjs=cYNm=?5_+MfRoWg$xo&Ls@9{9hQ0%Wm*0_3Lt2T@pq?L*#(LalmN`{Qzj
zjrzBTRj=yy1BoDOl)y8|*@tj!zD|07dwM9L<a0I)uXR%hfu%6yFo^Qac@x{kVmhgX
z>5c-mMFEhN<K@-PbLWF<PF@X<7c6Y^xXZpiFuc$LOKT3gWBYw|>IyC3$5#1rL<?Ik
zwnI{i385rVMOAe-At{Rt@vSg&S!(GwWVa}am@yX&(2E{6v@SPA{mgm|+iv)@rlar4
zphYZ8{6yCm*V5q*K1+X6EojM$VD-B9X(>qp>wRYc2ZS@^s6!a$*;-j_vgz`1%R+=B
z<;h4iWojWXSQ&BiGS39?Si!mWx7fm%$e`l~qINesR@^&iP#&dWZ5r=#_fu<SvSf_j
z7L8n<SM*4+VhBqV(_Rj&-@veZ4-eY7(jgo8Ew#2G14CV>C+8s~r|W``rOl816`S<z
zvTQ$Ltcq_YUN5oCDM)phd9Rq_Xcb^y`b)oE#`G@#S_arwi-byyat7F@^;gZA26pz0
z)UI;ejTVGtxU;OsrVC4KR}w0R+GJA$k1*~JafFX|rI}mOApI_yPRpu!SlTU|dwUI0
z%Y}H$s4h$-<nFnln#ctAaesK7EV*C`BwKUc+PAbSRT-vUZyOlT!1S-dpt~ZhIUg_d
z`#rF{qXg)B7)>bU4T_VQ<cu2qNxsy%ZQD>*%!qBhWc}W1HW$6<ep6qvov`O0QWU42
zck*(@d;HQBuW|tiKUaYGukRHb3BDI{WSr}{*ho^rzjEXT+`G>e-Pg>q5(~es?SYO2
zT6ZqUn-9O?o=OkiB|G;;I(^QoR`da!D)0(LSwwx6RgSiG-NF5t4aB%AXbn#-bd*i?
zL#J>IM;}V5SXkZI46zvDr&KFKNy?P*kHu=5c4<$kuGyHycG)P~9!)1hU%`ct-$LJB
zIlR}@yGP468jmVOx}9K`5GN3Kevk3@bspQvR!Kw$XM9F+V+X&8-jxo3L>m`h_x!?%
zO)QU--ByYA+r`_|xuo;IBJtQTEcQ(Nn~45@Qt}TY|A#Zi{y$h_k%123L})NDZld41
z#`U)~9;;n>)3v*_N96Q7N*9?niLl}yq#ryr{8nn?&7vd;)mPVaZ`S!e(F;s5l-?1C
z0-atqIVMi_85!L7R$w&_GT)&|by-%T-f)37H`7Dmq=siY3eqZl5?I!@G9ExEjldax
zbyr6;R9oQg%t`=uL`0gt8b*a6;IoNEaZ?F);137u+w9qXVah7jgd$biPZciWBkm*-
zEB5bp5r&Xc7F&X?4OI~y73prvit}O(yToQOyic#!Wc4|J^YvP`uZi8(46@98Nyg8e
z7peU9A+<X3%^aWAk$`!ld|xO+zuh@~Gsi6}3v#m_P*}XcFn0Ha3qbTXe;0-V3CB%_
z!gAVm8}JJ&CER;q_(fgr4^Q0&fAe*~t6xehbCo%dNFsSxSk(`S6rEcll-A>ouQ-!d
zZ+vyDP5Pcg2q0gkPjKe>_^|Loc0al{KJ8Sx$1P=0151Y64s6u4J#f?)-GEJs(8_Hm
zmcr~$?)h2y6>#Z_7<rzd452QAYRIuR*=a8iSm?I?aDgQF;Z{RKSkh#SP&q_4o4R<E
zXJ~y#%j5CeP)kz*{L76&Dpm#*HzG221^mujsd*)YP@D|>T}04ru@#9jQcj7efUGQE
z$gk;af=<kY=)at?c~wwTbp!;}HD(kOl_F${#3XA0z%?Jd-a;XiExJtHs-Zlg#6sj`
zEvcL=X#yd5Ni;%PT#D)44w5T_uCnkHiy{1AI=ey$r4k6^r-a8=+BqGn9{e|7`plL|
zu@nEnFIZvS%P)WAH#O71afcN%yfvgDD1?s6MlyN$)d1Qw>2VG(t&*ehmuvyoEr1tC
zwMWB9Xy@KJWvzj?n;4nI^R9l?aRcz2^+njLGKUCozLBtNq@1{Z@>%!vR$2IUdB6qW
z{mc6d+iNV+v2$7ZEKyb{z+AwY@aVwPrzZ>G$bI#kQ3Pe7!4Q#l)1X;h&?}~-Vob(D
zi*h3-s<P0G4psU&4aa0wT=Tox41o8{LVgYSv{uuhli(M?dwN7L9r6>_rE!LZ^sMAt
zPq$=7?=@EFi7GNI`DgFXq->h$IHhtY-Bk6mwHSJR`*xbXZ9EjT+EXQagXl@<;=Hg3
zB43INM7|is%xXr*hd;b9366T03a(M;s%#Ong@h2v%PBE})3K#OhKw+%fdJa=3-yud
zgVcC@&`>kg*B2wK|Aq;z0}~J^Dli7U1tZdHrQ{l&V)K@<(S)W=LYoL7;%7K!rb!o~
zph{1m;w6phBKL?(vm-jzL(efk#BCC9WZ!a&OC#p@)QlONu__sysn5-HMyeCRtOe)S
z%x(FxG?lUhpcFUd$*yA$GSw~|OiAHS4@$-T3NDJNIq3h_leq-_dI17b$k?=X-lDP6
z&|bXM$?`(E!$7&rPy5PU9BDHXhK!8OGJXBpd_Px2)dNHV->Fh}hW7V@M^0nrKS0FF
z9he7^4f(TzWGLluf(0>it7KcA^<45lYCW{qU;)s@2hm+51d>H4(@T6`D#}cr&6Pc5
z=<LWNqNsMFyvc)Q7-hdimp$lgS5!vhlBY*e>_opKucy4B&X&>*A1I<WjnXDhh@wG{
z8cdnXvCd}4LUzxNA~TKVA<vf~MVD7|Dzdpy=TWoy&w#Gj1e>+of!BcEdT9pA|G1Z{
zK*tw2LV+f+|A<UWlfIhMyV#>CgNP<Sh;}Ifk*uBa$erhE-_{<zZ?1I?_ir~C^})Cr
za|(I=!-Z-!Mdi1R(lbNyB+;`7Mv<-Ga#G8H5)jZ<Znq-k?%|9f6?uQoLy$uo84W44
z<TIIabCOY=qA4_{<Tk&AVj5Z7%1fzYKAT=wtjg=|TSpp!jdHxf)2TW@O7MM`<3S>f
zlgvnVp0_&1JxzjigU+Ml>o&z>0^Z_axi$89z$4m6#v-YePNiX}I8hWG?;#zY`MVxp
z99cleedoqKqi-em;d=K6Ox%98^p06iM0)9Re!m&g2jJ$$(PD~tGM><#lv?p}wbJ~F
zhWi~@i2JPJ2KIP#UMTIf?wXNGm3iq%u^;ODqAx)g0$x?s9Q`uI)>VloKPX@+Ykd9G
z@<JCC)UHshvxLdHJ4&Yo&<1mM0z9?8QSjr>kF2**`b-5-(B2HscEkkxZG#_KSf@1%
z;r!(-;e;A=qHZC;z#M=qRHOjlzq?8rXs)`yN~4`I9bOuViwSrkSOn^;dsCI?8*F=&
zl9ffxy=L~`a{!qpFHMUv_I-Rd*&)NbQRvoTbQX=whj&h6)2ZkOKI*O(__T-MgaNu;
z0Pan^PM7yVG&0uuqIT!2xokAHQ)}^x2uS%KhHcL-0n@cvc1#xyUNbRNs;O+_7UD2s
zj^`$sl_?O0afh6ofkE0TYsG<DV+(oKx@)1cYmObBlyQz|^(k_X;wYsOu}>cBxMN!I
zR~sM<=8=t^w5?}n_}TMKi5WHRU4?y+jAiKe)oUyJ;Wi$i<5dp=hH_?j4I)%mK}xG^
z%q#~3puLxP19ohp?xNj(W3-_FP_&C%7(aj+U&N@Owaa4gO+Fe6SqBY4LsWS10TkCT
z)OG%m1l{w@RxyyExCfqC4oi4~7Vfm2-4eslf2a-(oro5-P!(3XG2BX}vOtFFOpC+}
zws!&z`)ceF(_xwW+Dcxo!n^#bQB~OD=r}|P0Gj2@GoJ;e<|FrHuG1tPuN2d-o#s%n
zps;Mc^JJ;Sr?DB4#SQn*jZWX@9z5IqwQ&aRnWEpPwe_S1O-)tgwsBY-azBDktKF9$
zaw`)(z5O~-t@ZQmGxL8};ZtW$;X_F=gSPeHSd4EAYr^3P)}we#dR9g%PS3ms0jj>%
zDz6eo{z=8oD0U#FMBV6ZUaL?_sL?g~J7i4<*#@t^CY$q*6MSaBR8yWlxE6(7p3Z|(
zM4(}iJ!$51q4H7-zGtj)rC7?+7Q`Yho@4|t%3_V-P)3GH!u_Hcc>l|Sq})}5jF+<h
zTkI*fc%#uRTMBQAc*=;q-gjwp5&#9>Brj!j7xifFF-1dSa@+x7igvBrwPnZ4)&E<m
z^MBi7{qb}EzQg)$o*}_h<RRZ-0u`j-;QJ$2e+Inyw?zUVDJ?!gN8X8m8?_tf5I}mH
ze`q4x5;qb|CGV2M!)8w<nJ^c#VlimIEb}$NQPP!BrNCPETCJd!IHm<0y&&Uyu}yKh
zOyor@a!*cI*WI25`mq7IW$>Bj9B;z;vjHHJ)m<d4M8E+mt!f`EnR8D(TdHQuvB%4T
zNknbiJqk}*mLU^>j;Fur{rbhu+hB;24EDC%m8dmD`;+~b4=H4fW<F_Fl;Kts%;?oy
z<xz6t!6?Cw{c`iJ?W%IDijb*0??yK8;;m&=Hj~22s}A8S(;4=S&4w$+g7hj&3}LgD
z7VV=tf8BR)w-rsh((@}6Kz`E&zaaTJAmmd4F<AHEq5vxuFiu13pwF;<<By?EbHa3T
zoTc2&vdK!Qg`|vz^j(4Bi>{l8ozf;)rtYvu=eklWd$|=D52C0YgM6}OO{=EOq&2q4
z?J@hQQoCVo`#r>S8dX-#sb}ZLH;Z}+q>;cHrtQVA@hg_|eFRBEiy0Tqx0dCz&KJ<;
zm+@5oPR+{|fOc=l7(7Lgy!GjnAxpVXWE{`xZT})zx8}q)0;a@(zzk<Bn)1BP4s#z$
z+ab8Nld1MJ<|i8|DV{Ky^DT2y+J$cye0||Go^+}5^4TmYg#}r*P$ZVP4VI$kH?ikN
zJ>C(9J>I_8hY!=m9}|d$g?z7zOZV;9hHXBq-Tzpf1iU;wnq7k$oeTvZ`V%uGpFRuw
z{~VkA_$xR`5h2<z@Zq}kK-nnp=IDG=@UoX1iWqY{AR9H>5h&gS!k8Mz9Zt(d9DMM+
zHdDwW;(=Zvo-x4U?_E{@=HJHU67QufLUsx6q(IEma0W6hs))QA=Wz+*&u!$wF@#*3
z>w?EV11QgoOzJDdv<b$c_;>O1c<8U@nX4#Ti5@L=V_-;sGl)<Rr$|8$zvd949xPrD
zkR$3mn^JKem4;t0uZaJ&NEyKA3}I7mlS93g@k!arx!Eiy&#z)UOZk||b&H&gjP(k3
zlja-SE|P#M&t>w-ZWtC)AVSp=&ECt0FiCwadccLd1g`WtS@KNlK>}C+f>KX(p{LJ!
z_dr>6EM`nX0H>45MC=HY(Agl*Xczk|C$p(=k{PUOoM-#bKCiYMvPbTF@N$}tU&K@i
zoysilylnQOnqX0%&B+D2yF}Y^64kl-2S4D^o-36qSwEIFjML7Yjx`AKgSfU)uc#((
z>j9<#=x0XSMXeW<>r_ctW72aL;oa3!eSmJUkEvrgMj1GzxRuvM=UfKfu4UBz5_qdx
z>Qaufvyt*VaS?*roC6_Q#0QYN*aM4Q&tA^Sy_P8G`T-IeT_n4RzV|VB@)TLX1j?^%
zVZ1YpIpCRu#SX*T@D5#BV#mHh^%@eU&47VOPoGkr$cfR+Pn_&-=W)|9)Xtig3{~gD
z)uyF<)AqHtSQOcGNeQ-YX?-8ZD*JZUvZwY*5952kBH<o3pRoQimWo(8^OXA7F4@?v
zf-{PXB6~(rbJvhaKP=m*ofV!8lJ+eV7?Or^R8+hrYRemvmGo@JSQ_uTV>Zl2ngJlN
z91#koConQm)5uSoFjy~C-GGJj*g;X?=4{)g1TU64GP1>P1=4@$k9htxYTxg%Jzgfv
z$l%OidfHu5rw+qTHjkZNkA!qZivpgWO}05=R1>PhA_q=z<U}^~760A9CN<R3v<DMj
zyi7&R1aoyJKPu#iJl#iH^4r~QGJHT?^wPjqAW>z(WRhIlGb@(pYmgnkgTvJb(}d3q
z7RwmoM6Fr+5D!Z$pQKfm{Wf*t0ehyArY;5X;rFz4TCYG2g|9|Q^h2yKA3qgMVzeb_
z<JqEwIp;94&OEHZx2dp2(x(n_p1`$9a5O^}LDdFsv%SMwv&nf4Av}ELI$#7K0T);`
zDh2V$M7IUh>PgLUzNSQXt4C@4Y#d85=Ox2Xg>td593oDF?0@Wl=DJ=cF2EIA{{ZYl
z#9h0*cLN;YaYnfVuar;y?b$!J3QnyDJJg@NaNLqNi%qi-LMfgvTlrS{@OlsQ1x^$%
zXjW6jSKV-sG~?;hlRd*s2-5?^q!e(EA)6=f3R=j$Kbx#Fyi*Mr;v$y8k<lZDnc`WN
zBv8msb<qBKcH6&7+jcp@zs_qkGShnGK528}<MzG5&GmA8{oo-kPj9K-qki@#6F=qF
zJLbo&HP;+2D3jZ62*=R>0*vqr)Mr#`1=(Lew#f9!`cqAeWJG@I9c0%)5CDF+ousgE
z!NR3*u4sC{nKeh?G&dzc(P@<x@e9fl0JTY+S=Nl;IKvj^kng$$X$kt|ZC7YYSFfwX
z$vEuscr34@G>g1GwaSc?Bt!Nw|9-@Z=iF}W@VM&`zS}16!lKmf$f}u6@1m|>2d??h
zc5JhUAo+rYvRbojj-2@y$1}skVfS1*%V9oY*k|+Nnzept9U6B)BGls7tt4chV%ueb
zF;*&7Y=Ql&VAGd4vt+Mi^wxMwDXgas<wLi4ofH>Mh^DVOeg+qef^sQG7h`xXobGPD
zDmvB4+E`urE16Cd3&W+0W-jg+EzVWTNVo?scG1mT6|35K#cB{}daY`Br_M?6txzCm
z6S#>Zbn!eWC!1El?3F}sadx?!f3~@F?ocR6{%EYE0{#Uw5)WAh@d@;e4ywkOr^DK}
zA6PL1v!N$-eT{4pQ&HEX@6O`-<w{23P8l2)=864;Ih2kJOCa^Ex^(dEoEhkgQ#~<S
z*cZ}`3z^zn815;@Ry`01%9qZ8dohhmHE7$tHZkQ|;sZ8-i3z)EQB{sT=<WKU;M8HM
z>gBPfPz15hpT5AKk+*F51wzXkIqE7J;`%-DWC@p{PMlapnh#T5fAGvaDY<1;r$&t7
z;|F3(R=A-HjA+*?v_eL>rE=TRh<&Ds^;cR3&zR(An;Cxn(pJ}p|8aEYWc6pn5&${k
z9-e0<ZvV?Y03z?jqON0V;gOan43Ly;1pBeGA%qKiIE9X1C@+F4r@{TJ5=!|A#pzEU
zmVAilv)pTawN9Lj8DBv#&t%#LoXu^=<S>P$lMtM}`jw&kHk40bt5va+KMqfO9v|<8
zANl;QZgM^wUEEw$bS*5p`HnIZ{|~c7bW%0c_m+_(=s;UO3V^E9Ta^p6Nps4ty#dLh
zk^*LMOf6$j-t7Qvk@d^0q#ulE^YB--pP5qbaU+!D@37IMxkRU|Yv+Fr@Ol}N8B{EU
z?QS;l5L-;Mq#lz@&Q>LBghymNr40?44=?&`<dJ;QDWlD7L7hqBX_qamO_8lY%LfWD
zHKtaUFPuTD1~7ffvbGn4%4qT6@7a^o6}m*tU*Rl7BWgw8rNO^yyG`*_hbQ+D`<w%2
zw8#ZdF2^1eL*T>W9N<m(Zr~P(K1b_fg%fI%%--xH7F+o9<tvC@n6pS8*{cC-4c-+}
zji6JtLRW@?{v%%Bhg?)04Y~}m)s_!!iWMFNmIcit)qpXMLad3F8v**Q!fa4lCT*r7
z-6g`N0No@F8*QHqgJ`3>iAYTikDDZxO>>IxTg9MZ`AC(M;Z(<}QKjKx1n-Qpl@}ot
zBlke8_NT$gvy3e>K6Z<3Gp8$AS`$IK{6si$Qt1bJNGJ`lQp5gE&Q1f3x`kYSVYph_
z%k>}RoB-?K+lxvEp8Q6xe(c_QyrY$V@z^|CVuDJk$*~s9HbXN%-nECNu(NA8MunEz
zcMVfF_uJ-yd{C{HUY^({`XZq=oj@W?Xes@U(Ef=8HH)k<w=G;s%I^&u&>bF!84Tn4
zpJA@N<2JaMi5IJZH=6{Wg@pL$Xf{7>Tta^s=>bTELQYvJTn2<hn7U<ij9^RSSeCMS
zNJl0!Cq2@|_c188lX8(|St+NJ6e1}7`6<;2cU%UDL4_;4>2A^M`9@`5qqvOiF5h#V
zaIB3(9=HH&E+Nnunp(Zw4qQz-cb<w@3ZLbxX<!EVe$=1UZ-6=i)C{BEzg!TWwB%oF
z-wA1-n-MiHHuki?-n1a2MBOxjIzYZ8$D@E=J|5DAM&e`Q<ACaN7Lc#6uSj*GU$tNL
zU48yD(t_F3FWVQ$-=8H@4Lrv`zh(3-%3DSg0yvmG9qdms^z2qyF+rO!zmUcJM&xIB
zlxAaXVLB4J9YP)IB2{C4ueBiUnB(i&^AXo<PjfTwkyyy+eYb=czKoAZve{x7Pxxa$
zQfnCo?9<0WR#gid;29a$f9G}C@4T*h%WDm4owk!(&YpT{uB<utgagW6d%%SJrw8bT
zggFH+=2A7<#f~pV(IFWu0us(3b-vPFP17eav-~;H83#GpRGA|OjO~MkAyoxgXcf3B
zVhg$5QCUo+RU<i-;z3{lVwt4e+0RsLaiS-UsQK-g5oTV2Q4MF7Q6gHib*W<8!rW>x
z%XjLRV&)nyg(WTypPWCc9RN_i_2OpnezcNf(04fY(d<=!Kgm!(dw89iQg=BX`%?@(
z6usIN9&r}|!S{48oTYMWUM?#+T^?k4AbnZwj*~YoI!Z0QaZwO9woJ)4;K4jkzd0n7
zkc!q|0@SS0)MQH4bR2M{3gd)q_@HZ5P$>6_UEX`F@s|da9mPm%;~^j-WvYg`KPcJy
z<609Eh?juD030HxZxK0P_e9dgz5@R|!-4MN%kV9{H*7uK`VofqlD&bfdExQpta)H&
z$Yx_so;WVfl<`VnLuh}$46QKj$w`zQ$Uh%)l9Vi3<%VnOT)amY@4KI1CXLy3i_2)6
zWuAST(^A^4b(eU_C<rh_R0*&Nz6Yxvh@(XELur1h)DNU;ZYeZ&dWj=FXEVxXjHIWO
zDk9?*Bzcq<pv36nB}g&Tr`nAWM3HEhuRg$E73}%G^7SpN;Uu|EFa*9)oL|P)T*Biq
zDAak`o&PlH{L$aIQva-8tBgG2vF5h!_9gQ8R`o%;)(U3NcZL1$3vB_SkzRM3-JcUt
zMs8dtJtl96Ys+7QkB`9H^KbO;^v|~+@0Ht#cRxNqc>V6rF+AV+a`u+RMg)IlF_1+B
z2C(^-!f$Ic%rBC#+oX(=`t_#|ry`5ftqyfax2$q5_Eee-^22KS<$!0U^x38Xf>r#B
zSTS@uukz1lN+bm0{d@<oLSo7E{hey4Y37`X^R!m7bZKjcs^o8t;DkK^vr=D$CBKeG
zP^B(yx0SubzHd^Z1u-nqq)J%0i28J80QlqS#q5Yovf0<$;y;%9=yzQda4xun(Ri=P
z#2Hyv_kWf95g>{tPvNiUchNB)jANWCqf-_mI0Fmw^6e^OoLj0uwReIVdPYpLx)Jt3
z-0`_Ie&4ggCWr8(MBYbvI?r6}NDyW|*G6~c8)}JZ*i=iLx>UHMT-qJlg5`h&JU|q~
z;b)vPT%)+)B)oU!cONR-9Y6KBlYTFjuZ^;Z@Sid*)snnxC^v0<O=7<8x{eJ-&BjyD
z?(48JT}#@{Hjut#u?f_Duo}F>W}T#6wFE7J&@kv~qBkov6zbpaDxAv>&!onHudq&F
zj9TECP;>nZ$X#=iZ^xi}tiNz%Z2+v>QiNe*p1df|I`dc+R!msHk4Jgw3N<f9_16d6
zPBE>j82cifD(;8I$6u(Dj5f=8QbA#_5DVvhV<|UYAe>Z!jI+}bU-@nh)1FeCniU!{
zkziUIJ?dtRD?zKE;{nTg+&kG4l4*#dxsBVUsLxX?!(azk+>0UL6at$$yaGHT&-J_)
z@lWX(LJWWZK+Eu$cXSftT5I1o{DdWP#%dHp2`p7a^bzMYLLNs&Oh$16a7kd`4!9`h
zONVSwVY(wd8ybN0r>VD$>opekV*Ry!hRvv9cCkd#5$v1}hAsLN_LpOu!+*6JEBn6V
zcQ)p)>+E%`P16bq*q~TCK(@7q@x)mYb&GBBI)D9QHlVKIlt_*N(*oIvaB_tabb1GS
zI~RP~dj{d>3Q_lkxG$f#j%dliSlI}vrB(kH_IB;>K~IP}9Knz;e}~Wi2zKI=Xv{*i
zgnYtWhw}<?i~FLX_Y6dnz-mMhHztHj#d6j&cq<ZnG<qA6f&=_h<ll<KsA1<aPmv>c
zb3RjQ)s$E$Hxg9uL`fZoH_9~4_!|snp8VdHntif_s;p;$Mpv%YuYC7lX5pH=G16%8
zB5dQ~m9E$N69aVc)x|Nhx(p+#FxcZUSB-J+!)zc?q_~Z9!3exl8p7W5>hdn4#xdyy
zmw4_GeT|iOpaggXa$-5s>dnGCqf25i_Ngin4e=W<AxgN5xN?!^#DDSwv%<AaD##tW
zq#ZUggppbi#_`=BqzgaB-#kLSOC4V}BwW*DfHqrgDdsvgyMJN}G*}b9(GV~6Y&8vv
zSr4$aUGM+sN*A>0K2)4;tkrkJ@9t?PA<gl$n(DFuqSW~Y=FtoJ%DM~^pl3NcHQ<$&
zL?4=E?LeGuQ$SClHqY(0xCaHYf$fux_6+Wk$59jON|~i7$+&um2|zqLev7PXI)eYv
zgc+AoT8OUeun<rsx}m%&BzFn-mvRE_?2#LxP4yt*|NO-TXlJ1az6xx&U?TX_G7koZ
z`{pD6{A7XPmMpk`x9R?K3<(Wbh>8#7v_d8P({lJvP#ESL2n$$jiHG}lL(xA0JGj3A
zM{j_?`zZbipeOtdV6nm?{PX$tpP*Zk-yn}S(7&qYZ|&dzIhjuC-=NMn(4Rsg7#PN%
z>AsyE&;d`Zgdx&cflSuC5K5dtYinx4e{JsGViNoZeC(}S&4&x@x5gs;<G23_^5=g8
z5dl$cPznE7n16yig?>l+V1q^YM<M<RG8Fv{ViAD?7TM54gi8SXYzcrzHk=T-lD|)t
zCE$Sz7O21vw)o(Uz+PKIV5}`C#F+Hkpfbd+3=rK;0HR19sPcAbQ3ocx9pLnU-`)<8
m#z1&`a>9SH@=c)7{$r;5#y}Z+K?ol+V5B`cw58==gZ~BLO@}D}

delta 13177
zcmZX51yEaEw>A#J-Q9w_JEb@jcXxL$(BSUDwK&DCxLa{8?hY+(#rb*r-MR1F@6Swf
zcFua%K07nd-b>T923fNVSxXKNI%ywDbApC|&_RHJK!bpQ@V005a&omZb#k&}@piB;
z(ll@?<iqwivU!8P5`@Jfr)MDGHd>KusIYG0yjHbwH_P5)HeTRSET4J)(@7hcWJ%ga
zlDu@kOS<Xne<tfFLOj~%vm=M3QygOWFtCjS&M5PwfAH&e1rE%>l_Mj9bS{&H+!RF=
zKyRm(tn<54qtft@3^_?_=>nHR6Rw6qT}j+E&K8Us0Jw15!(&~SVfV>adL-71mnBAV
z^}UFMvS#DqQ%$z6gi3lJi#YNqJmYEqoM>GyYZ_l9eZa4}Zw9nugcLjl=m)K^G@nM2
z+6`!M>=BWl#HbFxUSREl(jJ#!)@57VUDWB*lXM*?a6b4-m$L4Zc1a)WAfl%6dlRW1
z(Z5okklY!5yU#Hk6xMu(l6o<cDwL^8WdzWk=UA(p>xq@{>Zx3?O8*j2UIdBTIF*xs
zf?Q$ZH1SSNVn^N2%~LEZ$>lESnGbjrq1k9Gz1HO?*!-ZD?n~_of{#zZ2oXnL&5?;&
z?Q_-$+;-C;*xRPBHrm&Z7M`NvSbG?Fa|qed)*fb!!O=4Y*xwK;oN<-QjBxmrr96jo
z7LxVIm?P@ngf7u3PoY8A_gPpxqqRXb@L#Y2q^0A}Yp{Y|%Zm#4R_Ne2&vD5vbzj4P
z9n?i1Jj=HQJ75JskDBZ$(@vX*g_lc4<EDQ;6Tmy-PyIME$C+B6Idx!*4&JBEohHMN
ziNd1IXMZ?;Wp<a7qSM5a2leIa<$sxtI2nvRe?4|_cq>0%__7pGHS3SzS7VfJ8a}yQ
zQO40KO&N+m=f;b;8^QPFd0a895ceBQH$?KIRkUSk|A-X@bd1jn#kmu?*-uN;w5M*|
zsW6h^6<z4X=%>G$?Q@8-W_OFTZtq)aYq#S`MHDRHSlm3p>!NNNuk58W_PTuli*Ft3
zk0OZugvY4fvB6L)OQ1VZ6~bVq;bYB7xD;c!y2e=6XkF~<tiZ<xV2v9QH)(6ZTtVOw
z0$eZo;&t+X^a?dbX6&s4ThEwRMn$6Gxq?razNjaCd26W(Vp&8+4C^y#=jq(F6Sck`
zMLQ_unIk%>vxMC^X389Yp5b=&mDF;Ssqy#XiI(8it(2NcXE_+Nwg5PS)WN}$Yh7xL
zKtsEUUam4K%?~RsrmOQdO?hjiz)ybsbHRIa1t@oDr3a)^VVKi{#~yEpC2x@6ZF=m6
zV?;+N!0nU8Ip_Pv3<Mjakb^Xr9Shm8{mE=zd#`*UvFaI$Y;`5)!V}AnnWo~{BQ<(r
zPMLCMo}cuyV1`w>>E)KrF8LSJc9C{h9dh%zwRDFmhX#<?c3{DgbhzxVm;RodJu*bP
z)LdZL_)Escd(W$fId7`@XhT;=bT({(A$haW#B;l*J^uMzCLhQ|bx?UWVT(GFQCig8
z0Lmz*C~{D4o|?*5lt~fGc7hE9J}ALH<+@LuI?C#&*z0F{xlXIGWg2&W4c3aQu#pa!
zsMg^SACDsZefyghx}t(><bzHc3=8XuExf+&hBTKyO3ALh$v-)`N!gw>CQIkY9wqYX
zGmqV0@n#Z+wfvY*ph5O!eeuN~+qkus5<Em?7Ai*UTYp;9B;miIc&?%lN3O$8IFdbu
z1I-S>KyCsef$p>nbs6RZ_S{}iG`<t1p)awrFiX9u5cIi2_i-^LMP?MJmueTWjyY6z
z8?DK#JQoPEu$u8!**ZiNh~eDbhdQ{s^hMPmK~=sg?$e&mFQ6>{gOk@H<vz8xQDo{@
zd1_07-9X+W&%%wbn}bhVysZy?gQhmd+r0)>#VS{iJ>1Yv5}Jlf-;)$96NFDxSWYNm
zoo1^99U0CNyj7M?(Ga#U^lYfEa~PImYoz)(J*c)@i1<%@sS!VWu4UVCV)3N85r0Vu
zVLXK9cIVJ^4F>hp2Q1KhM+u0jv^o*s&ww#^s`J7nHVX)h3CpYj_AVJoP%NXe3ri*t
zZU>uJAyGRR=nd3j`80KC!2X);N=vqgVvnf+BgQ8QNV-n~aSYUukwloJ+;eeM+dBqM
zecY$6Z)>(2w)&ZT4|PtXcTA}uTB<$Y%2rO7wmKFr2tmvci@N<(Rh}_ihC<GSfnGZ}
zxbW`G(uCJ-ux#J;6&o*?^S3<g5T{Wr%*J!J#BWf76AK4#CMg|jjAp&<f!^r;{fyID
zHT3jfvj|962VJ>m%=*=9BFfED+<UuKzoUmr5hxji!}?Z~!@`qh6X~pU*|iVr6CM>l
zRnYjgI<s4w1j4(1OSS|-zARsfxHT}S$jVpe4wBs)2f#C@1c;o8D4OcaSzzM~bxQn-
zFC77ax9D-;=wnkAo}j=K7*wFtb*C&g^qx)NDX8{2H^3%|mafQ_NaM7>$Y$k7gPDRU
z4qHPGtK!q97d&84o=iq1Xxs1kY_dx*eUl_A`=qE&E)8JOpMa@@m{abWGflRc2d5@i
zxoa`%W-uB;T3qZEJ(YLh>tJ{)B_|!Ayb|J=U9m%6Yf6RW7~YmdGY&Fe!CvdF(t&T{
zNGiRRkR<(R!UsUCmRmtLoeR%zPyn?4G}RlP<wmD+JL>_ptE33eq1vz-+TN3z3O~)c
zxGax13h}gE)F}5^S8(!Tw;>QNvJ_^+a<Wc-xB)b79eC)WK<JA*Go~<lRrVJbp<j-*
z{S2#dT}p%|auG*``vSsvb`n6nmY07qI4Gbos5WOsCW4hZ!!B7MVt8C)sCOBQ+A!uN
z3sKa?XGx1}W>7rJyKqM~Hupw^r|4HZyqCwjRu?}rQ6Cs)%KyXk#W)5o+9zj4G-Yjd
zlDk(<$L1vpi`5ti6Oz)3o!$`AUki@sn8JD#3O0Q*kl4ly5C_Sh80H--c{uoUQ>PGY
z9-OHhcB}v5QT}1i7d6zZh8uhAbC<OtD*t_hFCES%hsWi!z8vG%q@F4Yx?sh8H#P_>
zRXTe>g445d{1cX+Kun1lPcFG1O_i_eaM4Hms=2}CP^^mlp4@qCe5GXTwJc07mU+9n
zy*SnbsdE4|%@U{)=f^?6={RC#NOS@@?$1u?{YMPvg+AO9Tzeig&)V9ZqATZbn4wj2
z5K?+$##P-G>Yt7xk``PCUK^g{<ZikwWR_l2RS1}On|3IL++hrbw)8;^48A|(dI~>$
zxhXBEmp(-h`n5;+Kle5aB_u84S1`QQ{e15xmDNFaWh8;pgFG`?QKO{asV>TpH(s=o
zv21<U9v{F1@J3ptia4H0xS4i#<&eDob;qzcMq{jNl^Qm(EDdeoUeENB`A~@#&$hl-
z&*&F9!6TK!b`-yOnX+hSgwRBr(@9=swpM*N7a@OGbgW&*zK$Y~LbR2dyg|`-m+tKi
z8UlPskBX==Zw$HgZc?C7z$PdJAkKdrylk!2xE((9cgyldadpiT2@OdID>Wy^^};r5
z%i1tS%Wmv}Y31OS^=?~IbC5`|9O~x^TrkPV7U{J#qAqs1<#RZ3k}mxsFNP&hC(Lqd
zpKCvYUI>II5~*-Tu#&E_S!QD8Zd`81!Y}#;vxHR3iH(Q%F-cMwaH~O#yHpbq#Fs$_
zYz>E$REUm)y(ytzsoA6Jf8Gi^{jv<%FEiJUpdj|^wqjquxKEzZPoWMf`55?eHH@g!
zU;q~PizLUJ)?K}x98lz+yVxEH_9uvjP4lafVxL-@ZIn)#I&+h^!bQW##JKci$PH+L
za1l)9@ZfY!1P?|pN=$-amy4Y#wS*|@+;A=~w?WU7##b$K&S2R!HsQ~%=qR}@LAe>I
zRe~vRZiWNm{rj6EL!z=g4#SP@13LprG>>N*ASE`_mwjO$^x}i^+^0ZIlRomhMh5BW
z7fk=YXzQpo4Qt9L#zMS@vz3P-pF3%V8<p*2Y<>Y3&+KlSEk=<2Z^*(T`dAMHL6|U3
zgE=mP0tAV(7qhzvW3&^+gWmCzJs!QPz9)}P2Xogm-fe1krBB~<uWy`E%XTp*q_4@-
z6%HF9=E2?9V^*W6(y#};YFLWVdrkQqVKC6wQeVCf4hq43x%cEKVT28#XwAG-bv2_S
z5oFL9;KFg;0jq<ADtLx)BwO)6O*4hn8OBIsiKb7m+w&yD#!c4tXxE8V5nA^a=^~kZ
z=kY$w45T(LXNMSy+Vk0)f+kjY3rCDVkSL%kbS)BH!gKUtpe5w_n5aQ_KN#pfpQwgz
zka`J~-mZibr)DN<Q>z0752PTRVgwPlSs_%l1a-KJ(K>*J+n_NUDY8(98&Q@$#Vj2n
zMi9r=e4t%=A!G5#yTW|MgVPY8cFsbr#*pF@>a2Y7p@W!NZzP&t63W8FY@4_wc#sz>
zdmSdsHfPr$>w?VIA?vvR{71Z$xjh|bw_q&afHTTUg3g${NP@2CCThNl7HJz%MG!UV
zJV6Muw`syWFNQXy*=U0{wL1&9f1*%9spFQQA*$91W}f)NVS>0R4?iIH-~m6du!gpL
zC-F6yZo+{p?fkaB9H_k7os`^L_M{Iy0KPa5wGr4T^8eBlD@^407H~eyTI1uQYE+iI
z_mwE&YFiV5MHg+1#q5-|@mmf`;u{%CqT`jJ5`8aY^4n*L{MEH^$5o*+jnrff?Enj-
z2tKyyl-W)je#>DD3#FadE#1oxyp(Nr%Pz16lKx&Hn^@R}CL0KH{kd)$n*n%gpUE$)
ze(jLdO*7|uD=GkYc)ATW-=O~SB*7Q>FmS_A5D=dc{_c3L|Acg_aqp(kzvd^%xym0U
zY>kHqvVv>5>`D3CuiIsIT5Y8CmB}iNw@x#1>SaajHc0?}y4yjFzfdHodqTDGKn^Jb
ze#uQg#(UD7`pP1p+ds!k3g1TJI--Q*5OTvM_PoE@={Z(iNnm@8+<4d{0Lv_jog$M)
zsBi_82s^c{$fO~KvQ|zr6S~92Ag&|VRxL<PfWu#nKU9tq9%I&fuTZZs<ffKLm#Yf>
z=s}72?n-JCwZ)W|YRs5$Fb0@XsG4}v&`S{ziL$tTJE9D2Yc|Z|Xw<kk?rUyVQWF6s
z(qk#eGMnu}aDL<JW_pWT!x$wjNxN+>xO@bCy!>8^<o)$m*QX?2>y|e7`VGdBFqDiO
z7H-2xx}NU!$D^M424M5<-O2IMeXxWmW^G>+Ly<S}RwQpO8av=+R&*0~^Qe&hk;dXK
zRHa?rm?@@G7*ADXxkp;&H&o#Pz10F;M4!b(sgNS{GaF3!8<SfPI%WcG|D|`HgCQa4
z5&Quo{+NMVE#nnEQ5H(bRVr8<#Mef$^)fc)dL0;TVQOsGJ2%U;Yt6fx8R#pXGDpJ2
zJyDjQU$X}nh$V!sq1;~6!HyslkL+)0wz;~636;D;N5_p?x-dJsWKgiVuHe`v5-9;9
zslKxQ4dzZBvb}Fa>Co!uRx3TyA#w+GONlDIOv4G4^&VuL`4C;y!f|HswgJXnIA^x5
z7qEoudf(V*AetFrem9^#o%smYYOMUO>s`FjH{=&1+w(r$zj>34L+PSDEWm;AQ-*e|
zV*XsOT^qil^nSDCAgS3BLZmE^UTA~uI_lvgTe8Qe+EZqS=F0&%BP_qkDGLRxbjZho
zo%$cJ_aE7lzZA6C`f6Y?93_E@Z0TKNRY}uwe&JGSeCJpUtV;|%HELRBE@J0c-K4o7
zn#yNvMZITCojq}~w{v*q4gvj*lwc2aLQ&<x!o<`d&Hp8}*cYW97+1AnTkuKw|8V+~
zDYLPyV6*1SpE!zySQExJl|GtmeoiZ)T@z0)IjqB;q8-X;^fTXK69OFesV0oayW2Wm
zskm9{Xjw`Q+<Z>8Bgfa*^%J_ih86Y_660a@5KrAhwM7m0AF2h9lvV^n%OHddOq6EV
z-{txon0Qh)JCC`%?rMFU>KqCdc&K&D&ZC+^XKJa(NgwM&tRsw0y{j3T)Yc+dQ}+;b
zE4&RuS{U+v`tey0q`5MLViK*FQhJO{XEFfK_N1)j+7T2Z{WcDT$tT6vGJ*PeNUg)N
z$2T@v0NDjV<uCv~L;wgd0+D?Q#)AZxl%7#w>2h*&E=R(_SDV8-%Ketx?`DqXt89H{
zT;7T`@6Q}BlK=EaC~n&6{OY)uAS-zB8p_T4@*()3Azk8>KqJmP3t@ftxF5asT1RDh
zlV7kLZ*h{9VwN>Iy|}L*oxj9o;e!y-7FNq`fJ$A6idr9P;9MZH4k+p45v0WV3zMHv
z0NM%cuCuzgtDV;M@@<>zllILT@OZ#p;x5uKY7ajCmgU3lD<5+3wl|E6KIs?SA80X7
zUW9$YplI97p+g)sY(nI|-%t4%(vncm&hDy?o``f3u2JZF=i07MhM<tnksjx{&RNP%
zzCpa}S+*AZz6aBYOVGUYMDmwC+Eh?dKy+Wbtkp5elY#7-s%&B}&WyPTc;`qgGRyvO
zz}~|9(Uz7R7S7kHCloMHZVZd&24ZtVt0lgH(rFY6(i^S^0K%FRG?2{Qk?2PUqAc_!
z^p&;bLl%x%W&)^SeyhtPS2s%H;xZ>h($oO}!t$oD>^cLfpr(>mQ0SU<u-r-YM1a%x
z(9}aoq=Dwf?_!md+8PNw4vN={jlGoai<Uj9JQZi}>eI6eP<8SRL5((cncXEQx>fBw
z&`pgx>$TQ_+7|8h|9q5_JFH|L)YuTc#`8dT&GR%0Sj<4xje0JJATvEwAkLOl*=$T}
zqGr4P9+w?S^`~|n<J_mVR~wJUa5d*wHm2-3qj2Sv9M{%HLNUv4TkRbq5AOW=TY}0y
zWO^P&AuGfIAZ%pVOU%!Gts^oxNaylfYwE%1J#fMW6k%NC!9hSclKeFrJpY)DN}W~5
zT^@9Q)@cy>O}2K)jvlgL^!&tVKa^zix(Fi&2@+!l1I?l4&8A>7lTM}y&s9NhwGA05
zniXpJ&vbqN&m*e4kz<q?;onE3%;vVjwe7+EyDQJv@mrrw2qg(fnp7B2{T#Zk`kvkY
z{AR7NqUgr+os@@WB8~!mZ)F|xv(b_4IMh*(`8bI--GIdrEq6~xRz;a6w4!M=syap*
zZ>5P#1K-e8n=D};vCb0N3fwUpLbGuAc*Vug`(~QACuw$SAvnCDDK&*|x2yC+vgXX@
zv<`5|DsM`USjQ+`H(a^GO}Jm>u^P=89S)NzdG*BJ4cE^P^CtmRB>SKnYDcKk(E$=d
zAP_EMqhhtk92|-c0r5O{AM+8<B!;yU4_48X?~!xbo65c>m;DsTst0|N<Z*CVwld=t
zwXt6yRC68KXJ4+K01cIP+gWFPCNPRG{Z)MZr^cswg*L&LfkEPGedT#P_0Qc><C&(0
zElB>}ZN|r`1CR*vkEP<)ieaBqW;_&+IFIs1+l}ir?^PqMAm*X8&TSc+b9ZpGti6A&
zh6`+F3z(0a_TP-bRku6kwXuI-5hy8Uob66>Xf?`;o-q~38OOc;w(S1Ok}HubbS8g+
zIL#29x~dhA$IirbVQa25pHL^$E7C!zQ{agp(|uHIn`Z}FwHK(zQ{%~j-1Od(B%mHG
zu^ISLfW15d?OCs(7f4@ODdt<xkrTIM?iY!Y-Z~1`2qbDQ<R{wPC?kh|WanQFXo?Io
zu!jk&AP#I|cXg}7NIxD$Xtbg_lqJHSdNzk$0jFe*Seq&xB6;?n^rRX+K%I^tb`Tx)
zo*uVP^SFSPf`w~;nO7vOMUNgi(hN30wxgM=B})II?@u+>$*W8#j(XBbtf&zq#I0@0
zJv+%txrJe*6;UBZv6wAwN#<yH!d8+HhkhVLB&3^GE_n&$Asnq#0&pk$<hj-WKF_E!
zDYvB%rS3bo;`li9A=>We=1gdE?pV_`%Q(tq#3zE*P$V7XYztoCl+|K}A!(i1&}T86
z@Ks|}PdJN9@$EnK<EZU*nPTmANzl9PWmn{H!8p)WnWB7Et6U8qTCpDbkU?roM>lsa
zcy==oyhgy1OZjIot8oAQ3D&JGgovM-!)BuvAp|bje1ZcL;$&uHRX5~$y?-KWDI}`2
z<v)Nt=q%zq(X=Fbj~A462MNVWwVTmOvn-h!i$ng=juY&3SCgFI3Xh~LbeGUelb}aR
zmaiUqmt?AH_}rc?_(=l|=N9P<aH2*h;1G@@sMqQCJKcj_O<PQ(VrEfjqTV$++Eqr0
zOZ60w>E7Z0dB6|}mu8s@d?`pzLdD|qACvMgVrZXuC8c(@S3*^BBX(2~F>3M~?Y=|0
z?|3}bdM_!L5fIBu(w|$9V*B<Fs+m)h5-#T_!?HYcVc!H?63P|NBkV_Q{yx>h7N1*W
z(MJw(_pW!Q-xq$sCqC3y=40f+L>JX#jP{>|64mD>P69EPE+kM8edz{~eOvlH09f{H
z^39(l%T6MkY%egpvnyh0pVDh6oP=P=Ac8jt@`kwtV4|B#B-BRVc+z;P#2VGLPvM1=
zvxo}?-NOJ7>M3NFBp=1~)jc`&Q3YXK06Q21KpC+F6!AWmS+t#ktF-ODS&9gG^muX+
zai}1B1WeDwFX(t0MZZ8j-<qL_gIZuN5q7Ys$%n(pP4JH4T~sWFoiPs32WTRs(eub5
z#d%-=?J}Yfj&*QoD`=B)0vhr!dO0TlpLpiR1-%8%1z$@4|13H92$z)%K9CVI{J&J*
z=s;z9;(qur#KDE3is)(Nq~d--FfhyL)F+f4=(j#T{>d$&(1T7YnPL@ZQ8C)Y^Zd<K
zs5ia_QV*H^I{um1G#?xArv8LKYk%36#1Jm$d?vgp#;%87fWOV}tR$Y<U|14u(NN~K
z_`zlBb?|$>S6}D!FlGMF`P=Q{_q6$B<ZDpcycWeT+`W_Y`R#OJXP*@LOZ0hO4XxF}
zN_S0TmZGHRAME+$+mjUVt=%85*aJu89;(~G0TZ$dX}G%>9jubojtu2kAeBA`iEA4)
zyjbE%gdbGiN9Q4QS!l{`^zBAIKI5@Mqpi=G9le!(x>y{6S@EO4zkl;fRe>U?g@0yL
z6qEk7AK6x=vRG^Zg(pxtWz4Z3uk!r4jh$;x)c`p1<lu2^@}0NX>p%H(#JS}h^S1Y>
z#~=Q$`QrVX=ScbMAcmbY-jXnF!b1c20oteH6?X-y^QYS@{}c6?=-sXh4gq9`n?6Ki
z#&y3rg=Jt|PINKMq%ty_=WO*{v+r=fSFjcVuXndDTtD73=vIx%+)5Wni}$qtZ!8<K
zx0OhfC{W?Lz!tU=9n1oqUdD1^(DDUg?w+c~&z@pktaZ9i4W23*Nz;EyiwX*cCQZJD
zc@GQW#fZu@ee-t6j@+T)M(<|%SQzh!Lu9I@a@%1erJ#wTUHOrIT;489WiXX&@w>Xa
z#)pk3;|~NYjQ61zWMu58FN~@*@71e=A!w@|)2>e-aomvl$zy>xah~1~I#3r*Qvk5y
z!hr%-<q#)#e8`C_%hU4C<J-QJIpFn_tO<h{y8T-E_2=1SlPF&-CnQdj4#nU#%qHAe
zJ3Xz<+TD~_uh-)VF+tns9U7}JXsVv0wB&iL!Y21#`@lqp#gybJ&u!_1d$H+LDEY7B
z5zO6(78%B0$zK{Q1jvsgY(RbjqO)YHaBr2$LvC&K6P-+;h;eK{uJ1Kn=6ePMfsE`+
zVmq3CB@(ns?|PF@>hc)$B61r+!ByYz!SMSEYo$|a%>e5%e@<NvJNhri{RY4LU@u@%
zHO)wAh9d+8r?~}2HL*AgwZPPFy+ubx7x_b@{ha2~bPfNiW>V}}2FiIdw$U4m`=LRb
z5H!D`>k=dm$8PI=x#}(wJ%Vcp6qGqLOdJ;kKw3+~Wo?@?eqw9fjO;i{%rSSR&3Fz*
z=II~N^phox_caId1EqU^^9Nh<El2t8l}URhAzkfxaJU~k*W>&dei0BHfqkHEEo&g@
z?<KJu6n5>N{N(y6SLCrXd)>WV#OWv#RA9?&8!MOF+|S=H(|UX>%*R__M`d^DOM?=#
zepwzcMI1i#N0`oqrI_Ea#~OUrpYaO&PaKcf`L(-K9v%XMgcAHlPYBYLcOv3N@4`C*
zk>BNN*VFvqpoC9l8kBtwswx^z)KzJ3Y#oD@U7o3uiGxMNR*W`8yjhT#;SqznNqQx0
zreDDxHabesmz9r>sv43s@(l1M+eAK`Sk8Ca_?*n<t`JfrR0Y5PyBH#oYfm^+>UQf(
zslNk_nCf<8IDwL^Fg^&AU`^NetuxD8W00Bx;Yh7D>bH2MdQ)31xrBGjzDI5%_Qq|i
zOl^=J2J;P!2xJ#PPn#(v6%SKk4x*S3vZ?(=1O@QRI<<Et!Y8`dp}gn|I)}?w9jx#^
zV_#{D*aV4Xb0`usO*G>iQqkiCX$CQ5d05t8;v>&6NOrt<7q*~J13!^hpbS=@_uTo#
z&#h-F=?l)hXI$G=+w0g5oFEhwcXg{Y8#G9-IU$I{zJ5z79&VU5r010qLL&F&pJ>-9
zPi0`~>H*ov9q1O$lr6aGFsVV;r=`$8%h|NncO7k^4DVE^-OHJEF6(<<z?ol;somej
zuDA{k5-E>cf93*(E}%LzyPcWI<TWU=`oV57788qzi*v;k%CpWKe#^wM1BmPupg^FQ
zYtLYRGm(1a3!!t|DPp8URSr{Q*%r5{`^m!G`}SxwsIPHz<(;wxmE(^+`b}CFlUhQY
zQv>Y`YVj+C&-K+RFKb@m`#VSEi_7DE_#9B_6B)?KhaVJ@#TOeDQq%pq_lH&VW4G|z
zuiLq)XK$T>x94?f(U<#gWC#!3GjH>UQ(t^8Pv-H+-X15^(EJRx95FWX0g2Y(*<by0
z1009gs`rW*4+@V4(gm>j;5JAXm~jMqawsfC^n+W5e96By-Ul?phsRaI*WuIs(s?k<
z=Ie2o)CmS5)j*yde@DbUL!KS(Zo$WPdBP`qT5#O?>Ed9QAg&mikkjPPNgJ2$i(HB*
zO-+M%?DZUhj<zfzN7C*%r|dlPgeY887Wi$Enq<WT%H~6(EY`)?G?2VyrEcO5`sdbo
z=`9?;XJP_+x=+YC`a$|74=G*ty|$-wUliPCAHy98;KPP(g0Yx{C`5uWq*5VWbTe~;
zINiK;<BMOOv3R;_pWY&qf0P)1h<nfxv9={Wfe*~kvdkQkJR$Lcp@^Cx)e^bq@#6fr
z8BFFZFt?6I=d$Mig`!1?O1elk;D{G``*Po*Z{3f$UWB%bp@9Xnw6<CT<&UJa89IM{
zrz8TXAJ`&|&<6TJB_3=QHv#;x8Jk~jk7;DQKImBa>nU00DAt{mGp9sPU4X(@j8a4*
z_jEg>k8{8#Xpm`AaLWVeZsK&EI%|MSiW4Y^{aGXeuSUg_=7J8!%MN$36JM}BwJYt_
zn*{L(2ycXac&y*&8q<1sjZ1ln)jkW=auNyz4mdWBYVzC1l_CD4By-bcI}B&o)|J(W
z#Vf|U$Hfm&++4V1LqrvtyHu|cAYHA`IrX;GNS)yHpJZC=u5DlL33iX~3F9zmtIo`-
zIH0|s;^l2{eJC$(wY90%CES!=GMb%LXv4R5ro|%VJeReM^L<RNgXT6jm(Prh4?zQE
z84{6Zsw8Jf$Nd25AgAe)Y4`L>LJ<*Cw>$=v-;)m!Acv1ydJ~Lfdr+$k(TvB8n{i6f
zSb1ml(uF1_WVmSF-i)m%nSP(}y}ML(8RfS$D~cMY_Gp|=vUHvqOm*y7HML-zS7>N0
zPl>fIF>_WYs;E>Hs(;ny_U3D%^<xK7-A}?Xvi8^NjG?-5S=ackoGnOGH6InT&P#HH
zIrXhc>z><?rNn}z>Ntuc7W(0WbG37K3;S%jC+<iz+athHV;YuByjL6T-Zoq%Io`5|
zry1_%HthU5QJ;Q^SG8%QC`UAWo(_j()~6-#Hu!Z;GET|2kC0Mp0&?WYa?SyHHyogi
zag483t(RhZ{t}}tTIu0v$~Q6Gr|41T{KRv`s@YqaRvb`DXrK^PjDrx?t%z7#_*-}O
z9IT9!W8iIia157_rYfTRk$9jJ>!$^DU3~A0fVy}dX$qSxFOK+TPkusElg~cxq0n)Z
zG&`JCo>}T3A8WIKq;!V;b_qEs)`HDTQ<sYL=pbdCj<=hEcAr%?#T5T-N`2-D5Yv?Q
zc?FsxTNW|Kz2qk}i~4vzq6It3Dio@@!fAkYV6db&ia0E=@SvIe$7AyZ|5$$!w*uX~
zNNmE_TC#82==JaxBOf>U$44JFEOhF#6HnP_cSgsf)@go(c(|kJFZlX`0!$G+)D*8c
zjGe7^r#h!N%ZQF?5A#+YmdtaTh0{;UQ(uERWr~$=ko3eS+W%za+Ulcj)Ynb5&aXsq
z7hjs*_USpm!DBnKd?Tm*!0-W|D{}#_V{Q9MSq~#SbM+2UT-~T47Je3XHuVP!*QdLD
z6xxYVPeUIs-(v^t?5cwtl{0UH;aSnS<4Kzn7q<hJ!@Zr=_Nj-M9KEGFkGi?rbi$-N
z`>0=c);APp#AdG|h4(=04qt@0je!1NkZQA}5~4x7boZmRKh_hT7h1#IPtpFV7B>{K
z@gN}je02pHz59+3@ZdNRV(=RuJV?cEl@0w1*@-B^Ed}c*v#jq18iL9Kk4-2|y`<5m
z^mV#AR+Z7eZx+wZ`DF#0W}17`DH2`(9J_CN<;k*hrOy2DnMw0o<=?K*<(-`(QuEj{
ziM>W7oZ(IF_UEM7R+7=KQ7o3!kFyHX%AS!ZW;U7?hD9(M)6{DXX0A6pf&|G|r}VOu
zTuE`)`zM;n6)^b>;h9AaMiY(|#M8Je+FGiQoQEan3XacZHz0mh+YmZVI$ZY2q}}DU
zQZFpdf(~0Ttr}?CQ#3W}iq<Ad6IYs;QC#W`w)|t9{UfKiJXU;ta1#GJCSsi;Fl7c+
zy8J@f6&X)+E4!NQM&~Nd(-$N%g+RO{WVO@4C3?t2D&!;1%ENBgOT@JmqVH3y_fsoS
zovJ<7#jRqIZjwdmx&gPQ1pfnL{T#YC7ouHa@wttWUBkEdrP5XxIW$V5irb;AlsY)C
z0v<UVUNbv02Uy&;q`b=G8J#WJS!;VS4Ea_`K5M*bF0`OS>1V!Ka$=B8@fQ1md{=$_
z=dY%WMh=3C5;+3`GSA>_)NBz{$Po=j{H3B#yC&5N%GuPJ^3jwx+(?$DmC7ZX(Ea?Q
zS+Q8+jm!1g%ANx_(S`~ZyXWF-OKiJndy02H9!m4wvEUaNzK2kp;CqD5SI63}nT1Cg
zo{;V&DiB7Orv{3FeT-f)puI(*FA2ybzqGH@w=q+#jl9F3FlvC6r<dgA!&tRQWlR=H
z4Sxwsnr1T`G9ltEVQHfuqx6)3SWDG@KdVUlQYFs%S|EP?23MT%(ho->{H7bsa^&Rn
z0R7jwvy4%JD1N8h7mUA@s}+O?8N5>tJM0|(3GGAYaO@o4r`f7<5;~I3Rm4@}G|9ST
z?W#&Mra_O#-(KeYih@mbRl#Br6m<pHzmEKbQYmcoY=g+mbigwn4(Ss}I<A6D8iv84
z#Y6|n^(8-VPO}UF-D-WV-&{xf9GXjBg1PG_ic8$qfUuTGhO!g9wVp*FlX?$=ln<?w
zC#Q1B9iL4of&_DCLgEDQ3h_h^Y_Rm(6lo$c2|69jIvfo!fysLi-)n!V>^xrT5KFJ~
z(i|1`<#fXJQ{s3@N{z|D`;s&Ytm!VphYPk6p3$^ITU4%E-{>*`f1?xhDtA!Tm$xd<
zwy<clOqE>|9T>-q=frb@3c2H&9)#q#=(7$pmZ)<T=>bSehDCTm`0B<fYR33R+S7JG
zRzFpm8m&>}@A?8ff})C+DeXsx#On5bzksC4&agy3&Rz6nzB72P{K)6{&}6^nW0$Of
zRGhL}`vE{P5x|tPz;9i~`MG{TKF&W^oA6i;J)Oo#(7QZG%S==mlzVbxAn(6js|iIK
zBwUCzaqa8gRQu4Dy7aS4TZ!g1KN~J)jk}jZ1Ib#}C}*}$na;?CX(qw$ho3$)8WzPj
z<q6HB#e9P*<ZDTy@i`1#K?4C>?Q@=I$K{ke;nx8Bq^v=b&u=RpcUUdCcu7qsR9shC
zu7fE*SHjF+Zt$qDK)|^)TGp?qfr&03C#f){sUwG}DR;dvGIp3rzMr(Bb{d!pENZVm
za4%|l2AwJl*IDFkn~?fOli34n5>jSE9=@C?8aBzZ2)PLtcJH~>ON)T1@|GI<_ZuW$
zQ~Y)+jGYb-_l6I}uQS}?E&1f%Zl4Eu4>#_!Ngsr)MSyRClvCdEev@&2q+EG(i`V?Z
zn^bvL(jsV?Q_Gtd{C%PZ2b+J(FNimInje$inY)Yn&Rj~6mVwjCdy`sT<C}3ST89__
zQtmgcj7Djf1Gq`nw^G4ALhlD}lp^6fYe~XClsc@Wg@i=FxSR=hbserrP$!j4!-x*m
z3l#qzWy{*H-||#0tAs+5=N<lP?2&)G%1eUzB{XDT)gmc$)luo1$k!Y!t9a=B32OYH
zBRSK__*4coVZN)%5bEl%a+pHr{lFjr<0VvVCSWg*Y&4$^^q@QDDW;Yyy6K9%X!CNx
zR!f3Z?<=>;nDutYkbJqz!mL!;86gdJNK9%MDwZ?~_rs|bSR-8K4hOa(lG6A<{6pJH
zT{>XJG(oSXA4Pq6%R6mT8H7rHK0tG_Nm)+*0cFf0YUzpi!Wb@{aMtZw^Au9s3;9a2
za?LRIGL-hPb=7`d{AB2>3R|N&iro0l+D?J^(&F|rQ&c$~fj@yj8vRMB8JcRhZ<-IN
zPeRI4!(p*2w%_Y*^BY}>QW+_VjuVS#1)xzh<fc-+M-mSsXJxbk#8UlEw<HDQhkVF|
zKpWeuiJ!NX86lM@@@7L1mZ@;{mh25=iEV)k(F2<!Ha)}Z=)`TdkSNDL>6ws2uA&M*
zLRXVq``>nCAjpQ4H;Z+nUE5lu<5lr0g{3>VjMQ2hEMKhE>YZ9Wqjt)Qu#nb1rg_$r
z^>}=tknP`)VD4r}+mx$vAphN;H<Kue0NPAv3TNPwDiFyqxP4(h17nEj-;Ka%jntYv
zi7|Nv`<*FoJdiAEHU&6PC)h6my)&gncJWrXm@|vnC^oA6NBccb`r@lOkCjk9yNCr3
z0fTqO<CeZw9lp(cor5NH_v1KnA>wj0%Gh>zY<8L+fkYADC%9=r5qL#|BO#mqiarMx
zk6--?7G*y|e_NjH{trh*R8-?9-Z|<_@XzTA;$rc1us_LI|Hr2HY{Ea0g*+yt71gc#
zm3w}5zi}Y44;IG#p=$cQ-u$2Q^?z)7T@_?-oh;@4nR{D&9*#sG0~W%DDww0M{R3uB
z5mWRn-L7nY<4U?p0f(w3R5xf*yIfSKd?E}O)3V)I@)*rsDVN^ExCB)HSPF8P{Mex-
z7%cZbhly0{w|UJe8{EShu81MpxC@wgJG@;j5lxe0{Y?K!^iqtj8ljdm#VuHT8aW*G
zYb)ICo;;t8+SdUDzrEUq?}4?-q*+{j%gTKI+xDSE3=1&uf*Bu#Ic~nT{Iq24u^U2R
z)Hg@?P!~7d6hUt_`}14*4+IcVyGcQjd1Rf7R7=1zMIIi2xWOkQGWUK1o2vO@N8>X~
zgm^U<uf8m1&38<}{%QwR9ew*&Vn7qSd0Pi^No%X4xzzqZcViInlcJoT6XbRBrNifX
ze^M7&9m(evy%BCBx88}|qHI-Kc<cU9nPg?j5LO;%_RRx8u?ejjc0DM@ONzwKsw$t%
zm^Az3cj`)?cvabF*GGgJ$!wvZhVP#|VS`I}1_e^@OI_t?mcIeu%5*Cbx#ve32~`{!
zQE3PV`_PFNd-_6@Wqf?t8h?vTH76|<plEL6H-0zZW0Bom@UI!Jsug~PurHYU4S(@k
zb#?-@AT(6gPK+-%-)sfRI%S}*(AW(c1r8fUIGHDv_rZpBP-L=(dm@Alho^}W6#!zk
zTPRtXdpA+2-QjNy)(1Z{tZb3H=xFXmf7bQUGy}qKGfyYU0-X*tEhau4x3*X3-!zAO
zdT*vZ0)*yzWK0u2@?w6Tr_9$9&L;tKm7ytm<-VR1Afho<y^Ja=ef}W$>1dt+Ar#oL
zfRJ++;+Z`?LV$3qOYOdZwc_%j)QR?|9*Z8e`^G=ttgOi0XDIWdbEo651Vo|puW<ne
ze^my%(FMY~bu>l)damG)#>9Zkxg}LqNC*hqhG%1PNKi=9x_mzy_TVP$TX@nR9>uX|
z1P6EqUsrNn6+fgA1H1)lnn=%=^NF^qNwkYg*U6XtZ271M&F?et+Q3*A3|Up|zym3n
z_2i>g!vQTi*Ib)k-NO8j(0p7qlT(NrDvnG-+H--onEKc<D=Hf))~Ei`?4EuoQLWT`
z?&w-W44{L4YoaywE)$e3CPmg-s419h<nsU*UW^$?`umWR4F^Fem><(YP^TocfL<d*
zG|g5~ymw}OD<W%~K1TnO%zST{y`~2Qnr;4(&5LLM@FmD7?T?b#y!Nm4D!#9_UU{kf
z-N>c-vgo`Lx7OG#9)Tdn&@D5bQD!zYO2nEO|4&VYS4o>2^GNCz^sP~IV|+l5)ES||
zhUcmt(Yp#j8S2!VG&MD?PR`^}%=F1)h&N13s5inQzM%Rz5e<JK;~4hMkmHahxFfIY
zCGHp2-)@>m2h@fXGd*a;fB!TAeza7CtOm<iF%$i7!!`f{0{=aL`1eH&POxHydT&JS
zwZeu9MF%HZ;eb!92=M>Yi}&B?Cak{z4r>DZ|FrJ>H(&toFTnd9AphRY*Nyn!fIh;%
zfG%qs;=cc_<o-8Eg7hy4#s-J@-+cLhgYaqog5=&o|E3@zAh7=ZZSNxlI<T*eC{ze5
zxWk4YN|+P;$A*^p-#7{ag6RLi+jznJw#>x;dcFS*P!#~D3zC5IY|)ATbuj)L#3l?L
zd!Oas)BV2zPom(@;v`^pI|e9732^>L8nA^GI=EaC9r@qY;D3A<5P?VRxS&`i{|*W5
zxuB$_{toTlhep!ie0w1%QaSME`-MXlOyVE}6{-W)d%u7T!5I#e#Q(eP-faxVzf-@D
XOpL&b4#H5sOu)E~lyHmY{|Nsd&r5}?

diff --git a/tests/text_tox_timed_initialization.py b/tests/text_tox_timed_initialization.py
index 66782372..80e0ff48 100644
--- a/tests/text_tox_timed_initialization.py
+++ b/tests/text_tox_timed_initialization.py
@@ -18,6 +18,7 @@
 testdir = at.parent_dir()
 tmpdir = testdir / "temp"
 
+
 def test_timed_initialization():
 
     F = at.ProjectFramework(testdir / "timed_initialization_test_framework.xlsx")
@@ -32,17 +33,16 @@ def set_initialization_basic(F, D, year, y_factor=True):
         ps2 = at.ParameterSet(F, D, "basic")
 
         # Set initial compartment sizes for initialization quantities
-        for qty in itertools.chain(F.characs.index[F.characs['setup weight'] > 0], F.comps.index[F.comps['setup weight'] > 0]):
+        for qty in itertools.chain(F.characs.index[F.characs["setup weight"] > 0], F.comps.index[F.comps["setup weight"] > 0]):
             for pop in D.pops:
                 if y_factor:
                     ps2.pars[qty].meta_y_factor = 1
-                    ps2.pars[qty].y_factor[pop] = res.get_variable(qty, pop)[0].vals[-1]/ps2.pars[qty].ts[pop].interpolate(year)
+                    ps2.pars[qty].y_factor[pop] = res.get_variable(qty, pop)[0].vals[-1] / ps2.pars[qty].ts[pop].interpolate(year)
                 else:
                     ps2.pars[qty].ts[pop].insert(year, res.get_variable(qty, pop)[0].vals[-1])
 
         return ps2
 
-
     P = at.Project(framework=F, databook=D, do_run=False)
     P.settings.sim_dt = 1 / 12
     P.settings.sim_start = 2018
@@ -52,24 +52,25 @@ def set_initialization_basic(F, D, year, y_factor=True):
 
     # Basic steady state initialization
     ps = set_initialization_basic(F, D, 2018)
-    res2 = P.run_sim(ps, result_name='basic')
+    res2 = P.run_sim(ps, result_name="basic")
 
     # Advanced steady state initialization
-    ps2 = at.ParameterSet(F,D)
+    ps2 = at.ParameterSet(F, D)
     ps2.set_initialization(res1, 2021)
-    res3 = P.run_sim(ps2, result_name='advanced')
+    res3 = P.run_sim(ps2, result_name="advanced")
 
-    ps2.pars['alive'].y_factor[0] = 1.1
+    ps2.pars["alive"].y_factor[0] = 1.1
 
-    ps2.save_calibration(tmpdir/'test_cal.xlsx')
-    res3 = P.run_sim(ps2, result_name='advanced')
+    ps2.save_calibration(tmpdir / "test_cal.xlsx")
+    res3 = P.run_sim(ps2, result_name="advanced")
 
-    ps4 = P.make_parset('saved')
-    ps4.load_calibration(tmpdir/'test_cal.xlsx')
-    res4 = P.run_sim(ps4, result_name='loaded')
+    ps4 = P.make_parset("saved")
+    ps4.load_calibration(tmpdir / "test_cal.xlsx")
+    res4 = P.run_sim(ps4, result_name="loaded")
 
     d = at.PlotData([res1, res2, res3, res4], ["sus", "inf", "rec"])
     at.plot_series(d)
 
-if __name__ == '__main__':
-    test_timed_initialization()
\ No newline at end of file
+
+if __name__ == "__main__":
+    test_timed_initialization()

From 3e76bf789ea2463d106ecb03dd1ccd12e0bf252c Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Wed, 26 Jun 2024 16:50:33 +1000
Subject: [PATCH 26/33] Update to Python 3.12

---
 atomica/programs.py | 2 +-
 atomica/utils.py    | 2 +-
 azure-pipelines.yml | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/atomica/programs.py b/atomica/programs.py
index 8e3edad8..277eb3f5 100644
--- a/atomica/programs.py
+++ b/atomica/programs.py
@@ -1423,7 +1423,7 @@ def sample(self, rng_sampler=None) -> None:
             return
 
         if rng_sampler is None:
-            rng_sampler = np.random.default_rng()
+            rng_sampler = np.random
 
         for k, v in self.progs.items():
             self.progs[k] = v + self.sigma * rng_sampler.standard_normal(1)[0]
diff --git a/atomica/utils.py b/atomica/utils.py
index 87b711b1..5b532388 100644
--- a/atomica/utils.py
+++ b/atomica/utils.py
@@ -606,7 +606,7 @@ def sample(self, constant: bool = True, rng_sampler=None):
             raise Exception("Sampling has already been performed - can only sample once")
 
         if rng_sampler is None:
-            rng_sampler = np.random.default_rng()
+            rng_sampler = np.random  #
 
         new = self.copy()
         if self.sigma is not None:
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 4f7ed50a..92a02243 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -9,8 +9,8 @@ schedules:
 jobs:
 - job: 'tox'
   variables:
-    python.version: '3.11'
-    TOXENV: 'py311'
+    python.version: '3.12'
+    TOXENV: 'py312'
   strategy:
     matrix:
       linux:

From c2fe4ae8aa394903e5c7c0cf1645d587a62c4c56 Mon Sep 17 00:00:00 2001
From: Romesh Abeysuriya <romesh.abey@gmail.com>
Date: Mon, 19 Aug 2024 11:05:03 +1000
Subject: [PATCH 27/33] Rename stochastic to multinomial

---
 atomica/model.py              | 62 +++++++++++++++++------------------
 atomica/project.py            | 25 ++++++++------
 tests/test_tox_markovchain.py |  4 +--
 3 files changed, 48 insertions(+), 43 deletions(-)

diff --git a/atomica/model.py b/atomica/model.py
index f664664d..432bb7a3 100644
--- a/atomica/model.py
+++ b/atomica/model.py
@@ -230,14 +230,14 @@ def __setitem__(self, key, value) -> None:
 class Compartment(Variable):
     """A class to wrap up data for one compartment within a cascade network."""
 
-    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
         Variable.__init__(self, pop=pop, id=(pop.name, name))
         self.units = "Number of people"
         self.outlinks = []
         self.inlinks = []
-        self.stochastic = stochastic
+        self.multinomial = multinomial
         self.rng_sampler = rng_sampler
-        if self.stochastic and rng_sampler is None:
+        if self.multinomial and rng_sampler is None:
             self.rng_sampler = np.random.default_rng()
         self._cached_outflow = None
 
@@ -360,7 +360,7 @@ def resolve_outflows(self, ti: int) -> None:
 
         self._cached_outflow = 0
 
-        if self.stochastic:
+        if self.multinomial:
             # for every person in the compartment, we now have probabilities for each possible outflow or staying as the remainder
             stays = max(0.0, 1.0 - outflow)
             rescaled_frac_flows = [rescale * link._cache for link in self.outlinks] + [stays]
@@ -414,7 +414,7 @@ def connect(self, dest, par) -> None:
 
 
 class JunctionCompartment(Compartment):
-    def __init__(self, pop, name: str, duration_group: str = None, stochastic: bool = False, rng_sampler=None):
+    def __init__(self, pop, name: str, duration_group: str = None, multinomial: bool = False, rng_sampler=None):
         """
 
         A TimedCompartment has a duration group by virtue of having a `.parameter` attribute and a flush link.
@@ -433,7 +433,7 @@ def __init__(self, pop, name: str, duration_group: str = None, stochastic: bool
         :param duration_group: Optionally specify a duration group
 
         """
-        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
         self.duration_group = duration_group  #: Store the name of the duration group, if the junction belongs to one
 
     def connect(self, dest, par) -> None:
@@ -535,7 +535,7 @@ def balance(self, ti: int) -> None:
         outflow_fractions /= total_outflow  # proportionately accounting for the total outflow downscaling
 
         # assign outflow stochastically to compartments based on the probabilities
-        if self.stochastic and net_inflow > 0.0:
+        if self.multinomial and net_inflow > 0.0:
             if self.duration_group:
                 raise Exception("Need to implement code here to make sure outflows from EACH part of the duration group are integers")
             else:
@@ -568,7 +568,7 @@ def initial_flush(self) -> None:
             outflow_fractions /= np.sum(outflow_fractions)
 
             # assign outflow stochastically to compartments based on the probabilities
-            if self.stochastic:
+            if self.multinomial:
                 outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / self.vals[0]
 
             # Assign the inflow directly to the outflow compartments
@@ -621,7 +621,7 @@ def balance(self, ti: int) -> None:
             outflow_fractions /= total_outflow
 
         # assign outflow stochastically to compartments based on the probabilities
-        if self.stochastic and net_inflow > 0.0:
+        if self.multinomial and net_inflow > 0.0:
             if np.abs(np.round(net_inflow) - net_inflow).sum() > 1e-3:
                 raise Exception(f"Net inflow should be an integer for a stochastic run, was: {net_inflow}")
 
@@ -638,7 +638,7 @@ def balance(self, ti: int) -> None:
 
         # Finally, assign the inflow to the outflow proportionately accounting for the total outflow downscaling
         for frac, link in zip(outflow_fractions, self.outlinks):
-            if link.parameter is None and total_outflow < 1 and not self.stochastic:  # If we are stochastic then we have already taken the residual link into account
+            if link.parameter is None and total_outflow < 1 and not self.multinomial:  # If we are stochastic then we have already taken the residual link into account
                 flow = net_inflow - np.sum(outflow, axis=1)  # Sum after multiplying by outflow fractions to reduce numerical precision errors and enforce conserved quantities more accurately
             else:
                 flow = net_inflow * frac
@@ -672,7 +672,7 @@ def initial_flush(self) -> None:
                 has_residual = False
 
             # assign outflow stochastically to compartments based on the probabilities
-            if self.stochastic and total_outflow > 0.0:
+            if self.multinomial and total_outflow > 0.0:
                 outflow_fractions = self.rng_sampler.multinomial(self.vals[0], outflow_fractions, size=1)[0] / (self.vals[0] * sum(outflow_fractions))
 
             # Assign the inflow directly to the outflow compartments
@@ -705,8 +705,8 @@ class SourceCompartment(Compartment):
 
     """
 
-    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
-        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         super().preallocate(tvec, dt)
@@ -715,7 +715,7 @@ def preallocate(self, tvec: np.array, dt: float) -> None:
     def resolve_outflows(self, ti: int) -> None:
         # Unlike other Compartments, link._cache will still be in Number form and not converted to a proportion
         for link in self.outlinks:
-            if self.stochastic:
+            if self.multinomial:
                 link.vals[ti] = stochastic_rounding(link._cache, rng_sampler=self.rng_sampler)
             else:
                 link.vals[ti] = link._cache
@@ -725,8 +725,8 @@ def update(self, ti: int) -> None:
 
 
 class SinkCompartment(Compartment):
-    def __init__(self, pop, name, stochastic: bool = False, rng_sampler=None):
-        super().__init__(pop, name, stochastic=stochastic, rng_sampler=rng_sampler)
+    def __init__(self, pop, name, multinomial: bool = False, rng_sampler=None):
+        super().__init__(pop, name, multinomial=multinomial, rng_sampler=rng_sampler)
 
     def preallocate(self, tvec: np.array, dt: float) -> None:
         """
@@ -782,7 +782,7 @@ def update(self, ti: int) -> None:
 
 
 class TimedCompartment(Compartment):
-    def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic: bool = False, rng_sampler=None):
+    def __init__(self, pop, name: str, parameter: ParsetParameter, multinomial: bool = False, rng_sampler=None):
         """
         Instantiate the TimedCompartment
 
@@ -795,7 +795,7 @@ def __init__(self, pop, name: str, parameter: ParsetParameter, stochastic: bool
 
         """
 
-        Compartment.__init__(self, pop=pop, name=name, stochastic=stochastic, rng_sampler=rng_sampler)
+        Compartment.__init__(self, pop=pop, name=name, multinomial=multinomial, rng_sampler=rng_sampler)
         self._vals = None  #: Primary storage, a matrix of size (duration x timesteps). The first row is the one that people leave from, the last row is the one that people arrive in
         self.parameter = parameter  #: The parameter to read the duration from - this needs to be done after parameters are precomputed though, in case the duration is coming from a (constant) function
         self.flush_link = None  #: Reference to the timed outflow link that flushes the compartment. Note that this needs to be set externally because Links are instantiated after Compartments
@@ -900,7 +900,7 @@ def resolve_outflows(self, ti: int) -> None:
         """
 
         # TODO adjust for stochastic variables
-        assert self.stochastic == False, "resolve_outflows needs to be updated to work with stochastic transitions"
+        assert self.multinomial == False, "resolve_outflows needs to be updated to work with stochastic transitions"
 
         # First, work out the scale factors as usual
         self.flush_link._cache = 0.0  # At this stage, no outflow goes via the flush link
@@ -1587,7 +1587,7 @@ class Population:
     Each model population must contain a set of compartments with equivalent names.
     """
 
-    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, stochastic: bool = False, rng_sampler=None):
+    def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_type: str, multinomial: bool = False, rng_sampler=None):
         """
         Construct a Population
 
@@ -1601,7 +1601,7 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         :param label: The full name of the population
         :param progset: A ProgramSet instance
 
-        :param stochastic: whether to use stochasticity in determining population movement
+        :param multinomial: whether to use stochasticity in determining population movement
         :param rng_sampler: optional Generator for random numbers to give consistent results between model runs
 
         """
@@ -1609,7 +1609,7 @@ def __init__(self, framework, name: str, label: str, progset: ProgramSet, pop_ty
         self.name = name  #: The code name of the population
         self.label = label  #: The full name/label of the population
         self.type = pop_type  #: The population's type
-        self.stochastic = stochastic  #: Whether the population should be handled stochastically as discrete integer people instead of fractions
+        self.multinomial = multinomial  #: Whether the population should be handled stochastically as discrete integer people instead of fractions
         self.rng_sampler = rng_sampler
 
         self.comps = list()  #: List of Compartment objects
@@ -1821,17 +1821,17 @@ def build(self, framework, progset):
         for comp_name in list(comps.index):
             if comps.at[comp_name, "population type"] == self.type:
                 if comp_name in residual_junctions:
-                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(ResidualJunctionCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "is junction"] == "y":
-                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
+                    self.comps.append(JunctionCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, duration_group=comps.at[comp_name, "duration group"]))
                 elif comps.at[comp_name, "duration group"]:
-                    self.comps.append(TimedCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
+                    self.comps.append(TimedCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler, parameter=self.par_lookup[comps.at[comp_name, "duration group"]]))
                 elif comps.at[comp_name, "is source"] == "y":
-                    self.comps.append(SourceCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(SourceCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
                 elif comps.at[comp_name, "is sink"] == "y":
-                    self.comps.append(SinkCompartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(SinkCompartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
                 else:
-                    self.comps.append(Compartment(pop=self, name=comp_name, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
+                    self.comps.append(Compartment(pop=self, name=comp_name, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
 
         self.comp_lookup = {comp.name: comp for comp in self.comps}
 
@@ -2023,7 +2023,7 @@ def report_characteristic(charac, n_indent=0):
 
         # Otherwise, insert the values
         for i, c in enumerate(comps):
-            if self.stochastic:
+            if self.multinomial:
                 c[0] = stochastic_rounding(max(0.0, x[i]), rng_sampler=self.rng_sampler)  # integer values (still represented as floats)
             else:
                 c[0] = max(0.0, x[i])
@@ -2052,7 +2052,7 @@ def __init__(self, settings, framework, parset, progset=None, program_instructio
         self.program_instructions = sc.dcp(program_instructions)  # program instructions
         self.t = settings.tvec  #: Simulation time vector (this is a brand new instance from the `settings.tvec` property method)
         self.dt = settings.sim_dt  #: Simulation time step
-        self.stochastic = settings.stochastic  #: Run with discrete numbers of people per compartment instead of fractions.
+        self.multinomial = settings.multinomial  #: Run with discrete numbers of people per compartment instead of fractions.
 
         self.acceptance_criteria = sc.dcp(acceptance_criteria)  #: If the model fails to adhere to these criteria, raise a BadInitialization exception
 
@@ -2203,7 +2203,7 @@ def build(self, parset):
 
         # First construct populations
         for k, (pop_name, pop_label, pop_type) in enumerate(zip(parset.pop_names, parset.pop_labels, parset.pop_types)):
-            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, stochastic=self.stochastic, rng_sampler=self.rng_sampler))
+            self.pops.append(Population(framework=self.framework, name=pop_name, label=pop_label, progset=self.progset, pop_type=pop_type, multinomial=self.multinomial, rng_sampler=self.rng_sampler))
             self._pop_ids[pop_name] = k
 
         # Expand interactions into matrix form
diff --git a/atomica/project.py b/atomica/project.py
index 29c6f6ed..5f48d6e2 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -44,11 +44,24 @@
 
 
 class ProjectSettings:
-    def __init__(self, sim_start=2000, sim_end=2035, sim_dt=0.25, stochastic=False):
+    def __init__(self, sim_start:float=2000, sim_end:float=2035, sim_dt:float=0.25, multinomial:bool=False):
+        """
+        Project settings object
+
+        This contains various settings that are used by the Project when running simulations.
+        Most importantly, it contains the simulation start year, end year, and timestep, with
+        validation performed when setting these values.
+
+        :param sim_start: Simulation start year
+        :param sim_end: Simulation end year
+        :param sim_dt: Simulation time step
+        :param multinomial: Use multinomial mode with stochastic transitions
+        """
+
         self._sim_start = sim_start
         self._sim_dt = sim_dt
         self._sim_end = 0.0
-        self._stochastic = stochastic
+        self.multinomial = multinomial
         self.update_time_vector(end=sim_end)
 
     def __repr__(self):
@@ -68,10 +81,6 @@ def sim_end(self):
     def sim_dt(self):
         return self._sim_dt
 
-    @property
-    def stochastic(self):
-        return self._stochastic
-
     @sim_start.setter
     def sim_start(self, sim_start):
         self._sim_start = sim_start
@@ -88,10 +97,6 @@ def sim_dt(self, sim_dt):
         self._sim_dt = sim_dt
         self.sim_end = self.sim_end  # Call the setter function to change sim_end if it is no longer valid
 
-    @stochastic.setter
-    def stochastic(self, stochastic: bool = False):
-        self._stochastic = stochastic
-
     @property
     def tvec(self) -> np.ndarray:
         """
diff --git a/tests/test_tox_markovchain.py b/tests/test_tox_markovchain.py
index e5d723bd..c526cea2 100644
--- a/tests/test_tox_markovchain.py
+++ b/tests/test_tox_markovchain.py
@@ -14,7 +14,7 @@
 def test_functional_all_framework_databook(demo):
     # Test that all demo frameworks can be run with stochastic mode turned on
     P = at.demo(demo, do_run=False)
-    P.settings.stochastic = True
+    P.settings.multinomial = True
     P.run_sim(result_name="DTMC")
 
 
@@ -105,7 +105,7 @@ def run_simple_percent_test(P, baseline_results, mc_results):
 
     baseline_results = [res]
 
-    P.settings.stochastic = True
+    P.settings.multinomial = True
     mc_results = [None] * N
     with at.Quiet():  # Using at.Quiet should automatically reset the logging level if an exception occurs
         for i in range(N):

From d1c222309748cac89ae82f61500f08c4ea44f4ca Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 2 Sep 2024 23:47:25 +1000
Subject: [PATCH 28/33] Prefill interaction tables with Y/N more logically

Previous "if value" didn't make sense as an object - was always True, using (value.vals of value.assumption) will still generally return True based on defaults, but allows prefilling interaction tables when generating a databook if data has been specified already
---
 atomica/excel.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/atomica/excel.py b/atomica/excel.py
index 2a0119ce..b7801bc9 100644
--- a/atomica/excel.py
+++ b/atomica/excel.py
@@ -834,7 +834,7 @@ def _write_pop_matrix(self, worksheet, start_row, formats, references: dict = No
             from_idx = self.from_pops.index(from_pop)
             to_idx = self.to_pops.index(to_pop)
             if boolean_choice:
-                value = "Y" if value else "N"
+                value = "Y" if (value.vals or value.assumption) else "N"
             content[from_idx, to_idx] = value
 
         # Write the content

From c7520868b09d20bf9c3dd5de1d370d41cbec98c9 Mon Sep 17 00:00:00 2001
From: Rowanmh <rowanmh@gmail.com>
Date: Mon, 2 Sep 2024 23:48:13 +1000
Subject: [PATCH 29/33] More informative sc.prepr output (useful when
 debugging)

---
 atomica/excel.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/atomica/excel.py b/atomica/excel.py
index b7801bc9..fd2a43ae 100644
--- a/atomica/excel.py
+++ b/atomica/excel.py
@@ -424,7 +424,8 @@ def __init__(self, code_name: str, full_name: str, tvec: np.array, from_pops: li
             raise Exception('Unknown TimeDependentConnections type - must be "transfer" or "interaction"')
 
     def __repr__(self):
-        return '<TDC %s "%s">' % (self.type.title(), self.code_name)
+        # return '<TDC %s "%s">' % (self.type.title(), self.code_name)
+        return sc.prepr(self)
 
     @classmethod
     def from_tables(cls, tables: list, interaction_type):

From 59851fffddb40d7becc731419365c06a854a8643 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 28 Feb 2025 11:56:39 +1100
Subject: [PATCH 30/33] Fix save divide for pop_aggregation == "weighted" in
 PlotData

---
 atomica/plotting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/atomica/plotting.py b/atomica/plotting.py
index 9529ab8d..0cab8ca3 100644
--- a/atomica/plotting.py
+++ b/atomica/plotting.py
@@ -392,7 +392,7 @@ def placeholder_pop():
                         elif pop_aggregation == "weighted":
                             numerator = sum(aggregated_outputs[x][output_name] * popsize[x] for x in pop_labels)  # Add together all the outputs
                             denominator = sum([popsize[x] for x in pop_labels])
-                            vals = np.divide(numerator, denominator, out=np.full(numerator.shape, np.nan, dtype=float), where=numerator != 0)
+                            vals = np.divide(numerator, denominator, out=np.full(numerator.shape, np.nan, dtype=float), where=denominator != 0)
                         else:
                             raise Exception("Unknown pop aggregation method")
                         self.series.append(Series(tvecs[result_label], vals, result_label, pop_name, output_name, data_label[output_name], units=aggregated_units[output_name], timescale=aggregated_timescales[output_name], data_pop=pop_name))

From 5e42910a3cfeafcf8a5f538f98af4d560c04ebac Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 28 Feb 2025 11:58:09 +1100
Subject: [PATCH 31/33] Plotdata.set_colors add option for colors=tuple with
 (color, opacity)

---
 atomica/plotting.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/atomica/plotting.py b/atomica/plotting.py
index 0cab8ca3..2b7e9c05 100644
--- a/atomica/plotting.py
+++ b/atomica/plotting.py
@@ -977,7 +977,9 @@ def set_colors(self, colors=None, results="all", pops="all", outputs="all", over
         elif isinstance(colors, list):
             assert len(colors) == len(targets), "Number of colors must either be a string, or a list with as many elements as colors to set"
             colors = colors
-        elif colors.startswith("#") or colors not in [m for m in plt.cm.datad if not m.endswith("_r")]:
+        elif isinstance(colors, tuple) and colors[0].startswith("#"): # tuple with (color, opacity)
+            colors = [colors for _ in range(len(targets))]  # Apply color to all requested outputs
+        elif isinstance(colors, str) and (colors.startswith("#") or colors not in [m for m in plt.cm.datad if not m.endswith("_r")]):
             colors = [colors for _ in range(len(targets))]  # Apply color to all requested outputs
         else:
             color_norm = matplotlib_colors.Normalize(vmin=-1, vmax=len(targets))

From 86d19a66bc2b1c66abe564c2aeba9b8d92616e74 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 28 Feb 2025 11:59:50 +1100
Subject: [PATCH 32/33] plotting.py plot_series add optional kwarg colors=None,
 which gets passed to plotdata.set_colors

---
 atomica/plotting.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/atomica/plotting.py b/atomica/plotting.py
index 2b7e9c05..b77b9b6d 100644
--- a/atomica/plotting.py
+++ b/atomica/plotting.py
@@ -1441,7 +1441,8 @@ def process_input_stacks(input_stacks, available_items):
     return figs
 
 
-def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=None, lw=None, n_cols: int = None) -> list:
+def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=None, lw=None, n_cols: int = None,
+                colors=None) -> list:
     """
     Produce a time series plot
 
@@ -1454,6 +1455,7 @@ def plot_series(plotdata, plot_type="line", axis=None, data=None, legend_mode=No
     :param lw: override the default line width
     :param n_cols: If None (default), separate figures will be created for each axis. If provided, axes will be tiled as subplots in a single figure
                    window with the requested number of columns
+    :param colors: Colors to be passed to plotdata.set_colors
     :return: A list of newly created Figures
 
     """
@@ -1507,7 +1509,7 @@ def _prepare_figures(dim1, dim2, n_cols):
         logger.warning("At least one Series has only one timepoint. Series must have at least 2 time points to be rendered as a line - `plot_bars` may be more suitable for such data")
 
     if axis == "results":
-        plotdata.set_colors(results=plotdata.results.keys())
+        plotdata.set_colors(colors=colors, results=plotdata.results.keys())
 
         figs, axes = _prepare_figures(plotdata.pops, plotdata.outputs, n_cols)
 
@@ -1542,7 +1544,7 @@ def _prepare_figures(dim1, dim2, n_cols):
                 _render_legend(ax, plot_type)
 
     elif axis == "pops":
-        plotdata.set_colors(pops=plotdata.pops.keys())
+        plotdata.set_colors(colors=colors, pops=plotdata.pops.keys())
 
         figs, axes = _prepare_figures(plotdata.results, plotdata.outputs, n_cols)
 
@@ -1575,7 +1577,7 @@ def _prepare_figures(dim1, dim2, n_cols):
                 _render_legend(ax, plot_type)
 
     elif axis == "outputs":
-        plotdata.set_colors(outputs=plotdata.outputs.keys())
+        plotdata.set_colors(colors=colors, outputs=plotdata.outputs.keys())
 
         figs, axes = _prepare_figures(plotdata.results, plotdata.pops, n_cols)
 

From aab2c7c41b2db796ce3d1db4a920049fd9c83661 Mon Sep 17 00:00:00 2001
From: Kelvin Burke <kelv.burke@gmail.com>
Date: Fri, 28 Feb 2025 12:01:06 +1100
Subject: [PATCH 33/33] ProjectSettings.tvec update to latest function using
 arange - to match other branches

---
 atomica/project.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/atomica/project.py b/atomica/project.py
index 5f48d6e2..8b800381 100644
--- a/atomica/project.py
+++ b/atomica/project.py
@@ -108,8 +108,8 @@ def tvec(self) -> np.ndarray:
         :return: Array of simulation times
 
         """
-
-        return np.linspace(self.sim_start, self.sim_end, int((self.sim_end - self.sim_start) / self.sim_dt) + 1)
+        return np.arange(round((self.sim_end - self.sim_start) / self.sim_dt) + 1) * self.sim_dt + self.sim_start
+        # return np.linspace(self.sim_start, self.sim_end, int((self.sim_end - self.sim_start) / self.sim_dt) + 1)
 
     def update_time_vector(self, start: float = None, end: float = None, dt: float = None) -> None:
         """