Multivariate parameters

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(

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
 

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,

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