initial values instead of mean

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,16 +22,14 @@
 
 # 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), 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(
-    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),
 )
@@ -42,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,
@@ -117,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"})
@@ -126,4 +124,4 @@
     fig = result.plot_predictive(show=False)
     fig[0].show()
 
-    print_citations()
+    pybamm.print_citations()

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

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,

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,
         )

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):