From baeafb342d5caf4bab175cd70dd90819b4570605 Mon Sep 17 00:00:00 2001 From: Rowanmh Date: Tue, 12 Oct 2021 03:38:06 +1100 Subject: [PATCH 01/38] 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 Date: Tue, 12 Oct 2021 11:56:25 +1100 Subject: [PATCH 02/38] 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() Date: Sat, 6 Nov 2021 03:28:56 +1100 Subject: [PATCH 03/38] 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 Date: Sat, 6 Nov 2021 03:58:48 +1100 Subject: [PATCH 04/38] 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 Date: Sat, 6 Nov 2021 04:41:14 +1100 Subject: [PATCH 05/38] 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 Date: Sat, 6 Nov 2021 04:41:17 +1100 Subject: [PATCH 06/38] 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 Date: Mon, 8 Nov 2021 01:33:07 +1100 Subject: [PATCH 07/38] 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() Date: Mon, 8 Nov 2021 01:54:24 +1100 Subject: [PATCH 08/38] 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 Date: Mon, 8 Nov 2021 03:56:09 +1100 Subject: [PATCH 09/38] 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 Date: Mon, 8 Nov 2021 05:03:09 +1100 Subject: [PATCH 10/38] 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 Date: Mon, 8 Nov 2021 05:03:40 +1100 Subject: [PATCH 11/38] 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 Date: Mon, 8 Nov 2021 23:31:08 +1100 Subject: [PATCH 12/38] 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 Date: Wed, 1 Mar 2023 16:21:39 +1100 Subject: [PATCH 13/38] 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 Date: Fri, 3 Mar 2023 21:49:49 +1100 Subject: [PATCH 14/38] 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 Date: Wed, 1 Nov 2023 16:52:25 +1100 Subject: [PATCH 15/38] 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 Date: Fri, 3 Nov 2023 08:42:18 +1100 Subject: [PATCH 16/38] 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 Date: Fri, 3 Nov 2023 09:26:01 +1100 Subject: [PATCH 17/38] 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 Date: Mon, 6 Nov 2023 11:12:11 +1100 Subject: [PATCH 18/38] 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 Date: Fri, 24 Nov 2023 15:12:05 +1100 Subject: [PATCH 19/38] 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 Date: Fri, 24 Nov 2023 15:16:41 +1100 Subject: [PATCH 20/38] 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 Date: Thu, 11 Apr 2024 13:32:17 +1000 Subject: [PATCH 21/38] 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 Date: Fri, 7 Jun 2024 15:25:24 +1000 Subject: [PATCH 22/38] 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 Date: Fri, 7 Jun 2024 15:41:34 +1000 Subject: [PATCH 23/38] 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 Date: Fri, 7 Jun 2024 15:42:13 +1000 Subject: [PATCH 24/38] 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 Date: Fri, 7 Jun 2024 15:42:36 +1000 Subject: [PATCH 25/38] 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 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>NdSjHl3xgtIdIB@t-rrL0_d4$ezv!`v7S;tBZ85}UN_bNpXqQS+ zH5Bqjr*BfW^RVD%oE{U=tBSYQwUYP8Vbi41(ik2xx8>gp~y@iEjpgb_r5#0VhtUZ9oA zpHJKkL{&DBnZJhNq3Dlt*8g;!*I)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+rInjsLi0d!M`XI@kP%s%PcXw*pwQY;}m2tje;qK4S z>T>^|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&GvI3w9 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*YEhlrKsa5^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^i3***1d}!ZwNK=FA8JM4l zrb`^@p8At#jSz6;uSr+^B@@d~DzU98=bsEY#8 z;cP01apZi_9gXmVk6+*!B-GY`heZ`Pwq*WE*6Fuh=NC403p&V|~ z)%P;6@fVKDF?9jK$6E{a4 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 z@uMyl4ml6qCRiDB-j(}EfRj`?dE4!(jfb-NTjem03Z%uhqY zK5|@Cs|ic{wQ)#Pw4i0Ncj$FRs>_Yb=$RH zp;}{Ny(bR+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$>gHTaU#YkuU1fvldk$nXD8^Z0^%n5p+(qcieBqD^+ro5StWY6)tWf#5`Aae- zMrA$Y`|00YFolUy*4PyX<0FLWd&dGxZGH&*w-Hywp$Gc#P&u}5*0&D{H$&BN!U*-XS zSgoJ<#c_ZYG0aGEYi&+gRIfLwxe1Tme7ttJf$;;HK}>|%JDi~Mz}8XXS6R^Y zm{dvP65hZM4lU{)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*Kc0791fnWeVKiCa(|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=gHhzJ~m zhula6TH4TOO0!V%w0-Z7`(c!rs$n*0NW85~RJHQBBk;Y#6RugeQZ 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|07dwM9L5c-mMFEhNIyC3$5#1rLqp>wRYc2ZS@^s6!a$*;-j_vgz`1%R+=B z<;h4iWojWXSQ&BiGS39?Si!mWx7fm%$e`l~qINesR@^&iP#&dWZ5r=#_fuC+8s~r|W``rOl816`SbU4T_VQ*%!qBhWc}W1HW$6e-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+ouQ-!d zZ+vyDP5Pcg2q0gkPjKe>_^|Loc0al{KJ8Sx$1P=0151Y64s6u4J#f?)-GEJs(8_Hm zmcr~$?)h2y6>#Z_7T}04ru@#9jQcj7efUGQE z$gk;af=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_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=_opKucy4B&X&>*A1I+of!BcEdT9pA|G1Z{ zK*tw2LV+f+|ACgNPf 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@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_?O0afh6ofkE0TYsGcBxMN!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+Brj!j7xifFF-1dSa@+x7igvBrwPnZ4)&EBj9B;z;vjHHJ)mj;Fur{rbhu+hB;24EDC%m8dmD`;+~b4=H4fW{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@(zzkC(9J>I_8hY!=m9}|d$g?z7zOZV;9hHXBq-Tzpf1iU;wnq7k$oeTvZ`V%uGpFRuw z{~VkA_$xR`5h25h&gS!k8Mz9Zt(d9DMM+ zHdDwW;(=Zvo-x4U?_E{@=HJHU67QufLUsx6q(IEma0W6hs))QA=Wz+*&u!$wF@#*3 z>w?EV11QgoOzJDdvO1c<8U@nX4#Ti5@L=V_-;sGl)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*5952kBHZl2ngJlN z91#koConQm)5uSoFjy~C-GGJj*g;X?=4{)g1TU64GP1>P1=4@$k9htxYTxg%Jzgfv z$l%OidfHu5rw+qTHjkZNkA!qZivpgWO}05=R1>PhA_q=zz(WRhIlGb@(pYmgnkgTvJb(}d3q z7RwmoM6Fr+5D!Z$pQKfm{Wf*t0ehyArY;5X;rFz4TCYG2g|9|Q^h2yKA3qgMVzeb_ zPgLUzNSQXt4C@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$y8k0*vqr)Mr#`1=(Lew#f9!`cqAeWJG@I9c0%)5CDF+ousgE z!NR3*u4sC{nKeh?G&dzc(P@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^wxMwDXgasMh^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`-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-e0P$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=?&`W9NiAYTikDDZxO>>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)kbF!84Tn4 zpJA@N<2JaMi5IJZH=6{Wg@pL$Xf{7>Tta^s=>bTELQYvJTn22A^M`9@`5qqvOiF5h#V zaIB3(9=HH&E+Nnunp(Zw4qQz-cbx 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_CV6rF+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@{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^v0W}T#6wFE7J&@kv~qBkov6zbpaDxAv>&!onHudq&F zj9TECP;>nZ$X#=iZ^xi}tiNz%Z2+v>QiNe*p1df|I`dc+R!msHk4Jgw3Nj82cifD(;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}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=W0K2)4;tkrkJ@9t?PA8#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;l+V1q^YM9SH@=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@?KusIYG0yjHbwH_P5)HeTRSET4J)(@7hcWJ%ga zlDu@kOSl_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(l6oxq@{>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_lc40%__7pGHS3SzS7VfJ8a}yQ zQO40KO&N+m=f;b;8^QPFd0a895ceBQH$?KIRkUSk|A-X@bd1jn#kmu?*-uN;w5M*| zsW6h^6G$?Q@8-W_OFTZtq)aYq#S`MHDRHSlm3p>!NNNuk58W_PTuli*Ft3 zk0OZugvY4fvB6L)OQ1VZ6~bVq;bYB7xD;c!y2e=6XkF~vxMC^X389Yp5b=&mDF;Ssqy#XiI(8it(2NcXE_+Nwg5PS)WN}$Yh7xL zKtsEUUam4K%?~RsrmOQdO?hjiz)ybsbHRIa1t@oDr3a)^VVKi{#~yEpC2x@6ZF=m6 zV?;+N!0nU8Ip_Pv3>E)KrF8LSJc9C{h9dh%zwRDFmhX#CQIkY9wqYX zGmqV0@n#Z+wfvY*ph5O!eeuN~+qkus5KyCsef$p>nbs6RZ_S{}iG`{gOk@H#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}_ihCm%=*=9BFfED+l zRnYjgI$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_ptE33eq1vz-+TN3z3O~)c zxGax13h}gE)F}5^S8(!Tw;>QNvJ_^+a7rJyKqM~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 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{$ zxhXBEmp(-h`n5;+Kle5aB_u84S1`QQ{e15xmDNFaWh8;pgFG`?QKO{asV>TpH(s=o zv21Wy^^};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%V8Y3&+KlSEk=<2Z^*(T`dAMHL6|U3 zgE=mP0tAV(7qhzvW3&^+gWmCzJs!QPz9)}P2Xogm-fe1krBB~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 z=s}72?n-JCwZ)W|YRs5$Fb0@XsG4}v&`S{ziL$tTJE9D2Yc|Z|XwxO@bCy!>8^y|e7`VGdBFqDiO z7H-2xx}NU!$D^M424M5<-O2IMeXxWmW^G>+LyCo!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&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@v0NDjV2h*&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=i07MhMh($oO}!t$oD>^cLfpr(>mQ0SUeI6eP<8SRL5((cncXEQx>fBw z&`pgx>$TQ_+7|8h|9q5_JFH|L)YuTc#`8dT&GR%0Sj<4xje0JJATvEwAkLOl*=$T} zqGr4P9+w?S^`~|nO}2K)jvlgL^!&tVKa^zix(Fi&2@+!l1I?l4&8A>7lTM}y&s9NhwGA05 zniXpJ&vbqN&m*e4kz5P#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+8m;DsTst0|N}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$*W8#j(XBbtf&zq#I0@0 zJv+%txrJe*6;UBZv6wAwN#We=1gdE?pV_`%Q(tq#3zE*P$V7XYztoCl+|K}A!(i1&}T86 z@Ks|}PdJN9@$EnKlsa zcy==oyhgy1OZjIot8oAQ3D&JGgovM-!)BuvAp|bje1ZcL;$&uHRX5~$y?-KWDI}`2 zE7Z0dB6|}mu8s@d?`pzLdD|qACvMgVrZXuC8c(@S3*^BBX(2~F>3M~?Y=|0 z?|3}bdM_!L5fIBu(w|$9V*Bh7N1*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-d$&(1T7YnPL@ZQ8C)Y^ZdS6}D!FlGMF`P=Q{_q6$B9jubojtu2kAeBA`iEA4) zyjbE%gdbGiN9Q4QS!l{`^zBAIKI5@Mqpi=G9le!(x>y{6S@EO4zkl;fRe>U?g@0yL z6qEk7AK6x=vRG^Zg(pxtWz4Z3uk!r4jh$;x)c`p1p%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!8Xa z#)pk3;|~NYjQ61zWMu58FN~@*@71e=A!w@|)2>e-aomvl$zy>xah~1~I#3r*Qvk5y z!hr%-hc)$B61r+!ByYz!SMSEYo$|a%>e5%e@V}}2FiIdw$U4m`=LRb z5H!D`>k=dm$8PI=x#}(wJ%Vcp6qGqLOdJ;kKw3+~Wo?@?eqw9fjO;i{%rSSR&3Fz* z=II~N^phox_caId1EqU^^9Nh&dei0BHfqkHEEo&g@ z?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_?*nUQf( zslNk_nCf<8IDwL^Fg^&AU`^NetuxD8W00Bx;Yh7D>bH2MdQ)31xrBGjzDI5%_Qq|i zOl^=J2J;P!2xJ#PPn#(v6%SKk4x*S3vZ?(=1O@QRI<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(<cf93*(E}%LzyPcWIY`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(*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%?DZUhjnU00DAt{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^LLJ<*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%^z>Ntuc7W(0WbG37K3;S%jC+!$^DU3~A0fVy}dX$qSxFOK+TPkusElg~cxq0n)Z zG&`JCo>}T3A8WIKq;!V;b_qEs)`HDTQU$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<4Kzn7q0=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}iq1V!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{H7bsa^&Rn z0R7jwvy4%JD1N8h7mUA@s}+O?8N5>tJM0|(3GGAYaO@o4r`f7<5;~I3Rm4@}G|9ST z?W#&Mra_O#-(KeYih@mbRl#Br6m^xrT5KFJ~ z(i|1`<#fXJQ{s3@N{z|D`;s&Ytm!VphYPk6p3$^ITU4%E-{>*`f1?xhDtA!Tm$xd< zwy|9T>-q=frb@3c2H&9)#q#=(7$pmZ)bSehDCTm`0B}@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?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;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)xzh6ws2uA&M* zLRXVq``>nCAjpQ4H;Z+nUE5lu<5lr0g{3>VjMQ2hEMKhE>YZ9Wqjt)Qu#nb1rg_$r z^>}=tknP`)VD4r}+mx$vAphN;Hwj0%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^U1dt+Ar#oL zfRJ++;+Z`?LV$3qOYOdZwc_%j)QR?|9*Z8e`^G=ttgOi0XDIWdbEo651Vo|puWl)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<(YP^TocfLc-vgo`Lx7OG#9)Tdn&@D5bQD!zYO2nEO|4&VYS4o>2^GNCz^sP~IV|+l5)ES|| zhUcmt(Yp#j8S2!VG&MD?PR`^}%=F1)h&N13s5inQzM%Rz5em2h@fXGd*a;fB!TAeza7CtOmYi}&B?Cak{z4r>DZ|FrJ>H(&toFTnd9AphRY*Nyn!fIh;% zfG%qs;=cc_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 Date: Wed, 26 Jun 2024 16:50:33 +1000 Subject: [PATCH 26/38] 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 a8ef28dcef509adeaab9e99b3c2c2c042ff62849 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Fri, 28 Jun 2024 08:34:59 +1000 Subject: [PATCH 27/38] Enable auto calibration of transfers --- CHANGELOG.md | 4 ++++ atomica/calibration.py | 7 ++++--- atomica/version.py | 4 ++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a4c4fef..c1fefab4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,10 @@ This file records changes to the codebase grouped by version release. Unreleased changes are generally only present during development (relevant parts of the changelog can be written and saved in that section before a version number has been assigned) +## [1.28.2] - 2023-06-28 + +- Enable automated calibration of transfers and updated documentation to cover this feature + ## [1.28.1] - 2023-02-05 - Updated various Pandas operations to improve compatibility with Pandas 2.2.0 diff --git a/atomica/calibration.py b/atomica/calibration.py index e488330e..b8a728b6 100644 --- a/atomica/calibration.py +++ b/atomica/calibration.py @@ -43,7 +43,6 @@ def _update_parset(parset, y_factors, pars_to_adjust): tokens = par_name.split("_from_") par = parset.transfers[tokens[0]][tokens[1]] par.y_factor[pop_name] = y_factors[i] - raise NotImplementedError def _calculate_objective(y_factors, pars_to_adjust, output_quantities, parset, project): @@ -123,7 +122,10 @@ def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, :param pars_to_adjust: list of tuples, (par_name,pop_name,lower_limit,upper_limit) the pop name can be None, which will be expanded to all populations relevant to the parameter independently, or 'all' which will instead operate - on the meta y factor. + on the meta y factor. To calibrate a transfer, the parameter name should be set to + ``'_from_'`` and then the destination population can be specified + as the ``pop_name``. For example, to automatically calibrate an aging transfer 'age' from 5-14 + to 15-64, the tuple would contain ``pars_to_adjust=[('age_from_5-14','15-64',...)]`` :param output_quantities: list of tuples, (var_label,pop_name,weight,metric), for use in the objective function. pop_name=None will expand to all pops. pop_name='all' is not supported :param max_time: If using ASD, the maximum run time @@ -241,7 +243,6 @@ def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, tokens = par_name.split("_from_") par = args["parset"].transfers[tokens[0]][tokens[1]] logger.debug("parset.transfers['{}']['{}'].y_factor['{}']={:.2f}".format(tokens[0], tokens[1], pop_name, par.y_factor[pop_name])) - raise NotImplementedError # Transfers might be handled differently in Atomica args["parset"].name = "calibrated_" + args["parset"].name diff --git a/atomica/version.py b/atomica/version.py index 4e146119..88bb3c68 100644 --- a/atomica/version.py +++ b/atomica/version.py @@ -6,6 +6,6 @@ import sciris as sc -version = "1.28.1" -versiondate = "2024-02-05" +version = "1.28.2" +versiondate = "2024-06-28" gitinfo = sc.gitinfo(__file__) From c6dbf4b9ea8a14bac58da453f8d6a5c3acc51bc6 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Fri, 28 Jun 2024 08:35:08 +1000 Subject: [PATCH 28/38] Update test pipeline --- azure-pipelines.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 4f7ed50a..50fc1d9c 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: @@ -52,7 +52,7 @@ jobs: - task: UsePythonVersion@0 displayName: 'Select python version' inputs: - versionSpec: '3.11' + versionSpec: '3.12' - script: python -m pip install flake8 flake8-junit-report displayName: 'Install flake8' - script: flake8 atomica tests --exit-zero --output-file flake8.txt @@ -75,7 +75,7 @@ jobs: - task: UsePythonVersion@0 displayName: 'Select python version' inputs: - versionSpec: '3.11' + versionSpec: '3.12' - script: | sudo apt install pandoc -y displayName: 'Install pandoc' @@ -113,7 +113,7 @@ jobs: - task: UsePythonVersion@0 displayName: 'Select python version' inputs: - versionSpec: '3.11' + versionSpec: '3.12' - script: | python -m pip install wheel python -m pip install twine From 5af7e4be749302ccc327df609da9edce9bbcd328 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Fri, 28 Jun 2024 08:36:02 +1000 Subject: [PATCH 29/38] Update documentation --- docs/tutorial/T2-Calibration.ipynb | 232 ++++++++++++++++++++++++++++- 1 file changed, 224 insertions(+), 8 deletions(-) diff --git a/docs/tutorial/T2-Calibration.ipynb b/docs/tutorial/T2-Calibration.ipynb index df6bdbd5..9c0dbab2 100644 --- a/docs/tutorial/T2-Calibration.ipynb +++ b/docs/tutorial/T2-Calibration.ipynb @@ -83,7 +83,7 @@ "\n", "\n", " \n", - "An example of a suitable parameter for databook calibration is `relative seasonality of mosquito population size` - no hard data can exist, but it might be possible to calculate proxy values from other data such as rainfall patterns in different years, adjust annual values manually to match historical epidemic patterns, and then assume a median value for future years. Having this assumption in the databook allows for those calculations to be used as a starting point in the databook before manually editing, and allows for comparability with other projects that use the same parameter. \n", + "An example of a suitable parameter for databook calibration is `relative seasonality of mosquito population size` - it is unlikely that any existing data would be able to directly provide a value for this parameter, but it might be possible to calculate proxy values from other data such as rainfall patterns in different years, adjust annual values manually to match historical epidemic patterns, and then assume a median value for future years. Having this assumption in the databook allows for those calculations to be used as a starting point in the databook before manually editing, and allows for comparability with other projects that use the same parameter. \n", "\n", "
\n", "Suggestion: When designing a databook, it is recommended that all parameters intended for explicit databook calibration are placed on a single 'Calibration' sheet to provide clarity about what is data and what are calibrated assumptions.\n", @@ -91,14 +91,22 @@ "\n", "An example of a suitable parameter for scale factor calibration is `treatment initiation` used to determine model transitions from diagnosed to treated - a country has reported data for the number of people initiating treatment, and it is important to accurately maintain the official numbers in the databook. Nevertheless, there may be systematic under-reporting by an unknown degree and we want to account for those additional treatments in the model to ensure outcomes are representative, so it is appropriate to adjust the scale factor.\n", "\n", - "An example of a parameter that could be adjusted in either way depending on the circumstances or just personal preference would be `force of infection` - this is a clear calibration parameter that is not based on data, and could be adjusted in a databook if it's necessary to reflect changing circumstances external to the model over time, calibrated automatically with a scale factor via the `calibrate` function below in order to achieve the best fit to data, or even a mixture of the two methods.\n", + "An example of a parameter that could be adjusted in either way depending on the circumstances or just personal preference would be `force of infection` - this is a clear calibration parameter that is not based on data, and could be adjusted in a databook if it's necessary to reflect changing circumstances external to the model over time, calibrated automatically with a scale factor via the `calibrate` function below in order to achieve the best fit to data, or even a mixture of the two methods." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calibration scale factors\n", "\n", + "The values in the databook are intended to serve as a record of the original input data used to inform the model. Therefore, if a parameter needs to be scaled up or down to calibrate the model, it is typically preferable to do this by setting a 'scale factor' (also referred to as a 'y-factor') in the `ParameterSet`, rather than changing the data. \n", "\n", "
\n", "The web interfaces (such as the Cascade Analysis Tool) perform calibration using scale factors. The scale factors shown on the website correspond to the values being set here.\n", "
\n", "\n", - "To set a scale factor, create a `ParameterSet` either by copying an existing one, or creating a new one. Then, access the `pars` attribute to look up the parameter you wish to change, and set the `y_factor` for the population you want to change:" + "To set a scale factor, create a `ParameterSet` either by copying an existing one, or by creating a new one. Then, access the `pars` attribute to look up the parameter you wish to change, and set the `y_factor` for the population you want to change:" ] }, { @@ -107,7 +115,7 @@ "metadata": {}, "outputs": [], "source": [ - "p2 = P.parsets[0].copy()\n", + "p2 = P.parsets[0].copy() # Copy an existing ParameterSet\n", "p2.pars['b_rate'].y_factor['adults'] = 2" ] }, @@ -133,7 +141,206 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can see that we have considerably overshot the data, indicating that doubling the birth rate was much too big a change. This would typically be the first step in an iterative process, where you adjust the scale factor, inspect the data, and then make further adjustments. \n", + "We can see that we have considerably overshot the data, indicating that doubling the birth rate was much too big a change. This would typically be the first step in an iterative process, where you adjust the scale factor, inspect the data, and then make further adjustments. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calibrating compartment sizes\n", + "\n", + "The `ParameterSet` contains parameters for quantities entered on the 'Compartments' and 'Characteristics' sheet of the databook. These are used to calculate the initial compartment sizes used at the start of the simulation. Y-factors can be set for these quantities, the same as for any other parameter. This will not change the data points - just the calculated initial compartment size. For example, we can initialize the model with half as many susceptible people by setting a y-factor of 0.5 on the `sus` compartment size:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p2 = at.ParameterSet(P.framework, P.data) # Make a new ParameterSet\n", + "p2.pars['sus'].y_factor['adults'] = 0.5 \n", + "r2 = P.run_sim(parset=p2,result_name = 'Decreased initial sus')\n", + "d = at.PlotData([result,r2], outputs='sus',project=P)\n", + "at.plot_series(d,axis='results',data=P.data);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unlike setting y-factors for parameters, y-factors for compartment sizes may result in invalid initializations. This can occur if the scaled compartment sizes or characteristics are inconsistent with each other. For example, if the calibrated number of initial susceptible people results in needing more people than the total population size, an error will occur:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p2 = at.ParameterSet(P.framework, P.data) # Make a new ParameterSet\n", + "p2.pars['sus'].y_factor['adults'] = 1.1\n", + "try:\n", + " r2 = P.run_sim(parset=p2,result_name = 'Decreased initial sus')\n", + "except Exception as e:\n", + " print(e)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case, such an error could be resolved by also increasing the y-factor for the 'alive' characteristic, thereby increasing the total population size:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p2 = at.ParameterSet(P.framework, P.data) # Make a new ParameterSet\n", + "p2.pars['sus'].y_factor['adults'] = 1.1\n", + "p2.pars['alive'].y_factor['adults'] = 1.2\n", + "r2 = P.run_sim(parset=p2,result_name = 'Increased initial sus')\n", + "d = at.PlotData([result,r2], outputs=['sus','alive'],project=P)\n", + "at.plot_series(d,axis='results',data=P.data);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "See [https://atomica.tools/docs/master/general/Compartment-Initialization.html](https://atomica.tools/docs/master/general/Compartment-Initialization.html) for more information on how the initialization of compartment sizes is carried out by Atomica. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calibrating transfers\n", + "\n", + "It is also possible to set y-factors to calibrate transfers, that move people from one population to another. A transfer moves people from one population to another e.g., aging from 5-14 to 15-64. Therefore, a y-factor for a transfer must be identified by two populations rather than one. Further, the transfer parameter does not appear in the framework, because they are population-specific. Instead, the transfer is defined in the databook, and then a corresponding parameter is automatically created in the `ParameterSet`. For example, we can load the TB demo in Atomica and inspect the transfers:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2 = at.demo('tb')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the databook, there are two transfers, relating to aging and incarceration:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.data.transfers" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These both appear in the `ParameterSet` under the 'transfers' attribute:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.parsets[0].transfers.keys()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To set a scale factor, it is necessary to define both the source and destination population. The transfer is indexed first by the source population, and then by the 'to' population. The top level key in the transfer corresponds to the source population:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.parsets[0].transfers['age'].keys() # Source populations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we index the transfers using this code name, we can access a `Parameter` that contains `y_factor` entries for all of the destination populations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.parsets[0].transfers['age']['5-14']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice how this parameter only contains y-factors for the destination populations associated with the 'age' transfer out of the '5-14' population" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.parsets[0].transfers['age']['5-14'].y_factor" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This corresponds to the fact that the aging transfer can only move people in 5-14 population to the 15-64 population. To set a y-factor for this transfer, assign a value to this entry in the `y_factor` attribute:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P2.parsets[0].transfers['age']['5-14'].y_factor['15-64'] = 1.5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The example above will increase the transfer rate for aging from 5-14 to 15-64. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Automatic calibration\n", "\n", "Automated calibration is also available via the project's `calibrate` method. This will automatically adjust parameter values to match the data. To use this function, you need to specify which parameters to set scale factors for, and which variables in the databook to compute calibration quality from. The framework can provide defaults for which parameters to automatically calibrate, or you can pass a list of those parameters in directly. In this example, we will pass in `b_rate` because we want to adjust the birth rate, and we will use `sus` as a measurable since we want to match the number of susceptible people. The configuration therefore corresponds to the example shown above. " ] @@ -182,6 +389,15 @@ "at.plot_series(d,axis='results',data=P.data);" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Automated calibration of transfers\n", + "\n", + "As described above, transfers are identified by two populations rather than one. The 'adjustable' for a transfer parameter is therefore identified by concatenating the name of the transfer and the name of the source population, separated by `'_from_'`. For example, we could specify an adjustable of `'age_from_5-14'` to calibrate the aging transfer out of the 5-14 population. " + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -284,9 +500,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": { @@ -298,7 +514,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.12.2" }, "toc": { "base_numbering": 1, From 9a353e1e8fb5e77c9130daf348d96c687d2710e5 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Fri, 28 Jun 2024 08:37:11 +1000 Subject: [PATCH 30/38] Update setup.py --- setup.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 675bf903..58a58152 100644 --- a/setup.py +++ b/setup.py @@ -21,8 +21,13 @@ "Operating System :: OS Independent", "Programming Language :: Python", "Topic :: Software Development :: Libraries :: Python Modules", - "Development Status :: 3 - Alpha", + "Development Status :: 4 - Beta", "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", ] setup( @@ -41,7 +46,7 @@ include_package_data=True, install_requires=[ "matplotlib", - "numpy>=1.10.1", + "numpy>=1.10.1,<2", "scipy>=1.2.1", "pandas", "xlsxwriter", From 6a4970352564747b5c0e2e9f2bbd38377064eb35 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Fri, 28 Jun 2024 10:11:03 +1000 Subject: [PATCH 31/38] Update deprecated scipy function --- atomica/plotting.py | 4 ++-- setup.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/atomica/plotting.py b/atomica/plotting.py index 96362ca5..9529ab8d 100644 --- a/atomica/plotting.py +++ b/atomica/plotting.py @@ -447,9 +447,9 @@ def accumulate(self, accumulation_method) -> None: s.vals = np.cumsum(s.vals) elif accumulation_method == "integrate": if s.timescale: - s.vals = scipy.integrate.cumtrapz(s.vals, s.tvec / s.timescale) + s.vals = scipy.integrate.cumulative_trapezoid(s.vals, s.tvec / s.timescale) else: - s.vals = scipy.integrate.cumtrapz(s.vals, s.tvec) + s.vals = scipy.integrate.cumulative_trapezoid(s.vals, s.tvec) s.vals = np.insert(s.vals, 0, 0.0) # If integrating a quantity with a timescale, then lose the timescale factor diff --git a/setup.py b/setup.py index 58a58152..eb71024a 100644 --- a/setup.py +++ b/setup.py @@ -47,7 +47,7 @@ install_requires=[ "matplotlib", "numpy>=1.10.1,<2", - "scipy>=1.2.1", + "scipy>=1.6", "pandas", "xlsxwriter", "openpyxl", From c2fe4ae8aa394903e5c7c0cf1645d587a62c4c56 Mon Sep 17 00:00:00 2001 From: Romesh Abeysuriya Date: Mon, 19 Aug 2024 11:05:03 +1000 Subject: [PATCH 32/38] 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 Date: Mon, 2 Sep 2024 23:47:25 +1000 Subject: [PATCH 33/38] 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 Date: Mon, 2 Sep 2024 23:48:13 +1000 Subject: [PATCH 34/38] 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 '' % (self.type.title(), self.code_name) + # return '' % (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 Date: Fri, 28 Feb 2025 11:56:39 +1100 Subject: [PATCH 35/38] 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 Date: Fri, 28 Feb 2025 11:58:09 +1100 Subject: [PATCH 36/38] 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 Date: Fri, 28 Feb 2025 11:59:50 +1100 Subject: [PATCH 37/38] 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 Date: Fri, 28 Feb 2025 12:01:06 +1100 Subject: [PATCH 38/38] 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: """