diff --git a/examples/notebooks/comparing_cost_functions.ipynb b/examples/notebooks/comparing_cost_functions.ipynb index b3db5c192..ce54750cf 100644 --- a/examples/notebooks/comparing_cost_functions.ipynb +++ b/examples/notebooks/comparing_cost_functions.ipynb @@ -15,36 +15,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], + "outputs": [], "source": [ "%pip install --upgrade pip ipywidgets -q\n", "%pip install pybop -q\n", @@ -206,7 +177,7 @@ { "data": { "text/plain": [ - "2.0269652551213878e-10" + "1.1753460077019054e-09" ] }, "execution_count": null, @@ -240,7 +211,7 @@ { "data": { "text/plain": [ - "2.0269652551213878e-10" + "1.1753460077019054e-09" ] }, "execution_count": null, @@ -269,13 +240,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "[8.19241247e-05 5.09035035e-06]\n" + "[7.60957550e-05 5.48691392e-06]\n" ] }, { "data": { "text/plain": [ - "0.030472774516224807" + "0.014466627355628724" ] }, "execution_count": null, diff --git a/examples/notebooks/cost-compute-methods.ipynb b/examples/notebooks/cost-compute-methods.ipynb index f795f95c8..32d8099f1 100644 --- a/examples/notebooks/cost-compute-methods.ipynb +++ b/examples/notebooks/cost-compute-methods.ipynb @@ -19,27 +19,6 @@ "id": "1", "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, { "name": "stdout", "output_type": "stream", diff --git a/examples/notebooks/maximum_a_posteriori.ipynb b/examples/notebooks/maximum_a_posteriori.ipynb new file mode 100644 index 000000000..299ec9a79 --- /dev/null +++ b/examples/notebooks/maximum_a_posteriori.ipynb @@ -0,0 +1,878 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "expmkveO04pw" + }, + "source": [ + "## Maximum a Posteriori Parameter Inference\n", + "\n", + "In this notebook, we introduce the Maximum a Posteriori (MAP), which extends Maximum Likelihood Estimation (MLE) by inclusion of a prior $p(\\theta)$ into the cost function. To include this prior information, we construct a Bayesian Posterior with Bayesian's Theorem given as,\n", + "\n", + "$$\n", + "P(\\theta|D) = \\frac{P(D|\\theta)P(\\theta)}{P(D)}\n", + "$$\n", + "\n", + "where, \n", + "$~$ \n", + "$P(\\theta|D)$ represents the posterior and can be read as \"the probability of the parameters $(\\theta)$ given the data $(D)$\", \n", + "$P(D|\\theta)$ is the probability of the data given the parameters, commonly called the likelihood, \n", + "$P(\\theta)$ represents the probability of the parameters commonly called the prior, \n", + "$P(D)$ is the probability of the data and is commonly called the marginal probability. \n", + "\n", + "However, as the marginal probability is commonly difficult to compute and represents a normalisation constant, in the case of MAP this term is forgone and the proportional posterior is optmised instead. This is given as,\n", + "\n", + "$$\n", + "P(\\theta|D) \\propto P(D|\\theta)P(\\theta)\n", + "$$\n", + "\n", + "### Setting up the Environment\n", + "\n", + "Before we begin, we need to ensure that we have all the necessary tools. We will install and import PyBOP as well as upgrade dependencies. We also fix the random seed in order to generate consistent output during development, although this does not need to be done in practice." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "X87NUGPW04py", + "outputId": "0d785b07-7cff-4aeb-e60a-4ff5a669afbf" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Note: you may need to restart the kernel to use updated packages.\n", + "Note: you may need to restart the kernel to use updated packages.\n" + ] + } + ], + "source": [ + "%pip install --upgrade pip ipywidgets -q\n", + "%pip install pybop -q\n", + "\n", + "import time\n", + "\n", + "import numpy as np\n", + "\n", + "import pybop\n", + "\n", + "pybop.PlotlyManager().pio.renderers.default = \"notebook_connected\"\n", + "np.random.seed(8)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "5XU-dMtU04p2" + }, + "source": [ + "### Creating the model\n", + "\n", + "To demonstrate the MAP process, we will first need a forward model and data to run parameter inference on. As we are introducing this as a simple example, we will use the PyBOP forward model with white noise as the reference. This requires defining a parameter set and the model itself." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n", + "model = pybop.lithium_ion.SPM(parameter_set=parameter_set)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Simulating Forward Model\n", + "\n", + "We can now simulate the model using the `model.predict` method. This method is a light wrapper on the `Pybamm.Simulation` class and be used as such. For this example, we use the default current function for the `Chen2020` parameter set (5A) to generate the voltage data. As the goal is to investigate the MAP method, we will generate a range of observations from the forward model. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "sBasxv8U04p3" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of observations in each trajectory: [2, 4, 8, 16, 32, 64, 128, 256]\n" + ] + } + ], + "source": [ + "n = 8 # Number of time-series trajectories\n", + "observations = [\n", + " 2**j for j in range(1, n + 1)\n", + "] # Number of observations in each trajectory\n", + "values = []\n", + "for i in observations:\n", + " t_eval = np.linspace(0, 10, i)\n", + " values.append(model.predict(t_eval=t_eval)) # predict and store\n", + "\n", + "print(f\"Number of observations in each trajectory: {observations}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Adding noise to synthetic voltage data\n", + "\n", + "To make the parameter inference more realistic, we add gaussian noise with zero mean to the data. While this doesn't truly represent the challenge of parameter inference with experimental data, this does ensure the cost landscape curvature isn't perfect. For a more realistic representation of experimental data, a different noise function could be used. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sigma = 0.005\n", + "corrupt_values = values[1][\"Voltage [V]\"].data + np.random.normal(\n", + " 0, sigma, len(values[1][\"Voltage [V]\"].data)\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating the PyBOP dataset\n", + "\n", + "The reference trajectory needs to be included in the optimisation task, which is handed within the `Dataset` class. In this situation, this class is composed of the time, current, and the noisy voltage data; however, if we were performing parameter inference from a different measured signal, such as 'Cell Temperature', that would need to be included." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "zuvGHWID04p_" + }, + "outputs": [], + "source": [ + "dataset = pybop.Dataset(\n", + " {\n", + " \"Time [s]\": values[1][\"Time [s]\"].data,\n", + " \"Current function [A]\": values[1][\"Current [A]\"].data,\n", + " \"Voltage [V]\": corrupt_values,\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ffS3CF_704qA" + }, + "source": [ + "### Constructing Parameters Class\n", + "Next, we need to select the forward model parameters for inference. The PyBOP parameters class composes as many individual PyBOP parameters as the user wants (whether these parameters can be identified is left out of this example). This class requires the parameter name, which must resolve to a parameter within the `ParameterSet` defined above. Additionally, this class can accept an `initial_value` which will be used by the optimiser, as well as bounds. For this example, we will provide a `prior` to the parameter class, which will be used later by the MAP process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "WPCybXIJ04qA" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Default bounds applied based on prior distribution.\n", + "Default bounds applied based on prior distribution.\n" + ] + } + ], + "source": [ + "parameters = pybop.Parameters(\n", + " pybop.Parameter(\n", + " \"Negative particle radius [m]\",\n", + " prior=pybop.Gaussian(4e-6, 1e-6),\n", + " true_value=parameter_set[\"Negative particle radius [m]\"],\n", + " ),\n", + " pybop.Parameter(\n", + " \"Positive particle radius [m]\",\n", + " prior=pybop.Gaussian(5e-6, 1e-6),\n", + " true_value=parameter_set[\"Positive particle radius [m]\"],\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "n4OHa-aF04qA" + }, + "source": [ + "### Setting up the Fitting Problem, Likelihood, and Posterior\n", + "\n", + "With the datasets and parameters defined, we can now construct the `FittingProblem` which composes the model, parameters, and dataset providing a single class with the required information for simulating and assessing the forward model. \n", + "\n", + "As described in the introduction to this notebook, the MAP method uses the non-normalised posterior for optimisation. This is defined in PyBOP as the `LogPosterior` class, and has arguments for the likelihood and prior functions. If a prior is not provided, the parameter priors are used as default. In this example, we will use a `GaussianLogLikelihoodKnownSigma` likelihood function, and the default priors set above. For numerical reasons, we optimise the log of the posterior; however this doesn't affect the final results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "etMzRtx404qA" + }, + "outputs": [], + "source": [ + "problem = pybop.FittingProblem(model, parameters, dataset)\n", + "likelihood = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma)\n", + "posterior = pybop.LogPosterior(likelihood)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "caprp-bV04qB" + }, + "source": [ + "### Plotting the Posterior components\n", + "\n", + "Next, to investigate the individual components of the Posterior. The `LogPosterior` class provides attributes of the prior and likelihood. To investigate the contributions of each to the Posterior we plot the landscapes across a selected parameter range." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "-9OVt0EQ04qB" + }, + "outputs": [ + { + "data": { + "text/html": [ + " \n", + " " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "steps = 10 # Number of discretisation points\n", + "bounds = np.asarray([[1e-6, 9e-6], [1e-6, 9e-6]])\n", + "pybop.plot2d(posterior.prior, bounds=bounds, steps=steps, title=\"Log Prior\")\n", + "pybop.plot2d(posterior.likelihood, bounds=bounds, steps=steps, title=\"Log Likelihood\")\n", + "pybop.plot2d(posterior, bounds=bounds, steps=steps, title=\"Log Posterior\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected, the prior represents a two-dimensional gaussian distribution with a mode at $[4e-6,5e-6]$. The likelihood appears to have a banded shape with ridge of optimal points traversing the higher parameter values. Finally, the Posterior forms the combination of the two. This is the benefit of the MAP process, as it allows for previous information to be included in the parameter inference task. Previous knowledge is encapsulated within the prior function and influences the Posterior, depending on the magnitude of the likelihood function.\n", + "\n", + "To show how this is used within a PyBOP optimisation task, we select the Covariance Matrix Adaptation Evolution Strategy optimiser and run the optimisation. We can then plot the parameter trajectories to investigate how the optimiser performed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inferred Parameters: [5.35425483e-06 5.58415228e-06] in 8.549633979797363 seconds\n", + "True Parameters: [5.86e-06 5.22e-06]\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "optim = pybop.CMAES(\n", + " posterior,\n", + " min_iterations=20,\n", + " max_iterations=100,\n", + ")\n", + "start_time = time.time()\n", + "x, final_cost = optim.run()\n", + "print(f\"Inferred Parameters: {x} in {time.time() - start_time} seconds\")\n", + "print(f\"True Parameters: {parameters.true_value()}\")\n", + "pybop.plot_parameters(optim);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected, the optimisation process returns values close to the true optimal. In this case, as the synthetic data is drawn from the forward model (with additive noise), the underlying structure is close to identical. This gives the optimisation process a very well posed landscape, and as such it finds the correct parameter values. \n", + "\n", + "This is not always the case, especially in inference tasks with low-quality data, sloppy parameters, or poor excitation. In these cases, the prior influence can help 'steer' the optimisation process towards the combination of the likelihood and the user's prior knowledge of the parameters." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "-4pZsDmS04qC" + }, + "source": [ + "## Investigating how the number of observations effects the Posterior\n", + "\n", + "We've seen above that the proportional posterior can be represented from its components, the log-likelihood and log-prior. Next, to better understand when each of these terms can become dominating within the parameter inference problem we vary the number of measurement observations (i.e. the number of samples in the dataset) and inspect the construct posterior.\n", + "\n", + "This is completed below with an increasing series of $2^n$, where $n$ is set above in our original creation of the trajectories. This will give us a visual representation of how the posterior changes with increasing observations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "Hgz8SV4i04qC", + "outputId": "e1e42ae7-5075-4c47-dd68-1b22ecc170f6" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed Sigma for output 1: 0.005\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "for i, val in enumerate(values): # Loop through trajectories\n", + " dataset = pybop.Dataset(\n", + " {\n", + " \"Time [s]\": val[\"Time [s]\"].data,\n", + " \"Current function [A]\": val[\"Current [A]\"].data,\n", + " \"Voltage [V]\": val[\"Voltage [V]\"].data\n", + " + np.random.normal(0, sigma, len(val[\"Voltage [V]\"].data)),\n", + " }\n", + " )\n", + " problem = pybop.FittingProblem(model, parameters, dataset)\n", + " likelihood = pybop.GaussianLogLikelihood(problem, sigma0=sigma)\n", + " posterior = pybop.LogPosterior(likelihood)\n", + " pybop.plot2d(\n", + " posterior,\n", + " bounds=bounds,\n", + " steps=steps,\n", + " title=f\"Posterior with {observations[i]} observations\",\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above contour plots showcase the influence decay of the prior on the posterior as observations are increased. This is expected as prior knowledge on the posterior should be reduced as additional high fidelity observations are obtained. In the case of lower quality, or noisy data, the prior influence can maintain influence.\n", + "\n", + "The influence of the prior and likelihood on the posterior should be investigated during a parameter inference process, as presented above. This provides insight into how much influence prior knowledge has on the optimisation task, whether the likelihood function is well posed with smooth curvature, and finally the overall scale of the posterior." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Concluding Thoughts\n", + "\n", + "This notebook illustrates the process of parameter inference with the Maximum a Posteriori method. This process enables encapsulation of prior knowledge into the optimisation process with influence decay as observations of the system increase. This influence decay has been presented above across observations obtained from the set $({2^n \\mid n \\in \\mathbb{N}})$." + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/notebooks/multi_model_identification.ipynb b/examples/notebooks/multi_model_identification.ipynb index 6248e1fa2..cb838e321 100644 --- a/examples/notebooks/multi_model_identification.ipynb +++ b/examples/notebooks/multi_model_identification.ipynb @@ -25,36 +25,7 @@ "id": "X87NUGPW04py", "outputId": "0d785b07-7cff-4aeb-e60a-4ff5a669afbf" }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip\r\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], + "outputs": [], "source": [ "%pip install --upgrade pip ipywidgets pybamm -q\n", "%pip install pybop -q" diff --git a/examples/scripts/spm_MAP.py b/examples/scripts/spm_MAP.py index add1db448..dd8216674 100644 --- a/examples/scripts/spm_MAP.py +++ b/examples/scripts/spm_MAP.py @@ -60,8 +60,8 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset) -cost = pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=sigma) -optim = pybop.AdamW( +cost = pybop.LogPosterior(pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma)) +optim = pybop.IRPropMin( cost, sigma0=0.05, max_unchanged_iterations=20, diff --git a/pybop/__init__.py b/pybop/__init__.py index 99f30e1a6..6863ccf8d 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -111,7 +111,6 @@ GaussianLogLikelihood, GaussianLogLikelihoodKnownSigma, LogPosterior, - MAP, ) from .costs._weighted_cost import WeightedCost diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 6da737e5f..eb7610ef3 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -1,10 +1,12 @@ -from typing import Union +from typing import Optional, Union import numpy as np +import scipy.stats as stats +import pybop from pybop.costs.base_cost import BaseCost from pybop.parameters.parameter import Parameter, Parameters -from pybop.parameters.priors import JointLogPrior, Uniform +from pybop.parameters.priors import BasePrior, JointLogPrior, Uniform from pybop.problems.base_problem import BaseProblem @@ -102,7 +104,7 @@ class GaussianLogLikelihood(BaseLikelihood): def __init__( self, problem: BaseProblem, - sigma0: Union[float, list[float], list[Parameter]] = 0.002, + sigma0: Union[float, list[float], list[Parameter]] = 0.02, dsigma_scale: float = 1.0, ): super().__init__(problem) @@ -202,113 +204,42 @@ def compute( return e -class MAP(BaseLikelihood): - """ - Maximum a posteriori cost function. - - Computes the maximum a posteriori cost function, which is the sum of the - log likelihood and the log prior. The goal of maximising is achieved by - setting minimising = False in the optimiser settings. - - Inherits all parameters and attributes from ``BaseLikelihood``. - - """ - - def __init__(self, problem, likelihood, sigma0=None, gradient_step=1e-3): - super().__init__(problem) - self.sigma0 = sigma0 - self.gradient_step = gradient_step - if self.sigma0 is None: - self.sigma0 = [] - for param in self.problem.parameters: - self.sigma0.append(param.prior.sigma) - - try: - self.likelihood = likelihood(problem=self.problem, sigma0=self.sigma0) - except Exception as e: - raise ValueError( - f"An error occurred when constructing the Likelihood class: {e}" - ) from e - - if hasattr(self, "likelihood") and not isinstance( - self.likelihood, BaseLikelihood - ): - raise ValueError(f"{self.likelihood} must be a subclass of BaseLikelihood") - - self.parameters = self.likelihood.parameters - self._has_separable_problem = self.likelihood.has_separable_problem - - def compute( - self, - y: dict, - dy: np.ndarray = None, - calculate_grad: bool = False, - ) -> Union[float, tuple[float, np.ndarray]]: - """ - Compute the Maximum a Posteriori for the given parameters. - - If self._has_separable_problem is True, then this method only computes the - likelihood, without calling the problem.evaluate(). Else, problem.evaluate() - is called before computing the likelihood. - - Returns - ------- - float - The maximum a posteriori cost. - """ - # Verify we have dy if calculate_grad is True - self.verify_args(dy, calculate_grad) - - log_prior = sum(param.prior.logpdf(param.value) for param in self.parameters) - - 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.likelihood.compute(y, dy, calculate_grad=True) - - # Compute a finite difference approximation of the gradient of the log prior - delta = self.parameters.initial_value() * self.gradient_step - prior_gradient = [] - - for parameter, step_size in zip(self.problem.parameters, delta): - param_value = parameter.value - - log_prior_upper = parameter.prior.logpdf(param_value * (1 + step_size)) - log_prior_lower = parameter.prior.logpdf(param_value * (1 - step_size)) - - gradient = (log_prior_upper - log_prior_lower) / ( - 2 * step_size * param_value + np.finfo(float).eps - ) - prior_gradient.append(gradient) - - posterior = log_likelihood + log_prior - total_gradient = dl + prior_gradient - - return posterior, total_gradient - - log_likelihood = self.likelihood.compute(y) - posterior = log_likelihood + log_prior - return posterior - - -class LogPosterior(BaseCost): +class LogPosterior(BaseLikelihood): """ The Log Posterior for a given problem. - Computes the log posterior which is the sum of the log + Computes the log posterior which is proportional to the sum of the log likelihood and the log prior. - Inherits all parameters and attributes from ``BaseCost``. + Parameters + ---------- + log_likelihood : BaseLikelihood + The likelihood class of type ``BaseLikelihood``. + log_prior : Optional, Union[pybop.BasePrior, stats.rv_continuous] + The prior class of type ``BasePrior`` or ``stats.rv_continuous``. + If not provided, the prior class will be taken from the parameter priors + constructed in the `pybop.Parameters` class. + gradient_step : float, default: 1e-3 + The step size for the finite-difference gradient calculation + if the ``log_prior`` is not of type ``BasePrior``. """ - def __init__(self, log_likelihood, log_prior=None): + def __init__( + self, + log_likelihood: BaseLikelihood, + log_prior: Optional[Union[pybop.BasePrior, stats.rv_continuous]] = None, + gradient_step: float = 1e-3, + ): super().__init__(problem=log_likelihood.problem) + self.gradient_step = gradient_step # Store the likelihood and prior self._log_likelihood = log_likelihood + self.parameters = self._log_likelihood.parameters + self._has_separable_problem = self._log_likelihood.has_separable_problem + if log_prior is None: - self._prior = JointLogPrior(*log_likelihood.problem.parameters.priors()) + self._prior = JointLogPrior(*self.parameters.priors()) else: self._prior = log_prior @@ -319,26 +250,48 @@ def compute( calculate_grad: bool = False, ) -> Union[float, tuple[float, np.ndarray]]: """ - Calculate the posterior cost for a given set of parameters. + Calculate the posterior cost for a given forward model prediction. 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. + y : dict + The data for which to evaluate the cost. + dy : np.ndarray, optional + The correspond sensitivities in the data. + calculate_grad : bool, optional + Whether to calculate the gradient of the cost function. Returns ------- - float - The posterior cost. + Union[float, Tuple[float, np.ndarray]] + The posterior cost, and optionally the gradient. """ # 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()) + if isinstance(self._prior, BasePrior): + log_prior, dp = self._prior.logpdfS1(self.parameters.current_value()) + else: + # Compute log prior first + log_prior = self._prior.logpdf(self.parameters.current_value()) + + # Compute a finite difference approximation of the gradient of the log prior + delta = self.parameters.initial_value() * self.gradient_step + dp = [] + + for parameter, step_size in zip(self.parameters, delta): + param_value = parameter.value + upper_value = param_value * (1 + step_size) + lower_value = param_value * (1 - step_size) + + log_prior_upper = parameter.prior.logpdf(upper_value) + log_prior_lower = parameter.prior.logpdf(lower_value) + + gradient = (log_prior_upper - log_prior_lower) / ( + 2 * step_size * param_value + np.finfo(float).eps + ) + dp.append(gradient) else: log_prior = self._prior.logpdf(self.parameters.current_value()) @@ -359,24 +312,10 @@ def compute( posterior = log_likelihood + log_prior return posterior - def prior(self): - """ - Return the prior object. - - Returns - ------- - object - The prior object. - """ + @property + def prior(self) -> BasePrior: return self._prior - def likelihood(self): - """ - Returns the likelihood. - - Returns - ------- - object - The likelihood object. - """ + @property + def likelihood(self) -> BaseLikelihood: return self._log_likelihood diff --git a/pybop/costs/_weighted_cost.py b/pybop/costs/_weighted_cost.py index 1572fc408..10d7d5440 100644 --- a/pybop/costs/_weighted_cost.py +++ b/pybop/costs/_weighted_cost.py @@ -60,8 +60,9 @@ def __init__(self, *costs, weights: Optional[list[float]] = None): super().__init__(self.costs[0].problem) else: super().__init__() - for cost in self.costs: - self.parameters.join(cost.parameters) + + for cost in self.costs: + self.parameters.join(cost.parameters) # Weighted costs do not use this functionality self._has_separable_problem = False @@ -90,7 +91,7 @@ def compute( The weighted cost value. """ if self._has_identical_problems: - inputs = self.parameters.as_dict() + inputs = self.problem.parameters.as_dict() if calculate_grad: y, dy = self.problem.evaluateS1(inputs) else: diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 370abc276..f7555a815 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -559,3 +559,18 @@ def verify(self, inputs: Optional[Inputs] = None): raise TypeError( f"Inputs must be a dictionary or numeric. Received {type(inputs)}" ) + + def __repr__(self): + """ + Return a string representation of the Parameters instance. + + Returns + ------- + str + A string including the number of parameters and a summary of each parameter. + """ + param_summary = "\n".join( + f" {name}: prior= {param.prior}, value={param.value}, bounds={param.bounds}" + for name, param in self.param.items() + ) + return f"Parameters({len(self)}):\n{param_summary}" diff --git a/pybop/parameters/priors.py b/pybop/parameters/priors.py index 9e68cf008..5c40c42c1 100644 --- a/pybop/parameters/priors.py +++ b/pybop/parameters/priors.py @@ -255,7 +255,7 @@ def _logpdfS1(self, x): float The value(s) of the first derivative at x. """ - return self(x), -(x - self.loc) * self._multip + return self.__call__(x), -(x - self.loc) * self._multip class Uniform(BasePrior): @@ -417,17 +417,16 @@ def _logpdfS1(self, x: Union[float, np.ndarray]) -> tuple[float, np.ndarray]: f"Input x must have length {self._n_parameters}, got {len(x)}" ) - output = 0 - doutput = np.zeros(self._n_parameters) - index = 0 - - for prior in self._priors: - num_params = 1 - x_subset = x[index : index + num_params] - p, dp = prior.logpdfS1(x_subset) - output += p - doutput[index : index + num_params] = dp - index += num_params + log_probs = [] + derivatives = [] + + for prior, xi in zip(self._priors, x): + p, dp = prior.logpdfS1(xi) + log_probs.append(p) + derivatives.append(dp) + + output = sum(log_probs) + doutput = np.array(derivatives) return output, doutput diff --git a/pybop/plotting/plot2d.py b/pybop/plotting/plot2d.py index 5ea2a8519..09a49c0ad 100644 --- a/pybop/plotting/plot2d.py +++ b/pybop/plotting/plot2d.py @@ -4,7 +4,7 @@ import numpy as np from scipy.interpolate import griddata -from pybop import BaseOptimiser, Optimisation, PlotlyManager +from pybop import BaseCost, BaseOptimiser, Optimisation, PlotlyManager def plot2d( @@ -64,11 +64,11 @@ def plot2d( cost = cost_or_optim plot_optim = False - if len(cost.parameters) < 2: + if isinstance(cost, BaseCost) and len(cost.parameters) < 2: raise ValueError("This cost function takes fewer than 2 parameters.") additional_values = [] - if len(cost.parameters) > 2: + if isinstance(cost, BaseCost) and len(cost.parameters) > 2: warnings.warn( "This cost function requires more than 2 parameters. " "Plotting in 2d with fixed values for the additional parameters.", diff --git a/pybop/samplers/base_sampler.py b/pybop/samplers/base_sampler.py index 1766759bf..a46143e87 100644 --- a/pybop/samplers/base_sampler.py +++ b/pybop/samplers/base_sampler.py @@ -110,7 +110,7 @@ def _end_initial_phase(self): def _initialise_storage(self): self._prior = None if isinstance(self._log_pdf, LogPosterior): - self._prior = self._log_pdf.prior() + self._prior = self._log_pdf.prior # Storage of the received samples self._sampled_logpdf = np.zeros(self._n_chains) diff --git a/tests/integration/test_eis_parameterisation.py b/tests/integration/test_eis_parameterisation.py index 8022b67dc..e76f6d2cf 100644 --- a/tests/integration/test_eis_parameterisation.py +++ b/tests/integration/test_eis_parameterisation.py @@ -61,7 +61,7 @@ def init_soc(self, request): pybop.SumSquaredError, pybop.SumofPower, pybop.Minkowski, - pybop.MAP, + pybop.LogPosterior, ] ) def cost(self, request): @@ -112,13 +112,13 @@ def optim(self, optimiser, model, parameters, cost, init_soc): problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) # Construct the cost - if cost in [pybop.GaussianLogLikelihoodKnownSigma]: + if cost is pybop.GaussianLogLikelihoodKnownSigma: cost = cost(problem, sigma0=self.sigma0) - elif cost in [pybop.GaussianLogLikelihood]: + elif cost is pybop.GaussianLogLikelihood: cost = cost(problem, sigma0=self.sigma0 * 4) # Initial sigma0 guess - elif cost in [pybop.MAP]: + elif cost is pybop.LogPosterior: cost = cost( - problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0 + pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=self.sigma0) ) elif cost in [pybop.SumofPower, pybop.Minkowski]: cost = cost(problem, p=2) @@ -133,14 +133,14 @@ def optim(self, optimiser, model, parameters, cost, init_soc): "max_unchanged_iterations": 35, } - if isinstance(cost, pybop.MAP): + if isinstance(cost, pybop.LogPosterior): for i in cost.parameters.keys(): cost.parameters[i].prior = pybop.Uniform( 0.2, 2.0 ) # Increase range to avoid prior == np.inf # Set sigma0 and create optimiser - sigma0 = 0.05 if isinstance(cost, pybop.MAP) else None + sigma0 = 0.05 if isinstance(cost, pybop.LogPosterior) else None optim = optimiser(sigma0=sigma0, **common_args) return optim diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 1392ca83d..b6996b8c9 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -58,7 +58,7 @@ def init_soc(self, request): pybop.SumSquaredError, pybop.SumofPower, pybop.Minkowski, - pybop.MAP, + pybop.LogPosterior, ] ) def cost(self, request): @@ -103,13 +103,13 @@ def optim(self, optimiser, model, parameters, cost, init_soc): problem = pybop.FittingProblem(model, parameters, dataset) # Construct the cost - if cost in [pybop.GaussianLogLikelihoodKnownSigma]: + if cost is pybop.GaussianLogLikelihoodKnownSigma: cost = cost(problem, sigma0=self.sigma0) - elif cost in [pybop.GaussianLogLikelihood]: + elif cost is pybop.GaussianLogLikelihood: cost = cost(problem, sigma0=self.sigma0 * 4) # Initial sigma0 guess - elif cost in [pybop.MAP]: + elif cost is pybop.LogPosterior: cost = cost( - problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0 + pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=self.sigma0) ) elif cost in [pybop.SumofPower, pybop.Minkowski]: cost = cost(problem, p=2) @@ -124,14 +124,14 @@ def optim(self, optimiser, model, parameters, cost, init_soc): "max_unchanged_iterations": 55, } - if isinstance(cost, pybop.MAP): + if isinstance(cost, pybop.LogPosterior): for i in cost.parameters.keys(): cost.parameters[i].prior = pybop.Uniform( 0.2, 2.0 ) # Increase range to avoid prior == np.inf # Set sigma0 and create optimiser - sigma0 = 0.05 if isinstance(cost, pybop.MAP) else None + sigma0 = 0.05 if isinstance(cost, pybop.LogPosterior) else None optim = optimiser(sigma0=sigma0, **common_args) # AdamW will use lowest sigma0 for learning rate, so allow more iterations @@ -192,13 +192,13 @@ def spm_two_signal_cost(self, parameters, model, cost): signal = ["Voltage [V]", "Bulk open-circuit voltage [V]"] problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) - if cost in [pybop.GaussianLogLikelihoodKnownSigma]: + if cost is pybop.GaussianLogLikelihoodKnownSigma: return cost(problem, sigma0=self.sigma0) - elif cost in [pybop.GaussianLogLikelihood]: + elif cost is pybop.GaussianLogLikelihood: return cost(problem, sigma0=self.sigma0 * 4) # Initial sigma0 guess - elif cost in [pybop.MAP]: + elif cost is pybop.LogPosterior: return cost( - problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0 + pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=self.sigma0) ) else: return cost(problem) diff --git a/tests/integration/test_weighted_cost.py b/tests/integration/test_weighted_cost.py index 2436bd0c3..03d963fc3 100644 --- a/tests/integration/test_weighted_cost.py +++ b/tests/integration/test_weighted_cost.py @@ -55,7 +55,7 @@ def init_soc(self, request): pybop.GaussianLogLikelihoodKnownSigma, pybop.RootMeanSquaredError, pybop.SumSquaredError, - pybop.MAP, + pybop.LogPosterior, ) ] ) @@ -82,12 +82,12 @@ def weighted_fitting_cost(self, model, parameters, cost_class, init_soc): problem = pybop.FittingProblem(model, parameters, dataset) costs = [] for cost in cost_class: - if issubclass(cost, pybop.MAP): + if issubclass(cost, pybop.LogPosterior): costs.append( cost( - problem, - pybop.GaussianLogLikelihoodKnownSigma, - sigma0=self.sigma0, + pybop.GaussianLogLikelihoodKnownSigma( + problem, sigma0=self.sigma0 + ), ) ) elif issubclass(cost, pybop.BaseLikelihood): diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index c17b4a30d..9d11982be 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -75,18 +75,18 @@ def problem(self, model, parameters, dataset, signal, request): pybop.Minkowski, pybop.SumofPower, pybop.ObserverCost, - pybop.MAP, + pybop.LogPosterior, ] ) def cost(self, problem, request): cls = request.param if cls in [pybop.SumSquaredError, pybop.RootMeanSquaredError]: return cls(problem) - elif cls in [pybop.MAP]: - return cls(problem, pybop.GaussianLogLikelihoodKnownSigma) + elif cls is pybop.LogPosterior: + return cls(pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=0.002)) elif cls in [pybop.Minkowski, pybop.SumofPower]: return cls(problem, p=2) - elif cls in [pybop.ObserverCost]: + elif cls is pybop.ObserverCost: inputs = problem.parameters.initial_value() state = problem.model.reinit(inputs) n = len(state) @@ -130,41 +130,6 @@ def compute(self, y, dy=None, calculate_grad: bool = False): with pytest.raises(ValueError, match="Error in cost calculation: Error test."): cost([0.5], calculate_grad=True) - @pytest.mark.unit - def test_MAP(self, problem): - # Incorrect likelihood - with pytest.raises( - ValueError, - match="An error occurred when constructing the Likelihood class:", - ): - pybop.MAP(problem, pybop.SumSquaredError) - - # Incorrect construction of likelihood - with pytest.raises( - ValueError, - match="An error occurred when constructing the Likelihood class: could not convert string to float: 'string'", - ): - pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0="string") - - # Incorrect likelihood - with pytest.raises(ValueError, match="must be a subclass of BaseLikelihood"): - pybop.MAP(problem, self.InvalidLikelihood, sigma0=0.1) - - # Non finite prior - parameter = pybop.Parameter( - "Negative electrode active material volume fraction", - prior=pybop.Uniform(0.55, 0.6), - ) - dataset = pybop.Dataset(data_dictionary=problem.dataset) - problem_non_finite = pybop.FittingProblem( - problem.model, parameter, dataset, signal=problem.signal - ) - likelihood = pybop.MAP( - problem_non_finite, pybop.GaussianLogLikelihoodKnownSigma, sigma0=0.01 - ) - assert not np.isfinite(likelihood([0.7])) - assert not np.isfinite(likelihood([0.7], calculate_grad=True)[0]) - @pytest.mark.unit def test_costs(self, cost): if isinstance(cost, pybop.BaseLikelihood): @@ -190,7 +155,7 @@ def test_costs(self, cost): cost.set_fail_gradient(10) assert cost._de == 10 - if not isinstance(cost, (pybop.ObserverCost, pybop.MAP)): + if not isinstance(cost, (pybop.ObserverCost, pybop.LogPosterior)): e, de = cost([0.5], calculate_grad=True) assert np.isscalar(e) @@ -335,17 +300,17 @@ def test_weighted_fitting_cost(self, noisy_problem): TypeError, match="All costs must be instances of BaseCost.", ): - weighted_cost = pybop.WeightedCost(cost1.problem) + pybop.WeightedCost(cost1.problem) with pytest.raises( ValueError, match="Weights must be numeric values.", ): - weighted_cost = pybop.WeightedCost(cost1, cost2, weights="Invalid string") + pybop.WeightedCost(cost1, cost2, weights="Invalid string") with pytest.raises( ValueError, match="Number of weights must match number of costs.", ): - weighted_cost = pybop.WeightedCost(cost1, cost2, weights=[1]) + pybop.WeightedCost(cost1, cost2, weights=[1]) # Test with identical problems weight = 100 @@ -378,8 +343,8 @@ def test_weighted_fitting_cost(self, noisy_problem): np.testing.assert_allclose(errors_2, errors_3, atol=1e-5) np.testing.assert_allclose(sensitivities_2, sensitivities_3, atol=1e-5) - # Test MAP explicitly - cost4 = pybop.MAP(problem, pybop.GaussianLogLikelihood) + # Test LogPosterior explicitly + cost4 = pybop.LogPosterior(pybop.GaussianLogLikelihood(problem)) 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 diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py index 1a3806aa5..c467b787e 100644 --- a/tests/unit/test_parameters.py +++ b/tests/unit/test_parameters.py @@ -229,3 +229,11 @@ def test_initial_values_without_attributes(self): sample = parameter.initial_value() np.testing.assert_equal(sample, np.array([None])) + + @pytest.mark.unit + def test_parameters_repr(self, parameter): + params = pybop.Parameters(parameter) + assert ( + repr(params) + == "Parameters(1):\n Negative electrode active material volume fraction: prior= Gaussian, loc: 0.6, scale: 0.02, value=0.6, bounds=[0.375, 0.7]" + ) diff --git a/tests/unit/test_posterior.py b/tests/unit/test_posterior.py index 7073b94c3..b43e21ecd 100644 --- a/tests/unit/test_posterior.py +++ b/tests/unit/test_posterior.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import scipy.stats as st import pybop @@ -72,7 +73,7 @@ def test_log_posterior_construction(self, likelihood, prior): assert posterior._prior == prior # Test log posterior construction without parameters - likelihood.problem.parameters.priors = None + likelihood.parameters.priors = None with pytest.raises(TypeError, match="'NoneType' object is not callable"): pybop.LogPosterior(likelihood, log_prior=None) @@ -104,8 +105,8 @@ def test_log_posterior(self, posterior): assert np.allclose(dp, 2.0, atol=2e-2) # Get log likelihood and log prior - likelihood = posterior.likelihood() - prior = posterior.prior() + likelihood = posterior.likelihood + prior = posterior.prior assert likelihood == posterior._log_likelihood assert prior == posterior._prior @@ -119,3 +120,14 @@ def test_log_posterior_inf(self, posterior_uniform_prior): # Test prior np.inf assert not np.isfinite(posterior_uniform_prior([1])) assert not np.isfinite(posterior_uniform_prior([1], calculate_grad=True)[0]) + + @pytest.mark.unit + def test_non_logpdfS1_prior(self, likelihood): + # Scipy distribution + prior = st.norm(0.8, 0.01) + posterior = pybop.LogPosterior(likelihood, log_prior=prior) + p, dp = posterior([0.6], calculate_grad=True) + + # Assert to PyBOP.Gaussian + p2, dp2 = pybop.Gaussian(0.8, 0.01).logpdfS1(0.6) + np.testing.assert_allclose(dp, dp2, atol=2e-3)