Merge pull request #383 from pybop-team/minkowski-cost

Adds Minkowski Cost

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#393](https://github.com/pybop-team/PyBOP/pull/383) - Adds Minkowski and SumofPower cost classes, with an example and corresponding tests.
 - [#403](https://github.com/pybop-team/PyBOP/pull/403/) - Adds lychee link checking action.
 
 ## Bug Fixes

examples/scripts/spm_AdamW.py

@@ -53,7 +53,7 @@ def noise(sigma):
 problem = pybop.FittingProblem(
     model, parameters, dataset, signal=signal, init_soc=init_soc
 )
-cost = pybop.RootMeanSquaredError(problem)
+cost = pybop.Minkowski(problem, p=2)
 optim = pybop.AdamW(
     cost,
     verbose=True,

examples/scripts/spm_MAP.py

@@ -2,12 +2,16 @@
 
 import pybop
 
+# Set variables
+sigma = 0.002
+init_soc = 0.7
+
 # Construct and update initial parameter values
 parameter_set = pybop.ParameterSet.pybamm("Chen2020")
 parameter_set.update(
     {
-        "Negative electrode active material volume fraction": 0.63,
-        "Positive electrode active material volume fraction": 0.51,
+        "Negative electrode active material volume fraction": 0.43,
+        "Positive electrode active material volume fraction": 0.56,
     }
 )
 
@@ -18,38 +22,51 @@
 parameters = pybop.Parameters(
     pybop.Parameter(
         "Negative electrode active material volume fraction",
-        prior=pybop.Gaussian(0.6, 0.05),
-        bounds=[0.5, 0.8],
+        prior=pybop.Uniform(0.3, 0.8),
+        bounds=[0.3, 0.8],
+        initial_value=0.653,
         true_value=parameter_set["Negative electrode active material volume fraction"],
     ),
     pybop.Parameter(
         "Positive electrode active material volume fraction",
-        prior=pybop.Gaussian(0.48, 0.05),
+        prior=pybop.Uniform(0.3, 0.8),
         bounds=[0.4, 0.7],
+        initial_value=0.657,
         true_value=parameter_set["Positive electrode active material volume fraction"],
     ),
 )
 
-# Generate data
-sigma = 0.005
-t_eval = np.arange(0, 900, 3)
-values = model.predict(t_eval=t_eval)
-corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval))
+# Generate data and corrupt it with noise
+experiment = pybop.Experiment(
+    [
+        (
+            "Discharge at 0.5C for 3 minutes (4 second period)",
+            "Charge at 0.5C for 3 minutes (4 second period)",
+        ),
+    ]
+)
+values = model.predict(
+    init_soc=init_soc, experiment=experiment, inputs=parameters.true_value()
+)
+corrupt_values = values["Voltage [V]"].data + np.random.normal(
+    0, sigma, len(values["Voltage [V]"].data)
+)
 
 # Form dataset
 dataset = pybop.Dataset(
     {
-        "Time [s]": t_eval,
+        "Time [s]": values["Time [s]"].data,
         "Current function [A]": values["Current [A]"].data,
         "Voltage [V]": corrupt_values,
     }
 )
 
 # Generate problem, cost function, and optimisation class
