From c8640fffdfc552ea5c7ade0a20e4200ec53ca39d Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 5 Nov 2024 16:00:26 +0100 Subject: [PATCH 01/35] ESSOptimizer: Fix priority for local search startpoints (#1503) Fixes a bug in ESSOptimizer that lead to incorrect priorities for the selection of candidate solutions for local searches. A numpy array was incorrectly initialized from a generator, which led to the `balance` hyperparameter being practically ignored, and candidate solutions being ranked only by objective function value. --- pypesto/optimize/ess/ess.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index ca8ffa2f6..1f31df0f5 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -487,14 +487,19 @@ def _do_local_search( quality_order = np.argsort(fx_best_children) # compute minimal distance between the best children and all local # optima found so far - min_distances = np.array( - np.min( - np.linalg.norm( - y_i - optimizer_result.x[optimizer_result.free_indices] + min_distances = np.fromiter( + ( + min( + np.linalg.norm( + y_i + - optimizer_result.x[optimizer_result.free_indices] + ) + for optimizer_result in self.local_solutions ) - for optimizer_result in self.local_solutions - ) - for y_i in x_best_children + for y_i in x_best_children + ), + dtype=np.float64, + count=len(x_best_children), ) # sort by furthest distance to existing local optima diversity_order = np.argsort(min_distances)[::-1] From d5eedc205b7c7f6af941cbc0138c249d419a8a03 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Nov 2024 17:46:36 +0100 Subject: [PATCH 02/35] AmiciObjective/PEtab import: Fix plist (#1493) For further background, see #1448. We don't have to compute the full gradient in every model simulation. Some entries might not be required because of fixed parameter or some condition-specific objective parameter mapping. The former was supported (however, not tested), the latter was not supported. Now both are tested an supported. There was no good way to communicate the fixed parameters to AmiciCalculator where ExpData.plist gets set. Passing this information is currently only possible through the parameter mapping, based on which the required sensitivities are determined. Therefore, in addition to the general parameter mapping, there is now a parameter mapping that accounts for the fixed parameters. Not sure what could be a more elegant way to handle that. --- pypesto/objective/amici/amici.py | 41 +++++++++++++++- test/petab/test_petab_import.py | 84 +++++++++++++++++++++++--------- 2 files changed, 102 insertions(+), 23 deletions(-) diff --git a/pypesto/objective/amici/amici.py b/pypesto/objective/amici/amici.py index 3af3f7023..8edcf13fa 100644 --- a/pypesto/objective/amici/amici.py +++ b/pypesto/objective/amici/amici.py @@ -178,8 +178,12 @@ def __init__( parameter_mapping = create_identity_parameter_mapping( amici_model, len(edatas) ) + # parameter mapping where IDs of the currently fixed parameters + # have been replaced by their respective values + # (relevant for setting ``plist`` in ExpData later on) self.parameter_mapping = parameter_mapping - + # parameter mapping independent of fixed parameters + self._parameter_mapping_full = copy.deepcopy(parameter_mapping) # If supported, enable `guess_steadystate` by default. If not # supported, disable by default. If requested but unsupported, raise. if ( @@ -654,3 +658,38 @@ def check_gradients_match_finite_differences( return super().check_gradients_match_finite_differences( *args, x=x, x_free=x_free, **kwargs ) + + def update_from_problem( + self, + dim_full: int, + x_free_indices: Sequence[int], + x_fixed_indices: Sequence[int], + x_fixed_vals: Sequence[float], + ): + """Handle fixed parameters.""" + super().update_from_problem( + dim_full=dim_full, + x_free_indices=x_free_indices, + x_fixed_indices=x_fixed_indices, + x_fixed_vals=x_fixed_vals, + ) + + # To make amici aware of fixed parameters, and thus, avoid computing + # unnecessary sensitivities, we need to update the parameter mapping + # and replace the IDs of all fixed parameters by their respective + # values. + self.parameter_mapping = copy.deepcopy(self._parameter_mapping_full) + if not len(x_fixed_indices): + return + + id_to_val = { + self.x_ids[x_idx]: x_val + for x_idx, x_val in zip(x_fixed_indices, x_fixed_vals) + } + for condition_mapping in self.parameter_mapping: + for ( + model_par, + mapped_to_par, + ) in condition_mapping.map_sim_var.items(): + if (val := id_to_val.get(mapped_to_par)) is not None: + condition_mapping.map_sim_var[model_par] = val diff --git a/test/petab/test_petab_import.py b/test/petab/test_petab_import.py index aa4eb0067..1ce935880 100644 --- a/test/petab/test_petab_import.py +++ b/test/petab/test_petab_import.py @@ -5,6 +5,7 @@ import logging import os import unittest +from itertools import chain import amici import benchmark_models_petab as models @@ -18,8 +19,6 @@ import pypesto.petab from pypesto.petab import PetabImporter -from .test_sbml_conversion import ATOL, RTOL - # In CI, bionetgen is installed here BNGPATH = os.path.abspath( os.path.join(os.path.dirname(__file__), "..", "..", "BioNetGen-2.8.5") @@ -119,7 +118,7 @@ def test_check_gradients(self): # Check gradients of simple model (should always be a true positive) model_name = "Boehm_JProteomeRes2014" importer = pypesto.petab.PetabImporter.from_yaml( - os.path.join(models.MODELS_DIR, model_name, model_name + ".yaml") + models.get_problem_yaml_path(model_name) ) objective = importer.create_problem().objective @@ -137,35 +136,76 @@ def test_check_gradients(self): def test_plist_mapping(): - """Test that the AMICI objective created via PEtab correctly maps - gradient entries when some parameters are not estimated (realized via - edata.plist).""" - model_name = "Boehm_JProteomeRes2014" - petab_problem = pypesto.petab.PetabImporter.from_yaml( + """Test that the AMICI objective created via PEtab correctly uses + ExpData.plist. + + That means, ensure that + 1) it only computes gradient entries that are really required, and + 2) correctly maps model-gradient entries to the objective gradient + 3) with or without pypesto-fixed parameters. + + In Bruno_JExpBot2016, different parameters are estimated in different + conditions. + """ + model_name = "Bruno_JExpBot2016" + petab_importer = pypesto.petab.PetabImporter.from_yaml( os.path.join(models.MODELS_DIR, model_name, model_name + ".yaml") ) - - # define test parameter - par = np.asarray(petab_problem.petab_problem.x_nominal_scaled) - - problem = petab_problem.create_problem() + objective_creator = petab_importer.create_objective_creator() + problem = petab_importer.create_problem( + objective_creator.create_objective() + ) objective = problem.objective - # check that x_names are correctly subsetted - assert objective.x_names == [ - problem.x_names[ix] for ix in problem.x_free_indices - ] objective.amici_solver.setSensitivityMethod( - amici.SensitivityMethod_forward + amici.SensitivityMethod.forward ) - objective.amici_solver.setAbsoluteTolerance(1e-10) - objective.amici_solver.setRelativeTolerance(1e-12) + objective.amici_solver.setAbsoluteTolerance(1e-16) + objective.amici_solver.setRelativeTolerance(1e-15) + + # slightly perturb the parameters to avoid vanishing gradients + par = np.asarray(petab_importer.petab_problem.x_nominal_free_scaled) * 1.01 + + # call once to make sure ExpDatas are populated + fx1, grad1 = objective(par, sensi_orders=(0, 1)) + + plists1 = {edata.plist for edata in objective.edatas} + assert len(plists1) == 6, plists1 df = objective.check_grad_multi_eps( - par[problem.x_free_indices], multi_eps=[1e-3, 1e-4, 1e-5] + par[problem.x_free_indices], multi_eps=[1e-5, 1e-6, 1e-7, 1e-8] ) print("relative errors gradient: ", df.rel_err.values) print("absolute errors gradient: ", df.abs_err.values) - assert np.all((df.rel_err.values < RTOL) | (df.abs_err.values < ATOL)) + assert np.all((df.rel_err.values < 1e-8) | (df.abs_err.values < 1e-10)) + + # do the same after fixing some parameters + # we fix them to the previous values, so we can compare the results + fixed_ids = ["init_b10_1", "init_bcry_1"] + # the corresponding amici parameters (they are only ever mapped to the + # respective parameters in fixed_ids, and thus, should not occur in any + # `plist` later on) + fixed_model_par_ids = ["init_b10", "init_bcry"] + fixed_model_par_idxs = [ + objective.amici_model.getParameterIds().index(id) + for id in fixed_model_par_ids + ] + fixed_idxs = [problem.x_names.index(id) for id in fixed_ids] + problem.fix_parameters(fixed_idxs, par[fixed_idxs]) + assert objective is problem.objective + assert problem.x_fixed_indices == fixed_idxs + fx2, grad2 = objective(par[problem.x_free_indices], sensi_orders=(0, 1)) + assert np.isclose(fx1, fx2, rtol=1e-10, atol=1e-14) + assert np.allclose( + grad1[problem.x_free_indices], grad2, rtol=1e-10, atol=1e-14 + ) + plists2 = {edata.plist for edata in objective.edatas} + # the fixed parameters should have been in plist1, but not in plist2 + assert ( + set(fixed_model_par_idxs) - set(chain.from_iterable(plists1)) == set() + ) + assert ( + set(fixed_model_par_idxs) & set(chain.from_iterable(plists2)) == set() + ) def test_max_sensi_order(): From 9a69f3ce577dc584e8707f7af85f92169f252e72 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Tue, 12 Nov 2024 00:13:40 +0100 Subject: [PATCH 03/35] Select: support user-provided calibration results (#1338) * let users supply calibration results * handle calibrated model via petab select * custom petab select branch * functionality moved to petab-select * update method for next petab_select version * update for petab select * update for next petab-select version --------- Co-authored-by: Doresic <85789271+Doresic@users.noreply.github.com> --- pypesto/select/method.py | 80 ++++++++++++++++++++++++-------------- pypesto/select/problem.py | 40 +++++-------------- setup.cfg | 6 ++- test/select/test_select.py | 2 + 4 files changed, 67 insertions(+), 61 deletions(-) diff --git a/pypesto/select/method.py b/pypesto/select/method.py index 36a726be6..d1ff12443 100644 --- a/pypesto/select/method.py +++ b/pypesto/select/method.py @@ -8,6 +8,10 @@ import numpy as np import petab_select from petab_select import ( + CANDIDATE_SPACE, + MODELS, + PREDECESSOR_MODEL, + UNCALIBRATED_MODELS, VIRTUAL_INITIAL_MODEL, CandidateSpace, Criterion, @@ -213,6 +217,11 @@ class MethodCaller: Specify the predecessor (initial) model for the model selection algorithm. If ``None``, then the algorithm will generate an initial predecessor model if required. + user_calibrated_models: + Supply calibration results for models yourself, as a list of models. + If a model with the same hash is encountered in the current model + selection run, and the user-supplied calibrated model has the + `criterion` value set, the model will not be calibrated again. select_first_improvement: If ``True``, model selection will terminate as soon as a better model is found. If `False`, all candidate models will be tested. @@ -245,6 +254,7 @@ def __init__( # TODO deprecated model_to_pypesto_problem_method: Callable[[Any], Problem] = None, model_problem_options: dict = None, + user_calibrated_models: list[Model] = None, ): """Arguments are used in every `__call__`, unless overridden.""" self.petab_select_problem = petab_select_problem @@ -256,6 +266,12 @@ def __init__( self.select_first_improvement = select_first_improvement self.startpoint_latest_mle = startpoint_latest_mle + self.user_calibrated_models = {} + if user_calibrated_models is not None: + self.user_calibrated_models = { + model.get_hash(): model for model in user_calibrated_models + } + self.logger = MethodLogger() # TODO deprecated @@ -335,10 +351,7 @@ def __init__( # May have changed from `None` to `petab_select.VIRTUAL_INITIAL_MODEL` self.predecessor_model = self.candidate_space.get_predecessor_model() - def __call__( - self, - newly_calibrated_models: Optional[dict[str, Model]] = None, - ) -> tuple[list[Model], dict[str, Model]]: + def __call__(self) -> tuple[list[Model], dict[str, Model]]: """Run a single iteration of the model selection method. A single iteration here refers to calibration of all candidate models. @@ -347,14 +360,6 @@ def __call__( of all models that have both: the same 3 estimated parameters; and 1 additional estimated parameter. - The input `newly_calibrated_models` is from the previous iteration. The - output `newly_calibrated_models` is from the current iteration. - - Parameters - ---------- - newly_calibrated_models: - The newly calibrated models from the previous iteration. - Returns ------- A 2-tuple, with the following values: @@ -366,39 +371,56 @@ def __call__( # All calibrated models in this iteration (see second return value). self.logger.new_selection() - candidate_space = petab_select.ui.candidates( + iteration = petab_select.ui.start_iteration( problem=self.petab_select_problem, candidate_space=self.candidate_space, limit=self.limit, - calibrated_models=self.calibrated_models, - newly_calibrated_models=newly_calibrated_models, - excluded_model_hashes=self.calibrated_models.keys(), criterion=self.criterion, + user_calibrated_models=self.user_calibrated_models, ) - predecessor_model = self.candidate_space.predecessor_model - if not candidate_space.models: + if not iteration[UNCALIBRATED_MODELS]: raise StopIteration("No valid models found.") # TODO parallelize calibration (maybe not sensible if # `self.select_first_improvement`) - newly_calibrated_models = {} - for candidate_model in candidate_space.models: - # autoruns calibration - self.new_model_problem(model=candidate_model) - newly_calibrated_models[ - candidate_model.get_hash() - ] = candidate_model + calibrated_models = {} + for model in iteration[UNCALIBRATED_MODELS]: + if ( + model.get_criterion( + criterion=self.criterion, + compute=True, + raise_on_failure=False, + ) + is not None + ): + self.logger.log( + message=( + "Unexpected calibration result already available for " + f"model: `{model.get_hash()}`. Skipping " + "calibration." + ), + level="warning", + ) + else: + self.new_model_problem(model=model) + + calibrated_models[model.get_hash()] = model method_signal = self.handle_calibrated_model( - model=candidate_model, - predecessor_model=predecessor_model, + model=model, + predecessor_model=iteration[PREDECESSOR_MODEL], ) if method_signal.proceed == MethodSignalProceed.STOP: break - self.calibrated_models.update(newly_calibrated_models) + iteration_results = petab_select.ui.end_iteration( + candidate_space=iteration[CANDIDATE_SPACE], + calibrated_models=calibrated_models, + ) + + self.calibrated_models.update(iteration_results[MODELS]) - return predecessor_model, newly_calibrated_models + return iteration[PREDECESSOR_MODEL], iteration_results[MODELS] def handle_calibrated_model( self, diff --git a/pypesto/select/problem.py b/pypesto/select/problem.py index 9692890ef..e0bd08cf3 100644 --- a/pypesto/select/problem.py +++ b/pypesto/select/problem.py @@ -156,16 +156,13 @@ def select( self.handle_select_kwargs(kwargs) # TODO handle bidirectional method_caller = self.create_method_caller(**kwargs) - previous_best_model, newly_calibrated_models = method_caller( - # TODO add predecessor model to state - newly_calibrated_models=self.newly_calibrated_models, - ) + previous_best_model, newly_calibrated_models = method_caller() self.update_with_newly_calibrated_models( newly_calibrated_models=newly_calibrated_models, ) - best_model = petab_select.ui.best( + best_model = petab_select.ui.get_best( problem=self.petab_select_problem, models=self.newly_calibrated_models.values(), criterion=method_caller.criterion, @@ -198,9 +195,7 @@ def select_to_completion( while True: try: - previous_best_model, newly_calibrated_models = method_caller( - newly_calibrated_models=self.newly_calibrated_models, - ) + previous_best_model, newly_calibrated_models = method_caller() self.update_with_newly_calibrated_models( newly_calibrated_models=newly_calibrated_models, ) @@ -247,33 +242,18 @@ def multistart_select( """ self.handle_select_kwargs(kwargs) model_lists = [] - newly_calibrated_models_list = [ - self.newly_calibrated_models for _ in predecessor_models - ] - method_caller = self.create_method_caller(**kwargs) - for start_index, predecessor_model in enumerate(predecessor_models): - method_caller.candidate_space.previous_predecessor_model = ( - predecessor_model - ) - ( - best_model, - newly_calibrated_models_list[start_index], - ) = method_caller( - newly_calibrated_models=newly_calibrated_models_list[ - start_index - ], - ) - self.calibrated_models.update( - newly_calibrated_models_list[start_index] + for predecessor_model in predecessor_models: + method_caller = self.create_method_caller( + **(kwargs | {"predecessor_model": predecessor_model}) ) + (best_model, models) = method_caller() + self.calibrated_models |= models - model_lists.append( - newly_calibrated_models_list[start_index].values() - ) + model_lists.append(list(models.values())) method_caller.candidate_space.reset() - best_model = petab_select.ui.best( + best_model = petab_select.ui.get_best( problem=method_caller.petab_select_problem, models=[model for models in model_lists for model in models], criterion=method_caller.criterion, diff --git a/setup.cfg b/setup.cfg index 4fc1a4e1b..bc28ae610 100644 --- a/setup.cfg +++ b/setup.cfg @@ -163,10 +163,12 @@ example = ipywidgets >= 8.1.5 benchmark_models_petab @ git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master#subdirectory=src/python select = - # Remove when vis is moved to PEtab Select version + # Remove when vis is moved to PEtab Select networkx >= 2.5.1 # End remove - petab-select >= 0.1.12 + #petab-select >= 0.1.12 + # FIXME before merge + petab-select @ git+https://github.com/PEtab-dev/petab_select.git@develop test = pytest >= 5.4.3 pytest-cov >= 2.10.0 diff --git a/test/select/test_select.py b/test/select/test_select.py index b471365d6..5ba2d5a3e 100644 --- a/test/select/test_select.py +++ b/test/select/test_select.py @@ -279,6 +279,8 @@ def test_problem_multistart_select(pypesto_select_problem, initial_models): "M1_3": -4.705, # 'M1_7': -4.056, # skipped -- reproducibility requires many starts } + # As M1_7 criterion comparison is skipped, at least ensure it is present + assert {m.model_subspace_id for m in best_models} == {"M1_3", "M1_7"} test_best_models_criterion_values = { model.model_subspace_id: model.get_criterion(Criterion.AIC) for model in best_models From 80a7b8d09942be77800f27f8127b190202edd289 Mon Sep 17 00:00:00 2001 From: Doresic <85789271+Doresic@users.noreply.github.com> Date: Tue, 12 Nov 2024 09:47:42 +0100 Subject: [PATCH 04/35] Hierarchical: fix no error if inner observable parameter in noise formula & viceversa (#1504) * Search for observable inner pars in noise formulas & vice-versa * Fix if statements * Fix different censoring_types are allowed --- pypesto/hierarchical/petab.py | 135 +++++++++++++++++++++------------- 1 file changed, 82 insertions(+), 53 deletions(-) diff --git a/pypesto/hierarchical/petab.py b/pypesto/hierarchical/petab.py index c1b85b3be..dde1d0dad 100644 --- a/pypesto/hierarchical/petab.py +++ b/pypesto/hierarchical/petab.py @@ -76,63 +76,63 @@ def validate_hierarchical_petab_problem(petab_problem: petab.Problem) -> None: petab_problem: The PEtab problem. """ - if PARAMETER_TYPE not in petab_problem.parameter_df: - # not a hierarchical optimization problem - return - - # ensure we only have linear parameter scale - inner_parameter_table = petab_problem.parameter_df[ - petab_problem.parameter_df[PARAMETER_TYPE].isin( - [ - InnerParameterType.OFFSET, - InnerParameterType.SIGMA, - InnerParameterType.SCALING, - ] - ) - ] - if ( - petab.PARAMETER_SCALE in inner_parameter_table - and not ( - inner_parameter_table[petab.PARAMETER_SCALE].isna() - | (inner_parameter_table[petab.PARAMETER_SCALE] == petab.LIN) - | ( - inner_parameter_table[PARAMETER_TYPE] - != InnerParameterType.SIGMA + if PARAMETER_TYPE in petab_problem.parameter_df: + # ensure we only have linear parameter scale + inner_parameter_table = petab_problem.parameter_df[ + petab_problem.parameter_df[PARAMETER_TYPE].isin( + [ + InnerParameterType.OFFSET, + InnerParameterType.SIGMA, + InnerParameterType.SCALING, + ] ) - ).all() - ): - sub_df = inner_parameter_table.loc[ - :, [PARAMETER_TYPE, petab.PARAMETER_SCALE] ] - raise NotImplementedError( - "LOG and LOG10 parameter scale of inner parameters is not supported " - "for sigma parameters. Inner parameter table:\n" - f"{sub_df}" - ) - elif ( - petab.PARAMETER_SCALE in inner_parameter_table - and not ( - inner_parameter_table[petab.PARAMETER_SCALE].isna() - | (inner_parameter_table[petab.PARAMETER_SCALE] == petab.LIN) - ).all() - ): - sub_df = inner_parameter_table.loc[ - :, [PARAMETER_TYPE, petab.PARAMETER_SCALE] - ] - warnings.warn( - f"LOG and LOG10 parameter scale of inner parameters is used only " - f"for their visualization, and does not affect their optimization. " - f"Inner parameter table:\n{sub_df}", - stacklevel=1, - ) + if ( + petab.PARAMETER_SCALE in inner_parameter_table + and not ( + inner_parameter_table[petab.PARAMETER_SCALE].isna() + | (inner_parameter_table[petab.PARAMETER_SCALE] == petab.LIN) + | ( + inner_parameter_table[PARAMETER_TYPE] + != InnerParameterType.SIGMA + ) + ).all() + ): + sub_df = inner_parameter_table.loc[ + :, [PARAMETER_TYPE, petab.PARAMETER_SCALE] + ] + raise NotImplementedError( + "LOG and LOG10 parameter scale of inner parameters is not supported " + "for sigma parameters. Inner parameter table:\n" + f"{sub_df}" + ) + elif ( + petab.PARAMETER_SCALE in inner_parameter_table + and not ( + inner_parameter_table[petab.PARAMETER_SCALE].isna() + | (inner_parameter_table[petab.PARAMETER_SCALE] == petab.LIN) + ).all() + ): + sub_df = inner_parameter_table.loc[ + :, [PARAMETER_TYPE, petab.PARAMETER_SCALE] + ] + warnings.warn( + f"LOG and LOG10 parameter scale of inner parameters is used only " + f"for their visualization, and does not affect their optimization. " + f"Inner parameter table:\n{sub_df}", + stacklevel=1, + ) - inner_parameter_df = validate_measurement_formulae( - petab_problem=petab_problem - ) + inner_parameter_df = validate_measurement_formulae( + petab_problem=petab_problem + ) - validate_inner_parameter_pairings(inner_parameter_df=inner_parameter_df) + validate_inner_parameter_pairings( + inner_parameter_df=inner_parameter_df + ) - validate_observable_data_types(petab_problem=petab_problem) + if MEASUREMENT_TYPE in petab_problem.measurement_df: + validate_observable_data_types(petab_problem=petab_problem) def validate_inner_parameter_pairings( @@ -477,6 +477,18 @@ def _get_symbolic_formula_from_measurement( observable_id=observable_id, override_type=formula_type, ) + # Search for placeholders in the formula that are not allowed to be + # overridden by inner parameters -- noise inner parameters should not + # be in observable formulas, and vice versa. + disallowed_formula_type = ( + "noise" if formula_type == "observable" else "observable" + ) + disallowed_formula_placeholders = get_formula_placeholders( + formula_string=formula_string, + observable_id=observable_id, + override_type=disallowed_formula_type, + ) + if formula_placeholders: overrides = measurement[formula_type + "Parameters"] overrides = ( @@ -487,6 +499,20 @@ def _get_symbolic_formula_from_measurement( subs = dict(zip(formula_placeholders, overrides)) symbolic_formula = symbolic_formula.subs(subs) + if disallowed_formula_placeholders: + disallowed_overrides = measurement[ + disallowed_formula_type + "Parameters" + ] + disallowed_overrides = ( + disallowed_overrides.split(PARAMETER_SEPARATOR) + if isinstance(disallowed_overrides, str) + else [disallowed_overrides] + ) + disallowed_subs = dict( + zip(disallowed_formula_placeholders, disallowed_overrides) + ) + symbolic_formula = symbolic_formula.subs(disallowed_subs) + symbolic_formula_inner_parameters = { sp.Symbol(inner_parameter_id): inner_parameter_type for inner_parameter_id, inner_parameter_type in inner_parameters.items() @@ -614,7 +640,10 @@ def validate_observable_data_types(petab_problem: petab.Problem) -> None: other_data_type, other_observables, ) in observables_by_data_type.items(): - if data_type == other_data_type: + if data_type == other_data_type or ( + data_type in CENSORING_TYPES + and other_data_type in CENSORING_TYPES + ): continue if observables & other_observables: raise ValueError( From a23648320c52abca27f1d32d6d7a1045bf769ab9 Mon Sep 17 00:00:00 2001 From: Doresic <85789271+Doresic@users.noreply.github.com> Date: Tue, 12 Nov 2024 14:20:06 +0100 Subject: [PATCH 05/35] Hierarchical: remove inner datas from relative calculator (#1505) * Remove inner rdatas from relative calculator * Remove inner_rdatas from C.py --- pypesto/C.py | 1 - pypesto/hierarchical/relative/calculator.py | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/pypesto/C.py b/pypesto/C.py index 5a0e7438a..e79c90717 100644 --- a/pypesto/C.py +++ b/pypesto/C.py @@ -93,7 +93,6 @@ class EnsembleType(Enum): # HIERARCHICAL SCALING + OFFSET INNER_PARAMETERS = "inner_parameters" -INNER_RDATAS = "inner_rdatas" PARAMETER_TYPE = "parameterType" RELATIVE = "relative" diff --git a/pypesto/hierarchical/relative/calculator.py b/pypesto/hierarchical/relative/calculator.py index cfd74a42f..ceec8ef76 100644 --- a/pypesto/hierarchical/relative/calculator.py +++ b/pypesto/hierarchical/relative/calculator.py @@ -21,7 +21,6 @@ GRAD, HESS, INNER_PARAMETERS, - INNER_RDATAS, RDATAS, RES, SRES, @@ -365,6 +364,6 @@ def calculate_directly( rdatas=rdatas, inner_parameters=inner_parameters, ) - inner_result[INNER_RDATAS] = rdatas + inner_result[RDATAS] = rdatas return inner_result, inner_parameters From e3b3c8ab15f5302d9d6b266225554f1c7ef082b3 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Tue, 12 Nov 2024 15:19:14 +0100 Subject: [PATCH 06/35] Select: problem-specific minimize method for SaCeSS (#1339) * let users supply calibration results * minimize method maker for sacess * doc * user-supplied constructor args * doc how to specify e.g. `max_walltime_s` * allow saving of sacess histories * handle calibrated model via petab select * custom petab select branch * functionality moved to petab-select * change Iterable import * update method for next petab_select version * update for petab select * update for next petab-select version * fix tmpdir * handle no user tmpdir * test SacessMinimizeMethod partially * include fides dependency in select tests * Update pypesto/select/misc.py Co-authored-by: Daniel Weindl --------- Co-authored-by: Dilan Pathirana Co-authored-by: Dilan Pathirana Co-authored-by: Doresic <85789271+Doresic@users.noreply.github.com> Co-authored-by: Doresic Co-authored-by: Daniel Weindl --- pypesto/select/__init__.py | 2 +- pypesto/select/misc.py | 82 +++++++++++++++++++++++++++++++++ pypesto/select/model_problem.py | 17 ++++++- test/select/test_select.py | 30 +++++++++++- tox.ini | 2 +- 5 files changed, 128 insertions(+), 5 deletions(-) diff --git a/pypesto/select/__init__.py b/pypesto/select/__init__.py index eaff283d7..6fc388026 100644 --- a/pypesto/select/__init__.py +++ b/pypesto/select/__init__.py @@ -7,7 +7,7 @@ """ from . import postprocessors -from .misc import model_to_pypesto_problem +from .misc import SacessMinimizeMethod, model_to_pypesto_problem from .problem import Problem try: diff --git a/pypesto/select/misc.py b/pypesto/select/misc.py index 99fb5fae5..870e3339c 100644 --- a/pypesto/select/misc.py +++ b/pypesto/select/misc.py @@ -2,6 +2,7 @@ import logging from collections.abc import Iterable +from pathlib import Path import pandas as pd import petab.v1 as petab @@ -10,7 +11,13 @@ from petab_select import Model, parameter_string_to_value from petab_select.constants import PETAB_PROBLEM +from ..history import Hdf5History from ..objective import Objective +from ..optimize import Optimizer +from ..optimize.ess import ( + SacessOptimizer, + get_default_ess_options, +) from ..petab import PetabImporter from ..problem import Problem @@ -164,3 +171,78 @@ def correct_x_guesses( corrected_x_guess.append(corrected_value) corrected_x_guesses.append(corrected_x_guess) return corrected_x_guesses + + +class SacessMinimizeMethod: + """Create a minimize method for SaCeSS that adapts to each problem. + + When a pyPESTO SaCeSS optimizer is created, it takes the problem + dimension as input. Hence, an optimizer needs to be constructed for + each problem. Objects of this class act like a minimize method for model + selection, but a new problem-specific SaCeSS optimizer will be created + every time a model is minimized. + + Instance attributes correspond to pyPESTO's SaCeSS optimizer, and are + documented there. Extra keyword arguments supplied to the constructor + will be passed on to the constructor of the SaCeSS optimizer, for example, + `max_walltime_s` can be specified in this way. If specified, `tmpdir` will + be treated as a parent directory. + """ + + def __init__( + self, + num_workers: int, + local_optimizer: Optimizer = None, + tmpdir: str | Path | None = None, + save_history: bool = False, + **optimizer_kwargs, + ): + """Construct a minimize-like object.""" + self.num_workers = num_workers + self.local_optimizer = local_optimizer + self.optimizer_kwargs = optimizer_kwargs + self.save_history = save_history + + self.tmpdir = tmpdir + if self.tmpdir is not None: + self.tmpdir = Path(self.tmpdir) + + if self.save_history and self.tmpdir is None: + self.tmpdir = Path.cwd() / "sacess_tmpdir" + + def __call__(self, problem: Problem, model_hash: str, **minimize_options): + """Create then run a problem-specific sacess optimizer.""" + # create optimizer + ess_init_args = get_default_ess_options( + num_workers=self.num_workers, + dim=problem.dim, + ) + for x in ess_init_args: + x["local_optimizer"] = self.local_optimizer + model_tmpdir = None + if self.tmpdir is not None: + model_tmpdir = self.tmpdir / model_hash + model_tmpdir.mkdir(exist_ok=False, parents=True) + + ess = SacessOptimizer( + ess_init_args=ess_init_args, + tmpdir=model_tmpdir, + **self.optimizer_kwargs, + ) + + # optimize + result = ess.minimize( + problem=problem, + **minimize_options, + ) + + if self.save_history: + history_dir = model_tmpdir / "history" + history_dir.mkdir(exist_ok=False, parents=True) + for history_index, history in enumerate(ess.histories): + Hdf5History.from_history( + other=history, + file=history_dir / (str(history_index) + ".hdf5"), + id_=history_index, + ) + return result diff --git a/pypesto/select/model_problem.py b/pypesto/select/model_problem.py index cef94c0e3..fad3487b2 100644 --- a/pypesto/select/model_problem.py +++ b/pypesto/select/model_problem.py @@ -9,7 +9,7 @@ from ..optimize import minimize from ..problem import Problem from ..result import OptimizerResult, Result -from .misc import model_to_pypesto_problem +from .misc import SacessMinimizeMethod, model_to_pypesto_problem OBJECTIVE_CUSTOMIZER_TYPE = Callable[[ObjectiveBase], None] TYPE_POSTPROCESSOR = Callable[["ModelProblem"], None] # noqa: F821 @@ -142,9 +142,16 @@ def __init__( def minimize(self) -> Result: """Optimize the model. - Returns: + Returns + ------- The optimization result. """ + if isinstance(self.minimize_method, SacessMinimizeMethod): + return self.minimize_method( + self.pypesto_problem, + model_hash=self.model.get_hash(), + **self.minimize_options, + ) return self.minimize_method( self.pypesto_problem, **self.minimize_options, @@ -195,6 +202,12 @@ def create_fake_pypesto_result_from_fval( ---------- fval: The objective function value. + evaluation_time: + CPU time taken to compute the objective function value. + + Returns + ------- + The dummy result. """ result = Result() diff --git a/test/select/test_select.py b/test/select/test_select.py index 5ba2d5a3e..bdec69452 100644 --- a/test/select/test_select.py +++ b/test/select/test_select.py @@ -20,7 +20,7 @@ import pypesto.engine import pypesto.select import pypesto.visualize.select -from pypesto.select import model_problem +from pypesto.select import SacessMinimizeMethod, model_problem from pypesto.select.misc import correct_x_guesses model_problem_options = { @@ -496,3 +496,31 @@ def model_to_pypesto_problem( # The custom objective and gradient were returned. assert test_fun == expected_fun assert np.isclose(test_grad, expected_grad).all() + + +def test_sacess_minimize_method(pypesto_select_problem, initial_models): + """Test `SacessMinimizeMethod`. + + Only ensures that the pipeline runs. + """ + predecessor_model = initial_models[1] + + minimize_method = SacessMinimizeMethod( + num_workers=2, + max_walltime_s=1, + ) + + minimize_options = { + "startpoint_method": pypesto.startpoint.UniformStartpoints(), + } + + model_problem_options = { + "minimize_method": minimize_method, + "minimize_options": minimize_options, + } + + pypesto_select_problem.select_to_completion( + model_problem_options=model_problem_options, + method=petab_select.Method.FORWARD, + predecessor_model=predecessor_model, + ) diff --git a/tox.ini b/tox.ini index 6ad9f8719..e9193e0ad 100644 --- a/tox.ini +++ b/tox.ini @@ -120,7 +120,7 @@ description = Test hierarchical optimization module [testenv:select] -extras = test,amici,petab,select +extras = test,amici,petab,select,fides commands = pytest --cov=pypesto --cov-report=xml --cov-append \ test/select From cf56080831129f155514d10573bd848e1abf7c4f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Nov 2024 09:32:47 +0100 Subject: [PATCH 07/35] Fix ESSOptimizer min of empty sequence (#1510) Fixes `ValueError: min() arg is an empty sequence)` that may cause deadlocks. --- pypesto/optimize/ess/ess.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index 1f31df0f5..b38476830 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -487,19 +487,25 @@ def _do_local_search( quality_order = np.argsort(fx_best_children) # compute minimal distance between the best children and all local # optima found so far - min_distances = np.fromiter( - ( - min( - np.linalg.norm( - y_i - - optimizer_result.x[optimizer_result.free_indices] + min_distances = ( + np.fromiter( + ( + min( + np.linalg.norm( + y_i + - optimizer_result.x[ + optimizer_result.free_indices + ] + ) + for optimizer_result in self.local_solutions ) - for optimizer_result in self.local_solutions - ) - for y_i in x_best_children - ), - dtype=np.float64, - count=len(x_best_children), + for y_i in x_best_children + ), + dtype=np.float64, + count=len(x_best_children), + ) + if len(self.local_solutions) + else np.zeros(len(x_best_children)) ) # sort by furthest distance to existing local optima diversity_order = np.argsort(min_distances)[::-1] From e62c2cb94090309a564fda01d546e2ffb2902f12 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Nov 2024 13:47:52 +0100 Subject: [PATCH 08/35] PEtab: Fix warning from `fill_in_parameters` with fixed parameters (#1509) Fixes a warning `RuntimeWarning: The following problem parameters were not used: {[...]} amici.petab.conditions.fill_in_parameters(` in case of fixed parameters. The warning was issued since https://github.com/ICB-DCM/pyPESTO/pull/1493. This also fills in the fixed parameter values in the parameter mapping for the fixed parameters that was missing in https://github.com/ICB-DCM/pyPESTO/pull/1493. --- pypesto/objective/amici/amici.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/pypesto/objective/amici/amici.py b/pypesto/objective/amici/amici.py index 8edcf13fa..6edbe1dde 100644 --- a/pypesto/objective/amici/amici.py +++ b/pypesto/objective/amici/amici.py @@ -183,7 +183,12 @@ def __init__( # (relevant for setting ``plist`` in ExpData later on) self.parameter_mapping = parameter_mapping # parameter mapping independent of fixed parameters + # (i.e., all objective parameters are included as parameter IDs, + # not as their values) self._parameter_mapping_full = copy.deepcopy(parameter_mapping) + # IDs of fixed `Problem` parameters + self._fixed_parameter_ids = [] + # If supported, enable `guess_steadystate` by default. If not # supported, disable by default. If requested but unsupported, raise. if ( @@ -478,8 +483,16 @@ def call_unprocessed( edatas = self.edatas if parameter_mapping is None: parameter_mapping = self.parameter_mapping + # Some parameters may appear estimated in the original compiled model, + # but then are fixed during parameter estimation. These are removed + # from the parameter vector to avoid warnings about unused parameters. + x_dct_free = { + par_id: val + for par_id, val in x_dct.items() + if par_id not in self._fixed_parameter_ids + } ret = self.calculator( - x_dct=x_dct, + x_dct=x_dct_free, sensi_orders=sensi_orders, mode=mode, amici_model=self.amici_model, @@ -679,6 +692,7 @@ def update_from_problem( # and replace the IDs of all fixed parameters by their respective # values. self.parameter_mapping = copy.deepcopy(self._parameter_mapping_full) + self._fixed_parameter_ids = [self.x_ids[i] for i in x_fixed_indices] if not len(x_fixed_indices): return @@ -693,3 +707,9 @@ def update_from_problem( ) in condition_mapping.map_sim_var.items(): if (val := id_to_val.get(mapped_to_par)) is not None: condition_mapping.map_sim_var[model_par] = val + for ( + model_par, + mapped_to_par, + ) in condition_mapping.map_sim_fix.items(): + if (val := id_to_val.get(mapped_to_par)) is not None: + condition_mapping.map_sim_fix[model_par] = val From c5bfd72c22592ca53b1ab7c240067b0ae1bb1f5a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Nov 2024 15:23:58 +0100 Subject: [PATCH 09/35] Doc: extend scatter search documentation (#1506) Add example, background, cross-references, ... --------- Co-authored-by: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> --- pypesto/optimize/ess/ess.py | 120 ++++++++++++++++++++++++++------- pypesto/optimize/ess/refset.py | 14 ++-- pypesto/optimize/ess/sacess.py | 90 ++++++++++++++++++++----- pytest.ini | 1 + 4 files changed, 179 insertions(+), 46 deletions(-) diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index b38476830..2a86e18d9 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -3,11 +3,12 @@ See papers on ESS :footcite:p:`EgeaBal2009,EgeaMar2010`, CESS :footcite:p:`VillaverdeEge2012`, and saCeSS :footcite:p:`PenasGon2017`. """ +from __future__ import annotations import enum import logging import time -from typing import Callable, Optional, Union +from typing import Protocol from warnings import warn import numpy as np @@ -38,15 +39,80 @@ class ESSExitFlag(int, enum.Enum): MAX_TIME = -3 +class OptimizerFactory(Protocol): + def __call__( + self, max_eval: float, max_walltime_s: float + ) -> pypesto.optimize.Optimizer: + """Create a new optimizer instance. + + Parameters + ---------- + max_eval: + Maximum number of objective functions allowed. + max_walltime_s: + Maximum walltime in seconds. + """ + ... + + class ESSOptimizer: """Enhanced Scatter Search (ESS) global optimization. - See papers on ESS :footcite:p:`EgeaBal2009,EgeaMar2010`, - CESS :footcite:p:`VillaverdeEge2012`, and saCeSS :footcite:p:`PenasGon2017`. + Scatter search is a meta-heuristic for global optimization. A set of points + (the reference set, RefSet) is iteratively adapted to explore the parameter + space and to follow promising directions. - .. footbibliography:: + This implementation is based on :footcite:p:`EgeaBal2009,EgeaMar2010`, + but does not implement any constraint handling beyond box constraints. + + The basic steps of ESS are: + + * Initialization: Generate a diverse set of points (RefSet) in the + parameter space. + * Recombination: Generate new points by recombining the RefSet points. + * Improvement: Improve the RefSet by replacing points with better ones. + + The steps are repeated until a stopping criterion is met. + + ESS is gradient-free, unless a gradient-based local optimizer is used + (``local_optimizer``). - .. note: Does not implement any constraint handling beyond box constraints + Hyperparameters + --------------- + + Various hyperparameters control the behavior of ESS. + Initialization is controlled by ``dim_refset`` and ``n_diverse``. + Local optimizations are controlled by ``local_optimizer``, ``local_n1``, + ``local_n2``, and ``balance``. + + Exit criteria + ------------- + + The optimization stops if any of the following criteria are met: + + * The maximum number of iterations is reached (``max_iter``). + * The maximum number of objective function evaluations is reached + (``max_eval``). + * The maximum wall-time is reached (``max_walltime_s``). + + One of these criteria needs to be provided. + Note that the wall-time and function evaluation criteria are not checked + after every single function evaluation, and thus, the actual number of + function evaluations may slightly exceed the given value. + + Parallelization + --------------- + + Objective function evaluations inside :class:`ESSOptimizer` can be + parallelized using multiprocessing or multithreading by passing a value + >1 for ``n_procs`` or ``n_threads``, respectively. + + + .. seealso:: + + :class:`pypesto.optimize.ess.sacess.SacessOptimizer` + + .. footbibliography:: """ def __init__( @@ -57,10 +123,9 @@ def __init__( local_n1: int = 1, local_n2: int = 10, balance: float = 0.5, - local_optimizer: Union[ - "pypesto.optimize.Optimizer", - Callable[..., "pypesto.optimize.Optimizer"], - ] = None, + local_optimizer: pypesto.optimize.Optimizer + | OptimizerFactory + | None = None, max_eval=None, n_diverse: int = None, n_procs=None, @@ -68,7 +133,7 @@ def __init__( max_walltime_s=None, result_includes_refset: bool = False, ): - """Construct new ESS instance. + r"""Construct new ESS instance. For plausible values of hyperparameters, see :footcite:t:`VillaverdeEge2012`. @@ -81,10 +146,11 @@ def __init__( Maximum number of ESS iterations. local_n1: Minimum number of iterations before first local search. + Ignored if ``local_optimizer=None``. local_n2: Minimum number of iterations between consecutive local searches. Maximally one local search per performed in each - iteration. + iteration. Ignored if ``local_optimizer=None``. local_optimizer: Local optimizer for refinement, or a callable that creates an :class:`pypesto.optimize.Optimizer` or ``None`` to skip local searches. @@ -104,8 +170,14 @@ def __init__( optimizations and other simulations, and thus, may be exceeded by the duration of a local search. balance: - Quality vs diversity balancing factor [0, 1]; - 0 = only quality; 1 = only diversity + Quality vs. diversity balancing factor with + :math:`0 \leq balance \leq 1`; ``0`` = only quality, + ``1`` = only diversity. + Affects the choice of starting points for local searches. I.e., + whether local optimization should focus on improving the best + solutions found so far (quality), or on exploring new regions of + the parameter space (diversity). + Ignored if ``local_optimizer=None``. n_procs: Number of parallel processes to use for parallel function evaluation. Mutually exclusive with `n_threads`. @@ -144,8 +216,8 @@ def __init__( raise ValueError( "`n_procs` and `n_threads` are mutually exclusive." ) - self.n_procs: Optional[int] = n_procs - self.n_threads: Optional[int] = n_threads + self.n_procs: int | None = n_procs + self.n_threads: int | None = n_threads self.balance: float = balance # After how many iterations a stagnated solution is to be replaced by # a random one. Default value taken from [EgeaMar2010]_ @@ -162,9 +234,9 @@ def __init__( def _initialize(self): """(Re-)Initialize.""" # RefSet - self.refset: Optional[RefSet] = None + self.refset: RefSet | None = None # Overall best parameters found so far - self.x_best: Optional[np.array] = None + self.x_best: np.ndarray | None = None # Overall best function value found so far self.fx_best: float = np.inf # Results from local searches (only those with finite fval) @@ -177,15 +249,15 @@ def _initialize(self): # Whether self.x_best has changed in the current iteration self.x_best_has_changed: bool = False self.exit_flag: ESSExitFlag = ESSExitFlag.DID_NOT_RUN - self.evaluator: Optional[FunctionEvaluator] = None - self.starttime: Optional[float] = None + self.evaluator: FunctionEvaluator | None = None + self.starttime: float | None = None self.history: MemoryHistory = MemoryHistory() def _initialize_minimize( self, problem: Problem = None, startpoint_method: StartpointMethod = None, - refset: Optional[RefSet] = None, + refset: RefSet | None = None, ): """Initialize for optimizations. @@ -242,7 +314,7 @@ def minimize( self, problem: Problem = None, startpoint_method: StartpointMethod = None, - refset: Optional[RefSet] = None, + refset: RefSet | None = None, ) -> pypesto.Result: """Minimize the given objective. @@ -384,7 +456,7 @@ def _get_remaining_eval(self): return np.inf return self.max_eval - self.evaluator.n_eval - def _combine_solutions(self) -> tuple[np.array, np.array]: + def _combine_solutions(self) -> tuple[np.ndarray, np.ndarray]: """Combine solutions and evaluate. Creates the next generation from the RefSet by pair-wise combination @@ -418,7 +490,7 @@ def _combine_solutions(self) -> tuple[np.array, np.array]: break return y, fy - def _combine(self, i, j) -> np.array: + def _combine(self, i, j) -> np.ndarray: """Combine RefSet members ``i`` and ``j``. Samples a new point from a biased hyper-rectangle derived from the @@ -463,7 +535,7 @@ def _combine(self, i, j) -> np.array: ) def _do_local_search( - self, x_best_children: np.array, fx_best_children: np.array + self, x_best_children: np.ndarray, fx_best_children: np.ndarray ) -> None: """ Perform a local search to refine the next generation. diff --git a/pypesto/optimize/ess/refset.py b/pypesto/optimize/ess/refset.py index 0e3cff403..dfd8a4d42 100644 --- a/pypesto/optimize/ess/refset.py +++ b/pypesto/optimize/ess/refset.py @@ -29,8 +29,8 @@ def __init__( self, dim: int, evaluator: FunctionEvaluator, - x: Optional[np.array] = None, - fx: Optional[np.array] = None, + x: Optional[np.ndarray] = None, + fx: Optional[np.ndarray] = None, ): """Construct. @@ -65,7 +65,7 @@ def __init__( self.fx = fx self.n_stuck = np.zeros(shape=[dim]) - self.attributes: dict[Any, np.array] = {} + self.attributes: dict[Any, np.ndarray] = {} def __repr__(self): fx = ( @@ -97,7 +97,9 @@ def initialize_random( x_diverse, fx_diverse = self.evaluator.multiple_random(n_diverse) self.initialize_from_array(x_diverse=x_diverse, fx_diverse=fx_diverse) - def initialize_from_array(self, x_diverse: np.array, fx_diverse: np.array): + def initialize_from_array( + self, x_diverse: np.ndarray, fx_diverse: np.ndarray + ): """Create an initial reference set using the provided points. Populate half of the RefSet using the best given solutions and fill the @@ -168,7 +170,7 @@ def normalize(x): x[j], self.fx[j] = self.evaluator.single_random() self.sort() - def update(self, i: int, x: np.array, fx: float): + def update(self, i: int, x: np.ndarray, fx: float): """Update a RefSet entry.""" self.x[i] = x self.fx[i] = fx @@ -179,7 +181,7 @@ def replace_by_random(self, i: int): self.x[i], self.fx[i] = self.evaluator.single_random() self.n_stuck[i] = 0 - def add_attribute(self, name: str, values: np.array): + def add_attribute(self, name: str, values: np.ndarray): """ Add an attribute array to the refset members. diff --git a/pypesto/optimize/ess/sacess.py b/pypesto/optimize/ess/sacess.py index c54ef0ad6..113310f25 100644 --- a/pypesto/optimize/ess/sacess.py +++ b/pypesto/optimize/ess/sacess.py @@ -41,8 +41,10 @@ class SacessOptimizer: """SACESS optimizer. - A shared-memory-based implementation of the SaCeSS algorithm presented in - :footcite:t:`PenasGon2017`. Multiple processes (`workers`) run + A shared-memory-based implementation of the + `Self-Adaptive Cooperative Enhanced Scatter Search` (SaCeSS) algorithm + presented in :footcite:t:`PenasGon2017`. This is a meta-heuristic for + global optimization. Multiple processes (`workers`) run :class:`enhanced scatter searches (ESSs) ` in parallel. After each ESS iteration, depending on the outcome, there is a chance of exchanging good parameters, and changing ESS hyperparameters to those of @@ -51,6 +53,41 @@ class SacessOptimizer: :class:`SacessOptimizer` can be used with or without a local optimizer, but it is highly recommended to use one. + A basic example using :class:`SacessOptimizer` to minimize the Rosenbrock + function: + + >>> from pypesto.optimize import SacessOptimizer + >>> from pypesto.problem import Problem + >>> from pypesto.objective import Objective + >>> import scipy as sp + >>> import numpy as np + >>> import logging + >>> # Define some test Problem + >>> objective = Objective( + ... fun=sp.optimize.rosen, + ... grad=sp.optimize.rosen_der, + ... hess=sp.optimize.rosen_hess, + ... ) + >>> dim = 6 + >>> problem = Problem( + ... objective=objective, + ... lb=-5 * np.ones((dim, 1)), + ... ub=5 * np.ones((dim, 1)), + ... ) + >>> # Create and run the optimizer + >>> sacess = SacessOptimizer( + ... num_workers=2, + ... max_walltime_s=5, + ... sacess_loglevel=logging.WARNING + ... ) + >>> result = sacess.minimize(problem) + + .. seealso:: + + :class:`pypesto.optimize.ess.ess.ESSOptimizer` + + References + ---------- .. footbibliography:: Attributes @@ -83,11 +120,25 @@ def __init__( process. I.e., the length of this list is the number of ESSs. Ideally, this list contains some more conservative and some more aggressive configurations. - Resource limits such as ``max_eval`` apply to a single CESS + Resource limits such as ``max_eval`` apply to a single ESS iteration, not to the full search. Mutually exclusive with ``num_workers``. + Recommended default settings can be obtained from - :func:`get_default_ess_options`. + :func:`get_default_ess_options`. For example, to run + :class:`SacessOptimizer` without a local optimizer, use: + + >>> from pypesto.optimize.ess import get_default_ess_options + >>> ess_init_args = get_default_ess_options( + ... num_workers=12, + ... dim=10, # usually problem.dim + ... local_optimizer=False, + ... ) + >>> ess_init_args # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE + [{'dim_refset': 5, 'balance': 0.0, 'local_n1': 1, 'local_n2': 1}, + ... + {'dim_refset': 7, 'balance': 1.0, 'local_n1': 4, 'local_n2': 4}] + num_workers: Number of workers to be used. If this argument is given, (different) default ESS settings will be used for each worker. @@ -112,9 +163,10 @@ def __init__( parallel have a unique `tmpdir`. mp_start_method: The start method for the multiprocessing context. - See :mod:`multiprocessing` for details. + See :mod:`multiprocessing` for details. Running `SacessOptimizer` + under Jupyter may require ``mp_start_method="fork"``. options: - Further optimizer hyperparameters. + Further optimizer hyperparameters, see :class:`SacessOptions`. """ if (num_workers is None and ess_init_args is None) or ( num_workers is not None and ess_init_args is not None @@ -178,11 +230,12 @@ def minimize( Returns ------- - Result object with optimized parameters in - :attr:`pypesto.Result.optimize_result`. - Results are sorted by objective. At least the best parameters are - included. Additional results may be included - this is subject to - change. + _: + Result object with optimized parameters in + :attr:`pypesto.Result.optimize_result`. + Results are sorted by objective. At least the best parameters are + included. Additional results may be included - this is subject to + change. """ if startpoint_method is not None: warn( @@ -392,7 +445,7 @@ def __init__( self._logger = logging.getLogger() self._result_queue = shmem_manager.Queue() - def get_best_solution(self) -> tuple[np.array, float]: + def get_best_solution(self) -> tuple[np.ndarray, float]: """Get the best objective value and corresponding parameters.""" with self._lock: return np.array(self._best_known_x), self._best_known_fx.value @@ -416,7 +469,7 @@ def reconfigure_worker(self, worker_idx: int) -> dict: def submit_solution( self, - x: np.array, + x: np.ndarray, fx: float, sender_idx: int, elapsed_time_s: float, @@ -703,7 +756,7 @@ def _maybe_adapt(self, problem: Problem): f"neval: {self._neval} <= {problem.dim * self._options.adaptation_min_evals}." ) - def maybe_update_best(self, x: np.array, fx: float): + def maybe_update_best(self, x: np.ndarray, fx: float): """Maybe update the best known solution and send it to the manager.""" rel_change = ( abs((fx - self._best_known_fx) / fx) if fx != 0 else np.nan @@ -743,7 +796,7 @@ def maybe_update_best(self, x: np.array, fx: float): ) @staticmethod - def replace_solution(refset: RefSet, x: np.array, fx: float): + def replace_solution(refset: RefSet, x: np.ndarray, fx: float): """Replace the global refset member by the given solution.""" # [PenasGon2017]_ page 8, top if "cooperative_solution" not in refset.attributes: @@ -841,6 +894,11 @@ def get_default_ess_options( or a :obj:`Callable` returning an optimizer instance. The latter can be used to propagate walltime limits to the local optimizers. See :meth:`SacessFidesFactory.__call__` for an example. + The current default optimizer assumes that the optimized objective + function can provide its gradient. If this is not the case, the user + should provide a different local optimizer or consider using + :class:`pypesto.objective.finite_difference.FD` to approximate the + gradient using finite differences. """ min_dimrefset = 5 @@ -1086,7 +1144,7 @@ class SacessWorkerResult: Exit flag of the optimization process. """ - x: np.array + x: np.ndarray fx: float n_eval: int n_iter: int diff --git a/pytest.ini b/pytest.ini index 7571948eb..c852c729e 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,3 +1,4 @@ [pytest] +addopts = "--doctest-modules" filterwarnings = ignore:.*inspect.getargspec\(\) is deprecated.*:DeprecationWarning From 46115860f4ec4bcf88273f7d64d6950b0b1e2f1a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Nov 2024 10:50:05 +0100 Subject: [PATCH 10/35] Fix handling of PEtab fixed parameters (#1514) Follow-up to #1493 and #1509, which only filled in PEtab-fixed parameters for simulation but not preequilibration in the parameter mapping. Adding Weber_BMC2015 to the PEtab tests, where this issue popped up. --- pypesto/objective/amici/amici.py | 6 ++++++ test/petab/test_petab_import.py | 11 ++++++----- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/pypesto/objective/amici/amici.py b/pypesto/objective/amici/amici.py index 6edbe1dde..d9fd371a3 100644 --- a/pypesto/objective/amici/amici.py +++ b/pypesto/objective/amici/amici.py @@ -713,3 +713,9 @@ def update_from_problem( ) in condition_mapping.map_sim_fix.items(): if (val := id_to_val.get(mapped_to_par)) is not None: condition_mapping.map_sim_fix[model_par] = val + for ( + model_par, + mapped_to_par, + ) in condition_mapping.map_preeq_fix.items(): + if (val := id_to_val.get(mapped_to_par)) is not None: + condition_mapping.map_preeq_fix[model_par] = val diff --git a/test/petab/test_petab_import.py b/test/petab/test_petab_import.py index 1ce935880..d4a39f6d5 100644 --- a/test/petab/test_petab_import.py +++ b/test/petab/test_petab_import.py @@ -36,12 +36,13 @@ def setUpClass(cls): cls.obj_edatas = [] def test_0_import(self): - for model_name in ["Zheng_PNAS2012", "Boehm_JProteomeRes2014"]: + for model_name in [ + "Zheng_PNAS2012", + "Boehm_JProteomeRes2014", + "Weber_BMC2015", + ]: # test yaml import for one model: - yaml_config = os.path.join( - models.MODELS_DIR, model_name, model_name + ".yaml" - ) - petab_problem = petab.Problem.from_yaml(yaml_config) + petab_problem = models.get_problem(model_name) self.petab_problems.append(petab_problem) def test_1_compile(self): From fae775515e0b2f3f1d811d52b5732023482d0a02 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 26 Nov 2024 09:54:01 +0100 Subject: [PATCH 11/35] Exclude nlopt==2.9.0 (#1519) We want to avoid nlopt-python 2.9.0 because it does not match nlopt v2.9.0 and is affected by https://github.com/stevengj/nlopt/issues/575. Closes #1511. --- setup.cfg | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index bc28ae610..d9eb70919 100644 --- a/setup.cfg +++ b/setup.cfg @@ -105,7 +105,8 @@ ipopt = dlib = dlib >= 19.19.0 nlopt = - nlopt >= 2.6.2 + # != 2.9.0: https://github.com/stevengj/nlopt/issues/575 + nlopt >= 2.6.2, != 2.9.0 pyswarm = pyswarm >= 0.6 cma = From 6c049f599bbda2026ae1edfef76fb1df87d84358 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 26 Nov 2024 16:25:58 +0100 Subject: [PATCH 12/35] Fix test_deepcopy_objective, test_preeq_guesses (#1521) * Fix test_deepcopy_objective (Fixes #1520) * Fix test_preeq_guesses Failures were related to changes in amici's default steady-state settings. --- test/base/test_engine.py | 9 ++++++++- test/petab/test_amici_objective.py | 15 ++++++++++++--- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/test/base/test_engine.py b/test/base/test_engine.py index 6db8e79c3..1481c4ca1 100644 --- a/test/base/test_engine.py +++ b/test/base/test_engine.py @@ -86,7 +86,14 @@ def test_deepcopy_objective(): ) ) factory = petab_importer.create_objective_creator() - objective = factory.create_objective() + amici_model = factory.create_model() + amici_model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) + amici_model.setSteadyStateComputationMode( + amici.SteadyStateComputationMode.integrateIfNewtonFails + ) + objective = factory.create_objective(model=amici_model) objective.amici_solver.setSensitivityMethod( amici.SensitivityMethod_adjoint diff --git a/test/petab/test_amici_objective.py b/test/petab/test_amici_objective.py index 274962fc1..5ff4512b4 100644 --- a/test/petab/test_amici_objective.py +++ b/test/petab/test_amici_objective.py @@ -86,12 +86,21 @@ def test_preeq_guesses(): importer = pypesto.petab.PetabImporter.from_yaml( os.path.join(models.MODELS_DIR, model_name, model_name + ".yaml") ) - problem = importer.create_problem() - obj = problem.objective - obj.amici_solver.setNewtonMaxSteps(0) + obj_creator = importer.create_objective_creator() + amici_model = obj_creator.create_model() + amici_model.setSteadyStateComputationMode( + amici.SteadyStateComputationMode.integrateIfNewtonFails + ) + amici_model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) + obj = obj_creator.create_objective(model=amici_model) + problem = importer.create_problem(objective=obj) obj.amici_model.setSteadyStateSensitivityMode( amici.SteadyStateSensitivityMode.integrationOnly ) + obj = problem.objective + obj.amici_solver.setNewtonMaxSteps(0) obj.amici_solver.setAbsoluteTolerance(1e-12) obj.amici_solver.setRelativeTolerance(1e-12) From 6eb424520e039309d6087cadcefe6162bc2c174a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 26 Nov 2024 16:47:34 +0100 Subject: [PATCH 13/35] Update references (#1491) * ArmisteadHoe2024 * LakrisenkoPat2024 * .. --- doc/using_pypesto.bib | 161 ++++++++++++++++++++++++------------------ 1 file changed, 92 insertions(+), 69 deletions(-) diff --git a/doc/using_pypesto.bib b/doc/using_pypesto.bib index c4c85924d..40ca3247f 100644 --- a/doc/using_pypesto.bib +++ b/doc/using_pypesto.bib @@ -14,63 +14,67 @@ @Article{FalcoCoh2023 } @Article{LakrisenkoSta2023, - author = {Lakrisenko, Polina AND Stapor, Paul AND Grein, Stephan AND Paszkowski, Łukasz AND Pathirana, Dilan AND Fröhlich, Fabian AND Lines, Glenn Terje AND Weindl, Daniel AND Hasenauer, Jan}, - journal = {PLOS Computational Biology}, - title = {Efficient computation of adjoint sensitivities at steady-state in ODE models of biochemical reaction networks}, - year = {2023}, - month = {01}, - number = {1}, - pages = {1-19}, - volume = {19}, - abstract = {Dynamical models in the form of systems of ordinary differential equations have become a standard tool in systems biology. Many parameters of such models are usually unknown and have to be inferred from experimental data. Gradient-based optimization has proven to be effective for parameter estimation. However, computing gradients becomes increasingly costly for larger models, which are required for capturing the complex interactions of multiple biochemical pathways. Adjoint sensitivity analysis has been pivotal for working with such large models, but methods tailored for steady-state data are currently not available. We propose a new adjoint method for computing gradients, which is applicable if the experimental data include steady-state measurements. The method is based on a reformulation of the backward integration problem to a system of linear algebraic equations. The evaluation of the proposed method using real-world problems shows a speedup of total simulation time by a factor of up to 4.4. Our results demonstrate that the proposed approach can achieve a substantial improvement in computation time, in particular for large-scale models, where computational efficiency is critical.}, - creationdate = {2023-01-26T11:19:52}, - doi = {10.1371/journal.pcbi.1010783}, - publisher = {Public Library of Science}, + author = {Lakrisenko, Polina and Stapor, Paul and Grein, Stephan and Paszkowski, Łukasz and Pathirana, Dilan and Fröhlich, Fabian and Lines, Glenn Terje and Weindl, Daniel and Hasenauer, Jan}, + journal = {PLOS Computational Biology}, + title = {Efficient computation of adjoint sensitivities at steady-state in {ODE} models of biochemical reaction networks}, + year = {2023}, + month = {01}, + number = {1}, + pages = {1-19}, + volume = {19}, + abstract = {Dynamical models in the form of systems of ordinary differential equations have become a standard tool in systems biology. Many parameters of such models are usually unknown and have to be inferred from experimental data. Gradient-based optimization has proven to be effective for parameter estimation. However, computing gradients becomes increasingly costly for larger models, which are required for capturing the complex interactions of multiple biochemical pathways. Adjoint sensitivity analysis has been pivotal for working with such large models, but methods tailored for steady-state data are currently not available. We propose a new adjoint method for computing gradients, which is applicable if the experimental data include steady-state measurements. The method is based on a reformulation of the backward integration problem to a system of linear algebraic equations. The evaluation of the proposed method using real-world problems shows a speedup of total simulation time by a factor of up to 4.4. Our results demonstrate that the proposed approach can achieve a substantial improvement in computation time, in particular for large-scale models, where computational efficiency is critical.}, + creationdate = {2023-01-26T11:19:52}, + doi = {10.1371/journal.pcbi.1010783}, + modificationdate = {2024-11-08T18:28:07}, + publisher = {Public Library of Science}, } @Article{SchmiesterSch2021, - author = {Schmiester, Leonard AND Schälte, Yannik AND Bergmann, Frank T. AND Camba, Tacio AND Dudkin, Erika AND Egert, Janine AND Fröhlich, Fabian AND Fuhrmann, Lara AND Hauber, Adrian L. AND Kemmer, Svenja AND Lakrisenko, Polina AND Loos, Carolin AND Merkt, Simon AND Müller, Wolfgang AND Pathirana, Dilan AND Raimúndez, Elba AND Refisch, Lukas AND Rosenblatt, Marcus AND Stapor, Paul L. AND Städter, Philipp AND Wang, Dantong AND Wieland, Franz-Georg AND Banga, Julio R. AND Timmer, Jens AND Villaverde, Alejandro F. AND Sahle, Sven AND Kreutz, Clemens AND Hasenauer, Jan AND Weindl, Daniel}, - journal = {PLOS Computational Biology}, - title = {PEtab—Interoperable specification of parameter estimation problems in systems biology}, - year = {2021}, - month = {01}, - number = {1}, - pages = {1-10}, - volume = {17}, - abstract = {Author summary Parameter estimation is a common and crucial task in modeling, as many models depend on unknown parameters which need to be inferred from data. There exist various tools for tasks like model development, model simulation, optimization, or uncertainty analysis, each with different capabilities and strengths. In order to be able to easily combine tools in an interoperable manner, but also to make results accessible and reusable for other researchers, it is valuable to define parameter estimation problems in a standardized form. Here, we introduce PEtab, a parameter estimation problem definition format which integrates with established systems biology standards for model and data specification. As the novel format is already supported by eight software tools with hundreds of users in total, we expect it to be of great use and impact in the community, both for modeling and algorithm development.}, - creationdate = {2023-01-26T11:30:40}, - doi = {10.1371/journal.pcbi.1008646}, - publisher = {Public Library of Science}, - timestamp = {2021-01-30}, + author = {Schmiester, Leonard and Schälte, Yannik and Bergmann, Frank T. and Camba, Tacio and Dudkin, Erika and Egert, Janine and Fröhlich, Fabian and Fuhrmann, Lara and Hauber, Adrian L. and Kemmer, Svenja and Lakrisenko, Polina and Loos, Carolin and Merkt, Simon and Müller, Wolfgang and Pathirana, Dilan and Raimúndez, Elba and Refisch, Lukas and Rosenblatt, Marcus and Stapor, Paul L. and Städter, Philipp and Wang, Dantong and Wieland, Franz-Georg and Banga, Julio R. and Timmer, Jens 2and Villaverde, Alejandro F. and Sahle, Sven and Kreutz, Clemens and Hasenauer, Jan and Weindl, Daniel}, + journal = {PLOS Computational Biology}, + title = {{PEtab}—Interoperable specification of parameter estimation problems in systems biology}, + year = {2021}, + month = {01}, + number = {1}, + pages = {1-10}, + volume = {17}, + abstract = {Author summary Parameter estimation is a common and crucial task in modeling, as many models depend on unknown parameters which need to be inferred from data. There exist various tools for tasks like model development, model simulation, optimization, or uncertainty analysis, each with different capabilities and strengths. In order to be able to easily combine tools in an interoperable manner, but also to make results accessible and reusable for other researchers, it is valuable to define parameter estimation problems in a standardized form. Here, we introduce PEtab, a parameter estimation problem definition format which integrates with established systems biology standards for model and data specification. As the novel format is already supported by eight software tools with hundreds of users in total, we expect it to be of great use and impact in the community, both for modeling and algorithm development.}, + creationdate = {2023-01-26T11:30:40}, + doi = {10.1371/journal.pcbi.1008646}, + modificationdate = {2024-11-08T18:27:03}, + publisher = {Public Library of Science}, + timestamp = {2021-01-30}, } @Article{MishraWan2023, - author = {Shekhar Mishra and Ziyu Wang and Michael J. Volk and Huimin Zhao}, - journal = {Metabolic Engineering}, - title = {Design and application of a kinetic model of lipid metabolism in Saccharomyces cerevisiae}, - year = {2023}, - issn = {1096-7176}, - pages = {12-18}, - volume = {75}, - abstract = {Lipid biosynthesis plays a vital role in living cells and has been increasingly engineered to overproduce various lipid-based chemicals. However, owing to the tightly constrained and interconnected nature of lipid biosynthesis, both understanding and engineering of lipid metabolism remain challenging, even with the help of mathematical models. Here we report the development of a kinetic metabolic model of lipid metabolism in Saccharomyces cerevisiae that integrates fatty acid biosynthesis, glycerophospholipid metabolism, sphingolipid metabolism, storage lipids, lumped sterol synthesis, and the synthesis and transport of relevant target-chemicals, such as fatty acids and fatty alcohols. The model was trained on lipidomic data of a reference S. cerevisiae strain, single knockout mutants, and lipid overproduction strains reported in literature. The model was used to design mutants for fatty alcohol overproduction and the lipidomic analysis of the resultant mutant strains coupled with model-guided hypothesis led to discovery of a futile cycle in the triacylglycerol biosynthesis pathway. In addition, the model was used to explain successful and unsuccessful mutant designs in metabolic engineering literature. Thus, this kinetic model of lipid metabolism can not only enable the discovery of new phenomenon in lipid metabolism but also the engineering of mutant strains for overproduction of lipids.}, - creationdate = {2023-01-26T11:31:17}, - doi = {https://doi.org/10.1016/j.ymben.2022.11.003}, - keywords = {Lipid metabolism, Kinetic model, Free fatty acid, Fatty alcohol}, + author = {Shekhar Mishra and Ziyu Wang and Michael J. Volk and Huimin Zhao}, + journal = {Metabolic Engineering}, + title = {Design and application of a kinetic model of lipid metabolism in Saccharomyces cerevisiae}, + year = {2023}, + issn = {1096-7176}, + pages = {12-18}, + volume = {75}, + abstract = {Lipid biosynthesis plays a vital role in living cells and has been increasingly engineered to overproduce various lipid-based chemicals. However, owing to the tightly constrained and interconnected nature of lipid biosynthesis, both understanding and engineering of lipid metabolism remain challenging, even with the help of mathematical models. Here we report the development of a kinetic metabolic model of lipid metabolism in Saccharomyces cerevisiae that integrates fatty acid biosynthesis, glycerophospholipid metabolism, sphingolipid metabolism, storage lipids, lumped sterol synthesis, and the synthesis and transport of relevant target-chemicals, such as fatty acids and fatty alcohols. The model was trained on lipidomic data of a reference S. cerevisiae strain, single knockout mutants, and lipid overproduction strains reported in literature. The model was used to design mutants for fatty alcohol overproduction and the lipidomic analysis of the resultant mutant strains coupled with model-guided hypothesis led to discovery of a futile cycle in the triacylglycerol biosynthesis pathway. In addition, the model was used to explain successful and unsuccessful mutant designs in metabolic engineering literature. Thus, this kinetic model of lipid metabolism can not only enable the discovery of new phenomenon in lipid metabolism but also the engineering of mutant strains for overproduction of lipids.}, + creationdate = {2023-01-26T11:31:17}, + doi = {10.1016/j.ymben.2022.11.003}, + keywords = {Lipid metabolism, Kinetic model, Free fatty acid, Fatty alcohol}, + modificationdate = {2024-11-08T18:25:04}, } @Article{FroehlichSor2022, - author = {Fröhlich, Fabian AND Sorger, Peter K.}, - journal = {PLOS Computational Biology}, - title = {Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models}, - year = {2022}, - month = {07}, - number = {7}, - pages = {1-28}, - volume = {18}, - abstract = {Ordinary differential equation (ODE) models are widely used to study biochemical reactions in cellular networks since they effectively describe the temporal evolution of these networks using mass action kinetics. The parameters of these models are rarely known a priori and must instead be estimated by calibration using experimental data. Optimization-based calibration of ODE models on is often challenging, even for low-dimensional problems. Multiple hypotheses have been advanced to explain why biochemical model calibration is challenging, including non-identifiability of model parameters, but there are few comprehensive studies that test these hypotheses, likely because tools for performing such studies are also lacking. Nonetheless, reliable model calibration is essential for uncertainty analysis, model comparison, and biological interpretation. We implemented an established trust-region method as a modular Python framework (fides) to enable systematic comparison of different approaches to ODE model calibration involving a variety of Hessian approximation schemes. We evaluated fides on a recently developed corpus of biologically realistic benchmark problems for which real experimental data are available. Unexpectedly, we observed high variability in optimizer performance among different implementations of the same mathematical instructions (algorithms). Analysis of possible sources of poor optimizer performance identified limitations in the widely used Gauss-Newton, BFGS and SR1 Hessian approximation schemes. We addressed these drawbacks with a novel hybrid Hessian approximation scheme that enhances optimizer performance and outperforms existing hybrid approaches. When applied to the corpus of test models, we found that fides was on average more reliable and efficient than existing methods using a variety of criteria. We expect fides to be broadly useful for ODE constrained optimization problems in biochemical models and to be a foundation for future methods development.}, - creationdate = {2023-01-26T11:31:44}, - doi = {10.1371/journal.pcbi.1010322}, - publisher = {Public Library of Science}, + author = {Fröhlich, Fabian and Sorger, Peter K.}, + journal = {PLOS Computational Biology}, + title = {Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models}, + year = {2022}, + month = {07}, + number = {7}, + pages = {1-28}, + volume = {18}, + abstract = {Ordinary differential equation (ODE) models are widely used to study biochemical reactions in cellular networks since they effectively describe the temporal evolution of these networks using mass action kinetics. The parameters of these models are rarely known a priori and must instead be estimated by calibration using experimental data. Optimization-based calibration of ODE models on is often challenging, even for low-dimensional problems. Multiple hypotheses have been advanced to explain why biochemical model calibration is challenging, including non-identifiability of model parameters, but there are few comprehensive studies that test these hypotheses, likely because tools for performing such studies are also lacking. Nonetheless, reliable model calibration is essential for uncertainty analysis, model comparison, and biological interpretation. We implemented an established trust-region method as a modular Python framework (fides) to enable systematic comparison of different approaches to ODE model calibration involving a variety of Hessian approximation schemes. We evaluated fides on a recently developed corpus of biologically realistic benchmark problems for which real experimental data are available. Unexpectedly, we observed high variability in optimizer performance among different implementations of the same mathematical instructions (algorithms). Analysis of possible sources of poor optimizer performance identified limitations in the widely used Gauss-Newton, BFGS and SR1 Hessian approximation schemes. We addressed these drawbacks with a novel hybrid Hessian approximation scheme that enhances optimizer performance and outperforms existing hybrid approaches. When applied to the corpus of test models, we found that fides was on average more reliable and efficient than existing methods using a variety of criteria. We expect fides to be broadly useful for ODE constrained optimization problems in biochemical models and to be a foundation for future methods development.}, + creationdate = {2023-01-26T11:31:44}, + doi = {10.1371/journal.pcbi.1010322}, + modificationdate = {2024-11-08T18:20:34}, + publisher = {Public Library of Science}, } @Article{FroehlichGer2022, @@ -242,7 +246,7 @@ @Article{ArrudaSch2023 @Article{MerktAli2024, author = {Merkt, Simon and Ali, Solomon and Gudina, Esayas Kebede and Adissu, Wondimagegn and Gize, Addisu and Muenchhoff, Maximilian and Graf, Alexander and Krebs, Stefan and Elsbernd, Kira and Kisch, Rebecca and Betizazu, Sisay Sirgu and Fantahun, Bereket and Bekele, Delayehu and Rubio-Acero, Raquel and Gashaw, Mulatu and Girma, Eyob and Yilma, Daniel and Zeynudin, Ahmed and Paunovic, Ivana and Hoelscher, Michael and Blum, Helmut and Hasenauer, Jan and Kroidl, Arne and Wieser, Andreas}, journal = {Nature Communications}, - title = {Long-term monitoring of SARS-CoV-2 seroprevalence and variants in Ethiopia provides prediction for immunity and cross-immunity}, + title = {Long-term monitoring of {SARS-CoV-2} seroprevalence and variants in {Ethiopia} provides prediction for immunity and cross-immunity}, year = {2024}, issn = {2041-1723}, month = apr, @@ -250,7 +254,7 @@ @Article{MerktAli2024 volume = {15}, creationdate = {2024-04-29T08:32:16}, doi = {10.1038/s41467-024-47556-2}, - modificationdate = {2024-04-29T08:32:16}, + modificationdate = {2024-11-08T18:27:48}, publisher = {Springer Science and Business Media LLC}, } @@ -282,29 +286,18 @@ @Article{HoepflAlb2024 modificationdate = {2024-05-16T07:58:55}, } -@Misc{LakrisenkoPat2024, - author = {Polina Lakrisenko and Dilan Pathirana and Daniel Weindl and Jan Hasenauer}, - title = {Exploration of methods for computing sensitivities in ODE models at dynamic and steady states}, - year = {2024}, - archiveprefix = {arXiv}, - creationdate = {2024-05-30T09:47:51}, - eprint = {2405.16524}, - modificationdate = {2024-05-30T09:47:51}, - primaryclass = {q-bio.QM}, -} - @Article{PhilippsKoe2024, author = {Maren Philipps and Antonia Körner and Jakob Vanhoefer and Dilan Pathirana and Jan Hasenauer}, + journal = {IFAC-PapersOnLine}, title = {Non-Negative Universal Differential Equations With Applications in Systems Biology}, year = {2024}, - journal = {IFAC-PapersOnLine}, - volume = {58}, + issn = {2405-8963}, number = {23}, pages = {25-30}, - issn = {2405-8963}, - doi = {https://doi.org/10.1016/j.ifacol.2024.10.005}, - url = {https://www.sciencedirect.com/science/article/pii/S2405896324017518}, - abstract = {Universal differential equations (UDEs) leverage the respective advantages of mechanistic models and artificial neural networks and combine them into one dynamic model. However, these hybrid models can suffer from unrealistic solutions, such as negative values for biochemical quantities. We present non-negative UDE (nUDEs), a constrained UDE variant that guarantees non-negative values. Furthermore, we explore regularisation techniques to improve generalisation and interpretability of UDEs.} + volume = {58}, + abstract = {Universal differential equations (UDEs) leverage the respective advantages of mechanistic models and artificial neural networks and combine them into one dynamic model. However, these hybrid models can suffer from unrealistic solutions, such as negative values for biochemical quantities. We present non-negative UDE (nUDEs), a constrained UDE variant that guarantees non-negative values. Furthermore, we explore regularisation techniques to improve generalisation and interpretability of UDEs.}, + doi = {10.1016/j.ifacol.2024.10.005}, + modificationdate = {2024-11-08T18:25:20}, } @Article{SchmiesterBra2024, @@ -319,8 +312,7 @@ @Article{SchmiesterBra2024 creationdate = {2024-08-01T09:44:04}, doi = {10.1158/1078-0432.CCR-24-0244}, eprint = {https://aacrjournals.org/clincancerres/article-pdf/doi/10.1158/1078-0432.CCR-24-0244/3478451/ccr-24-0244.pdf}, - modificationdate = {2024-08-01T09:44:04}, - url = {https://doi.org/10.1158/1078-0432.CCR-24-0244}, + modificationdate = {2024-11-08T18:17:14}, } @InProceedings{JacksonCha2023, @@ -339,4 +331,35 @@ @InProceedings{JacksonCha2023 modificationdate = {2024-09-06T15:49:47}, } +@Article{ArmisteadHoe2024, + author = {Armistead, Joy and Höpfl, Sebastian and Goldhausen, Pierre and Müller-Hartmann, Andrea and Fahle, Evelin and Hatzold, Julia and Franzen, Rainer and Brodesser, Susanne and Radde, Nicole E. and Hammerschmidt, Matthias}, + journal = {Cell Death & Disease}, + title = {A sphingolipid rheostat controls apoptosis versus apical cell extrusion as alternative tumour-suppressive mechanisms}, + year = {2024}, + issn = {2041-4889}, + month = oct, + number = {10}, + volume = {15}, + creationdate = {2024-10-17T16:30:11}, + doi = {10.1038/s41419-024-07134-2}, + modificationdate = {2024-10-17T16:30:11}, + publisher = {Springer Science and Business Media LLC}, +} + +@Article{LakrisenkoPat2024, + author = {Lakrisenko, Polina and Pathirana, Dilan and Weindl, Daniel and Hasenauer, Jan}, + journal = {PLOS ONE}, + title = {Benchmarking methods for computing local sensitivities in ordinary differential equation models at dynamic and steady states}, + year = {2024}, + month = {10}, + number = {10}, + pages = {1-19}, + volume = {19}, + abstract = {Estimating parameters of dynamic models from experimental data is a challenging, and often computationally-demanding task. It requires a large number of model simulations and objective function gradient computations, if gradient-based optimization is used. In many cases, steady-state computation is a part of model simulation, either due to steady-state data or an assumption that the system is at steady state at the initial time point. Various methods are available for steady-state and gradient computation. Yet, the most efficient pair of methods (one for steady states, one for gradients) for a particular model is often not clear. In order to facilitate the selection of methods, we explore six method pairs for computing the steady state and sensitivities at steady state using six real-world problems. The method pairs involve numerical integration or Newton’s method to compute the steady-state, and—for both forward and adjoint sensitivity analysis—numerical integration or a tailored method to compute the sensitivities at steady-state. Our evaluation shows that all method pairs provide accurate steady-state and gradient values, and that the two method pairs that combine numerical integration for the steady-state with a tailored method for the sensitivities at steady-state were the most robust, and amongst the most computationally-efficient. We also observed that while Newton’s method for steady-state computation yields a substantial speedup compared to numerical integration, it may lead to a large number of simulation failures. Overall, our study provides a concise overview across current methods for computing sensitivities at steady state. While our study shows that there is no universally-best method pair, it also provides guidance to modelers in choosing the right methods for a problem at hand.}, + creationdate = {2024-11-08T18:16:35}, + doi = {10.1371/journal.pone.0312148}, + modificationdate = {2024-11-08T18:17:06}, + publisher = {Public Library of Science}, +} + @Comment{jabref-meta: databaseType:bibtex;} From bd98004511b8b93a20e77f453c849806b939174a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 27 Nov 2024 11:20:09 +0100 Subject: [PATCH 14/35] Fix `NLoptOptimizer.__repr__` (#1518) An incorrect value was used for `local_options`. --- pypesto/optimize/optimizer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypesto/optimize/optimizer.py b/pypesto/optimize/optimizer.py index b05570e73..c445a2776 100644 --- a/pypesto/optimize/optimizer.py +++ b/pypesto/optimize/optimizer.py @@ -1169,7 +1169,7 @@ def __repr__(self) -> str: if self.options is not None: rep += f" options={self.options}" if self.local_options is not None: - rep += f" local_options={self.local_methods}" + rep += f" local_options={self.local_options}" return rep + ">" @minimize_decorator_collection From 8f1bf803fd00a504be1a5979d3679059108d7c9b Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 27 Nov 2024 11:28:04 +0100 Subject: [PATCH 15/35] Improve exception-handling in SacessOptimizer (#1517) Previously, uncaught exceptions in SacessOptimizer worker processes would have resulted in deadlocks in `SacessOptimizer.minimize`. These changes will (in most cases) prevent deadlocks and other errors due to missing results from individual workers. Furthermore, this will lead to termination of minimize() relatively soon after an error occurred on some worker - not only after some other exit criterion is met. Closes #1512. --- pypesto/optimize/ess/ess.py | 2 + pypesto/optimize/ess/sacess.py | 158 +++++++++++++++++++++++++-------- test/optimize/test_optimize.py | 35 ++++++++ 3 files changed, 157 insertions(+), 38 deletions(-) diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index 2a86e18d9..74d88f641 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -37,6 +37,8 @@ class ESSExitFlag(int, enum.Enum): MAX_EVAL = -2 # Exited after exhausting wall-time budget MAX_TIME = -3 + # Termination because for other reason than exit criteria + ERROR = -99 class OptimizerFactory(Protocol): diff --git a/pypesto/optimize/ess/sacess.py b/pypesto/optimize/ess/sacess.py index 113310f25..9340b23da 100644 --- a/pypesto/optimize/ess/sacess.py +++ b/pypesto/optimize/ess/sacess.py @@ -7,6 +7,7 @@ import multiprocessing import os import time +from contextlib import suppress from dataclasses import dataclass from math import ceil, sqrt from multiprocessing import get_context @@ -20,6 +21,7 @@ import pypesto +from ... import MemoryHistory from ...startpoint import StartpointMethod from ...store.read_from_hdf5 import read_result from ...store.save_to_hdf5 import write_result @@ -331,12 +333,18 @@ def minimize( n_eval_total = sum( worker_result.n_eval for worker_result in self.worker_results ) - logger.info( - f"{self.__class__.__name__} stopped after {walltime:3g}s " - f"and {n_eval_total} objective evaluations " - f"with global best {result.optimize_result[0].fval}." - ) - + if len(result.optimize_result): + logger.info( + f"{self.__class__.__name__} stopped after {walltime:3g}s " + f"and {n_eval_total} objective evaluations " + f"with global best {result.optimize_result[0].fval}." + ) + else: + logger.error( + f"{self.__class__.__name__} stopped after {walltime:3g}s " + f"and {n_eval_total} objective evaluations without producing " + "a result." + ) return result def _create_result(self, problem: Problem) -> pypesto.Result: @@ -345,25 +353,40 @@ def _create_result(self, problem: Problem) -> pypesto.Result: Creates an overall Result object from the results saved by the workers. """ # gather results from workers and delete temporary result files - result = None + result = pypesto.Result() + retry_after_sleep = True for worker_idx in range(self.num_workers): tmp_result_filename = SacessWorker.get_temp_result_filename( worker_idx, self._tmpdir ) + tmp_result = None try: tmp_result = read_result( tmp_result_filename, problem=False, optimize=True ) except FileNotFoundError: # wait and retry, maybe the file wasn't found due to some filesystem latency issues - time.sleep(5) - tmp_result = read_result( - tmp_result_filename, problem=False, optimize=True - ) + if retry_after_sleep: + time.sleep(5) + # waiting once is enough - don't wait for every worker + retry_after_sleep = False + + try: + tmp_result = read_result( + tmp_result_filename, problem=False, optimize=True + ) + except FileNotFoundError: + logger.error( + f"Worker {worker_idx} did not produce a result." + ) + continue + else: + logger.error( + f"Worker {worker_idx} did not produce a result." + ) + continue - if result is None: - result = tmp_result - else: + if tmp_result: result.optimize_result.append( tmp_result.optimize_result, sort=False, @@ -375,7 +398,8 @@ def _create_result(self, problem: Problem) -> pypesto.Result: filename = SacessWorker.get_temp_result_filename( worker_idx, self._tmpdir ) - os.remove(filename) + with suppress(FileNotFoundError): + os.remove(filename) # delete tmpdir if empty try: self._tmpdir.rmdir() @@ -397,6 +421,7 @@ class SacessManager: Attributes ---------- + _dim: Dimension of the optimization problem _num_workers: Number of workers _ess_options: ESS options for each worker _best_known_fx: Best objective value encountered so far @@ -410,6 +435,7 @@ class SacessManager: _rejection_threshold: Threshold for relative objective improvements that incoming solutions have to pass to be accepted _lock: Lock for accessing shared state. + _terminate: Flag to signal termination of the SACESS run to workers _logger: A logger instance _options: Further optimizer hyperparameters. """ @@ -421,6 +447,7 @@ def __init__( dim: int, options: SacessOptions = None, ): + self._dim = dim self._options = options or SacessOptions() self._num_workers = len(ess_options) self._ess_options = [shmem_manager.dict(o) for o in ess_options] @@ -440,6 +467,7 @@ def __init__( self._worker_scores = shmem_manager.Array( "d", range(self._num_workers) ) + self._terminate = shmem_manager.Value("b", False) self._worker_comms = shmem_manager.Array("i", [0] * self._num_workers) self._lock = shmem_manager.RLock() self._logger = logging.getLogger() @@ -550,6 +578,16 @@ def submit_solution( ) self._rejections.value = 0 + def abort(self): + """Abort the SACESS run.""" + with self._lock: + self._terminate.value = True + + def aborted(self) -> bool: + """Whether this run was aborted.""" + with self._lock: + return self._terminate.value + class SacessWorker: """A SACESS worker. @@ -641,7 +679,7 @@ def run( ess = self._setup_ess(startpoint_method) # run ESS until exit criteria are met, but start at least one iteration - while self._keep_going() or ess.n_iter == 0: + while self._keep_going(ess) or ess.n_iter == 0: # perform one ESS iteration ess._do_iteration() @@ -667,19 +705,42 @@ def run( f"(best: {self._best_known_fx}, " f"n_eval: {ess.evaluator.n_eval})." ) - - ess.history.finalize(exitflag=ess.exit_flag.name) - worker_result = SacessWorkerResult( - x=ess.x_best, - fx=ess.fx_best, - history=ess.history, - n_eval=ess.evaluator.n_eval, - n_iter=ess.n_iter, - exit_flag=ess.exit_flag, - ) + self._finalize(ess) + + def _finalize(self, ess: ESSOptimizer = None): + """Finalize the worker.""" + # Whatever happens here, we need to put something to the queue before + # returning to avoid deadlocks. + worker_result = None + if ess is not None: + try: + ess.history.finalize(exitflag=ess.exit_flag.name) + ess._report_final() + worker_result = SacessWorkerResult( + x=ess.x_best, + fx=ess.fx_best, + history=ess.history, + n_eval=ess.evaluator.n_eval, + n_iter=ess.n_iter, + exit_flag=ess.exit_flag, + ) + except Exception as e: + self._logger.exception( + f"Worker {self._worker_idx} failed to finalize: {e}" + ) + if worker_result is None: + # Create some dummy result + worker_result = SacessWorkerResult( + x=np.full(self._manager._dim, np.nan), + fx=np.nan, + history=MemoryHistory(), + n_eval=0, + n_iter=0, + exit_flag=ESSExitFlag.ERROR, + ) self._manager._result_queue.put(worker_result) + self._logger.debug(f"Final configuration: {self._ess_kwargs}") - ess._report_final() def _setup_ess(self, startpoint_method: StartpointMethod) -> ESSOptimizer: """Run ESS.""" @@ -821,7 +882,7 @@ def replace_solution(refset: RefSet, x: np.ndarray, fx: float): fx=fx, ) - def _keep_going(self): + def _keep_going(self, ess: ESSOptimizer) -> bool: """Check exit criteria. Returns @@ -830,14 +891,26 @@ def _keep_going(self): """ # elapsed time if time.time() - self._start_time >= self._max_walltime_s: - self.exit_flag = ESSExitFlag.MAX_TIME + ess.exit_flag = ESSExitFlag.MAX_TIME self._logger.debug( f"Max walltime ({self._max_walltime_s}s) exceeded." ) return False - + # other reasons for termination (some worker failed, ...) + if self._manager.aborted(): + ess.exit_flag = ESSExitFlag.ERROR + self._logger.debug("Manager requested termination.") + return False return True + def abort(self): + """Send signal to abort.""" + self._logger.error(f"Worker {self._worker_idx} aborting.") + # signal to manager + self._manager.abort() + + self._finalize(None) + @staticmethod def get_temp_result_filename(worker_idx: int, tmpdir: str | Path) -> str: return str(Path(tmpdir, f"sacess-{worker_idx:02d}_tmp.h5").absolute()) @@ -853,15 +926,24 @@ def _run_worker( Helper function as entrypoint for sacess worker processes. """ - # different random seeds per process - np.random.seed((os.getpid() * int(time.time() * 1000)) % 2**32) - - # Forward log messages to the logging process - h = logging.handlers.QueueHandler(log_process_queue) - worker._logger = logging.getLogger(multiprocessing.current_process().name) - worker._logger.addHandler(h) + try: + # different random seeds per process + np.random.seed((os.getpid() * int(time.time() * 1000)) % 2**32) + + # Forward log messages to the logging process + h = logging.handlers.QueueHandler(log_process_queue) + worker._logger = logging.getLogger( + multiprocessing.current_process().name + ) + worker._logger.addHandler(h) - return worker.run(problem=problem, startpoint_method=startpoint_method) + return worker.run(problem=problem, startpoint_method=startpoint_method) + except Exception as e: + with suppress(Exception): + worker._logger.exception( + f"Worker {worker._worker_idx} failed: {e}" + ) + worker.abort() def get_default_ess_options( diff --git a/test/optimize/test_optimize.py b/test/optimize/test_optimize.py index 48ebdea55..3e6a5b4ed 100644 --- a/test/optimize/test_optimize.py +++ b/test/optimize/test_optimize.py @@ -18,6 +18,7 @@ import pypesto import pypesto.optimize as optimize +from pypesto import Objective from pypesto.optimize.ess import ( ESSOptimizer, FunctionEvaluatorMP, @@ -577,6 +578,40 @@ def test_ess_refset_repr(): ) +class FunctionOrError: + """Callable that raises an error every nth invocation.""" + + def __init__(self, fun, error_frequency=100): + self.counter = 0 + self.error_frequency = error_frequency + self.fun = fun + + def __call__(self, *args, **kwargs): + self.counter += 1 + if self.counter % self.error_frequency == 0: + raise RuntimeError("Intentional error.") + return self.fun(*args, **kwargs) + + +def test_sacess_worker_error(capsys): + """Check that SacessOptimizer does not hang if an error occurs on a worker.""" + objective = Objective( + fun=FunctionOrError(sp.optimize.rosen), grad=sp.optimize.rosen_der + ) + problem = pypesto.Problem( + objective=objective, lb=0 * np.ones((1, 2)), ub=1 * np.ones((1, 2)) + ) + sacess = SacessOptimizer( + num_workers=2, + max_walltime_s=2, + sacess_loglevel=logging.DEBUG, + ess_loglevel=logging.DEBUG, + ) + res = sacess.minimize(problem) + assert isinstance(res, pypesto.Result) + assert "Intentional error." in capsys.readouterr().err + + def test_scipy_integrated_grad(): integrated = True obj = rosen_for_sensi(max_sensi_order=2, integrated=integrated)["obj"] From d39b959fef3a08746b9718ef4de1a7a37a471b09 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 27 Nov 2024 11:28:25 +0100 Subject: [PATCH 16/35] CI: don't cache amici; pytest ignore amici_models; passenv BNGPATH (#1523) * Don't cache the amici installation. We are using amici/develop on purpose to detect problems early on. E.g., #1520 should have been detected much earlier. * Ignore `amici_models` when collecting tests * Use external BNGPATH if set --- .github/workflows/ci.yml | 2 +- pytest.ini | 1 + tox.ini | 10 +++++++--- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4db56ab55..f55f8f6b5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -157,7 +157,7 @@ jobs: - name: Run tests timeout-minutes: 35 - run: tox -e petab + run: tox -e petab && tox e -e petab -- pip uninstall -y amici env: CC: clang CXX: clang++ diff --git a/pytest.ini b/pytest.ini index c852c729e..82ae58578 100644 --- a/pytest.ini +++ b/pytest.ini @@ -2,3 +2,4 @@ addopts = "--doctest-modules" filterwarnings = ignore:.*inspect.getargspec\(\) is deprecated.*:DeprecationWarning +norecursedirs = amici_models diff --git a/tox.ini b/tox.ini index e9193e0ad..d1937e4ee 100644 --- a/tox.ini +++ b/tox.ini @@ -29,7 +29,7 @@ envlist = # Base-environment [testenv] -passenv = AMICI_PARALLEL_COMPILE,CC,CXX,MPLBACKEND +passenv = AMICI_PARALLEL_COMPILE,CC,CXX,MPLBACKEND,BNGPATH # Sub-environments # inherit settings defined in the base @@ -75,10 +75,14 @@ description = Test basic functionality on Windows [testenv:petab] -extras = test,amici,petab,pyswarm,roadrunner +extras = test,petab,pyswarm,roadrunner deps = git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master\#subdirectory=src/python - git+https://github.com/AMICI-dev/amici.git@develop\#egg=amici&subdirectory=python/sdist +# always install amici from develop branch, avoid caching +# to skip re-installation, run `tox -e petab --override testenv:petab.commands_pre=` +commands_pre = + python3 -m pip uninstall -y amici + python3 -m pip install git+https://github.com/AMICI-dev/amici.git@develop\#egg=amici&subdirectory=python/sdist commands = python3 -m pip install git+https://github.com/PEtab-dev/petab_test_suite@main python3 -m pip install git+https://github.com/pysb/pysb@master From 1eb61d798686c4f3bb907dc20d3f7abbfc9279b7 Mon Sep 17 00:00:00 2001 From: Stephan Grein Date: Wed, 27 Nov 2024 11:49:45 +0100 Subject: [PATCH 17/35] Update CODEOWNERS (#1516) Co-authored-by: Vincent Wieland --- .github/CODEOWNERS | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 7d10fd310..fe396968f 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -42,6 +42,7 @@ /pypesto/select/ @dilpath /pypesto/startpoint/ @PaulJonasJost /pypesto/store/ @PaulJonasJost +/pypesto/visualize/ @stephanmg # Tests /test/base/ @PaulJonasJost @vwiela From 2ed5f8f621b3eba275d7ee4678cbf7b383ee8c22 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 27 Nov 2024 16:09:39 +0100 Subject: [PATCH 18/35] Don't modify sys.path for amici model imports (#1522) Not necessary when using `amici.import_model_module`. --- pypesto/petab/objective_creator.py | 5 ----- test/petab/test_amici_predictor.py | 6 ------ test/petab/test_sbml_conversion.py | 3 --- 3 files changed, 14 deletions(-) diff --git a/pypesto/petab/objective_creator.py b/pypesto/petab/objective_creator.py index 72f98cf03..697c0f9a3 100644 --- a/pypesto/petab/objective_creator.py +++ b/pypesto/petab/objective_creator.py @@ -7,7 +7,6 @@ import os import re import shutil -import sys import warnings from abc import ABC, abstractmethod from collections.abc import Iterable, Sequence @@ -145,10 +144,6 @@ def create_model( f"compilation: Not a folder." ) - # add module to path - if self.output_folder not in sys.path: - sys.path.insert(0, self.output_folder) - # compile if self._must_compile(force_compile): logger.info( diff --git a/test/petab/test_amici_predictor.py b/test/petab/test_amici_predictor.py index 2d23620a4..d99fa7d52 100644 --- a/test/petab/test_amici_predictor.py +++ b/test/petab/test_amici_predictor.py @@ -2,7 +2,6 @@ import os import shutil -import sys import amici import libsbml @@ -34,13 +33,10 @@ def conversion_reaction_model(): # try to import the exisiting model, if possible try: - sys.path.insert(0, os.path.abspath(model_output_dir)) model_module = amici.import_model_module(model_name, model_output_dir) model = model_module.getModel() except ValueError: # read in and adapt the sbml slightly - if os.path.abspath(model_output_dir) in sys.path: - sys.path.remove(os.path.abspath(model_output_dir)) sbml_importer = amici.SbmlImporter(sbml_file) # add observables to sbml model @@ -95,7 +91,6 @@ def create_intial_assignment(sbml_model, spec_id): ) # Importing the module and loading the model - sys.path.insert(0, os.path.abspath(model_output_dir)) model_module = amici.import_model_module(model_name, model_output_dir) model = model_module.getModel() except RuntimeError as err: @@ -107,7 +102,6 @@ def create_intial_assignment(sbml_model, spec_id): "Delete the conversion_reaction_enhanced model from your python " "path and retry. Your python path is currently:" ) - print(sys.path) print("Original error message:") raise err diff --git a/test/petab/test_sbml_conversion.py b/test/petab/test_sbml_conversion.py index 22653a0a6..7284fb8a7 100644 --- a/test/petab/test_sbml_conversion.py +++ b/test/petab/test_sbml_conversion.py @@ -1,6 +1,5 @@ import os import re -import sys import unittest import warnings @@ -11,8 +10,6 @@ from ..util import load_amici_objective -sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) - optimizers = { "scipy": [ "Nelder-Mead", From 861318a1fb11c5e6bf3345bd2413dc96510c4c71 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Nov 2024 12:38:54 +0100 Subject: [PATCH 19/35] Previously, only writing `Result` was supported. (#1528) Changes: * `OptimizeResult`, `OptimizerResult`(s) can be written by `OptimizationResultHDF5Writer.write()` (Closes #1526) * `OptimizerResult`s can be written incrementally by `OptimizationResultHDF5Writer.write_optimizer_result()` (Closes #1527) * Fixed an `AttributeError` in `pypesto.store.save_to_hdf5.check_overwrite` with `h5py.Group`s --- pypesto/store/save_to_hdf5.py | 105 ++++++++++++++++++++++++---------- test/base/test_store.py | 27 ++++++++- 2 files changed, 99 insertions(+), 33 deletions(-) diff --git a/pypesto/store/save_to_hdf5.py b/pypesto/store/save_to_hdf5.py index 3804bfbbb..a4ac4e703 100644 --- a/pypesto/store/save_to_hdf5.py +++ b/pypesto/store/save_to_hdf5.py @@ -1,23 +1,22 @@ """Include functions for saving various results to hdf5.""" +from __future__ import annotations import logging import os from numbers import Integral from pathlib import Path -from typing import Union import h5py import numpy as np +from .. import OptimizeResult, OptimizerResult from ..result import ProfilerResult, Result, SampleResult from .hdf5 import write_array, write_float_array logger = logging.getLogger(__name__) -def check_overwrite( - f: Union[h5py.File, h5py.Group], overwrite: bool, target: str -): +def check_overwrite(f: h5py.File | h5py.Group, overwrite: bool, target: str): """ Check whether target already exists. @@ -36,7 +35,7 @@ def check_overwrite( del f[target] else: raise RuntimeError( - f"File `{f.filename}` already exists and contains " + f"File `{f.file.filename}` already exists and contains " f"information about {target} result. " f"If you wish to overwrite the file, set " f"`overwrite=True`." @@ -53,7 +52,7 @@ class ProblemHDF5Writer: HDF5 result file name """ - def __init__(self, storage_filename: Union[str, Path]): + def __init__(self, storage_filename: str | Path): """ Initialize writer. @@ -106,7 +105,7 @@ class OptimizationResultHDF5Writer: HDF5 result file name """ - def __init__(self, storage_filename: Union[str, Path]): + def __init__(self, storage_filename: str | Path): """ Initialize Writer. @@ -117,32 +116,76 @@ def __init__(self, storage_filename: Union[str, Path]): """ self.storage_filename = str(storage_filename) - def write(self, result: Result, overwrite=False): - """Write HDF5 result file from pyPESTO result object.""" - # Create destination directory - if isinstance(self.storage_filename, str): - basedir = os.path.dirname(self.storage_filename) - if basedir: - os.makedirs(basedir, exist_ok=True) + def write( + self, + result: Result + | OptimizeResult + | OptimizerResult + | list[OptimizerResult], + overwrite=False, + ): + """Write HDF5 result file from pyPESTO result object. + + Parameters + ---------- + result: Result to be saved. + overwrite: Boolean, whether already existing results should be + overwritten. This applies to the whole list of results, not only to + individual results. See :meth:`write_optimizer_result` for + incrementally writing a sequence of `OptimizerResult`. + """ + Path(self.storage_filename).parent.mkdir(parents=True, exist_ok=True) + + if isinstance(result, Result): + results = result.optimize_result.list + elif isinstance(result, OptimizeResult): + results = result.list + elif isinstance(result, list): + results = result + elif isinstance(result, OptimizerResult): + results = [result] + else: + raise ValueError(f"Unsupported type for `result`: {type(result)}.") with h5py.File(self.storage_filename, "a") as f: check_overwrite(f, overwrite, "optimization") optimization_grp = f.require_group("optimization") - # settings = - # optimization_grp.create_dataset("settings", settings, dtype=) results_grp = optimization_grp.require_group("results") - for start in result.optimize_result.list: - start_id = start["id"] - start_grp = results_grp.require_group(start_id) - for key in start.keys(): - if key == "history": - continue - if isinstance(start[key], np.ndarray): - write_array(start_grp, key, start[key]) - elif start[key] is not None: - start_grp.attrs[key] = start[key] - f.flush() + for start in results: + self._do_write_optimizer_result(start, results_grp, overwrite) + + def write_optimizer_result( + self, result: OptimizerResult, overwrite: bool = False + ): + """Write HDF5 result file from pyPESTO result object. + + Parameters + ---------- + result: Result to be saved. + overwrite: Boolean, whether already existing results with the same ID + should be overwritten.s + """ + Path(self.storage_filename).parent.mkdir(parents=True, exist_ok=True) + + with h5py.File(self.storage_filename, "a") as f: + results_grp = f.require_group("optimization/results") + self._do_write_optimizer_result(result, results_grp, overwrite) + + def _do_write_optimizer_result( + self, result: OptimizerResult, g: h5py.Group = None, overwrite=False + ): + """Write an OptimizerResult to the given group.""" + sub_group_id = result["id"] + check_overwrite(g, overwrite, sub_group_id) + start_grp = g.require_group(sub_group_id) + for key in result.keys(): + if key == "history": + continue + if isinstance(result[key], np.ndarray): + write_array(start_grp, key, result[key]) + elif result[key] is not None: + start_grp.attrs[key] = result[key] class SamplingResultHDF5Writer: @@ -155,7 +198,7 @@ class SamplingResultHDF5Writer: HDF5 result file name """ - def __init__(self, storage_filename: Union[str, Path]): + def __init__(self, storage_filename: str | Path): """ Initialize Writer. @@ -208,7 +251,7 @@ class ProfileResultHDF5Writer: HDF5 result file name """ - def __init__(self, storage_filename: Union[str, Path]): + def __init__(self, storage_filename: str | Path): """ Initialize Writer. @@ -241,7 +284,7 @@ def write(self, result: Result, overwrite: bool = False): @staticmethod def _write_profiler_result( - parameter_profile: Union[ProfilerResult, None], result_grp: h5py.Group + parameter_profile: ProfilerResult | None, result_grp: h5py.Group ) -> None: """Write a single ProfilerResult to hdf5. @@ -267,7 +310,7 @@ def _write_profiler_result( def write_result( result: Result, - filename: Union[str, Path], + filename: str | Path, overwrite: bool = False, problem: bool = True, optimize: bool = False, diff --git a/test/base/test_store.py b/test/base/test_store.py index d90a8030d..840440c70 100644 --- a/test/base/test_store.py +++ b/test/base/test_store.py @@ -1,9 +1,11 @@ """Test the `pypesto.store` module.""" import os -import tempfile +from pathlib import Path +from tempfile import TemporaryDirectory import numpy as np +import pytest import scipy.optimize as so import pypesto @@ -52,7 +54,7 @@ def test_storage_opt_result(): minimize_result = create_optimization_result() - with tempfile.TemporaryDirectory(dir=".") as tmpdirname: + with TemporaryDirectory(dir=".") as tmpdirname: result_file_name = os.path.join(tmpdirname, "a", "b", "result.h5") opt_result_writer = OptimizationResultHDF5Writer(result_file_name) opt_result_writer.write(minimize_result) @@ -89,6 +91,27 @@ def test_storage_opt_result_update(hdf5_file): assert opt_res[key] == read_result.optimize_result[i][key] +def test_write_optimizer_results_incrementally(): + """Test writing optimizer results incrementally to the same file.""" + res = create_optimization_result() + res1, res2 = res.optimize_result.list[:2] + + with TemporaryDirectory() as tmp_dir: + result_path = Path(tmp_dir, "result.h5") + writer = OptimizationResultHDF5Writer(result_path) + writer.write_optimizer_result(res1) + writer.write_optimizer_result(res2) + reader = OptimizationResultHDF5Reader(result_path) + read_res = reader.read() + assert len(read_res.optimize_result) == 2 + + # overwriting works + writer.write_optimizer_result(res1, overwrite=True) + # overwriting attempt fails without overwrite=True + with pytest.raises(RuntimeError): + writer.write_optimizer_result(res1) + + def test_storage_problem(hdf5_file): problem = create_problem() problem_writer = ProblemHDF5Writer(hdf5_file) From d1c9e527aeb106618519b39e003be34a3c1bff35 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Nov 2024 12:39:15 +0100 Subject: [PATCH 20/35] Set `OptimizerResult.optimizer` in `Optimizer.minimize` (#1525) Fixes #1524. --- pypesto/optimize/ess/ess.py | 1 - pypesto/optimize/optimizer.py | 37 +++++++++++++++++++++++----------- pypesto/optimize/task.py | 1 - test/optimize/test_optimize.py | 1 + 4 files changed, 26 insertions(+), 14 deletions(-) diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index 74d88f641..cdcab2f12 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -403,7 +403,6 @@ def _create_result(self) -> pypesto.Result: for i, optimizer_result in enumerate(self.local_solutions): i_result += 1 optimizer_result.id = f"Local solution {i}" - optimizer_result.optimizer = str(self.local_optimizer) result.optimize_result.append(optimizer_result) if self._result_includes_refset: diff --git a/pypesto/optimize/optimizer.py b/pypesto/optimize/optimizer.py index c445a2776..659303aed 100644 --- a/pypesto/optimize/optimizer.py +++ b/pypesto/optimize/optimizer.py @@ -39,7 +39,7 @@ def __init__(self, optimizer: str): def hierarchical_decorator(minimize): """Add inner parameters to the optimizer result. - Default decorator for the minimize() method. + Default decorator for the :meth:`Optimizer.minimize` method. """ @wraps(minimize) @@ -81,7 +81,7 @@ def wrapped_minimize( def history_decorator(minimize): """Initialize and extract information stored in the history. - Default decorator for the minimize() method. + Default decorator for the :meth:`Optimizer.minimize` method. """ @wraps(minimize) @@ -140,7 +140,11 @@ def wrapped_minimize( logger.error(f"start {id} failed:\n{trace}") result = OptimizerResult( - x0=x0, exitflag=-1, message=str(err), id=id + x0=x0, + exitflag=-1, + message=str(err), + id=id, + optimizer=str(self), ) else: raise @@ -163,7 +167,7 @@ def wrapped_minimize( def time_decorator(minimize): """Measure time of optimization. - Default decorator for the minimize() method to take time. + Default decorator for the :meth:`Optimizer.minimize` method to take time. Currently, the method time.time() is used, which measures the wall-clock time. """ @@ -196,8 +200,8 @@ def wrapped_minimize( def fix_decorator(minimize): """Include also fixed parameters in the result arrays of minimize(). - Default decorator for the minimize() method (nans will be inserted in the - derivatives). + Default decorator for the :meth:`Optimizer.minimize` method (nans will be + inserted in the derivatives). """ @wraps(minimize) @@ -523,6 +527,7 @@ def fun(x): hess=getattr(res, "hess", None), exitflag=res.status, message=res.message, + optimizer=str(self), ) return optimizer_result @@ -612,7 +617,10 @@ def minimize( # the ipopt return object is a scipy.optimize.OptimizeResult return OptimizerResult( - x=ret.x, exitflag=ret.status, message=ret.message + x=ret.x, + exitflag=ret.status, + message=ret.message, + optimizer=str(self), ) def is_least_squares(self): @@ -630,7 +638,7 @@ def __init__(self, options: dict = None): if self.options is None: self.options = DlibOptimizer.get_default_options(self) elif "maxiter" not in self.options: - raise KeyError("Dlib options are missing the key word " "maxiter.") + raise KeyError("Dlib options are missing the keyword maxiter.") def __repr__(self) -> str: rep = f"<{self.__class__.__name__}" @@ -677,7 +685,7 @@ def get_fval_vararg(*x): 0.002, ) - optimizer_result = OptimizerResult() + optimizer_result = OptimizerResult(optimizer=str(self)) return optimizer_result @@ -737,7 +745,9 @@ def minimize( problem.objective.get_fval, lb, ub, **self.options ) - optimizer_result = OptimizerResult(x=np.array(xopt), fval=fopt) + optimizer_result = OptimizerResult( + x=np.array(xopt), fval=fopt, optimizer=str(self) + ) return optimizer_result @@ -821,7 +831,7 @@ def minimize( ) optimizer_result = OptimizerResult( - x=np.array(result[0]), fval=result[1] + x=np.array(result[0]), fval=result[1], optimizer=str(self) ) return optimizer_result @@ -901,7 +911,7 @@ def minimize( ) optimizer_result = OptimizerResult( - x=np.array(result.x), fval=result.fun + x=np.array(result.x), fval=result.fun, optimizer=str(self) ) return optimizer_result @@ -1019,6 +1029,7 @@ def successively_working_fval(swarm: np.ndarray) -> np.ndarray: optimizer_result = OptimizerResult( x=pos, fval=float(cost), + optimizer=str(self), ) return optimizer_result @@ -1249,6 +1260,7 @@ def nlopt_objective(x, grad): fval=opt.last_optimum_value(), message=msg, exitflag=opt.last_optimize_result(), + optimizer=str(self), ) return optimizer_result @@ -1433,6 +1445,7 @@ def minimize( hess=opt.hess, message=msg, exitflag=opt.exitflag, + optimizer=str(self), ) return optimizer_result diff --git a/pypesto/optimize/task.py b/pypesto/optimize/task.py index 10ae83dd8..7482097e1 100644 --- a/pypesto/optimize/task.py +++ b/pypesto/optimize/task.py @@ -63,7 +63,6 @@ def execute(self) -> OptimizerResult: history_options=self.history_options, optimize_options=self.optimize_options, ) - optimizer_result.optimizer = str(self.optimizer) if not self.optimize_options.report_hess: optimizer_result.hess = None diff --git a/test/optimize/test_optimize.py b/test/optimize/test_optimize.py index 3e6a5b4ed..6f8260b90 100644 --- a/test/optimize/test_optimize.py +++ b/test/optimize/test_optimize.py @@ -309,6 +309,7 @@ def check_minimize(problem, library, solver, allow_failed_starts=False): ]: assert np.isfinite(result.optimize_result.list[0]["fval"]) assert result.optimize_result.list[0]["x"] is not None + assert result.optimize_result.list[0]["optimizer"] is not None def test_trim_results(problem): From 7448b92e243c6ea0ddeb56f7a16174e93e71402f Mon Sep 17 00:00:00 2001 From: Stephan Grein Date: Thu, 28 Nov 2024 18:17:36 +0100 Subject: [PATCH 21/35] Push latest pyPESTO (develop) build as image with tag `latest` to docker.io (#1083) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Adding proposal to deploy to docker hub. * Fix path to Dockerfile. * Update publish_dockerhub.yml * Update publish_dockerhub.yml * Update publish_dockerhub.yml --------- Co-authored-by: Yannik Schälte <31767307+yannikschaelte@users.noreply.github.com> Co-authored-by: Paul Jonas Jost <70631928+PaulJonasJost@users.noreply.github.com> --- .github/workflows/publish_dockerhub.yml | 32 +++++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 .github/workflows/publish_dockerhub.yml diff --git a/.github/workflows/publish_dockerhub.yml b/.github/workflows/publish_dockerhub.yml new file mode 100644 index 000000000..55d970055 --- /dev/null +++ b/.github/workflows/publish_dockerhub.yml @@ -0,0 +1,32 @@ +name: Build and Push Docker Image + +on: + push: + branches: + - main + paths: + - 'docker/**' + - '.github/workflows/docker-publish.yml' + workflow_dispatch: + +jobs: + build-and-push: + runs-on: ubuntu-latest + + steps: + - name: Check out the repository + uses: actions/checkout@v2 + + - name: Log in to Docker Hub + uses: docker/login-action@v2 + with: + username: ${{ secrets.DOCKER_USERNAME }} + password: ${{ secrets.DOCKER_TOKEN }} + + - name: Build and tag the Docker image + run: | + docker build -t ICB_DCM/pypesto:latest -f docker/Dockerfile . + + - name: Push the Docker image to Docker Hub + run: | + docker push ICB_DCM/pypesto:latest From 010e93a877bab40eaaaffecc9f14af714c455c77 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 2 Dec 2024 11:50:31 +0100 Subject: [PATCH 22/35] Bump actions/checkout from 2 to 4 (#1532) Bumps [actions/checkout](https://github.com/actions/checkout) from 2 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v2...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/publish_dockerhub.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/publish_dockerhub.yml b/.github/workflows/publish_dockerhub.yml index 55d970055..4b05ad1fc 100644 --- a/.github/workflows/publish_dockerhub.yml +++ b/.github/workflows/publish_dockerhub.yml @@ -15,7 +15,7 @@ jobs: steps: - name: Check out the repository - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: Log in to Docker Hub uses: docker/login-action@v2 From 9a20cd880c8841a6f876aa52b2fe503321aa736b Mon Sep 17 00:00:00 2001 From: Doresic <85789271+Doresic@users.noreply.github.com> Date: Mon, 2 Dec 2024 14:57:59 +0100 Subject: [PATCH 23/35] Visualize: fix flatten of observable mapping with one observable (#1515) * Fix flatten of axes if only one obs * Same for other, implement tests of both * Decrease maxiter of visualize tests --- pypesto/visualize/observable_mapping.py | 8 ++-- pypesto/visualize/ordinal_categories.py | 26 ++++------- test/visualize/test_visualize.py | 57 +++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 22 deletions(-) diff --git a/pypesto/visualize/observable_mapping.py b/pypesto/visualize/observable_mapping.py index ec378fd1c..250a82dd5 100644 --- a/pypesto/visualize/observable_mapping.py +++ b/pypesto/visualize/observable_mapping.py @@ -129,7 +129,7 @@ def visualize_estimated_observable_mapping( n_axes = n_relative_observables + n_semiquant_observables n_rows = int(np.ceil(np.sqrt(n_axes))) n_cols = int(np.ceil(n_axes / n_rows)) - _, axes = plt.subplots(n_rows, n_cols, **kwargs) + _, axes = plt.subplots(n_rows, n_cols, squeeze=False, **kwargs) axes = axes.flatten() # Plot the estimated observable mapping for relative observables. @@ -246,8 +246,7 @@ def plot_linear_observable_mappings_from_pypesto_result( n_cols = int(np.ceil(n_relative_observables / n_rows)) # Make as many subplots as there are relative observables - _, axes = plt.subplots(n_rows, n_cols, **kwargs) - + _, axes = plt.subplots(n_rows, n_cols, squeeze=False, **kwargs) # Flatten the axes array axes = axes.flatten() @@ -590,8 +589,7 @@ def plot_splines_from_inner_result( n_cols = int(np.ceil(n_groups / n_rows)) # Make as many subplots as there are groups - _, axes = plt.subplots(n_rows, n_cols, **kwargs) - + _, axes = plt.subplots(n_rows, n_cols, squeeze=False, **kwargs) # Flatten the axes array axes = axes.flatten() diff --git a/pypesto/visualize/ordinal_categories.py b/pypesto/visualize/ordinal_categories.py index b93e97f22..c8e81def5 100644 --- a/pypesto/visualize/ordinal_categories.py +++ b/pypesto/visualize/ordinal_categories.py @@ -612,26 +612,18 @@ def _get_data_for_plotting( def _get_default_axes(n_groups, **kwargs): """Return a list of axes with the default layout.""" - # If there is only one group, make a figure with only one plot - if n_groups == 1: - # Make figure with only one plot - fig, ax = plt.subplots(1, 1, **kwargs) + # Choose number of rows and columns to be used for the subplots + n_rows = int(np.ceil(np.sqrt(n_groups))) + n_cols = int(np.ceil(n_groups / n_rows)) - axes = [ax] - # If there are multiple groups, make a figure with multiple plots - else: - # Choose number of rows and columns to be used for the subplots - n_rows = int(np.ceil(np.sqrt(n_groups))) - n_cols = int(np.ceil(n_groups / n_rows)) - - # Make as many subplots as there are groups - fig, axes = plt.subplots(n_rows, n_cols, **kwargs) + # Make as many subplots as there are groups + fig, axes = plt.subplots(n_rows, n_cols, squeeze=False, **kwargs) - # Increase the spacing between the subplots - fig.subplots_adjust(hspace=0.35, wspace=0.25) + # Increase the spacing between the subplots + fig.subplots_adjust(hspace=0.35, wspace=0.25) - # Flatten the axes array - axes = axes.flatten() + # Flatten the axes array + axes = axes.flatten() return axes diff --git a/test/visualize/test_visualize.py b/test/visualize/test_visualize.py index 149c65719..ecdc38eab 100644 --- a/test/visualize/test_visualize.py +++ b/test/visualize/test_visualize.py @@ -3,6 +3,7 @@ import os from collections.abc import Sequence from functools import wraps +from pathlib import Path import matplotlib.pyplot as plt import numpy as np @@ -1240,3 +1241,59 @@ def test_parameters_correlation_matrix(result_creation): result = result_creation() visualize.parameters_correlation_matrix(result) + + +@close_fig +def test_plot_ordinal_categories(): + example_ordinal_yaml = ( + Path(__file__).parent + / ".." + / ".." + / "doc" + / "example" + / "example_ordinal" + / "example_ordinal.yaml" + ) + petab_problem = petab.Problem.from_yaml(example_ordinal_yaml) + # Set seed for reproducibility. + np.random.seed(0) + optimizer = pypesto.optimize.ScipyOptimizer( + method="L-BFGS-B", options={"maxiter": 1} + ) + importer = pypesto.petab.PetabImporter(petab_problem, hierarchical=True) + problem = importer.create_problem() + result = pypesto.optimize.minimize( + problem=problem, n_starts=1, optimizer=optimizer + ) + visualize.plot_categories_from_pypesto_result(result) + + +@close_fig +def test_visualize_estimated_observable_mapping(): + example_semiquantitative_yaml = ( + Path(__file__).parent + / ".." + / ".." + / "doc" + / "example" + / "example_semiquantitative" + / "example_semiquantitative_linear.yaml" + ) + petab_problem = petab.Problem.from_yaml(example_semiquantitative_yaml) + # Set seed for reproducibility. + np.random.seed(0) + optimizer = pypesto.optimize.ScipyOptimizer( + method="L-BFGS-B", + options={ + "disp": None, + "ftol": 2.220446049250313e-09, + "gtol": 1e-5, + "maxiter": 1, + }, + ) + importer = pypesto.petab.PetabImporter(petab_problem, hierarchical=True) + problem = importer.create_problem() + result = pypesto.optimize.minimize( + problem=problem, n_starts=1, optimizer=optimizer + ) + visualize.visualize_estimated_observable_mapping(result, problem) From 25c2f841c4ddbf4143a2123a471efc112785c6d3 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 9 Dec 2024 10:30:45 +0100 Subject: [PATCH 24/35] Pin petab-select (#1536) Tests fail because of incompatible un-released changes in petab-select. --- setup.cfg | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/setup.cfg b/setup.cfg index d9eb70919..9019fe480 100644 --- a/setup.cfg +++ b/setup.cfg @@ -167,9 +167,9 @@ select = # Remove when vis is moved to PEtab Select networkx >= 2.5.1 # End remove - #petab-select >= 0.1.12 - # FIXME before merge - petab-select @ git+https://github.com/PEtab-dev/petab_select.git@develop + # TODO: pinned until https://github.com/ICB-DCM/pyPESTO/pull/1530 is done + # petab-select @ git+https://github.com/PEtab-dev/petab_select.git@develop + petab-select == 0.2.1 test = pytest >= 5.4.3 pytest-cov >= 2.10.0 From 2e5b8de7f1153a3381d45a573c64e7a0d07d91b7 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 10 Dec 2024 10:16:00 +0100 Subject: [PATCH 25/35] SacessOptimizer: More efficient saving of intermediate results (#1529) Don't write old local optima again and again in every iteration to save time and to avoid huge fragmented files. Some other smaller modifications/fixes: * make sacess_history plotting more robust. * Set SacessOptimizer.exit_flag --- pypesto/optimize/ess/__init__.py | 2 +- pypesto/optimize/ess/ess.py | 12 +++-- pypesto/optimize/ess/sacess.py | 70 +++++++++++++++++++++----- pypesto/visualize/optimizer_history.py | 6 +++ test/optimize/test_optimize.py | 3 ++ test/visualize/test_visualize.py | 3 +- 6 files changed, 79 insertions(+), 17 deletions(-) diff --git a/pypesto/optimize/ess/__init__.py b/pypesto/optimize/ess/__init__.py index c5f2d0df4..d4c53756a 100644 --- a/pypesto/optimize/ess/__init__.py +++ b/pypesto/optimize/ess/__init__.py @@ -1,6 +1,6 @@ """Enhanced Scatter Search.""" -from .ess import ESSOptimizer +from .ess import ESSExitFlag, ESSOptimizer from .function_evaluator import ( FunctionEvaluator, FunctionEvaluatorMP, diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index cdcab2f12..2d6f0775e 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -27,17 +27,21 @@ class ESSExitFlag(int, enum.Enum): - """Exit flags used by :class:`ESSOptimizer`.""" + """Scatter search exit flags. + + Exit flags used by :class:`pypesto.ess.ESSOptimizer` and + :class:`pypesto.ess.SacessOptimizer`. + """ # ESS did not run/finish yet DID_NOT_RUN = 0 - # Exited after reaching maximum number of iterations + # Exited after reaching the maximum number of iterations MAX_ITER = -1 # Exited after exhausting function evaluation budget MAX_EVAL = -2 # Exited after exhausting wall-time budget MAX_TIME = -3 - # Termination because for other reason than exit criteria + # Termination because of other reasons than exit criteria ERROR = -99 @@ -242,6 +246,8 @@ def _initialize(self): # Overall best function value found so far self.fx_best: float = np.inf # Results from local searches (only those with finite fval) + # (there is potential to save memory here by only keeping the + # parameters in memory and not the full result) self.local_solutions: list[OptimizerResult] = [] # Index of current iteration self.n_iter: int = 0 diff --git a/pypesto/optimize/ess/sacess.py b/pypesto/optimize/ess/sacess.py index 9340b23da..04428d9e1 100644 --- a/pypesto/optimize/ess/sacess.py +++ b/pypesto/optimize/ess/sacess.py @@ -24,7 +24,10 @@ from ... import MemoryHistory from ...startpoint import StartpointMethod from ...store.read_from_hdf5 import read_result -from ...store.save_to_hdf5 import write_result +from ...store.save_to_hdf5 import ( + OptimizationResultHDF5Writer, + ProblemHDF5Writer, +) from ..optimize import Problem from .ess import ESSExitFlag, ESSOptimizer from .function_evaluator import create_function_evaluator @@ -162,7 +165,7 @@ def __init__( Directory for temporary files. This defaults to a directory in the current working directory named ``SacessOptimizerTemp-{random suffix}``. When setting this option, make sure any optimizers running in - parallel have a unique `tmpdir`. + parallel have a unique `tmpdir`. Expected to be empty. mp_start_method: The start method for the multiprocessing context. See :mod:`multiprocessing` for details. Running `SacessOptimizer` @@ -326,7 +329,9 @@ def minimize( self.histories = [ worker_result.history for worker_result in self.worker_results ] - + self.exit_flag = min( + worker_result.exit_flag for worker_result in self.worker_results + ) result = self._create_result(problem) walltime = time.time() - start_time @@ -649,6 +654,10 @@ def run( ): self._start_time = time.time() + # index of the local solution in ESSOptimizer.local_solutions + # that was most recently saved by _autosave + last_saved_local_solution = -1 + self._logger.setLevel(self._loglevel) # Set the manager logger to one created within the current process self._manager._logger = self._logger @@ -683,20 +692,14 @@ def run( # perform one ESS iteration ess._do_iteration() - if self._tmp_result_file: - # TODO maybe not in every iteration? - ess_results = ess._create_result() - write_result( - ess_results, - self._tmp_result_file, - overwrite=True, - optimize=True, - ) # check if the best solution of the last local ESS is sufficiently # better than the sacess-wide best solution self.maybe_update_best(ess.x_best, ess.fx_best) self._best_known_fx = min(ess.fx_best, self._best_known_fx) + self._autosave(ess, last_saved_local_solution) + last_saved_local_solution = len(ess.local_solutions) - 1 + self._cooperate() self._maybe_adapt(problem) @@ -911,6 +914,49 @@ def abort(self): self._finalize(None) + def _autosave(self, ess: ESSOptimizer, last_saved_local_solution: int): + """Save intermediate results. + + If a temporary result file is set, save the (part of) the current state + of the ESS to that file. + + We save the current best solution and the local optimizer results. + """ + if not self._tmp_result_file: + return + + # save problem in first iteration + if ess.n_iter == 1: + pypesto_problem_writer = ProblemHDF5Writer(self._tmp_result_file) + pypesto_problem_writer.write( + ess.evaluator.problem, overwrite=False + ) + + opt_res_writer = OptimizationResultHDF5Writer(self._tmp_result_file) + for i in range( + last_saved_local_solution + 1, len(ess.local_solutions) + ): + optimizer_result = ess.local_solutions[i] + optimizer_result.id = str(i + ess.n_iter) + opt_res_writer.write_optimizer_result( + optimizer_result, overwrite=False + ) + + # save the current best solution + optimizer_result = pypesto.OptimizerResult( + id=str(len(ess.local_solutions) + ess.n_iter), + x=ess.x_best, + fval=ess.fx_best, + message=f"Global best (iteration {ess.n_iter})", + time=time.time() - ess.starttime, + n_fval=ess.evaluator.n_eval, + optimizer=str(ess), + ) + optimizer_result.update_to_full(ess.evaluator.problem) + opt_res_writer.write_optimizer_result( + optimizer_result, overwrite=False + ) + @staticmethod def get_temp_result_filename(worker_idx: int, tmpdir: str | Path) -> str: return str(Path(tmpdir, f"sacess-{worker_idx:02d}_tmp.h5").absolute()) diff --git a/pypesto/visualize/optimizer_history.py b/pypesto/visualize/optimizer_history.py index 474c030ff..dba5b0200 100644 --- a/pypesto/visualize/optimizer_history.py +++ b/pypesto/visualize/optimizer_history.py @@ -1,4 +1,5 @@ import logging +import warnings from collections.abc import Iterable from typing import Optional, Union @@ -501,6 +502,8 @@ def sacess_history( The plot axes. `ax` or a new axes if `ax` was `None`. """ ax = ax or plt.subplot() + if len(histories) == 0: + warnings.warn("No histories to plot.", stacklevel=2) # plot overall minimum # merge results @@ -530,6 +533,9 @@ def sacess_history( # plot steps of individual workers for worker_idx, history in enumerate(histories): x, y = history.get_time_trace(), history.get_fval_trace() + if len(x) == 0: + warnings.warn(f"No trace for worker #{worker_idx}.", stacklevel=2) + continue # extend from last decrease to last timepoint x = np.append(x, [np.max(t)]) y = np.append(y, [np.min(y)]) diff --git a/test/optimize/test_optimize.py b/test/optimize/test_optimize.py index 6f8260b90..ce83bc069 100644 --- a/test/optimize/test_optimize.py +++ b/test/optimize/test_optimize.py @@ -20,6 +20,7 @@ import pypesto.optimize as optimize from pypesto import Objective from pypesto.optimize.ess import ( + ESSExitFlag, ESSOptimizer, FunctionEvaluatorMP, RefSet, @@ -501,12 +502,14 @@ def test_ess(problem, local_optimizer, ess_type, request): adaptation_sent_coeff=5, ), ) + else: raise ValueError(f"Unsupported ESS type {ess_type}.") res = ess.minimize( problem=problem, ) + assert ess.exit_flag in (ESSExitFlag.MAX_TIME, ESSExitFlag.MAX_ITER) print("ESS result: ", res.summary()) # best values roughly: cr: 4.701; rosen 7.592e-10 diff --git a/test/visualize/test_visualize.py b/test/visualize/test_visualize.py index ecdc38eab..57a4d8ddd 100644 --- a/test/visualize/test_visualize.py +++ b/test/visualize/test_visualize.py @@ -1220,7 +1220,7 @@ def test_time_trajectory_model(): @close_fig def test_sacess_history(): """Test pypesto.visualize.optimizer_history.sacess_history""" - from pypesto.optimize.ess.sacess import SacessOptimizer + from pypesto.optimize.ess.sacess import ESSExitFlag, SacessOptimizer from pypesto.visualize.optimizer_history import sacess_history problem = create_problem() @@ -1228,6 +1228,7 @@ def test_sacess_history(): max_walltime_s=1, num_workers=2, sacess_loglevel=logging.WARNING ) sacess.minimize(problem) + assert sacess.exit_flag == ESSExitFlag.MAX_TIME sacess_history(sacess.histories) From 50f830f87955f9245855c0109ca58b831d5cd6d6 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 10 Dec 2024 18:42:56 +0100 Subject: [PATCH 26/35] Fix get_parameter_prior_dict docstring (#1537) Documentation did not match the code. --- pypesto/objective/priors.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pypesto/objective/priors.py b/pypesto/objective/priors.py index 4ffcdaf6a..9164c70b2 100644 --- a/pypesto/objective/priors.py +++ b/pypesto/objective/priors.py @@ -244,11 +244,13 @@ def get_parameter_prior_dict( "parameterScaleUniform", "parameterScaleNormal", "parameterScaleLaplace"} prior_parameters: - Parameters of the priors. Parameters are defined in linear scale. + Parameters of the priors. Parameters are defined in the parameter + scale if the `prior_type` starts with ``parameterScale``, otherwise in + the linear scale. parameter_scale: scale in which the parameter is defined (since a parameter can be log-transformed, while the prior is always defined in the linear - space, unless it starts with "parameterScale") + space, unless it starts with ``parameterScale``). """ log_f, d_log_f_dx, dd_log_f_ddx, res, d_res_dx = _prior_densities( prior_type, prior_parameters From f2fb5cdc4b32ffcdbec4fbf1b7e48e76181d189a Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 10 Dec 2024 17:48:36 +0000 Subject: [PATCH 27/35] Bump codecov/codecov-action from 4 to 5 (#1508) * Bump codecov/codecov-action from 4 to 5 Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v4...v5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] * Update ci.yml --------- Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Daniel Weindl --- .github/workflows/ci.yml | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f55f8f6b5..72d27d7f3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -55,10 +55,10 @@ jobs: CXX: clang++ - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml mac: runs-on: macos-13 # TODO: change to macos-latest after the next release @@ -91,10 +91,10 @@ jobs: run: ulimit -n 65536 65536 && tox -e base - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml windows: runs-on: windows-latest @@ -163,10 +163,10 @@ jobs: CXX: clang++ - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml julia: runs-on: ubuntu-latest @@ -214,10 +214,10 @@ jobs: run: tox -e julia - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml optimize: runs-on: ubuntu-latest @@ -250,10 +250,10 @@ jobs: run: tox -e optimize - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml hierarchical: runs-on: ubuntu-latest @@ -286,10 +286,10 @@ jobs: run: tox -e hierarchical - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml select: runs-on: ubuntu-latest @@ -322,10 +322,10 @@ jobs: run: tox -e select - name: Coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + files: ./coverage.xml quality: runs-on: ubuntu-latest From 20fd9dec68ff796f5be0aa38de182be4b1d00563 Mon Sep 17 00:00:00 2001 From: Stephan Grein Date: Wed, 18 Dec 2024 16:02:18 +0100 Subject: [PATCH 28/35] Update publish_dockerhub.yml (#1538) * Update publish_dockerhub.yml * Update publish_dockerhub.yml Corrected workflow. --- .github/workflows/publish_dockerhub.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/publish_dockerhub.yml b/.github/workflows/publish_dockerhub.yml index 4b05ad1fc..c0e8fd6c8 100644 --- a/.github/workflows/publish_dockerhub.yml +++ b/.github/workflows/publish_dockerhub.yml @@ -3,7 +3,7 @@ name: Build and Push Docker Image on: push: branches: - - main + - develop paths: - 'docker/**' - '.github/workflows/docker-publish.yml' @@ -25,8 +25,8 @@ jobs: - name: Build and tag the Docker image run: | - docker build -t ICB_DCM/pypesto:latest -f docker/Dockerfile . + docker build -t stephanmg/pypesto:latest -f docker/Dockerfile . - name: Push the Docker image to Docker Hub run: | - docker push ICB_DCM/pypesto:latest + docker push stephanmg/pypesto:latest From e31852b049dac17891a96d38f51a526ba528115a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 8 Jan 2025 11:25:20 +0100 Subject: [PATCH 29/35] SciPy 1.15 compatibility (#1544) Some tests were incompatible with the recently release `scipy==1.15`. Fixed here. --- test/util.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/util.py b/test/util.py index 696db87e4..10bbe4371 100644 --- a/test/util.py +++ b/test/util.py @@ -40,6 +40,8 @@ def obj_for_sensi(fun, grad, hess, max_sensi_order, integrated, x): ret: dict With fields obj, max_sensi_order, x, fval, grad, hess. """ + x = np.asarray(x) + if integrated: if max_sensi_order == 2: From 7c89cd27c4d471b78a70a6d0cff09a695471c483 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 8 Jan 2025 11:26:23 +0100 Subject: [PATCH 30/35] Doc: IpoptOptimizer (#1543) Fix incorrect package URL. Link available options. --- pypesto/optimize/optimizer.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pypesto/optimize/optimizer.py b/pypesto/optimize/optimizer.py index 659303aed..7a6d83488 100644 --- a/pypesto/optimize/optimizer.py +++ b/pypesto/optimize/optimizer.py @@ -549,7 +549,7 @@ def get_default_options(self): class IpoptOptimizer(Optimizer): - """Use IpOpt (https://pypi.org/project/ipopt/) for optimization.""" + """Use Ipopt (https://pypi.org/project/cyipopt/) for optimization.""" def __init__(self, options: dict = None): """ @@ -560,6 +560,9 @@ def __init__(self, options: dict = None): options: Options are directly passed on to `cyipopt.minimize_ipopt`, except for the `approx_grad` option, which is handled separately. + + For a list of available options, see the Ipopt documentation + (https://coin-or.github.io/Ipopt/OPTIONS.html). """ super().__init__() self.approx_grad = False @@ -599,7 +602,7 @@ def minimize( jac = objective.get_grad else: raise ValueError( - "For IPOPT, the objective must either be able to return " + "For Ipopt, the objective must either be able to return " "gradients or the `approx_grad` must be set to True." ) From 11caccd777b4a672ea9eee897dcc688f1b66abf6 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 8 Jan 2025 11:53:19 +0100 Subject: [PATCH 31/35] Bump docker/login-action from 2 to 3 (#1531) Bumps [docker/login-action](https://github.com/docker/login-action) from 2 to 3. - [Release notes](https://github.com/docker/login-action/releases) - [Commits](https://github.com/docker/login-action/compare/v2...v3) --- updated-dependencies: - dependency-name: docker/login-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Paul Jonas Jost <70631928+PaulJonasJost@users.noreply.github.com> --- .github/workflows/publish_dockerhub.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/publish_dockerhub.yml b/.github/workflows/publish_dockerhub.yml index c0e8fd6c8..a98ad41fd 100644 --- a/.github/workflows/publish_dockerhub.yml +++ b/.github/workflows/publish_dockerhub.yml @@ -18,7 +18,7 @@ jobs: uses: actions/checkout@v4 - name: Log in to Docker Hub - uses: docker/login-action@v2 + uses: docker/login-action@v3 with: username: ${{ secrets.DOCKER_USERNAME }} password: ${{ secrets.DOCKER_TOKEN }} From 1121bcc55495a70418ae291857b093e1cf37f038 Mon Sep 17 00:00:00 2001 From: Stephan Grein Date: Wed, 8 Jan 2025 18:18:50 +0100 Subject: [PATCH 32/35] Update parameters.py (#1542) * Update parameters.py Fix issue #1535 * Update parameters.py * Update parameters.py Linting * Update parameters.py Revert changes. * Update read_from_hdf5.py Decode byte str to str. --- pypesto/store/read_from_hdf5.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pypesto/store/read_from_hdf5.py b/pypesto/store/read_from_hdf5.py index 8f011d114..6c1603bfe 100644 --- a/pypesto/store/read_from_hdf5.py +++ b/pypesto/store/read_from_hdf5.py @@ -144,6 +144,7 @@ def read(self, objective: ObjectiveBase = None) -> Problem: problem.x_fixed_vals = [float(val) for val in problem.x_fixed_vals] problem.x_fixed_indices = [int(ix) for ix in problem.x_fixed_indices] problem.x_names = [name.decode() for name in problem.x_names] + problem.x_scales = [scale.decode() for scale in problem.x_scales] return problem From e9a969e5adae4ee5943d364fe1815d7e78dfb4b8 Mon Sep 17 00:00:00 2001 From: Doresic <85789271+Doresic@users.noreply.github.com> Date: Wed, 8 Jan 2025 21:58:41 +0100 Subject: [PATCH 33/35] Fix not scaling inner pars when applying to rdatas (#1534) --- pypesto/hierarchical/relative/solver.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pypesto/hierarchical/relative/solver.py b/pypesto/hierarchical/relative/solver.py index 930fdeba2..6ee1e2612 100644 --- a/pypesto/hierarchical/relative/solver.py +++ b/pypesto/hierarchical/relative/solver.py @@ -238,6 +238,8 @@ def apply_inner_parameters_to_rdatas( """ sim = [rdata["y"] for rdata in rdatas] sigma = [rdata["sigmay"] for rdata in rdatas] + inner_parameters = copy.deepcopy(inner_parameters) + inner_parameters = scale_back_value_dict(inner_parameters, problem) # apply offsets, scalings and sigmas for x in problem.get_xs_for_type(InnerParameterType.SCALING): From 9f0b6f1498367533a674fa7f365dcbacd7030362 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Sat, 11 Jan 2025 03:46:33 +0100 Subject: [PATCH 34/35] Select: use the new "standardized" `petab_select` classes (#1530) * switch to the new `Models` class * use feature branch * fixme: use petab select PR branch * deprecate model_id_binary_postprocessor * update method.py * add `ModelProblem` to `pypesto.select` interface * update for new `end_iteration` and `Problem.state` * fix tests * temp fix reqs * update notebook * update temp fix petab-select * deprecate `pypesto.visualize.select` in favor of `petab_select.plot` * install petab_select[plot] * remove tests for deprecated viz methods * move test cases to pypesto * unfix petab-select * update notebook; clear cell outputs * disable other tests * print debug * increase logging * print debug * update petab-select version * try explicit return * check only famos_cli * fail fast * pass env var; undo debug code * add other tests * review: clarify reloading the problem * review: numpy class attribute docstring typehints * review: typo * fix doc type hint refs, bump petab-select version * add save/load info and postprocessors to notebook --- doc/example/model_selection.ipynb | 283 ++++++++++++-------- doc/example/model_selection/model_space.tsv | 2 +- pypesto/select/__init__.py | 1 + pypesto/select/method.py | 37 ++- pypesto/select/model_problem.py | 5 +- pypesto/select/postprocessors.py | 20 +- pypesto/select/problem.py | 91 +++---- pypesto/visualize/select.py | 201 +++----------- setup.cfg | 7 +- test/select/get_test_cases.sh | 8 + test/select/test_select.py | 107 +++----- test/select/test_test_cases.py | 218 +++++++++++++++ tox.ini | 4 +- 13 files changed, 544 insertions(+), 440 deletions(-) create mode 100755 test/select/get_test_cases.sh create mode 100644 test/select/test_test_cases.py diff --git a/doc/example/model_selection.ipynb b/doc/example/model_selection.ipynb index f31338dc7..5d1ef6800 100644 --- a/doc/example/model_selection.ipynb +++ b/doc/example/model_selection.ipynb @@ -105,28 +105,6 @@ "metadata": {}, "outputs": [], "source": [ - "from petab_select import VIRTUAL_INITIAL_MODEL\n", - "\n", - "\n", - "# Helpers for plotting etc.\n", - "def get_labels(models):\n", - " labels = {\n", - " model.get_hash(): str(model.model_subspace_id)\n", - " for model in models\n", - " if model != VIRTUAL_INITIAL_MODEL\n", - " }\n", - " return labels\n", - "\n", - "\n", - "def get_digraph_labels(models, criterion):\n", - " zero = min(model.get_criterion(criterion) for model in models)\n", - " labels = {\n", - " model.get_hash(): f\"{model.model_subspace_id}\\n{model.get_criterion(criterion) - zero:.2f}\"\n", - " for model in models\n", - " }\n", - " return labels\n", - "\n", - "\n", "# Disable some logged messages that make the model selection log more\n", "# difficult to read.\n", "import tqdm\n", @@ -190,16 +168,18 @@ "outputs": [], "source": [ "# Reduce notebook runtime\n", - "minimize_options = {\n", - " \"n_starts\": 10,\n", - " \"filename\": None,\n", - " \"progress_bar\": False,\n", + "model_problem_options = {\n", + " \"minimize_options\": {\n", + " \"n_starts\": 10,\n", + " \"filename\": None,\n", + " \"progress_bar\": False,\n", + " },\n", "}\n", "\n", "best_model_1, _ = pypesto_select_problem_1.select(\n", " method=Method.FORWARD,\n", " criterion=Criterion.AIC,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", ")" ] }, @@ -221,7 +201,7 @@ "best_model_2, _ = pypesto_select_problem_1.select(\n", " method=Method.FORWARD,\n", " criterion=Criterion.AIC,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", " predecessor_model=best_model_1,\n", ")" ] @@ -239,21 +219,19 @@ "metadata": {}, "outputs": [], "source": [ - "import pypesto.visualize.select as pvs\n", + "import petab_select.plot as plot\n", "\n", "selected_models = [best_model_1, best_model_2]\n", - "ax = pvs.plot_selected_models(\n", - " [best_model_1, best_model_2],\n", - " criterion=Criterion.AIC,\n", - " relative=False,\n", - " labels=get_labels(selected_models),\n", - ")\n", - "ax = pvs.plot_selected_models(\n", - " [best_model_1, best_model_2],\n", + "\n", + "plot_data = plot.PlotData(\n", + " models=petab_select_problem.state.models,\n", " criterion=Criterion.AIC,\n", - " labels=get_labels(selected_models),\n", + " relative_criterion=False,\n", ")\n", - "ax.plot();" + "plot.line_best_by_iteration(plot_data=plot_data)\n", + "\n", + "plot_data.relative_criterion = True\n", + "plot.line_best_by_iteration(plot_data=plot_data);" ] }, { @@ -262,13 +240,12 @@ "metadata": {}, "outputs": [], "source": [ - "pvs.plot_calibrated_models_digraph(\n", - " problem=pypesto_select_problem_1,\n", - " labels=get_digraph_labels(\n", - " pypesto_select_problem_1.calibrated_models.values(),\n", - " criterion=Criterion.AIC,\n", - " ),\n", - ");" + "plot_data = plot.PlotData(\n", + " models=petab_select_problem.state.models,\n", + " criterion=Criterion.AIC,\n", + ")\n", + "plot_data.augment_labels(criterion=True)\n", + "plot.graph_history(plot_data=plot_data);" ] }, { @@ -290,15 +267,22 @@ "import numpy as np\n", "from petab_select import Model\n", "\n", - "petab_select_problem.model_space.reset_exclusions()\n", + "# The PEtab Select problem contains information about calibrated models,\n", + "# so that they are not recalibrated. In order to do a completely new\n", + "# model selection, one can reload the original PEtab Select problem.\n", + "petab_select_problem = petab_select.Problem.from_yaml(\n", + " \"model_selection/petab_select_problem.yaml\"\n", + ")\n", "pypesto_select_problem_2 = pypesto.select.Problem(\n", " petab_select_problem=petab_select_problem\n", ")\n", "\n", - "petab_yaml = \"model_selection/example_modelSelection.yaml\"\n", + "model_subspace_petab_yaml = \"model_selection/example_modelSelection.yaml\"\n", "initial_model = Model(\n", " model_id=\"myModel\",\n", - " petab_yaml=petab_yaml,\n", + " model_subspace_petab_yaml=model_subspace_petab_yaml,\n", + " model_subspace_id=\"dummy_myModel\",\n", + " model_subspace_indices=[0,0,0],\n", " parameters={\n", " \"k1\": 0.1,\n", " \"k2\": ESTIMATE,\n", @@ -321,7 +305,7 @@ " method=Method.BACKWARD,\n", " criterion=Criterion.AIC,\n", " predecessor_model=initial_model,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", ");" ] }, @@ -331,18 +315,12 @@ "metadata": {}, "outputs": [], "source": [ - "initial_model_label = {initial_model.get_hash(): initial_model.model_id}\n", - "\n", - "pvs.plot_calibrated_models_digraph(\n", - " problem=pypesto_select_problem_2,\n", - " labels={\n", - " **get_digraph_labels(\n", - " pypesto_select_problem_2.calibrated_models.values(),\n", - " criterion=Criterion.AIC,\n", - " ),\n", - " **initial_model_label,\n", - " },\n", - ");" + "plot_data = plot.PlotData(\n", + " models=petab_select_problem.state.models,\n", + " criterion=Criterion.AIC,\n", + ")\n", + "plot_data.augment_labels(criterion=True)\n", + "plot.graph_history(plot_data=plot_data);" ] }, { @@ -374,7 +352,9 @@ "metadata": {}, "outputs": [], "source": [ - "petab_select_problem.model_space.reset_exclusions()\n", + "petab_select_problem = petab_select.Problem.from_yaml(\n", + " \"model_selection/petab_select_problem.yaml\"\n", + ")\n", "pypesto_select_problem_3 = pypesto.select.Problem(\n", " petab_select_problem=petab_select_problem\n", ")\n", @@ -383,7 +363,7 @@ " criterion=Criterion.BIC,\n", " select_first_improvement=True,\n", " startpoint_latest_mle=True,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", ")" ] }, @@ -393,11 +373,11 @@ "metadata": {}, "outputs": [], "source": [ - "pvs.plot_selected_models(\n", - " selected_models=best_models,\n", + "plot_data = plot.PlotData(\n", + " petab_select_problem.state.models,\n", " criterion=Criterion.BIC,\n", - " labels=get_labels(best_models),\n", - ");" + ")\n", + "plot.line_best_by_iteration(plot_data=plot_data);" ] }, { @@ -406,15 +386,18 @@ "metadata": {}, "outputs": [], "source": [ - "pvs.plot_calibrated_models_digraph(\n", - " problem=pypesto_select_problem_3,\n", - " criterion=Criterion.BIC,\n", - " relative=False,\n", - " labels=get_digraph_labels(\n", - " pypesto_select_problem_3.calibrated_models.values(),\n", - " criterion=Criterion.BIC,\n", - " ),\n", - ");" + "plot_data.relative_criterion = False\n", + "plot_data.augment_labels(criterion=True)\n", + "plot.graph_history(plot_data=plot_data);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Postprocessors\n", + "\n", + "You can optionally provide postprocessing methods that operate on calibrated models. For example, there are postprocessors to save the pyPESTO result to disk, or save a summarized report of models to a TSV file. You can also write your own method, which should take only a single `pypesto.select.ModelProblem` object as input." ] }, { @@ -424,17 +407,42 @@ "outputs": [], "source": [ "# Repeat with AICc and criterion_threshold == 10\n", - "petab_select_problem.model_space.reset_exclusions()\n", + "petab_select_problem = petab_select.Problem.from_yaml(\n", + " \"model_selection/petab_select_problem.yaml\"\n", + ")\n", "pypesto_select_problem_4 = pypesto.select.Problem(\n", " petab_select_problem=petab_select_problem\n", - ")\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Setup a postprocessor that saves the pyPESTO result objects to disk, as well as a summarized report.\n", + "from functools import partial\n", + "from pypesto.select.postprocessors import save_postprocessor, report_postprocessor, multi_postprocessor\n", + "from pathlib import Path\n", + "import shutil\n", + "output_path = Path() / \"output_select\"\n", + "if output_path.exists():\n", + " shutil.rmtree(str(output_path))\n", + "output_path.mkdir(exist_ok=True, parents=True)\n", + "model_postprocessor = partial(multi_postprocessor, postprocessors=[\n", + " partial(save_postprocessor, output_path=output_path),\n", + " partial(report_postprocessor, output_filepath=output_path / \"report.tsv\"),\n", + "])\n", + "\n", "best_models = pypesto_select_problem_4.select_to_completion(\n", " method=Method.FORWARD,\n", " criterion=Criterion.AICC,\n", " select_first_improvement=True,\n", " startpoint_latest_mle=True,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", " criterion_threshold=10,\n", + " model_postprocessor=model_postprocessor,\n", ")" ] }, @@ -444,11 +452,41 @@ "metadata": {}, "outputs": [], "source": [ - "pvs.plot_selected_models(\n", - " selected_models=best_models,\n", + "df = pd.read_csv(output_path / \"report.tsv\", sep=\"\\t\")\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_data = plot.PlotData(\n", + " models=petab_select_problem.state.models,\n", " criterion=Criterion.AICC,\n", - " labels=get_labels(best_models),\n", - ");" + ")\n", + "plot.line_best_by_iteration(plot_data=plot_data);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_data.relative_criterion = False\n", + "plot_data.augment_labels(criterion=True)\n", + "plot.graph_history(plot_data=plot_data);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Saving results\n", + "\n", + "Although pyPESTO results and summaries can be saved to disk via the postprocessors, most analysis methods will expect a `petab_select.Models` object, which is created during model selection and can also be saved to disk for later analysis." ] }, { @@ -457,14 +495,42 @@ "metadata": {}, "outputs": [], "source": [ - "pvs.plot_calibrated_models_digraph(\n", - " problem=pypesto_select_problem_4,\n", + "petab_select_problem.state.models.to_yaml(output_path / \"my_models.yaml\")\n", + "loaded_models = petab_select.Models.from_yaml(output_path / \"my_models.yaml\")\n", + "\n", + "plot_data = plot.PlotData(\n", + " models=loaded_models,\n", " criterion=Criterion.AICC,\n", - " relative=False,\n", - " labels=get_digraph_labels(\n", - " pypesto_select_problem_4.calibrated_models.values(),\n", - " criterion=Criterion.AICC,\n", - " ),\n", + " relative_criterion=False\n", + ")\n", + "plot_data.augment_labels(criterion=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Additional plotting methods\n", + "\n", + "There are additional plotting methods provided by the PEtab Select library, which are demonstrated in another notebook: https://petab-select.readthedocs.io/en/stable/examples/visualization.html\n", + "\n", + "An example is this ordering of models according to the iteration in which they were calibrated. N.B. The predecessor of M1_7 is actually M1_1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot.graph_iteration_layers(\n", + " plot_data=plot_data,\n", + " draw_networkx_kwargs={\n", + " \"arrowstyle\": \"-|>\",\n", + " \"node_shape\": \"s\",\n", + " \"node_size\": 4000,\n", + " \"edgecolors\": \"k\",\n", + " },\n", ");" ] }, @@ -482,14 +548,18 @@ "metadata": {}, "outputs": [], "source": [ - "petab_select_problem.model_space.reset_exclusions()\n", + "petab_select_problem = petab_select.Problem.from_yaml(\n", + " \"model_selection/petab_select_problem.yaml\"\n", + ")\n", "pypesto_select_problem_5 = pypesto.select.Problem(\n", " petab_select_problem=petab_select_problem\n", ")\n", "\n", "initial_model_1 = Model(\n", " model_id=\"myModel1\",\n", - " petab_yaml=petab_yaml,\n", + " model_subspace_petab_yaml=model_subspace_petab_yaml,\n", + " model_subspace_id=\"dummy_myModel1\",\n", + " model_subspace_indices=[0,0,0],\n", " parameters={\n", " \"k1\": 0,\n", " \"k2\": 0,\n", @@ -500,7 +570,9 @@ "\n", "initial_model_2 = Model(\n", " model_id=\"myModel2\",\n", - " petab_yaml=petab_yaml,\n", + " model_subspace_petab_yaml=model_subspace_petab_yaml,\n", + " model_subspace_id=\"dummy_myModel2\",\n", + " model_subspace_indices=[0,0,0],\n", " parameters={\n", " \"k1\": ESTIMATE,\n", " \"k2\": ESTIMATE,\n", @@ -514,7 +586,7 @@ " method=Method.FORWARD,\n", " criterion=Criterion.AIC,\n", " predecessor_models=initial_models,\n", - " minimize_options=minimize_options,\n", + " model_problem_options=model_problem_options,\n", ")" ] }, @@ -524,23 +596,12 @@ "metadata": {}, "outputs": [], "source": [ - "initial_model_labels = {\n", - " initial_model.get_hash(): initial_model.model_id\n", - " for initial_model in initial_models\n", - "}\n", - "\n", - "pvs.plot_calibrated_models_digraph(\n", - " problem=pypesto_select_problem_5,\n", + "plot_data = plot.PlotData(\n", + " models=petab_select_problem.state.models,\n", " criterion=Criterion.AIC,\n", - " relative=False,\n", - " labels={\n", - " **get_digraph_labels(\n", - " pypesto_select_problem_5.calibrated_models.values(),\n", - " criterion=Criterion.AIC,\n", - " ),\n", - " **initial_model_labels,\n", - " },\n", - ");" + " relative_criterion=False,\n", + ")\n", + "plot.graph_history(plot_data);" ] } ], diff --git a/doc/example/model_selection/model_space.tsv b/doc/example/model_selection/model_space.tsv index 42e3f1f54..d4eb7ec42 100644 --- a/doc/example/model_selection/model_space.tsv +++ b/doc/example/model_selection/model_space.tsv @@ -1,4 +1,4 @@ -model_subspace_id petab_yaml k1 k2 k3 +model_subspace_id model_subspace_petab_yaml k1 k2 k3 M1_0 example_modelSelection.yaml 0 0 0 M1_1 example_modelSelection.yaml 0.2 0.1 estimate M1_2 example_modelSelection.yaml 0.2 estimate 0 diff --git a/pypesto/select/__init__.py b/pypesto/select/__init__.py index 6fc388026..aa927bc39 100644 --- a/pypesto/select/__init__.py +++ b/pypesto/select/__init__.py @@ -8,6 +8,7 @@ from . import postprocessors from .misc import SacessMinimizeMethod, model_to_pypesto_problem +from .model_problem import ModelProblem from .problem import Problem try: diff --git a/pypesto/select/method.py b/pypesto/select/method.py index d1ff12443..ec30eea97 100644 --- a/pypesto/select/method.py +++ b/pypesto/select/method.py @@ -3,7 +3,7 @@ import logging from dataclasses import dataclass from enum import Enum -from typing import Any, Callable, Optional +from typing import Any, Callable import numpy as np import petab_select @@ -17,6 +17,7 @@ Criterion, Method, Model, + Models, ) from ..problem import Problem @@ -206,8 +207,7 @@ class MethodCaller: example, in `ForwardSelector`, test models are compared to the previously selected model. calibrated_models: - The calibrated models of the model selection, as a `dict` where keys - are model hashes and values are models. + All calibrated models of the model selection. limit: Limit the number of calibrated models. NB: the number of accepted models may (likely) be fewer. @@ -233,7 +233,7 @@ class MethodCaller: def __init__( self, petab_select_problem: petab_select.Problem, - calibrated_models: dict[str, Model], + calibrated_models: Models, # Arguments/attributes that can simply take the default value here. criterion_threshold: float = 0.0, limit: int = np.inf, @@ -266,11 +266,9 @@ def __init__( self.select_first_improvement = select_first_improvement self.startpoint_latest_mle = startpoint_latest_mle - self.user_calibrated_models = {} + self.user_calibrated_models = Models() if user_calibrated_models is not None: - self.user_calibrated_models = { - model.get_hash(): model for model in user_calibrated_models - } + self.user_calibrated_models = user_calibrated_models self.logger = MethodLogger() @@ -351,7 +349,7 @@ def __init__( # May have changed from `None` to `petab_select.VIRTUAL_INITIAL_MODEL` self.predecessor_model = self.candidate_space.get_predecessor_model() - def __call__(self) -> tuple[list[Model], dict[str, Model]]: + def __call__(self) -> tuple[Model, Models]: """Run a single iteration of the model selection method. A single iteration here refers to calibration of all candidate models. @@ -365,8 +363,7 @@ def __call__(self) -> tuple[list[Model], dict[str, Model]]: A 2-tuple, with the following values: 1. the predecessor model for the newly calibrated models; and - 2. the newly calibrated models, as a `dict` where keys are model - hashes and values are models. + 2. the newly calibrated models. """ # All calibrated models in this iteration (see second return value). self.logger.new_selection() @@ -384,7 +381,7 @@ def __call__(self) -> tuple[list[Model], dict[str, Model]]: # TODO parallelize calibration (maybe not sensible if # `self.select_first_improvement`) - calibrated_models = {} + calibrated_models = Models() for model in iteration[UNCALIBRATED_MODELS]: if ( model.get_criterion( @@ -405,7 +402,7 @@ def __call__(self) -> tuple[list[Model], dict[str, Model]]: else: self.new_model_problem(model=model) - calibrated_models[model.get_hash()] = model + calibrated_models.append(model) method_signal = self.handle_calibrated_model( model=model, predecessor_model=iteration[PREDECESSOR_MODEL], @@ -414,18 +411,19 @@ def __call__(self) -> tuple[list[Model], dict[str, Model]]: break iteration_results = petab_select.ui.end_iteration( + problem=self.petab_select_problem, candidate_space=iteration[CANDIDATE_SPACE], calibrated_models=calibrated_models, ) - self.calibrated_models.update(iteration_results[MODELS]) + self.calibrated_models += iteration_results[MODELS] - return iteration[PREDECESSOR_MODEL], iteration_results[MODELS] + return iteration_results[MODELS] def handle_calibrated_model( self, model: Model, - predecessor_model: Optional[Model], + predecessor_model: Model, ) -> MethodSignal: """Handle the model selection method, given a new calibrated model. @@ -454,8 +452,7 @@ def handle_calibrated_model( # Reject the model if it doesn't improve on the predecessor model. if ( - predecessor_model is not None - and predecessor_model != VIRTUAL_INITIAL_MODEL + predecessor_model.hash != VIRTUAL_INITIAL_MODEL.hash and not self.model1_gt_model0( model1=model, model0=predecessor_model ) @@ -549,7 +546,9 @@ def new_model_problem( predecessor_model = self.calibrated_models[ model.predecessor_model_hash ] - if str(model.petab_yaml) != str(predecessor_model.petab_yaml): + if str(model.model_subspace_petab_yaml) != str( + predecessor_model.model_subspace_petab_yaml + ): raise NotImplementedError( "The PEtab YAML files differ between the model and its " "predecessor model. This may imply different (fixed union " diff --git a/pypesto/select/model_problem.py b/pypesto/select/model_problem.py index fad3487b2..3872cd9d5 100644 --- a/pypesto/select/model_problem.py +++ b/pypesto/select/model_problem.py @@ -15,6 +15,9 @@ TYPE_POSTPROCESSOR = Callable[["ModelProblem"], None] # noqa: F821 +__all__ = ["ModelProblem"] + + class ModelProblem: """Handles all required calibration tasks on a model. @@ -149,7 +152,7 @@ def minimize(self) -> Result: if isinstance(self.minimize_method, SacessMinimizeMethod): return self.minimize_method( self.pypesto_problem, - model_hash=self.model.get_hash(), + model_hash=self.model.hash, **self.minimize_options, ) return self.minimize_method( diff --git a/pypesto/select/postprocessors.py b/pypesto/select/postprocessors.py index ae8a9d587..b975065a5 100644 --- a/pypesto/select/postprocessors.py +++ b/pypesto/select/postprocessors.py @@ -1,10 +1,11 @@ """Process a model selection :class:`ModelProblem` after calibration.""" +import warnings from pathlib import Path import matplotlib.pyplot as plt import numpy as np -from petab_select.constants import ESTIMATE, TYPE_PATH, Criterion +from petab_select.constants import TYPE_PATH, Criterion from .. import store, visualize from .model_problem import TYPE_POSTPROCESSOR, ModelProblem @@ -48,7 +49,7 @@ def waterfall_plot_postprocessor( See :meth:`save_postprocessor` for usage hints and argument documentation. """ visualize.waterfall(problem.minimize_result) - plot_output_path = Path(output_path) / (problem.model.model_hash + ".png") + plot_output_path = Path(output_path) / (str(problem.model.hash) + ".png") plt.savefig(str(plot_output_path)) @@ -85,7 +86,7 @@ def save_postprocessor( """ stem = problem.model.model_id if use_model_hash: - stem = problem.model.get_hash() + stem = str(problem.model.hash) store.write_result( problem.minimize_result, Path(output_path) / (stem + ".hdf5"), @@ -109,10 +110,15 @@ def model_id_binary_postprocessor(problem: ModelProblem): problem: A model selection :class:`ModelProblem` that has been optimized. """ - model_id = "M_" - for parameter_value in problem.model.parameters.values(): - model_id += "1" if parameter_value == ESTIMATE else "0" - problem.model.model_id = model_id + warnings.warn( + ( + "This is obsolete. Model IDs are by default the model hash, which " + "is now similar to the binary string." + ), + DeprecationWarning, + stacklevel=2, + ) + problem.model.model_id = str(problem.model.hash) def report_postprocessor( diff --git a/pypesto/select/problem.py b/pypesto/select/problem.py index e0bd08cf3..a0034e218 100644 --- a/pypesto/select/problem.py +++ b/pypesto/select/problem.py @@ -1,11 +1,10 @@ """Manage all components of a pyPESTO model selection problem.""" import warnings -from collections.abc import Iterable from typing import Any, Optional import petab_select -from petab_select import Model +from petab_select import Model, Models from .method import MethodCaller from .model_problem import TYPE_POSTPROCESSOR, ModelProblem # noqa: F401 @@ -20,18 +19,18 @@ class Problem: Attributes ---------- - calibrated_models: - Storage for all calibrated models. A dictionary, where keys are - model hashes, and values are :class:`petab_select.Model` objects. - newly_calibrated_models: - Storage for models that were calibrated in the previous iteration of - model selection. Same type as ``calibrated_models``. - method_caller: - A :class:`MethodCaller`, used to run a single iteration of a model + calibrated_models : petab_select.Models + All calibrated models. + newly_calibrated_models : petab_select.Models + All models that were calibrated in the latest iteration of model + selection. + method_caller : pypesto.select.method.MethodCaller + Used to run a single iteration of a model selection method. - model_problem_options: - Passed to the constructor of :class:``ModelProblem``. - petab_select_problem: + model_problem_options : dict[str, typing.Any] + Keyword arguments, passed to the constructor of + :class:`ModelProblem`. + petab_select_problem : petab_select.Problem A PEtab Select problem. """ @@ -60,8 +59,8 @@ def __init__( self.model_problem_options["postprocessor"] = model_postprocessor self.set_state( - calibrated_models={}, - newly_calibrated_models={}, + calibrated_models=Models(), + newly_calibrated_models=Models(), ) # TODO default caller, based on petab_select.Problem @@ -90,8 +89,8 @@ def create_method_caller(self, **kwargs) -> MethodCaller: def set_state( self, - calibrated_models: dict[str, Model], - newly_calibrated_models: dict[str, Model], + calibrated_models: Models, + newly_calibrated_models: Models, ) -> None: """Set the state of the problem. @@ -102,7 +101,7 @@ def set_state( def update_with_newly_calibrated_models( self, - newly_calibrated_models: Optional[dict[str, Model]] = None, + newly_calibrated_models: Optional[Models] = None, ) -> None: """Update the state of the problem with newly calibrated models. @@ -111,7 +110,7 @@ def update_with_newly_calibrated_models( See attributes of :class:`Problem`. """ self.newly_calibrated_models = newly_calibrated_models - self.calibrated_models.update(self.newly_calibrated_models) + self.calibrated_models += self.newly_calibrated_models def handle_select_kwargs( self, @@ -132,7 +131,7 @@ def handle_select_kwargs( def select( self, **kwargs, - ) -> tuple[Model, dict[str, Model], dict[str, Model]]: + ) -> tuple[Model, Models]: """Run a single iteration of a model selection algorithm. The result is the selected model for the current run, independent of @@ -142,21 +141,14 @@ def select( Returns ------- - A 3-tuple, with the following values: + A 2-tuple, with the following values from this iteration: - 1. the best model; - 2. all candidate models in this iteration, as a `dict` with - model hashes as keys and models as values; and - 3. all candidate models from all iterations, as a `dict` with - model hashes as keys and models as values. + 1. the best model; and + 2. all models. """ - # TODO move some options to PEtab Select? e.g.: - # - startpoint_latest_mle - # - select_first_improvement self.handle_select_kwargs(kwargs) - # TODO handle bidirectional method_caller = self.create_method_caller(**kwargs) - previous_best_model, newly_calibrated_models = method_caller() + newly_calibrated_models = method_caller() self.update_with_newly_calibrated_models( newly_calibrated_models=newly_calibrated_models, @@ -164,56 +156,45 @@ def select( best_model = petab_select.ui.get_best( problem=self.petab_select_problem, - models=self.newly_calibrated_models.values(), + models=self.newly_calibrated_models, criterion=method_caller.criterion, ) - # TODO: Reconsider return value. `result` could be stored in attribute, - # then no values need to be returned, and users can request values - # manually. return best_model, newly_calibrated_models def select_to_completion( self, **kwargs, - ) -> list[Model]: - """Run an algorithm until an exception `StopIteration` is raised. + ) -> Models: + """Perform model selection until the method terminates. ``kwargs`` are passed to the :class:`MethodCaller` constructor. - An exception ``StopIteration`` is raised by - :meth:`pypesto.select.method.MethodCaller.__call__` when no candidate models - are found. - Returns ------- - The best models (the best model at each iteration). + All models. """ - best_models = [] + calibrated_models = Models(problem=self.petab_select_problem) self.handle_select_kwargs(kwargs) method_caller = self.create_method_caller(**kwargs) while True: try: - previous_best_model, newly_calibrated_models = method_caller() + iteration_calibrated_models = method_caller() self.update_with_newly_calibrated_models( - newly_calibrated_models=newly_calibrated_models, + newly_calibrated_models=iteration_calibrated_models, ) - best_models.append(previous_best_model) + calibrated_models += iteration_calibrated_models except StopIteration: - previous_best_model = ( - method_caller.candidate_space.predecessor_model - ) - best_models.append(previous_best_model) break - return best_models + return calibrated_models # TODO method that automatically generates initial models, for a specific # number of starts. TODO parallelise? def multistart_select( self, - predecessor_models: Iterable[Model] = None, + predecessor_models: Models = None, **kwargs, ) -> tuple[Model, list[Model]]: """Run an algorithm multiple times, with different predecessor models. @@ -247,10 +228,10 @@ def multistart_select( method_caller = self.create_method_caller( **(kwargs | {"predecessor_model": predecessor_model}) ) - (best_model, models) = method_caller() - self.calibrated_models |= models + models = method_caller() + self.calibrated_models += models - model_lists.append(list(models.values())) + model_lists.append(models) method_caller.candidate_space.reset() best_model = petab_select.ui.get_best( diff --git a/pypesto/visualize/select.py b/pypesto/visualize/select.py index 0a8efde34..e77126184 100644 --- a/pypesto/visualize/select.py +++ b/pypesto/visualize/select.py @@ -1,20 +1,18 @@ -"""Visualization routines for model selection with pyPESTO.""" +"""Deprecated. Use ``petab_select.plot`` methods instead.""" +import warnings + import matplotlib.pyplot as plt -import networkx as nx -from petab_select import Model -from petab_select.constants import VIRTUAL_INITIAL_MODEL, Criterion +import petab_select.plot +from petab_select import Criterion, Model, Models from .. import select as pypesto_select -# TODO move methods to petab_select -RELATIVE_LABEL_FONTSIZE = -2 - def default_label_maker(model: Model) -> str: """Create a model label, for plotting.""" - return model.model_hash[:4] + return str(model.hash)[:4] # FIXME supply the problem to automatically detect the correct criterion? @@ -27,89 +25,19 @@ def plot_selected_models( labels: dict[str, str] = None, ax: plt.Axes = None, ) -> plt.Axes: - """Plot criterion for calibrated models. - - Parameters - ---------- - selected_models: - First result in tuple returned by `ModelSelector.select()`. - criterion: - Key of criterion value in `selected_models` to be plotted. - relative: - If `True`, criterion values are plotted relative to the lowest - criterion value. TODO is the lowest value, always the best? May not - be for different criterion. - fz: - fontsize - size: - Figure size in inches. - labels: - A dictionary of model labels, where keys are model hashes, and - values are model labels, for plotting. If a model label is not - provided, it will be generated from its model ID. - ax: - The axis to use for plotting. - - Returns - ------- - matplotlib.pyplot.Axes - The plot axis. - """ - zero = 0 - if relative: - zero = selected_models[-1].get_criterion(criterion) - - if labels is None: - labels = {} - - # FIGURE - if ax is None: - _, ax = plt.subplots(figsize=size) - linewidth = 3 - - selected_models = [ - model for model in selected_models if model != VIRTUAL_INITIAL_MODEL - ] - - criterion_values = { - labels.get( - model.get_hash(), default_label_maker(model) - ): model.get_criterion(criterion) - zero - for model in selected_models - } - - ax.plot( - criterion_values.keys(), - criterion_values.values(), - linewidth=linewidth, - color="lightgrey", - # edgecolor='k' - ) - - ax.get_xticks() - ax.set_xticks(list(range(len(criterion_values)))) - ax.set_ylabel( - criterion + ("(relative)" if relative else "(absolute)"), fontsize=fz - ) - # could change to compared_model_ids, if all models are plotted - ax.set_xticklabels( - criterion_values.keys(), - fontsize=fz + RELATIVE_LABEL_FONTSIZE, + """Use `petab_select.plot.line_best_by_iteration`` instead. Deprecated.""" + warnings.warn( + "Use `petab_select.plot.line_best_by_iteration` instead.", + DeprecationWarning, + stacklevel=2, ) - for tick in ax.yaxis.get_major_ticks(): - tick.label1.set_fontsize(fz + RELATIVE_LABEL_FONTSIZE) - ytl = ax.get_yticks() - ax.set_ylim([min(ytl), max(ytl)]) - # removing top and right borders - ax.spines["top"].set_visible(False) - ax.spines["right"].set_visible(False) - return ax - - -def plot_history_digraph(*args, **kwargs): - """Ignore me. Renamed to `plot_calibrated_models_digraph`.""" - raise NotImplementedError( - "This method is now `plot_calibrated_models_digraph`." + return petab_select.plot.line_best_by_iteration( + models=Models(selected_models), + criterion=criterion, + relative=relative, + fz=fz, + labels=labels, + ax=ax, ) @@ -123,82 +51,17 @@ def plot_calibrated_models_digraph( labels: dict[str, str] = None, ax: plt.Axes = None, ) -> plt.Axes: - """Plot all calibrated models in the model space, as a directed graph. - - TODO replace magic numbers with options/constants - - Parameters - ---------- - problem: - The pyPESTO Select problem. - calibrated_models: - The models calibrated during model selection, in the format of - `pypesto.select.Problem.calibrated_models`. - criterion: - The criterion. - optimal_distance: - See docs for argument `k` in `networkx.spring_layout`. - relative: - If `True`, criterion values are offset by the minimum criterion - value. - options: - Additional keyword arguments for `networkx.draw_networkx`. - labels: - A dictionary of model labels, where keys are model hashes, and - values are model labels, for plotting. If a model label is not - provided, it will be generated from its model ID. - ax: - The axis to use for plotting. - - Returns - ------- - matplotlib.pyplot.Axes - The plot axis. - """ - if criterion is None: - criterion = problem.petab_select_problem.criterion - if calibrated_models is None: - calibrated_models = problem.calibrated_models - if labels is None: - labels = {} - - G = nx.DiGraph() - edges = [] - for model in calibrated_models.values(): - predecessor_model_hash = model.predecessor_model_hash - if predecessor_model_hash is not None: - from_ = labels.get(predecessor_model_hash, predecessor_model_hash) - # may only not be the case for - # COMPARED_MODEL_ID == INITIAL_VIRTUAL_MODEL - if predecessor_model_hash in calibrated_models: - predecessor_model = calibrated_models[predecessor_model_hash] - from_ = labels.get( - predecessor_model.get_hash(), - default_label_maker(predecessor_model), - ) - else: - raise NotImplementedError( - "Plots for models with `None` as their predecessor model are " - "not yet implemented." - ) - from_ = "None" - to = labels.get(model.get_hash(), default_label_maker(model)) - edges.append((from_, to)) - - G.add_edges_from(edges) - default_options = { - "node_color": "lightgrey", - "arrowstyle": "-|>", - "node_shape": "s", - "node_size": 2500, - } - if options is not None: - default_options.update(options) - - if ax is None: - _, ax = plt.subplots(figsize=(12, 12)) - - pos = nx.spring_layout(G, k=optimal_distance, iterations=20) - nx.draw_networkx(G, pos, ax=ax, **default_options) - - return ax + """Use `petab_select.plot.graph_history`` instead. Deprecated.""" + warnings.warn( + "Use `petab_select.plot.graph_history` instead.", + DeprecationWarning, + stacklevel=2, + ) + return petab_select.plot.graph_history( + models=calibrated_models.values(), + criterion=criterion, + draw_networkx_kwargs=options, + relative=relative, + labels=labels, + ax=ax, + ) diff --git a/setup.cfg b/setup.cfg index 9019fe480..071c44248 100644 --- a/setup.cfg +++ b/setup.cfg @@ -164,12 +164,7 @@ example = ipywidgets >= 8.1.5 benchmark_models_petab @ git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master#subdirectory=src/python select = - # Remove when vis is moved to PEtab Select - networkx >= 2.5.1 - # End remove - # TODO: pinned until https://github.com/ICB-DCM/pyPESTO/pull/1530 is done - # petab-select @ git+https://github.com/PEtab-dev/petab_select.git@develop - petab-select == 0.2.1 + petab-select[plot] >= 0.3.3 test = pytest >= 5.4.3 pytest-cov >= 2.10.0 diff --git a/test/select/get_test_cases.sh b/test/select/get_test_cases.sh new file mode 100755 index 000000000..639f8ed94 --- /dev/null +++ b/test/select/get_test_cases.sh @@ -0,0 +1,8 @@ +#!/bin/bash +if [ ! -d petab_select ] +then + git clone -n --depth=1 --filter=tree:0 https://github.com/PEtab-dev/petab_select.git + cd petab_select + git sparse-checkout set --no-cone /test_cases + git checkout +fi diff --git a/test/select/test_select.py b/test/select/test_select.py index bdec69452..40553f7ed 100644 --- a/test/select/test_select.py +++ b/test/select/test_select.py @@ -15,6 +15,7 @@ Criterion, Method, Model, + get_best_by_iteration, ) import pypesto.engine @@ -80,7 +81,9 @@ def initial_models(petab_problem_yaml) -> list[Model]: """Models that can be used to initialize a search.""" initial_model_1 = Model( model_id="myModel1", - petab_yaml=petab_problem_yaml, + model_subspace_id="dummy", + model_subspace_indices=[], + model_subspace_petab_yaml=petab_problem_yaml, parameters={ "k1": 0, "k2": 0, @@ -90,7 +93,9 @@ def initial_models(petab_problem_yaml) -> list[Model]: ) initial_model_2 = Model( model_id="myModel2", - petab_yaml=petab_problem_yaml, + model_subspace_id="dummy", + model_subspace_indices=[], + model_subspace_petab_yaml=petab_problem_yaml, parameters={ "k1": ESTIMATE, "k2": ESTIMATE, @@ -183,7 +188,7 @@ def test_problem_select_to_completion(pypesto_select_problem): method=Method.FORWARD ) ) - best_models = pypesto_select_problem.select_to_completion( + models = pypesto_select_problem.select_to_completion( criterion=Criterion.BIC, select_first_improvement=True, startpoint_latest_mle=True, @@ -191,63 +196,53 @@ def test_problem_select_to_completion(pypesto_select_problem): candidate_space=candidate_space, ) - expected_calibrated_models_subspace_ids = { + expected_calibrated_ids = { "M1_0", "M1_1", "M1_4", "M1_5", "M1_7", } - test_calibrated_models_subspace_ids = { - model.model_subspace_id - for model in pypesto_select_problem.calibrated_models.values() - } + test_calibrated_ids = {model.model_subspace_id for model in models} # Expected models were calibrated during the search. - assert ( - test_calibrated_models_subspace_ids - == expected_calibrated_models_subspace_ids + assert test_calibrated_ids == expected_calibrated_ids + + best_by_iteration = get_best_by_iteration( + models=models, criterion=Criterion.BIC ) - expected_best_model_subspace_ids = [ - # The first iteration is from a virtual model, which will appear - # as the best model for that iteration. - VIRTUAL_INITIAL_MODEL, + expected_best_by_iteration = [ "M1_0", "M1_1", - # This iteration with models `{'M1_4', 'M1_5'}` didn't have a better - # model than the previous iteration. - "M1_1", + "M1_5", "M1_7", ] - test_best_model_subspace_ids = [ - (model.model_subspace_id if model != VIRTUAL_INITIAL_MODEL else model) - for model in best_models + test_best_by_iteration = [ + model.model_subspace_id for model in best_by_iteration.values() ] # Expected best models were found. - assert test_best_model_subspace_ids == expected_best_model_subspace_ids + assert test_best_by_iteration == expected_best_by_iteration - expected_best_model_criterion_values = [ - # VIRTUAL_INITIAL_MODEL - np.inf, + expected_best_values = [ 36.767, -4.592, # This iteration with models `{'M1_4', 'M1_5'}` didn't have a better # model than the previous iteration. - -4.592, + -3.330, -4.889, ] - test_best_model_criterion_values = [ + test_best_values = [ ( model.get_criterion(Criterion.BIC) if model != VIRTUAL_INITIAL_MODEL else np.inf ) - for model in best_models + for model in best_by_iteration.values() ] # The best models have the expected criterion values. assert np.isclose( - test_best_model_criterion_values, - expected_best_model_criterion_values, + test_best_values, + expected_best_values, **tolerances, ).all() @@ -268,7 +263,8 @@ def test_problem_multistart_select(pypesto_select_problem, initial_models): assert test_best_model_subspace_id == expected_best_model_subspace_id best_models = [ - pypesto_select_problem.petab_select_problem.get_best( + petab_select.ui.get_best( + problem=pypesto_select_problem.petab_select_problem, models=models, criterion=criterion, ) @@ -294,7 +290,7 @@ def test_problem_multistart_select(pypesto_select_problem, initial_models): ) initial_model_id_hash_map = { - initial_model.model_id: initial_model.get_hash() + initial_model.model_id: initial_model.hash for initial_model in initial_models } expected_predecessor_model_hashes = { @@ -305,7 +301,7 @@ def test_problem_multistart_select(pypesto_select_problem, initial_models): } test_predecessor_model_hashes = { model.model_subspace_id: model.predecessor_model_hash - for model in pypesto_select_problem.calibrated_models.values() + for model in pypesto_select_problem.calibrated_models } # All calibrated models have the expected predecessor model. assert test_predecessor_model_hashes == expected_predecessor_model_hashes @@ -344,7 +340,7 @@ def test_postprocessors(petab_select_problem): expected_newly_calibrated_models_subspace_ids = ["M1_0"] test_newly_calibrated_models_subspace_ids = [ - model.model_subspace_id for model in newly_calibrated_models_1.values() + model.model_subspace_id for model in newly_calibrated_models_1 ] # The expected "forward" models were found. assert ( @@ -363,8 +359,8 @@ def test_postprocessors(petab_select_problem): ) # End Iteration 1 ######################################################### - expected_png_file = output_path / (best_model_1.model_hash + ".png") - expected_hdf5_file = output_path / (best_model_1.model_hash + ".hdf5") + expected_png_file = output_path / f"{best_model_1.hash}.png" + expected_hdf5_file = output_path / f"{best_model_1.hash}.hdf5" # The expected files exist. assert expected_png_file.is_file() @@ -400,33 +396,6 @@ def test_model_problem_fake_result(): assert test_fval == expected_fval -def test_vis(pypesto_select_problem): - """Test plotting routines.""" - criterion = Criterion.AIC - best_models = pypesto_select_problem.select_to_completion( - method=Method.FORWARD, - criterion=criterion, - model_problem_options=model_problem_options, - ) - labels = { - model.get_hash(): model.model_subspace_id - for model in pypesto_select_problem.calibrated_models.values() - } - pypesto.visualize.select.plot_selected_models( - selected_models=best_models, - criterion=criterion, - labels=labels, - ) - # import matplotlib.pyplot as plt - # plt.savefig('output/selected.png') - pypesto.visualize.select.plot_calibrated_models_digraph( - problem=pypesto_select_problem, - criterion=criterion, - labels=labels, - ) - # plt.savefig('output/digraph.png') - - def test_custom_objective(petab_problem_yaml): parameters = { "k2": 0.333, @@ -445,7 +414,12 @@ def fun_grad(x, fun=expected_fun, grad=expected_grad): grad=True, ) - model = Model(petab_yaml=petab_problem_yaml) + model = Model( + model_subspace_id="dummy", + model_subspace_indices=[], + model_subspace_petab_yaml=petab_problem_yaml, + parameters={}, + ) def get_pypesto_problem(model, objective, petab_problem, x_guesses): corrected_x_guesses = correct_x_guesses( @@ -510,13 +484,8 @@ def test_sacess_minimize_method(pypesto_select_problem, initial_models): max_walltime_s=1, ) - minimize_options = { - "startpoint_method": pypesto.startpoint.UniformStartpoints(), - } - model_problem_options = { "minimize_method": minimize_method, - "minimize_options": minimize_options, } pypesto_select_problem.select_to_completion( diff --git a/test/select/test_test_cases.py b/test/select/test_test_cases.py new file mode 100644 index 000000000..8bafdfd87 --- /dev/null +++ b/test/select/test_test_cases.py @@ -0,0 +1,218 @@ +import os +import shlex +import shutil +import subprocess +from pathlib import Path + +import numpy as np +import pandas as pd +import petab_select +import pytest +import yaml +from petab_select import Model +from petab_select.constants import ( + CRITERIA, + ESTIMATED_PARAMETERS, + TERMINATE, +) + +import pypesto.engine +import pypesto.optimize +import pypesto.select + +# Set to `[]` to test all +test_cases = [ + # "0009", +] + +# Do not use computationally-expensive test cases in CI +skip_test_cases = [ + "0009", +] + +# Download test cases from GitHub +test_cases_path = Path("petab_select") / "test_cases" +if not test_cases_path.exists(): + subprocess.run([Path(__file__).parent / "get_test_cases.sh"]) # noqa: S603 + +# Reduce runtime but with high reproducibility +minimize_options = { + "n_starts": 10, + "engine": pypesto.engine.MultiProcessEngine(), + "filename": None, + "progress_bar": False, +} + + +def objective_customizer(obj): + obj.amici_solver.setRelativeTolerance(1e-12) + + +model_problem_options = { + "minimize_options": minimize_options, + "objective_customizer": objective_customizer, +} + + +@pytest.mark.parametrize( + "test_case_path_stem", + sorted( + [test_case_path.stem for test_case_path in test_cases_path.glob("*")] + ), +) +def test_pypesto(test_case_path_stem): + """Run all test cases with pyPESTO.""" + if test_cases and test_case_path_stem not in test_cases: + pytest.skip("Test excluded from subset selected for debugging.") + return + + if test_case_path_stem in skip_test_cases: + pytest.skip("Test marked to be skipped.") + return + + test_case_path = test_cases_path / test_case_path_stem + expected_model_yaml = test_case_path / "expected.yaml" + # Setup the pyPESTO model selector instance. + petab_select_problem = petab_select.Problem.from_yaml( + test_case_path / "petab_select_problem.yaml", + ) + pypesto_select_problem = pypesto.select.Problem( + petab_select_problem=petab_select_problem + ) + + # Run the selection process until "exhausted". + pypesto_select_problem.select_to_completion( + model_problem_options=model_problem_options, + ) + + # Get the best model + best_model = petab_select.analyze.get_best( + models=pypesto_select_problem.calibrated_models, + criterion=petab_select_problem.criterion, + compare=petab_select_problem.compare, + ) + + # Load the expected model. + expected_model = Model.from_yaml(expected_model_yaml) + + def get_series(model, dict_attribute) -> pd.Series: + return pd.Series( + getattr(model, dict_attribute), + dtype=np.float64, + ).sort_index() + + # The estimated parameters and criteria values are as expected. + for dict_attribute in [CRITERIA, ESTIMATED_PARAMETERS]: + pd.testing.assert_series_equal( + get_series(expected_model, dict_attribute), + get_series(best_model, dict_attribute), + rtol=1e-2, + ) + # FIXME ensure `current model criterion` trajectory also matches, in summary.tsv file, + # for test case 0009, after summary format is revised + + +@pytest.mark.skipif( + os.getenv("GITHUB_ACTIONS") == "true", + reason="Too CPU heavy for CI.", +) +def test_famos_cli(): + """Run test case 0009 with pyPESTO and the CLI interface.""" + test_case_path = test_cases_path / "0009" + expected_model_yaml = test_case_path / "expected.yaml" + problem_yaml = test_case_path / "petab_select_problem.yaml" + + problem = petab_select.Problem.from_yaml(problem_yaml) + + # Setup working directory for intermediate files + work_dir = Path("output_famos_cli") + work_dir_str = str(work_dir) + if work_dir.exists(): + shutil.rmtree(work_dir_str) + work_dir.mkdir(exist_ok=True, parents=True) + + models_yamls = [] + metadata_yaml = work_dir / "metadata.yaml" + state_dill = work_dir / "state.dill" + iteration = 0 + while True: + iteration += 1 + uncalibrated_models_yaml = ( + work_dir / f"uncalibrated_models_{iteration}.yaml" + ) + calibrated_models_yaml = ( + work_dir / f"calibrated_models_{iteration}.yaml" + ) + models_yaml = work_dir / f"models_{iteration}.yaml" + models_yamls.append(models_yaml) + # Start iteration + subprocess.run( + shlex.split( # noqa: S603 + f"""petab_select start_iteration + --problem {problem_yaml} + --state {state_dill} + --output-uncalibrated-models {uncalibrated_models_yaml} + --relative-paths + """ + ) + ) + # Calibrate models + models = petab_select.Models.from_yaml(uncalibrated_models_yaml) + for model in models: + pypesto.select.ModelProblem( + model=model, + criterion=problem.criterion, + **model_problem_options, + ) + models.to_yaml(filename=calibrated_models_yaml) + # End iteration + subprocess.run( + shlex.split( # noqa: S603 + f"""petab_select end_iteration + --output-models {models_yaml} + --output-metadata {metadata_yaml} + --state {state_dill} + --calibrated-models {calibrated_models_yaml} + --relative-paths + """ + ) + ) + with open(metadata_yaml) as f: + metadata = yaml.safe_load(f) + if metadata[TERMINATE]: + break + + # Get the best model + models_yamls_arg = " ".join( + f"--models {models_yaml}" for models_yaml in models_yamls + ) + subprocess.run( + shlex.split( # noqa: S603 + f"""petab_select get_best + --problem {problem_yaml} + {models_yamls_arg} + --output {work_dir / "best_model.yaml"} + --relative-paths + """ + ) + ) + best_model = petab_select.Model.from_yaml(work_dir / "best_model.yaml") + + # Load the expected model. + expected_model = Model.from_yaml(expected_model_yaml) + + def get_series(model, dict_attribute) -> pd.Series: + return pd.Series( + getattr(model, dict_attribute), + dtype=np.float64, + ).sort_index() + + # The estimated parameters and criteria values are as expected. + for dict_attribute in [CRITERIA, ESTIMATED_PARAMETERS]: + pd.testing.assert_series_equal( + get_series(expected_model, dict_attribute), + get_series(best_model, dict_attribute), + rtol=1e-2, + ) + # FIXME ensure `current model criterion` trajectory also matches, in summary.tsv file, + # after summary format is revised diff --git a/tox.ini b/tox.ini index d1937e4ee..98f96b7bf 100644 --- a/tox.ini +++ b/tox.ini @@ -29,7 +29,7 @@ envlist = # Base-environment [testenv] -passenv = AMICI_PARALLEL_COMPILE,CC,CXX,MPLBACKEND,BNGPATH +passenv = AMICI_PARALLEL_COMPILE,BNGPATH,CC,CXX,GITHUB_ACTIONS,MPLBACKEND # Sub-environments # inherit settings defined in the base @@ -126,7 +126,7 @@ description = [testenv:select] extras = test,amici,petab,select,fides commands = - pytest --cov=pypesto --cov-report=xml --cov-append \ + pytest -svvv --cov=pypesto --cov-report=xml --cov-append \ test/select description = Test model selection functionality From c0339bb982192c5217f2edabd7ba22bb4e63e585 Mon Sep 17 00:00:00 2001 From: Paul Jonas Jost <70631928+PaulJonasJost@users.noreply.github.com> Date: Mon, 13 Jan 2025 14:50:53 +0100 Subject: [PATCH 35/35] Prepare release 0.5.5 (#1546) * Updated Changelog and Version * Correct intendation * Added #1544 to changelog preemtively as it needs to be merged before release anyways * Updated recent pushs * Updated changelod for recent pushes * Forgot to mention one PR with # * Correct intendation --- CHANGELOG.rst | 50 ++++++++++++++++++++++++++++++++++++++++++++++ pypesto/version.py | 2 +- 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index c56d4b29e..a46a808e0 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -6,6 +6,56 @@ Release notes .......... +0.5.5 (2025-01-10) +------------------- + +- **Breaking Changes** + - **PETab select**: There are some deprecated features that will show up as warnings. In addition: + + - The plotting methods ignore some arguments. You will need to reimplement these with the newer approach, which uses + plotting methods from the PEtab Select library instead. See the model selection notebook for examples. + + - All objects containing multiple models (e.g., dictionaries or lists) are now replaced by `petab_select.Models`, + which supports dictionary and list methods. + + To convert your old list of models: + + ```python + petab_select.Models(list_of_Model) + ``` +- General + - Exclude nlopt==2.9.0 from setup (#1519) + - Improve CI (#1521, #1523, #1532, #1536, #1508, #1544, #1531) + - Update references/documentation (#1506, #1491, #1516, #1543) + - **Docker Image** (#1083, #1538) +- Hierarchical + - Fix no error if inner observable parameter in noise formula & viceversa (#1504) + - Remove inner datas from relative calculator (#1505) + - Fix not scaling inner pars when applying to rdatas (#1534) +- Optimization + - ESSOptimizer: Fix priority for local search startpoints (#1503) + - Fix NLoptOptimizer.__repr__ (#1518) + - Improve exception-handling in SacessOptimizer (#1517) + - Fix ESSOptimizer min of empty sequence (#1510) + - Don't modify sys.path for amici model imports (#1522) + - Set OptimizerResult.optimizer in Optimizer.minimize (#1525) + - SacessOptimizer: More efficient saving of intermediate results (#1529) +- Objective + - AmiciObjective/PEtab import: Fix plist (#1493) + - PEtab: Fix warning from fill_in_parameters with fixed parameters (#1509) + - Amici: Fix handling of PEtab fixed parameters (#1514) + - Fix get_parameter_prior_dict docstring (#1537) +- Select + - Support user-provided calibration results (#1338) + - Problem-specific minimize method for SaCeSS (#1339) + - Update for the latest PEtab Select version; see example notebook or the PEtab Select repo (#1530) +- Storage + - Enable writing Optimize(r)Result directly in Writer (#1528) + - Update parameter scale storage (#1542 +- Visualize + - Fix flatten of observable mapping with one observable (#1515) + + 0.5.4 (2024-10-19) ------------------- diff --git a/pypesto/version.py b/pypesto/version.py index 6b27eeebf..86716a713 100644 --- a/pypesto/version.py +++ b/pypesto/version.py @@ -1 +1 @@ -__version__ = "0.5.4" +__version__ = "0.5.5"