Make likelihoods a type of cost

examples/standalone/cost.py

@@ -18,7 +18,7 @@ class StandaloneCost(pybop.BaseCost):
         BaseCost interface.
     x0 : array-like
         The initial guess for the optimization problem, set to [4.2].
-    n_parameters : int
+    _n_parameters : int
         The number of parameters in the model, which is 1 in this case.
     bounds : dict
         A dictionary containing the lower and upper bounds for the parameter,
@@ -40,7 +40,7 @@ def __init__(self, problem=None):
         super().__init__(problem)
 
         self.x0 = np.array([4.2])
-        self.n_parameters = len(self.x0)
+        self._n_parameters = len(self.x0)
 
         self.bounds = dict(
             lower=[-1],

pybop/__init__.py

@@ -28,11 +28,6 @@
 #
 from ._problem import BaseProblem, FittingProblem, DesignProblem
 
-#
-# Likelihood classes
-#
-from ._likelihoods import BaseLikelihood, GaussianLogLikelihood, GaussianLogLikelihoodKnownSigma
-
 #
 # Cost function class
 #
@@ -41,13 +36,17 @@
     RootMeanSquaredError,
     SumSquaredError,
     ObserverCost,
-    ProbabilityCost,
 )
 from .costs.design_costs import (
     DesignCost,
     GravimetricEnergyDensity,
     VolumetricEnergyDensity,
 )
+from .costs._likelihoods import (
+    BaseLikelihood,
+    GaussianLogLikelihood,
+    GaussianLogLikelihoodKnownSigma,
+)
 
 #
 # Dataset class

pybop/_optimisation.py

@@ -29,7 +29,7 @@ class Optimisation:
         Initial parameter values for the optimization.
     bounds : dict
         Dictionary containing the parameter bounds with keys 'lower' and 'upper'.
-    n_parameters : int
+    _n_parameters : int
         Number of parameters in the optimization problem.
     sigma0 : float or sequence
         Initial step size or standard deviation for the optimiser.
@@ -40,27 +40,24 @@ class Optimisation:
     def __init__(
         self,
         cost,
-        x0=None,
         optimiser=None,
         sigma0=None,
         verbose=False,
         physical_viability=True,
         allow_infeasible_solutions=True,
     ):
         self.cost = cost
-        self.x0 = x0
+        self.x0 = cost.x0
         self.optimiser = optimiser
         self.verbose = verbose
         self.bounds = cost.bounds
         self.sigma0 = sigma0 or cost.sigma0
-        self.n_parameters = cost.n_parameters
+        self._n_parameters = cost._n_parameters
         self.physical_viability = physical_viability
         self.allow_infeasible_solutions = allow_infeasible_solutions
         self.log = []
 
-        # Catch x0, and convert to pints vector
-        if x0 is None:
-            self.x0 = cost.x0
+        # Convert x0 to pints vector
         self._x0 = pints.vector(self.x0)
 
         # Set whether to allow infeasible locations