-problem = pybop.FittingProblem(model, parameters, dataset)
+problem = pybop.FittingProblem(model, parameters, dataset, init_soc=init_soc)
 cost = pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=sigma)
 optim = pybop.AdamW(
     cost,
+    sigma0=0.05,
     max_unchanged_iterations=20,
     min_iterations=20,
     max_iterations=100,
@@ -61,7 +78,7 @@
 print("Estimated parameters:", x)
 
 # Plot the timeseries output
-pybop.quick_plot(problem, problem_inputs=x[0:2], title="Optimised Comparison")
+pybop.quick_plot(problem, problem_inputs=x, title="Optimised Comparison")
 
 # Plot convergence
 pybop.plot_convergence(optim)
@@ -73,5 +90,5 @@
 pybop.plot2d(cost, steps=15)
 
 # Plot the cost landscape with optimisation path
-bounds = np.asarray([[0.55, 0.77], [0.48, 0.68]])
+bounds = np.asarray([[0.35, 0.7], [0.45, 0.625]])
 pybop.plot2d(optim, bounds=bounds, steps=15)

examples/scripts/spm_SNES.py

@@ -35,7 +35,7 @@
 
 # Generate problem, cost function, and optimisation class
 problem = pybop.FittingProblem(model, parameters, dataset)
-cost = pybop.SumSquaredError(problem)
+cost = pybop.SumofPower(problem, p=2)
 optim = pybop.SNES(cost, max_iterations=100)
 
 x, final_cost = optim.run()

pybop/__init__.py

@@ -85,6 +85,8 @@
 from .costs.fitting_costs import (
     RootMeanSquaredError,
     SumSquaredError,
+    Minkowski,
+    SumofPower,
     ObserverCost,
 )
 from .costs.design_costs import (

pybop/costs/_likelihoods.py

@@ -38,19 +38,17 @@ def __init__(self, problem: BaseProblem, sigma0: Union[list[float], float]):
         self.sigma2 = sigma0**2.0
         self._offset = -0.5 * self.n_time_data * np.log(2 * np.pi * self.sigma2)
         self._multip = -1 / (2.0 * self.sigma2)
-        self._dl = np.ones(self.n_parameters)
 
     def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> float:
         """
         Evaluates the Gaussian log-likelihood for the given parameters with known sigma.
         """
         y = self.problem.evaluate(inputs)
-        if any(
-            len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
-        ):
-            return -np.inf  # prediction length doesn't match target
 
-        e = np.sum(
+        if not self.verify_prediction(y):
+            return -np.inf
+
+        e = np.asarray(
             [
                 np.sum(
                     self._offset
@@ -60,18 +58,16 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
             ]
         )
 
-        return e if self.n_outputs != 1 else e.item()
+        return e.item() if self.n_outputs == 1 else np.sum(e)
 
     def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
         """
         Calls the problem.evaluateS1 method and calculates the log-likelihood and gradient.
         """
         y, dy = self.problem.evaluateS1(inputs)
 
-        if any(
-            len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
-        ):
-            return -np.inf, -self._dl
+        if not self.verify_prediction(y):
+            return -np.inf, -self._de * np.ones(self.n_parameters)
 
         likelihood = self._evaluate(inputs)
 
@@ -125,7 +121,6 @@ def __init__(
         self.sigma = Parameters()
         self._add_sigma_parameters(sigma0)
         self.parameters.join(self.sigma)
-        self._dl = np.ones(self.n_parameters)
 
     def _add_sigma_parameters(self, sigma0):
         sigma0 = [sigma0] if not isinstance(sigma0, list) else sigma0
@@ -195,12 +190,10 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
             return -np.inf
 
         y = self.problem.evaluate(self.problem.parameters.as_dict())
-        if any(
-            len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
-        ):
-            return -np.inf  # prediction length doesn't match target
+        if not self.verify_prediction(y):
+            return -np.inf
 
-        e = np.sum(
+        e = np.asarray(
             [
                 np.sum(
                     self._logpi
@@ -212,7 +205,7 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
             ]
         )
 
-        return e if self.n_outputs != 1 else e.item()
+        return e.item() if self.n_outputs == 1 else np.sum(e)
 
     def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
         """
@@ -232,13 +225,11 @@ def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
 
         sigma = self.sigma.current_value()
         if np.any(sigma <= 0):
-            return -np.inf, -self._dl
+            return -np.inf, -self._de * np.ones(self.n_parameters)
 
         y, dy = self.problem.evaluateS1(self.problem.parameters.as_dict())
-        if any(
-            len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
-        ):
-            return -np.inf, -self._dl
+        if not self.verify_prediction(y):
+            return -np.inf, -self._de * np.ones(self.n_parameters)
 
         likelihood = self._evaluate(inputs)
 
@@ -302,11 +293,14 @@ def _evaluate(self, inputs: Inputs, grad=None) -> float:
         float
             The maximum a posteriori cost.
         """
-        log_likelihood = self.likelihood._evaluate(inputs)
         log_prior = sum(
             self.parameters[key].prior.logpdf(value) for key, value in inputs.items()
         )
 
+        if not np.isfinite(log_prior).any():
+            return -np.inf
+
+        log_likelihood = self.likelihood._evaluate(inputs)
         posterior = log_likelihood + log_prior
         return posterior
 
@@ -331,10 +325,13 @@ def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
         ValueError
             If an error occurs during the calculation of the cost or gradient.
         """
-        log_likelihood, dl = self.likelihood._evaluateS1(inputs)
         log_prior = sum(
             self.parameters[key].prior.logpdf(value) for key, value in inputs.items()
         )
+        if not np.isfinite(log_prior).any():
+            return -np.inf, -self._de * np.ones(self.n_parameters)
+
+        log_likelihood, dl = self.likelihood._evaluateS1(inputs)
 
         # Compute a finite difference approximation of the gradient of the log prior
         delta = self.parameters.initial_value() * self.gradient_step

pybop/costs/base_cost.py

@@ -1,3 +1,5 @@
+from typing import Union
+
 from pybop import BaseProblem
 from pybop.parameters.parameter import Inputs, Parameters
 
@@ -25,6 +27,7 @@ class BaseCost:
     def __init__(self, problem=None):
         self.parameters = Parameters()
         self.problem = problem
+        self.set_fail_gradient()
         if isinstance(self.problem, BaseProblem):
             self._target = self.problem._target
             self.parameters.join(self.problem.parameters)
@@ -35,20 +38,20 @@ def __init__(self, problem=None):
     def n_parameters(self):
         return len(self.parameters)
 
-    def __call__(self, x, grad=None):
+    def __call__(self, inputs: Union[Inputs, list], grad=None):
         """
         Call the evaluate function for a given set of parameters.
         """
-        return self.evaluate(x, grad)
+        return self.evaluate(inputs, grad)
 
-    def evaluate(self, x, grad=None):
+    def evaluate(self, inputs: Union[Inputs, list], grad=None):
         """
         Call the evaluate function for a given set of parameters.
 
         Parameters
         ----------
-        x : array-like
-            The parameters for which to evaluate the cost.
+        inputs : Inputs or array-like
+            The parameters for which to compute the cost and gradient.
         grad : array-like, optional
             An array to store the gradient of the cost function with respect
             to the parameters.
@@ -63,7 +66,7 @@ def evaluate(self, x, grad=None):
         ValueError
             If an error occurs during the calculation of the cost.
         """
-        inputs = self.parameters.verify(x)
+        inputs = self.parameters.verify(inputs)
 
         try:
             return self._evaluate(inputs, grad)
@@ -100,27 +103,27 @@ def _evaluate(self, inputs: Inputs, grad=None):
         """
         raise NotImplementedError
 
-    def evaluateS1(self, x):
+    def evaluateS1(self, inputs: Union[Inputs, list]):
         """
         Call _evaluateS1 for a given set of parameters.
 
         Parameters
         ----------
-        x : array-like
+        inputs : Inputs or array-like
             The parameters for which to compute the cost and gradient.
 
         Returns
         -------
         tuple
             A tuple containing the cost and the gradient. The cost is a float,
-            and the gradient is an array-like of the same length as `x`.
+            and the gradient is an array-like of the same length as `inputs`.
 
         Raises
         ------
         ValueError
             If an error occurs during the calculation of the cost or gradient.
         """
-        inputs = self.parameters.verify(x)
+        inputs = self.parameters.verify(inputs)
 
         try:
             return self._evaluateS1(inputs)
@@ -144,11 +147,48 @@ def _evaluateS1(self, inputs: Inputs):
         -------
         tuple
             A tuple containing the cost and the gradient. The cost is a float,
-            and the gradient is an array-like of the same length as `x`.
+            and the gradient is an array-like of the same length as `inputs`.
 
         Raises
         ------
         NotImplementedError
             If the method has not been implemented by the subclass.
         """
         raise NotImplementedError
+
+    def set_fail_gradient(self, de: float = 1.0):
+        """
+        Set the fail gradient to a specified value.
+
+        The fail gradient is used if an error occurs during the calculation
+        of the gradient. This method allows updating the default gradient value.
+
+        Parameters
+        ----------
+        de : float
+            The new fail gradient value to be used.
+        """
+        if not isinstance(de, float):
+            de = float(de)
+        self._de = de
+
+    def verify_prediction(self, y):
+        """
+        Verify that the prediction matches the target data.
+
+        Parameters
+        ----------
+        y : dict
+            The model predictions.
+
+        Returns
+        -------
+        bool
+            True if the prediction matches the target data, otherwise False.
+        """
+        if any(
+            len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
+        ):
+            return False
+
+        return True