diff --git a/pypesto/C.py b/pypesto/C.py index 7b8ff0fd2..25f1468cb 100644 --- a/pypesto/C.py +++ b/pypesto/C.py @@ -232,6 +232,12 @@ class InnerParameterType(str, Enum): SUFFIXES_HDF5 = ["hdf5", "h5"] SUFFIXES = SUFFIXES_CSV + SUFFIXES_HDF5 +CPU_TIME_TOTAL = 'cpu_time_total' +PREEQ_CPU_TIME = 'preeq_cpu_time' +PREEQ_CPU_TIME_BACKWARD = 'preeq_cpu_timeB' +POSTEQ_CPU_TIME = 'posteq_cpu_time' +POSTEQ_CPU_TIME_BACKWARD = 'posteq_cpu_timeB' + ############################################################################### # PRIOR diff --git a/pypesto/__init__.py b/pypesto/__init__.py index d1f31749f..7bcd97165 100644 --- a/pypesto/__init__.py +++ b/pypesto/__init__.py @@ -16,7 +16,9 @@ CountHistory, CountHistoryBase, CsvHistory, + CsvAmiciHistory, Hdf5History, + Hdf5AmiciHistory, NoHistory, HistoryBase, HistoryOptions, diff --git a/pypesto/hierarchical/base_problem.py b/pypesto/hierarchical/base_problem.py index 28279f65f..5fa58f725 100644 --- a/pypesto/hierarchical/base_problem.py +++ b/pypesto/hierarchical/base_problem.py @@ -74,9 +74,9 @@ def get_interpretable_x_ids(self) -> list[str]: Interpretable parameters need to be easily interpretable by the user. Examples are scaling factors, offsets, or noise parameters. An example - for a non-interpretable inner parameters are spline heights of spline - approximation for semiquantitative data: it is hard to interpret what - the spline heights are just by looking at the parameter value. + of non-interpretable inner parameters is the spline heights of spline + approximation for semiquantitative data. It is challenging to interpret + the meaning of these parameters based solely on their value. """ return list(self.xs.keys()) diff --git a/pypesto/hierarchical/inner_calculator_collector.py b/pypesto/hierarchical/inner_calculator_collector.py index f5e4c794c..843dd0ae3 100644 --- a/pypesto/hierarchical/inner_calculator_collector.py +++ b/pypesto/hierarchical/inner_calculator_collector.py @@ -248,7 +248,10 @@ def get_inner_par_ids(self) -> list[str]: ] def get_interpretable_inner_par_ids(self) -> list[str]: - """Return the ids of interpretable inner parameters of all inner problems.""" + """Return the ids of interpretable inner parameters of all inner problems. + + See :func:`InnerProblem.get_interpretable_x_ids`. + """ return [ parameter_id for inner_calculator in self.inner_calculators diff --git a/pypesto/history/__init__.py b/pypesto/history/__init__.py index 574418a48..37fc18377 100644 --- a/pypesto/history/__init__.py +++ b/pypesto/history/__init__.py @@ -8,6 +8,7 @@ and evaluate performance. """ +from .amici import CsvAmiciHistory, Hdf5AmiciHistory from .base import CountHistory, CountHistoryBase, HistoryBase, NoHistory from .csv import CsvHistory from .generate import create_history diff --git a/pypesto/history/amici.py b/pypesto/history/amici.py new file mode 100644 index 000000000..cfd330cd5 --- /dev/null +++ b/pypesto/history/amici.py @@ -0,0 +1,236 @@ +from pathlib import Path +from typing import Sequence, Union + +import numpy as np + +from ..C import ( + CPU_TIME_TOTAL, + POSTEQ_CPU_TIME, + POSTEQ_CPU_TIME_BACKWARD, + PREEQ_CPU_TIME, + PREEQ_CPU_TIME_BACKWARD, + RDATAS, +) +from .csv import CsvHistory +from .hdf5 import Hdf5History +from .options import HistoryOptions +from .util import trace_wrap + + +class Hdf5AmiciHistory(Hdf5History): + """ + Stores history extended by AMICI-specific time traces in an HDF5 file. + + Stores AMICI-specific traces of total simulation time, pre-equilibration + time and post-equilibration time. + + Parameters + ---------- + id: + Id of the history + file: + HDF5 file name. + options: + History options. Defaults to ``None``. + """ + + def __init__( + self, + id: str, + file: Union[str, Path], + options: Union[HistoryOptions, dict, None] = None, + ): + super().__init__(id, file, options=options) + + @staticmethod + def _simulation_to_values(x, result, used_time): + values = Hdf5History._simulation_to_values(x, result, used_time) + # default unit for time in amici is [ms], converted to [s] + values |= { + key: sum([rdata[key] for rdata in result[RDATAS]]) * 0.001 + for key in ( + CPU_TIME_TOTAL, + PREEQ_CPU_TIME, + PREEQ_CPU_TIME_BACKWARD, + POSTEQ_CPU_TIME, + POSTEQ_CPU_TIME_BACKWARD, + ) + } + return values + + @trace_wrap + def get_cpu_time_total_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative simulation CPU time [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return self._get_hdf5_entries(CPU_TIME_TOTAL, ix) + + @trace_wrap + def get_preeq_time_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative pre-equilibration time, [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return self._get_hdf5_entries(PREEQ_CPU_TIME, ix) + + @trace_wrap + def get_preeq_timeB_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative pre-equilibration time of the backward problem, [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return self._get_hdf5_entries(PREEQ_CPU_TIME_BACKWARD, ix) + + @trace_wrap + def get_posteq_time_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative post-equilibration time [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return self._get_hdf5_entries(POSTEQ_CPU_TIME, ix) + + @trace_wrap + def get_posteq_timeB_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative post-equilibration time of the backward problem [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return self._get_hdf5_entries(POSTEQ_CPU_TIME_BACKWARD, ix) + + +class CsvAmiciHistory(CsvHistory): + """ + Stores history extended by AMICI-specific time traces in a CSV file. + + Stores AMICI-specific traces of total simulation time, pre-equilibration + time and post-equilibration time. + + Parameters + ---------- + file: + CSV file name. + x_names: + Parameter names. + options: + History options. + load_from_file: + If True, history will be initialized from data in the specified file. + """ + + def __init__( + self, + file: str, + x_names: Sequence[str] = None, + options: Union[HistoryOptions, dict] = None, + load_from_file: bool = False, + ): + super().__init__(file, x_names, options, load_from_file=load_from_file) + + def _trace_columns(self) -> list[tuple]: + columns = super()._trace_columns() + return columns + [ + (c, np.nan) + for c in [ + CPU_TIME_TOTAL, + PREEQ_CPU_TIME, + PREEQ_CPU_TIME_BACKWARD, + POSTEQ_CPU_TIME, + POSTEQ_CPU_TIME_BACKWARD, + ] + ] + + def _simulation_to_values(self, result, used_time): + values = super()._simulation_to_values(result, used_time) + # default unit for time in amici is [ms], converted to [s] + values |= { + key: sum([rdata[key] for rdata in result[RDATAS]]) * 0.001 + for key in ( + CPU_TIME_TOTAL, + PREEQ_CPU_TIME, + PREEQ_CPU_TIME_BACKWARD, + POSTEQ_CPU_TIME, + POSTEQ_CPU_TIME_BACKWARD, + ) + } + return values + + @trace_wrap + def get_cpu_time_total_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative simulation CPU time [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return list(self._trace[(CPU_TIME_TOTAL, np.nan)].values[ix]) + + @trace_wrap + def get_preeq_time_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative pre-equilibration time [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return list(self._trace[(PREEQ_CPU_TIME, np.nan)].values[ix]) + + @trace_wrap + def get_preeq_timeB_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative pre-equilibration time of the backward problem [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return list(self._trace[(PREEQ_CPU_TIME_BACKWARD, np.nan)].values[ix]) + + @trace_wrap + def get_posteq_time_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative post-equilibration time [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return list(self._trace[(POSTEQ_CPU_TIME, np.nan)].values[ix]) + + @trace_wrap + def get_posteq_timeB_trace( + self, ix: Union[int, Sequence[int], None] = None, trim: bool = False + ) -> Union[Sequence[float], float]: + """ + Cumulative post-equilibration time of the backward problem [s]. + + Takes as parameter an index or indices and returns corresponding trace + values. If only a single value is requested, the list is flattened. + """ + return list(self._trace[(POSTEQ_CPU_TIME_BACKWARD, np.nan)].values[ix]) diff --git a/pypesto/history/base.py b/pypesto/history/base.py index 27f5ef455..8870830ff 100644 --- a/pypesto/history/base.py +++ b/pypesto/history/base.py @@ -249,7 +249,7 @@ def get_time_trace( trim: bool = False, ) -> Union[Sequence[float], float]: """ - Cumulative execution times. + Cumulative execution times [s]. Takes as parameter an index or indices and returns corresponding trace values. If only a single value is requested, the list is flattened. diff --git a/pypesto/history/csv.py b/pypesto/history/csv.py index 07df30fbf..c672fe8fe 100644 --- a/pypesto/history/csv.py +++ b/pypesto/history/csv.py @@ -100,6 +100,20 @@ def finalize(self, message: str = None, exitflag: str = None): super().finalize(message=message, exitflag=exitflag) self._save_trace(finalize=True) + def _simulation_to_values(self, result, used_time): + values = { + TIME: used_time, + N_FVAL: self._n_fval, + N_GRAD: self._n_grad, + N_HESS: self._n_hess, + N_RES: self._n_res, + N_SRES: self._n_sres, + FVAL: result[FVAL], + RES: result[RES], + HESS: result[HESS], + } + return values + def _update_trace( self, x: np.ndarray, @@ -127,17 +141,7 @@ def _update_trace( name=len(self._trace), index=self._trace.columns, dtype='object' ) - values = { - TIME: used_time, - N_FVAL: self._n_fval, - N_GRAD: self._n_grad, - N_HESS: self._n_hess, - N_RES: self._n_res, - N_SRES: self._n_sres, - FVAL: result[FVAL], - RES: result[RES], - HESS: result[HESS], - } + values = self._simulation_to_values(result, used_time) for var, val in values.items(): row[(var, np.nan)] = val @@ -162,12 +166,8 @@ def _update_trace( # save trace to file self._save_trace() - def _init_trace(self, x: np.ndarray): - """Initialize the trace.""" - if self.x_names is None: - self.x_names = [f'x{i}' for i, _ in enumerate(x)] - - columns: list[tuple] = [ + def _trace_columns(self) -> list[tuple]: + return [ (c, np.nan) for c in [ TIME, @@ -183,6 +183,13 @@ def _init_trace(self, x: np.ndarray): ] ] + def _init_trace(self, x: np.ndarray): + """Initialize the trace.""" + if self.x_names is None: + self.x_names = [f'x{i}' for i, _ in enumerate(x)] + + columns = self._trace_columns() + for var in [X, GRAD]: if var == X or self.options[f'trace_record_{var}']: columns.extend([(var, x_name) for x_name in self.x_names]) diff --git a/pypesto/history/generate.py b/pypesto/history/generate.py index 5fdbbfe4c..3e291c7c7 100644 --- a/pypesto/history/generate.py +++ b/pypesto/history/generate.py @@ -13,9 +13,7 @@ def create_history( - id: str, - x_names: Sequence[str], - options: HistoryOptions, + id: str, x_names: Sequence[str], options: HistoryOptions ) -> HistoryBase: """Create a :class:`HistoryBase` object; Factory method. diff --git a/pypesto/history/hdf5.py b/pypesto/history/hdf5.py index a56366ea0..4d05cbd97 100644 --- a/pypesto/history/hdf5.py +++ b/pypesto/history/hdf5.py @@ -302,6 +302,19 @@ def exitflag(self) -> str: except KeyError: return None + @staticmethod + def _simulation_to_values(x, result, used_time): + values = { + X: x, + FVAL: result[FVAL], + GRAD: result[GRAD], + RES: result[RES], + SRES: result[SRES], + HESS: result[HESS], + TIME: used_time, + } + return values + @with_h5_file("a") def _update_trace( self, @@ -322,15 +335,7 @@ def _update_trace( used_time = time.time() - self.start_time - values = { - X: x, - FVAL: result[FVAL], - GRAD: result[GRAD], - RES: result[RES], - SRES: result[SRES], - HESS: result[HESS], - TIME: used_time, - } + values = self._simulation_to_values(x, result, used_time) iteration = self._require_group().attrs[N_ITERATIONS] diff --git a/pypesto/objective/amici/amici.py b/pypesto/objective/amici/amici.py index 521dba3b8..13e2d87e8 100644 --- a/pypesto/objective/amici/amici.py +++ b/pypesto/objective/amici/amici.py @@ -3,11 +3,29 @@ import os import tempfile from collections import OrderedDict +from pathlib import Path from typing import TYPE_CHECKING, Dict, Optional, Sequence, Tuple, Union import numpy as np -from ...C import FVAL, INNER_PARAMETERS, MODE_FUN, MODE_RES, RDATAS, ModeType +from ...C import ( + FVAL, + INNER_PARAMETERS, + MODE_FUN, + MODE_RES, + RDATAS, + SUFFIXES_CSV, + SUFFIXES_HDF5, + ModeType, +) +from ...history import ( + CountHistory, + CsvAmiciHistory, + Hdf5AmiciHistory, + HistoryOptions, + HistoryTypeError, + MemoryHistory, +) from ..base import ObjectiveBase, ResultDict from .amici_calculator import AmiciCalculator from .amici_util import ( @@ -227,6 +245,33 @@ def get_config(self) -> dict: return info + def create_history( + self, id: str, x_names: Sequence[str], options: HistoryOptions + ): + """See `history.generate.create_history` documentation.""" + # create different history types based on the inputs + if options.storage_file is None: + if options.trace_record: + return MemoryHistory(options=options) + else: + return CountHistory(options=options) + + # replace id template in storage file + storage_file = options.storage_file.replace("{id}", id) + + # evaluate type + suffix = Path(storage_file).suffix[1:] + + # create history type based on storage type + if suffix in SUFFIXES_CSV: + return CsvAmiciHistory( + x_names=x_names, file=storage_file, options=options + ) + elif suffix in SUFFIXES_HDF5: + return Hdf5AmiciHistory(id=id, file=storage_file, options=options) + else: + raise HistoryTypeError(suffix) + def initialize(self): """See `ObjectiveBase` documentation.""" super().initialize() diff --git a/pypesto/objective/base.py b/pypesto/objective/base.py index 3e47f0ef7..a38114bd4 100644 --- a/pypesto/objective/base.py +++ b/pypesto/objective/base.py @@ -7,7 +7,7 @@ import pandas as pd from ..C import FVAL, GRAD, HESS, MODE_FUN, MODE_RES, RES, SRES, ModeType -from ..history import NoHistory +from ..history import NoHistory, create_history from .pre_post_process import FixedParametersProcessor, PrePostProcessor ResultDict = Dict[str, Union[float, np.ndarray, Dict]] @@ -119,6 +119,10 @@ def initialize(self): By default does nothing. """ + def create_history(self, id, x_names, options): + """See `history.generate.create_history` documentation.""" + return create_history(id, x_names, options) + def __call__( self, x: np.ndarray, diff --git a/pypesto/optimize/ess/ess.py b/pypesto/optimize/ess/ess.py index 80f25cad0..68de3f4d6 100644 --- a/pypesto/optimize/ess/ess.py +++ b/pypesto/optimize/ess/ess.py @@ -6,7 +6,7 @@ import enum import logging import time -from typing import Callable, List, Optional, Tuple, Union +from typing import Callable, Optional, Union from warnings import warn import numpy as np @@ -65,6 +65,7 @@ def __init__( n_procs=None, n_threads=None, max_walltime_s=None, + result_includes_refset: bool = False, ): """Construct new ESS instance. @@ -113,6 +114,10 @@ def __init__( history: History of the best values/parameters found so far. (Monotonously decreasing objective values.) + result_includes_refset: + Whether the :meth:`minimize` result should include the final + RefSet, or just the local search results and the overall best + parameters. """ if max_eval is None and max_walltime_s is None and max_iter is None: # in this case, we'd run forever @@ -151,6 +156,7 @@ def __init__( self.logger = logging.getLogger( f"{self.__class__.__name__}-{id(self)}" ) + self._result_includes_refset = result_includes_refset def _initialize(self): """(Re-)Initialize.""" @@ -160,8 +166,8 @@ def _initialize(self): self.x_best: Optional[np.array] = None # Overall best function value found so far self.fx_best: float = np.inf - # Final parameters from local searches - self.local_solutions: List[np.array] = [] + # Results from local searches + self.local_solutions: list[OptimizerResult] = [] # Index of current iteration self.n_iter: int = 0 # ESS iteration at which the last local search took place @@ -313,25 +319,29 @@ def _create_result(self) -> pypesto.Result: **common_result_fields, ) optimizer_result.update_to_full(result.problem) - # TODO DW: Create a single History with the global best? result.optimize_result.append(optimizer_result) - # save refset - for i in range(self.refset.dim): + # save local solutions + for i, optimizer_result in enumerate(self.local_solutions): i_result += 1 - result.optimize_result.append( - pypesto.OptimizerResult( - id=str(i_result), - x=self.refset.x[i], - fval=self.refset.fx[i], - message=f"RefSet[{i}]", - **common_result_fields, + 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: + # save refset + for i in range(self.refset.dim): + i_result += 1 + result.optimize_result.append( + pypesto.OptimizerResult( + id=str(i_result), + x=self.refset.x[i], + fval=self.refset.fx[i], + message=f"RefSet[{i}]", + **common_result_fields, + ) ) - ) - result.optimize_result[-1].update_to_full(result.problem) - - # TODO DW: also save local solutions? - # (need to track fvals or re-evaluate) + result.optimize_result[-1].update_to_full(result.problem) return result @@ -370,20 +380,20 @@ 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.array, np.array]: """Combine solutions and evaluate. - Creates the next generation from the RefSet by pair-wise combinations + Creates the next generation from the RefSet by pair-wise combination of all RefSet members. Creates ``RefSet.dim ** 2 - RefSet.dim`` new parameter vectors, tests them, and keeps the best child of each parent. Returns ------- y: - Contains the best of all pairwise combinations of all RefSet - members, for each RefSet members. + The next generation of parameter vectors + (`dim_refset` x `dim_problem`). fy: - The corresponding objective values. + The objective values corresponding to the parameters in `y`. """ y = np.zeros(shape=(self.refset.dim, self.evaluator.problem.dim)) fy = np.full(shape=self.refset.dim, fill_value=np.inf) @@ -471,8 +481,10 @@ def _do_local_search( # optima found so far min_distances = np.array( np.min( - np.linalg.norm(y_i - local_solution) - for local_solution in self.local_solutions + np.linalg.norm( + y_i - optimizer_result.x[optimizer_result.free_indices] + ) + for optimizer_result in self.local_solutions ) for y_i in x_best_children ) @@ -521,15 +533,11 @@ def _do_local_search( f"{optimizer_result.exitflag}: {optimizer_result.message}" ) if np.isfinite(optimizer_result.fval): - local_solution_x = optimizer_result.x[ - optimizer_result.free_indices - ] - local_solution_fx = optimizer_result.fval - - self.local_solutions.append(local_solution_x) + self.local_solutions.append(optimizer_result) self._maybe_update_global_best( - local_solution_x, local_solution_fx + optimizer_result.x[optimizer_result.free_indices], + optimizer_result.fval, ) break diff --git a/pypesto/optimize/ess/sacess.py b/pypesto/optimize/ess/sacess.py index 7879821af..195087b7f 100644 --- a/pypesto/optimize/ess/sacess.py +++ b/pypesto/optimize/ess/sacess.py @@ -543,7 +543,6 @@ def run( ) ) - ess_results = pypesto.Result(problem=problem) ess = self._setup_ess(startpoint_method) # run ESS until exit criteria are met, but start at least one iteration @@ -551,17 +550,9 @@ def run( # perform one ESS iteration ess._do_iteration() - # TODO maybe not in every iteration? - # drop all but the 50 best results - cur_ess_results = ess._create_result() - ess_results.optimize_result.append( - cur_ess_results.optimize_result, - prefix=f"{self._worker_idx}_{ess.n_iter}_", - ) - ess_results.optimize_result.list = ( - ess_results.optimize_result.list[:50] - ) if self._tmp_result_file: + # TODO maybe not in every iteration? + ess_results = ess._create_result() write_result( ess_results, self._tmp_result_file, @@ -580,6 +571,7 @@ def run( f"sacess worker {self._worker_idx} iteration {ess.n_iter} " f"(best: {self._best_known_fx})." ) + ess.history.finalize(exitflag=ess.exit_flag.name) self._manager._result_queue.put(ess.history) ess._report_final() diff --git a/pypesto/optimize/optimizer.py b/pypesto/optimize/optimizer.py index 52e83a269..166020958 100644 --- a/pypesto/optimize/optimizer.py +++ b/pypesto/optimize/optimizer.py @@ -12,12 +12,7 @@ import scipy.optimize from ..C import FVAL, GRAD, INNER_PARAMETERS, MODE_FUN, MODE_RES -from ..history import ( - HistoryOptions, - NoHistory, - OptimizerHistory, - create_history, -) +from ..history import HistoryOptions, NoHistory, OptimizerHistory from ..objective import Objective from ..problem import Problem from ..result import OptimizerResult @@ -102,7 +97,7 @@ def wrapped_minimize( objective.initialize() # initialize the history - history = create_history( + history = objective.create_history( id=id, x_names=[problem.x_names[ix] for ix in problem.x_free_indices], options=history_options, diff --git a/pypesto/petab/importer.py b/pypesto/petab/importer.py index bb04ffc23..be4a32aee 100644 --- a/pypesto/petab/importer.py +++ b/pypesto/petab/importer.py @@ -10,6 +10,7 @@ import warnings from dataclasses import dataclass from functools import partial +from importlib.metadata import version from typing import ( Any, Callable, @@ -45,7 +46,7 @@ from ..objective.amici import AmiciObjectBuilder from ..objective.priors import NegLogParameterPriors, get_parameter_prior_dict from ..predict import AmiciPredictor -from ..problem import Problem +from ..problem import HierarchicalProblem, Problem from ..result import PredictionResult from ..startpoint import CheckedStartpoints, StartpointMethod @@ -64,7 +65,7 @@ ) from petab.models import MODEL_TYPE_SBML except ImportError: - pass + amici = None logger = logging.getLogger(__name__) @@ -81,7 +82,7 @@ class PetabImporter(AmiciObjectBuilder): `PEtab documentation `_. """ # noqa - MODEL_BASE_DIR = "amici_models" + MODEL_BASE_DIR = f"amici_models/{version('amici') if amici else ''}" def __init__( self, @@ -780,7 +781,12 @@ def create_problem( ) objective = AggregatedObjective([objective, prior]) - problem = Problem( + if self._hierarchical: + problem_class = HierarchicalProblem + else: + problem_class = Problem + + problem = problem_class( objective=objective, lb=lb, ub=ub, diff --git a/pypesto/problem/__init__.py b/pypesto/problem/__init__.py index 63512a23a..1fc5f7ca6 100644 --- a/pypesto/problem/__init__.py +++ b/pypesto/problem/__init__.py @@ -8,3 +8,4 @@ """ from .base import Problem +from .hierarchical import HierarchicalProblem diff --git a/pypesto/problem/hierarchical.py b/pypesto/problem/hierarchical.py new file mode 100644 index 000000000..26061d691 --- /dev/null +++ b/pypesto/problem/hierarchical.py @@ -0,0 +1,72 @@ +import logging +from typing import Iterable, List, Optional, SupportsFloat, SupportsInt, Union + +import numpy as np + +from .base import Problem + +SupportsFloatIterableOrValue = Union[Iterable[SupportsFloat], SupportsFloat] +SupportsIntIterableOrValue = Union[Iterable[SupportsInt], SupportsInt] + +logger = logging.getLogger(__name__) + + +class HierarchicalProblem(Problem): + """ + The Hierarchical Problem. + + A hierarchical problem is a problem with a nested structure: One or + multiple inner problems are nested inside the outer problem. The inner + problems are optimized for each evaluation of the outer problem. The + objective's calculator is used to collect the inner problems' objective + values. + + Parameters + ---------- + hierarchical: + A flag indicating the problem is hierarchical. + inner_x_names: + Names of the inner optimization parameters. Only relevant if + hierarchical is True. Contains the names of easily interpretable + inner parameters only, e.g. noise parameters, scaling factors, offsets. + inner_lb, inner_ub: + The lower and upper bounds for the inner optimization parameters. + Only relevant if hierarchical is True. Contains the bounds of easily + interpretable inner parameters only, e.g. noise parameters, scaling + factors, offsets. + """ + + def __init__( + self, + inner_x_names: Optional[Iterable[str]] = None, + inner_lb: Optional[Union[np.ndarray, List[float]]] = None, + inner_ub: Optional[Union[np.ndarray, List[float]]] = None, + **problem_kwargs: dict, + ): + super().__init__(**problem_kwargs) + + if inner_x_names is None: + inner_x_names = ( + self.objective.calculator.get_interpretable_inner_par_ids() + ) + if len(set(inner_x_names)) != len(inner_x_names): + raise ValueError("Parameter names inner_x_names must be unique") + self.inner_x_names = inner_x_names + + if inner_lb is None or inner_ub is None: + ( + default_inner_lb, + default_inner_ub, + ) = self.objective.calculator.get_interpretable_inner_par_bounds() + inner_lb = default_inner_lb if inner_lb is None else inner_lb + inner_ub = default_inner_ub if inner_ub is None else inner_ub + + if len(inner_lb) != len(inner_ub): + raise ValueError("Parameter bounds must have same length") + if len(inner_lb) != len(inner_x_names): + raise ValueError( + "Parameter bounds must have same length as parameter names" + ) + + self.inner_lb = np.array(inner_lb) + self.inner_ub = np.array(inner_ub) diff --git a/pypesto/visualize/parameters.py b/pypesto/visualize/parameters.py index 9d0088d43..2e5ae162d 100644 --- a/pypesto/visualize/parameters.py +++ b/pypesto/visualize/parameters.py @@ -462,26 +462,22 @@ def _handle_inner_inputs( inner_lb = None inner_ub = None - if any(inner_x is not None for inner_x in inner_xs): - from ..hierarchical import InnerCalculatorCollector - - if hasattr(result.problem.objective, 'calculator') and isinstance( - inner_calculator := result.problem.objective.calculator, - InnerCalculatorCollector, - ): - inner_xs_names = inner_calculator.get_interpretable_inner_par_ids() - # replace None with a list of nans - inner_xs = [ - np.full(len(inner_xs_names), np.nan) - if inner_xs_idx is None - else np.asarray(inner_xs_idx) - for inner_xs_idx in inner_xs - ] - # set bounds for inner parameters - ( - inner_lb, - inner_ub, - ) = inner_calculator.get_interpretable_inner_par_bounds() + from ..problem import HierarchicalProblem + + if any(inner_x is not None for inner_x in inner_xs) and isinstance( + result.problem, HierarchicalProblem + ): + inner_xs_names = result.problem.inner_x_names + # replace None with a list of nans + inner_xs = [ + np.full(len(inner_xs_names), np.nan) + if inner_xs_idx is None + else np.asarray(inner_xs_idx) + for inner_xs_idx in inner_xs + ] + # set bounds for inner parameters + inner_lb = result.problem.inner_lb + inner_ub = result.problem.inner_ub if inner_xs_names is None: inner_xs = None diff --git a/test/base/test_history.py b/test/base/test_history.py index 8338206aa..070231526 100644 --- a/test/base/test_history.py +++ b/test/base/test_history.py @@ -14,7 +14,9 @@ import pypesto import pypesto.optimize as optimize from pypesto import ( + CsvAmiciHistory, CsvHistory, + Hdf5AmiciHistory, Hdf5History, HistoryOptions, MemoryHistory, @@ -689,6 +691,70 @@ def test_hdf5_history_mp(): ) +def test_hdf5_amici_history(): + objective1 = pypesto.Objective( + fun=so.rosen, grad=so.rosen_der, hess=so.rosen_hess + ) + objective2 = load_amici_objective('conversion_reaction')[0] + lb = -2 * np.ones((1, 2)) + ub = 2 * np.ones((1, 2)) + problem1 = pypesto.Problem(objective=objective1, lb=lb, ub=ub) + problem2 = pypesto.Problem(objective=objective2, lb=lb, ub=ub) + + optimizer = pypesto.optimize.ScipyOptimizer(options={'maxiter': 10}) + + with tempfile.TemporaryDirectory(dir=".") as tmpdirname: + for f_ext, amici_history_class in zip( + [".csv", ".hdf5"], [CsvAmiciHistory, Hdf5AmiciHistory] + ): + _, fn = tempfile.mkstemp(f_ext, '{id}', dir=f"{tmpdirname}") + + history_options = pypesto.HistoryOptions( + trace_record=True, storage_file=fn + ) + + result1 = pypesto.optimize.minimize( + problem=problem1, + optimizer=optimizer, + n_starts=1, + history_options=history_options, + progress_bar=False, + ) + assert not isinstance( + result1.optimize_result.list[0].history, amici_history_class + ) + os.remove(fn) + os.remove(fn.replace('{id}', '0')) + + # optimizing with amici history saved in hdf5 + result2 = pypesto.optimize.minimize( + problem=problem2, + optimizer=optimizer, + n_starts=1, + history_options=history_options, + progress_bar=False, + ) + history = result2.optimize_result.list[0].history + assert isinstance(history, amici_history_class) + + assert np.all( + history.get_cpu_time_total_trace() + >= history.get_preeq_time_trace() + ) + assert np.all( + history.get_cpu_time_total_trace() + >= history.get_preeq_timeB_trace() + ) + assert np.all( + history.get_cpu_time_total_trace() + >= history.get_posteq_time_trace() + ) + assert np.all( + history.get_cpu_time_total_trace() + >= history.get_posteq_timeB_trace() + ) + + def test_trim_history(): """ Test whether the history gets correctly trimmed to be monotonically