From a18bd959443b090db84f03d1cbccefd51c17caed Mon Sep 17 00:00:00 2001 From: u2370093 Date: Thu, 19 Mar 2026 13:30:10 +0000 Subject: [PATCH 1/5] initial values instead of mean --- pybop/optimisers/ep_bolfi_optimiser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pybop/optimisers/ep_bolfi_optimiser.py b/pybop/optimisers/ep_bolfi_optimiser.py index 42541d819..ce5cb8f14 100644 --- a/pybop/optimisers/ep_bolfi_optimiser.py +++ b/pybop/optimisers/ep_bolfi_optimiser.py @@ -264,7 +264,7 @@ def _set_up_optimiser(self): feature_extractors, fixed_parameters={}, # probably baked into self.problem.model free_parameters={ - name: par.get_mean(transformed=True) + name: par.get_initial_value(transformed=True) for name, par in self.problem.parameters.items() }, initial_covariance=self.problem.parameters.get_covariance(transformed=True), From c7023817fda455d99d73ce918c8e505f7d0bd89f Mon Sep 17 00:00:00 2001 From: u2370093 Date: Thu, 19 Mar 2026 14:29:51 +0000 Subject: [PATCH 2/5] fix mean --- .../battery_parameterisation/bayesian_feature_fitting.py | 7 ++++--- pybop/optimisers/ep_bolfi_optimiser.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py index 743183ba7..358119b9c 100644 --- a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py +++ b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py @@ -23,16 +23,17 @@ # Set multivariate parameters (defined in model space) distribution = pybop.MultivariateLogNormal( - mean_log_x=[np.log(original_D_n), np.log(original_D_p)], + mean_log_x=[ + np.log(0.9 * original_D_n) - 0.5 * np.log(2), + np.log(1.1 * original_D_p) - 0.5 * np.log(2), + ], covariance_log_x=[[np.log(2), 0.0], [0.0, np.log(2)]], ) parameter_values["Negative particle diffusivity [m2.s-1]"] = pybop.Parameter( - initial_value=0.9 * original_D_n, transformation=pybop.LogTransformation(), distribution=pybop.MarginalDistribution(distribution, 0), ) parameter_values["Positive particle diffusivity [m2.s-1]"] = pybop.Parameter( - initial_value=1.1 * original_D_p, transformation=pybop.LogTransformation(), distribution=pybop.MarginalDistribution(distribution, 1), ) diff --git a/pybop/optimisers/ep_bolfi_optimiser.py b/pybop/optimisers/ep_bolfi_optimiser.py index ce5cb8f14..42541d819 100644 --- a/pybop/optimisers/ep_bolfi_optimiser.py +++ b/pybop/optimisers/ep_bolfi_optimiser.py @@ -264,7 +264,7 @@ def _set_up_optimiser(self): feature_extractors, fixed_parameters={}, # probably baked into self.problem.model free_parameters={ - name: par.get_initial_value(transformed=True) + name: par.get_mean(transformed=True) for name, par in self.problem.parameters.items() }, initial_covariance=self.problem.parameters.get_covariance(transformed=True), From 557610808856802f73ba5cb1577c30c00a3fcd51 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 2 Apr 2026 13:14:03 +0100 Subject: [PATCH 3/5] Update cost failures --- pybop/costs/base_cost.py | 4 +--- pybop/costs/feature_distances.py | 21 +++++++++++++-------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index ee4a221d5..c140c84f8 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -183,9 +183,7 @@ def evaluate( l = self.parameters.distribution.logpdf(input_values) if not np.isfinite(l).any(): - return self.failure( - inputs=inputs, calculate_sensitivities=calculate_sensitivities - ) + return self.failure(self.parameters.names, calculate_sensitivities) if calculate_sensitivities: return l, dl diff --git a/pybop/costs/feature_distances.py b/pybop/costs/feature_distances.py index 918e6f486..fb6716bb8 100644 --- a/pybop/costs/feature_distances.py +++ b/pybop/costs/feature_distances.py @@ -90,14 +90,19 @@ def _fit(self, t, y): involves applying the fit function to data and comparing to identity. """ t = t - t[0] - fit_guess = self._fit_guess(t, y) - return self._feature_selection( - minimize( - lambda x: np.sum((t - self._inverse_fit_function(y, *x)) ** 2) ** 0.5, - x0=fit_guess, - method="trust-constr", - ).x - ) + try: + fit_guess = self._fit_guess(t, y) + return self._feature_selection( + minimize( + lambda x: ( + np.sum((t - self._inverse_fit_function(y, *x)) ** 2) ** 0.5 + ), + x0=fit_guess, + method="trust-constr", + ).x + ) + except IndexError: + return self.failure(self.parameters.names, calculate_sensitivities=False) def __call__( self, From c4a8f71e47e3b537222419a7c3029e41699bde23 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 2 Apr 2026 14:11:38 +0100 Subject: [PATCH 4/5] Update EP-BOLFI posterior --- .../bayesian_feature_fitting.py | 13 +- pybop/optimisers/ep_bolfi_optimiser.py | 160 +++++++++--------- 2 files changed, 83 insertions(+), 90 deletions(-) diff --git a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py index 358119b9c..bb59d52a7 100644 --- a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py +++ b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py @@ -1,7 +1,6 @@ import numpy as np import pybamm from ep_bolfi.models.solversetup import spectral_mesh_pts_and_method -from pybamm import CasadiSolver, Experiment, print_citations import pybop @@ -23,10 +22,7 @@ # Set multivariate parameters (defined in model space) distribution = pybop.MultivariateLogNormal( - mean_log_x=[ - np.log(0.9 * original_D_n) - 0.5 * np.log(2), - np.log(1.1 * original_D_p) - 0.5 * np.log(2), - ], + mean_log_x=[np.log(0.9 * original_D_n), np.log(1.1 * original_D_p)], covariance_log_x=[[np.log(2), 0.0], [0.0, np.log(2)]], ) parameter_values["Negative particle diffusivity [m2.s-1]"] = pybop.Parameter( @@ -43,13 +39,13 @@ simulator = pybop.pybamm.Simulator( model=model, parameter_values=parameter_values, - protocol=Experiment( + protocol=pybamm.Experiment( [ "Discharge at 1.0 C for 15 minutes (1 second period)", "Rest for 15 minutes (1 second period)", ] ), - solver=CasadiSolver( + solver=pybamm.CasadiSolver( rtol=1e-5, atol=1e-5, root_tol=1e-3, @@ -118,6 +114,7 @@ ) optim = pybop.EP_BOLFI(problem, options=options) result = optim.run() + print("True values:", [original_D_n, original_D_p]) # Plot the optimisation result result.plot_convergence(yaxis={"type": "log"}) @@ -127,4 +124,4 @@ fig = result.plot_predictive(show=False) fig[0].show() - print_citations() + pybamm.print_citations() diff --git a/pybop/optimisers/ep_bolfi_optimiser.py b/pybop/optimisers/ep_bolfi_optimiser.py index 42541d819..70dd6a51e 100644 --- a/pybop/optimisers/ep_bolfi_optimiser.py +++ b/pybop/optimisers/ep_bolfi_optimiser.py @@ -8,27 +8,25 @@ import numpy as np from pybamm import citations -import pybop from pybop import plot from pybop._logging import Logger -from pybop.optimisers.base_optimiser import BaseOptimiser, OptimisationResult -from pybop.parameters.multivariate_distributions import MultivariateGaussian -from pybop.parameters.parameter import Parameters - - -def ep_bolfi_problem_processing(y, problem): - if isinstance(y, dict): - evaluation = problem.cost(y[problem.simulator.output_variables[0]]) - else: - evaluation = problem.cost(y) - if isinstance(evaluation, pybop.costs.evaluation.Evaluation): - return [evaluation.values] - else: - return [evaluation] +from pybop.optimisers.base_optimiser import ( + BaseOptimiser, + OptimisationResult, + OptimiserOptions, +) +from pybop.parameters.multivariate_distributions import ( + MarginalDistribution, + MultivariateGaussian, + MultivariateLogNormal, +) +from pybop.parameters.parameter import Parameter, Parameters +from pybop.problems.meta_problem import MetaProblem +from pybop.problems.problem import Problem @dataclass -class EPBOLFIOptions(pybop.OptimiserOptions): +class EPBOLFIOptions(OptimiserOptions): """ A class to hold EP-BOLFI options for the optimisation process. @@ -187,17 +185,17 @@ class EP_BOLFI(BaseOptimiser): objects, but will be converted to an ep_bolfi.EP_BOLFI instance upon instantiation of this class. To change attributes, re-init. - Only compatible with MultivariateParameters with a - MultivariateGaussian distribution. + Only compatible with MultivariateParameters with a MultivariateGaussian + distribution, or a MultivariateLogNormal with a LogTransformation. """ def __init__( self, - problem: pybop.Problem, + problem: Problem, options: EPBOLFIOptions | None = None, ): - if type(problem) is not pybop.MetaProblem: - problem = pybop.MetaProblem(problem) + if type(problem) is not MetaProblem: + problem = MetaProblem(problem) super().__init__(problem, options) # citations.register("""@article{ # Minka2013, @@ -243,26 +241,25 @@ def __init__( def _set_up_optimiser(self): import ep_bolfi - # Use the first output variable to pass to EP-BOLFI; define separate simulators - # for multiple output variables. + # Define separate simulators for multiple target variables. simulators = [ - lambda inputs, sim=problem.simulator: ( - sim.solve(inputs)[sim.output_variables[0]].data + lambda inputs, problem=problem: ( + problem.simulate(inputs=inputs)[problem.target[0]].data ) for problem in self.problem.problems ] experimental_datasets = [ - problem.target_data for problem in self.problem.problems + problem.target_data[problem.target[0]] for problem in self.problem.problems ] feature_extractors = [ - lambda y, prob=problem: ep_bolfi_problem_processing(y, prob) + lambda y, problem=problem: [problem.cost(y=y).values] for problem in self.problem.problems ] self.optimiser = ep_bolfi.EP_BOLFI( simulators, experimental_datasets, feature_extractors, - fixed_parameters={}, # probably baked into self.problem.model + fixed_parameters={}, # probably baked into each problem.simulator free_parameters={ name: par.get_mean(transformed=True) for name, par in self.problem.parameters.items() @@ -282,11 +279,8 @@ def _set_up_optimiser(self): display_current_feature=None, # ToDo: costs with names fixed_parameter_order=list(enumerate(self.problem.parameters.keys())), ) - self._logger = Logger( - minimising=True, - verbose=self.verbose, - verbose_print_rate=self.verbose_print_rate, - ) + assert self.problem.minimising is True + self._logger = Logger(minimising=self.problem.minimising, verbose=False) def _run(self) -> "BayesianOptimisationResult": verbose_log_target = stdout if self._options.verbose else None @@ -305,31 +299,33 @@ def _run(self) -> "BayesianOptimisationResult": + self._options.bolfi_optimally_acquired_samples ) self.bolfi_posterior = self.optimiser.run( - self._options.bolfi_initial_sobol_samples, - total_samples, - self._options.bolfi_posterior_effective_sample_size, - self._options.ep_iterations, - self._options.ep_stepwise_dampener, - self._options.ep_total_dampening, - -1, # ep_dampener_reduction_steps; better re-init with another dampening factor - self._options.posterior_gelman_rubin_threshold, - self._options.posterior_ess_ratio_threshold_resampling, - self._options.posterior_ess_ratio_threshold_evaluation_at_centre, - self._options.posterior_ess_ratio_threshold_skip_feature, - self._options.max_posterior_sampling_retries, - self._options.posterior_actual_sample_size_increase, - self._options.posterior_model_resample_size_increase, - 4, # independent_mcmc_chains; 4 generally works well - self._options.ep_randomise_feature_order, - False, # normalize_features; does not work when features assume 0, normalise within PyBOP - False, # show_trials; use the PyBOP visualization tools instead - self._options.verbose, - self._options.seed, + bolfi_initial_evidence=self._options.bolfi_initial_sobol_samples, + bolfi_total_evidence=total_samples, + bolfi_posterior_samples=self._options.bolfi_posterior_effective_sample_size, + ep_iterations=self._options.ep_iterations, + ep_dampener=self._options.ep_stepwise_dampener, + final_dampening=self._options.ep_total_dampening, + ep_dampener_reduction_steps=-1, # better re-init with another dampening factor + gelman_rubin_threshold=self._options.posterior_gelman_rubin_threshold, + ess_ratio_resample=self._options.posterior_ess_ratio_threshold_resampling, + ess_ratio_sampling_from_zero=self._options.posterior_ess_ratio_threshold_evaluation_at_centre, + ess_ratio_abort=self._options.posterior_ess_ratio_threshold_skip_feature, + max_heuristic_steps=self._options.max_posterior_sampling_retries, + posterior_sampling_increase=self._options.posterior_actual_sample_size_increase, + model_resampling_increase=self._options.posterior_model_resample_size_increase, + independent_mcmc_chains=4, # 4 generally works well + scramble_ep_feature_order=self._options.ep_randomise_feature_order, + normalize_features=False, # does not work when features assume 0, normalise within PyBOP + show_trials=False, # use the PyBOP visualization tools instead + verbose=self._options.verbose, + seed=self._options.seed, ) end = time.time() ep_bolfi_result = json.loads( self.optimiser.result_to_json(seed=self._options.seed) ) + + # Get the optimiser log ep_bolfi_log = json.loads(self.optimiser.log_to_json()) x_list = np.array(list(ep_bolfi_log["tried parameters"].values())).T # Collect all features into one cost. Note: they are logarithms, @@ -371,46 +367,46 @@ def _run(self) -> "BayesianOptimisationResult": ] self._logger.x_search_best = x_search_best_over_time[-1] self._logger.cost_best = cost_best[-1] - model_mean_dict = { - key: value[0] - for key, value in ep_bolfi_result["inferred parameters"].items() - } - model_mean_array = np.array(list(model_mean_dict.values())) - search_mean_array = [ - par.transformation.to_search(entry)[0] - for entry, par in zip( - model_mean_array, - self.problem.parameters.values(), - strict=False, - ) - ] + + # Get the mean and the 95% confidence error bounds + model_mean = np.array( + [val[0] for val in ep_bolfi_result["inferred parameters"].values()] + ) lower_bounds = np.array( [bounds[0][0] for bounds in ep_bolfi_result["error bounds"].values()] ) upper_bounds = np.array( [bounds[1][0] for bounds in ep_bolfi_result["error bounds"].values()] ) - # The re-use of `parameters` makes transformations easily usable. - posterior = deepcopy(self.problem.parameters) - posterior._distribution = MultivariateGaussian( # noqa: SLF001 - search_mean_array, np.array(ep_bolfi_result["covariance"]) - ) - self._logger.iteration = { - "EP iterations": self._options.ep_iterations, - "total feature iterations": self._options.ep_iterations - * len(self.problem.problems), - } - self._logger.evaluations = { - "model evaluations": len( - list(ep_bolfi_log["tried parameters"].values())[0] - ), - # "surrogate evaluations" are not directly accessible + + # Create the posterior distribution, using the existing parameter transformations + n = len(self.problem.parameters) + if isinstance(self.problem.parameters.distribution, MultivariateLogNormal): + covariance_log_x = np.array(ep_bolfi_result["covariance"]) + mean_log_x = np.zeros(n) + for i in range(n): + mean_log_x[i] = np.log(model_mean[i]) - 0.5 * covariance_log_x[i, i] + posterior_distribution = MultivariateLogNormal( + mean_log_x=mean_log_x, covariance_log_x=covariance_log_x + ) + else: + posterior_distribution = MultivariateGaussian( + mean=model_mean, covariance=np.array(ep_bolfi_result["covariance"]) + ) + posterior_parameters = { + key: Parameter( + distribution=MarginalDistribution(posterior_distribution, i), + initial_value=p.initial_value, + transformation=p.transformation, + ) + for i, (key, p) in enumerate(self.problem.parameters.items()) } + return BayesianOptimisationResult( optim=self, time=end - start, method_name="EP-BOLFI", - posterior=posterior, + posterior=Parameters(posterior_parameters), lower_bounds=lower_bounds, upper_bounds=upper_bounds, ) From a8665e14dda8fe44b08311fdd515d9c3541668b1 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 2 Apr 2026 14:39:37 +0100 Subject: [PATCH 5/5] EP-BOLFI only works with feature distances --- tests/unit/test_optimisation.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index 7d21f414e..d57a1fe9a 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -141,9 +141,8 @@ def multivariate_simulator(self, model, dataset): @pytest.fixture def multivariate_problem(self, multivariate_simulator, dataset): - cost = pybop.SumSquaredError(dataset) - problem = pybop.Problem(multivariate_simulator, cost) - return problem + cost = pybop.SquareRootFeatureDistance(dataset, feature="offset") + return pybop.Problem(multivariate_simulator, cost) @pytest.fixture def gitt_like_problem(self, multivariate_simulator, dataset):