diff --git a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py index d851937d4..7b9ae5516 100644 --- a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py +++ b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py @@ -11,9 +11,27 @@ original_D_n = parameter_values["Negative particle diffusivity [m2.s-1]"] original_D_p = parameter_values["Positive particle diffusivity [m2.s-1]"] -# Put empty Parameter slots as placeholders -parameter_values["Negative particle diffusivity [m2.s-1]"] = pybop.Parameter() -parameter_values["Positive particle diffusivity [m2.s-1]"] = pybop.Parameter() +# Set multivariate parameters +parameter_values["Negative particle diffusivity [m2.s-1]"] = ( + pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + initial_value=0.9 * original_D_n, + bounds=[original_D_n / 2, original_D_n * 2], + transformation=pybop.LogTransformation(), + distribution=pybop.MultivariateGaussian( + [np.log(original_D_n), np.log(original_D_p)], + [[np.log(2), 0.0], [0.0, np.log(2)]], + ), + ) +) +parameter_values["Positive particle diffusivity [m2.s-1]"] = ( + pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + initial_value=1.1 * original_D_p, + bounds=[original_D_p / 2, original_D_p * 2], + transformation=pybop.LogTransformation(), + ) +) # Set up simulator with custom settings submesh_types, var_pts, spatial_methods = spectral_mesh_pts_and_method(10, 10, 10) @@ -58,26 +76,6 @@ } ) -# Override the forced univariate Parameters -simulator.parameters = pybop.MultivariateParameters( - { - "Negative particle diffusivity [m2.s-1]": pybop.Parameter( - initial_value=0.9 * original_D_n, - bounds=[original_D_n / 2, original_D_n * 2], - transformation=pybop.LogTransformation(), - ), - "Positive particle diffusivity [m2.s-1]": pybop.Parameter( - initial_value=1.1 * original_D_p, - bounds=[original_D_p / 2, original_D_p * 2], - transformation=pybop.LogTransformation(), - ), - }, - distribution=pybop.MultivariateGaussian( - [np.log(original_D_n), np.log(original_D_p)], - [[np.log(2), 0.0], [0.0, np.log(2)]], - ), -) - ICI_cost = pybop.SquareRootFeatureDistance( dataset["Time [s]"], dataset["Voltage [V]"], @@ -98,9 +96,6 @@ GITT_problem = pybop.Problem(simulator, GITT_cost) problem = pybop.MetaProblem(ICI_problem, GITT_problem) - # Copy the MultivariateParameters to the meta-problem - problem.parameters = simulator.parameters - # Set up and run the optimiser, increase the number of iterations # and samples to improve accuracy options = pybop.EPBOLFIOptions( diff --git a/pybop/__init__.py b/pybop/__init__.py index 68dab1a03..e2ca82b22 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -47,8 +47,7 @@ # # Parameter classes # -from .parameters.parameter import Parameter, Parameters -from .parameters.multivariate_parameters import MultivariateParameters +from .parameters.parameter import Parameter, Parameters, MultivariateParameter from .parameters.distributions import Distribution, Gaussian, Uniform, Exponential, JointDistribution from .parameters.multivariate_distributions import MultivariateNonparametric, MultivariateUniform, MultivariateGaussian diff --git a/pybop/_result.py b/pybop/_result.py index 5df606f59..894df88c9 100644 --- a/pybop/_result.py +++ b/pybop/_result.py @@ -5,7 +5,7 @@ if TYPE_CHECKING: from pybop import BaseOptimiser, BaseSampler from pybop import Logger, plot -from pybop.parameters.multivariate_parameters import MultivariateParameters +from pybop.parameters.parameter import Parameters class OptimisationResult: @@ -390,7 +390,7 @@ class BayesianOptimisationResult(OptimisationResult): The lower confidence parameter boundaries. upper_bounds: ndarray The upper confidence parameter boundaries. - posterior : MultivariateParameters + posterior : Parameters The probability distribution of the optimisation. maximum_a_posteriori : Inputs or ndarray Complementing the best observed value in `x`, this is the @@ -414,7 +414,7 @@ def __init__( message: str | None = None, lower_bounds: np.ndarray | None = None, upper_bounds: np.ndarray | None = None, - posterior: MultivariateParameters | None = None, + posterior: Parameters | None = None, maximum_a_posteriori: np.ndarray | None = None, log_evidence_mean: float | None = None, log_evidence_variance: float | None = None, diff --git a/pybop/parameters/multivariate_parameters.py b/pybop/parameters/multivariate_parameters.py deleted file mode 100644 index e89fc8cf1..000000000 --- a/pybop/parameters/multivariate_parameters.py +++ /dev/null @@ -1,84 +0,0 @@ -import numpy as np - -from pybop import Parameters - - -class MultivariateParameters(Parameters): - """ - Represents a correlated set of uncertain parameters within the PyBOP - framework. - - This class encapsulates the definition of each of its parameters, - including their names and bounds to ensure the parameter - stays within feasible limits during optimisation or sampling. - - Parameters - ---------- - parameters : pybop.Parameters or dict - The optimisation parameters. - distribution : pybop.BaseMultivariateDistribution - The joint multivariate distribution. - """ - - def __init__(self, parameters: dict | Parameters, distribution=None): - self.distribution = distribution - super().__init__(parameters) - for param in self._parameters.values(): - # Ensure that no individual distributions are mixed with the joint - # one. They may have been used for setting boundaries. - param._distribution = None # noqa: SLF001 - - def pdf(self, x: np.ndarray) -> np.ndarray: - """ - Evaluate the probability density function in parameter space. - Just use the pdf method of the distribution directly for search space. - - Parameters - ---------- - x : np.ndarray - The parameter vector(s) to evaluate the pdf at. - - Returns - ------- - np.ndarray - An array of the pdf values at each parameter vector. - """ - if len(x.shape) == 1: # one-dimensional - x = np.atleast_2d(x) - x = np.asarray([self.transformation.to_search(m) for m in x]) - return self.distribution.pdf(x) - - def sample_from_distribution( - self, n_samples: int = 1, random_state=None, apply_transform: bool = False - ) -> np.ndarray: - return self.rvs(n_samples, random_state, apply_transform) - - def rvs( - self, n_samples: int = 1, random_state=None, apply_transform: bool = False - ) -> np.ndarray: - """ - Draw random samples from the joint parameters distribution. - - Parameters - ---------- - n_samples : int - The number of samples to draw (default: 1). - random_state : int, optional - The random state seed for reproducibility (default: None). - apply_transform : bool - If True, the transformation is applied to the output - (default: False). - - Returns - ------- - array-like - A 2D array of samples drawn from the joint distribution. - """ - samples = self.distribution.rvs(n_samples, random_state=random_state) - if samples.ndim < 2: - samples = np.atleast_2d(samples) - - if apply_transform: - samples = np.asarray([self.transformation.to_model(s) for s in samples]) - - return samples diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 0c18389b9..2659fa6e6 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -11,6 +11,7 @@ from numpy.typing import NDArray from pybop.parameters.distributions import Distribution +from pybop.parameters.multivariate_distributions import BaseMultivariateDistribution from pybop.transformation.base_transformation import Transformation from pybop.transformation.transformations import ( ComposedTransformation, @@ -238,6 +239,32 @@ def transformation(self) -> Transformation: return self._transformation +class MultivariateParameter(Parameter): + def __init__( + self, + distribution_param, + distribution: stats.distributions.rv_frozen | Distribution | None = None, + bounds: BoundsPair | None = None, + initial_value: float = None, + transformation: Transformation | None = None, + ): + if isinstance(distribution, BaseMultivariateDistribution): + super().__init__( + bounds=bounds, + initial_value=initial_value, + transformation=transformation, + ) + self._distribution = distribution + else: + super().__init__( + distribution=distribution, + bounds=bounds, + initial_value=initial_value, + transformation=transformation, + ) + self.distribution_param = distribution_param + + class Parameters: """ Container for managing multiple Parameter objects with additional functionality. @@ -253,12 +280,12 @@ def __init__(self, parameters: dict | Parameters = None) -> None: raise TypeError( "parameters must be either a dictionary or a pybop.Parameters instance" ) - self._parameters = OrderedDict() for name, param in parameters.items(): - self._add(name, param, update_transform=False) + self._add(name, param, update_transform=False, check_multivariate=False) self._transform = self.construct_transformation() + self.check_multivariate() def __getitem__(self, name: str) -> Parameter: return self.get(name) @@ -280,12 +307,37 @@ def names(self) -> list[str]: def __iter__(self) -> Iterator[Parameter]: return iter(self._parameters.values()) - def add(self, name: str, parameter: Parameter) -> None: + def add(self, name: str, parameter: Parameter, check_multivariate=True) -> None: """Add a parameter to the collection.""" - self._add(name, parameter) + self._add(name, parameter, check_multivariate=check_multivariate) + + def check_multivariate(self): + self._multivariate = any( + isinstance(param, MultivariateParameter) for param in self + ) + if self._multivariate: + if not all(isinstance(param, MultivariateParameter) for param in self): + raise TypeError( + "A MultivariateParameter cannot be combined with other types of parameters" + ) + dist_param = next(iter(self._parameters.values())).distribution_param + if not all(param.distribution_param == dist_param for param in self): + raise ValueError( + "All MultivariateParameters must share the same distribution." + ) + self.distribution = self._parameters[dist_param].distribution + + if not isinstance(self.distribution, BaseMultivariateDistribution): + raise TypeError( + "The distribution provided for MultivariateParameters must be a BaseMultivariateDistribution." + ) def _add( - self, name: str, parameter: Parameter, update_transform: bool = True + self, + name: str, + parameter: Parameter, + update_transform: bool = True, + check_multivariate=True, ) -> None: """ Internal method to add a parameter to the collection. @@ -307,6 +359,8 @@ def _add( if update_transform: self._transform = self.construct_transformation() + if check_multivariate: + self.check_multivariate() def remove(self, name: str) -> Parameter: """Remove parameter and return it.""" @@ -326,23 +380,27 @@ def join(self, parameters=None): """ for name, param in parameters.items(): if name not in self._parameters.keys(): - self.add(name, param) + self.add(name, param, check_multivariate=False) else: print(f"Discarding duplicate {name}.") + self.check_multivariate() + def get(self, name: str) -> Parameter: """Get a parameter by name.""" if name not in self._parameters: raise ParameterNotFoundError(f"Parameter for '{name}' not found") return self._parameters[name] - def set(self, name: str, param: Parameter) -> None: + def set(self, name: str, param: Parameter, check_multivariate=True) -> None: """Get a parameter by name.""" if name not in self._parameters: raise ParameterNotFoundError(f"Parameter for '{name}' not found") if not isinstance(param, Parameter): - raise TypeError({"Paremeter must be of type pybop.ParemterInfo"}) + raise TypeError({"Paremeter must be of type pybop.Parameter"}) self._parameters[name] = param + if check_multivariate: + self.check_multivariate() def get_bounds(self, transformed: bool = False) -> dict: """ @@ -449,22 +507,47 @@ def sample_from_distribution( """ Sample from each parameter distribution. + or + + Draw random samples from the joint parameters distribution for multivariate parameters. + + Parameters + ---------- + n_samples : int + The number of samples to draw (default: 1). + random_state : int, optional + The random state seed for reproducibility (default: None). + transformed: bool + If True, the transformation is applied to the output + (default: False). + Returns ------- NDArray[np.floating] or None Array of shape (n_samples, n_parameters) or None if any distribution is missing """ - all_samples = [] - for param in self._parameters.values(): - samples = param.sample_from_distribution( - n_samples, random_state=random_state, transformed=transformed - ) - if samples is None: - return None - all_samples.append(samples) + if self._multivariate: + samples = self.distribution.rvs(n_samples, random_state=random_state) + if samples.ndim < 2: + samples = np.atleast_2d(samples) + + if transformed: + samples = np.asarray([self.transformation.to_model(s) for s in samples]) + + return samples + else: + all_samples = [] + + for param in self._parameters.values(): + samples = param.sample_from_distribution( + n_samples, random_state=random_state, transformed=transformed + ) + if samples is None: + return None + all_samples.append(samples) - return np.column_stack(all_samples) + return np.column_stack(all_samples) def get_sigma0(self, transformed: bool = False) -> list: """ diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index dcfa35ee8..c12932a96 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -138,18 +138,22 @@ def multivariate_simulator(self, model, dataset): model, parameter_values=parameter_values, protocol=dataset ) # Override the forced univariate Parameters - simulator.parameters = pybop.MultivariateParameters( + simulator.parameters = pybop.Parameters( { - "Negative electrode active material volume fraction": pybop.Parameter( - initial_value=0.6, bounds=[0.001, 0.999] + "Negative electrode active material volume fraction": pybop.MultivariateParameter( + distribution_param="Negative electrode active material volume fraction", + distribution=pybop.MultivariateGaussian( + [0.6, 0.5], [[0.02, 0.0], [0.0, 0.05]] + ), + initial_value=0.6, + bounds=[0.001, 0.999], ), - "Positive electrode active material volume fraction": pybop.Parameter( - initial_value=0.5, bounds=[0.001, 0.999] + "Positive electrode active material volume fraction": pybop.MultivariateParameter( + distribution_param="Negative electrode active material volume fraction", + initial_value=0.5, + bounds=[0.001, 0.999], ), }, - distribution=pybop.MultivariateGaussian( - [0.6, 0.5], [[0.02, 0.0], [0.0, 0.05]] - ), ) return simulator @@ -157,8 +161,6 @@ def multivariate_simulator(self, model, dataset): def multivariate_problem(self, multivariate_simulator, dataset): cost = pybop.SumSquaredError(dataset) problem = pybop.Problem(multivariate_simulator, cost) - # Copy the MultivariateParameters to the problem - problem.parameters = multivariate_simulator.parameters return problem @pytest.fixture @@ -181,8 +183,6 @@ def gitt_like_problem(self, multivariate_simulator, dataset): pybop.Problem(multivariate_simulator, sqrt_cost_1), pybop.Problem(multivariate_simulator, sqrt_cost_2), ) - # Copy the MultivariateParameters to the problem - problem.parameters = multivariate_simulator.parameters return problem @pytest.mark.parametrize( diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py index 0c8c56016..025eb5456 100644 --- a/tests/unit/test_parameters.py +++ b/tests/unit/test_parameters.py @@ -146,7 +146,7 @@ def test_parameters_construction(self, name, parameter): with pytest.raises(ParameterNotFoundError, match="not found"): params["not a parameter"] = pybop.Parameter(initial_value=0.8) with pytest.raises( - TypeError, match="Paremeter must be of type pybop.ParemterInfo" + TypeError, match="Paremeter must be of type pybop.Parameter" ): params[name] = pybop.Gaussian(0.5, 0.02) @@ -274,7 +274,7 @@ def test_parameters_repr(self, name, parameter): ) -class TestMultivariateParameters: +class TestMultivariateParameter: """ A class to test the multivariate parameters class. """ @@ -283,31 +283,134 @@ class TestMultivariateParameters: @pytest.fixture def multivariate_parameters(self): - return pybop.MultivariateParameters( + return pybop.Parameters( { - "Negative particle diffusivity [m2.s-1]": pybop.Parameter( + "Negative particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + distribution=pybop.MultivariateGaussian( + [np.log(3.9e-14), np.log(1e-15)], + [[np.log(10), 0.0], [0.0, np.log(10)]], + ), initial_value=3.9e-14, bounds=[3.9e-15, 3.9e-13], transformation=pybop.LogTransformation(), ), - "Positive particle diffusivity [m2.s-1]": pybop.Parameter( + "Positive particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", initial_value=1e-15, bounds=[1e-16, 1e-14], transformation=pybop.LogTransformation(), ), }, - distribution=pybop.MultivariateGaussian( - [np.log(3.9e-14), np.log(1e-15)], - [[np.log(10), 0.0], [0.0, np.log(10)]], - ), ) def test_rvs(self, multivariate_parameters): - samples = multivariate_parameters.sample_from_distribution( - 1, apply_transform=True - ) + samples = multivariate_parameters.sample_from_distribution(1, transformed=True) assert samples.shape == (1, 2) assert samples.T[1].min() >= 1e-16 assert samples.T[1].max() <= 1e-14 - assert multivariate_parameters.pdf(np.asarray([3.9e-14, 1e-15])) > 0 + # assert multivariate_parameters.pdf(np.asarray([3.9e-14, 1e-15])) > 0 assert multivariate_parameters.distribution is not None + + def test_input_checks_multivariate_parameters(self): + with pytest.raises( + TypeError, + match="A MultivariateParameter cannot be combined with other types of parameters", + ): + pybop.Parameters( + { + "Negative particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + distribution=pybop.MultivariateGaussian( + [np.log(3.9e-14), np.log(1e-15)], + [[np.log(10), 0.0], [0.0, np.log(10)]], + ), + initial_value=3.9e-14, + bounds=[3.9e-15, 3.9e-13], + transformation=pybop.LogTransformation(), + ), + "Positive particle diffusivity [m2.s-1]": pybop.Parameter( + initial_value=1e-15, + bounds=[1e-16, 1e-14], + transformation=pybop.LogTransformation(), + ), + }, + ) + + params = pybop.Parameters( + { + "Negative particle diffusivity [m2.s-1]": pybop.Parameter( + initial_value=3.9e-14, + bounds=[3.9e-15, 3.9e-13], + transformation=pybop.LogTransformation(), + ), + "Positive particle diffusivity [m2.s-1]": pybop.Parameter( + initial_value=1e-15, + bounds=[1e-16, 1e-14], + transformation=pybop.LogTransformation(), + ), + }, + ) + + with pytest.raises( + TypeError, + match="A MultivariateParameter cannot be combined with other types of parameters", + ): + params["Negative particle diffusivity [m2.s-1]"] = ( + pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + distribution=pybop.MultivariateGaussian( + [np.log(3.9e-14), np.log(1e-15)], + [[np.log(10), 0.0], [0.0, np.log(10)]], + ), + initial_value=3.9e-14, + bounds=[3.9e-15, 3.9e-13], + transformation=pybop.LogTransformation(), + ) + ) + + with pytest.raises( + ValueError, + match="All MultivariateParameters must share the same distribution.", + ): + pybop.Parameters( + { + "Negative particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + distribution=pybop.MultivariateGaussian( + [np.log(3.9e-14), np.log(1e-15)], + [[np.log(10), 0.0], [0.0, np.log(10)]], + ), + initial_value=3.9e-14, + bounds=[3.9e-15, 3.9e-13], + transformation=pybop.LogTransformation(), + ), + "Positive particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Positive particle diffusivity [m2.s-1]", + initial_value=1e-15, + bounds=[1e-16, 1e-14], + transformation=pybop.LogTransformation(), + ), + }, + ) + + with pytest.raises( + TypeError, + match="The distribution provided for MultivariateParameters must be a BaseMultivariateDistribution.", + ): + pybop.Parameters( + { + "Negative particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + distribution=pybop.Gaussian(0.0, 1.0), + initial_value=3.9e-14, + transformation=pybop.LogTransformation(), + ), + "Positive particle diffusivity [m2.s-1]": pybop.MultivariateParameter( + distribution_param="Negative particle diffusivity [m2.s-1]", + initial_value=1e-15, + bounds=[1e-16, 1e-14], + transformation=pybop.LogTransformation(), + ), + }, + )