@@ -77,12 +74,10 @@ def __init__(
         self._transformation = None
 
         # Check if minimising or maximising
-        self._minimising = not isinstance(cost, pybop.BaseLikelihood)
-        if self._minimising:
-            self._function = self.cost
-        else:
-            self._function = pybop.ProbabilityCost(cost)
-        del cost
+        if isinstance(cost, pybop.BaseLikelihood):
+            self.cost._minimising = False
+        self._minimising = self.cost._minimising
+        self._function = self.cost
 
         # Construct Optimiser
         self.pints = True

pybop/_likelihoods.py → pybop/costs/_likelihoods.py

@@ -1,27 +1,19 @@
 import numpy as np
+from pybop.costs.base_cost import BaseCost
 
 
-class BaseLikelihood:
+class BaseLikelihood(BaseCost):
     """
     Base class for likelihoods
     """
 
     def __init__(self, problem, sigma=None):
-        self.problem = problem
+        super(BaseLikelihood, self).__init__(problem)
         self._n_output = problem.n_outputs
         self._n_times = problem.n_time_data
         self.sigma0 = sigma or np.zeros(self._n_output)
-        self.x0 = problem.x0
-        self.bounds = problem.bounds
         self._n_parameters = problem.n_parameters
-        self._target = problem._target
-
-    def __call__(self, x):
-        """
-        Calls the problem.evaluate method and calculates
-        the log-likelihood
-        """
-        raise NotImplementedError
+        self.log_likelihood = problem
 
     def set_sigma(self, sigma):
         """
@@ -44,10 +36,6 @@ def get_n_parameters(self):
         """
         return self._n_parameters
 
-    @property
-    def n_parameters(self):
-        return self._n_parameters
-
 
 class GaussianLogLikelihoodKnownSigma(BaseLikelihood):
     """
@@ -68,15 +56,15 @@ def __init__(self, problem, sigma=None):
         self.sigma2 = self.sigma0**-2
         self._dl = np.ones(self._n_parameters)
 
-    def __call__(self, x):
+    def _evaluate(self, x, grad=None):
         """
         Calls the problem.evaluate method and calculates
         the log-likelihood
         """
         e = self._target - self.problem.evaluate(x)
         return np.sum(self._offset + self._multip * np.sum(e**2, axis=0))
 
-    def _evaluateS1(self, x):
+    def _evaluateS1(self, x, grad=None):
         """
         Calls the problem.evaluateS1 method and calculates
         the log-likelihood
@@ -116,7 +104,7 @@ def __init__(self, problem):
         self._logpi = -0.5 * self._n_times * np.log(2 * np.pi)
         self._dl = np.ones(self._n_parameters + self._n_output)
 
-    def __call__(self, x):
+    def _evaluate(self, x, grad=None):
         """
         Evaluates the Gaussian log-likelihood for the given parameters.
 
@@ -140,7 +128,7 @@ def __call__(self, x):
             - np.sum(e**2, axis=0) / (2.0 * sigma**2)
         )
 
-    def _evaluateS1(self, x):
+    def _evaluateS1(self, x, grad=None):
         """
         Calls the problem.evaluateS1 method and calculates
         the log-likelihood
@@ -163,7 +151,7 @@ def _evaluateS1(self, x):
                 )
             )
             e = self._target - y
-            likelihood = self.__call__(x)
+            likelihood = self._evaluate(x)
             dl = np.sum((sigma**-(2.0) * np.sum((e.T * dy.T), axis=2)), axis=1)
 
             # Add sigma gradient to dl

pybop/costs/base_cost.py

@@ -1,5 +1,4 @@
 from pybop import BaseProblem
-from pybop import BaseLikelihood
 
 
 class BaseCost:
@@ -26,7 +25,7 @@ class BaseCost:
         Initial standard deviation around ``x0``. Either a scalar value (one
         standard deviation for all coordinates) or an array with one entry
         per dimension. Not all methods will use this information.
-    n_parameters : int
+    _n_parameters : int
         The number of parameters in the model.
     n_outputs : int
         The number of outputs in the model.
@@ -37,19 +36,28 @@ def __init__(self, problem=None):
         self.x0 = None
         self.bounds = None
         self.sigma0 = None
+        self._minimising = True
         if isinstance(self.problem, BaseProblem):
             self._target = problem._target
             self.x0 = problem.x0
             self.bounds = problem.bounds
             self.sigma0 = problem.sigma0
-            self.n_parameters = problem.n_parameters
+            self._n_parameters = problem.n_parameters
             self.n_outputs = problem.n_outputs
-        elif isinstance(self.problem, BaseLikelihood):
-            self.log_likelihood = problem
+
+    @property
+    def n_parameters(self):
+        return self._n_parameters
 
     def __call__(self, x, grad=None):
         """
         Call the evaluate function for a given set of parameters.
+        """
+        return self.evaluate(x, grad)
+
+    def evaluate(self, x, grad=None):
+        """
+        Call the evaluate function for a given set of parameters.
 
         Parameters
         ----------
@@ -70,7 +78,10 @@ def __call__(self, x, grad=None):
             If an error occurs during the calculation of the cost.
         """
         try:
-            return self._evaluate(x, grad)
+            if self._minimising:
+                return self._evaluate(x, grad)
+            else:  # minimise the negative cost
+                return -self._evaluate(x, grad)
 
         except NotImplementedError as e:
             raise e
