Adds Monte Carlo Samplers

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#6](https://github.com/pybop-team/PyBOP/issues/6) - Adds Monte Carlo functionality, with methods based on Pints' algorithms. A base class is added `BaseSampler`, in addition to `PintsBaseSampler`.
 - [#405](https://github.com/pybop-team/PyBOP/pull/405) - Adds frequency-domain based EIS prediction methods via `model.simulateEIS` and updates to `problem.evaluate` with examples and tests.
 - [#460](https://github.com/pybop-team/PyBOP/pull/460) - Notebook example files added for ECM and folder structure updated.
 - [#450](https://github.com/pybop-team/PyBOP/pull/450) - Adds support for IDAKLU with output variables, and corresponding examples, tests.
@@ -34,7 +35,6 @@
 
 ## Bug Fixes
 
-
 ## Breaking Changes
 
 # [v24.6](https://github.com/pybop-team/PyBOP/tree/v24.6) - 2024-07-08

examples/scripts/mcmc_example.py

@@ -0,0 +1,99 @@
+import numpy as np
+import plotly.graph_objects as go
+import pybamm
+
+import pybop
+
+# Parameter set and model definition
+solver = pybamm.IDAKLUSolver()
+parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+parameter_set.update(
+    {
+        "Negative electrode active material volume fraction": 0.63,
+        "Positive electrode active material volume fraction": 0.71,
+    }
+)
+synth_model = pybop.lithium_ion.DFN(parameter_set=parameter_set, solver=solver)
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Negative electrode active material volume fraction",
+        prior=pybop.Gaussian(0.68, 0.02),
+        transformation=pybop.LogTransformation(),
+    ),
+    pybop.Parameter(
+        "Positive electrode active material volume fraction",
+        prior=pybop.Gaussian(0.68, 0.02),
+        transformation=pybop.LogTransformation(),
+    ),
+)
+
+# Generate data
+init_soc = 1.0
+sigma = 0.001
+experiment = pybop.Experiment(
+    [
+        (
+            "Discharge at 0.5C until 3.5V (10 second period)",
+            "Charge at 0.5C until 4.0V (10 second period)",
+        ),
+    ]
+    # * 2
+)
+values = synth_model.predict(initial_state={"Initial SoC": 1.0}, experiment=experiment)
+
+
+def noise(sigma):
+    return np.random.normal(0, sigma, len(values["Voltage [V]"].data))
+
+
+# Form dataset
+dataset = pybop.Dataset(
+    {
+        "Time [s]": values["Time [s]"].data,
+        "Current function [A]": values["Current [A]"].data,
+        "Voltage [V]": values["Voltage [V]"].data + noise(sigma),
+        "Bulk open-circuit voltage [V]": values["Bulk open-circuit voltage [V]"].data
+        + noise(sigma),
+    }
+)
+
+model = pybop.lithium_ion.SPM(parameter_set=parameter_set, solver=pybamm.IDAKLUSolver())
+model.build(initial_state={"Initial SoC": 1.0})
+signal = ["Voltage [V]", "Bulk open-circuit voltage [V]"]
+
+# Generate problem, likelihood, and sampler
+problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+likelihood = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=0.002)
+prior1 = pybop.Gaussian(0.59, 0.05)
+prior2 = pybop.Gaussian(0.65, 0.05)
+composed_prior = pybop.JointLogPrior(prior1, prior2)
+posterior = pybop.LogPosterior(likelihood, composed_prior)
+
+optim = pybop.DREAM(
+    posterior,
+    chains=3,
+    max_iterations=300,
+    burn_in=100,
+    verbose=True,
+    # parallel=True,  # uncomment to enable parallelisation (MacOS/WSL/Linux only)
+)
+result = optim.run()
+
+# Create a histogram
+fig = go.Figure()
+for _i, data in enumerate(result):
+    fig.add_trace(go.Histogram(x=data[:, 0], name="Neg", opacity=0.75))
+    fig.add_trace(go.Histogram(x=data[:, 1], name="Pos", opacity=0.75))
+
+# Update layout for better visualization
+fig.update_layout(
+    title="Posterior distribution of volume fractions",
+    xaxis_title="Value",
+    yaxis_title="Count",
+    barmode="overlay",
+)
+
+# Show the plot
+fig.show()

pybop/__init__.py

@@ -71,7 +71,7 @@
 #
 from .parameters.parameter import Parameter, Parameters
 from .parameters.parameter_set import ParameterSet
-from .parameters.priors import BasePrior, Gaussian, Uniform, Exponential
+from .parameters.priors import BasePrior, Gaussian, Uniform, Exponential, JointLogPrior
 
 #
 # Model classes
@@ -83,15 +83,15 @@
 from .models.base_model import Inputs
 
 #
-# Problem class
+# Problem classes
 #
 from .problems.base_problem import BaseProblem
 from .problems.fitting_problem import FittingProblem
 from .problems.multi_fitting_problem import MultiFittingProblem
 from .problems.design_problem import DesignProblem
 
 #
-# Cost function class
+# Cost classes
 #
 from .costs.base_cost import BaseCost
 from .costs.fitting_costs import (
@@ -110,12 +110,13 @@
     BaseLikelihood,
     GaussianLogLikelihood,
     GaussianLogLikelihoodKnownSigma,
+    LogPosterior,
     MAP,
 )
 from .costs._weighted_cost import WeightedCost
 
 #
-# Optimiser class
+# Optimiser classes
 #
 
 from .optimisers._cuckoo import CuckooSearchImpl
@@ -141,14 +142,32 @@
 )
 from .optimisers.optimisation import Optimisation
 
