diff --git a/CHANGELOG.md b/CHANGELOG.md index 657e7ec4..a305f917 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ ## Optimisations +- [#588](https://github.com/pybop-team/PyBOP/pull/588) - Makes `minimising` a property of `BaseOptimiser` set by the cost class. - [#512](https://github.com/pybop-team/PyBOP/pull/513) - Refactors `LogPosterior` with attributes pointing to composed likelihood object. - [#551](https://github.com/pybop-team/PyBOP/pull/551) - Refactors Optimiser arguments, `population_size` and `max_iterations` as default args, improves optimiser docstrings diff --git a/examples/standalone/optimiser.py b/examples/standalone/optimiser.py index 2e577958..93802ebe 100644 --- a/examples/standalone/optimiser.py +++ b/examples/standalone/optimiser.py @@ -67,12 +67,10 @@ def callback(x): ) return OptimisationResult( + optim=self, x=result.x, - cost=self.cost, - final_cost=self.cost(result.x), n_iterations=result.nit, scipy_result=result, - optim=self, ) def name(self): diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 2472f351..8fc23809 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -18,6 +18,7 @@ class BaseLikelihood(BaseCost): def __init__(self, problem: BaseProblem): super().__init__(problem) self.n_data = problem.n_data + self.minimising = False class BaseMetaLikelihood(BaseLikelihood): diff --git a/pybop/costs/_weighted_cost.py b/pybop/costs/_weighted_cost.py index 39c42ef6..9ead437f 100644 --- a/pybop/costs/_weighted_cost.py +++ b/pybop/costs/_weighted_cost.py @@ -2,7 +2,7 @@ import numpy as np -from pybop import BaseCost, BaseLikelihood, DesignCost +from pybop import BaseCost class WeightedCost(BaseCost): @@ -32,9 +32,6 @@ def __init__(self, *costs, weights: Optional[list[float]] = None): self.costs = [cost for cost in costs] if len(set(type(cost.problem) for cost in self.costs)) > 1: raise TypeError("All problems must be of the same class type.") - self.minimising = not any( - isinstance(cost, (BaseLikelihood, DesignCost)) for cost in self.costs - ) # Check if weights are provided if weights is not None: @@ -62,6 +59,14 @@ def __init__(self, *costs, weights: Optional[list[float]] = None): for cost in self.costs: self.join_parameters(cost.parameters) + # Apply the minimising property from each cost + for i, cost in enumerate(self.costs): + self.weights[i] = self.weights[i] * (1 if cost.minimising else -1) + if all(not cost.minimising for cost in self.costs): + # If all costs are maximising, convert the weighted cost to maximising + self.weights = -self.weights + self.minimising = False + # Weighted costs do not use this functionality self._has_separable_problem = False diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 7ec08f0e..2edf1faa 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -32,6 +32,9 @@ class BaseCost: _de : float The gradient of the cost function to use if an error occurs during evaluation. Defaults to 1.0. + minimising : bool, optional, default=True + If False, switches the sign of the cost and gradient to perform maximisation + instead of minimisation. """ def __init__(self, problem: Optional[BaseProblem] = None): @@ -43,6 +46,7 @@ def __init__(self, problem: Optional[BaseProblem] = None): self.y = None self.dy = None self._de = 1.0 + self.minimising = True if isinstance(self.problem, BaseProblem): self._target = self.problem.target self._parameters.join(self.problem.parameters) @@ -58,6 +62,7 @@ def __call__( inputs: Union[Inputs, list], calculate_grad: bool = False, apply_transform: bool = False, + for_optimiser: bool = False, ) -> Union[float, tuple[float, np.ndarray]]: """ This method calls the forward model via problem.evaluate(inputs), @@ -73,6 +78,9 @@ def __call__( cost is computed. apply_transform : bool, optional, default=False If True, applies a transformation to the inputs before evaluating the model. + for_optimiser : bool, optional, default=False + If True, returns the cost value if self.minimising=True and the negative of + the cost value if self.minimising=False (i.e. the cost is being maximised). Returns ------- @@ -96,6 +104,15 @@ def __call__( else: model_inputs = inputs + # Check whether we are maximising or minimising via: + # | `minimising` | `self.minimising` | `for_optimiser` | + # |--------------|-------------------|-----------------| + # | `True` | `True` | `True` | + # | `True` | `True` | `False` | + # | `False` | `False` | `True` | + # | `True` | `False` | `False` | + minimising = self.minimising or not for_optimiser + # Validate inputs, update parameters model_inputs = self.parameters.verify(model_inputs) self.parameters.update(values=list(model_inputs.values())) @@ -110,10 +127,15 @@ def __call__( jac = self.transformation.jacobian(inputs) grad = np.matmul(grad, jac) - return cost, grad + return cost * (1 if minimising else -1), grad * ( + 1 if minimising else -1 + ) y = self.problem.evaluate(self.problem.parameters.as_dict()) - return self.compute(y, dy=dy, calculate_grad=calculate_grad) + + return self.compute(y, dy=dy, calculate_grad=calculate_grad) * ( + 1 if minimising else -1 + ) def compute(self, y: dict, dy: ndarray, calculate_grad: bool = False): """ diff --git a/pybop/costs/design_costs.py b/pybop/costs/design_costs.py index 49714cff..c9651974 100644 --- a/pybop/costs/design_costs.py +++ b/pybop/costs/design_costs.py @@ -19,7 +19,7 @@ class DesignCost(BaseCost): def __init__(self, problem): """ - Initialises the gravimetric energy density calculator with a problem. + Initialises the design cost calculator with a problem. Parameters ---------- @@ -27,15 +27,14 @@ def __init__(self, problem): The problem instance containing the model and data. """ super().__init__(problem) - self.problem = problem + self.minimising = False class GravimetricEnergyDensity(DesignCost): """ Calculates the gravimetric energy density (specific energy) of a battery cell, when applied to a normalised discharge from upper to lower voltage limits. The - goal of maximising the energy density is achieved by setting minimising = False - in the optimiser settings. + goal of maximising the energy density is achieved with self.minimising=False. The gravimetric energy density [Wh.kg-1] is calculated as @@ -92,8 +91,7 @@ class VolumetricEnergyDensity(DesignCost): """ Calculates the (volumetric) energy density of a battery cell, when applied to a normalised discharge from upper to lower voltage limits. The goal of maximising - the energy density is achieved by setting minimising = False in the optimiser - settings. + the energy density is achieved with self.minimising = False. The volumetric energy density [Wh.m-3] is calculated as @@ -150,8 +148,7 @@ class GravimetricPowerDensity(DesignCost): """ Calculates the gravimetric power density (specific power) of a battery cell, when applied to a discharge from upper to lower voltage limits. The goal of - maximising the power density is achieved by setting minimising = False in the - optimiser settings. + maximising the power density is achieved with self.minimising=False. The time-averaged gravimetric power density [W.kg-1] is calculated as @@ -214,7 +211,7 @@ class VolumetricPowerDensity(DesignCost): """ Calculates the (volumetric) power density of a battery cell, when applied to a discharge from upper to lower voltage limits. The goal of maximising the power - density is achieved by setting minimising = False in the optimiser settings. + density is achieved with self.minimising=False. The time-averaged volumetric power density [W.m-3] is calculated as diff --git a/pybop/optimisers/base_optimiser.py b/pybop/optimisers/base_optimiser.py index fd8375db..88f2e3e6 100644 --- a/pybop/optimisers/base_optimiser.py +++ b/pybop/optimisers/base_optimiser.py @@ -4,15 +4,7 @@ import numpy as np from scipy.optimize import OptimizeResult -from pybop import ( - BaseCost, - BaseLikelihood, - DesignCost, - Inputs, - Parameter, - Parameters, - WeightedCost, -) +from pybop import BaseCost, Inputs, Parameter, Parameters class BaseOptimiser: @@ -43,9 +35,6 @@ class BaseOptimiser: Not all methods will use this information. verbose : bool, optional If True, the optimisation progress is printed (default: False). - minimising : bool, optional - If True, the target is to minimise the cost, else target is to maximise by minimising - the negative cost (default: True). physical_viability : bool, optional If True, the feasibility of the optimised parameters is checked (default: False). allow_infeasible_solutions : bool, optional @@ -62,13 +51,13 @@ def __init__( # First set attributes to default values self.parameters = Parameters() self.x0 = optimiser_kwargs.get("x0", []) - self.log = dict(x=[], x_best=[], cost=[], cost_best=[], x0=[]) + self.log = dict(x=[], x_best=[], x_search=[], x0=[], cost=[], cost_best=[]) self.bounds = None self.sigma0 = 0.02 self.verbose = True - self.minimising = True self._transformation = None self._needs_sensitivities = False + self._minimising = True self.physical_viability = False self.allow_infeasible_solutions = False self.default_max_iterations = 1000 @@ -79,10 +68,7 @@ def __init__( self.parameters = self.cost.parameters self._transformation = self.cost.transformation self.set_allow_infeasible_solutions() - if isinstance(cost, WeightedCost): - self.minimising = cost.minimising - if isinstance(cost, (BaseLikelihood, DesignCost)): - self.minimising = False + self._minimising = self.cost.minimising else: try: @@ -98,7 +84,6 @@ def __init__( self.parameters.add( Parameter(name=f"Parameter {i}", initial_value=value) ) - self.minimising = True except Exception as e: raise Exception( @@ -147,7 +132,6 @@ def set_base_options(self): # Set other options self.verbose = self.unset_options.pop("verbose", self.verbose) - self.minimising = self.unset_options.pop("minimising", self.minimising) if "allow_infeasible_solutions" in self.unset_options.keys(): self.set_allow_infeasible_solutions( self.unset_options.pop("allow_infeasible_solutions") @@ -169,6 +153,38 @@ def _set_up_optimiser(self): """ raise NotImplementedError + def cost_call( + self, + x: Union[Inputs, list], + calculate_grad: bool = False, + ) -> Union[float, tuple[float, np.ndarray]]: + """ + Call the cost function to minimise, applying any given transformation to the + input parameters. + + Parameters + ---------- + x : Inputs or list-like + The input parameters for which the cost and optionally the gradient + will be computed. + calculate_grad : bool, optional, default=False + If True, both the cost and gradient will be computed. Otherwise, only the + cost is computed. + + Returns + ------- + float or tuple + - If `calculate_grad` is False, returns the computed cost (float). + - If `calculate_grad` is True, returns a tuple containing the cost (float) + and the gradient (np.ndarray). + """ + return self.cost( + x, + calculate_grad=calculate_grad, + apply_transform=True, + for_optimiser=True, + ) + def run(self): """ Run the optimisation and return the optimised parameters and final cost. @@ -243,6 +259,7 @@ def apply_transformation(values): if x is not None: x = convert_to_list(x) + self.log["x_search"].extend(x) x = apply_transformation(x) self.log["x"].extend(x) @@ -252,10 +269,17 @@ def apply_transformation(values): if cost is not None: cost = convert_to_list(cost) + cost = [ + internal_cost * (1 if self.minimising else -1) for internal_cost in cost + ] self.log["cost"].extend(cost) if cost_best is not None: cost_best = convert_to_list(cost_best) + cost_best = [ + internal_cost * (1 if self.minimising else -1) + for internal_cost in cost_best + ] self.log["cost_best"].extend(cost_best) if x0 is not None: @@ -302,6 +326,10 @@ def set_allow_infeasible_solutions(self, allow: bool = True): def needs_sensitivities(self): return self._needs_sensitivities + @property + def minimising(self): + return self._minimising + class OptimisationResult: """ @@ -321,22 +349,26 @@ class OptimisationResult: def __init__( self, + optim: BaseOptimiser, x: Union[Inputs, np.ndarray] = None, - cost: Union[BaseCost, None] = None, final_cost: Optional[float] = None, n_iterations: Optional[int] = None, - optim: Optional[BaseOptimiser] = None, time: Optional[float] = None, scipy_result=None, ): - self.x = x - self.cost = cost + self.optim = optim + self.cost = self.optim.cost + self.minimising = self.optim.minimising + self._transformation = self.optim._transformation # noqa: SLF001 + + self.x = self._transformation.to_model(x) if self._transformation else x self.final_cost = ( - final_cost if final_cost is not None else self._calculate_final_cost() + final_cost * (1 if self.minimising else -1) + if final_cost is not None + else self._calculate_final_cost() ) self.n_iterations = n_iterations self.scipy_result = scipy_result - self.optim = optim self.time = time if isinstance(self.optim, BaseOptimiser): self.x0 = self.optim.parameters.initial_value() @@ -458,7 +490,7 @@ def add_run(self, result: OptimisationResult): def best_run(self) -> Optional[OptimisationResult]: """Returns the result with the best final cost.""" valid_results = [res for res in self.results if res.final_cost is not None] - if self.results[0].optim.minimising is True: + if self.results[0].minimising is True: return min(valid_results, key=lambda res: res.final_cost) return max(valid_results, key=lambda res: res.final_cost) diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index bdf5b828..70de351d 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -198,13 +198,7 @@ def _run(self): # Choose method to evaluate def fun(x): - if self._needs_sensitivities: - L, dl = self.cost(x, calculate_grad=True, apply_transform=True) - else: - L = self.cost(x, apply_transform=True) - dl = None - sign = -1 if not self.minimising else 1 - return (sign * L, sign * dl) if dl is not None else sign * L + return self.cost_call(x, calculate_grad=self._needs_sensitivities) # Create evaluator object if self._parallel: @@ -222,9 +216,6 @@ def fun(x): # Keep track of current best and best-guess scores. fb = fg = np.inf - # Internally we always minimise! Keep a 2nd value to show the user. - fg_user = (fb, fg) if self.minimising else (-fb, -fg) - # Keep track of the last significant change f_sig = np.inf @@ -244,7 +235,6 @@ def fun(x): # Update the scores fb = self.optimiser.f_best() fg = self.optimiser.f_guessed() - fg_user = (fb, fg) if self.minimising else (-fb, -fg) # Check for significant changes against the absolute and relative tolerance f_new = fg if self._use_f_guessed else fb @@ -263,8 +253,8 @@ def fun(x): self.log_update( x=xs, x_best=self.optimiser.x_best(), - cost=_fs if self.minimising else [-x for x in _fs], - cost_best=fb if self.minimising else -fb, + cost=_fs, + cost_best=fb, ) # Check stopping criteria: @@ -327,14 +317,11 @@ def fun(x): # Show last result and exit print("\n" + "-" * 40) print("Unexpected termination.") - print("Current score: " + str(fg_user)) + print("Current score: " + str((fb, fg))) print("Current position:") - # Show current parameters - x_user = self.optimiser.x_guessed() - if self._transformation: - x_user = self._transformation.to_model(x_user) - for p in x_user: + # Show current parameters (with any transformation applied) + for p in self.optimiser.x_guessed(): print(PintsStrFloat(p)) print("-" * 40) raise @@ -355,16 +342,11 @@ def fun(x): x = self.optimiser.x_best() f = self.optimiser.f_best() - # Inverse transform search parameters - if self._transformation: - x = self._transformation.to_model(x) - return OptimisationResult( + optim=self, x=x, - cost=self.cost, - final_cost=f if self.minimising else -f, + final_cost=f, n_iterations=self._iterations, - optim=self, time=total_time, ) diff --git a/pybop/optimisers/scipy_optimisers.py b/pybop/optimisers/scipy_optimisers.py index 838b9eb8..2547c84b 100644 --- a/pybop/optimisers/scipy_optimisers.py +++ b/pybop/optimisers/scipy_optimisers.py @@ -88,14 +88,10 @@ def _run(self): nit = -1 return OptimisationResult( - x=self._transformation.to_model(result.x) - if self._transformation - else result.x, - cost=self.cost, - final_cost=self.cost(result.x, apply_transform=True), + optim=self, + x=result.x, n_iterations=nit, scipy_result=result, - optim=self, time=total_time, ) @@ -200,7 +196,7 @@ def cost_wrapper(self, x): Scale the cost function, preserving the sign convention, and eliminate nan values """ if not self._options["jac"]: - cost = self.cost(x, apply_transform=True) + cost = self.cost_call(x) self.log_update(x=[x], cost=cost) scaled_cost = cost / self._cost0 if np.isinf(scaled_cost): @@ -208,13 +204,11 @@ def cost_wrapper(self, x): scaled_cost = np.sign(cost) * ( 1 + 0.9**self.inf_count ) # for fake finite gradient - return scaled_cost * (1 if self.minimising else -1) + return scaled_cost - L, dl = self.cost(x, calculate_grad=True, apply_transform=True) + L, dl = self.cost_call(x, calculate_grad=True) self.log_update(x=[x], cost=L) - scaled_L = L / self._cost0 - scaled_dl = dl / self._cost0 - return (scaled_L, scaled_dl) if self.minimising else (-scaled_L, -scaled_dl) + return (L / self._cost0, dl / self._cost0) def _run_optimiser(self): """ @@ -237,14 +231,10 @@ def base_callback(intermediate_result: Union[OptimizeResult, np.ndarray]): """ if isinstance(intermediate_result, OptimizeResult): x_best = intermediate_result.x - cost_best = ( - intermediate_result.fun - * self._cost0 - * (1 if self.minimising else -1) - ) + cost_best = intermediate_result.fun * self._cost0 else: x_best = intermediate_result - cost_best = self.cost(x_best, apply_transform=True) + cost_best = self.cost_call(x_best) self.log_update(x_best=x_best, cost_best=cost_best) @@ -255,7 +245,7 @@ def base_callback(intermediate_result: Union[OptimizeResult, np.ndarray]): ) # Compute the absolute initial cost and resample if required - self._cost0 = np.abs(self.cost(self.x0, apply_transform=True)) + self._cost0 = np.abs(self.cost_call(self.x0)) if np.isinf(self._cost0): for _i in range(1, self.num_resamples): try: @@ -267,7 +257,7 @@ def base_callback(intermediate_result: Union[OptimizeResult, np.ndarray]): stacklevel=2, ) break - self._cost0 = np.abs(self.cost(self.x0, apply_transform=True)) + self._cost0 = np.abs(self.cost_call(self.x0)) if not np.isinf(self._cost0): break if np.isinf(self._cost0): @@ -433,13 +423,13 @@ def _run_optimiser(self): def callback(intermediate_result: OptimizeResult): self.log_update( x_best=intermediate_result.x, - cost_best=intermediate_result.fun * (1 if self.minimising else -1), + cost_best=intermediate_result.fun, ) def cost_wrapper(x): - cost = self.cost(x, apply_transform=True) + cost = self.cost_call(x) self.log_update(x=[x], cost=cost) - return cost * (1 if self.minimising else -1) + return cost return differential_evolution( cost_wrapper, diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index 134a7d36..ea5a5792 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -357,7 +357,7 @@ def test_weighted_fitting_cost(self, noisy_problem): # Test LogPosterior explicitly cost4 = pybop.LogPosterior(pybop.GaussianLogLikelihood(problem)) - weighted_cost_4 = pybop.WeightedCost(cost1, cost4, weights=[1, -1 / weight]) + weighted_cost_4 = pybop.WeightedCost(cost1, cost4, weights=[1, 1 / weight]) assert weighted_cost_4.has_identical_problems is True assert weighted_cost_4.has_separable_problem is False sigma = 0.01 diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index b9b5c71f..2f865f46 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -142,9 +142,15 @@ def test_no_optimisation_parameters(self, model, dataset): @pytest.mark.unit def test_optimiser_kwargs(self, cost, optimiser): def assert_log_update(optim): - optim.log_update(x=[0.7], cost=[0.01]) - assert optim.log["x"][-1] == 0.7 - assert optim.log["cost"][-1] == 0.01 + x_search = 0.7 + optim.log_update(x=[x_search], cost=[0.01]) + assert optim.log["x_search"][-1] == x_search + assert optim.log["cost"][-1] == 0.01 * (1 if optim.minimising else -1) + assert ( + optim.log["x"][-1] == optim._transformation.to_model(x_search) + if optim._transformation + else x_search + ) def check_max_iterations(optim): results = optim.run() @@ -621,7 +627,8 @@ def test_unphysical_result(self, cost): @pytest.mark.unit def test_optimisation_results(self, cost): # Construct OptimisationResult - results = pybop.OptimisationResult(x=[1e-3], cost=cost, n_iterations=1) + optim = pybop.Optimisation(cost=cost) + results = pybop.OptimisationResult(optim=optim, x=[1e-3], n_iterations=1) # Asserts assert results.x[0] == 1e-3 @@ -632,4 +639,4 @@ def test_optimisation_results(self, cost): with pytest.raises( ValueError, match="Optimised parameters do not produce a finite cost value" ): - pybop.OptimisationResult(x=[1e-5], cost=cost, n_iterations=1) + pybop.OptimisationResult(optim=optim, x=[1e-5], n_iterations=1)