diff --git a/CHANGELOG.md b/CHANGELOG.md index 53b1bd35..44ecf4d9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,24 +2,92 @@ 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.5] - 2023-06-28 +## [1.31.5] - 2026-05-26 + +- Update optimization function for compatibility with latest `pyswarm` release +- Update minimum Python version to 3.10 as previous versions are deprecated + +## [1.31.4] - 2026-01-28 + +- Updated internal functions for Python 3.14 and Pandas 3.0.0 + +## [1.31.3] - 2025-08-11 + +- Changed initialization method to increase prioritization for characteristics with zero value. If a characteristic has a zero value, this will now automatically force included compartments to have a zero value. In practice, this means that it is less likely that frameworks will cause negative initial popsizes, so some frameworks that previously did not work will now initialize correctly. +- Added `at.TimeSeries.clear()` to reset a time series while preserving the units. This function can be useful when updating the value of databook quantities programatically. +- Added a separate numerical tolerance used for initialization (`at.model.model_settings['initialization_tolerance']`) which permits more approximate initializations while still maintaining the same numerical tolerance for the rest of the integration. +- Added another logging level (`at.VERBOSE`) which enables more targeted additional output + +- *Backwards-compatibility notes* + +- Some initializations might show numerical (e.g., `1e-10`) differences in their values due to the new algorithm. In a small number of cases (depending on the framework), it is possible that the updated initialization method could result in a slightly different initialization. +- The default initialization tolerance is now `1e-3` instead of `1e-6`, some models that previously raised a `BadInitialization` error will now run without error. Users should note that if it's necessary to guarantee an exact initialization, this tolerance should be reduced. + +## [1.31.2] - 2025-06-26 + +- Fix bug in creating databook if framework 'Databook pages' sheet contains multiple code names mapping to the same full name. + +## [1.31.1] - 2025-05-30 + +- Switched from `setup.py` to `pyproject.toml`-based installation. + +## [1.31.0] - 2025-02-18 + +- Framework definitions of compartments, characteristics, and parameters support a new column 'databook default all'. If set to 'y', then when a databook is produced, the data entry table will only contain a record for 'All' instead of having population-specific rows. Further manual editing of the tables is supported as normal. + +## [1.30.0] - 2025-02-17 + +- Automatic calibration can now selectively weight parts of the time series to select or prioritise a subset of time points. +- Added YAML-based calibration support to Atomica, covered in Tutorial 7 in the online documentation. + +## [1.29.0] - 2024-10-30 + +- `ProjectSettings` now computes the simulation time vector in a more robust way to reduce edge cases where the reported `sim_dt` doesn't match the input. +- `ParameterSet.load_calibration()` now clears any existing initialization if the calibration being loaded does not contain an initialization. Previously, the absence of an 'initialization' sheet in the calibration would be treated as not making any change to the initialization. This could cause calibrations to become mixed if a calibration without an initialization was loaded after a calibration with an initialization. Now, a missing initialization sheet is treated as meaning 'no initialization' and any existing initialization will be cleared when the calibration is loaded. +- `Project.load_databook()` will no longer populate `Project.databook` when a `ProjectData` instance is supplied rather than a spreadsheet. The intention of `Project.databook` as opposed to `Project.data.to_spreadsheet()` is that the original databook may contain comments or other content that is not preserved when the databook is loaded into a `ProjectData` instance. Therefore, `Project.databook` serves as a record of the original inputs. However, in previous versions of Atomica, if a `ProjectData` instance was provided rather than a spreadsheet, `ProjectData.to_spreadsheet()` would be used to populate `Project.databook`. For large databooks, this can be computationally expensive and particularly affect the use case of passing in a preloaded databook to improve performance. Since the conversion of the `ProjectData` to a spreadsheet upon loading offers no functional difference to creating the spreadsheet from `Project.data` when required, Atomica no longer performs this conversion upfront. +- `PlotData.time_aggregate()` now only interpolates if necessary, using simulated time points as much as possible. + +*Backwards-compatibility notes* + +- In some edge cases, the simulation time points in the output may be different. In those cases, the difference between simulation time points in the model output would not have matched the model input, although the correct time step would have been used to calculate parameter values. In these cases, there may be an extra time point in the model output. Re-running the model should produce results that are close to the original results. +- If accessing `Project.databook`, in some cases this may now be `None` rather than an `sc.Spreadsheet()`. If that occurs, `Project.data.to_spreadsheet()` should be used to produce an equivalent spreadsheet. +- Time aggregation of `PlotData` may produce slightly different results due to the more accurate selection of time points in this version. +- +## [1.28.7] - 2024-10-29 + +- If aggregating characteristics with a denominator, weighted aggregations use the denominator of the quantity rather than the total population size to perform the weighting. This can be useful for quantities that are proportions of things other than the population size e.g., proportion of active infections that are diagnosed +- If aggregating characteristics with a denominator, weighted aggregation will be used by default (rather than 'average') +- Output aggregation of durations now uses 'sum' instead of 'average' by default +- If no aggregation method is specified, the aggregation method will now be selected separately for each output +- Population aggregation of non-number units now uses 'weighted' by default rather than 'average' + +*Backwards-compatibility notes* + +- Results obtained when aggregating characteristics with a denominator using the 'weighted' method will change. To reproduce the previous results, it is necessary to perform the population size weighting manually by removing the aggregation, then extracting the population sizes, and using those to aggregate the outputs. This is considered to be a rare use case because the updated result is a more useful weighting compared to the previous result. +- Results obtained when aggregating characteristics with a denominator _without_ explicitly specifying a method will change, because 'weighted' aggregation is now used by default. + +## [1.28.6] - 2024-07-12 + +- Support entering `'total'` as the population name in auto-calibration measurables to calibrate aggregated values across populations in the model to aggregate values entered in the databook under a 'Total' population + +## [1.28.5] - 2024-06-28 - Enable automated calibration of transfers and updated documentation to cover this feature -## [1.28.4] - 2023-05-27 +## [1.28.4] - 2024-05-27 - Added an option to save initial compartment sizes inside a `ParameterSet`. Importantly, this saved representation allows setting the initial _subcompartment_ sizes for a `TimedCompartment`. It therefore offers the possibility of initializing the model in a steady state computed from a previous simulation run, that would not be possible to initialize conventionally because standard initialization uniformly distributes people into the subcompartments of a timed compartment. - Added `ParameterSet.make_constant` to facilitate constructing copies of `ParameterSet` instances that are constant over time. -## [1.28.3] - 2023-05-16 +## [1.28.3] - 2024-05-16 - Updated `at.Project()` to explicitly take in the settings arguments for `sim_start`, `sim_end`, and `sim_dt`. These are now applied after databooks are loaded, fixing a bug where these arguments would get overwritten when loading the databook. -## [1.28.2] - 2023-04-05 +## [1.28.2] - 2024-04-05 - `at.calibrate` now supports passing any additional arguments into the optimization function e.g., `sc.asd` allowing additional options for customizing the optimization. -## [1.28.1] - 2023-02-05 +## [1.28.1] - 2024-02-05 - Updated various Pandas operations to improve compatibility with Pandas 2.2.0 - Replaced 'probability' units with 'rate' units in many of the library example frameworks diff --git a/LICENSE b/LICENSE index 43382f75..c75ac9ae 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2020 Atomica Development Team +Copyright (c) 2017 - 2025 Atomica Development Team Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index 2316c6cb..59573fd3 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![Build Status](https://dev.azure.com/AtomicaTeam/Atomica/_apis/build/status/atomicateam.atomica?branchName=master)](https://dev.azure.com/AtomicaTeam/Atomica/_build/latest?definitionId=1&branchName=master) -[![PyPi version](https://badgen.net/pypi/v/atomica/)](https://pypi.com/project/atomica) +[![PyPi version](https://badgen.net/pypi/v/atomica/)](https://pypi.org/project/atomica) Atomica is a simulation engine for compartmental models. It can be used to simulate disease epidemics, health care cascades, and many other things. diff --git a/atomica/__init__.py b/atomica/__init__.py index fa6a254f..efcf58d0 100644 --- a/atomica/__init__.py +++ b/atomica/__init__.py @@ -26,7 +26,10 @@ import sys import logging +VERBOSE = 15 +logging.addLevelName(VERBOSE, "VERBOSE") logger = logging.getLogger("atomica") +logging.Logger.verbose = lambda self, msg, *args, **kwargs: self.log(VERBOSE, msg, *args, **kwargs) if not logger.handlers: # Only add handlers if they don't already exist in the module-level logger @@ -34,16 +37,16 @@ # prior to importing Atomica, and the user's custom logger won't be overwritten as long as it has # at least one handler already added. The use case was originally to suppress messages on import, but since # importing is silent now, it doesn't matter so much. - debug_handler = logging.StreamHandler(sys.stdout) # info_handler will handle all messages below WARNING sending them to STDOUT - info_handler = logging.StreamHandler(sys.stdout) # info_handler will handle all messages below WARNING sending them to STDOUT + debug_handler = logging.StreamHandler(sys.stdout) # handle all messages below WARNING sending them to STDOUT + info_handler = logging.StreamHandler(sys.stdout) # handle all messages below WARNING sending them to STDOUT warning_handler = logging.StreamHandler(sys.stderr) # warning_handler will send all messages at or above WARNING to STDERR - debug_handler.setLevel(0) # Handle all lower levels - the output should be filtered further by setting the logger level, not the handler level - info_handler.setLevel(logging.INFO) # Handle all lower levels - the output should be filtered further by setting the logger level, not the handler level + debug_handler.setLevel(0) # Handle all levels - the output is then filtered further by setting the logger level, not the handler level + info_handler.setLevel(logging.INFO) # Handle levels INFO or higher warning_handler.setLevel(logging.WARNING) - debug_handler.addFilter(type("ThresholdFilter", (object,), {"filter": lambda x, logRecord: logRecord.levelno < logging.INFO})()) # Display anything INFO or higher - info_handler.addFilter(type("ThresholdFilter", (object,), {"filter": lambda x, logRecord: logRecord.levelno < logging.WARNING})()) # Don't display WARNING or higher + debug_handler.addFilter(type("ThresholdFilter", (object,), {"filter": lambda x, logRecord: logRecord.levelno < logging.INFO})()) # Display only messages below the INFO level + info_handler.addFilter(type("ThresholdFilter", (object,), {"filter": lambda x, logRecord: logRecord.levelno < logging.WARNING})()) # Display only messages below WARNING (the info_handler level already rejects anything below INFO) warning_formatter = logging.Formatter("%(levelname)s {%(filename)s:%(lineno)d} - %(message)s") warning_handler.setFormatter(warning_formatter) @@ -56,12 +59,16 @@ logger.addHandler(warning_handler) logger.setLevel("INFO") # Set the overall log level -from .version import version as __version__, versiondate as __versiondate__, gitinfo as __gitinfo__ -# Check scipy version -import scipy +# Get version information +from .version import version as __version__, versiondate as __versiondate__ import sciris as sc +__gitinfo__ = sc.gitinfo(__file__) +version.gitinfo = __gitinfo__ # Add back (so version.py does not require imports) + +# Check scipy version +import scipy if sc.compareversions(scipy.__version__, "1.2.1") < 0: raise Exception(f"Atomica requires Scipy 1.2.1 or later - installed version is {scipy.__version__}") @@ -95,3 +102,7 @@ from .scenarios import * from .system import * from .utils import * + +import pathlib + +rootdir = pathlib.Path(__file__).parent.parent diff --git a/atomica/calibration.py b/atomica/calibration.py index 16046035..c100a8b5 100644 --- a/atomica/calibration.py +++ b/atomica/calibration.py @@ -12,6 +12,7 @@ from .system import logger from .parameters import ParameterSet import logging +import atomica __all__ = ["calibrate"] @@ -20,35 +21,38 @@ calibration_settings["tolerance"] = 1e-6 -def _update_parset(parset, y_factors, pars_to_adjust): - # Insert updated y-values into the parset - # - parset : a ParameterSet object - # - y_factors : Array with as many elements as pars_to_adjust - # - pars_to_adjust : Array of tuples (par_name,pop_name,...) with special value pop='all' supported to set meta factor - # Must have as many elements as y_factors. pop=None is not allowed - it must be converted - # to a full list of pops previously (in perform_autofit) +def _update_parset(parset, y_factors, pars_to_adjust) -> None: + """ - for i, x in enumerate(pars_to_adjust): - par_name = x[0] - pop_name = x[1] + :param parset: a ParameterSet object + :param y_factors: Array with as many elements as pars_to_adjust + :param pars_to_adjust: Array of tuples (par_name,pop_name,...) with special value pop='all' supported to set meta factor + Must have as many elements as y_factors. pop=None is not allowed - it must be converted + to a full list of pops previously (in perform_autofit) + :return: None (parset is modified in-place) + """ - if par_name in parset.pars: - if pop_name == "all": - par = parset.pars[par_name] - par.meta_y_factor = y_factors[i] - else: - parset.pars[par_name].y_factor[pop_name] = y_factors[i] + for i, (par_name, pop_name, *_) in enumerate(pars_to_adjust): + par = parset.get_par(par_name) + if pop_name.lower() == "all": + par.meta_y_factor = y_factors[i] else: - # Handle transfers here - tokens = par_name.split("_from_") - par = parset.transfers[tokens[0]][tokens[1]] par.y_factor[pop_name] = y_factors[i] -def _calculate_objective(y_factors, pars_to_adjust, output_quantities, parset, project): - # y-factors, array of y-factors to apply to specified output_quantities - # pars_to_adjust - list of tuples (par_name,pop_name,...) recognized by parset.update() - # output_quantities - a tuple like (pop,var,weight,metric) understood by model.get_pop[pop].getVar +def _calculate_objective(y_factors, pars_to_adjust, output_quantities, parset, project, *args, **kwargs) -> float: + """ + Run the model for a given set of y-factors and return the objective/goodness-of-fit + + Additional extra arguments will be ignored but will not raise an error. + + :param y_factors: array of y-factors to apply to specified output_quantities + :param pars_to_adjust: list of tuples (par_name,pop_name,...) recognized by parset.update() + :param output_quantities: a tuple containing (pop, var, weight, metric, start_year, end_year) - start year and end year are inclusive + :param parset: + :param project: + :return: The value of the objective function defined by the output_quantities + """ _update_parset(parset, y_factors, pars_to_adjust) @@ -59,20 +63,36 @@ def _calculate_objective(y_factors, pars_to_adjust, output_quantities, parset, p objective = 0.0 - for var_label, pop_name, weight, metric in output_quantities: + for var_label, pop_name, weight, metric, start_year, end_year in output_quantities: + target = project.data.get_ts(var_label, pop_name) # This is the TimeSeries with the data for the requested quantity if target is None: continue if not target.has_time_data: # Only use this output quantity if the user entered time-specific data continue - var = result.model.get_pop(pop_name).get_variable(var_label) + + if pop_name.lower() == "total": + var = atomica.PlotData(result, outputs=var_label, pops="total", project=project) + else: + var = result.model.get_pop(pop_name).get_variable(var_label) + data_t, data_v = target.get_arrays() + inds = (data_t >= start_year) & (data_t <= end_year) + if np.count_nonzero(inds) == 0: + # If no time points remain after filtering down to the time points the user requested + logger.info(f"No data points remaining after filtering down to requested time period. Skipping...") + continue + data_t = data_t[inds] + data_v = data_v[inds] # Interpolate the model outputs onto the data times # If there is data outside the range when the model was simulated, don't # extrapolate the model outputs y = data_v - y2 = np.interp(data_t, var[0].t, var[0].vals, left=np.nan, right=np.nan) + if pop_name.lower() == "total": + y2 = np.interp(data_t, var.series[0].tvec, var.series[0].vals, left=np.nan, right=np.nan) + else: + y2 = np.interp(data_t, var[0].t, var[0].vals, left=np.nan, right=np.nan) idx = ~np.isnan(y) & ~np.isnan(y2) objective += weight * sum(_calculate_fitscore(y[idx], y2[idx], metric)) @@ -113,13 +133,13 @@ def _calc_wape(y_obs, y_fit): return abs(y_fit - y_obs) / (y_obs.mean() + calibration_settings["tolerance"]) -def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, max_time=60, method="asd", **kwargs) -> ParameterSet: +def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, max_time=None, method="asd", time_period=(-np.inf, np.inf), **kwargs) -> ParameterSet: """ Run automated calibration :param project: A project instance to provide data and sim settings :param parset: A :class:`ParameterSet` instance to calibrate - :param pars_to_adjust: list of tuples, (par_name,pop_name,lower_limit,upper_limit) + :param pars_to_adjust: list of tuples, `(par_name, pop_name, lower_bound, upper_bound, initial_value)` 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. To calibrate a transfer, the parameter name should be set to @@ -127,8 +147,18 @@ def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, 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 + function. pop_name=None will expand to all pops. pop_name='all' is not supported. The + output can optionally contain `(var_label, pop_name, weight, metric, start_year, end_year)` + to select a subset of the data for evaluating the objective. The start year and end year + specified here will take precedence over the time_period argument. In some cases, it may be desirable + to fit to an aggregated total value across populations. In that case, the databook should have + an extra row in the TDVE table for a population called "Total". The measurable + can then be given pop_name="Total" which will cause the Atomica model outputs to be aggregated over + all populations, and the aggregate value compared to the "Total" data. The aggregation methods will + be automatically selected depending on units of the quantity (sum for "number" units, average for others). :param max_time: If using ASD, the maximum run time + :param time_period: Tuple of start and end years to use for the objective function. Applies to all outputs unless + the output has an explicitly specified start and end year :param method: 'asd' or 'pso'. If using 'pso' all upper and lower limits must be finite :param kwargs: Dictionary of additional arguments to be passed to the optimization function, e.g. stepsize or pinitial :return: A calibrated :class:`ParameterSet` @@ -138,95 +168,134 @@ def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, # Expand out pop=None in pars_to_adjust p2 = [] for par_tuple in pars_to_adjust: + if len(par_tuple) == 5: + initial_value = par_tuple[4] + else: + initial_value = None + if par_tuple[1] is None: # If the pop name is None - par = parset.pars[par_tuple[0]] + par = parset.get_par(par_tuple[0]) for pop_name in par.pops: - p2.append((par_tuple[0], pop_name, par_tuple[2], par_tuple[3])) + p2.append((par_tuple[0], pop_name, par_tuple[2], par_tuple[3], initial_value)) else: - p2.append(par_tuple) + p2.append((par_tuple[0], par_tuple[1], par_tuple[2], par_tuple[3], initial_value)) + pars_to_adjust = p2 # Expand out pop=None in output_quantities - o2 = [] + outputs = [] for output_tuple in output_quantities: - if output_tuple[1] is None: # If the pop name is None - pops = project.data.pops.keys() - for pop_name in pops: - o2.append((output_tuple[0], pop_name, output_tuple[2], output_tuple[3])) + var_label = output_tuple[0] + pop_name = output_tuple[1] + weight = output_tuple[2] + metric = output_tuple[3] + if len(output_tuple) == 6: + start_year = output_tuple[4] + end_year = output_tuple[5] else: - o2.append(output_tuple) - output_quantities = o2 + start_year = time_period[0] + end_year = time_period[1] - args = { - "project": project, - "parset": parset.copy(), - "pars_to_adjust": pars_to_adjust, - "output_quantities": output_quantities, - } + if pop_name is None: # If the pop name is None + for pop in project.data.pops.keys(): + outputs.append((var_label, pop, weight, metric, start_year, end_year)) + else: + outputs.append((var_label, pop_name, weight, metric, start_year, end_year)) + + output_quantities = outputs x0 = [] xmin = [] xmax = [] + filtered_pars_to_adjust = [] + parset = parset.copy() + for i, x in enumerate(pars_to_adjust): - par_name, pop_name, scale_min, scale_max = x - if par_name in parset.pars: - par = parset.pars[par_name] + par_name, pop_name, scale_min, scale_max, initial_value = x + + par = parset.get_par(par_name) + + # if initial_value has not been explicitly set in the tuple: use y_factor in parset + if initial_value is None: if pop_name == "all": - x0.append(par.meta_y_factor) + initial_value = par.meta_y_factor else: - x0.append(par.y_factor[pop_name]) + initial_value = par.y_factor[pop_name] + # if this value is outside the min and max bounds, make it equal to min or max (whichever is closest) + # if min == max, this should make the initial value equal to min and max + initial_value = np.clip(initial_value, scale_min, scale_max) else: - tokens = par_name.split("_from_") - par = parset.transfers[tokens[0]][tokens[1]] - x0.append(par.y_factor[pop_name]) - xmin.append(scale_min) - xmax.append(scale_max) + assert (initial_value >= scale_min) and (initial_value <= scale_max), "Initial value is not consistent with the lower/upper bounds" - original_sim_end = project.settings.sim_end - project.settings.sim_end = min(project.data.tvec[-1], original_sim_end) + # update y_factors in parset + if pop_name == "all": + par.meta_y_factor = initial_value + else: + par.y_factor[pop_name] = initial_value - try: - if method == "asd": - optim_args = { - "stepsize": 0.1, - "maxiters": 2000, - "sinc": 1.5, - "sdec": 2.0, - "reltol": 1e-3, - "abstol": 1e-6, - "xmin": xmin, - "xmax": xmax, - } - optim_args.update(kwargs) - - if max_time is not None: - optim_args["maxtime"] = max_time - - log_level = logger.getEffectiveLevel() - if log_level < logging.WARNING: - optim_args["verbose"] = 2 - else: - optim_args["verbose"] = 0 + if scale_min != scale_max: + # Only include the y-factor in the adjustments made in the optimization function if a range + # of y-factor values is permitted + filtered_pars_to_adjust.append(x) + x0.append(initial_value) + xmin.append(scale_min) + xmax.append(scale_max) - opt_result = sc.asd(_calculate_objective, x0, args, **optim_args) - x1 = opt_result["x"] - elif method == "pso": - import pyswarm + args = { + "project": project, + "parset": parset, + "pars_to_adjust": filtered_pars_to_adjust, + "output_quantities": output_quantities, + } - optim_args = {"maxiter": 3, "lb": xmin, "ub": xmax, "minstep": 1e-3, "debug": True} - if np.any(~np.isfinite(xmin)) or np.any(~np.isfinite(xmax)): - errormsg = "PSO optimization requires finite upper and lower bounds to specify the search domain (i.e. every parameter being adjusted needs to have finite bounds)" - raise Exception(errormsg) + original_sim_end = project.settings.sim_end + project.settings.sim_end = min(project.data.tvec[-1], original_sim_end) - x1, _ = pyswarm.pso(_calculate_objective, kwargs=args, **optim_args) - else: - raise Exception("Unrecognized method") - except Exception as e: - raise e - finally: - project.settings.sim_end = original_sim_end # Restore the simulation end year + if len(filtered_pars_to_adjust) > 0: + try: + if method == "asd": + optim_args = { + "stepsize": 0.1, + "maxiters": 2000, + "sinc": 1.5, + "sdec": 2.0, + "reltol": 1e-3, + "abstol": 1e-6, + "xmin": xmin, + "xmax": xmax, + "maxtime": 60 if max_time is None else max_time, + "minval": 0, + } + optim_args.update(kwargs) + + if "verbose" not in optim_args: + log_level = logger.getEffectiveLevel() + if log_level < logging.WARNING: + optim_args["verbose"] = 2 + else: + optim_args["verbose"] = 0 + + opt_result = sc.asd(_calculate_objective, x0, args, **optim_args) + x1 = opt_result["x"] + elif method == "pso": + import pyswarm + + optim_args = {"maxiter": 3, "lb": xmin, "ub": xmax, "minstep": 1e-3, "debug": True} + if np.any(~np.isfinite(xmin)) or np.any(~np.isfinite(xmax)): + errormsg = "PSO optimization requires finite upper and lower bounds to specify the search domain (i.e. every parameter being adjusted needs to have finite bounds)" + raise Exception(errormsg) + + x1, _ = pyswarm.pso(_calculate_objective, kwargs=args, **optim_args) + else: + raise Exception("Unrecognized method") + except Exception as e: + raise e + finally: + project.settings.sim_end = original_sim_end # Restore the simulation end year - _update_parset(args["parset"], x1, pars_to_adjust) + _update_parset(args["parset"], x1, args["pars_to_adjust"]) + else: + logger.info("No parameters to adjust provided to the optimisation function. Skipping optimisation...") # Log out the commands required for equivalent manual calibration if desired for i, x in enumerate(pars_to_adjust): @@ -236,7 +305,7 @@ def calibrate(project, parset: ParameterSet, pars_to_adjust, output_quantities, if par_name in parset.pars: par = args["parset"].pars[par_name] - if pop_name == "all": + if pop_name.lower() == "all": logger.debug("parset.get_par('{}').meta_y_factor={:.2f}".format(par_name, par.meta_y_factor)) else: logger.debug("parset.get_par('{}').y_factor['{}']={:.2f}".format(par_name, pop_name, par.y_factor[pop_name])) diff --git a/atomica/data.py b/atomica/data.py index 4422e715..2ecf0713 100644 --- a/atomica/data.py +++ b/atomica/data.py @@ -21,8 +21,7 @@ from collections import defaultdict import pandas as pd import itertools - -_DEFAULT_PROVENANCE = "Framework-supplied default" +from .version import version, gitinfo __all__ = ["InvalidDatabook", "ProjectData"] @@ -57,14 +56,24 @@ def __init__(self, framework): self.interpops = list() #: This stores a list of :class:`TimeDependentConnections` instances for interactions self.tvec = None #: This is the data's tvec used when instantiating new tables. Not _guaranteed_ to be the same for every TDVE/TDC table self.tdve = sc.odict() #: This is an odict storing :class:`TimeDependentValuesEntry` instances keyed by the code name of the TDVE - self.tdve_pages = sc.odict() #: This is an odict mapping worksheet name to an (ordered) list of TDVE code names appearing on that sheet + self.tdve_pages = sc.odict(defaultdict=list) #: This is an odict mapping worksheet name to an (ordered) list of TDVE code names appearing on that sheet # Internal storage used with methods while writing - self._pop_types = list(framework.pop_types.keys()) # : Store set of valid population types from framework + self._pop_types = list(framework.pop_types.keys()) #: Store set of valid population types from framework self._formats = None #: Temporary storage for the Excel formatting while writing a databook self._book = None #: Temporary storage for the workbook while writing a databook self._references = None #: Temporary storage for cell references while writing a databook + self.version = version #: Current Atomica version + self.gitinfo = sc.dcp(gitinfo) #: Atomica Git version information, if being run in a Git repository + + def __setstate__(self, d): + from .migration import migrate + + self.__dict__ = d + projectdata = migrate(self) + self.__dict__ = projectdata.__dict__ + def tables(self): """ Return iterator over all TDVE and TDC tables @@ -161,7 +170,12 @@ def get_ts(self, name: str, key=None): is therefore encountered on the :class:`Result` and plotting side. If retrieving values for a comp/charac/par and the databook contains an entry for 'all' rather - than specific populations, then the 'all' time series will be returned regardless of the key + than specific populations, then the 'all' time series will be returned if the key does not match + any population names. + + Note that the TDVE can contain time series for population names that don't correspond to actual populations. + For example, users could add an extra row for 'Total', or anything else. `get_ts()` will retrieve these + rows without an error, so population names are not validated against the populations in the databook. """ # Exit immediately if the name is not specified @@ -255,7 +269,7 @@ class instance (e.g. if creating a new databook). data = ProjectData(framework=framework) data.tvec = sc.promotetoarray(tvec) pages = defaultdict(list) # This will store {sheet_name:(code_name,databook_order)} which will then get sorted further - + page_names = framework.sheets["databook pages"][0].set_index('datasheet code name')['datasheet title'].to_dict() for obj_type, df in zip(["comps", "characs", "pars"], [framework.comps, framework.characs, framework.pars]): for _, spec in df.iterrows(): databook_page = spec.get("databook page") @@ -263,13 +277,15 @@ class instance (e.g. if creating a new databook). pop_type = spec.get("population type") databook_order = spec.get("databook order") full_name = spec["display name"] + default_all = spec["databook default all"] == "y" + allowed_units = [framework.get_databook_units(full_name)] if pd.isna(databook_order): order = np.inf else: order = databook_order - pages[databook_page].append((spec.name, order)) - data.tdve[spec.name] = TimeDependentValuesEntry(full_name, data.tvec, allowed_units=[framework.get_databook_units(full_name)], comment=spec["guidance"], pop_type=pop_type) + pages[page_names[databook_page]].append((spec.name, order)) + data.tdve[spec.name] = TimeDependentValuesEntry(full_name, data.tvec, allowed_units=allowed_units, comment=spec["guidance"], pop_type=pop_type, default_all=default_all) data.tdve[spec.name].write_units = True data.tdve[spec.name].write_uncertainty = True if obj_type == "pars": @@ -279,14 +295,14 @@ class instance (e.g. if creating a new databook). data.tdve[spec.name].write_uncertainty = False # Don't show uncertainty for timed parameters. In theory users could manually add the column and sample over it, but because the duration is rounded to the timestep, it's likely to have confusing stepped effects data.tdve[spec.name].pop_type = pop_type - # Now convert pages to full names and sort them into the correct order - for _, spec in framework.sheets["databook pages"][0].iterrows(): + if default_all: + # add_pop normally adds TDVE rows, but it won't operate on any TDVEs that default to 'All' so we need to add the 'All' rows here + data.tdve[spec.name].ts["All"] = TimeSeries(units=allowed_units[0]) - if spec["datasheet code name"] in pages: - pages[spec["datasheet code name"]].sort(key=lambda x: x[1]) - data.tdve_pages[spec["datasheet title"]] = [x[0] for x in pages[spec["datasheet code name"]]] - else: - data.tdve_pages[spec["datasheet title"]] = list() + # Now sort them into the correct order + for page, tables in pages.items(): + tables.sort(key=lambda x: x[1]) + data.tdve_pages[page] = [x[0] for x in tables] # Now, proceed to add pops, transfers, and interactions for code_name, spec in new_pops.items(): @@ -303,7 +319,7 @@ class instance (e.g. if creating a new databook). ts = TimeSeries(units=interpop.allowed_units[0]) ts.insert(None, spec["default value"]) interpop.ts[(from_pop, to_pop)] = ts - interpop.ts_attributes["Provenance"][(from_pop, to_pop)] = _DEFAULT_PROVENANCE + interpop.ts_attributes["Provenance"][(from_pop, to_pop)] = spec["provenance"] if "provenance" in spec else FS.DEFAULT_PROVENANCE # Finally, insert parameter and characteristic default values for df in [framework.comps, framework.characs, framework.pars]: @@ -315,7 +331,7 @@ class instance (e.g. if creating a new databook). tdve = data.tdve[spec.name] for key, ts in tdve.ts.items(): ts.insert(None, spec["default value"]) - tdve.ts_attributes["Provenance"][key] = _DEFAULT_PROVENANCE + tdve.ts_attributes["Provenance"][key] = spec["provenance"] if "provenance" in spec else FS.DEFAULT_PROVENANCE return data @@ -413,7 +429,9 @@ def from_spreadsheet(spreadsheet, framework): ts.units = tdve.allowed_units[0] if not spec["databook page"]: - logger.warning('A TDVE table for "%s" (%s) was read in and will be used, but the Framework did not mark this quantity as appearing in the databook', tdve.name, code_name) + # Note that if the parameter doesn't have a databook page and the framework is valid, then the parameter must have a function. Therefore, + # if data is also read in, it will not change the simulation outputs and would only be used for calibration/validation + logger.warning('A TDVE table for "%s" (%s) was read in and data will be available for calibration, but the Framework did not mark this quantity as appearing in the databook', tdve.name, code_name) tdve.comment = spec["guidance"] if code_name in self.tdve: @@ -489,8 +507,10 @@ def validate(self, framework) -> bool: for pop in self.pops.keys(): self.tdve[spec_name].ts[pop] = TimeSeries(assumption=spec["default value"], units=units) tdve_page = framework.sheets["databook pages"][0][framework.sheets["databook pages"][0]["datasheet code name"] == spec["databook page"]]["datasheet title"].values[0] - if tdve_page in self.tdve_pages: - self.tdve_pages[tdve_page].append(spec_name) + for existing in self.tdve_pages.keys(): + if existing.lower() == tdve_page.lower(): + self.tdve_pages[existing].append(spec_name) + break else: self.tdve_pages[tdve_page] = [spec_name] else: @@ -625,7 +645,7 @@ def add_pop(self, code_name: str, full_name: str, pop_type: str = None) -> None: :param code_name: The code name for the new population :param full_name: The full name/label for the new population - :param pop_type: String with the population type code name + :param pop_type: String with the population type code name (optional) - default is the type of the first population """ @@ -651,8 +671,11 @@ def add_pop(self, code_name: str, full_name: str, pop_type: str = None) -> None: for tdve in self.tdve.values(): # Since TDVEs in databooks must have the unit set in the framework, all ts objects must share the same units # And, there is only supposed to be one type of unit allowed for TDVE tables (if the unit is empty, it will be 'N.A.') - # so can just pick the first of the allowed units - if tdve.pop_type == pop_type: + # so can just pick the first of the allowed units. We will add the population row if the pop type matches and if + # the TDVE is either not a 'default_all' or if the user has removed the 'All' row from the TDVE despite it being default_all + if tdve.pop_type != pop_type or (tdve.default_all and ("All" in tdve.ts or "all" in tdve.ts)): + continue + else: tdve.ts[code_name] = TimeSeries(units=tdve.allowed_units[0]) def rename_pop(self, existing_code_name: str, new_code_name: str, new_full_name: str) -> None: diff --git a/atomica/excel.py b/atomica/excel.py index 2a0119ce..c40f46dc 100644 --- a/atomica/excel.py +++ b/atomica/excel.py @@ -403,7 +403,7 @@ def __init__(self, code_name: str, full_name: str, tvec: np.array, from_pops: li self.ts = ts if ts is not None else sc.odict() self.attributes = {} #: Attributes associated with the table - self.ts_attributes = {} #: Attributes associated with each TimeSeries row + self.ts_attributes = {} #: Attributes associated with each TimeSeries row. This is a dict-of-dicts with `ts_attributes[attribute]={(from_pop, to_pop):value}` e.g., `ts_attributes["Provenance"] = {('0-4', '5-14'): 'value'}` self.ts_attributes["Provenance"] = {} # Include provenance attribute by default self.assumption_heading = "Constant" #: Heading to use for assumption column @@ -884,7 +884,7 @@ class TimeDependentValuesEntry: """ - def __init__(self, name, tvec: np.array = None, ts=None, allowed_units: list = None, comment: str = None, pop_type: str = None): + def __init__(self, name, tvec: np.array = None, ts=None, allowed_units: list = None, comment: str = None, pop_type: str = None, default_all: bool = False): if ts is None: ts = sc.odict() @@ -907,6 +907,8 @@ def __init__(self, name, tvec: np.array = None, ts=None, allowed_units: list = N self.write_uncertainty = None #: Write a column for uncertainty (if None, uncertainty will be written if any of the TimeSeries have uncertainty) self.write_assumption = None #: Write a column for assumption/constant (if None, assumption will be written if any of the TimeSeries have an assumption) + self.default_all = default_all #: Record whether the framework specifies that this TDVE should default to having an 'All' row instead of population-specific rows (the user can manually modify further) + def __repr__(self): output = sc.prepr(self) return output diff --git a/atomica/framework.py b/atomica/framework.py index 7e686f49..cd378003 100644 --- a/atomica/framework.py +++ b/atomica/framework.py @@ -627,12 +627,24 @@ def _sanitize_compartments(self) -> None: self.sheets["compartments"] = [pd.DataFrame(columns=["code name", "display name"])] required_columns = ["display name"] - defaults = {"is sink": "n", "is source": "n", "is junction": "n", "databook page": None, "default value": None, "databook order": None, "guidance": None, "population type": None} # Default is for it to be randomly ordered if the databook page is not None + defaults = { + "is sink": "n", + "is source": "n", + "is junction": "n", + "databook page": None, + "databook default all": "n", + "default value": None, + "population type": None, + "databook order": None, # Default is for it to be randomly ordered if the databook page is not None + "guidance": None, + "provenance": FS.DEFAULT_PROVENANCE, + } valid_content = { "display name": None, # Valid content being `None` means that it just cannot be empty "is sink": {"y", "n"}, "is source": {"y", "n"}, "is junction": {"y", "n"}, + "databook default all": {"y", "n"}, } numeric_columns = ["databook order", "default value"] @@ -655,9 +667,9 @@ def _sanitize_compartments(self) -> None: self.comps["setup weight"] = ((~self.comps["databook page"].isna() | ~self.comps["default value"].isna()) & (self.comps["is source"] == "n") & (self.comps["is sink"] == "n")).astype(float) else: fill_ones = self.comps["setup weight"].isna() & (~self.comps["databook page"].isna() | ~self.comps["default value"].isna()) & (self.comps["is source"] == "n") & (self.comps["is sink"] == "n") + self.comps["setup weight"] = self.comps["setup weight"].astype(float) self.comps.loc[fill_ones, "setup weight"] = 1 self.comps.loc[self.comps["setup weight"].isna(), "setup weight"] = 0 - self.comps["setup weight"] = self.comps["setup weight"].astype(float) if "calibrate" not in self.comps: # If calibration column is not present, then it calibrate if in the databook @@ -707,10 +719,21 @@ def _sanitize_characteristics(self) -> None: self.sheets["characteristics"] = [pd.DataFrame(columns=["code name", "display name"])] required_columns = ["display name"] - defaults = {"components": None, "denominator": None, "default value": None, "databook page": None, "databook order": None, "guidance": None, "population type": None} + defaults = { + "components": None, + "denominator": None, + "default value": None, + "databook page": None, + "databook order": None, + "databook default all": "n", + "guidance": None, + "population type": None, + "provenance": FS.DEFAULT_PROVENANCE, + } valid_content = { "display name": None, "components": None, + "databook default all": {"y", "n"}, } numeric_columns = ["databook order", "default value"] @@ -727,9 +750,9 @@ def _sanitize_characteristics(self) -> None: self.characs["setup weight"] = (~self.characs["databook page"].isna() | ~self.characs["default value"].isna()).astype(float) else: fill_ones = self.characs["setup weight"].isna() & (~self.characs["databook page"].isna() | ~self.characs["default value"].isna()) + self.characs["setup weight"] = self.characs["setup weight"].astype(float) self.characs.loc[fill_ones, "setup weight"] = 1 self.characs.loc[self.characs["setup weight"].isna(), "setup weight"] = 0 - self.characs["setup weight"] = self.characs["setup weight"].astype(float) if "calibrate" not in self.characs: # If calibration column is not present, then it calibrate if in the databook @@ -801,7 +824,12 @@ def _sanitize_interactions(self) -> None: self.sheets["interactions"] = [pd.DataFrame(columns=["code name", "display name", "to population type", "from population type"])] required_columns = ["display name"] - defaults = {"default value": None, "from population type": None, "to population type": None} + defaults = { + "default value": None, + "from population type": None, + "to population type": None, + "provenance": FS.DEFAULT_PROVENANCE, + } valid_content = { "display name": None, } @@ -837,18 +865,21 @@ def _sanitize_parameters(self) -> None: "function": None, "databook page": None, "databook order": None, + "databook default all": "n", "targetable": "n", "guidance": None, "timescale": None, "population type": None, "is derivative": "n", "timed": "n", + "provenance": FS.DEFAULT_PROVENANCE, } valid_content = { "display name": None, "targetable": {"y", "n"}, "is derivative": {"y", "n"}, "timed": {"y", "n"}, + "databook default all": {"y", "n"}, } numeric_columns = ["databook order", "default value", "minimum value", "maximum value", "timescale"] @@ -872,7 +903,7 @@ def _sanitize_parameters(self) -> None: # If framework has units that case-insensitively match the standard units, then correct the case if not self.pars.empty: lower_idx = self.pars["format"].str.lower().isin(FS.STANDARD_UNITS) - self.pars["format"][lower_idx] = self.pars["format"][lower_idx].str.lower() + self.pars.loc[lower_idx, "format"] = self.pars.loc[lower_idx, "format"].str.lower() def _process_transitions(self) -> None: """ diff --git a/atomica/library/framework_template_advanced.xlsx b/atomica/library/framework_template_advanced.xlsx index 639bfb54..73cc4298 100644 Binary files a/atomica/library/framework_template_advanced.xlsx and b/atomica/library/framework_template_advanced.xlsx differ diff --git a/atomica/migration.py b/atomica/migration.py index 7c893a48..724855a7 100644 --- a/atomica/migration.py +++ b/atomica/migration.py @@ -13,7 +13,7 @@ import sys import io -from distutils.version import LooseVersion +from packaging.version import Version from .system import logger from .version import version, gitinfo import sciris as sc @@ -198,7 +198,7 @@ def migrate(obj, registry=migrations, version=version, gitinfo=gitinfo): print("Skipping migration") return obj # If migration is disabled then don't make any changes EXCEPT to add in version and gitinfo which may otherwise be hard to catch - migrations_to_run = sorted(registry[type(obj).__name__], key=lambda m: LooseVersion(m.original_version)) + migrations_to_run = sorted(registry[type(obj).__name__], key=lambda m: Version(m.original_version)) if sc.compareversions(obj.version, version) >= 0: return obj else: @@ -807,9 +807,49 @@ def _convert_framework_columns(framework): framework._validate() return framework -@migration("ParameterSet", "1.28.3", "1.28.4", "Add initialization atttribute") + +@migration("ParameterSet", "1.28.3", "1.28.4", "Add initialization attribute") def _parset_add_initialization(parset): parset.initialization = None return parset +@migration("ProjectData", "1.30.0", "1.31.0", "Add pop types attribute and default_all") +def _projectdata_add_types_default(D): + + if not hasattr(D, "_pop_types"): + # We have to check for existence because + # Are any population types present? + pop_types = list(x["type"] for x in D.pops.values() if "type" in x) + if not pop_types: + pop_types = [FS.DEFAULT_POP_TYPE] + + for pop, spec in D.pops.items(): + if "type" not in spec: + spec["type"] = pop_types[0] + D._pop_types = list({x["type"]: None for x in D.pops.values()}.keys()) # Using a dictionary here allows for order-preserving unique + + for tdve in D.tdve.values(): + + # Fix the '_add_pop_type' migration which added tdve.type instead of tdve.pop_type + if hasattr(tdve, "type"): + tdve.pop_type = tdve.type + delattr(tdve, "type") + elif not hasattr(tdve, "pop_type"): + tdve.pop_type = D._pop_types[0] # The majority of the time, if the pop_type is missing, the default needs to be added + + # Also add the default_all attribute + tdve.default_all = False + + # Some old TDVE tables have a lowercase 'n.a.' instead of the correct default value + for ts in tdve.ts.values(): + if ts.units == FS.DEFAULT_SYMBOL_INAPPLICABLE.lower(): + ts.units = FS.DEFAULT_SYMBOL_INAPPLICABLE + + # A TDVE should always have some allowed units, if the allowed units have not been populated, then draw + # them from the ts entries. There are some older saved files that may have no allowed units even though the ts + # entries themselves have units. + if not hasattr(tdve, "allowed_units") or not tdve.allowed_units: + tdve.allowed_units = list({x.units: None for x in tdve.ts.values()}.keys()) + + return D diff --git a/atomica/model.py b/atomica/model.py index 173de76d..59b2d145 100644 --- a/atomica/model.py +++ b/atomica/model.py @@ -26,6 +26,7 @@ model_settings = dict() model_settings["tolerance"] = 1e-6 +model_settings["initialization_tolerance"] = 1e-3 __all__ = [ "BadInitialization", @@ -602,7 +603,7 @@ def balance(self, ti: int) -> None: if self.duration_group: link._vals[:, ti] = flow else: - link.vals[ti] = flow + link.vals[ti] = flow[0] def initial_flush(self) -> None: """ @@ -1848,7 +1849,6 @@ def initialize_compartments(self, parset: ParameterSet, framework, t_init: float """ - if not self.comps: # If this population has no compartments, then nothing needs to be done return @@ -1862,34 +1862,66 @@ def initialize_compartments(self, parset: ParameterSet, framework, t_init: float # values for the compartments by solving the set of characteristics simultaneously # Build up the comps and characs containing the setup values in the databook - the `b` in `x=A*b` - characs_to_use = framework.characs.index[(framework.characs["setup weight"] > 0) & (framework.characs["population type"] == self.type)] - comps_to_use = framework.comps.index[(framework.comps["setup weight"] > 0) & (framework.comps["population type"] == self.type)] - b_objs = [self.charac_lookup[x] for x in characs_to_use] + [self.comp_lookup[x] for x in comps_to_use] + characs_to_use = set(framework.characs.index[(framework.characs["setup weight"] > 0) & (framework.characs["population type"] == self.type)]) + comps_to_use = set(framework.comps.index[(framework.comps["setup weight"] > 0) & (framework.comps["population type"] == self.type)]) + + b_vals = {} + b_objs = {} + + # First add all of the compartments that will be used from the databook + zero_compartments = set() # Keep track of compartments that will be initialized as zero + for comp in comps_to_use: + par = parset.pars[comp] + val = par.interpolate(t_init, pop_name=self.name)[0] * par.y_factor[self.name] * par.meta_y_factor + obj = self.comp_lookup[comp] + if val == 0: + zero_compartments.add(obj) + else: + b_objs[comp] = obj + b_vals[comp] = val + + # Then add all of the characteristics + for charac in characs_to_use: + par = parset.pars[charac] + val = par.interpolate(t_init, pop_name=self.name)[0] * par.y_factor[self.name] * par.meta_y_factor + obj = self.charac_lookup[charac] + if val == 0: + # If the characteristic is zero, then all of the compartments included in it must also be zero + for comp in obj.get_included_comps(): + if comp.name in b_vals: + # If a separate databook entry just for this compartment says the compartment should be non-zero, then we have two essentially equally + # direct specifications for the compartment. In that case, we should raise an error as the user has explicitly specified contradictory values + raise BadInitialization(f'Compartment {comp.name} was explicitly specified as having a non-zero value, but characteristic {charac} has a zero value - input data not consistent') + else: + zero_compartments.add(comp) + else: + b_objs[charac] = obj + b_vals[charac] = val + if obj.denominator is not None: + denom_par = parset.pars[obj.denominator.name] + b_vals[charac] *= denom_par.interpolate(t_init, pop_name=self.name)[0] * denom_par.y_factor[self.name] * denom_par.meta_y_factor + # Build up the comps corresponding to the `x` values in `x=A*b` i.e. the compartments being solved for - comps = [c for c in self.comps if not (isinstance(c, SourceCompartment) or isinstance(c, SinkCompartment))] - charac_indices = {c.name: i for i, c in enumerate(b_objs)} # Make lookup dict for characteristic indices + comps = [c for c in self.comps if not (isinstance(c, SourceCompartment) or isinstance(c, SinkCompartment) or c in zero_compartments)] + b_indices = {c.name: i for i, c in enumerate(b_objs.values())} # Make lookup dict for characteristic indices comp_indices = {c.name: i for i, c in enumerate(comps)} # Make lookup dict for compartment indices - b = np.zeros((len(b_objs), 1)) + b = np.fromiter(b_vals.values(),dtype=float).reshape(-1,1) A = np.zeros((len(b_objs), len(comps))) - # Construct the characteristic value vector (b) and the includes matrix (A) - for i, obj in enumerate(b_objs): - # Look up the characteristic value - par = parset.pars[obj.name] - b[i] = par.interpolate(t_init, pop_name=self.name)[0] * par.y_factor[self.name] * par.meta_y_factor + # Fill out the includes matrix (A) + for i, obj in enumerate(b_objs.values()): if isinstance(obj, Characteristic): - if obj.denominator is not None: - denom_par = parset.pars[obj.denominator.name] - b[i] *= denom_par.interpolate(t_init, pop_name=self.name)[0] * denom_par.y_factor[self.name] * denom_par.meta_y_factor for inc in obj.get_included_comps(): - A[i, comp_indices[inc.name]] = 1.0 + if inc not in zero_compartments: + A[i, comp_indices[inc.name]] = 1.0 else: A[i, comp_indices[obj.name]] = 1.0 # Solve the linear system (nb. lstsq returns the minimum norm solution x = np.linalg.lstsq(A, b.ravel(), rcond=None)[0].reshape(-1, 1) + proposed = np.matmul(A, x) residual = np.sum((proposed.ravel() - b.ravel()) ** 2) @@ -1899,10 +1931,10 @@ def initialize_compartments(self, parset: ParameterSet, framework, t_init: float characteristic_tolerence_failed = False # Print warning for characteristics that are not well matched by the compartment size solution - for i in range(0, len(b_objs)): - if abs(proposed[i] - b[i]) > model_settings["tolerance"]: + for i, obj in enumerate(b_objs.values()): + if abs(proposed[i] - b[i]) > model_settings["initialization_tolerance"]: characteristic_tolerence_failed = True - error_msg += "Characteristic '{0}' '{1}' - Requested {2}, Calculated {3}\n".format(self.name, b_objs[i].name, b[i], proposed[i]) + error_msg += f"{obj.__class__.__name__} '{obj.name}' ({self.name})- Requested {b[i]}, Calculated {proposed[i]}\n" # Print expanded diagnostic for negative compartments showing parent characteristics def report_characteristic(charac, n_indent=0): @@ -1923,8 +1955,8 @@ def report_characteristic(charac, n_indent=0): """ msg = "" - if charac.name in charac_indices: - msg += n_indent * "\t" + "Characteristic '{0}': Target value = {1}\n".format(charac.name, b[charac_indices[charac.name]]) + if charac.name in b_indices: + msg += n_indent * "\t" + "Characteristic '{0}': Target value = {1}\n".format(charac.name, b[b_indices[charac.name]]) else: msg += n_indent * "\t" + "Characteristic '{0}' not in databook: Target value = N/A (0.0)\n".format(charac.name) @@ -1933,13 +1965,15 @@ def report_characteristic(charac, n_indent=0): for inc in charac.includes: if isinstance(inc, Characteristic): msg += report_characteristic(inc, n_indent) + elif inc in zero_compartments: + msg += n_indent * "\t" + "Compartment %s: Preassigned value = 0.0\n" % (inc.name) else: msg += n_indent * "\t" + "Compartment %s: Computed value = %f\n" % (inc.name, x[comp_indices[inc.name]]) return msg for i in range(0, len(comps)): - if x[i] < -model_settings["tolerance"]: - error_msg += "Compartment %s %s - Calculated %f\n" % (self.name, comps[i].name, x[i]) + if x[i] < -model_settings["initialization_tolerance"]: + error_msg += "Compartment %s %s - Calculated %f\n" % (self.name, comps[i].name, x[i,0]) for charac in b_objs: try: if comps[i] in charac.get_included_comps(): @@ -1948,10 +1982,10 @@ def report_characteristic(charac, n_indent=0): if comps[i] == charac: error_msg += report_characteristic(charac) - if residual > model_settings["tolerance"]: + if residual > model_settings["initialization_tolerance"]: # Halt for an unsatisfactory overall solution - raise BadInitialization("Global residual was %g which is unacceptably large (should be < %g)\n%s" % (residual, model_settings["tolerance"], error_msg)) - elif np.any(np.less(x, -model_settings["tolerance"])): + raise BadInitialization("Global residual was %g which is unacceptably large (should be < %g)\n%s" % (residual, model_settings["initialization_tolerance"], error_msg)) + elif np.any(np.less(x, -model_settings["initialization_tolerance"])): # Halt for any negative popsizes raise BadInitialization(f"Negative initial popsizes:\n{error_msg}") elif characteristic_tolerence_failed: @@ -1962,9 +1996,13 @@ def report_characteristic(charac, n_indent=0): # (but it exists as a fallback to ensure that any inconsistencies result in the error being raised) raise BadInitialization(f"Initialization error\n{error_msg}") - # Otherwise, insert the values + # Initialize any compartments that were specified as zero + for c in zero_compartments: + c[0] = 0.0 + + # Insert the calculated initial values for i, c in enumerate(comps): - c[0] = max(0.0, x[i]) + c[0] = max(0.0, x[i,0]) class Model: @@ -2640,7 +2678,7 @@ def update_pars(self) -> None: for par, val in zip(pars, par_vals): if par.skip_function is None or (self.t[ti] < par.skip_function[0]) or (self.t[ti] > par.skip_function[1]): # Careful - note how the < here matches >= in Parameter.update() - par[ti] = par.scale_factor * val + par[ti] = par.scale_factor * val[0] # Restrict the parameter's value if a limiting range was defined for par in pars: diff --git a/atomica/optimization.py b/atomica/optimization.py index f43319d0..7620a87f 100644 --- a/atomica/optimization.py +++ b/atomica/optimization.py @@ -998,7 +998,7 @@ def get_hard_constraint(self, optimization, instructions: ProgramInstructions) - # Lastly, apply the budget factor if len(self.budget_factor) == 1: - hard_constraints["initial_total_spend"][t] = total_spend * self.budget_factor + hard_constraints["initial_total_spend"][t] = total_spend * self.budget_factor.item() else: hard_constraints["initial_total_spend"][t] = total_spend * self.budget_factor[idx] @@ -1452,7 +1452,12 @@ def optimize(project, optimization, parset: ParameterSet, progset: ProgramSet, i errormsg = "PSO optimization requires finite upper and lower bounds to specify the search domain (i.e. every Adjustable needs to have finite bounds)" raise Exception(errormsg) - x_opt, _ = pyswarm.pso(_objective_fcn, kwargs=args, **optim_args) + if sc.compareversions(pyswarm, ">=1.0.0"): + x_opt = pyswarm.pso(_objective_fcn, kwargs=args, **optim_args).x + else: + # On Mac OS, Pyswarm 1.0.0 is not installing yet. This can be revisited and hopefully + # removed once the incompatibility is resolved (possibly by switching to more recent UV pipeline) + x_opt, _ = pyswarm.pso(_objective_fcn, kwargs=args, **optim_args) elif optimization.method == "hyperopt": diff --git a/atomica/parameters.py b/atomica/parameters.py index 5f668300..4f88b15f 100644 --- a/atomica/parameters.py +++ b/atomica/parameters.py @@ -299,7 +299,20 @@ def apply(self, pop, framework=None, parset=None) -> None: if (comp.name, pop.name) not in self.values: comp._vals[:, 0] = 0 else: - comp._vals[:, 0] = self.values[(comp.name, pop.name)] + vals = sc.promotetoarray(self.values[(comp.name, pop.name)]) + if comp._vals.shape[0] != vals.shape[0]: + if np.all(self.values[(comp.name, pop.name)] == 0): + # If there is a mismatch between the saved initialization duration and the + # duration for the current simulation AND if the values were all zero it is probably + # safe to assume the values can remain zero - this considerably increases flexibility in usage + comp._vals[:, 0] = 0 + else: + # Otherwise, although we could try and guess how to assign the values, there is a risk that the time points + # are different because the step size was changed rather than the duration being changed. Therefore, + # if the sizes don't match and any values are nonzero, simply raise an error + raise Exception(f'The saved initialization for "{comp.name}" ({pop.name}) has {len(self.values[(comp.name, pop.name)])} time points, but the current parameters lead to a timed compartment duration with {comp._vals.shape[0]} time points. As nonzero values are present, the initialization cannot be applied.') + else: + comp._vals[:, 0] = self.values[(comp.name, pop.name)] else: if (comp.name, pop.name) not in self.values: comp.vals[0] = 0 @@ -341,16 +354,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 @@ -490,7 +503,10 @@ def get_par(self, name: str, pop: str = None) -> Parameter: The parameter values for interactions and transfers are stored keyed by the source/from population. Thus, if the quantity name is an interaction or transfer, it is also necessary to specify the source population in order - to return a :class:`Parameter` instance. + to return a :class:`Parameter` instance. However, transfer parameters can + also be identified by a parameter name "_from_" such + that ``ParameterSet.get_par('age','5-14')`` is equivalent to + ``ParameterSet.get_par('age_from_5-14')``. :param name: The code name of a parameter, interaction, or transfer :param pop: @@ -507,7 +523,12 @@ def get_par(self, name: str, pop: str = None) -> Parameter: assert not pd.isna(pop), f'"{name}" is an interaction, so the ``pop`` must be specified' return self.interactions[name][pop] else: - raise KeyError(f'Parameter "{name}" not found') + for transfer in self.transfers.values(): + for par in transfer.values(): + if par.name == name: + return par + + raise KeyError(f'Parameter "{name}" not found') def sample(self, constant=True): """ @@ -534,7 +555,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)) @@ -677,5 +698,7 @@ def load_calibration(self, spreadsheet: sc.Spreadsheet) -> None: if "Initialization" in excelfile.sheet_names: self.initialization = Initialization.from_excel(excelfile) + else: + self.initialization = None logger.debug("Loaded calibration from file") diff --git a/atomica/plotting.py b/atomica/plotting.py index 9529ab8d..5749feab 100644 --- a/atomica/plotting.py +++ b/atomica/plotting.py @@ -43,6 +43,28 @@ settings["dpi"] = 150 # average quality settings["transparent"] = False +# Define default aggregation methods for standard units +DEFAULT_OUTPUT_AGGREGATIONS = { + FS.QUANTITY_TYPE_FRACTION: "average", + FS.QUANTITY_TYPE_PROPORTION: "average", + FS.QUANTITY_TYPE_PROBABILITY: "average", + FS.QUANTITY_TYPE_RATE: "average", + FS.QUANTITY_TYPE_DURATION: "sum", # e.g., diagnosis time + treatment time should be summed within a population + FS.QUANTITY_TYPE_NUMBER: "sum", + "Number of people": "sum", # This unit is used by compartments and characteristics without denominators +} + +DEFAULT_POP_AGGREGATIONS = { + FS.QUANTITY_TYPE_FRACTION: "weighted", + FS.QUANTITY_TYPE_PROPORTION: "weighted", + FS.QUANTITY_TYPE_PROBABILITY: "weighted", + FS.QUANTITY_TYPE_RATE: "weighted", + FS.QUANTITY_TYPE_DURATION: "weighted", # e.g., diagnosis time should be averaged across populations + FS.QUANTITY_TYPE_NUMBER: "sum", + "Number of people": "sum", # This unit is used by compartments and characteristics without denominators + "": "weighted", # This unit is used by characteristics with denominators +} + def save_figs(figs, path=".", prefix="", fnames=None, file_format="png") -> None: """ @@ -139,20 +161,41 @@ class PlotData: :param output_aggregation: If an output aggregation is requested, combine the outputs listed using one of - - 'sum' - just add values together - - 'average' - unweighted average of quantities + - 'sum' - just add values together. If summing together characteristics with a denominator, then + all of the items being included in the sum must have the same denominator. This denominator + will then continue to be tracked for the sum. + - 'average' - unweighted average of quantities. If averaging characteristics with a denominator, then + the denominator will be dropped afterwards (so subsequent population aggregation will be + performed based on total population size) - 'weighted' - weighted average where the weight is the - compartment size, characteristic value, or link source - compartment size (summed over duplicate links). 'weighted' - method cannot be used with non-transition parameters and a - KeyError will result in that case - - :param pop_aggregation: Same as output_aggregation, except that 'weighted' - uses population sizes. Note that output aggregation is performed - before population aggregation. This also means that population - aggregation can be used to combine already aggregated outputs (e.g. - can first sum 'sus'+'vac' within populations, and then take weighted - average across populations) + - Link source compartment size, or + - Characteristic denominator size, or + - Compartment size + depending on the quantities being aggregated. If one quantity has a denominator, then + all quantities must have denominators. This method cannot be used with non-transition parameters and a + KeyError will result in that case. Note that if using this method to aggregate characteristics, there + should generally not be any overlap in the compartments contributing to the denominators. + + For aggregating characteristics, the behaviour for combining characteristics ``a/b`` and ``c/d`` is as follows: + + - 'sum' will give ``(a+c)/b`` and raise an error if ``b!=d``. The denominator ``b`` can be used for weighted population aggregation. + - 'average' will give ``((a/b)+(c/d))/2``. No denominator is tracked for weighted population aggregation, so such aggregation would use the total population size. + - 'weighted' will give ``(a+c)/(b+d)``. The denominator ``b+d`` can be used for weighted population aggregation. + + The default is 'weighted' for dimensionless units, probabilities, rates, and fractions, and 'sum' for everything else, with the exception of characteristics + with denominators, where + - 'sum' is used if the characteristics all have the same denominator + - 'weighted' is used if the characteristics all have different denominators + - 'average' is used otherwise + + :param pop_aggregation: Same options as output_aggregation, except + - 'weighted' uses population sizes OR characteristic denominator size. Note that output aggregation + is performed before population aggregation. This also means that population + aggregation can be used to combine already aggregated outputs (e.g. + can first sum 'sus'+'vac' within populations, and then take weighted + average across populations). This weighting is principally intended for aggregating + nondimensional quantities (fractions, proportions etc.). + - By default, quantities in 'duration' units are averaged across populations rather than summed. :param project: Optionally provide a :class:`Project` object, which will be used to convert names to labels in the outputs for plotting. @@ -203,6 +246,7 @@ def __init__(self, results, outputs=None, pops=None, output_aggregation=None, po assert output_aggregation in [None, "sum", "average", "weighted"] assert pop_aggregation in [None, "sum", "average", "weighted"] + assert time_aggregation in [None, "integrate", "average"] # First, get all of the pops and outputs requested by flattening the lists pops_required = _extract_labels(pops) @@ -221,9 +265,12 @@ def __init__(self, results, outputs=None, pops=None, output_aggregation=None, po aggregated_outputs = defaultdict(dict) # Dict with aggregated_outputs[pop_label][aggregated_output_label] aggregated_units = dict() # Dict with aggregated_units[aggregated_output_label] aggregated_timescales = dict() + aggregated_denominators = defaultdict(dict) # Store aggregated denominators where available at the same dimensionality as outputs output_units = dict() output_timescales = dict() compsize = dict() + denoms = sc.odict() # Record denominators associated with input characteristics + popsize = dict() # Defaultdict won't throw key error when checking outputs. data_label = defaultdict(str) # Label used to identify which data to plot, maps output label to data label. @@ -241,7 +288,10 @@ def __init__(self, results, outputs=None, pops=None, output_aggregation=None, po vars = pop.get_variable(output_label) except NotFoundError as e: in_pops = [x.name for x in result.model.pops if output_label in x] - message = f'Variable "{output_label}" was requested in population "{pop.name}" but it is only defined in these populations: {in_pops}' + if not in_pops: + message = f'Variable "{output_label}" was requested in population "{pop.name}" but it is not defined in any populations' + else: + message = f'Variable "{output_label}" was requested in population "{pop.name}" but it is only defined in these populations: {in_pops}' raise NotFoundError(message) from e if vars[0].vals is None: @@ -280,7 +330,8 @@ def __init__(self, results, outputs=None, pops=None, output_aggregation=None, po output_units[output_label] = vars[0].units output_timescales[output_label] = None data_label[output_label] = vars[0].name - + if isinstance(vars[0], Characteristic) and vars[0].denominator is not None: + denoms[output_label] = vars[0].denominator # Record the denominator for this quantity else: raise Exception("Unknown type") @@ -303,10 +354,9 @@ def placeholder_pop(): par = Parameter(pop=placeholder_pop, name=output_label) fcn, dep_labels = parse_function(f_stack_str) deps = {} - displayed_annualization_warning = False for dep_label in dep_labels: vars = pop.get_variable(dep_label) - if t_bins is not None and (isinstance(vars[0], Link) or isinstance(vars[0], Parameter)) and time_aggregation == "sum" and not displayed_annualization_warning: + if t_bins is not None and (isinstance(vars[0], Link) or isinstance(vars[0], Parameter)) and time_aggregation == "integrate": raise Exception("Function includes Parameter/Link so annualized rates are being used. Aggregation should therefore use 'average' rather than 'sum'.") deps[dep_label] = vars par._fcn = fcn @@ -320,7 +370,7 @@ def placeholder_pop(): # Third pass, aggregate them according to any aggregations present for output in outputs: # For each final output if isinstance(output, dict): - output_name = list(output.keys())[0] + output_name = list(output.keys())[0] # nb. there should be only one item in the output dict at this point labels = output[output_name] # If this was a function, aggregation over outputs doesn't apply so just put it straight in. @@ -328,24 +378,44 @@ def placeholder_pop(): aggregated_outputs[pop_label][output_name] = data_dict[output_name] aggregated_units[output_name] = "unknown" # Also, we don't know what the units of a function are aggregated_timescales[output_name] = None # Timescale is lost + aggregated_denominators[pop_label][output_name] = None # Denominators are not tracked continue units = list(set([output_units[x] for x in labels])) timescales = list(set([np.nan if isna(output_timescales[x]) else output_timescales[x] for x in labels])) # Ensure that None and nan don't appear as different timescales - # Set default aggregation method depending on the units of the quantity - if output_aggregation is None: - if units[0] in ["", FS.QUANTITY_TYPE_FRACTION, FS.QUANTITY_TYPE_PROPORTION, FS.QUANTITY_TYPE_PROBABILITY, FS.QUANTITY_TYPE_RATE]: - output_aggregation = "average" + # Set default aggregation method depending on the units of the quantity. If the aggregation method has not been specified + # then select a default aggregation method separately for each output + if output_aggregation is not None: + aggregation = output_aggregation + elif denoms: + try: + denominators = set([denoms[x].name for x in labels]) + except KeyError as e: + raise Exception("If an aggregation includes characteristics with denominators, then all of the included quantities must be characteristics with denominators") from e + + if len(denominators) == 1: + # If all characteristics share the same denominator e.g., SP and SN prevalence, then we want to add the characteristics together + aggregation = "sum" + elif len(denominators) == len(labels): + # If all characteristics have different denominators e.g., proportion of SP/SN on treatment aggregated over SP and SN, then use weighted averaging + aggregation = "weighted" else: - output_aggregation = "sum" + # Otherwise if mixing denominators, use a direct average. In general this should not happen because the result is not readily interpretable + aggregation = "average" + elif units[0] in DEFAULT_OUTPUT_AGGREGATIONS: + aggregation = DEFAULT_OUTPUT_AGGREGATIONS[units[0]] + else: + logger.warning(f'Unit quantity "{units[0]}" was not recognized and an output aggregation method was not specified, using "sum" by default') + aggregation = "sum" if len(units) > 1: logger.warning("Aggregation for output '%s' is mixing units, this is almost certainly not desired.", output_name) aggregated_units[output_name] = "unknown" else: - if units[0] in ["", FS.QUANTITY_TYPE_FRACTION, FS.QUANTITY_TYPE_PROPORTION, FS.QUANTITY_TYPE_PROBABILITY, FS.QUANTITY_TYPE_RATE] and output_aggregation == "sum" and len(labels) > 1: # Dimensionless, like prevalance - logger.warning("Output '%s' is not in number units, so output aggregation probably should not be 'sum'.", output_name) + # If the user has requested a sum for a unit that doesn't use summation by default, there's a decent chance that this is a mistake + if aggregation == "sum" and DEFAULT_OUTPUT_AGGREGATIONS.get(units[0], "sum") != "sum" and len(labels) > 1 and not labels[0] in denoms: # Dimensionless, like prevalance + logger.warning(f"Output '{output_name}' is in '{units[0]}' units, so output aggregation probably should not be 'sum'.") aggregated_units[output_name] = output_units[labels[0]] if len(timescales) > 1: @@ -354,18 +424,56 @@ def placeholder_pop(): else: aggregated_timescales[output_name] = output_timescales[labels[0]] - if output_aggregation == "sum": - aggregated_outputs[pop_label][output_name] = sum(data_dict[x] for x in labels) # Add together all the outputs - elif output_aggregation == "average": + # Perform additional validation on the denominators + if denoms: + unique_denoms = set([denoms[x].name for x in labels]) + + if len(unique_denoms) > 1 and len(unique_denoms) < len(labels): + logger.warning("When aggregating characteristics with denominators, they should generally either all have the same denominator, or all have different denominators. Partially duplicate denominators is likely not desired.") + + if len(unique_denoms) > 1: + # Check that the denominators are non-overlapping + comps = [set(result.framework.get_charac_includes(x)) for x in labels] + if sum(len(x) for x in comps) != len(set().union(*comps)): + logger.warning("When aggregating characteristics with different denominators, generally there should be no overlap in the denominators. However, some compartments are shared between the requested characteristic denominators. Unless this is intentional, this is likely a mistake that should be investigated.") + + if aggregation == "sum": aggregated_outputs[pop_label][output_name] = sum(data_dict[x] for x in labels) # Add together all the outputs + + # If summing characteristics that have denominators, then all quantities should have denominators, and they should all be equal + # For example, summing SP-TB and SN-TB prevalance where they both have the 'alive' denominator + has_denominator = [x in denoms for x in labels] + if all(has_denominator): + assert len(set([denoms[x].name for x in labels])) == 1, 'When aggregating characteristics with denominators, if summing them then all quantities must have the same denominator. The "weighted" aggregation can be used to combine characteristics with different denominators, on the assumption that their denominators do not overlap.' + aggregated_denominators[pop_label][output_name] = denoms[0].vals + elif not any(has_denominator): + aggregated_denominators[pop_label][output_name] = None + else: + raise Exception(f"When aggregating outputs, if any quantities being aggregated have denominators, then they must all have denominators. These quantities have denominators {[x for x in labels if x in denoms]} while these do not {[x for x in labels if x not in denoms]}") + + elif aggregation == "average": + aggregated_outputs[pop_label][output_name] = sum(data_dict[x] for x in labels) aggregated_outputs[pop_label][output_name] /= len(labels) - elif output_aggregation == "weighted": - aggregated_outputs[pop_label][output_name] = sum(data_dict[x] * compsize[x] for x in labels) # Add together all the outputs - aggregated_outputs[pop_label][output_name] /= sum([compsize[x] for x in labels]) + aggregated_denominators[pop_label][output_name] = None # If taking a direct average then drop the denominator + elif aggregation == "weighted": + has_denominator = [x in denoms for x in labels] + assert all(has_denominator) or not any(has_denominator) # An output aggregation must either have no quantities with denominators, or only contain quantities with denominators + if all(has_denominator): + # Take weighted average based on denominator size, and then store the combined denominator further population aggregation + aggregated_outputs[pop_label][output_name] = sum(data_dict[x] * denoms[x].vals for x in labels) + aggregated_denominators[pop_label][output_name] = sum([denoms[x].vals for x in labels]) + aggregated_outputs[pop_label][output_name] /= aggregated_denominators[pop_label][output_name] + elif not any(has_denominator): + aggregated_outputs[pop_label][output_name] = sum(data_dict[x] * compsize[x] for x in labels) + aggregated_outputs[pop_label][output_name] /= sum([compsize[x] for x in labels]) + aggregated_denominators[pop_label][output_name] = None # If taking a direct average then drop the denominator + else: + raise Exception(f"When aggregating outputs, if any quantities being aggregated have denominators, then they must all have denominators. These quantities have denominators {[x for x in labels if x in denoms]} while these do not {[x for x in labels if x not in denoms]}") else: aggregated_outputs[pop_label][output] = data_dict[output] aggregated_units[output] = output_units[output] aggregated_timescales[output] = output_timescales[output] + aggregated_denominators[pop_label][output] = denoms[output].vals if output in denoms else None # Now aggregate over populations # If we have requested a reduction over populations, this is done for every output present @@ -375,26 +483,39 @@ def placeholder_pop(): pop_name = list(pop.keys())[0] pop_labels = pop[pop_name] - # Set population aggregation method depending on - if pop_aggregation is None: - if aggregated_units[output_name] in ["", FS.QUANTITY_TYPE_FRACTION, FS.QUANTITY_TYPE_PROPORTION, FS.QUANTITY_TYPE_PROBABILITY, FS.QUANTITY_TYPE_RATE]: - pop_aggregation = "average" - else: - pop_aggregation = "sum" + # Set population aggregation method depending on the quantity being aggregated + if pop_aggregation is not None: + aggregation = pop_aggregation + elif aggregated_denominators[pop_labels[0]][output_name] is not None: + # Outputs with denominators use weighted by default + aggregation = "weighted" + elif aggregated_units[output_name] in DEFAULT_POP_AGGREGATIONS: + aggregation = DEFAULT_POP_AGGREGATIONS[aggregated_units[output_name]] + else: + logger.warning(f'Unit quantity "{aggregated_units[output_name]}" was not recognized and an output aggregation method was not specified, using "sum" by default') + aggregation = "sum" - if pop_aggregation == "sum": - if aggregated_units[output_name] in ["", FS.QUANTITY_TYPE_FRACTION, FS.QUANTITY_TYPE_PROPORTION, FS.QUANTITY_TYPE_PROBABILITY, FS.QUANTITY_TYPE_RATE] and len(pop_labels) > 1: - logger.warning("Output '%s' is not in number units, so population aggregation probably should not be 'sum'", output_name) + # Perform aggregation + if aggregation == "sum": + if aggregation == "sum" and DEFAULT_POP_AGGREGATIONS.get(aggregated_units[output_name], "sum") != "sum" and len(pop_labels) > 1: + logger.warning(f"Output '{output_name}' is in '{aggregated_units[output_name]}' units, so output aggregation should probably be '{DEFAULT_POP_AGGREGATIONS.get(aggregated_units[output_name])}', not 'sum'.") vals = sum(aggregated_outputs[x][output_name] for x in pop_labels) # Add together all the outputs - elif pop_aggregation == "average": + elif aggregation == "average": + if aggregated_denominators[pop_labels[0]][output_name] is not None: + logger.warning(f"Output '{output_name}' has a denominator, so output aggregation should probably be 'weighted' rather than 'average'") vals = sum(aggregated_outputs[x][output_name] for x in pop_labels) # Add together all the outputs vals /= len(pop_labels) - 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]) + elif aggregation == "weighted": + if aggregated_denominators[pop_labels[0]][output_name] is not None: + weights = {x: aggregated_denominators[x][output_name] for x in pop_labels} + else: + weights = {x: popsize[x] for x in pop_labels} + numerator = sum(aggregated_outputs[x][output_name] * weights[x] for x in pop_labels) # Add together all the outputs + denominator = sum([weights[x] for x in pop_labels]) vals = np.divide(numerator, denominator, out=np.full(numerator.shape, np.nan, dtype=float), where=numerator != 0) else: - raise Exception("Unknown pop aggregation method") + raise Exception(f'Unknown population aggregation method "{aggregation}"') + 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)) else: vals = aggregated_outputs[pop][output_name] @@ -501,12 +622,16 @@ def time_aggregate(self, t_bins, time_aggregation=None, interpolation_method=Non interpolation_method = "linear" if not hasattr(t_bins, "__len__"): - # If a scalar bin is provided, then it is + # If a scalar bin is provided, then it is assumed to correspond to the bin width, in which case we need + # to create uniformly spaced bins up to the last complete bin. if t_bins > (self.series[0].tvec[-1] - self.series[0].tvec[0]): # If bin width is greater than the sim duration, treat it the same as aggregating over all times t_bins = "all" else: if not (self.series[0].tvec[-1] - self.series[0].tvec[0]) % t_bins: + # If the simulation ends in a multiple of the t_bins, then we can have a final bin extending up to the last + # simulation year. If we just use tvec[-1] in np.arange then that final bin will be excluded. So we add + # t_bins on to that final bin to ensure it's included. upper = self.series[0].tvec[-1] + t_bins else: upper = self.series[0].tvec[-1] @@ -539,19 +664,46 @@ def time_aggregate(self, t_bins, time_aggregation=None, interpolation_method=Non else: scale = 1.0 - # We interpolate in time-aggregation because the time bins are independent of the step size. In contrast, - # accumulation preserves the same time bins, so we don't need the interpolation step and instead go straight - # to summation or trapezoidal integration - max_step = 0.5 * min(np.diff(s.tvec)) # Subdivide for trapezoidal integration with at least 2 divisions per timestep. Could be a lot of memory for integrating daily timesteps over a full simulation, but unlikely to be prohibitive vals = np.full(lower.shape, fill_value=np.nan) - for i, (l, u) in enumerate(zip(lower, upper)): - n = np.ceil((u - l) / max_step) + 1 # Add 1 so that in most cases, we can use the actual timestep values - t2 = np.linspace(l, u, int(n)) + lower_idx = np.searchsorted(s.tvec, lower, side="left") # Time index for bin start + upper_idx = np.searchsorted(s.tvec, upper, side="right") # Time index for bin end + + for i, (l, u, l_idx, u_idx) in enumerate(zip(lower, upper, lower_idx, upper_idx)): + + # For bins that partially extend out of bounds, return NaN as the value immediately + if l < s.tvec[0] or u > s.tvec[-1]: + vals[i] = np.nan + continue + elif l == u: + vals[i] = 0 + continue + + # The bins will consist of the actual simulation time points, plus + # partial bins that are interpolated before and after if the requested + # bins don't line up with the simulation timepoints + idx = np.arange(l_idx, u_idx) + t2 = s.tvec[idx] + interpolate = False + if t2[0] > l or t2[-1] < u: + interpolate = True + t2 = list(t2) + if t2[0] > l: + t2.insert(0, l) + if t2[-1] < u: + t2.append(u) + t2 = np.array(t2, dtype=float) + if interpolation_method == "linear": - v2 = np.interp(t2, s.tvec, s.vals, left=np.nan, right=np.nan) # Return NaN outside bounds - it should never be valid to use extrapolated output values in time aggregation - vals[i] = np.trapz(y=v2 / scale, x=t2) # Note division by timescale here, which annualizes it + if interpolate: + v2 = np.interp(t2, s.tvec, s.vals, left=np.nan, right=np.nan) # Return NaN outside bounds - it should never be valid to use extrapolated output values in time aggregation + else: + v2 = s.vals[idx] + vals[i] = np.trapezoid(y=v2 / scale, x=t2) # Note division by timescale here, which annualizes it elif interpolation_method == "previous": - v2 = scipy.interpolate.interp1d(s.tvec, s.vals, kind="previous", copy=False, assume_sorted=True, bounds_error=False, fill_value=(np.nan, np.nan))(t2) + if interpolate: + v2 = scipy.interpolate.interp1d(s.tvec, s.vals, kind="previous", copy=False, assume_sorted=True, bounds_error=False, fill_value=(np.nan, np.nan))(t2) + else: + v2 = s.vals[idx] vals[i] = sum(v2[:-1] / scale * np.diff(t2)) s.tvec = (lower + upper) / 2.0 @@ -590,7 +742,7 @@ def time_aggregate(self, t_bins, time_aggregation=None, interpolation_method=Non if sc.isstring(t_bins) and t_bins == "all": s.t_labels = ["All"] else: - s.t_labels = ["%d-%d" % (low, high) for low, high in zip(lower, upper)] + s.t_labels = [f"{np.format_float_positional(low,trim='-')}-{np.format_float_positional(high,trim='-')}" for low, high in zip(lower, upper)] return self @@ -1088,6 +1240,70 @@ def interpolate(self, new_tvec): return np.interp(sc.promotetoarray(new_tvec), self.tvec, self.vals, left=np.nan, right=np.nan) +# Temporary copy of function from Sciris to remove after Sciris update +def _get_legend_handles(ax, handles, labels): + """ + Construct handle and label list, from one of: + + - A list of handles and a list of labels + - A list of handles, where each handle contains the label + - An axis object, containing the objects that should appear in the legend + - A figure object, from which the first axis will be used + """ + if handles is None: + if ax is None: + ax = plt.gca() + elif isinstance(ax, plt.Figure): # Allows an argument of a figure instead of an axes # pragma: no cover + ax = ax.axes[-1] + handles, labels = ax.get_legend_handles_labels() + else: # pragma: no cover + if labels is None: + labels = [h.get_label() for h in handles] + else: + assert len(handles) == len(labels), f"Number of handles ({len(handles)}) and labels ({len(labels)}) must match" + return ax, handles, labels + + +# Temporary copy of function from Sciris to remove after Sciris update +def separatelegend(ax=None, handles=None, labels=None, reverse=False, figsettings=None, legendsettings=None): + """Allows the legend of a figure to be rendered in a separate window instead""" + + # Handle settings + f_settings = sc.mergedicts({"figsize": (4.0, 4.8)}, figsettings) # (6.4,4.8) is the default, so make it a bit narrower + l_settings = sc.mergedicts({"loc": "center", "bbox_to_anchor": None, "frameon": False}, legendsettings) + + # Get handles and labels + _, handles, labels = _get_legend_handles(ax, handles, labels) + + # Set up new plot + fig = plt.figure(**f_settings) + ax = fig.add_subplot(111) + ax.set_position([-0.05, -0.05, 1.1, 1.1]) # This cuts off the axis labels, ha-ha + ax.set_axis_off() # Hide axis lines + + # A legend renders the line/patch based on the object handle. However, an object + # can only appear in one figure. Thus, if the legend is in a different figure, the + # object cannot be shown in both the original figure and in the legend. Thus we need + # to copy the handles, and use the copies to render the legend + handles2 = [] + for h in handles: + h2 = sc.cp(h) + h2.axes = None + h2._parent_figure = None + h2.figure = None + handles2.append(h2) + + # Reverse order, e.g. for stacked plots + if reverse: # pragma: no cover + handles2 = handles2[::-1] + labels = labels[::-1] + + # Plot the new legend + ax.legend(handles=handles2, labels=labels, **l_settings) + + return fig + + def plot_bars(plotdata, stack_pops=None, stack_outputs=None, outer=None, legend_mode=None, show_all_labels=False, orientation="vertical") -> list: """ Produce a bar plot @@ -1434,7 +1650,7 @@ def process_input_stacks(input_stacks, available_items): if legend_mode == "together": _render_legend(ax, plot_type="bar", handles=legend_patches) elif legend_mode == "separate": - figs.append(sc.separatelegend(handles=legend_patches, reverse=True)) + figs.append(separatelegend(handles=legend_patches, reverse=True)) return figs @@ -1615,7 +1831,7 @@ def _prepare_figures(dim1, dim2, n_cols): if not subplots: # Replace the last figure with a legend figure plt.close(figs[-1]) # TODO - update Sciris to allow passing in an existing figure - figs[-1] = sc.separatelegend(ax, reverse=reverse_legend) + figs[-1] = separatelegend(ax, reverse=reverse_legend) else: legend_ax = axes[-1] handles, labels = ax.get_legend_handles_labels() @@ -1735,7 +1951,7 @@ def plot_legend(entries: dict, plot_type=None, fig=None, legendsettings: dict = raise Exception(f'Unknown plot type "{p_type}"') if fig is None: # Draw in a new figure - fig = sc.separatelegend(handles=h, legendsettings=legendsettings) + fig = separatelegend(handles=h, legendsettings=legendsettings) else: existing_legend = fig.findobj(Legend) if existing_legend and existing_legend[0].parent is fig: # If existing legend and this is a separate legend fig @@ -1992,7 +2208,7 @@ def _extract_labels(input_arrays) -> set: for x in input_arrays: if isinstance(x, dict): k = list(x.keys()) - assert len(k) == 1, "Aggregation dict can only have one key" + assert len(k) == 1, "Aggregation dict can only have one key, for multiple aggregations, use a list of dicts" if sc.isstring(x[k[0]]): continue else: diff --git a/atomica/programs.py b/atomica/programs.py index f3f60f1f..02640a22 100644 --- a/atomica/programs.py +++ b/atomica/programs.py @@ -1552,7 +1552,7 @@ def get_outcome(self, prop_covered): # RANDOM CALCULATION elif self.cov_interaction == "random": # Outcome += c1(1-c2)* delta_out1 + c2(1-c1)*delta_out2 + c1c2* max(delta_out1,delta_out2) - combination_coverage = np.product(self.combinations * cov + (self.combinations ^ 1) * (1 - cov), axis=1) + combination_coverage = np.prod(self.combinations * cov + (self.combinations ^ 1) * (1 - cov), axis=1) outcome += np.sum(combination_coverage.ravel() * self._combination_outcomes.ravel()) else: raise Exception('Unknown reachability type "%s"', self.cov_interaction) diff --git a/atomica/project.py b/atomica/project.py index e6881957..a7911dd4 100644 --- a/atomica/project.py +++ b/atomica/project.py @@ -88,14 +88,13 @@ def tvec(self) -> np.ndarray: """ Return simulation time vector - This method uses `linspace` rather than `arange` to avoid accumulating numerical errors that prevent - integer years aligning exactly. + This method uses `arange` with a step size of 1, then multiplies with the correct step size and add sim start time. :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 def update_time_vector(self, start: float = None, end: float = None, dt: float = None) -> None: """ @@ -260,7 +259,7 @@ def load_databook(self, databook_path, make_default_parset: bool = True, do_run: if isinstance(databook_path, ProjectData): # If the user passed in a ProjectData instance, use it directly data = databook_path - databook = data.to_spreadsheet() + databook = None else: if isinstance(databook_path, sc.Spreadsheet): databook = databook_path @@ -577,7 +576,7 @@ def run_sampled_sims(self, parset, progset=None, progset_instructions=None, resu return results - def calibrate(self, parset=None, adjustables=None, measurables=None, max_time=60, save_to_project=False, new_name=None, default_min_scale=0.0, default_max_scale=2.0, default_weight=1.0, default_metric="fractional", **kwargs) -> ParameterSet: + def calibrate(self, parset=None, adjustables=None, measurables=None, save_to_project=False, new_name=None, default_min_scale=0.0, default_max_scale=2.0, default_weight=1.0, default_metric="fractional", yaml=None, **kwargs) -> ParameterSet: """ Method to perform automatic calibration. @@ -603,6 +602,23 @@ def calibrate(self, parset=None, adjustables=None, measurables=None, max_time=60 To calibrate a project-attached parameter set in place, provide its key as the new name argument to this method. Current fitting metrics are: "fractional", "meansquare", "wape" Note that scaling limits are absolute, not relative. + + Additional keyword arguments will be passed to `at.calibrate` so additional arguments such as `max_time` + or `methods` can be provided to this function. + + :param parset: + :param adjustables: + :param measurables: + :param save_to_project: + :param new_name: + :param default_min_scale: + :param default_max_scale: + :param default_weight: + :param default_metric: + :param yaml: Optionally specify a YAML file or YAML calibration dictionary to use for calibration + :param kwargs: + :return: + """ if parset is None: @@ -610,23 +626,37 @@ def calibrate(self, parset=None, adjustables=None, measurables=None, max_time=60 elif not isinstance(parset, ParameterSet): parset = self.parsets[parset] + if yaml is not None: + import atomica.yaml_calibration # Avoid circular import + + assert adjustables is None, "If a YAML file is specified, adjustables should not be set" + assert measurables is None, "If a YAML file is specified, measurables should not be set" + assert "time_period" not in kwargs, "If a YAML file is specified, time_period should not be set - instead, set cal_start and cal_end in the YAML file" + new_parset = atomica.yaml_calibration.run(yaml, self, parset, **kwargs) + else: + if adjustables is None: + adjustables = list(self.framework.pars.index[~self.framework.pars["calibrate"].isnull()]) + adjustables += list(self.framework.comps.index[~self.framework.comps["calibrate"].isnull()]) + adjustables += list(self.framework.characs.index[~self.framework.characs["calibrate"].isnull()]) + + if measurables is None: + measurables = list(self.framework.comps.index) + measurables += list(self.framework.characs.index) + + for index, adjustable in enumerate(adjustables): + if sc.isstring(adjustable): # Assume that a parameter name was passed in if not a tuple. + adjustables[index] = (adjustable, None, default_min_scale, default_max_scale) + + for index, measurable in enumerate(measurables): + if sc.isstring(measurable): # Assume that a parameter name was passed in if not a tuple. + measurables[index] = (measurable, None, default_weight, default_metric) + + new_parset = calibrate(project=self, parset=parset, pars_to_adjust=adjustables, output_quantities=measurables, **kwargs) + if new_name is None: new_name = parset.name + " (auto-calibrated)" - if adjustables is None: - adjustables = list(self.framework.pars.index[~self.framework.pars["calibrate"].isnull()]) - adjustables += list(self.framework.comps.index[~self.framework.comps["calibrate"].isnull()]) - adjustables += list(self.framework.characs.index[~self.framework.characs["calibrate"].isnull()]) - if measurables is None: - measurables = list(self.framework.comps.index) - measurables += list(self.framework.characs.index) - for index, adjustable in enumerate(adjustables): - if sc.isstring(adjustable): # Assume that a parameter name was passed in if not a tuple. - adjustables[index] = (adjustable, None, default_min_scale, default_max_scale) - for index, measurable in enumerate(measurables): - if sc.isstring(measurable): # Assume that a parameter name was passed in if not a tuple. - measurables[index] = (measurable, None, default_weight, default_metric) - new_parset = calibrate(project=self, parset=parset, pars_to_adjust=adjustables, output_quantities=measurables, max_time=max_time, **kwargs) new_parset.name = new_name # The new parset is a calibrated copy of the old, so change id. + if save_to_project: self.parsets.append(new_parset) diff --git a/atomica/reconciliation.py b/atomica/reconciliation.py index e87f6e8e..e27627ea 100644 --- a/atomica/reconciliation.py +++ b/atomica/reconciliation.py @@ -161,7 +161,7 @@ def _objective(x, mapping, progset, eval_years, target_vals, num_eligible, dt): return obj -def _convert_to_single_year(progset, reconciliation_year): +def _convert_to_single_year(progset, reconciliation_year: float): # Take in a progset # Return a progset with values only in the reconciliation year # This is then what actually gets reconciled @@ -242,11 +242,15 @@ def reconcile(project, parset, progset, reconciliation_year: float, max_time=10, logger.warning("Reconcilation when parameter is in number units not fully tested") - reconciliation_year = sc.promotetoarray(reconciliation_year) - assert len(reconciliation_year) == 1, "Reconciliation year must be a scalar" + # Convert reconciliation year to a scalar + if hasattr(reconciliation_year, "__len__"): + if len(reconciliation_year) == 1: + reconciliation_year = reconciliation_year[0] + else: + raise ValueError("Reconciliation year must be a scalar") if eval_range is None: - eval_range = [reconciliation_year[0], reconciliation_year[0] + project.settings.sim_dt] + eval_range = [reconciliation_year, reconciliation_year + project.settings.sim_dt] # Do a prerun to get the baseline values and coverage denominator parset_results = project.run_sim(parset=parset, progset=progset, store_results=False) @@ -255,7 +259,7 @@ def reconcile(project, parset, progset, reconciliation_year: float, max_time=10, target_vals, num_eligible = _extract_targets(parset_results, progset, ti, eval_pars) # Prepare ASD inputs - new_progset = _convert_to_single_year(progset, reconciliation_year[0]) + new_progset = _convert_to_single_year(progset, reconciliation_year) bounds = _prepare_bounds(new_progset, unit_cost_bounds, baseline_bounds, capacity_bounds, outcome_bounds) # Convert reconcile() inputs into full detailed bounds x0, xmin, xmax, mapping = _prepare_asd_inputs(new_progset, bounds) # Assemble ASD variables diff --git a/atomica/system.py b/atomica/system.py index c72ab22a..9e3c149c 100644 --- a/atomica/system.py +++ b/atomica/system.py @@ -89,3 +89,5 @@ class FrameworkSettings: RESERVED_KEYWORDS += supported_functions.keys() RESERVED_SYMBOLS = set(":,;/+-*'\" @") # A code_name in the framework (for characs, comps, pars) cannot contain any of these characters + + DEFAULT_PROVENANCE = "Framework-supplied default" diff --git a/atomica/utils.py b/atomica/utils.py index ee342e54..c20dfd7c 100644 --- a/atomica/utils.py +++ b/atomica/utils.py @@ -36,6 +36,7 @@ "parallel_progress", "start_logging", "stop_logging", + "get_sigfigs_necessary", ] @@ -376,8 +377,15 @@ def insert(self, t, v) -> None: :param v: Value to insert. If ``None``, this function will return immediately without doing anything """ + # TODO - could potentially incorporate iterability check above for greater efficiency + # TODO - add handling series inputs here too? + + # Convert 1-element lists/arrays to scalars + if hasattr(t, "__len__") and len(t) == 1: + t = t[0] + if hasattr(v, "__len__") and len(v) == 1: + v = v[0] - # Check if inputs are iterable iterable_input = True try: assert len(t) == len(v), "Cannot insert non-matching lengths or types of time and values %s and %s" % (t, v) @@ -489,7 +497,7 @@ def remove_after(self, t_remove) -> None: """ Remove times from start - :param tval: Remove times up to but not including this time + :param tval: Remove times after but not including this time """ @@ -623,6 +631,22 @@ def sample(self, constant=True): return new + def clear(self) -> None: + """ + Clear the TimeSeries + + This method resets TimeSeries but keeps the units present, so it is useful when + wanting to replace the values associated with a particular TimeSeries instance. + + :return: None + """ + + self.t = [] + self.vals = [] + self.assumption = None + self.sigma = None + self._sampled = False + return None def evaluate_plot_string(plot_string: str): """ @@ -656,7 +680,7 @@ def evaluate_plot_string(plot_string: str): fcn_ast = ast.parse(plot_string, mode="eval") for node in ast.walk(fcn_ast): if not (node is fcn_ast): - assert isinstance(node, ast.Dict) or isinstance(node, ast.Str) or isinstance(node, ast.List) or isinstance(node, ast.Load), "Only allowed to initialize lists and dicts of strings here" + assert isinstance(node, ast.Dict) or isinstance(node, ast.List) or isinstance(node, ast.Load) or isinstance(node, ast.Constant), "Only allowed to initialize lists and dicts of strings here" compiled_code = compile(fcn_ast, filename="", mode="eval") return eval(compiled_code) else: @@ -982,3 +1006,21 @@ 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 get_sigfigs_necessary(x, y, min_sigfigs: int = 2) -> int: + """ + Get how many significant figures are necessary to tell the difference between two numbers + + :param x, y: numbers to compare + :param min_sigfigs: minimum number of sigfigs to use if no difference + :return: Number of significant figures required + """ + msf = min_sigfigs + assert sc.isnumber(x) and sc.isnumber(y), f"Cannot compare sigfigs as {x} and {y} are not both numbers" + if x == y or np.isnan(x) or np.isnan(y): + return msf + else: + while sc.sigfig(x=x, sigfigs=msf) == sc.sigfig(x=y, sigfigs=msf): + msf += 1 + return msf diff --git a/atomica/version.py b/atomica/version.py index 235b678a..fa7668d7 100644 --- a/atomica/version.py +++ b/atomica/version.py @@ -3,9 +3,5 @@ Standard location for module version number and date. """ - -import sciris as sc - -version = "1.28.5" -versiondate = "2024-06-28" -gitinfo = sc.gitinfo(__file__) +version = "1.31.5" +versiondate = "2026-05-26" diff --git a/atomica/yaml_calibration.py b/atomica/yaml_calibration.py new file mode 100644 index 00000000..25ebd4ef --- /dev/null +++ b/atomica/yaml_calibration.py @@ -0,0 +1,664 @@ +""" +This file implements YAML calibration, a mechanism of programming multiple automated calibration steps +(and related operations). Essentially this scheme coordinates repeated calls to at.calibrate with different +adjustables and measurables in a pre-defined sequence of automated calibration steps. +""" + +import sciris as sc +from pathlib import Path +import shutil +import atomica as at +import numpy as np +import yaml +import time +import re + +__all__ = ["build", "run"] + +from atomica import ParameterSet + + +def _get_named_nodes(): + """ + Return dictionary with all named Node subclasses + """ + return {x._name: x for x in BaseNode.__subclasses__() if x._name is not None} + + +def build(instructions=None, context=None, name="calibration"): + """ + Construct nodes representing a calibration + + :param instructions: A dictionary of attributes/settings defined for this node OR a string filename + containing a YAML file that can be loaded to provide instructions + :param context: A dictionary of attributes/settings inherited from parent nodes + :param name: The name to assign this node + :param fname: Optionally read the instructions from a file + :return: A Node subclass instance, the type of which depends on the instructions + """ + + if (sc.isstring(instructions) or isinstance(instructions, Path)) and Path(instructions).exists(): + with open(instructions) as file: + instructions = yaml.load(file, Loader=yaml.FullLoader) + + named_nodes = _get_named_nodes() + if isinstance(instructions, dict) and ("adjustables" in instructions or (context is not None and "adjustables" in context)) and ("measurables" in instructions or (context is not None and "measurables" in context)): + return CalibrationNode(instructions, context, name) + elif name in named_nodes: + return named_nodes[name](instructions, context, name) + else: + return SectionNode(instructions, context, name) + + +def run(node, project, parset, savedir=None, save_intermediate=False, log_output: bool = False, *args, **kwargs): + """ + Run YAML calibration + + This will execute the YAML calibration using the passed-in node (or instructions to build a node), and any associated children + + :param node: Calibration node to execute. If not a node (i.e., a YAML file, or node instructions), it will be converted into a node + :param P: Project to which to apply these instructions + :param parset: An `at.ParameterSet` instance to calibrate + :param savedir: Optionally specify a directory to save the results. Defaults to the current working directory + :param save_intermediate: Set whether to save intermediate calibrations (defaults to False) + :return new_parset: A calibrated `at.ParameterSet` instance + """ + + if savedir is None: + savedir = Path(".") + else: + savedir = Path(savedir) + savedir.mkdir(exist_ok=True, parents=True) + + if not isinstance(node, BaseNode): + # Save a copy of the yaml-file if saving log output + if isinstance(node, Path) and log_output: + shutil.copyfile(node, savedir / node.name) + node = build(node) + + parset = sc.dcp(project.parset(parset)) + + nodes = list(node.walk()) # Make a flat list of all nodes to execute in order + n_steps = len([x for x in nodes if not isinstance(x[1], SectionNode)]) + n = 1 + + if log_output: + at.start_logging(savedir / "calibration_log.txt") + + at.logger.info(f"\nStarting calibration ({n_steps} steps)") + + for n_reps, node in nodes: + + if isinstance(node, SectionNode): + at.logger.info(f'\nSection "{node.name}" (repeat {n_reps} of {node.repeats})') + else: + at.logger.info(f'\nStep {n} of {n_steps}: "{node.name}" (repeat {n_reps} of {node.repeats})') + parset = node.apply(project, parset, savedir, save_intermediate, *args, **kwargs) + n += 1 + + if save_intermediate and not isinstance(node, SaveCalibrationNode): + output = savedir / f'intermediate_calibration_{n:0{len(str(n_steps))}}_{node.name.replace(" ", "_")}' + at.logger.info(f"Saving intermediate calibration...") + parset.save_calibration(output) + + t = time.process_time() + at.logger.info(f"\nCalibration completed. Total time elapsed: {round(t, 2)} seconds ({round(t/60, 2)} minutes)") + + if log_output: + at.stop_logging() + + return parset + + +class BaseNode: + """ + Node base class + + The base node class implements basic node features. Typically there should not be any + instances of this class, only instances of subclasses + """ + + _name = None # If specified, this key can be used as the name of the step to create a node of this type + + def __init__(self, instructions, context, name): + self.name = name + self.instructions = sc.dcp(instructions) + self.context = context # Attributes inherited from parent nodes + self.children = [] + self.validate() + + def walk(self): + n_reps = 0 + for repeat in range(self.repeats): + n_reps += 1 + yield (n_reps, self) + for child in self.children: + yield from child.walk() + + @property + def n_steps(self): + if type(self) == BaseNode: + return self.repeats * sum(child.n_steps for child in self.children) + else: + return self.repeats + + def __repr__(self): + return f'<{self.__class__.__name__} "{self.name}" x{self.repeats}>' + + def __str__(self, indent=0): + """ + Print a tree representation of this node and all children + + :param indent: Recursively increase the indent for child nodes + :return: + """ + s = "\t" * indent + self.__repr__() + for child in self.children: + s += "\n" + child.__str__(indent=indent + 1) + return s + + @property + def attributes(self): + return sc.mergedicts(self.context, self.instructions) + + def __getitem__(self, item): + # Directly index the Node to extract attributes without merging the dictionaries every time + if item in self.instructions: + return self.instructions[item] + elif item in self.context: + return self.context[item] + else: + raise KeyError(item) + + def __setitem__(self, key, value): + self.instructions[key] = value + + def __contains__(self, item): + return item in self.instructions or item in self.context + + @property + def repeats(self): + # Although repeats may be part of the context, we only repeat a node if the instructions requested a repeat + # i.e., repeats are not inherited + if isinstance(self.instructions, dict) and "repeats" in self.instructions: + return self.instructions["repeats"] + else: + return 1 + + def validate(self): + """ + Validate/sanitize contents of this node + + If the node isn't valid, an error should be raised + """ + return + + def apply(self, project: at.Project, parset: at.ParameterSet, savedir, *args, **kwargs) -> ParameterSet: + """ + Perform the action associated with this node + """ + return parset + + +class SectionNode(BaseNode): + """ + A section node is a special kind of node, that contains other nodes + """ + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.children = self.make_children() + self.validate() + + def make_children(self): + + children = [] + + # Remove any keys in the instructions that correspond to named nodes + # These should be used as instructions for the child node, rather than + # forming part of the context that is passed down to all children + named_nodes = _get_named_nodes() + step_instructions = {k: v for k, v in self.instructions.items() if isinstance(v, dict) or k in named_nodes} + for k in step_instructions: + del self.instructions[k] + + # Create the child nodes + for name, instructions in step_instructions.items(): + children.append(build(instructions, self.attributes, name)) + + return children + + +class CalibrationNode(BaseNode): + + # Order for list of adjustable parameters and default values + adj_defaults = { + "lower_bound": 0.1, + "upper_bound": 10.0, + "starting_y_factor": None, + } + + # Order for list of measurable parameters and default values + meas_defaults = { + "weight": 1.0, + "metric": "fractional", + "cal_start": -np.inf, + "cal_end": np.inf, + } + + @staticmethod + def parse_list(l, defaults): + # Routine to parse list of arguments into a dictionary of values + d = {} + # convert number strings back to numerical values + for i, e in enumerate(l.copy()): + try: + l[i] = float(e) + except ValueError: + pass + + for k, v in zip(list(defaults.keys())[: len(l)], l): + d[k] = v + return d + + def validate(self): + """ + Pre-parse calibration inputs + """ + + def separate_keys(keys_str: str) -> list: + """ + Separate inputs that kave been defined together as one key in the YAML file but actually represent multiple + parameters. + :param str keys_str: Unprocessed input key from the YAML file. + :return list : A list of strings, each of which represents a single key. + """ + in_brackets = False + brackets_str = "" + nobrackets_str = "" + separated_keys = [] + + for ch in keys_str: + if ch == "(": + in_brackets = True + continue + elif ch == ")": + in_brackets = False + separated_keys.append(brackets_str) + brackets_str = "" + continue + + if in_brackets: + brackets_str += ch + else: + if ch == ",": + if nobrackets_str == " " or nobrackets_str == "": + nobrackets_str = "" + continue + else: + separated_keys.append(nobrackets_str) + nobrackets_str = "" + else: + nobrackets_str += ch + + if nobrackets_str != "" and nobrackets_str != " ": + separated_keys.append(nobrackets_str) + + return [x.strip() for x in separated_keys] + + def process_key(key: str) -> tuple: + """ + Sanitize the key name, separating the parameter codename from the optional population name/s. + :par str key: Key representing an adjustable or measurable parameter. It can also contain one or two + population names (two in the case of a transfer), separated from the parameter name by a comma. + :returns: A tuple of the parameter codename and population name/s. Population defaults to None + if not specified. + + EXAMPLES: + INPUT: 'b_rate' + OUTPUT: (b_rate, None) + + INPUT: 'b_rate, 0-4' + OUTPUT: (b_rate, 0-4) + + INPUT: 'aging, 0-4, 5-14' + OUTPUT: (aging, 0-4, 5-14) + """ + if "," in key: + return tuple([x.strip() for x in key.strip("() ").split(",") if x]) + else: + return (key.strip(), None) + + def process_list(l: list) -> (tuple, list): + """ + Process list-format inputs and separate them into a key (par_name, pop_name tuple) and a value (list of + parameter settings). + @param l: List representing the settings for one parameter. + @return: tuple key: Tuple of the parameter codename and the population (default None). + list value: List of calibration settings for this parameter. + + EXAMPLES: + INPUT: [b_rate, 0.1, 10] + OUTPUT: (b_rate, None) ; [0.1, 10] + + INPUT: [(b_rate, 0-4), 0.1, 10, 1.5] + OUTPUT: (b_rate, 0-4) ; [0.1, 10, 1.5] + """ + if len(l) == 1: + # if the list is already just one string, return that string as key with None pop and vals + return (l[0].strip("() "), None), None + elif "(" in str(l): + # separate out the parenthesis contents as the par/pop/s, + # then output the key (par, pop tuple) and value + + # process keys + s = str(l).strip("[] ").replace("'", "") + s1 = re.findall(r"\(.*?\)", s) + key = process_key(s1[0].replace("(", "").replace(")", "")) + + # process values/settings + value = s.replace(s1[0], "").strip(", ").split(",") + value = [x.strip(", ") for x in value if x] + return key, value + else: + key = process_key(l[0]) + value = l[1:] + return key, value + + def process_inputs(inputs, defaults: dict) -> dict: + """ + Process adjustables and measurables, which can be specified as a string, list or nested dict representation. + * In string representation, only the parameter name is specified, and the default settings are used. + * In list representation, the input is a list of lists, where the first item in each list is the parameter + (with optional population) and the remaining items are the supported arguments for the input type, in the + order defined by the defaults dictionary. + *In dict representation, the key is the quantity with optional population, and the value can either be + a list (in the order defined by the dictionary) or a dictionary explicitly naming the inputs. + + This function returns a flat dictionary with {(quantity, pop_name):{argument:value}} e.g., {('b_rate','0-5'):{'lower_bound':0.5}}. + In the dict representation, the key can be a comma separated list of quantities with optional values e.g., + 'b_rate 0-5, d_rate'. + In the list representation, multiple quantities are not supported (as a comma is already used to separate + the arguments), but multiple lists (one for each quantity) can be provided. + #TODO could add examples + """ + out = {} + + if sc.isstring(inputs): + # Support a comma separated string with "quantity pop" specifications of adjustables and measurables + # In this case, default values should be used for all other items. Proceed by splitting into a list + inputs = inputs.split(",") + + if isinstance(inputs, (tuple, list)): + for l in inputs: + l = sc.promotetolist(l) + keyspops, v = process_list(l) + + # process key + if len(keyspops) == 2: + key, pop_name = keyspops + else: + assert len(keyspops) == 3, f"Number of populations must be 0, 1 or 2." + key = f"{keyspops[0]}_from_{keyspops[1]}" + pop_name = keyspops[2] + + # process value + if v is None: + value = defaults + else: + value = self.parse_list(v, defaults) + + out[key, pop_name] = sc.mergedicts(out.get((key, pop_name), {}), value) + + elif isinstance(inputs, dict): + for keys, v in inputs.items(): + + separated_keys = separate_keys(keys) + for key in separated_keys: + # separate par name from pop name + keyspops = process_key(key.strip()) + if len(keyspops) == 2: + key, pop_name = keyspops + else: + assert len(keyspops) == 3, f"Number of populations must be 0, 1 or 2." + key = f"{keyspops[0]}_from_{keyspops[1]}" + pop_name = keyspops[2] + + # process values + if isinstance(v, (tuple, list)): + value = self.parse_list(v, defaults) + elif v is None: + value = defaults + else: + value = v.copy() + + # add keys and values to outputs dict + out[key, pop_name] = sc.mergedicts(out.get((key, pop_name), {}), value) + return out + + self["adjustables"] = process_inputs(self["adjustables"], self.adj_defaults) + self["measurables"] = process_inputs(self["measurables"], self.meas_defaults) + + def check_optional_number(key, v, defaults): + if key in v and v[key] is not None: + if not sc.isnumber(v[key], isnan=False): + raise TypeError(f"Adjustable argument {key} needs to be a number or None (defaults to {defaults[key]}). Provided value: {v[key]} ") + + # Validate adjustables + assert len(self["adjustables"]) > 0, f"Cannot calibrate with no adjustables for calibration section {self.name}" + for (quantity, pop_name), v in self["adjustables"].items(): + assert "pop_name" not in v, f'Setting the population name through "pop_name: {v["pop_name"]}" is not supported. Instead, the name of the adjustable quantity should include the population name ("{quantity} {v["pop_name"]}")' + assert isinstance(quantity, str), f"Adjustable codename {quantity} needs to be a string" + assert pop_name is None or isinstance(pop_name, str), f"Adjustable population {pop_name} needs to be a string or None (defaults to all populations for that parameter)" + check_optional_number("lower_bound", v, self.adj_defaults) + check_optional_number("upper_bound", v, self.adj_defaults) + check_optional_number("starting_y_factor", v, self.adj_defaults) + + # Validate measurables + assert len(self["measurables"]) > 0, f"Cannot calibrate with no measurables for calibration section {self.name}" + for (quantity, pop_name), v in self["measurables"].items(): + assert isinstance(quantity, str), f"Measurable codename {quantity} needs to be a string" + assert pop_name is None or isinstance(pop_name, str), f"Adjustable population {pop_name} needs to be a string or None (defaults to all populations for that parameter)" + assert "metric" not in v or v["metric"] is None or isinstance(v["metric"], str), f"Measurable metric {v['metric']} needs to be a number or None (defaults to 'fractional')" + + check_optional_number("weight", v, self.meas_defaults) + check_optional_number("cal_start", v, self.meas_defaults) + check_optional_number("cal_end", v, self.meas_defaults) + + def apply(self, project: at.Project, parset: at.ParameterSet, n: int, *args, quiet=False, compare_results=False, **kwargs) -> ParameterSet: + + step_name = self.name + attributes = self.attributes + + at.logger.info(f"Calibrating adjustable(s) {set([adj[0] for adj in attributes['adjustables']])} to match measurable(s) {set([mea[0] for mea in attributes['measurables']])}...") + + # Expand adjustables + adjustables = {} + par_names = {x[0] for x in attributes["adjustables"]}.intersection(x.name for x in parset.all_pars()) + pop_names = {x[1] for x in attributes["adjustables"]}.intersection({*parset.pop_names} | {"all", None}) + + adj_defaults = {k: self.attributes[k] if k in self.attributes else self.adj_defaults[k] for k in self.adj_defaults} + + for par_name, pop_name in attributes["adjustables"]: + + if par_name not in par_names: + at.logger.warning(f"Extra YAML adjustable parameter '{par_name}' does not exist in this project's framework/parset and will be ignored") + continue + elif pop_name not in pop_names: + at.logger.warning(f"Extra YAML adjustable population '{pop_name}' does not exist in this project's databook and will be ignored") + continue + + if pop_name is None: + pops = parset.pop_names + else: + pops = sc.promotetolist(pop_name) + + for pop in pops: + d = sc.mergedicts(adj_defaults, attributes["adjustables"].get((par_name, None), None), attributes["adjustables"].get((par_name, pop), None)) + adjustables[(par_name, pop)] = (d["lower_bound"], d["upper_bound"], d["starting_y_factor"]) + adjustables = [(*k, *v) for k, v in adjustables.items()] + + # Expand measurables + measurables = {} + par_names = {x[0] for x in attributes["measurables"]}.intersection(x.name for x in parset.all_pars()) # TODO: This is probably OK for now but will need to validate that pars have databook entries in the future + pop_names = {x[1] for x in attributes["measurables"]}.intersection({*parset.pop_names} | {None}) + + meas_defaults = {k: self.attributes[k] if k in self.attributes else self.meas_defaults[k] for k in self.meas_defaults} + + for par_name, pop_name in attributes["measurables"]: + + if par_name not in par_names: + at.logger.warning(f"Extra YAML measurable variable '{par_name}' does not exist in this project's framework and will be ignored") + continue + elif pop_name not in pop_names: + if not (pop_name.lower() == "total" and pop_name.lower() in {x.lower() for x in project.data.tdve[par_name].ts.keys()}): + at.logger.warning(f"Extra YAML measurable population '{pop_name}' does not exist in this project's databook and will be ignored") + continue + + if pop_name is None: + pops = parset.pop_names + else: + pops = sc.promotetolist(pop_name) + + for pop in pops: + d = sc.mergedicts(meas_defaults, attributes["measurables"].get((par_name, None), None), attributes["measurables"].get((par_name, pop), None)) + measurables[(par_name, pop)] = (d["weight"], d["metric"], d["cal_start"], d["cal_end"]) + measurables = [(*k, *v) for k, v in measurables.items()] + + # Calibration + if len(adjustables): + # note: attributes = instructions + context + kwargs = sc.mergedicts(self.attributes, kwargs) + + del kwargs["adjustables"] # supplied via the adjustables variable + del kwargs["measurables"] # supplied via the measurables variable + + if "repeats" in kwargs: + del kwargs["repeats"] + + if quiet: + with at.Quiet(show_warnings=False): + new_cal_parset = at.calibrate(project, parset, adjustables, measurables, **kwargs) + else: + new_cal_parset = at.calibrate(project, parset, adjustables, measurables, **kwargs) + else: + new_cal_parset = parset + + at.logger.info(f'Completed "{step_name}"...') + made_changes = False + + for par, pop, *_ in adjustables: + + if pop == "all": + old = parset.get_par(par).meta_y_factor + new = new_cal_parset.get_par(par).meta_y_factor + else: + old = parset.get_par(par).y_factor[pop] + new = new_cal_parset.get_par(par).y_factor[pop] + + if new != old: + at.logger.info(f"...adjusted the y-factor for {par} in {pop} from {old} to {new}") + made_changes = True + else: + at.logger.debug(f"...did NOT adjust the y-factor for {par} in {pop} from {old} to {new}") + + if not made_changes: + at.logger.info(f"...made no changes!") + + if compare_results: + base_res = project.run_sim(parset=parset) + cal_res = project.run_sim(parset=new_cal_parset) + for par_name in [par_measure[0] for par_measure in measurables]: + base_rms_error = 0 + cal_rms_error = 0 + for pop in parset.pars[par_name].ts.keys(): + for time_par_ind, time_value in enumerate(parset.pars[par_name].ts[pop].t): + data_time_val = parset.pars[par_name].ts[pop].vals[time_par_ind] + base_res_time_ind = list(base_res.get_variable(par_name, pop)[0].t).index(time_value) + base_time_val = base_res.get_variable(par_name, pop)[0].vals[base_res_time_ind] + cal_res_time_ind = list(cal_res.get_variable(par_name, pop)[0].t).index(time_value) # probably redundant as they *should* be the same + cal_time_val = cal_res.get_variable(par_name, pop)[0].vals[cal_res_time_ind] + + base_rms_error += (data_time_val - base_time_val) ** 2 + cal_rms_error += (data_time_val - cal_time_val) ** 2 + + sf = at.get_sigfigs_necessary(base_time_val, cal_time_val) + at.logger.info(f"...for parameter {par_name} and population {pop} at time {time_value} the data value was {sc.sigfig(data_time_val, sf)}, the baseline value was {sc.sigfig(base_time_val, sf)}, and the calibrated value was {sc.sigfig(cal_time_val, sf)}.") + + base_rms_error = base_rms_error**0.5 + cal_rms_error = cal_rms_error**0.5 + sf = at.get_sigfigs_necessary(base_rms_error, cal_rms_error) + at.logger.info(f"...RMS error for parameter {par_name} has changed from baseline {sc.sigfig(base_rms_error, sf)} to calibrated {sc.sigfig(cal_rms_error, sf)}") + + return new_cal_parset + + +class InitializationNode(BaseNode): + _name = "set_initialization" + + def __init__(self, instructions, context, name): + new_instructions = {} + + if isinstance(instructions, dict): + new_instructions.update(instructions) + elif type(instructions) is int: + new_instructions.update({"init_year": instructions}) + elif isinstance(instructions, (tuple, list)): + sc.promotetolist(instructions) + new_instructions.update({"init_year": instructions[0]}) + if len(instructions) > 1: + new_instructions.update({"constant_parset": instructions[1]}) + + super().__init__(new_instructions, context, name) + + def validate(self): + assert "init_year" in self, f"Initialization year must be specified" + assert sc.isnumber(self["init_year"]), f'Initialization year {self["init_year"]} must be numeric.' + if "constant_parset" in self: + assert isinstance(self["constant_parset"], int), f'Constant parset (optional) {self["constant_parset"]} must be numeric (boolean, or int to specify constant parset year).' + + def apply(self, project: at.Project, parset: at.ParameterSet, n: int, *args, **kwargs) -> ParameterSet: + p2 = sc.dcp(parset) + if "constant_parset" in self: + if self["constant_parset"] == False: + pass + elif self["constant_parset"] == True: + p2 = parset.make_constant(year=project.settings.sim_start) + elif sc.isnumber(self["constant_parset"]): # constant parset year was provided + p2 = parset.make_constant(year=self["constant_parset"]) + new_settings = sc.dcp(project.settings) + new_settings.update_time_vector(end=self["init_year"]) + res = at.run_model(settings=new_settings, framework=project.framework, parset=p2) + parset.set_initialization(res, self["init_year"]) + return parset + + +class ClearInitializationNode(BaseNode): + _name = "clear_initialization" + + def __init__(self, instructions, context, name): + super().__init__(instructions=None, context=context, name=name) + + def apply(self, project: at.Project, parset: at.ParameterSet, n: int, *args, **kwargs) -> ParameterSet: + parset.initialization = None + return parset + + +class SaveCalibrationNode(BaseNode): + """ + Block in YAML file with "save calibration: " + """ + + _name = "save_calibration" + + def __init__(self, instructions, context, name): + if not isinstance(instructions, dict): + instructions = {"fname": instructions} + super().__init__(instructions, context, name) + + def validate(self): + assert self["fname"] is not None, 'A "save calibration" node must have a file name explicitly specified' + + def apply(self, project: at.Project, parset: at.ParameterSet, savedir=None, *args, **kwargs) -> ParameterSet: + parset.save_calibration(savedir / self["fname"]) + return parset diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 50fc1d9c..ed2b192b 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -3,14 +3,14 @@ schedules: displayName: Nightly build branches: include: - - master + - main - develop jobs: - job: 'tox' variables: - python.version: '3.12' - TOXENV: 'py312' + python.version: '3.14' + TOXENV: 'py314' strategy: matrix: linux: @@ -39,34 +39,12 @@ jobs: testRunTitle: 'Publish test results for Python $(python.version)' condition: succeededOrFailed() - - task: PublishCodeCoverageResults@1 + - task: PublishCodeCoverageResults@2 inputs: codeCoverageTool: Cobertura summaryFileLocation: '$(System.DefaultWorkingDirectory)/**/coverage.xml' # reportDirectory: '$(System.DefaultWorkingDirectory)/**/htmlcov' -- job: 'flake8' - pool: - vmImage: 'ubuntu-latest' - steps: - - task: UsePythonVersion@0 - displayName: 'Select python version' - inputs: - 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 - # Note - using --exit-zero will cause flake8 'errors' to show up as 'errors' in Azure Pipelines - # but the test will have been considered to pass. Therefore, the overall tests suite will still - # be considered to be passing as long as tox runs, and flake8 errors are just guidance for style - displayName: 'Run flake8' - - script: flake8_junit flake8.txt flake8_junit.xml - displayName: 'Generate flake8 junit output' - - task: PublishTestResults@2 - inputs: - testResultsFiles: '**/flake8_junit.xml' - testRunTitle: 'Publish flake8 results' - - job: 'docs' condition: and(ne(variables['Build.Reason'], 'PullRequest'), ne(variables['Build.Reason'], 'Schedule')) # Skip if it's a PR build or a scheduled build pool: @@ -106,7 +84,7 @@ jobs: - job: 'deploy' dependsOn: 'tox' - condition: and(ne(variables['Build.Reason'], 'Schedule'),and(succeeded(), eq(variables['Build.SourceBranch'], 'refs/heads/master'))) + condition: and(ne(variables['Build.Reason'], 'Schedule'),and(succeeded(), eq(variables['Build.SourceBranch'], 'refs/heads/main'))) pool: vmImage: 'ubuntu-latest' steps: @@ -115,11 +93,9 @@ jobs: inputs: versionSpec: '3.12' - script: | - python -m pip install wheel + python -m pip install build python -m pip install twine - python setup.py build - python setup.py sdist - python setup.py bdist_wheel + python -m build displayName: 'Build Atomica' - task: TwineAuthenticate@1 displayName: 'Twine Authenticate' @@ -127,4 +103,5 @@ jobs: pythonUploadServiceConnection: pypi_atomica - script: | python -m twine upload -r "pypi_atomica" --config-file $(PYPIRC_PATH) dist/* - displayName: 'Upload to PyPI' \ No newline at end of file + displayName: 'Upload to PyPI' + diff --git a/docs/conf.py b/docs/conf.py index edda0d44..699d1c5b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -24,7 +24,7 @@ # -- Project information ----------------------------------------------------- project = "Atomica" -copyright = "2023, Atomica Team" +copyright = "2025, Atomica Team" author = "Atomica Team" import atomica diff --git a/docs/examples/Plotting.ipynb b/docs/examples/Plotting.ipynb index 7c8d3180..dd4c50b0 100644 --- a/docs/examples/Plotting.ipynb +++ b/docs/examples/Plotting.ipynb @@ -335,7 +335,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Some aggregations do not make sense. For instance, if you aggregate prevalances, they should be averaged, not summed. In cases where the code is able to identify this, a warning will be displayed:" + "Some aggregations do not make sense. For instance, if you aggregate prevalances, a weighted average should be performed. In cases where the code is able to identify this, a warning will be displayed:" ] }, { @@ -348,11 +348,21 @@ "figs = at.plot_series(d)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(par_results,pop_aggregation='average',pops='total',outputs=['ac_prev'])\n", + "figs = at.plot_series(d)" + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "If you do not specify an aggregation method, the aggregation method will automatically be selected based on the units of the quantity being plotted. Non-dimensional quantities, proportions, fractions, and prevalances will be averaged instead of summed. " + "If you do not specify an aggregation method, the aggregation method will automatically be selected based on the units of the quantity being plotted. Non-dimensional quantities, proportions, fractions, and prevalances will use a population-weighted average instead of summation." ] }, { @@ -365,6 +375,57 @@ "figs = at.plot_series(d)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Characteristic aggregations\n", + "\n", + "The default aggregation for characteristics with denominators incorporates information about the denominators to select the aggregation method. For output aggregations involving characteristics, if the characteristics share a common denominator, they will be summed. This is useful for summing characteristics like prevalance of mutually exclusive strains e.g., SP and SN prevalence for TB:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(par_results,pops='total',outputs=['sp_prev',\n", + " 'sn_prev',\n", + " 'ac_prev',\n", + " {'Aggregated':['sp_prev','sn_prev']}])\n", + "figs = at.plot_series(d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note how the use of summation means that aggregating the SP and SN prevalences returns a value equal to the framework's definition of total prevalence. Alternatively, if aggregating characteristics with different denominators e.g., if there was a characteristic for diagnosed infections as a proportion of all infections, for SP and SN TB separately, a weighted average will be used to combine the characteristics. If neither of these statements are true, a raw average will be used, and a warning will be displayed. If aggregating characteristics that have different denominators, these should generally be no overlap in the compartments for those denominators in order for the weighted average to make sense. If any common compartments are detected, a warning will be displayed. \n", + "\n", + "For population aggregations of characteristics with denominators, the characteristic denominator will be used to perform weighted averages, rather than the total population size. This is useful when the denominator is something other than the population size. For example, the diagnosis sufficiency is given by people diagnosis divided by all infections. We can calculate this at the population level explicitly by summing the number of people diagnosed and the number of known infections and dividing the totals, and also by using aggregation of a characteristic that uses the number of infections as a denominator. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d1 = at.PlotData(par_results, 'known_inf',pops='total')\n", + "d2 = at.PlotData(par_results, 'ac_inf' ,pops='total')\n", + "d = at.PlotData(par_results, 'diag_suff' ,pops='total', pop_aggregation='weighted')\n", + "figs = at.plot_series(d)\n", + "figs[0].axes[0].plot(d1.series[0].tvec, d1.series[0].vals/d2.series[0].vals, 'r--');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice how the `'weighted'` average method produces the expected result, using the characteristic's denominator (`ac_inf`) as the weighting factor rather than the total population size. " + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -429,9 +490,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This is obviously incorrect because the value is greater than 1. Instead, we need to specify that we are averaging across populations. However, note that because the aggregation is an arbitrary function, the units of the output quantity are unknown. Therefore, a warning is not displayed (because the code does not know what the units are).\n", - "\n", - "Because the populations have very different sizes, a popsize-weighted average would be most appropriate:" + "This is obviously incorrect because the value is greater than 1. Because the aggregation is an arbitrary function, the units of the output quantity are unknown, and Atomica defaults to using summation for the aggregation. We need to instead explicitly specify the aggregation method. As the populations have very different sizes, a popsize-weighted average would be most appropriate." ] }, { @@ -444,6 +503,15 @@ "figs = at.plot_series(d)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "In general, averaging across populations should use the 'weighted' method rather than 'average'. For quantities that use standard units recognized by Atomica, the 'weighted' method will be selected by default. \n", + "
" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -486,17 +554,6 @@ "figs = at.plot_series(d)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Characteristic prevalance units (dimensionless)\n", - "d = at.PlotData(par_results,pops='total',output_aggregation='average',pop_aggregation='average',outputs=[{'Test':['lt_prev','spd_prev']}])\n", - "figs = at.plot_series(d)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1536,19 +1593,29 @@ "- `np.arange(2005,2040,5) = array([2005, 2010, 2015, 2020, 2025, 2030, 2035])` so the last bin is correct\n", "- `np.arange(2005,2045,5) = array([2005, 2010, 2015, 2020, 2025, 2030, 2035, 2040])` so the final bin is out of bounds\n", "\n", - "If you simply specify the bin size, e.g. `t_bins = 5` then the middle example will automatically be used. \n", - "\n", - "
\n", - "In general, it is easiest to only specify the bin width, as long as you are happy for the bins to start from the first simulation year\n", - "
\n", - "\n", - "Finally, note that because time aggregation is implemented by `PlotData`, you can apply time aggregation to `plot_series` as well as `plot_bars` e.g. if you want to plot actual annual values." + "If you simply specify the bin size, e.g. `t_bins = 5` then the middle example will automatically be used. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData([par_results,scen_results],pops=['0-4','5-14'],outputs=['sus','vac'],t_bins=5)\n", + "d.series[0].t_labels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "
\n", + "In general, it is easiest to only specify the bin width, as long as you are happy for the bins to start from the first simulation year\n", + "
\n", + "\n", + "Finally, note that because time aggregation is implemented by `PlotData`, you can apply time aggregation to `plot_series` as well as `plot_bars` e.g. if you want to plot actual annual values.\n", + "\n", "So, to make our bar graph, we first select the results, pops, and outputs using the `PlotData` constructor, then perform any necessary time aggregation, and finally render the plot using `plotBars`:" ] }, @@ -2455,9 +2522,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:atomica311] *", + "display_name": "Python [conda env:atomica312]", "language": "python", - "name": "conda-env-atomica311-py" + "name": "conda-env-atomica312-py" }, "language_info": { "codemirror_mode": { @@ -2469,7 +2536,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.12.2" }, "toc": { "base_numbering": 1, diff --git a/docs/general/YAML_calibration.ipynb b/docs/general/YAML_calibration.ipynb new file mode 100644 index 00000000..0029dfa3 --- /dev/null +++ b/docs/general/YAML_calibration.ipynb @@ -0,0 +1,493 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# YAML calibration documentation" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "This page contains additional information about features and functionality supported by the YAML calibration system\n" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## YAML internal structure\n", + "\n", + "The native representation of the YAML calibration is as a graph of 'nodes' where each node corresponds to either an action or a section containing other nodes. Each node has associated with it a collection of attributes/settings defined for that node (`instructions`), and a collection of attributes inherited from all parent nodes (`context`). The `atomica.yaml_calibration.build()` function takes in the content of the YAML file and loads it into a tree of `at.yaml_calibration.BaseNode` instances with a root node called 'calibration' at the top level. For example, consider the following YAML file from Tutorial 7:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica as at\n", + "import pathlib\n", + "yaml_file = pathlib.Path.cwd()/'..'/'tutorial'/'T7'/'calibrations'/'T7_YAML_3_repeats.yaml'\n", + "print(open(yaml_file).read())" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "This file would be loaded into the following tree of nodes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica.yaml_calibration \n", + "nodes = at.yaml_calibration.build(yaml_file) \n", + "print(nodes)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "Section nodes have a `children` attribute that in turn contains others nodes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "nodes.children[0].children[0]" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "The nodes corresponding to actions each have their own type, with methods that implement the action performed by the node:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "nodes.children[0].children[0].children[0]" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "The context for a node consists of all settings defined in the node's parents. For example, the 'Match population sizes' node inherits the 'repeats' context from the parent 'calibration' node:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "nodes.children[0].children[0].children[0].context" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "The 'instructions' for the 'Match population sizes' node contains all of the settings that are defined within the node itself - in this case, the adjustables and measurables:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "nodes.children[0].children[0].children[0].instructions" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "The YAML file is executed by sequentially traversing the tree of nodes, and calling the `apply()` method on each node in turn. The order of execution can be obtained using the `walk()` method e.g., " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "list(nodes.walk())" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "This returns a flat list of tuples, where the first item corresponds to the number of times the node has been repeated (which is used when printing progress during execution) and the second item is the node itself. " + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## `CalibrationNode` functionality\n", + "\n", + "### Adjustables and measurables settings\n", + "\n", + "A calibration node contains adjustables and measurables. Each adjustable and measurable in turn has its own settings. Each adjustable has:\n", + "\n", + "- `adj_label` (required): Adjustable codename (can be found in the framework)\n", + "- `pop_name`: Population to calibrate (default: all populations)\n", + "- `lower_bound`: Lowest value the y-factor will be allowed to take (default: 0.1)\n", + "- `upper_bound`: Highest value the y-factor will be allowed to take (default: 1)\n", + "- `starting_y_factor`: Y-factor value the autocalibration will start from when running the optimisation algorithm (default: the adjustable's current `y_factor` in the parset)\n", + "\n", + "\n", + "Each measurable has:\n", + "\n", + "- `meas_label` (required): Measurable codename (can be found in the framework)\n", + "- `pop_name`: Population to use for calibration (default: all populations)\n", + "- `weight`: Weight for a particular population (default: weight = 1. This implies that, by default, all populations are weighted equally regardless of size. See the [section on setting weights](#measurable-weights \"Measurable weights\") for further details)\n", + "- `metric`: Metric to be used by the optimisation algorithm (default: fractional)\n", + "- `cal_start`: Starting year that the calibration will be evaluated for (default: `sim_start`)\n", + "- `cal_end`: End year that the calibration will be evaluated for (default: `sim_end`)\n", + "\n", + "
\n", + "Note that sim_start and sim_end are governed by the project settings and are not set as part of the YAML calibration routine. \n", + "
\n", + "\n", + "When creating a calibration node in the YAML file, it is possible to create the adjustables and measurables using their labels only, using the default values for all other settings. Alternatively, users can specify some or all of the settings. The YAML calibration framework for Atomica supports four ways of setting adjustables and measurables, so as to give users a high level of flexibility. These are: \n", + "\n", + "- String format\n", + "- List format \n", + "- Dictionary format\n", + "- Combined format \n", + "\n", + "In general, we can only use one notation within any particular block of adjustables or measurables (with the exception of **combined format**). However, notation does not necessarily have to be consistent between different calibration blocks, or even between the adjustables and the measurables of the same calibration block. \n", + "\n", + "Below, we will describe each notation and how to use them in a YAML calibration file." + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "#### String format\n", + "\n", + "The simplest notation is string format, where only parameter names are passed to the adjustables and measurables, like so: \n", + "\n", + " Calibration: \n", + " match population sizes:\n", + " adjustables: births, mig_rate\n", + " measurables: alive\n", + "\n", + "Multiple parameters can be provided, separated by commas. When we use string notation, the optimisation algorithm will perform an autocalibration run using the default settings for adjustables and measurables. If we want to the optimisation algorithm to use specific settings, rather than just the defaults, it is necessary to use one of the other formats." + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "#### Dictionary format\n", + "\n", + "Dictionary format allows us to explicitly set calibration settings for each adjustable and measurable. We do this by writing the setting names and their values under the relevant parameter name. Each adjustable and measurable is placed on a new line, and their respective settings are also specified on separate indented lines, like so: \n", + "\n", + " Calibration: \n", + " Match population sizes:\n", + " adjustables: \n", + " \tb_rate: \n", + " \t\tstarting_y_factor: 1.2 \n", + " \tmig_rate: \n", + " lower_bound: 0.5\n", + " \t\tupper_bound: 20\t\t\n", + " measurables: \n", + " alive:\n", + " \tcal_start: 2000\n", + " \tcal_end: 2040\n", + "\n", + "We can also specify the same settings for multiple adjustables or measurables at once, by placing them together, separated by commas: \n", + "\n", + " Calibration: \n", + " match population sizes:\n", + " adjustables: \n", + " \tbirths, mig_rate: \n", + " \t\tlower_bound: 0.5\n", + " \t\tupper_bound: 20\t\t\n", + " measurables: alive \n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "#### List format\n", + "\n", + "In list format, as the name suggests, we specify the adjustables and measurables settings in a list. It can be useful as a shorthand of dictionary form, since the labels for each setting don't need to be explicitly written. Instead, we simply write the value of each setting, following the same order as in the [Adjustables and Measurables Settings](#Adjustables-and-measurables-settings) section.\n", + "\n", + "To use list format, place the parameter name and ordered settings values in a list (that is, in square brackets, separated by commas) after the adjustables and/or measurables keyword. The general structure and order to follow are shown below. \n", + "\n", + " adjustables: [adj_label, lower_bound, upper_bound, starting_y_factor]\n", + " measurables: [meas_label, weight, metric, start_year, end_year]\n", + "\n", + "Although this might seem like a lot of information for each adjustable and measurable, it is not necessary to include each item every time we use list format – only up to the point where the last setting we want to change is. For example, if we just want to set the `b_rate` adjustable's `lower bound` to, say, 0.2, we only need to list the `par_label` and `lower_bound` values in order. Any subsequent settings will retain their default values.\n", + "\n", + " adjustables: [b_rate, 0.2]\n", + " measurables: [alive]" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "However, if we wanted to set values that are at the end of the list order, we need to explicitly specify the default values of all the previous settings. For example, to set the `starting_y_factor` to 1.2 and the `end_year` to 2020, assuming the simulation start year `sim_start` was 2000, we would write:\n", + "\n", + " adjustables: [b_rate, 0.1, 10, 1.2]\n", + " measurables: [alive, 1.0, fractional, 2000, 2020]\n", + "\n", + "We can also specify settings for multiple adjustables/measurables at once, by writing a list for each adjustable or measurable, and placing them together in a list of lists. The first example from the previous section on Dictionary format would be written as\n", + "\n", + " adjustables: [[b_rate, 0.1, 10, 1.2], [mig_rate, 0.5, 20]]\n", + " measurables: [alive]\n", + "\n", + "In this particular case, it might be most practical to use the dictionary format for the `mig_rate`, while the `b_rate` is more concise in list format. List notation may also become convoluted and hard to read if there are parameters to calibrate in the same block. In case like these, we can use the combined format instead, described below." + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "#### Combined format\n", + "\n", + "Combined format uses dictionary keys, while the values are in list form. This has two main benefits: Firstly, it separates out the parameters in a clear and organised way, which avoids ending up with a dense list of lists containing long series of numbers. Secondly, it allows us to use both the list and dictionary formats under one same adjustables or measurables block.\n", + "\n", + " calibration: \n", + " match population sizes:\n", + " adjustables: \n", + " b_rate: \t\t\n", + " \t\t\tstarting_y_factor: \t1.2 \t------> dictionary format\n", + " \t\tmig_rate: [0.5, 20]\t ------> combined format\n", + " measurables: alive \n", + "\n", + "In the above example, the `b_rate` adjustable settings are in dictionary format, while the `mig_rate` is now in combined format. When using the combined format, the list of settings is defined in the same order as in the [Adjustables and Measurables Settings](#Adjustables-and-measurables-settings) section. In other words, the order is the same as when using the list format, except we don't specify the first entry (corresponding to the parameter code name) inside the list, as it is already specified before the colon. " + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "### Calibrating populations\n", + "\n", + "The YAML calibration framework allows us to indicate specific populations to calibrate. This can be useful if we wish to calibrate some populations separately, or use different calibration settings for different populations. By default, if only the code name of the adjustable or measurable is provided, a separate copy will be created for every population. To specify the population in any format (except for string format, which does not support populations), the `par_name` and `pop_name` must be placed in a tuple, i.e. in round brackets and separated by a comma, like so: `(births, 0-4)`. Calibrating populations with spaces in the `pop_label` is supported, and follows the same syntax: `(births, 0-4 HIV+)`. The following are examples of this feature’s usage in all supported formats: \n", + "\n", + "Dictionary format:\n", + "\n", + " adjustables:\n", + " (births, 0-4), mig_rate: \n", + " lower_bound: 0.5\n", + " upper_bound: 20\n", + " measurables: \n", + " (alive, 0-4):\n", + " weight: 0.1\n", + "\n", + "List format:\n", + "\n", + " adjustables: [ [(births, 0-4), 0.5, 20], [mig_rate, 0.5, 20] ]\n", + " measurables: [(alive, 0-4), 1.0]\n", + "\n", + "Combined format:\n", + "\n", + " adjustables:\n", + " (births, 0-4), mig_rate: [0.5, 20] \n", + " measurables: \n", + " (alive, 0-4): [1.0]" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "#### Meta Y-factors\n", + "\n", + "For each parameter, the meta Y-factor is applied to all populations. To calibrate the meta Y-factor and apply the same changes to every population, set the population name to `all` using the syntax above e.g., `(births, all)` would set the meta Y-factor for the `births` parameter" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "#### Overriding population specific settings\n", + "\n", + "For each adjustable that has been created, the population-specific settings will take precedence over non-population-specific settings. Recall that if no population is specified, this is equivalent to defined adjustables and measurables for each population separately. For example:\n", + "\n", + " adjustables:\n", + " births:\n", + " lower_bound: 0.5\n", + " upper_bound: 20\n", + " (births, 0-4) \n", + " upper_bound: 10\n", + "\n", + "In this case, `births` will be adjusted in every population with a lower bound of 0.5 and upper bound of 20, _except_ for the 0-4 population, in which the lower bound will be 0.5 but the upper bound will be 10." + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "### Calibrating transfers and interactions\n", + "\n", + "In the case of transfers and interactions, there are two populations involved: the _to_ population, and the _from_ population. The approach is the same as with regular populations, except that we now have two population names in the tuple instead of one:\n", + "\n", + " adjustables:\n", + " (aging, 5-14, 15-64): \n", + " lower_bound: 0.5\n", + " \tupper_bound: 20\n", + " measurables: alive\n", + "\n", + "And in list form:\n", + "\n", + "\tadjustables: [(aging, 5-14, 15-64), 0.1, 10]\n", + " measurables: alive\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "### Calibrating to total population data\n", + "\n", + "Usually, our data is structured like below, with each parameter containing several populations, such as age groups. \n", + "\n", + "\"Example\n", + "\n", + "In the YAML file, we could write \n", + "\n", + " calibration\n", + " adjustables: b_rate\n", + " measurables: alive\n", + "\n", + "Since no population has been specified, all populations will be calibrated. \n", + "\n", + "However, for some parameters, our source data might not be broken down into populations or age groups. In that case, the above YAML file will not work, since there is no data available at the individual population level. What we can do in that situation is to add an extra row to the databook with a population called `Total` (which is a reserved keyword in Atomica), and explicitly set the population name to `Total` in the YAML file. If it was our 'alive' data that was not broken down by populations, the databook would look like so\n", + "\n", + "\"Example\n", + "\n", + "And we would adjust the previous YAML file like so.\n", + "\n", + "\tcalibration\n", + " adjustables: b_rate\n", + " measurables: (alive, Total)" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "## Measurable Weights \n", + "\n", + "In one calibration block, we can include several measurables at the same time, or multiple populations of the same measurable. But say we trusted the data from one measurable's data source more than another, or wanted to prioritise fitting a particular population – how would we indicate this to the optimisation algorithm?\n", + "\n", + "In the measurable settings, we can set weights for this purpose. The default value for the `weight` setting is 1.0 which is used for all measurables and populations, which corresponds to giving them each an equal weight, regardless of size. That might be desirable if, for example, we have a key population that is smaller than the other populations – if they were weighted proportionally to size, the small key population might be effectively ignored in the optimisation. However, there might be cases where we want to do things differently. For example, we could give a key population a higher weight, or we could weight different age bins according to their size. Another reason to use measurable weights could simply be that we trust one data source more than another. \n", + "\n", + "In the following example, we set the `0-4 HIV+` population of the `alive` measurable to have double the weight than the `0-4` population. \n", + "\n", + " calibration: \n", + " match population sizes:\n", + " adjustables: \n", + " \tb_rate: \n", + " lower_bound: 0.1\n", + " upper_bound: 10\n", + " measurables: \n", + " (alive, 0-4 HIV+): \n", + " \tweight: 2\n", + " (alive, 0-4): \n", + " \t\t\t\tweight: 1\n", + "\n", + "Note that the important factor is the proportion between the different weights, not the weight values themselves. That is to say, if we instead set the above weight values to `4` and `2` respectively, the result would be the same. \n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:atomica312]", + "language": "python", + "name": "conda-env-atomica312-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "352px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/general/databook_pops.png b/docs/general/databook_pops.png new file mode 100644 index 00000000..be026587 Binary files /dev/null and b/docs/general/databook_pops.png differ diff --git a/docs/general/databook_pops_Total.png b/docs/general/databook_pops_Total.png new file mode 100644 index 00000000..f7b3fea0 Binary files /dev/null and b/docs/general/databook_pops_Total.png differ diff --git a/docs/index.rst b/docs/index.rst index 16beee87..3397fb48 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -5,7 +5,25 @@ :align: center -Atomica is a simulation engine for compartmental models. It can be used to simulate disease epidemics, health care cascades, and many other things. +`Atomica `_ is a simulation engine for compartmental models. It can be used to simulate disease epidemics, health care cascades, and many other things. + +To install via PyPI: + +.. code-block:: python + + pip install atomica + +You can run and plot an Atomica demo with: + +.. code-block:: python + + import atomica as at + import matplotlib.pyplot as plt + + P = at.demo("sir") + d = at.PlotData(P.results[0], project=P) + figs = at.plot_series(d) + plt.show() .. toctree:: :maxdepth: 1 diff --git a/docs/tutorial/T7/T7_YAML_autocalibration.ipynb b/docs/tutorial/T7/T7_YAML_autocalibration.ipynb new file mode 100644 index 00000000..f557fd60 --- /dev/null +++ b/docs/tutorial/T7/T7_YAML_autocalibration.ipynb @@ -0,0 +1,779 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# T7 - YAML calibration\n", + "\n", + "As we saw in the earlier calibration tutorial, most autocalibrations involve multiple steps, since the optimization algorithm often gets stuck in local minima if we try to optimize too many parameters at once. We could run each of these steps ourselves through separate calls to `P.calibrate()` – however, there is significant value in being able to explicitly frame the overall calibration process as an algorithm, as this makes it easier to modify the calibration steps and apply the calibration algorithm across a collection of projects. This is implemented in Atomica through the 'YAML calibration' feature, in which the calibration steps are specified in a file which Atomica can then read and use to execute the calibration.\n", + "\n", + "## YAML files\n", + "\n", + "The calibration algorithm files used by Atomica are written in YAML. YAML is a plain-text, human-readable data serialization language used to make configuration files. Essentially, a YAML file can be read into Python variables (dictionaries, lists, strings) which in turn can be used as arguments to Python functions. Here is an example of how variables can be specified in a YAML file:\n", + "\n", + "```\n", + " foo: a string\n", + " bar: 1\n", + " baz: [a,b,c]\n", + " list:\n", + " - i\n", + " - j\n", + " nested:\n", + " x: 1\n", + " y: 2\n", + " ```\n", + "\n", + "When parsed into Python, this becomes\n", + "\n", + "```\n", + "{'foo': 'a string',\n", + " 'bar': 1,\n", + " 'baz': ['a', 'b', 'c'],\n", + " 'list': ['i', 'j'],\n", + " 'nested': {'x': 1, 'y': 2}}\n", + "```\n", + "\n", + "Using YAML files provides a simple way to define a calibration algorithm in a format that is easy to work with and that Atomica can directly execute. This can cut down the time we spend manually calibrating, or even running autocalibrations. It allows us to conduct reproducible calibration runs, and is also highly scalable, since it allows us to apply the same calibration algorithm in multiple countries or settings.\n", + "\n", + "The following tutorial outlines how to use the YAML framework that has been developed for Atomica calibration. Specifically, it will cover how to write a YAML configuration file with calibration instructions for Atomica, and how to use this file to execute a calibration. Bear in mind that YAML calibration is not intended to be a standalone tool that will perfectly calibrate any model – rather, it is one part of the calibration toolbox. It can be used to reduce the time spent on calibration by autocalibrating Atomica models to a reasonable level, but additional tweaking may be required to obtain a consistently high calibration quality across all parameters, populations and/or countries." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Basic calibration example\n", + "\n", + "In this tutorial, we will work with a simple version of a typhoid model. This model captures typhoid infections, as well as asymptomatic carriers and vaccination. Firstly, we need to create an Atomica `Project` by loading in the Framework and Databook files, just like we did in [the first Atomica tutorial](https://atomica.tools/docs/master/tutorial/T1-Defining-a-model.html). The Framework and Databook for this project can be found in the Atomica repository under `T7/assets`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica as at\n", + "F = at.ProjectFramework('assets/T7_framework.xlsx')\n", + "D = at.ProjectData.from_spreadsheet('assets/T7_databook.xlsx', framework=F)\n", + "P = at.Project(framework=F, databook=D, do_run=False)\n", + "P.settings.update_time_vector(start=2000, end=2040, dt=1/52)" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "In the example above, no calibration has been loaded, so all the calibration Y-factors are equal to 1, and the model is uncalibrated. We can run a simulation and plot it, to see what our model looks like before calibration:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "uncal = P.make_parset()\n", + "res = P.run_sim(parset=uncal, result_name='Uncalibrated')\n", + "d = at.PlotData(res, outputs=['alive','deaths', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d, axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(10,7)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "Two issues with these simulation outputs immediately stand out. First, for the variables and years in which data is available, the model output doesn't match the data very well at all. Second, there is a large sudden change in the values of several variables right at the start of the simulation. In the absence of associated changes in interventions or disease transmission, disease burden is typically much more stable over time, so we would not expect to see such sudden changes at the start of the simulation.\n", + "\n", + "Tutorial 2 covers some of the detail around how to approach calibration. Although the parameters used for calibration vary from model to model, as a general rule, calibration proceeds by\n", + "\n", + "1. Demographic calibration: Adjusting birth and background death rates to match quantities like the total population size.\n", + "2. Epidemiological calibration: Adjusting parameters such as force of infection, diagnosis rate and mortality rate, to match quantities like incidence, prevalence and deaths.\n", + "\n", + "Each of these steps involves calibrating multiple parameters, which may be adjusted sequentially and repeatedly. They may also be interspersed with adjusting the model's initialization, to help minimize sudden transients at the start of the simulation. The purpose of the YAML file is to specify the sequence of steps to follow in carrying out a calibration run, in terms of which parameters to adjust, the order in which to adjust them (or whether some should be adjusted simultaneously), and how to assess the quality of the calibration at each step of the process. In this tutorial, we will use YAML files to specify a sequence of steps to calibrate the simple typhoid model shown above, starting with a minimal example and then introducing key features provided by Atomica's YAML calibration system." + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Minimal YAML file\n", + "\n", + "To illustrate how YAML calibration works, let’s start with a simple example of what a YAML file could look like:\n", + "\n", + " calibration:\n", + " \tadjustables: b_rate, mig_rate\n", + " \tmeasurables: alive\n", + "\n", + "The above YAML file represents a single call to the autocalibration optimization function, where the `b_rate` and `mig_rate` y-factors are adjusted to match `alive`. It is equivalent to running a calibration with the following command:\n", + "\n", + " calibrated_parset = P.calibrate(parset = uncal, adjustables = ['b_rate', 'mig_rate'], measurables = 'alive')\n", + "\n", + "Running the YAML calibration is very similar to performing a standalone auto-calibration. After saving the YAML file to disk (e.g., `T7_YAML_1_minimal.yaml`), the calibration can be run using\n", + "\n", + " calibrated_parset = P.calibrate(parset = uncal, yaml='T7_YAML_1_minimal.yaml')\n", + "\n", + "As the YAML calibration framework allows us to specify the adjustables and measurables, it is not necessary to provide them to `Project.calibrate()` – simply providing the YAML file is sufficient. However, the YAML file can contain multiple calibration commands, and therefore a single call to `Project.calibrate()` with a YAML file might be equivalent to multiple explicit calibration steps. The resulting simulation after running this simple calibration is like so: \n", + "\n", + "\"Simple\n", + "\"Simple\n", + "\n", + "As you can see, these plots don't look very different to the uncalibrated simulation results. If we plot the uncalibrated simulation in the same plot as our simple calibration, we can see that there has been a slight change, but not nearly enough to be able to consider this a good calibration. Next, we will illustrate the YAML calibration features that we can use to improve on this initial result. " + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Sections\n", + "\n", + "At the most basic level, a YAML calibration file defines a sequence of individual steps, where each step incrementally modifies the calibration. The YAML file therefore defines an overall _algorithm_ for performing an automatic calibration. This algorithm is defined using two structures in the YAML file:\n", + "\n", + "- _actions_, which are associated with particular operations, like running a gradient-descent calibration step with a particular set of parameters and data. A calibration action contains `adjustables` and `measurables`. Other examples of actions are detailed below.\n", + "- _sections_, which are containers for actions. Sections can contain attributes, such as how many times to repeat the contents of the section, or they can define settings that are applied to any relevant actions within that section.\n", + "\n", + "The original YAML file above consisted of a single calibration action. If we wanted to extend the algorithm by adding a step to calibrate the death rate, we could update our YAML file to the following: \n", + "\n", + "\n", + "```\n", + "calibration:\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + " Match deaths:\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + "```\n", + "\n", + "To organize the YAML file further, we could group these into an additional section. Both of these actions affect the overall population calibration, so we could logically group them as follows:\n", + "\n", + "```\n", + "calibration: \n", + " Population calibration:\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + " Match deaths:\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + "```\n", + "\n", + "The overall structure of this YAML file is thus:\n", + "\n", + "1. A top-level section titled `calibration`, which has one sub-section (`Population calibration`)\n", + "2. A sub-section called `Population calibration`, which in turn contains two actions\n", + "3. An action called `Match population sizes`, corresponding to the original calibration step we used to adjust the birth rate and migration rate\n", + "4. An action called `Match deaths`, corresponding to the new step of adjusting the death rate\n", + "\n", + "An action is differentiated from a section in two possible ways:\n", + "\n", + "- By its contents (e.g. if it contains `adjustables` and `measurables`, then it will be interpreted as a calibration action rather than a section)\n", + "- By the title (if the name corresponds to the name of a supported operation, as described below)\n", + "\n", + "Additionally, actions can never contain any sub-sections, whereas sections always do.\n", + "Apart from the names of supported operations, sections can be freely named." + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "### Repeating a section\n", + "\n", + "The *repeats* keyword can be used to loop over any part of the calibration multiple times. We do this by writing `repeats: n` inside a particular section, where *n* is the number of times we would like to loop over that section. All subsections contained in it will also be looped over *n* times. \n", + "\n", + " calibration: \n", + " repeats: 2\n", + " Population calibration:\n", + " Match population sizes:\n", + " repeats: 2\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + " Match deaths:\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + "\n", + "In the above example, we have set `repeats: 2` inside the `calibration` section, so the entire YAML calibration will be repeated twice. Then, the `Match population sizes` section also has `repeats: 2`, so the calibration step defined in that section will also be repeated twice each time. In total, there will be four calls to the optimization algorithm to match `alive`, and two to match `deaths, so the YAML file above would be equivalent to\n", + "\n", + " parset = P.calibrate(parset = parset, adjustables = [b_rate, mig_rate], measurables=alive)\n", + " parset = P.calibrate(parset = parset, adjustables = [b_rate, mig_rate], measurables=alive)\n", + " parset = P.calibrate(parset = parset, adjustables = [d_rate], measurables=deaths)\n", + " parset = P.calibrate(parset = parset, adjustables = [b_rate, mig_rate], measurables=alive)\n", + " parset = P.calibrate(parset = parset, adjustables = [b_rate, mig_rate], measurables=alive)\n", + " parset = P.calibrate(parset = parset, adjustables = [d_rate], measurables=deaths)\n", + "\n", + "In this way, it is possible for even a very compact YAML file to correspond to a large number of individual autocalibration steps. " + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Sections vs Actions\n", + "\n", + "As shown above, sections can help us to structure the calibration in a way that is practical and intuitive. They can be used to group blocks of YAML code that are conceptually related, that we want to repeat together several times, or that we want to apply similar settings to (we will cover which settings are supported in [the corresponding section](#specifying-settings-in-outer-sections)).\n", + "\n", + "Importantly though, _sections_ do not modify the calibration itself – they are merely wrappers for the innermost blocks that actually correspond to specific actions. It is in these _actions_ that operations are performed on the calibration, such as modifying the calibration or saving it.\n", + "\n", + "
\n", + "We can tell action blocks apart from sections because action blocks contain keywords indicating what kind of block they are, and they don't contain any sub-sections.\n", + "
\n", + "\n", + "It is possible to load and inspect a YAML file in Atomica without executing it. This can help confirm that the YAML file has been parsed correctly. After loading in the YAML file, it can be printed to show a summary of the sections and actions that are present, how they are nested, and how many times they are repeated: " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica.yaml_calibration \n", + "calibration_tree = at.yaml_calibration.build('calibrations/T7_YAML_3_repeats.yaml') \n", + "print(calibration_tree)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Types of Action Blocks\n", + "\n", + "YAML calibration files can contain the following types of actions, or action blocks: \n", + "\n", + "- Calibration block\n", + "- Initialization block\n", + "- Intialization clearing block\n", + "- Saving block" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### Calibration blocks\n", + "\n", + "In all of the YAML examples shown above, we have worked with 'calibration blocks', which are the main type of action. Calibration blocks are defined by the fact that they contain the keywords `adjustables` and `measurables`. Under `adjustables`, we list the parameters for Atomica’s optimization algorithm to adjust, and under `measurables`, we list the parameters to calibrate to. Each calibration block provides instructions for one optimization run, and is equivalent to making one call to `P.calibrate()` with the same adjustables and measurables.\n", + "\n", + "#### Adjustable and measurable settings\n", + "\n", + "So far, we have only specified the names of the adjustables and measurables, with no further information – in that case, the optimization algorithm will use the default calibration settings for adjustables and measurables. For more flexibility, we can customise the settings that will be used for the optimization, which directly map to the options supported by `P.calibrate()`.\n", + "\n", + "Each adjustable has:\n", + "\n", + "- `adj_label` (required): Adjustable parameter codename (can be found in the framework)\n", + "- `pop_name`: Population to calibrate (default: all populations)\n", + "- `lower_bound`: Lowest value the y-factor will be allowed to take (default: 0.1)\n", + "- `upper_bound`: Highest value the y-factor will be allowed to take (default: 10)\n", + "- `starting_y_factor`: Y-factor value the optimization algorithm will start from (default: the adjustable's current `y_factor` in the parset)\n", + "\n", + "Each measurable has:\n", + "\n", + "- `meas_label` (required): Measurable parameter codename (can be found in the framework)\n", + "- `pop_name`: Population to use for calibration (default: all populations)\n", + "- `weight`: Weight for a particular population (default: 1). By default, all populations are weighted equally regardless of size. See [the documentation on weights](https://atomica.tools/docs/master/general/YAML_calibration.html#Measurable-Weights \"Measurable weights\") for further details.\n", + "- `metric`: Metric to be used by the optimization algorithm (default: fractional)\n", + "- `cal_start`: Starting year that the calibration will be evaluated from (default: `sim_start`)\n", + "- `cal_end`: End year until which the calibration will be evaluated (default: `sim_end`)\n", + "\n", + "Note that `sim_start` and `sim_end` are the start and end years that the simulation will run for (the simulation timespan). These are distinct from `cal_start` and `cal_end`, which specify the time period for which we want to calibrate the model, i.e. a subset of the simulation timespan. For more information, see the [section on outer settings](#specifying-settings-in-outer-sections \"Specifying settings in outer sections\"). The `cal_start` and `cal_end` years can be set independently for each measurable, so it is possible to prioritize different years for different variables or for different steps of the calibration.\n", + "To specify these adjustables and measurables settings in the YAML file, we simply write the setting names and their values under the relevant parameter name. Each adjustable and measurable is placed on a new line, and their respective settings are also specified on separate indented lines, like so:\n", + "\n", + " calibration:\n", + " Match population sizes:\n", + " adjustables: \n", + " \tb_rate: \n", + " \t\tlower_bound: 0.5\n", + " \t\tupper_bound: 20\t\t\n", + " \tmig_rate: \n", + " \t\tstarting_y_factor: 1.2 \n", + " measurables: \n", + " alive:\n", + " \tcal_start: 2000\n", + " \tcal_end: 2040\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "#### Specifying Populations\n", + "\n", + "In some cases, we may only want to set adjustables or measurables for a specific population, or we may wish to use different settings for one population compared to another. A population can optionally be specified after the parameter name, as the second element of a tuple. Thus, if we only wish to calibrate some populations, we can rewrite our previous YAML file like so:\n", + "\n", + " calibration: \n", + " Match population sizes:\n", + " adjustables: \n", + " \t(b_rate, 0-4): \n", + " \t\tlower_bound: 0.5\n", + " \t\tupper_bound: 20\t\t\n", + " \t(mig_rate, 5-14): \n", + " \t\tstarting_y_factor: 1.2 \n", + " measurables: \n", + " (alive, 0-4), (alive, 5-14):\n", + " \tcal_start: 2000\n", + " \tcal_end: 2040\n", + " \n", + "Note that we can specify the same settings for more than one adjustable/measurable at once by placing several parameter names before the colon, separated by commas – this is applicable to any set of adjustables/measurables, not only different populations of the same parameter. For example, we could write: \n", + "\n", + " calibration: \n", + " Match population sizes:\n", + " adjustables: \n", + " \t(b_rate, 0-4), mig_rate: \n", + " \t\tlower_bound: 0.5\n", + " \t\tupper_bound: 20\n", + " starting_y_factor: 1.2\t\t \n", + " measurables: \n", + " (alive, 0-4), (alive, 5-14):\n", + " \tcal_start: 2000\n", + " \tcal_end: 2040\n", + "\n", + "Finally, the same syntax we use to calibrate populations can be used to calibrate transfers and interactions, but in such cases the tuple should have three elements - the parameter name, the _from_ population and the _to_ population. For more information on how to calibrate in these cases, see the documentation on [Calibrating transfers and interactions](https://atomica.tools/docs/master/general/YAML_calibration.html#Calibrating-transfers-and-interactions \"Calibrating transfers and interactions\")." + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "### Setting initializations\n", + "\n", + "In some cases, the model may exhibit an unrealistically large transient at the start of the simulation. This can occur if the initial compartment sizes calculated by Atomica are very different to the equilibrium compartment sizes associated with the model's parameters. Two common reasons for this are:\n", + "\n", + "- The initial conditions are underdetermined, for example, if there are two strains of a disease with different levels of infectiousness, but data is only available on the total prevalence. In the absence of any other information, Atomica will automatically split up the total prevalence equally between the two strains, when in fact maybe one or the other strain may be dominant, which would change the overall transmission.\n", + "- The model parameters may not give rise to an equilibrium solution that matches the data. This could be due to data sources being mixed, combined across different years, collected in different ways or with different definitions, or because the simplified dynamics in the model simply don't capture all the processes in the real world.\n", + "\n", + "Regardless of the cause, there can often be an initial transience period at the beginning of the simulation, where we can observe abrupt spikes in some parameters until the system reaches equilibrium. \n", + "\n", + "\"Pre-initialization\n", + "\n", + "In some models, this can be treated as a 'burn-in' period, and the initial part of the simulation is simply discarded. However, there is a risk that being too far from the correct initialization can contaminate the parameter estimates during the calibration period. For example, if the model is initialized with the incorrect prevalence, the calibration applied to the force of infection in order to match the observed incidence would be impacted, which might subsequently affect the model's sensitivity to interventions that change the prevalence later on. Therefore, we wish to minimize the effect of transients on the calibration.\n", + "\n", + "We can sometimes achieve this by setting the calibration start year `cal_start` to be a few years after the simulation start year `sim_start`, so that the simulation has a few years to reach equilibrium before the calibration process itself begins. However, this extends the duration of the simulation, and it might limit the amount of data that can be used for calibration.\n", + "\n", + "An alternative approach is to override the initial values of our simulation and replace them with new initialization values. Setting an initialization works by running a normal Atomica simulation for a few years (past the initial transient), taking the compartment sizes of that future year, and setting the initial compartment sizes to those stabilized values. In cases where the model parameters are approximately constant, these future compartment sizes will be roughly equivalent to what the initial compartment sizes should be, which will prevent the appearance of an initial transient in the model.\n", + "\n", + "\"Post-initialization\n", + "\n", + "If the model parameters are not roughly constant, the equilibrium that the model converges to in the future might not correspond to the equilibrium solution for the model's initial parameter values. In that case, an initial transient can still occur. To address this, we can remove any time variation in the model's parameter values using the `ParameterSet.make_constant()` method. This will return a copy of the parset in which all parameters are constant over time, thus ensuring that the future compartment sizes are computed based on the same parameter values as the initial simulation year. This often provides a suitable solution, although changes to the total population size due to births and deaths can still take place, so in some cases a small initial transient may still be present. In such cases, repeatedly setting initialization based on a shorter simulation can help minimize the discrepancy. For more information on Atomica initializations, see the [documentation on this topic](https://atomica.tools/docs/master/general/Compartment-Initialization.html \"Initializing compartments\").\n", + "\n", + "In the YAML file, we indicate that we want to set a new initialization by making a YAML block with the title `set_initialization`. Under this title, we can specify further settings:\n", + "\n", + "- `init_year` (required): The future year to take the compartment sizes from. The simulation will be run up to this year.\n", + "- `constant_parset` (default: False): Whether to make the parameters constant while the initialization is being calculated, and which year to extract these constant parameter values from, i.e. which year to use in `parset.make_constant()`. It can be a numerical value (representing the year from which to draw the constant values for the parset) or Boolean (if `True`, the default year `sim_start` will be used).\n", + "\n", + "There are thus several valid ways to set an initialization. For example, only setting the initialization year:\n", + "\n", + "```\n", + "calibration:\n", + " set_initialization:\n", + " init_year: 2030\n", + "```\n", + "\n", + "Setting the initialization year with a constant parset, using the parameter values from the default `sim_start` year:\n", + "\n", + "```\n", + "calibration:\n", + " set_initialization:\n", + " init_year: 2030\n", + " constant_parset: True\n", + "```\n", + "\n", + "Setting the initialization year with a constant parset, using constant parameter values from a specific year:\n", + "\n", + "\n", + "```\n", + "calibration:\n", + " set_initialization:\n", + " init_year: 2030\n", + " constant_parset: 2005\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "### Clearing initializations\n", + "\n", + "If we have previously set an initialization in our YAML calibration algorithm, and we then set another initialization later in the YAML file, it uses information from the previous initialization to calculate the next, since the new simulation will start from the initial compartment sizes calculated in that previous initialization.\n", + "\n", + "Sometimes we might wish to calculate a new initialization from scratch, without using our previous initializations as a starting point. This could be useful if we have done some calibration steps between the previous initialization and now, in which case the y-factors will have changed, which in turn could have caused the equilibrium values to vary.\n", + "\n", + "To do this in the YAML file, we can make a section titled `clear_initialization`, followed by a Boolean value. If it is set to `True`, any existing initialization will be cleared; if `False`, nothing will happen. The most simple example is as follows:\n", + "\n", + "```\n", + "calibration:\n", + " clear_initialization: True\n", + "```\n", + "\n", + "Below is a more practical example of how we would use the `clear_initialization` feature in a real YAML file, where we first set an initialization, perform some calibration steps, and then clear our previous initialization before making a new one.\n", + "\n", + " calibration:\n", + " First initialization:\n", + " set_initialization:\n", + " init_year: 2030\n", + "\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " \tmeasurables: alive\n", + " \n", + " \tclear_initialization: True\n", + "\n", + " Second initialization:\n", + " set_initialization:\n", + " init_year: 2030\n", + "\n", + "\n", + "Note that we need to wrap `set_initialization` in an outer section, since YAML does not allow duplicate keys (i.e. identical section/action titles at the same level)." + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "### Saving calibrations\n", + "\n", + "Throughout a YAML calibration, we might wish to save the calibration state at specific points in the process. For example, if our YAML file has a population section and an epidemiological section, we might want to save the calibration after the population calibration section so we can see the progress made up to that point, or otherwise identify the effect that different parts of the algorithm are having on the calibration. To save a calibration, we simply make a section titled `save_calibration`. Under the title, we indicate the filename we wish to save the calibration to, either by providing it directly, or by using the keyword `fname` - both examples are shown below:\n", + "\n", + " calibration:\n", + " \tPopulation section: \n", + " \t\t[…]\n", + " \n", + " \tsave_calibration: \n", + " \t\tfname: pop_calibration.xlsx\n", + " \n", + " \tEpidemiology section: \n", + " \t\t[…]\n", + "\n", + " \tsave_calibration: epi_calibration.xlsx\n", + " \n", + "Note that when we save a calibration, if the initial compartment sizes have been explicitly specified using `set_initialization`, these compartment sizes will be saved along with the y-factors in the same Excel file. Loading the calibration from this file will thus also set the initial compartment sizes to those values.\n", + "\n", + "Another option that can be useful for debugging is to save the calibration state at every intermediate step of the YAML file. In that case, we can use the `save_intermediate_calibrations` option when loading in and running the calibration, e.g.\n", + "\n", + " calibrated_parset = P.calibrate(parset = cal, yaml='T7_YAML_1_minimal.yaml', save_intermediate=True)\n", + "\n", + "More options for running the calibration can be found [here](#running-the-yaml-calibration \"Running the YAML calibration\")." + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Specifying settings in outer sections\n", + "\n", + "We saw previously that sections (i.e. YAML blocks that don't correspond to actions) don't directly modify the calibration, and are mainly used to structure the YAML file. However, it is possible to determine calibration settings inside a section, such that they are passed down to any subsections or action blocks contained within. Another way to think of this is that an action block will inherit any settings that are defined in any of its parent sections.\n", + "\n", + "This can simplify the process of writing YAML files where we want to override the default settings in several action blocks, as it allows us to specify those settings once in a parent section, rather than repeatedly in every action block. Additionally, this feature is hierarchical, so settings that are specified further in (e.g. in the action block itself) will always override those set in a section that is further out, allowing for more flexibility.\n", + "\n", + "Some settings that this feature is commonly used for include: \n", + "\n", + "- `max_time`: Maximum amount of time each call to the optimization algorithm will run for (default: 60 seconds)\n", + "- `stepsize`: Initial step size, i.e. how much the y-factors will be increased/decreased by the optimization algorithm (default: 0.1)\n", + "\n", + "Additionally, any adjustable or measurable setting can also be set outside of its calibration block, although some settings lend themselves more to this feature than others. For example, changing the measurable `weight` for an entire calibration block has no effect, as what matters is the proportion between the weights of different measurables in the same block. The following settings can be useful to set in parent sections:\n", + "\n", + "- `cal_start` and `cal_end`: These can be used to change the calibration timespan for multiple calibration steps at the same time.\n", + "- `metric`: To change the metric used to assess calibration quality, i.e. how well the simulation fits the measurables data.\n", + "\n", + "We will now provide some examples of how this functionality can be used.\n", + "\n", + "### Inheritance\n", + "\n", + "Settings are automatically inherited by all subsections. To set a `max_time` of 120 for every calibration step in the entire calibration, we would write the following. This will result in every calibration step running for at most 120 seconds, rather than the default 60 seconds.\n", + "\n", + " calibration:\n", + " max_time: 120\n", + " Population calibration:\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + " Match deaths:\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + " \n", + " Epidemiological calibration\n", + " [...]\n", + "\n", + "### Settings hierarchy\n", + "\n", + "Inherited settings can be overwritten inside subsections (in both sections and actions). Say we wanted all calibration blocks to run for a `max_time` of 120 seconds, except for `Match_deaths`, which we want to run for only 60 seconds. In that case, we can override the setting determined in the parent section by specifying the updated value inside `Match_deaths` subsection (which happens to be an action block):\n", + "\n", + " calibration:\n", + " max_time: 120\n", + " Population calibration:\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + " Match deaths:\n", + " max_time: 60\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + " \n", + " Epidemiological calibration\n", + " [...]\n", + "\n", + "### Outer adjustables/measurables settings\n", + "\n", + "Finally, here is an example of how to use the adjustables and measurables settings outside a calibration block:\n", + "\n", + " calibration:\n", + " Population calibration:\n", + " upper_bound: 2\n", + " lower_bound: 0.5\n", + "\n", + " Match population sizes:\n", + " adjustables: b_rate, mig_rate\n", + " measurables: alive\n", + "\n", + " Match deaths:\n", + " adjustables: d_rate\n", + " measurables: deaths\n", + "\n", + " Epidemiological calibration:\n", + " [...]\n", + "\n", + "In this case, the upper and lower bounds have been updated to 2 and 0.5 respectively, for all the adjustables in the `Population calibration` section (that is, `b_rate`, `mig_rate` and `d_rate`). The upper and lower bounds in the Epidemiological section will remain unchanged, using the default values of 0.1 and 10." + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "## Running the YAML calibration\n", + "\n", + "Now that we have finished writing our YAML calibration file, we can proceed to running the calibration. Having loaded a project, `P.calibrate()` can then be called directly, passing in the uncalibrated parset and the name of the YAML file:\n", + "\n", + " calibrated_parset = P.calibrate(parset=uncal, yaml='T7_YAML_1_minimal.yaml')\n", + "\n", + "This function also supports several optional arguments:\n", + "\n", + "- `savedir`: Any saved calibrations and logs will be saved into this folder (default: current working directory).\n", + "- `save_intermediate`: If `True`, this will save a calibration file at every step (default: `False`).\n", + "- `max_time`: Maximum amount of time each calibration step will run for (default: 60 seconds).\n", + "- `log_output`: If `True`, this will save a text file into the output directory containing all console output, e.g. objective function values (default: `False`).\n", + "\n", + "The YAML file can also be supplied as a dictionary, which would normally be obtained via\n", + "\n", + " import yaml\n", + " yaml_instructions = yaml.load(file, Loader=yaml.FullLoader)\n", + " \n", + "The `yaml_instructions` variable above can be passed to `P.calibrate(..., yaml=yaml_instructions)`. This enables changes to be made to the YAML content programmatically prior to running the calibration." + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "## Exercise: Worked Example\n", + "\n", + "Now that we understand what all the working parts of a YAML file are, let's put it all together. For each question, write a YAML file to calibrate the model as described. Each question will consist of incremental additions to the previous solution.\n", + "\n", + "_**Question 1.** We want to do a basic population calibration, where we calibrate the birth rate and migration rate to match the data corresponding to the total number of people alive. We also want to calibrate the death rate. What should this YAML file look like?_\n", + "\n", + "**Hint**: Open the framework and look at the `Compartments`, `Parameters` and `Characteristics` pages. The `Code Names` and `Display Names` show us how we have to write the parameter names in the YAML file, and what quantities they correspond to. Those with a value in the `Databook Page` column have data values supplied, and can therefore be used as measurables.\n", + "\n", + "**Hint 2**: Some of the parameter `Code Names` and `Display Names` in the databook are as follows: \n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Framework SheetCode NameDisplay NameDatabook Page
CompartmentsbirthBirthnone
deathDeathnone
Parametersb_rateBirth ratedemographic
deathsAll-cause deathsdemographic
d_rateBackground mortality ratedemographic
mig_rateMigration Ratedemographic
CharacteristicsaliveTotal populationdemographic
" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "_**Question 2.** A single calibration run may not be enough to get good results, so let's loop over our simple population calibration ten times. How would we make that change to the YAML file?_\n", + "\n", + "**Hint**: Use the `repeats` keyword.\n", + "\n", + "_**Question 3.** Say that, for this particular project, we are only interested in calibrating results from 2005 to 2040. How would we specify this in the YAML file to reduce calibration time?_\n", + "\n", + "**Hint**: Using the `cal_start` and `cal_end` keywords.\n", + "\n", + "_**Question 4.** When calibrating the birth rate, it only really makes sense to calibrate the 0-4 population. Modify the YAML file to reflect this._\n", + "\n", + "**Hint**: We place the parameter code name and population name in a tuple.\n", + "\n", + "_**Question 5.** Say we know that our data source underestimates the birth rate. Let’s set the starting y_factor for this parameter to 1.2 to speed up the optimization._\n", + "\n", + "**Hint**: Use `starting_y_factor`.\n", + "\n", + "_**Question 6.** We want to avoid the presence of transients in our calibrated simulation. Let’s initialize the calibration in order to eliminate any jumps._\n", + "\n", + "**Hint**: Use `set_initialization` after the calibration blocks, and make sure the initialization gets re-calculated in every loop! \n", + "\n", + "_**Question 7.** We also want to clear the previous initialization every time we make a new one, instead of using information from the previous initialization. How can we update the YAML file to reflect this?_\n", + "\n", + "**Hint**: Use `clear_initialization`, and make sure the initialization gets cleared and re-calculated in every loop!\n", + "\n", + "_**Question 8.** We are almost ready to run the YAML calibration! Now, what instructions do we need to add to automatically save our population calibration once it is done?_\n", + "\n", + "**Hint**: Use `save_calibration` and specify a `fname`.\n", + "\n", + "All solutions to the worked example can be found in the `docs/tutorial/T7/calibrations` folder of the Atomica repository. The calibration obtained after running the last YAML file in this exercise yields the following simulation result:\n", + "\n", + "\"Worked\n", + "\n", + "The above calibration result is an improvement on the simple example shown earlier in the tutorial - while the simulation result for the 15-64 age group (in green) may not look as close to the data, the remaining populations have an excellent fit, leading to a very good overall population calibration. Additionally, all-cause deaths are now calibrated. Note that we do not expect the calibrated simulation to match the data points exactly, but to follow the same overall trends." + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "## Further resources\n", + "\n", + "### Epidemiological calibration\n", + "\n", + "In this tutorial we have demonstrated the key functionality of Atomica's YAML calibration system applied to a population calibration, focussing on adjusting births, migration and deaths to match population size data, without considering disease burden. The same functionality can be used for epidemiological calibration, just with different model parameters and data. A complete YAML file of the population and epidemiological calibration instructions can be found in the Atomica code repository, at `docs/tutorial/T7/calibrations/typ_calibration_instructions.yaml`.\n", + "\n", + "In this YAML file, we have set the `max_time` to 120, to give the calibration a bit more time to reach an optimal solution in each calibration step. We then calibrate the population following the same principles illustrated previously, repeating the population calibration ten times, and setting a new initialization at the end of each loop, with `constant_parset=True` and the `init_year` set to 2030.\n", + "\n", + "In the epidemiological calibration section, the typhoid incidence, prevalence, and typhoid deaths are calibrated. Since we don't have a lot of information on the order of magnitude of the susceptibility and infectiousness, the lower and upper bounds are expanded to leave more room for variation. However, we don't want the disease duration `typ_gen_dur` to vary a lot, since we have a pretty clear idea of its magnitude and don't want it to vary too much between calibrations and settings, so we set stricter bounds. \n", + "\n", + "At the end of the typhoid calibration, we set an initialization in the same way we did for the population calibration. We then repeat the epidemiological calibration ten times, and finally, repeat the whole YAML calibration process (except silencing the epi y-factors) twice.\n", + "\n", + "The resulting calibration from running this YAML file is like so:\n", + "\n", + "\"Complete\n", + "\n", + "As you can see, this is an improvement on the previous calibrations shown, since the epidemiological parameters (like typhoid prevalence and typhoid deaths) are now also calibrated to match the data.\n", + "\n", + "__Alternative approach__\n", + "\n", + "If the epidemiological parameters have a very large impact on the population parameters (for example, if the uncalibrated disease prevalence or disease mortality is much larger than the data would indicate), this can make it hard to properly calibrate the population in the first instance. We might observe that the population y-factors relating to births and migration increase (or, in the case of deaths, decrease) by a lot, to compensate for the large amount of people dying from disease.\n", + "\n", + "To prevent this behaviour, we can do the following: At the beginning of the YAML file, we set the parameters relating to the typhoid disease (`typ_active_inf` and `typ_car`) to zero. This essentially \"turns off\" of the typhoid disease, giving us the opportunity to reasonably calibrate the population without interference from the disease components. Then, just after the population calibration section, we \"reactivate\" the disease by setting the `typ_active_inf` and `typ_car` y-factors back to 1, before proceeding to the epidemiological calibration. Looping over this whole algorithm a couple of times (except for the first step where we \"turn off\" the disease) allows the population to be fine-tuned while taking the epidemiological parameters into account.\n", + "\n", + "### Documentation\n", + "\n", + "For more information on using Atomica's YAML calibration functionality, see [the documentation](https://atomica.tools/docs/master/general/YAML_calibration.html \"YAML documentation page\").\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.11" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "343.4px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorial/T7/assets/T7_databook.xlsx b/docs/tutorial/T7/assets/T7_databook.xlsx new file mode 100644 index 00000000..e41e519a Binary files /dev/null and b/docs/tutorial/T7/assets/T7_databook.xlsx differ diff --git a/docs/tutorial/T7/assets/T7_framework.xlsx b/docs/tutorial/T7/assets/T7_framework.xlsx new file mode 100644 index 00000000..4d819763 Binary files /dev/null and b/docs/tutorial/T7/assets/T7_framework.xlsx differ diff --git a/docs/tutorial/T7/assets/T7_plot_1_simple_compare.png b/docs/tutorial/T7/assets/T7_plot_1_simple_compare.png new file mode 100644 index 00000000..5d73a052 Binary files /dev/null and b/docs/tutorial/T7/assets/T7_plot_1_simple_compare.png differ diff --git a/docs/tutorial/T7/assets/T7_plot_2_simple.png b/docs/tutorial/T7/assets/T7_plot_2_simple.png new file mode 100644 index 00000000..a75ba620 Binary files /dev/null and b/docs/tutorial/T7/assets/T7_plot_2_simple.png differ diff --git a/docs/tutorial/T7/assets/T7_plot_3_worked_example.png b/docs/tutorial/T7/assets/T7_plot_3_worked_example.png new file mode 100644 index 00000000..5980d53c Binary files /dev/null and b/docs/tutorial/T7/assets/T7_plot_3_worked_example.png differ diff --git a/docs/tutorial/T7/assets/T7_plot_4_complete.png b/docs/tutorial/T7/assets/T7_plot_4_complete.png new file mode 100644 index 00000000..3fe308bd Binary files /dev/null and b/docs/tutorial/T7/assets/T7_plot_4_complete.png differ diff --git a/docs/tutorial/T7/assets/fig_post_init.png b/docs/tutorial/T7/assets/fig_post_init.png new file mode 100644 index 00000000..5cc4a4e7 Binary files /dev/null and b/docs/tutorial/T7/assets/fig_post_init.png differ diff --git a/docs/tutorial/T7/assets/fig_pre_init.png b/docs/tutorial/T7/assets/fig_pre_init.png new file mode 100644 index 00000000..5ce00e95 Binary files /dev/null and b/docs/tutorial/T7/assets/fig_pre_init.png differ diff --git a/docs/tutorial/T7/calibrations/Plot_YAML_tutorial_plots.ipynb b/docs/tutorial/T7/calibrations/Plot_YAML_tutorial_plots.ipynb new file mode 100644 index 00000000..b8995e53 --- /dev/null +++ b/docs/tutorial/T7/calibrations/Plot_YAML_tutorial_plots.ipynb @@ -0,0 +1,516 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Calibration plots for the YAML files shown in the YAML tutorial\n", + "\n", + "This plots the calibration reslts from `run_tutorial_autocals` (which includes the final worked example as well as the YAML examples in the tutorial)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "# Make project" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica as at\n", + "import os\n", + "from os.path import isfile\n", + "import re\n", + "import sciris as sc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "inputs = '../assets'\n", + "F = at.ProjectFramework(f'{inputs}/T7_framework.xlsx')\n", + "D = at.ProjectData.from_spreadsheet(f'{inputs}/T7_databook.xlsx', framework=F)\n", + "P = at.Project(framework=F,databook=D)\n", + "P.settings.update_time_vector(start=2000, end=2040, dt=1/52)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "# Pre-calibraton" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "cal0 = P.make_parset()\n", + "res0 = P.run_sim(parset=cal0, result_name = 'Uncalibrated')\n", + "d = at.PlotData(res0, outputs=['alive','deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(8,4)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "cal0 = P.make_parset()\n", + "res0 = P.run_sim(parset=cal0, result_name = 'Uncalibrated')\n", + "d = at.PlotData(res0, outputs=['alive', 'typ_incidence', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(8,6)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "# Import calibrations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "#calibration files to import\n", + "calpaths = ['cal_1', '', 'cal_3', 'cal_4', 'cal_5', 'cal_6', 'cal_7', 'cal_8', 'cal_9', 'cal_WE8']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "resdict = sc.odict()\n", + "for i, file in enumerate(calpaths):\n", + " print(i+1)\n", + " print(file)\n", + " print(f'res_{i+1}')\n", + " if file == '':\n", + " continue\n", + " cal = P.make_parset() \n", + " cal = cal.load_calibration(f'{file}.xlsx')\n", + " resdict[f'res_{i+1}'] = P.run_sim(parset=cal, result_name = file)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "# Plot calibrations" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## 1 - Minimal yaml file (shown in tutorial)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "# specific calibration\n", + "file = 'cal_1.xlsx'\n", + "cal1 = P.make_parset() \n", + "cal1 = cal1.load_calibration(file)\n", + "res1 = P.run_sim(parset=cal1, result_name = 'Simple calibration')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "d = at.PlotData([res0, res1], outputs='alive', project=P, pops='65+')\n", + "at.plot_series(d, axis='results', data=P.data);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(res1, outputs=['alive','typ_incidence', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(12,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(res1, outputs=['alive', 'deaths', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(12,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 3 - Repeats (not shown at the moment)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "d = at.PlotData(resdict['res_3'], outputs=['alive','deaths', 'typ_incidence', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(10,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(resdict['res_3'], outputs=['alive','typ_incidence', 'typ_prev', 'typ_num_deaths', ], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(10,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 4 - cal dict format (not shown)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_4'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 5 - calibrate specific pops (not shown)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_5'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 6 - cal with mult params (not shown)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_6'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 7 - init (showing a different example)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_7'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 8 - clear init (showing a different example)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_8'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "# cal 9 - outer settings (not shown)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "#plot total pop: cals vs no calibration\n", + "d = at.PlotData(resdict['res_9'], outputs='alive', project=P)\n", + "at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "33", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "### add -cal with outer params (not shown)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "# #plot total pop: cals vs no calibration\n", + "# d = at.PlotData(resdict['res_9'], outputs='alive', project=P)\n", + "# at.plot_series(d, axis='pops', data=P.data);" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "# WE_8: Worked example calibration (shown in tutorial)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(resdict['res_10'], outputs=['alive', 'deaths', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(12,8)\n", + "fig.tight_layout()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "AtomicaTest", + "language": "python", + "name": "atomicatest" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorial/T7/calibrations/Plot_final_typhoid_calibration.ipynb b/docs/tutorial/T7/calibrations/Plot_final_typhoid_calibration.ipynb new file mode 100644 index 00000000..6c4ca419 --- /dev/null +++ b/docs/tutorial/T7/calibrations/Plot_final_typhoid_calibration.ipynb @@ -0,0 +1,202 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Calibration plots for the pop + epi calibration\n", + "\n", + "This plots the results from `run_model_autocals.py` (whole calibration)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import atomica as at\n", + "import os\n", + "from os.path import isfile\n", + "import re\n", + "import sciris as sc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "#load typhoid project\n", + "inputs = '../assets'\n", + "F = at.ProjectFramework(f'{inputs}/T7_framework.xlsx')\n", + "D = at.ProjectData.from_spreadsheet(f'{inputs}/T7_databook.xlsx', framework=F)\n", + "P = at.Project(framework=F, databook=D, do_run=False)\n", + "P.settings.update_time_vector(start=2000, end=2040, dt=1/52)" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "# Plot uncalibrated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "cal0 = P.make_parset()\n", + "res0 = P.run_sim(parset=cal0, result_name = 'Uncalibrated')\n", + "d = at.PlotData(res0, outputs=['alive', 'deaths'], project=P)\n", + "at.plot_series(d, axis='pops', data=P.data, n_cols=2, legend_mode='none');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(res0, outputs=['alive','typ_incidence', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(12,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "# Load in calibrations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "calfile= \"typ_calibration.xlsx\"\n", + "print(calfile)\n", + "cal = P.make_parset() \n", + "cal = cal.load_calibration(f'{calfile}')\n", + "res1 = P.run_sim(parset=cal, result_name = 'Complete YAML calibration')" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "# Plot calibrations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(res1, outputs=['alive', 'deaths'], project=P)\n", + "at.plot_series(d, axis='pops', data=P.data, n_cols=2, legend_mode='none');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "d = at.PlotData(res1, outputs=['alive','deaths', 'typ_prev', 'typ_num_deaths'], project=P)\n", + "fig = at.plot_series(d,axis='pops', data=P.data, n_cols=2, legend_mode='none')[0]\n", + "fig.set_size_inches(12,8)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "#display deaths\n", + "d = at.PlotData(res1, outputs='d_rate:flow', project=P) \n", + "for s in d.series:\n", + " s.data_label = 'deaths'\n", + "at.plot_series(d,axis='pops',data=P.data);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "#display other typh vars\n", + "d = at.PlotData(res1, outputs=['life_exp'], project=P)\n", + "at.plot_series(d,axis='pops',data=P.data);" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "AtomicaTest", + "language": "python", + "name": "atomicatest" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorial/T7/calibrations/README.txt b/docs/tutorial/T7/calibrations/README.txt new file mode 100644 index 00000000..db78fa2e --- /dev/null +++ b/docs/tutorial/T7/calibrations/README.txt @@ -0,0 +1,10 @@ +Contains example YAML files and scripts to run calibrations with them + +- `T7_YAML_*` files correspond to the examples of YAML calibrations used in Tutorial 7 +- `WE_*.yaml` are the solutions to the 'worked example' section in Tutorial 7 +- `typ_calibration_instructions.yaml` corresponds to a full finished YAML calibration for the model + +This folder also contains + +- `*.py` files to run the YAML files and produce `*.xlsx` calibration outputs +- `*.ipynb` notebooks to generate plots using those calibration outputs \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/T7_YAML_1_minimal.yaml b/docs/tutorial/T7/calibrations/T7_YAML_1_minimal.yaml new file mode 100644 index 00000000..dac077d3 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_1_minimal.yaml @@ -0,0 +1,7 @@ + +calibration: + adjustables: b_rate, mig_rate + measurables: alive + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_2_sections.yaml b/docs/tutorial/T7/calibrations/T7_YAML_2_sections.yaml new file mode 100644 index 00000000..227d32d5 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_2_sections.yaml @@ -0,0 +1,12 @@ + +calibration: + + Population calibration: + + Match population sizes: + adjustables: b_rate, mig_rate + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/T7_YAML_3_repeats.yaml b/docs/tutorial/T7/calibrations/T7_YAML_3_repeats.yaml new file mode 100644 index 00000000..f5b526f0 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_3_repeats.yaml @@ -0,0 +1,11 @@ + +calibration: + repeats: 2 + Population calibration: + Match population sizes: + repeats: 5 + adjustables: b_rate, mig_rate + measurables: alive + Match deaths: + adjustables: d_rate + measurables: deaths diff --git a/docs/tutorial/T7/calibrations/T7_YAML_4_cal.yaml b/docs/tutorial/T7/calibrations/T7_YAML_4_cal.yaml new file mode 100644 index 00000000..5fca7b66 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_4_cal.yaml @@ -0,0 +1,16 @@ + +calibration: + Match population sizes: + adjustables: + b_rate: + lower_bound: 0.5 + upper_bound: 20 + mig_rate: + starting_y_factor: 1.2 + measurables: + alive: + cal_start: 2000 + cal_end: 2040 + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_5_cal_pop.yaml b/docs/tutorial/T7/calibrations/T7_YAML_5_cal_pop.yaml new file mode 100644 index 00000000..caef4bbc --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_5_cal_pop.yaml @@ -0,0 +1,17 @@ + +calibration: + Match population sizes: + adjustables: + (b_rate, 0-4): + lower_bound: 0.5 + upper_bound: 20 + (mig_rate, 5-14): + starting_y_factor: 1.2 + measurables: + (alive, 0-4), (alive, 5-14): + cal_start: 2000 + cal_end: 2040 + + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_6_cal_multiple_params.yaml b/docs/tutorial/T7/calibrations/T7_YAML_6_cal_multiple_params.yaml new file mode 100644 index 00000000..c97d9248 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_6_cal_multiple_params.yaml @@ -0,0 +1,17 @@ + +calibration: + Match population sizes: + adjustables: + b_rate, mig_rate: + lower_bound: 0.5 + upper_bound: 20 + starting_y_factor: 1.2 + measurables: + alive: + cal_start: 2000 + cal_end: 2040 + + + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_7_transience_initialisation.yaml b/docs/tutorial/T7/calibrations/T7_YAML_7_transience_initialisation.yaml new file mode 100644 index 00000000..508c0749 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_7_transience_initialisation.yaml @@ -0,0 +1,14 @@ + + + + + + set_initialization: + init_year: 2010 + constant_parset: True + + + + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_8_clearing_init.yaml b/docs/tutorial/T7/calibrations/T7_YAML_8_clearing_init.yaml new file mode 100644 index 00000000..314f0069 --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_8_clearing_init.yaml @@ -0,0 +1,14 @@ + +calibration: + clear_initialization: True + + + + + + save_calibration: + fname: Pop_calibration.xlsx + + + + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings.yaml b/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings.yaml new file mode 100644 index 00000000..b6e3622c --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings.yaml @@ -0,0 +1,15 @@ + + +calibration: + Population calibration: + upper_bound: 2 + lower_bound: 0.5 + + Match population sizes: + adjustables: b_rate, mig_rate + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths + diff --git a/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings_maxtime.yaml b/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings_maxtime.yaml new file mode 100644 index 00000000..f742fedc --- /dev/null +++ b/docs/tutorial/T7/calibrations/T7_YAML_9_outer_settings_maxtime.yaml @@ -0,0 +1,12 @@ + +calibration: + Population calibration: + max_time: 120 + Match population sizes: + adjustables: b_rate, mig_rate + measurables: alive + Match deaths: + adjustables: d_rate + measurables: deaths + + diff --git a/docs/tutorial/T7/calibrations/WE_1.yaml b/docs/tutorial/T7/calibrations/WE_1.yaml new file mode 100644 index 00000000..e93ca5e4 --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_1.yaml @@ -0,0 +1,12 @@ +calibration: + Population calibration: + + Match population sizes: + adjustables: b_rate, mig_rate + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths + + \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/WE_2.yaml b/docs/tutorial/T7/calibrations/WE_2.yaml new file mode 100644 index 00000000..00af8391 --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_2.yaml @@ -0,0 +1,14 @@ +calibration: + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: + b_rate, mig_rate: + lower_bound: 0.1 + upper_bound: 10 + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths diff --git a/docs/tutorial/T7/calibrations/WE_3.yaml b/docs/tutorial/T7/calibrations/WE_3.yaml new file mode 100644 index 00000000..6e90a45e --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_3.yaml @@ -0,0 +1,14 @@ +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: b_rate, mig_rate + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/WE_4.yaml b/docs/tutorial/T7/calibrations/WE_4.yaml new file mode 100644 index 00000000..a422925a --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_4.yaml @@ -0,0 +1,14 @@ +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: (b_rate, 0-4), mig_rate + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/WE_5.yaml b/docs/tutorial/T7/calibrations/WE_5.yaml new file mode 100644 index 00000000..da9fa653 --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_5.yaml @@ -0,0 +1,18 @@ +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: + (b_rate, 0-4): + starting_y_factor: 1.2 + mig_rate: + lower_bound: 0.1 + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths diff --git a/docs/tutorial/T7/calibrations/WE_6.yaml b/docs/tutorial/T7/calibrations/WE_6.yaml new file mode 100644 index 00000000..02a6554c --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_6.yaml @@ -0,0 +1,21 @@ +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: + (b_rate, 0-4): + starting_y_factor: 1.2 + mig_rate: + lower_bound: 0.1 + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths + + set_initialization: + init_year: 2025 \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/WE_7.yaml b/docs/tutorial/T7/calibrations/WE_7.yaml new file mode 100644 index 00000000..1628e563 --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_7.yaml @@ -0,0 +1,23 @@ +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: + (b_rate, 0-4): + starting_y_factor: 1.2 + mig_rate: + lower_bound: 0.1 + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths + + clear_initialization: True + + set_initialization: + init_year: 2025 \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/WE_8.yaml b/docs/tutorial/T7/calibrations/WE_8.yaml new file mode 100644 index 00000000..849d5750 --- /dev/null +++ b/docs/tutorial/T7/calibrations/WE_8.yaml @@ -0,0 +1,28 @@ + + +calibration: + cal_start: 2005 + cal_end: 2040 + + Population calibration: + repeats: 10 + + Match population sizes: + adjustables: + (b_rate, 0-4): + starting_y_factor: 1.2 + mig_rate: + lower_bound: 0.1 + measurables: alive + + Match deaths: + adjustables: d_rate + measurables: deaths + + clear_initialization: True + + set_initialization: + init_year: 2010 + + save_calibration: + fname: Pop_calibration.xlsx \ No newline at end of file diff --git a/docs/tutorial/T7/calibrations/cal_1.xlsx b/docs/tutorial/T7/calibrations/cal_1.xlsx new file mode 100644 index 00000000..90de0236 Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_1.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_3.xlsx b/docs/tutorial/T7/calibrations/cal_3.xlsx new file mode 100644 index 00000000..fa332aba Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_3.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_4.xlsx b/docs/tutorial/T7/calibrations/cal_4.xlsx new file mode 100644 index 00000000..2e775583 Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_4.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_5.xlsx b/docs/tutorial/T7/calibrations/cal_5.xlsx new file mode 100644 index 00000000..6f8723fa Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_5.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_6.xlsx b/docs/tutorial/T7/calibrations/cal_6.xlsx new file mode 100644 index 00000000..5362c5a5 Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_6.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_7.xlsx b/docs/tutorial/T7/calibrations/cal_7.xlsx new file mode 100644 index 00000000..6a9a4b1e Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_7.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_8.xlsx b/docs/tutorial/T7/calibrations/cal_8.xlsx new file mode 100644 index 00000000..6a9a4b1e Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_8.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_9.xlsx b/docs/tutorial/T7/calibrations/cal_9.xlsx new file mode 100644 index 00000000..b28f538d Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_9.xlsx differ diff --git a/docs/tutorial/T7/calibrations/cal_WE8.xlsx b/docs/tutorial/T7/calibrations/cal_WE8.xlsx new file mode 100644 index 00000000..816c756e Binary files /dev/null and b/docs/tutorial/T7/calibrations/cal_WE8.xlsx differ diff --git a/docs/tutorial/T7/calibrations/run_model_autocals.py b/docs/tutorial/T7/calibrations/run_model_autocals.py new file mode 100644 index 00000000..0dc315d6 --- /dev/null +++ b/docs/tutorial/T7/calibrations/run_model_autocals.py @@ -0,0 +1,80 @@ +# This file produces a calibrations of the whole typhoid model. +# Output files will be saved locally in the outputs folder. +# Desirable results must be manually moved to the calibrations folder so that they can be plotted by the +# Plot_final_typhoid_calibration.ipynb notebook. This is to avoid crowding the final results folder (calibrations) when +# trying out new ways of calibrating the model. + + +import atomica as at +import atomica.yaml_calibration +import os +import time +import shutil +import winsound + +print(os.getcwd()) +# run settings +yaml_dir = "../calibrations" +yaml_fname = "typ_calibration_instructions" +yaml_path = f"{yaml_dir}/{yaml_fname}.yaml" +testing = False # <------------ +savecals = True +dis_model = "typh" +out_path = "outputs" + +# make project ------------------------------------------------ +print(os.getcwd()) +inputs = "../assets" +F = at.ProjectFramework(f"{inputs}/T7_framework.xlsx") +D = at.ProjectData.from_spreadsheet(f"{inputs}/T7_databook.xlsx", framework=F) + +P = at.Project(framework=F, databook=D, do_run=False) +P.settings.update_time_vector(start=2000, end=2050, dt=1 / 52) +cal = P.make_parset() +# starting calibration (if needed) +# cal.load_calibration(out_path + "/PAK_typh_YAML_autocalibrations_2024-04-25_1941_reinitialising/PAK_typh_YAML-autocalibration_2024-04-25_1959.xlsx") +# print('Existing calibration loaded') + +# new folder to put calibrations in +date = time.strftime("%Y-%m-%d_%H%M") +calspath = f"{out_path}/cal_{yaml_fname}" +if testing: + calspath += "_TESTING" + savecals = False +if savecals: + # save everything in folder + try: + os.makedirs(calspath) + except FileExistsError: + calspath += f"_{date}" + os.makedirs(f"{calspath}") + + # save yaml file for reference + new_fname = f"{yaml_fname}_{date}.yaml" + shutil.copy(yaml_path, calspath) + shutil.move(f"{calspath}/{yaml_fname}.yaml", f"{calspath}/{new_fname}.yaml") # renaming file + + # save logfile + at.start_logging(f"{calspath}/logfile_{date}.txt") + + # save current runfile for reference + runfile_name = f"{os.path.basename(__file__)}_{date}" + shutil.copy(__file__, f"{calspath}/{runfile_name}") + +# display tree +cal_tree = at.yaml_calibration.build(yaml_path) +print(cal_tree) + +# run calibration w yaml instructions # <-------------------------- +newcal = P.calibrate(parset=cal, yaml=yaml_path, quiet=False, savedir=calspath, save_intermediate=False, verbose=0) + +# save final calibration +date = time.strftime("%Y-%m-%d_%H%M") +newcal.save_calibration(f"{calspath}/typ_calibration.xlsx") +at.stop_logging() + + +# Notification - calibration finished +winsound.Beep(frequency=2500, duration=200) +winsound.Beep(frequency=2750, duration=200) +winsound.Beep(frequency=3050, duration=200) diff --git a/docs/tutorial/T7/calibrations/run_tutorial_autocals.py b/docs/tutorial/T7/calibrations/run_tutorial_autocals.py new file mode 100644 index 00000000..5510fe43 --- /dev/null +++ b/docs/tutorial/T7/calibrations/run_tutorial_autocals.py @@ -0,0 +1,81 @@ +# This file produces calibrations of each YAML file shown in T7: YAML autocalibration tutorial. +# Output files will be saved locally in the outputs folder. +# Desirable results must be manually moved to the calibrations folder so that they can be plotted by the +# YAML_tutorial_plots_final.ipynb notebook. This is to avoid crowding the final results folder (calibrations) when +# re-running the calibrations or trying out new versions of the model. + +import sys + +in_path = "../inputs" +sys.path.append(in_path) +import atomica as at +import atomica.yaml_calibration +import os +from os.path import isfile +import time +import shutil +import winsound + +yaml_dir = os.getcwd() +fnames = [f for f in os.listdir(yaml_dir) if isfile(f"{yaml_dir}/{f}")] +yaml_files = [f for f in fnames if ".yaml" in f and "T7_YAML" in f] +yaml_files += "WE_8" + +testing = False # <------------ +savecals = True +out_path = "outputs" + +# make project ------------------------------------------------ +inputs = "../assets" +F = at.ProjectFramework(f"{inputs}/T7_framework.xlsx") +D = at.ProjectData.from_spreadsheet(f"{inputs}/T7_databook.xlsx", framework=F) + +P = at.Project(framework=F, databook=D, do_run=False) +P.settings.update_time_vector(start=2000, end=2050, dt=1 / 52) +cal = P.make_parset() +# starting calibration (if needed) +# cal.load_calibration(out_path + "/PAK_typh_YAML_autocalibrations_2024-04-25_1941_reinitialising/PAK_typh_YAML-autocalibration_2024-04-25_1959.xlsx") +# print('Existing calibration loaded') + +for yaml_fname in yaml_files: + yaml_path = f"{yaml_dir}/{yaml_fname}" + + # new folder to put calibrations in + date = time.strftime("%Y-%m-%d_%H%M") + calspath = f"{out_path}/cals_{yaml_fname[8]}" + if savecals: + # save everything in folder + try: + os.makedirs(calspath) + except FileExistsError: + calspath += f"_{date}" + os.makedirs(f"{calspath}") + + # save yaml file for reference + new_fname = f"{yaml_fname}_{date}.yaml" + shutil.copy(yaml_path, calspath) + shutil.move(f"{calspath}/{yaml_fname}", f"{calspath}/{new_fname}") # renaming file + + # save logfile + at.start_logging(f"{calspath}/logfile_{date}.txt") + + # save current runfile for reference + runfile_name = f"{os.path.basename(__file__)}_{date}" + shutil.copy(__file__, f"{calspath}/{runfile_name}") + + # display tree + cal_tree = at.yaml_calibration.build(yaml_path) + print(cal_tree) + + # run calibration w yaml instructions # <-------------------------- + newcal = P.calibrate(parset=cal, yaml=yaml_path, quiet=False, savedir=calspath, save_intermediate=False, verbose=0) + + # save final calibration + newcal.save_calibration(f"{calspath}/cal_{yaml_fname[8]}") # newcal.save_calibration(f'{calspath}/cal_{yaml_fname}') + at.stop_logging() + + +# Notification - calibration finished +winsound.Beep(frequency=2500, duration=200) +winsound.Beep(frequency=2750, duration=200) +winsound.Beep(frequency=3050, duration=200) diff --git a/docs/tutorial/T7/calibrations/typ_calibration.xlsx b/docs/tutorial/T7/calibrations/typ_calibration.xlsx new file mode 100644 index 00000000..db1b00f8 Binary files /dev/null and b/docs/tutorial/T7/calibrations/typ_calibration.xlsx differ diff --git a/docs/tutorial/T7/calibrations/typ_calibration_instructions.yaml b/docs/tutorial/T7/calibrations/typ_calibration_instructions.yaml new file mode 100644 index 00000000..9fc49bdc --- /dev/null +++ b/docs/tutorial/T7/calibrations/typ_calibration_instructions.yaml @@ -0,0 +1,55 @@ +--- +#Timestep for the model +#sim_dt: 0.019 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + repeats: 3 + max_time: 120 + + Population calibration: #-------------------------------------------------- + repeats: 5 + + match population sizes: + adjustables: + (b_rate, 0-4), mig_rate: [ 0.1, 10 ] + (b_rate, 0-4): + initial_value: 1.2 + measurables: alive + + match all-cause deaths: + adjustables: d_rate + measurables: deaths + + pop initialization: + clear_initialization: True + set_initialization: + init_year: 2010 + constant_parset: True + + + Epidemiology calibration: # ----------------------------------------------------- + repeats: 10 + + Typhoid calibration: + + match typhoid incidence: + adjustables: + typ_susceptibility: [0.01, 100 ] + measurables: [ typ_incidence ] + + match typhoid prevalence: + adjustables: [ [ typ_gen_dur, 0.5, 2 ], [ typ_infx, 0.01, 100 ] ] + measurables: [ typ_prev ] + + match typhoid deaths: + adjustables: [ typ_d_rate ] + measurables: [ typ_num_deaths ] + + reinitialization typhoid: + clear_initialization: True + set_initialization: + init_year: 2010 + constant_parset: True diff --git a/docs/tutorial/index.rst b/docs/tutorial/index.rst index 8cf8a11c..ec817ec7 100644 --- a/docs/tutorial/index.rst +++ b/docs/tutorial/index.rst @@ -12,3 +12,4 @@ These tutorials walk through getting starting with building models in Atomica T4-Characteristics T5-Junctions T6-Programs + T7/T7_YAML_autocalibration diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..e70c3c0f --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,56 @@ +[build-system] +requires = ["setuptools>=61.0", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "atomica" +dynamic = ["version"] +description = "Toolbox for compartment-based dynamic systems with costing and optimization" +readme = "README.md" +authors = [ + { name="Atomica Team", email="info@atomica.tools" }, +] +license = { text="MIT" } +keywords = ["dynamic", "compartment", "optimization", "disease"] +classifiers = [ + "Environment :: Console", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Programming Language :: Python", + "Topic :: Software Development :: Libraries :: Python Modules", + "Development Status :: 5 - Production/Stable", + "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", + "Programming Language :: Python :: 3.13" +] +requires-python = ">=3.10" + +dependencies = [ + "matplotlib", + "numpy>=1.10.1", + "scipy>=1.6", + "pandas", + "xlsxwriter", + "openpyxl", + "pyswarm", + "hyperopt", + "sciris", + "tqdm", +] + +[project.urls] +Homepage = "https://atomica.tools" + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +where = ["."] + +[tool.setuptools.dynamic] +version = {attr = "atomica.version.version"} diff --git a/setup.py b/setup.py deleted file mode 100644 index eb71024a..00000000 --- a/setup.py +++ /dev/null @@ -1,59 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -from setuptools import setup, find_packages -import os - -cwd = os.path.abspath(os.path.dirname(__file__)) - -# Read version -with open(os.path.join(cwd, "atomica", "version.py"), "r") as f: - version = [x.split("=")[1].replace('"', "").strip() for x in f if x.startswith("version =")][0] - -# Read README.md for description -with open(os.path.join(cwd, "README.md"), "r") as f: - long_description = f.read() - -CLASSIFIERS = [ - "Environment :: Console", - "Intended Audience :: Science/Research", - "License :: OSI Approved :: MIT License", - "Operating System :: OS Independent", - "Programming Language :: Python", - "Topic :: Software Development :: Libraries :: Python Modules", - "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( - name="atomica", - version=version, - author="Atomica Team", - author_email="info@atomica.tools", - description="Toolbox for compartment-based dynamic systems with costing and optimization", - long_description=long_description, - long_description_content_type="text/markdown", - url="https://atomica.tools", - keywords=["dynamic", "compartment", "optimization", "disease"], - platforms=["OS Independent"], - classifiers=CLASSIFIERS, - packages=find_packages(), - include_package_data=True, - install_requires=[ - "matplotlib", - "numpy>=1.10.1,<2", - "scipy>=1.6", - "pandas", - "xlsxwriter", - "openpyxl", - "pyswarm", - "hyperopt", - "sciris", - "tqdm", - ], -) diff --git a/tests/framework_init_zero.xlsx b/tests/framework_init_zero.xlsx new file mode 100644 index 00000000..b7409836 Binary files /dev/null and b/tests/framework_init_zero.xlsx differ diff --git a/tests/simple.py b/tests/simple.py new file mode 100644 index 00000000..26d6a963 --- /dev/null +++ b/tests/simple.py @@ -0,0 +1,10 @@ +""" +Minimal example of running Atomica. +""" +import atomica as at +import matplotlib.pyplot as plt + +P = at.demo("sir") +d = at.PlotData(P.results[0], project=P) +figs = at.plot_series(d) +plt.show() \ No newline at end of file diff --git a/tests/sir_framework_default_all.xlsx b/tests/sir_framework_default_all.xlsx new file mode 100644 index 00000000..99b7b429 Binary files /dev/null and b/tests/sir_framework_default_all.xlsx differ diff --git a/tests/test_time_vector.py b/tests/test_time_vector.py new file mode 100644 index 00000000..5fc3fae2 --- /dev/null +++ b/tests/test_time_vector.py @@ -0,0 +1,25 @@ +# Test time vector +# The problem: +# In atomica, we want to include the end year: Say we're running a simulation from 2000 to 2050 -> +# we want to be able to get the results for year 2000, 2001, ... and 2050. Hence, we need to run +# more time steps after 2050, in particular up to 2051-1/12 + +# Currently, at.ProjectSettings.tvec creates 611 entries -> if taken the integer,the last year doesn't have enough entries + +import atomica as at +import numpy as np +import pandas as pd + + +def test_time_vector(): + # Test 1 + # Creating a setting object + analysis_years = (2000, 2050) + dt = 1 / 12 + settings = at.ProjectSettings(sim_start=analysis_years[0], sim_end=analysis_years[1] + 1 - dt, sim_dt=dt) #: Atomica project settings + tvec = settings.tvec #: Simulation/result time points + assert set(np.unique(tvec.astype(int), return_counts=True)[1]) == {1 / dt} # or pd.Series(tvec).astype(int).value_counts() + + +if __name__ == "__main__": + test_time_vector() diff --git a/tests/test_tox_databooks.py b/tests/test_tox_databooks.py index 7958424f..03110d42 100644 --- a/tests/test_tox_databooks.py +++ b/tests/test_tox_databooks.py @@ -181,9 +181,27 @@ def test_databook_all(): at.plot_series(d, axis="pops", data=P.data) +def test_databook_default_all(): + F = at.ProjectFramework(testdir / "sir_framework_default_all.xlsx") + D = at.ProjectData.new(F, np.arange(2000, 2005), pops=4, transfers=1) + D.save(tmpdir / "databook_default_all_test.xlsx") # Test saving it back + + D.add_pop("new", "new") + assert list(D.tdve["sus"].ts.keys()) == ["pop_0", "pop_1", "pop_2", "pop_3", "new"] + assert list(D.tdve["ch_prev"].ts.keys()) == ["All"] + + D.tdve["contacts"].ts["pop_0"] = D.tdve["contacts"].ts["All"] + del D.tdve["contacts"].ts["All"] + D.add_pop("new2", "new2") + assert list(D.tdve["ch_prev"].ts.keys()) == ["All"] + assert list(D.tdve["contacts"].ts.keys()) == ["pop_0", "new2"] + D.save(tmpdir / "databook_default_all_test_2.xlsx") # Test saving it back + + if __name__ == "__main__": # test_mixed_years_2() # test_mixed_years_1() # test_databooks() # test_databook_comments() - test_databook_all() + # test_databook_all() + test_databook_default_all() diff --git a/tests/test_tox_frameworks.py b/tests/test_tox_frameworks.py index 303ee535..ccded5d8 100644 --- a/tests/test_tox_frameworks.py +++ b/tests/test_tox_frameworks.py @@ -56,7 +56,7 @@ def test_framework_single_char(): P.settings = P1.settings r2 = P.run_sim() - test_equal = lambda var1, var2: np.testing.assert_array_equal(r1.get_variable(var1)[0].vals, r2.get_variable(var2)[0].vals) + test_equal = lambda var1, var2: np.testing.assert_allclose(r1.get_variable(var1)[0].vals, r2.get_variable(var2)[0].vals) test_equal("sus", "s") test_equal("inf", "i") diff --git a/tests/test_tox_library.py b/tests/test_tox_library.py index 704742dc..110dd246 100644 --- a/tests/test_tox_library.py +++ b/tests/test_tox_library.py @@ -239,6 +239,6 @@ def test_model(model): if __name__ == "__main__": np.seterr(all="raise") - # models = ['tb'] + models = ["combined"] for m in models: test_model(m) diff --git a/tests/test_tox_migration.py b/tests/test_tox_migration.py index 63393827..76c63e4e 100644 --- a/tests/test_tox_migration.py +++ b/tests/test_tox_migration.py @@ -33,6 +33,11 @@ def test_migration(): P.databook.save(tmpdir / "migration_test_databook_save") # Save original databook P.data.save(tmpdir / "migration_test_data_save") # Re-convert data to spreadsheet and save + # Test databook operations + P.data.add_pop("test_pop", "test_pop") # Requires migration + P.data.add_transfer("test_transfer", "test_transfer") + P.data.add_interaction("test_interaction", "test_interaction") + if __name__ == "__main__": test_migration() diff --git a/tests/test_tox_timed_compartments.py b/tests/test_tox_timed_compartments.py index c8db579f..cf306ba1 100644 --- a/tests/test_tox_timed_compartments.py +++ b/tests/test_tox_timed_compartments.py @@ -2,19 +2,8 @@ import sciris as sc import os import pytest -import sys - -# # P = at.Project(framework='dummy_framework.xlsx',databook='dummy_databook.xlsx') -# # # P.run_sim() -# # P = at.demo('tb',do_run=True) -# # P.results[0].plot() -# # P = at.Project.load('asdf2.prj') -# -# # P = at.Project(framework=at.LIBRARY_PATH+'tb_framework.xlsx') -# # P.load_databook(at.LIBRARY_PATH+'tb_databook.xlsx') -# P = at.demo('tb',do_run=True) -# -# +import numpy as np + testdir = at.parent_dir() tmpdir = testdir / "temp" @@ -55,14 +44,15 @@ def test_zero_duration(): P = get_project() ps = P.parsets[0].copy() ps.pars["nat_rec"].ts[0].insert(None, 0) # Sub-timestep size + P.settings.sim_dt = 1/24 # One of the tests had a numerical precision issue on just one platform, hard to reproduce but seeing if this avoids it res2 = P.run_sim(ps) d = at.PlotData(res2, [":inf", "inf:", "inf"]) at.plot_series(d) pop = res2.model.pops[0] - assert pop.get_variable("foi")[0].vals[0] == 24 # Flow into the compartment - assert pop.get_variable("inf")[0].vals[1] == 24 * res2.dt # Compartment contents should now equal the inflow - assert pop.get_variable("inf")[0].vals[2] == 24 * res2.dt # Same again, contents equals the inflow because it was flushed entirely + assert np.isclose(pop.get_variable("foi")[0].vals[0], 24) # Flow into the compartment + assert np.isclose(pop.get_variable("inf")[0].vals[1], 24 * res2.dt) # Compartment contents should now equal the inflow + assert np.isclose(pop.get_variable("inf")[0].vals[2], 24 * res2.dt) # Same again, contents equals the inflow because it was flushed entirely # Check formally that total inflows equal total outflows assert pop.get_variable("inf")[0].vals[0] == sum(link.vals[0] for link in pop.get_variable("inf")[0].outlinks) diff --git a/tests/test_tox_yaml.py b/tests/test_tox_yaml.py new file mode 100644 index 00000000..53c90bf9 --- /dev/null +++ b/tests/test_tox_yaml.py @@ -0,0 +1,54 @@ +# import unittest +import atomica as at +import os +import pytest +import numpy as np +import sciris as sc + +# Basic minimum working product test to verify all yaml features run +testdir = at.parent_dir() +tmpdir = testdir / "temp" +yaml_dir = testdir / "yaml_tests" + +# List yaml files in directory +yaml_files = list(yaml_dir.glob("*.yaml")) + + +@pytest.fixture() +def project(): + P = at.demo("udt_dyn", do_run=False) + P.settings.update_time_vector(start=2000, end=2050, dt=1 / 52) + return P + + +def test_save_intermediate(project): + project.calibrate(yaml=yaml_files[0], parset=project.make_parset(), quiet=False, savedir=tmpdir / "yaml_save_test", save_intermediate=True, verbose=0, max_time=0.1) + + +@pytest.mark.parametrize("yaml_fname", yaml_files, ids=[x.stem for x in yaml_files]) +def test_yaml_calibration(project, yaml_fname): + + # TODO - this demo project has only one population. It would be good + # to have a test with multiple populations and where there is a space in the population name + + print(f"\n\nTESTING yaml config {yaml_fname}") + + # run calibration with yaml instructions + # config_file_path = yaml_dir/yaml_fname + # print(config_file_path) + project.calibrate(yaml=yaml_fname, parset=project.make_parset(), quiet=False, savedir=tmpdir, save_intermediate=False, verbose=0, max_time=0.1) + + # or can import in all the diff ways - make function + # also run in several diff ways? Or have a separate test for that? + + # Test saving the calibration + # save() + + print("Test complete") + + +if __name__ == "__main__": + + np.seterr(all="raise") + for f in yaml_files: + test_yaml_calibration(f) diff --git a/tests/test_tox_zero_init.py b/tests/test_tox_zero_init.py new file mode 100644 index 00000000..c9ca9956 --- /dev/null +++ b/tests/test_tox_zero_init.py @@ -0,0 +1,23 @@ +# Test frameworks where characteristics with zero values get extra initialization conditions for constituent compartments + +import atomica as at +import numpy as np +from numpy.testing import assert_allclose +import pytest + +testdir = at.parent_dir() +tmpdir = testdir / "temp" + + +def test_zero_initialization(): + F = at.ProjectFramework(testdir / "framework_init_zero.xlsx") + D = at.ProjectData.new(F, [2020], 1, 0) + P = at.Project(framework=F, databook=D, do_run=False) + P.settings.update_time_vector(start=2020 ,end=2021, dt=1) + res = P.run_sim() + return res + + + +if __name__ == "__main__": + test_zero_initialization() 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() diff --git a/tests/timed_initialization_test_framework.xlsx b/tests/timed_initialization_test_framework.xlsx new file mode 100644 index 00000000..796297da Binary files /dev/null and b/tests/timed_initialization_test_framework.xlsx differ diff --git a/tests/yaml_tests/config_bounds_test.yaml b/tests/yaml_tests/config_bounds_test.yaml new file mode 100644 index 00000000..13c4a025 --- /dev/null +++ b/tests/yaml_tests/config_bounds_test.yaml @@ -0,0 +1,13 @@ + +calibration: + pop cal: + + match population sizes: + adjustables: + b_rate, inf_death: [ 0.1, 10 ] + measurables: all_people + + turn off epi y-factors: + adjustables: + typ_active_inf, typ_car, par_active_inf, par_car: [0, 0] #set y-factor w/o going through cal algorithm + measurables: deaths \ No newline at end of file diff --git a/tests/yaml_tests/config_cal_formats.yaml b/tests/yaml_tests/config_cal_formats.yaml new file mode 100644 index 00000000..85eb63f0 --- /dev/null +++ b/tests/yaml_tests/config_cal_formats.yaml @@ -0,0 +1,47 @@ +--- +#Timestep for the model +sim_dt: 0.019 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + + Population calibration: + + #dict format + match population sizes 1: + adjustables: + b_rate: + lower_bound: 0.1 + upper_bound: 10 + inf_death: + starting_y_factor: 1.2 + measurables: + all_people: + metric: fractional + + #combined format + match population sizes 2: + version 1: + adjustables: + b_rate: [ 0.1,10 ] + inf_death: [ 0.1,10 ] + measurables: + all_people: [1.0] + + version 2: + adjustables: + b_rate, inf_death: [ 0.1,10 ] + measurables: + all_people: [1.0] + + #list format + match population sizes 3: + adjustables: [[b_rate, 0.1, 10], [inf_death, 0.1, 10]] + measurables: [all_people] + + #string format + match population sizes 4: + adjustables: [b_rate, inf_death] + measurables: all_people diff --git a/tests/yaml_tests/config_caltime_test.yaml b/tests/yaml_tests/config_caltime_test.yaml new file mode 100644 index 00000000..deb40fb0 --- /dev/null +++ b/tests/yaml_tests/config_caltime_test.yaml @@ -0,0 +1,39 @@ +--- +#Timestep for the model +sim_dt: 0.019 +cal_start: 2000 +cal_end: 2040 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + cal_start: 2000 + cal_end: 2040 + + Population calibration: + + round 1: + cal_start: 2000 + cal_end: 2040 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + cal_start: 2000 #is this meant to work? + cal_end: 2040 + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 3: + match population sizes 3: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: + all_people: + cal_start: 2000 + cal_end: 2040 diff --git a/tests/yaml_tests/config_hierarchical_syntax_test.yaml b/tests/yaml_tests/config_hierarchical_syntax_test.yaml new file mode 100644 index 00000000..9dcc1ba3 --- /dev/null +++ b/tests/yaml_tests/config_hierarchical_syntax_test.yaml @@ -0,0 +1,23 @@ +calibration: + pop cal: + repeats: 1 + + match population sizes: + adjustables: + #general, then pop-specific + b_rate, mig_rate: [ 0.1, 10 ] + (b_rate, 0-4): + initial_value: 2 + upper_bound: 10 + lower_bound: 0.1 + + typ_active_inf, typ_car, par_active_inf, par_car: [0, 0] + measurables: alive + + + other formats, hierarchical tests: + adjustables: #[ [ b_rate, 0.1, 10 ], [ mig_rate ] ] + b_rate, mig_rate: [ 0.1, 10 ] + (b_rate, 0-4), mig_rate: [ 0.5, 20 ] + (b_rate, 0-4): + initial_value: 1.2 \ No newline at end of file diff --git a/tests/yaml_tests/config_maxtime_test.yaml b/tests/yaml_tests/config_maxtime_test.yaml new file mode 100644 index 00000000..95d8f9a8 --- /dev/null +++ b/tests/yaml_tests/config_maxtime_test.yaml @@ -0,0 +1,36 @@ +--- +#Timestep for the model +sim_dt: 0.019 +cal_start: 2000 +cal_end: 2040 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + max_time: 1 + + Population calibration: + + round 1: + max_time: 1 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + max_time: 1 + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + +# round 3: +# match population sizes 3: +# adjustables: +# b_rate,inf_death: [ 0.1,10 ] +# b_rate: +# cal_start: 2000 #is this meant to work? +# cal_end: 2040 +# measurables: [ all_people ] diff --git a/tests/yaml_tests/config_nonexistent_setting_test.yaml b/tests/yaml_tests/config_nonexistent_setting_test.yaml new file mode 100644 index 00000000..cd7e6b44 --- /dev/null +++ b/tests/yaml_tests/config_nonexistent_setting_test.yaml @@ -0,0 +1,31 @@ +--- +#Timestep for the model +sim_dt: 0.019 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + strawberry: 2 + Population calibration: + + round 1: + raddish: 3 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + kiwi: 4 + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + +# round 3: +# match population sizes 3: +# adjustables: +# stepsize: 0.5 +# b_rate,inf_death: [ 0.1,10 ] +# measurables: [ all_people ] diff --git a/tests/yaml_tests/config_repeats_test.yaml b/tests/yaml_tests/config_repeats_test.yaml new file mode 100644 index 00000000..dcd60944 --- /dev/null +++ b/tests/yaml_tests/config_repeats_test.yaml @@ -0,0 +1,36 @@ +--- +#Timestep for the model +sim_dt: 0.019 +cal_start: 2000 +cal_end: 2040 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + repeats: 1 + + Population calibration: + repeats: 1 + round 1: + repeats: 1 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + repeats: 1 #should error? + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + +# round 3: +# match population sizes 3: +# adjustables: +# b_rate,inf_death: [ 0.1,10 ] +# b_rate: +# cal_start: 2000 #is this meant to work? +# cal_end: 2040 +# measurables: [ all_people ] diff --git a/tests/yaml_tests/config_simtime_test.yaml b/tests/yaml_tests/config_simtime_test.yaml new file mode 100644 index 00000000..2800f5e4 --- /dev/null +++ b/tests/yaml_tests/config_simtime_test.yaml @@ -0,0 +1,30 @@ +--- +#Timestep for the model +sim_dt: 0.019 +sim_start: 2000 +sim_end: 2040 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + sim_start: 2000 + sim_end: 2040 + + Population calibration: + + round 1: + sim_start: 2000 + sim_end: 2040 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + sim_start: 2000 #is this meant to work? + sim_end: 2040 + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] diff --git a/tests/yaml_tests/config_stepsize_test.yaml b/tests/yaml_tests/config_stepsize_test.yaml new file mode 100644 index 00000000..d94db8f3 --- /dev/null +++ b/tests/yaml_tests/config_stepsize_test.yaml @@ -0,0 +1,32 @@ +--- +#Timestep for the model +sim_dt: 0.019 + +#Ordered instructions for how to calibrate the model + #adj order: (adjustable var_label, populations to adjust for (default all), minimum value (default 0.1), maximum value (default 10), initial y_factor value (optional, default is parset value)) + #meas order: #(measurable var_label, populations to measure for (default all), weight (default 1.0), metric (default "fractional"), start_year (optional), end_year (optional)) +calibration: + stepsize: 1 + + Population calibration: + + round 1: + stepsize: 0.5 + match population sizes 1: + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + + round 2: + match population sizes 2: + stepsize: 0.5 + adjustables: + b_rate,inf_death: [ 0.1,10 ] + measurables: [all_people] + +# round 3: +# match population sizes 3: +# adjustables: +# stepsize: 0.5 +# b_rate,inf_death: [ 0.1,10 ] +# measurables: [ all_people ] diff --git a/tox.ini b/tox.ini index e81a142c..889520bb 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist = py311 +envlist = py312 [testenv] description = Run basic usage tests @@ -12,7 +12,7 @@ commands = pytest --nbsmoke-run {toxinidir} {posargs} [pytest] -addopts = -ra -v -n 2 --junitxml=junit/test-results.xml --color=no --cov=atomica --cov-report=xml --cov-report=html +addopts = -ra -v -n 2 --junitxml=junit/test-results.xml --color=no --cov=atomica --cov-report=xml --cov-report=html --log-level=DEBUG junit_family=xunit1 nbsmoke_cell_timeout = 600 console_output_style = count @@ -29,20 +29,3 @@ filterwarnings = ignore::PendingDeprecationWarning ignore::UserWarning:.*openpyxl -[flake8] -# Use these commands to call autopep8 to tidy up whitespace/trivial formatting -# -# autopep8 --in-place -j 0 --recursive atomica -# autopep8 --in-place -j 0 --recursive tests -# -# Alternatively, use Black via -# -# pip install black -# black -l 999 . # Use `-l 999` to more or less ignore line length - -exclude = .git,__pycache__,docs,build,dist,junit,.tox -ignore = E501,F401,E401,E402,C901,F403 - -# flake8 does not check with E203 properly with list slices -extend-ignore = E203 -