From 16b1ad0d950a071e88d8aa8e539a665cd50fb4bd Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sat, 6 Apr 2024 20:07:16 +0100 Subject: [PATCH 1/6] Initial maximum a posteriori implementation & example --- examples/scripts/spm_MAP.py | 71 +++++++++++++++++++++++++++++++++ pybop/__init__.py | 1 + pybop/costs/fitting_costs.py | 76 ++++++++++++++++++++++++++++++++++++ 3 files changed, 148 insertions(+) create mode 100644 examples/scripts/spm_MAP.py diff --git a/examples/scripts/spm_MAP.py b/examples/scripts/spm_MAP.py new file mode 100644 index 000000000..dc34c2def --- /dev/null +++ b/examples/scripts/spm_MAP.py @@ -0,0 +1,71 @@ +import numpy as np + +import pybop + +# Define model +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +model = pybop.lithium_ion.SPM(parameter_set=parameter_set) + +# Fitting parameters +parameters = [ + pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.05), + bounds=[0.5, 0.8], + ), + pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Gaussian(0.48, 0.05), + bounds=[0.4, 0.7], + ), +] + +# Set initial parameter values +parameter_set.update( + { + "Negative electrode active material volume fraction": 0.63, + "Positive electrode active material volume fraction": 0.51, + } +) +# Generate data +sigma = 0.005 +t_eval = np.arange(0, 900, 2) +values = model.predict(t_eval=t_eval) +corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval)) + +# Form dataset +dataset = pybop.Dataset( + { + "Time [s]": t_eval, + "Current function [A]": values["Current [A]"].data, + "Voltage [V]": corrupt_values, + } +) + +# Generate problem, cost function, and optimisation class +problem = pybop.FittingProblem(model, parameters, dataset) +cost = pybop.MAP(pybop.GaussianLogLikelihoodKnownSigma, problem, sigma=[0.01, 0.01]) +optim = pybop.Optimisation(cost, optimiser=pybop.CMAES) +optim.set_max_unchanged_iterations(20) +optim.set_min_iterations(20) +optim.set_max_iterations(100) + +# Run the optimisation +x, final_cost = optim.run() +print("Estimated parameters:", x) + +# Plot the timeseries output +pybop.quick_plot(problem, parameter_values=x[0:2], title="Optimised Comparison") + +# Plot convergence +pybop.plot_convergence(optim) + +# Plot the parameter traces +pybop.plot_parameters(optim) + +# Plot the cost landscape +pybop.plot2d(cost, steps=15) + +# Plot the cost landscape with optimisation path +bounds = np.array([[0.55, 0.77], [0.48, 0.68]]) +pybop.plot2d(optim, bounds=bounds, steps=15) diff --git a/pybop/__init__.py b/pybop/__init__.py index 5ae5c4d23..bd6d8f9cd 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -42,6 +42,7 @@ RootMeanSquaredError, SumSquaredError, ObserverCost, + MAP, ) from .costs.design_costs import ( DesignCost, diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index 930715576..80538d543 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -1,5 +1,6 @@ import numpy as np +from pybop.costs._likelihoods import BaseLikelihood from pybop.costs.base_cost import BaseCost from pybop.observers.observer import Observer @@ -279,3 +280,78 @@ def evaluateS1(self, x): If an error occurs during the calculation of the cost or gradient. """ raise NotImplementedError + + +class MAP(BaseLikelihood): + """ + Maximum a posteriori cost function. + + Computes the maximum a posteriori cost function, which is the sum of the + negative log likelihood and the log prior. + + Inherits all parameters and attributes from ``BaseLikelihood``. + + """ + + def __init__(self, likelihood, problem, sigma): + super(MAP, self).__init__(problem) + self.likelihood = likelihood + self.sigma = sigma + + def _evaluate(self, x, grad=None): + """ + Calculate the maximum a posteriori 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 maximum a posteriori cost. + """ + log_likelihood = self.likelihood( + problem=self.problem, sigma=self.sigma + ).evaluate(x) + log_prior = sum( + param.prior.logpdf(x_i) for x_i, param in zip(x, self.problem.parameters) + ) + + posterior = log_likelihood + log_prior + return posterior + + def _evaluateS1(self, x): + """ + Compute the maximum a posteriori with respect to the parameters. + The method passes the likelihood gradient to the optimiser without modification. + + 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. + """ + log_likelihood, dl = self.likelihood( + problem=self.problem, sigma=self.sigma + ).evaluateS1(x) + log_prior = sum( + param.prior.logpdf(x_i) for x_i, param in zip(x, self.problem.parameters) + ) + + posterior = log_likelihood + log_prior + return posterior, dl From 097864b9ad568b7956e6647134bbac9c02242779 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sun, 7 Apr 2024 13:34:17 +0100 Subject: [PATCH 2/6] Use parameter.prior sigma for MAP if not provided, add tests --- examples/scripts/spm_MAP.py | 2 +- pybop/costs/fitting_costs.py | 10 +++++++--- tests/unit/test_cost.py | 11 +++++++++-- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/examples/scripts/spm_MAP.py b/examples/scripts/spm_MAP.py index dc34c2def..cc22c315c 100644 --- a/examples/scripts/spm_MAP.py +++ b/examples/scripts/spm_MAP.py @@ -44,7 +44,7 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset) -cost = pybop.MAP(pybop.GaussianLogLikelihoodKnownSigma, problem, sigma=[0.01, 0.01]) +cost = pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma) optim = pybop.Optimisation(cost, optimiser=pybop.CMAES) optim.set_max_unchanged_iterations(20) optim.set_min_iterations(20) diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index 80538d543..28581431d 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -293,10 +293,14 @@ class MAP(BaseLikelihood): """ - def __init__(self, likelihood, problem, sigma): + def __init__(self, problem, likelihood, sigma=None): super(MAP, self).__init__(problem) self.likelihood = likelihood - self.sigma = sigma + self.sigma0 = sigma + if self.sigma0 is None: + self.sigma0 = [] + for param in self.problem.parameters: + self.sigma0.append(param.prior.sigma) def _evaluate(self, x, grad=None): """ @@ -316,7 +320,7 @@ def _evaluate(self, x, grad=None): The maximum a posteriori cost. """ log_likelihood = self.likelihood( - problem=self.problem, sigma=self.sigma + problem=self.problem, sigma=self.sigma0 ).evaluate(x) log_prior = sum( param.prior.logpdf(x_i) for x_i, param in zip(x, self.problem.parameters) diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index 9dd772c89..a4c65ecf9 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -73,12 +73,15 @@ def problem(self, model, parameters, dataset, signal, x0, request): pybop.RootMeanSquaredError, pybop.SumSquaredError, pybop.ObserverCost, + pybop.MAP, ] ) 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 in [pybop.ObserverCost]: inputs = {p.name: problem.x0[i] for i, p in enumerate(problem.parameters)} state = problem._model.reinit(inputs) @@ -121,8 +124,12 @@ def test_design_base(self, problem): @pytest.mark.unit def test_costs(self, cost): - higher_cost = cost([0.55]) - lower_cost = cost([0.52]) + if isinstance(cost, pybop.BaseLikelihood): + higher_cost = cost([0.52]) + lower_cost = cost([0.55]) + else: + higher_cost = cost([0.55]) + lower_cost = cost([0.52]) assert higher_cost > lower_cost or ( higher_cost == lower_cost and higher_cost == np.inf ) From 69e545780081a8cca6f7138a6ad49be3356b609e Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sun, 7 Apr 2024 14:12:41 +0100 Subject: [PATCH 3/6] Updt changelog, construct likelihoods on init, add cost to integration tests --- CHANGELOG.md | 1 + pybop/costs/fitting_costs.py | 10 +++------- tests/integration/test_parameterisations.py | 11 ++++++++++- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1453728d5..8228cd079 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#275](https://github.com/pybop-team/PyBOP/pull/275) - Adds Maximum a Posteriori (MAP) cost function with corresponding tests. - [#241](https://github.com/pybop-team/PyBOP/pull/241) - Adds experimental circuit model fitting notebook with LG M50 data. - [#268](https://github.com/pybop-team/PyBOP/pull/268) - Fixes the GitHub Release artifact uploads, allowing verification of codesigned binaries and source distributions via `sigstore-python`. diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index 28581431d..33e729fde 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -295,12 +295,12 @@ class MAP(BaseLikelihood): def __init__(self, problem, likelihood, sigma=None): super(MAP, self).__init__(problem) - self.likelihood = likelihood self.sigma0 = sigma if self.sigma0 is None: self.sigma0 = [] for param in self.problem.parameters: self.sigma0.append(param.prior.sigma) + self.likelihood = likelihood(problem=self.problem, sigma=self.sigma0) def _evaluate(self, x, grad=None): """ @@ -319,9 +319,7 @@ def _evaluate(self, x, grad=None): float The maximum a posteriori cost. """ - log_likelihood = self.likelihood( - problem=self.problem, sigma=self.sigma0 - ).evaluate(x) + log_likelihood = self.likelihood.evaluate(x) log_prior = sum( param.prior.logpdf(x_i) for x_i, param in zip(x, self.problem.parameters) ) @@ -350,9 +348,7 @@ def _evaluateS1(self, x): ValueError If an error occurs during the calculation of the cost or gradient. """ - log_likelihood, dl = self.likelihood( - problem=self.problem, sigma=self.sigma - ).evaluateS1(x) + log_likelihood, dl = self.likelihood.evaluateS1(x) log_prior = sum( param.prior.logpdf(x_i) for x_i, param in zip(x, self.problem.parameters) ) diff --git a/tests/integration/test_parameterisations.py b/tests/integration/test_parameterisations.py index 4da45fd4d..1f7d8b28c 100644 --- a/tests/integration/test_parameterisations.py +++ b/tests/integration/test_parameterisations.py @@ -44,6 +44,7 @@ def init_soc(self, request): pybop.GaussianLogLikelihoodKnownSigma, pybop.RootMeanSquaredError, pybop.SumSquaredError, + pybop.MAP, ] ) def cost_class(self, request): @@ -72,6 +73,10 @@ def spm_costs(self, model, parameters, cost_class, init_soc): ) if cost_class in [pybop.GaussianLogLikelihoodKnownSigma]: return cost_class(problem, sigma=[0.03, 0.03]) + elif cost_class in [pybop.MAP]: + return cost_class( + problem, pybop.GaussianLogLikelihoodKnownSigma, sigma=[0.03, 0.03] + ) else: return cost_class(problem) @@ -125,7 +130,9 @@ def test_spm_optimisers(self, optimiser, spm_costs): assert parameterisation._max_iterations == 125 elif optimiser in [pybop.GradientDescent]: - if isinstance(spm_costs, pybop.GaussianLogLikelihoodKnownSigma): + if isinstance( + spm_costs, (pybop.GaussianLogLikelihoodKnownSigma, pybop.MAP) + ): parameterisation.optimiser.set_learning_rate(1.8e-5) parameterisation.set_min_iterations(150) else: @@ -177,6 +184,8 @@ def spm_two_signal_cost(self, parameters, model, cost_class): if cost_class in [pybop.GaussianLogLikelihoodKnownSigma]: return cost_class(problem, sigma=[0.05, 0.05]) + elif cost_class in [pybop.MAP]: + return cost_class(problem, pybop.GaussianLogLikelihoodKnownSigma) else: return cost_class(problem) From 9c1b714c8c97b439d98eac6034da3754b9af885b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Mon, 8 Apr 2024 10:00:48 +0100 Subject: [PATCH 4/6] Update codecov options to debug upload failure --- .github/workflows/test_on_push.yaml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test_on_push.yaml b/.github/workflows/test_on_push.yaml index af98e8ad9..f4de4fe8b 100644 --- a/.github/workflows/test_on_push.yaml +++ b/.github/workflows/test_on_push.yaml @@ -144,5 +144,7 @@ jobs: - name: Upload coverage report uses: codecov/codecov-action@v3 - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + with: + fail_ci_if_error: true + verbose: true + token: ${{ secrets.CODECOV_TOKEN }} From 40ac8e1dd66c3bf97525f11c5cf105d2d4afc85d Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Mon, 8 Apr 2024 13:46:53 +0100 Subject: [PATCH 5/6] bump codecov to v4 --- .github/workflows/test_on_push.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test_on_push.yaml b/.github/workflows/test_on_push.yaml index f4de4fe8b..882aa151c 100644 --- a/.github/workflows/test_on_push.yaml +++ b/.github/workflows/test_on_push.yaml @@ -143,7 +143,7 @@ jobs: run: nox -s coverage - name: Upload coverage report - uses: codecov/codecov-action@v3 + uses: codecov/codecov-action@v4 with: fail_ci_if_error: true verbose: true From d28029ce54e47f76dc5036e4d54e8473436de440 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 11 Apr 2024 17:50:16 +0100 Subject: [PATCH 6/6] Add catches on MAP likelihood construction --- pybop/costs/fitting_costs.py | 13 ++++++++++++- tests/unit/test_cost.py | 10 ++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index 33e729fde..0665454b0 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -300,7 +300,18 @@ def __init__(self, problem, likelihood, sigma=None): self.sigma0 = [] for param in self.problem.parameters: self.sigma0.append(param.prior.sigma) - self.likelihood = likelihood(problem=self.problem, sigma=self.sigma0) + + try: + self.likelihood = likelihood(problem=self.problem, sigma=self.sigma0) + except Exception as e: + raise ValueError( + f"An error occurred when constructing the Likelihood class: {e}" + ) + + if hasattr(self, "likelihood") and not isinstance( + self.likelihood, BaseLikelihood + ): + raise ValueError(f"{self.likelihood} must be a subclass of BaseLikelihood") def _evaluate(self, x, grad=None): """ diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index a4c65ecf9..cfde8709e 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -122,6 +122,16 @@ def test_design_base(self, problem): with pytest.raises(NotImplementedError): design_cost([0.5]) + @pytest.mark.unit + def test_MAP(self, problem): + # Incorrect likelihood + with pytest.raises(ValueError): + pybop.MAP(problem, pybop.SumSquaredError) + + # Incorrect construction of likelihood + with pytest.raises(ValueError): + pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma, sigma="string") + @pytest.mark.unit def test_costs(self, cost): if isinstance(cost, pybop.BaseLikelihood):