@@ -125,7 +136,10 @@ def evaluateS1(self, x):
             If an error occurs during the calculation of the cost or gradient.
         """
         try:
-            return self._evaluateS1(x)
+            if self._minimising:
+                return self._evaluateS1(x)
+            else:  # minimise the negative cost
+                return -self._evaluateS1(x)
 
         except NotImplementedError as e:
             raise e

pybop/costs/fitting_costs.py

@@ -70,13 +70,13 @@ def _evaluateS1(self, x):
         y, dy = self.problem.evaluateS1(x)
         if len(y) < len(self._target):
             e = np.float64(np.inf)
-            de = self._de * np.ones(self.n_parameters)
+            de = self._de * np.ones(self._n_parameters)
         else:
             dy = dy.reshape(
                 (
                     self.problem.n_time_data,
                     self.n_outputs,
-                    self.n_parameters,
+                    self._n_parameters,
                 )
             )
             r = y - self._target
@@ -177,13 +177,13 @@ def _evaluateS1(self, x):
         y, dy = self.problem.evaluateS1(x)
         if len(y) < len(self._target):
             e = np.float64(np.inf)
-            de = self._de * np.ones(self.n_parameters)
+            de = self._de * np.ones(self._n_parameters)
         else:
             dy = dy.reshape(
                 (
                     self.problem.n_time_data,
                     self.n_outputs,
-                    self.n_parameters,
+                    self._n_parameters,
                 )
             )
             r = y - self._target
@@ -208,61 +208,6 @@ def set_fail_gradient(self, de):
         self._de = de
 
 
-class ProbabilityCost(BaseCost):
-    """
-    Probability based cost function.
-
-    Changes the sign of the log likelihood to make it a cost function.
-
-    Inherits all parameters and attributes from ``BaseCost``.
-    """
-
-    def __init__(self, log_likelihood):
-        super(ProbabilityCost, self).__init__(log_likelihood)
-
-    def _evaluate(self, x, grad=None):
-        """
-        Calculate the probability based cost for a given set of parameters.
-
-        Parameters
-        ----------
-        x : array-like
-            The parameters for which to evaluate the cost.
-        grad : array-like, optional
-            An array to store the gradient of the cost function with respect
-            to the parameters.
-
-        Returns
-        -------
-        float
-            The probability based cost.
-        """
-        return -self.log_likelihood(x)
-
-    def _evaluateS1(self, x):
-        """
-        Compute the cost and its gradient with respect to the parameters.
-
-        Parameters
-        ----------
-        x : 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`.
-
-        Raises
-        ------
-        ValueError
-            If an error occurs during the calculation of the cost or gradient.
-        """
-        likelihood, dl = self.log_likelihood._evaluateS1(x)
-        return -likelihood, -dl
-
-
 class ObserverCost(BaseCost):
     """
     Observer cost function.

tests/unit/test_cost.py

@@ -70,16 +70,12 @@ def problem(self, model, parameters, dataset, signal, x0, request):
             pybop.RootMeanSquaredError,
             pybop.SumSquaredError,
             pybop.ObserverCost,
-            pybop.ProbabilityCost,
         ]
     )
     def cost(self, problem, request):
         cls = request.param
         if cls in [pybop.SumSquaredError, pybop.RootMeanSquaredError]:
             return cls(problem)
-        elif cls in [pybop.ProbabilityCost]:
-            likelihood = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma=0.01)
-            return cls(likelihood)
         elif cls in [pybop.ObserverCost]:
             inputs = {p.name: problem.x0[i] for i, p in enumerate(problem.parameters)}
             state = problem._model.reinit(inputs)
@@ -144,7 +140,7 @@ def test_costs(self, cost):
             # Test option setting
             cost.set_fail_gradient(1)
 
-        if isinstance(cost, (pybop.SumSquaredError, pybop.ProbabilityCost)):
+        if isinstance(cost, pybop.SumSquaredError):
             e, de = cost.evaluateS1([0.5])
 
             assert type(e) == np.float64

tests/unit/test_likelihoods.py

@@ -75,12 +75,6 @@ def test_base_likelihood_init(self, problem):
         assert likelihood._n_parameters == 1
         assert np.array_equal(likelihood._target, problem._target)
 
-    @pytest.mark.unit
-    def test_base_likelihood_call_raises_not_implemented_error(self, problem):
-        likelihood = pybop.BaseLikelihood(problem)
-        with pytest.raises(NotImplementedError):
-            likelihood(np.array([0.5, 0.5]))
-
     @pytest.mark.unit
     def test_base_likelihood_set_get_sigma(self, problem):
         likelihood = pybop.BaseLikelihood(problem)