+#
+# Monte Carlo classes
+#
+from .samplers.base_sampler import BaseSampler
+from .samplers.base_pints_sampler import BasePintsSampler
+from .samplers.pints_samplers import (
+    NUTS, DREAM, AdaptiveCovarianceMCMC,
+    DifferentialEvolutionMCMC, DramACMC,
+    EmceeHammerMCMC,
+    HaarioACMC, HaarioBardenetACMC,
+    HamiltonianMCMC, MALAMCMC,
+    MetropolisRandomWalkMCMC, MonomialGammaHamiltonianMCMC,
+    PopulationMCMC, RaoBlackwellACMC,
+    RelativisticMCMC, SliceDoublingMCMC,
+    SliceRankShrinkingMCMC, SliceStepoutMCMC,
+)
+from .samplers.mcmc_sampler import MCMCSampler
+
 #
 # Observer classes
 #
 from .observers.unscented_kalman import UnscentedKalmanFilterObserver
 from .observers.observer import Observer
 
 #
-# Plotting class
+# Plotting classes
 #
 from .plotting.plotly_manager import PlotlyManager
 from .plotting.quick_plot import StandardPlot, StandardSubplot, plot_trajectories

pybop/costs/_likelihoods.py

@@ -290,3 +290,97 @@ def compute(
         log_likelihood = self.likelihood.compute(y)
         posterior = log_likelihood + log_prior
         return posterior
+
+
+class LogPosterior(BaseCost):
+    """
+    The Log Posterior for a given problem.
+
+    Computes the log posterior which is the sum of the log
+    likelihood and the log prior.
+
+    Inherits all parameters and attributes from ``BaseCost``.
+    """
+
+    def __init__(self, log_likelihood, log_prior=None):
+        super().__init__(problem=log_likelihood.problem)
+
+        # Store the likelihood and prior
+        self._log_likelihood = log_likelihood
+        self._prior = log_prior
+        if self._prior is None:
+            try:
+                self._prior = log_likelihood.problem.parameters.priors()
+            except Exception as e:
+                raise ValueError(
+                    f"An error occurred when constructing the Prior class: {e}"
+                ) from e
+
+    def compute(
+        self,
+        y: dict,
+        dy: np.ndarray = None,
+        calculate_grad: bool = False,
+    ) -> Union[float, tuple[float, np.ndarray]]:
+        """
+        Calculate the posterior 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 posterior cost.
+        """
+        # Verify we have dy if calculate_grad is True
+        self.verify_args(dy, calculate_grad)
+
+        if calculate_grad:
+            log_prior, dp = self._prior.logpdfS1(self.parameters.current_value())
+        else:
+            log_prior = self._prior.logpdf(self.parameters.current_value())
+
+        if not np.isfinite(log_prior).any():
+            return (-np.inf, -self.grad_fail) if calculate_grad else -np.inf
+
+        if calculate_grad:
+            log_likelihood, dl = self._log_likelihood.compute(
+                y, dy, calculate_grad=True
+            )
+
+            posterior = log_likelihood + log_prior
+            total_gradient = dl + dp
+
+            return posterior, total_gradient
+
+        log_likelihood = self._log_likelihood.compute(y)
+        posterior = log_likelihood + log_prior
+        return posterior
+
+    def prior(self):
+        """
+        Return the prior object.
+
+        Returns
+        -------
+        object
+            The prior object.
+        """
+        return self._prior
+
+    def likelihood(self):
+        """
+        Returns the likelihood.
+
+        Returns
+        -------
+        object
+            The likelihood object.
+        """
+        return self._log_likelihood

pybop/optimisers/base_optimiser.py

@@ -115,7 +115,7 @@ def __init__(
         self.set_base_options()
         self._set_up_optimiser()
 
-        # Throw an warning if any options remain
+        # Throw a warning if any options remain
         if self.unset_options:
             warnings.warn(
                 f"Unrecognised keyword arguments: {self.unset_options} will not be used.",

pybop/parameters/parameter.py

@@ -56,6 +56,9 @@ def __init__(
         self.value = initial_value
         self.transformation = transformation
         self.applied_prior_bounds = False
+        self.bounds = None
+        self.lower_bounds = None
+        self.upper_bounds = None
         self.set_bounds(bounds)
         self.margin = 1e-4
 
@@ -162,17 +165,9 @@ def set_bounds(self, bounds=None, boundary_multiplier=6):
         if bounds is not None:
             if bounds[0] >= bounds[1]:
                 raise ValueError("Lower bound must be less than upper bound")
-            elif self.transformation is not None:
-                self.lower_bound = np.ndarray.item(
-                    self.transformation.to_search(bounds[0])
-                )
-                self.upper_bound = np.ndarray.item(
-                    self.transformation.to_search(bounds[1])
-                )
             else:
                 self.lower_bound = bounds[0]
                 self.upper_bound = bounds[1]
-
         elif self.prior is not None:
             self.applied_prior_bounds = True
             self.lower_bound = self.prior.mean - boundary_multiplier * self.prior.sigma
@@ -181,6 +176,13 @@ def set_bounds(self, bounds=None, boundary_multiplier=6):
         else:
             self.bounds = None
             return
+        if self.transformation is not None:
+            self.lower_bound = np.ndarray.item(
+                self.transformation.to_search(self.lower_bound)
+            )
+            self.upper_bound = np.ndarray.item(
+                self.transformation.to_search(self.upper_bound)
+            )
 
         self.bounds = [self.lower_bound, self.upper_bound]
 
@@ -409,6 +411,12 @@ def get_sigma0(self) -> list:
             sigma0 = None
         return sigma0
 
+    def priors(self) -> list:
+        """
+        Return the prior distribution of each parameter.
+        """
+        return [param.prior for param in self.param.values()]
+
     def initial_value(self) -> np.ndarray:
         """
         Return the initial value of each parameter.