From 7d8f2e2612874f079dd8868b67e641b22d6350f0 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 20 Dec 2024 14:43:40 +0000 Subject: [PATCH] Add Hessian check for 2-parameter problems (#363) * Add Hessian check * Update test_classification.py * Update to implicitly concatenated strings * Update test_classification * Add tests on insensitivity * Update CHANGELOG.md * Fix typo Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> * Update input to OptimisationResult * Improve insensitivity tests * Improve correlation checks * Update in line with OptimisationResult * Make one test a unit test * Update function name * Check bounds proximity if infinite cost --------- Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- CHANGELOG.md | 1 + pybop/__init__.py | 5 + pybop/_classification.py | 139 ++++++++++++++++++ tests/integration/test_classification.py | 174 +++++++++++++++++++++++ tests/unit/test_classifier.py | 50 +++++++ 5 files changed, 369 insertions(+) create mode 100644 pybop/_classification.py create mode 100644 tests/integration/test_classification.py create mode 100644 tests/unit/test_classifier.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a305f9172..a5c166b45 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#362](https://github.com/pybop-team/PyBOP/issues/362) - Adds the `classify_using_Hessian` functionality to classify the optimised result. - [#584](https://github.com/pybop-team/PyBOP/pull/584) - Adds the `GroupedSPMe` model for parameter identification. - [#571](https://github.com/pybop-team/PyBOP/pull/571) - Adds Multistart functionality to optimisers via initialisation arg `multistart`. - [#582](https://github.com/pybop-team/PyBOP/pull/582) - Fixes `population_size` arg for Pints' based optimisers, reshapes `parameters.rvs` to be parameter instances. diff --git a/pybop/__init__.py b/pybop/__init__.py index 486b52953..d453109e9 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -170,6 +170,11 @@ from .observers.unscented_kalman import UnscentedKalmanFilterObserver from .observers.observer import Observer +# +# Classification classes +# +from ._classification import classify_using_hessian + # # Plotting classes # diff --git a/pybop/_classification.py b/pybop/_classification.py new file mode 100644 index 000000000..0a4c2432e --- /dev/null +++ b/pybop/_classification.py @@ -0,0 +1,139 @@ +from typing import Optional + +import numpy as np + +from pybop import OptimisationResult + + +def classify_using_hessian( + result: OptimisationResult, dx=None, cost_tolerance: Optional[float] = 1e-5 +): + """ + A simple check for parameter correlations based on numerical approximation + of the Hessian matrix at the optimal point using central finite differences. + + Parameters + --------- + result : OptimisationResult + The PyBOP optimisation results. + dx : array-like, optional + An array of small positive values used to check proximity to the parameter + bounds and as the perturbation distance in the finite difference calculations. + cost_tolerance : float, optional + A small positive tolerance used for cost value comparisons (default: 1e-5). + """ + x = result.x + dx = np.asarray(dx) if dx is not None else np.maximum(x, 1e-40) * 1e-2 + final_cost = result.final_cost + cost = result.cost + parameters = cost.parameters + minimising = result.minimising + + n = len(x) + if n != 2 or len(dx) != n: + raise ValueError( + "The function classify_using_hessian currently only works in the case " + "of 2 parameters, and dx must have the same length as x." + ) + + # Get a list of parameter names for use in the output message + names = list(parameters.keys()) + + # Evaluate the cost for a grid of surrounding points + costs = np.zeros((3, 3, 2)) + for i in np.arange(0, 3): + for j in np.arange(0, 3): + if i == j == 1: + costs[1, 1, 0] = final_cost + costs[1, 1, 1] = final_cost + else: + costs[i, j, 0] = cost(x + np.multiply([i - 1, j - 1], dx)) + costs[i, j, 1] = cost(x + np.multiply([i - 1, j - 1], 2 * dx)) + + def check_proximity_to_bounds(parameters, x, dx, names) -> str: + bounds = parameters.get_bounds() + if bounds is not None: + for i, value in enumerate(x): + if value > bounds["upper"][i] - dx[i]: + return f" The result is near the upper bound of {names[i]}." + + if value < bounds["lower"][i] + dx[i]: + return f" The result is near the lower bound of {names[i]}." + return "" + + # Classify the result + if (minimising and np.any(costs < final_cost)) or ( + not minimising and np.any(costs > final_cost) + ): + message = "The optimiser has not converged to a stationary point." + message += check_proximity_to_bounds(parameters, x, dx, names) + + elif not np.all([np.isfinite(cost) for cost in costs]): + message = "Classification cannot proceed due to infinite cost value(s)." + message += check_proximity_to_bounds(parameters, x, dx, names) + + else: + # Estimate the Hessian using second-order accurate central finite differences + # cfd_hessian = np.zeros((2, 2)) + # cfd_hessian[0, 0] = costs[2,1,0] - 2 * costs[1,1,0] + costs[0,1,0] + # cfd_hessian[0, 1] = (costs[2,2,0] - costs[2,0,0] + costs[0,0,0] - costs[0,2,0]) / 4 + # cfd_hessian[1, 0] = cfd_hessian[0, 1] + # cfd_hessian[1, 1] = costs[1,2,0] - 2 * costs[1,1,0] + costs[1,0,0] + + # Estimate the Hessian using fourth-order accurate central finite differences + cfd_hessian = np.zeros((2, 2)) + cfd_hessian[0, 0] = ( + -costs[2, 1, 1] + + 16 * costs[2, 1, 0] + - 30 * costs[1, 1, 0] + + 16 * costs[0, 1, 0] + - costs[0, 1, 1] + ) / 12 + cfd_hessian[0, 1] = ( + -(costs[2, 2, 1] - costs[2, 0, 1] + costs[0, 0, 1] - costs[0, 2, 1]) + + 16 * (costs[2, 2, 0] - costs[2, 0, 0] + costs[0, 0, 0] - costs[0, 2, 0]) + ) / 48 + cfd_hessian[1, 0] = cfd_hessian[0, 1] + cfd_hessian[1, 1] = ( + -costs[1, 2, 1] + + 16 * costs[1, 2, 0] + - 30 * costs[1, 1, 0] + + 16 * costs[1, 0, 0] + - costs[1, 0, 1] + ) / 12 + + # Compute the eigenvalues and sort into ascending order + eigenvalues, eigenvectors = np.linalg.eig(cfd_hessian) + idx = eigenvalues.argsort() + eigenvalues = eigenvalues[idx] + eigenvectors = eigenvectors[:, idx] + + # Classify the result + if np.all(eigenvalues > cost_tolerance): + message = "The optimiser has located a minimum." + elif np.all(eigenvalues < -cost_tolerance): + message = "The optimiser has located a maximum." + elif np.all(np.abs(eigenvalues) > cost_tolerance): + message = "The optimiser has located a saddle point." + elif np.all(np.abs(eigenvalues) < cost_tolerance): + message = f"The cost variation is smaller than the cost tolerance: {cost_tolerance}." + else: + # One eigenvalue is too small to classify with certainty + message = "The cost variation is too small to classify with certainty." + + # Check for parameter correlations + if np.any(np.abs(eigenvalues) > cost_tolerance): + if np.allclose(eigenvectors[0], np.array([1, 0])): + message += f" The cost is insensitive to a change of {dx[0]:.2g} in {names[0]}." + elif np.allclose(eigenvectors[0], np.array([0, 1])): + message += f" The cost is insensitive to a change of {dx[1]:.2g} in {names[1]}." + else: + diagonal_costs = [ + cost(x - np.multiply(eigenvectors[:, 0], dx)), + cost(x + np.multiply(eigenvectors[:, 0], dx)), + ] + if np.allclose(final_cost, diagonal_costs, atol=cost_tolerance, rtol=0): + message += " There may be a correlation between these parameters." + + print(message) + return message diff --git a/tests/integration/test_classification.py b/tests/integration/test_classification.py new file mode 100644 index 000000000..f6fc3ab4a --- /dev/null +++ b/tests/integration/test_classification.py @@ -0,0 +1,174 @@ +import numpy as np +import pytest +from pybamm import Parameter + +import pybop + + +class TestClassification: + """ + A class to test the classification of different optimisation results. + """ + + @pytest.fixture( + params=[ + np.asarray([0.05, 0.05]), + np.asarray([0.1, 0.05]), + np.asarray([0.05, 0.01]), + ] + ) + def parameters(self, request): + ground_truth = request.param + return pybop.Parameters( + pybop.Parameter( + "R0 [Ohm]", + prior=pybop.Gaussian(0.05, 0.01), + bounds=[0.02, 0.08], + true_value=ground_truth[0], + ), + pybop.Parameter( + "R1 [Ohm]", + prior=pybop.Gaussian(0.05, 0.01), + bounds=[0.02, 0.08], + true_value=ground_truth[1], + ), + ) + + @pytest.fixture + def parameter_set(self): + parameter_set = pybop.ParameterSet( + json_path="examples/parameters/initial_ecm_parameters.json" + ) + parameter_set.import_parameters() + parameter_set.params.update({"C1 [F]": 1000}) + return parameter_set + + @pytest.fixture + def model(self, parameter_set, parameters): + parameter_set.params.update(parameters.as_dict(parameters.true_value())) + return pybop.empirical.Thevenin(parameter_set=parameter_set) + + @pytest.fixture + def dataset(self, model): + experiment = pybop.Experiment( + [ + "Discharge at 0.5C for 2 minutes (4 seconds period)", + "Charge at 0.5C for 2 minutes (4 seconds period)", + ] + ) + solution = model.predict(experiment=experiment) + return pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": solution["Voltage [V]"].data, + } + ) + + @pytest.fixture + def problem(self, model, parameters, dataset): + return pybop.FittingProblem(model, parameters, dataset) + + @pytest.mark.integration + def test_classify_using_hessian(self, problem): + cost = pybop.RootMeanSquaredError(problem) + x = cost.parameters.true_value() + bounds = cost.parameters.get_bounds() + x0 = np.clip(x, bounds["lower"], bounds["upper"]) + optim = pybop.Optimisation(cost=cost) + results = pybop.OptimisationResult(x=x0, optim=optim) + + if np.all(x == np.asarray([0.05, 0.05])): + message = pybop.classify_using_hessian(results) + assert message == "The optimiser has located a minimum." + elif np.all(x == np.asarray([0.1, 0.05])): + message = pybop.classify_using_hessian(results) + assert message == ( + "The optimiser has not converged to a stationary point." + " The result is near the upper bound of R0 [Ohm]." + ) + elif np.all(x == np.asarray([0.05, 0.01])): + message = pybop.classify_using_hessian(results) + assert message == ( + "The optimiser has not converged to a stationary point." + " The result is near the lower bound of R1 [Ohm]." + ) + else: + raise Exception(f"Please add a check for these values: {x}") + + if np.all(x == np.asarray([0.05, 0.05])): + cost = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=0.002) + optim = pybop.Optimisation(cost=cost) + results = pybop.OptimisationResult(x=x, optim=optim) + + message = pybop.classify_using_hessian(results) + assert message == "The optimiser has located a maximum." + + # message = pybop.classify_using_hessian(results) + # assert message == "The optimiser has located a saddle point." + + @pytest.mark.integration + def test_insensitive_classify_using_hessian(self, parameter_set): + param_R0_a = pybop.Parameter( + "R0_a [Ohm]", + bounds=[0, 0.002], + true_value=0.001, + ) + param_R0_b = pybop.Parameter( + "R0_b [Ohm]", + bounds=[-1e-4, 1e-4], + true_value=0, + ) + parameter_set.params.update( + {"R0_a [Ohm]": 0.001, "R0_b [Ohm]": 0}, + check_already_exists=False, + ) + parameter_set.params.update( + {"R0 [Ohm]": Parameter("R0_a [Ohm]") + Parameter("R0_b [Ohm]")}, + ) + model = pybop.empirical.Thevenin(parameter_set=parameter_set) + + experiment = pybop.Experiment( + ["Discharge at 0.5C for 2 minutes (4 seconds period)"] + ) + solution = model.predict(experiment=experiment) + dataset = pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": solution["Voltage [V]"].data, + } + ) + + for parameters in [ + pybop.Parameters(param_R0_b, param_R0_a), + pybop.Parameters(param_R0_a, param_R0_b), + ]: + problem = pybop.FittingProblem(model, parameters, dataset) + cost = pybop.SumofPower(problem, p=1) + x = cost.parameters.true_value() + optim = pybop.Optimisation(cost=cost) + results = pybop.OptimisationResult(x=x, optim=optim) + + message = pybop.classify_using_hessian(results) + assert message == ( + "The cost variation is too small to classify with certainty." + " The cost is insensitive to a change of 1e-42 in R0_b [Ohm]." + ) + + message = pybop.classify_using_hessian(results, dx=[0.0001, 0.0001]) + assert message == ( + "The optimiser has located a minimum." + " There may be a correlation between these parameters." + ) + + message = pybop.classify_using_hessian(results, cost_tolerance=1e-2) + assert message == ( + "The cost variation is smaller than the cost tolerance: 0.01." + ) + + message = pybop.classify_using_hessian(results, dx=[1, 1]) + assert message == ( + "Classification cannot proceed due to infinite cost value(s)." + " The result is near the upper bound of R0_a [Ohm]." + ) diff --git a/tests/unit/test_classifier.py b/tests/unit/test_classifier.py new file mode 100644 index 000000000..3079ab429 --- /dev/null +++ b/tests/unit/test_classifier.py @@ -0,0 +1,50 @@ +import numpy as np +import pytest + +import pybop + + +class TestClassifier: + """ + A class to test the classification of different optimisation results. + """ + + @pytest.fixture + def problem(self): + model = pybop.empirical.Thevenin() + experiment = pybop.Experiment( + [ + "Discharge at 0.5C for 2 minutes (4 seconds period)", + "Charge at 0.5C for 2 minutes (4 seconds period)", + ] + ) + solution = model.predict(experiment=experiment) + dataset = pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": solution["Voltage [V]"].data, + } + ) + parameters = pybop.Parameters( + pybop.Parameter( + "R0 [Ohm]", + prior=pybop.Uniform(0.001, 0.1), + bounds=[1e-4, 0.1], + ), + ) + return pybop.FittingProblem(model, parameters, dataset) + + @pytest.mark.unit + def test_classify_using_hessian_invalid(self, problem): + cost = pybop.SumSquaredError(problem) + optim = pybop.Optimisation(cost=cost) + x = np.asarray([0.001]) + results = pybop.OptimisationResult(x=x, optim=optim) + + with pytest.raises( + ValueError, + match="The function classify_using_hessian currently only works" + " in the case of 2 parameters, and dx must have the same length as x.", + ): + pybop.classify_using_hessian(results)