diff --git a/CHANGELOG.md b/CHANGELOG.md index b81be05cb..999f919bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ ## Optimisations +- [#925](https://github.com/pybop-team/PyBOP/pull/925) - Add `UnboundedDistribution` and the `get_transformed_distribution` functionality. + ## Bug Fixes - [#915](https://github.com/pybop-team/PyBOP/pull/915) - Fixes axis labels for non-standard domain names, adds `Dataset` length property and adds `kind` property to `Interpolant`. diff --git a/examples/notebooks/getting_started/using_transformations.ipynb b/examples/notebooks/getting_started/using_transformations.ipynb index 236dc8bba..b252dc04d 100644 --- a/examples/notebooks/getting_started/using_transformations.ipynb +++ b/examples/notebooks/getting_started/using_transformations.ipynb @@ -318,7 +318,9 @@ "problem.parameters[\"Positive electrode active material volume fraction\"] = (\n", " pybop.Parameter(\n", " initial_value=old_param.initial_value,\n", - " bounds=old_param.bounds,\n", + " distribution=pybop.LogUniform(\n", + " lower=old_param.bounds[0], upper=old_param.bounds[1]\n", + " ),\n", " transformation=pybop.LogTransformation(),\n", " )\n", ")\n", diff --git a/examples/scripts/getting_started/maximum_a_posteriori.py b/examples/scripts/getting_started/maximum_a_posteriori.py index 5ad09eaf5..29e423d64 100644 --- a/examples/scripts/getting_started/maximum_a_posteriori.py +++ b/examples/scripts/getting_started/maximum_a_posteriori.py @@ -41,12 +41,10 @@ "Negative electrode active material volume fraction": pybop.Parameter( distribution=pybop.Uniform(0.3, 0.8), initial_value=0.653, - transformation=pybop.LogTransformation(), ), "Positive electrode active material volume fraction": pybop.Parameter( distribution=pybop.Uniform(0.3, 0.8), initial_value=0.657, - transformation=pybop.LogTransformation(), ), } ) diff --git a/examples/scripts/getting_started/monte_carlo_sampling.py b/examples/scripts/getting_started/monte_carlo_sampling.py index 7e557843c..f588bb6af 100644 --- a/examples/scripts/getting_started/monte_carlo_sampling.py +++ b/examples/scripts/getting_started/monte_carlo_sampling.py @@ -1,3 +1,4 @@ +import numpy as np import pybamm import pybop @@ -32,11 +33,11 @@ parameter_values.update( { "Negative electrode active material volume fraction": pybop.Parameter( - distribution=pybop.Gaussian(0.68, 0.02), + distribution=pybop.LogNormal(np.log(0.68), 0.02), transformation=pybop.LogTransformation(), ), "Positive electrode active material volume fraction": pybop.Parameter( - distribution=pybop.Gaussian(0.65, 0.02), + distribution=pybop.LogNormal(np.log(0.65), 0.02), transformation=pybop.LogTransformation(), ), } diff --git a/pybop/__init__.py b/pybop/__init__.py index 7094fdfe2..11436f11c 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -49,7 +49,7 @@ # Parameter classes # from .parameters.parameter import Parameter, Parameters -from .parameters.distributions import Distribution, Gaussian, Uniform, Exponential, JointDistribution +from .parameters.distributions import Distribution, Exponential, Gaussian, JointDistribution, LogNormal, LogUniform, Unbounded, Uniform from .parameters.multivariate_distributions import MultivariateNonparametric, MultivariateUniform, MultivariateGaussian, MultivariateLogNormal, MarginalDistribution # diff --git a/pybop/optimisers/base_optimiser.py b/pybop/optimisers/base_optimiser.py index cad70cd58..a8914bc21 100644 --- a/pybop/optimisers/base_optimiser.py +++ b/pybop/optimisers/base_optimiser.py @@ -144,8 +144,6 @@ def run(self) -> "OptimisationResult": results = [] for i in range(self._multistart): if i >= 1: - if self.problem.parameters.distribution is None: - raise RuntimeError("Distributions must be provided for multi-start") initial_values = self.problem.parameters.sample_from_distribution(1)[0] self.problem.parameters.update(initial_values=initial_values) self._set_up_optimiser() diff --git a/pybop/parameters/distributions.py b/pybop/parameters/distributions.py index 4409be2f8..326224a62 100644 --- a/pybop/parameters/distributions.py +++ b/pybop/parameters/distributions.py @@ -1,6 +1,13 @@ import numpy as np import scipy.stats as stats +from pybop.transformation.base_transformation import Transformation +from pybop.transformation.transformations import ( + IdentityTransformation, + LogTransformation, + ScaledTransformation, +) + class Distribution: """ @@ -15,16 +22,25 @@ class Distribution: ---------- distribution : scipy.stats.distributions.rv_frozen The underlying continuous random variable distribution. + n_parameters : int + The number of dimensions (default: 1). """ def __init__( self, distribution: stats.distributions.rv_frozen | None = None, + properties: dict | None = None, + n_parameters: int = 1, ): self.distribution = distribution + self.properties = properties or {} + self._n_parameters = n_parameters + + def support(self) -> tuple[float]: + if self.distribution is None: + return (-np.inf, np.inf) - def support(self): - return self.distribution.support() + return tuple(float(x) for x in self.distribution.support()) def pdf(self, x): """ @@ -42,8 +58,8 @@ def pdf(self, x): """ if self.distribution is None: raise NotImplementedError - else: - return self.distribution.pdf(x) + + return self.distribution.pdf(x) def logpdf(self, x): """ @@ -61,8 +77,8 @@ def logpdf(self, x): """ if self.distribution is None: raise NotImplementedError - else: - return self.distribution.logpdf(x) + + return self.distribution.logpdf(x) def icdf(self, q): """ @@ -80,8 +96,8 @@ def icdf(self, q): """ if self.distribution is None: raise NotImplementedError - else: - return self.distribution.ppf(q) + + return self.distribution.ppf(q) def cdf(self, x): """ @@ -99,10 +115,10 @@ def cdf(self, x): """ if self.distribution is None: raise NotImplementedError - else: - return self.distribution.cdf(x) - def rvs(self, size=1, random_state=None): + return self.distribution.cdf(x) + + def rvs(self, size: int = 1, random_state: int | None = None): """ Generates random variates from the distribution. @@ -123,19 +139,10 @@ def rvs(self, size=1, random_state=None): ValueError If the size parameter is negative. """ - if not isinstance(size, int | tuple): - raise ValueError( - "size must be a positive integer or tuple of positive integers" - ) - if isinstance(size, int) and size < 1: - raise ValueError("size must be a positive integer") - if isinstance(size, tuple) and any(s < 1 for s in size): - raise ValueError("size must be a tuple of positive integers") - if self.distribution is None: raise NotImplementedError - else: - return self.distribution.rvs(size=size, random_state=random_state) + + return self.distribution.rvs(size=size, random_state=random_state) def logpdfS1(self, x): """ @@ -171,115 +178,149 @@ def _dlogpdf_dx(self, x): float The value(s) of the first derivative at x. """ - if self.distribution is None: - raise NotImplementedError - else: - # Use a finite difference approximation of the gradient - delta = max(abs(x) * 1e-3, np.finfo(float).eps) - log_distribution_upper = self.logpdf(x + delta) - log_distribution_lower = self.logpdf(x - delta) + # Use a finite difference approximation of the gradient + delta = max(abs(x) * 1e-3, np.finfo(float).eps) + log_distribution_upper = self.logpdf(x + delta) + log_distribution_lower = self.logpdf(x - delta) - return (log_distribution_upper - log_distribution_lower) / (2 * delta) + return (log_distribution_upper - log_distribution_lower) / (2 * delta) - def verify(self, x): - """ - Verifies that the input is a numpy array and converts it if necessary. - """ + def verify(self, x) -> np.ndarray: + """Verifies that the input is a numpy array and converts it if necessary.""" if isinstance(x, dict): x = np.asarray(list(x.values())) elif not isinstance(x, np.ndarray): x = np.asarray(x) return x - def __repr__(self): - """ - Returns a string representation of the object. - """ - return f"{self.__class__.__name__}, mean: {self.mean()}, standard deviation: {self.std()}" + def __repr__(self) -> str: + return f"{self.name}, properties: {self.properties}" - def mean(self): - """ - Get the mean of the distribution. + def mean(self) -> float: + """Get the mean of the distribution.""" + if self.distribution is None: + raise NotImplementedError - Returns - ------- - float - The mean of the distribution. - """ return self.distribution.mean() - def std(self): - """ - Get the standard deviation of the distribution. + def std(self) -> float: + """Get the standard deviation of the distribution.""" + if self.distribution is None: + raise NotImplementedError - Returns - ------- - float - The standard deviation of the distribution. - """ return self.distribution.std() + def get_transformed_distribution(self, transform: Transformation): + """Get the transformed distribution in the search space.""" + if isinstance(transform, IdentityTransformation): + return self + + transformed_distribution = self._transform(transform) + + if transformed_distribution is None: + raise TypeError( + f"The transformation of a {self.name} distribution by a " + f"{transform.name} is undefined or not yet implemented." + ) + return transformed_distribution + + def _transform(self, transform: Transformation): + """Get the transformed distribution in the search space.""" + return NotImplementedError + + @property + def name(self): + return self.__class__.__name__ + class Gaussian(Distribution): """ Represents a Gaussian (normal) distribution with a given mean and standard deviation. - This class provides methods to calculate the probability density function (pdf), - the logarithm of the pdf, and to generate random variates (rvs) from the distribution. - Parameters ---------- mean : float The mean (mu) of the Gaussian distribution. sigma : float The standard deviation (sigma) of the Gaussian distribution. + truncated_at : list[float], optional + If provided, the distribution becomes a truncated normal distribution. """ def __init__( - self, - mean, - sigma, - truncated_at: list[float] = None, + self, mean: float, sigma: float, truncated_at: list[float] | None = None ): if truncated_at is not None: - distribution = stats.truncnorm( - (truncated_at[0] - mean) / sigma, - (truncated_at[1] - mean) / sigma, - loc=mean, - scale=sigma, - ) + properties = { + "a": float((truncated_at[0] - mean) / sigma), + "b": float((truncated_at[1] - mean) / sigma), + "loc": float(mean), + "scale": float(sigma), + } + distribution = stats.truncnorm(**properties) else: - distribution = stats.norm(loc=mean, scale=sigma) - super().__init__(distribution) - self.name = "Gaussian" - self._n_parameters = 1 - self.loc = mean - self.scale = sigma + properties = {"loc": float(mean), "scale": float(sigma)} + distribution = stats.norm(**properties) + super().__init__(distribution, properties=properties) def _dlogpdf_dx(self, x): - """ - Evaluates the first derivative of the log Gaussian distribution at x. + """Evaluates the first derivative of the log Gaussian distribution at x.""" + if "a" in self.properties.keys(): + return super()._dlogpdf_dx(x) # use estimate for a truncated normal + + return (self.properties["loc"] - x) / self.properties["scale"] ** 2 + + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if isinstance(transform, ScaledTransformation): + truncated_at = None + if "a" in self.properties.keys(): + original_truncation = self.properties["loc"] + ( + self.properties["scale"] + * np.asarray([self.properties["a"], self.properties["b"]]) + ) + truncated_at = np.sort( + [transform.to_search(x).item() for x in original_truncation] + ) + return Gaussian( + mean=transform.to_search(self.properties["loc"]).item(), + sigma=self.properties["scale"] * transform.coefficient.item(), + truncated_at=truncated_at, + ) + return None - Parameters - ---------- - x : float - The point(s) at which to evaluate the first derivative. - Returns - ------- - float - The value(s) of the first derivative at x. - """ - return (self.loc - x) / self.scale**2 +class LogNormal(Distribution): + """ + Represents a log-normal distribution corresponding to a Gaussian distribution for + log(x) with a given mean and standard deviation. + + Parameters + ---------- + mean_log_x : float + The mean (mu) of the Gaussian distribution of log(x). + sigma : float + The standard deviation (sigma) of the Gaussian distribution of log(x). + """ + + def __init__(self, mean_log_x: float, sigma: float): + properties = {"scale": float(np.exp(mean_log_x)), "s": float(sigma)} + super().__init__(stats.lognorm(**properties), properties=properties) + + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if isinstance(transform, LogTransformation): + return Gaussian( + mean=np.log(self.properties["scale"]), + sigma=self.properties["s"], + ) + return None class Uniform(Distribution): """ Represents a uniform distribution over a specified interval. - This class provides methods to calculate the pdf, the log pdf, and to generate - random variates from the distribution. - Parameters ---------- lower : float @@ -288,91 +329,136 @@ class Uniform(Distribution): The upper bound of the distribution. """ - def __init__( - self, - lower, - upper, - ): - super().__init__(stats.uniform(loc=lower, scale=upper - lower)) - self.name = "Uniform" - self.lower = lower - self.upper = upper - self._n_parameters = 1 + def __init__(self, lower: float, upper: float): + properties = {"loc": float(lower), "scale": float(upper - lower)} + super().__init__(stats.uniform(**properties), properties=properties) def _dlogpdf_dx(self, x): - """ - Evaluates the first derivative of the log uniform distribution at x. + """Evaluates the first derivative of the log uniform distribution at x.""" + return np.zeros_like(x) - Parameters - ---------- - x : float - The point(s) at which to evaluate the first derivative. + def mean(self) -> float: + """Returns the mean of the distribution.""" + return self.properties["loc"] + self.properties["scale"] / 2 - Returns - ------- - float - The value(s) of the first derivative at x. - """ - return np.zeros_like(x) + def __repr__(self) -> str: + return f"{self.name}, bounds: {self.support()}" - def mean(self): - """ - Returns the mean of the distribution. - """ - return (self.lower + self.upper) / 2 + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if isinstance(transform, ScaledTransformation): + bounds = [transform.to_search(x) for x in self.support()] + return Uniform(lower=np.min(bounds), upper=np.max(bounds)) + return None - def __repr__(self): - """ - Returns a string representation of the object. - """ - return f"{self.__class__.__name__}, lower: {self.lower}, upper: {self.upper}" +class LogUniform(Distribution): + """ + Represents a log-uniform distribution over a specified interval. -class Exponential(Distribution): + Parameters + ---------- + lower : float + The lower bound of the distribution. + upper : float + The upper bound of the distribution. """ - Represents an exponential distribution with a specified scale parameter. - This class provides methods to calculate the pdf, the log pdf, and to generate random - variates from the distribution. + def __init__(self, lower: float, upper: float): + # Validate that upper > lower for all elements + if not np.all(0 < lower < upper): + raise ValueError( + "All elements of upper bounds must be greater than lower bounds." + ) + properties = {"a": float(lower), "b": float(upper)} + super().__init__(stats.loguniform(**properties), properties=properties) + + def __repr__(self) -> str: + return f"{self.name}, bounds: {self.support()}" + + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if isinstance(transform, ScaledTransformation): + bounds = [transform.to_search(x) for x in self.support()] + return LogUniform(lower=np.min(bounds), upper=np.max(bounds)) + + elif isinstance(transform, LogTransformation): + bounds = [transform.to_search(x) for x in self.support()] + return Uniform(lower=np.min(bounds), upper=np.max(bounds)) + return None + + +class Unbounded(Distribution): + """ + Represents an unbounded distribution with either zero or one finite bound. Parameters ---------- - scale : float - The scale parameter (lambda) of the exponential distribution. + lower : float + The lower bound of the distribution. + upper : float + The upper bound of the distribution. """ def __init__( - self, - scale: float, - loc: float = 0, + self, initial_value: float = 1.0, lower: float = -np.inf, upper: float = np.inf ): - super().__init__(stats.expon(loc=loc, scale=scale)) - self.name = "Exponential" - self._n_parameters = 1 - self.loc = loc - self.scale = scale + super().__init__() + self.lower = float(lower) + self.upper = float(upper) + self.initial_value = ( + None + if initial_value is None + else float(np.minimum(np.maximum(initial_value, lower), upper)) + ) + + def support(self) -> tuple[float]: + return (self.lower, self.upper) def _dlogpdf_dx(self, x): - """ - Evaluates the first derivative of the log exponential distribution at x. + """Evaluates the first derivative of the log uniform distribution at x.""" + return np.zeros_like(x) - Parameters - ---------- - x : float - The point(s) at which to evaluate the first derivative. + def mean(self) -> float: + """Returns an estimate for the mean of the distribution.""" + return self.initial_value - Returns - ------- - float - The value(s) of the first derivative at x. - """ - return -1 / self.scale * np.ones_like(x) + def std(self) -> float: + """Returns an estimate for the standard deivation of the distribution.""" + return 0.05 if self.initial_value == 0 else 0.05 * self.initial_value - def __repr__(self): - """ - Returns a string representation of the object. - """ - return f"{self.__class__.__name__}, loc: {self.loc}, scale: {self.scale}" + def __repr__(self) -> str: + return f"{self.name}, initial value: {self.initial_value}, bounds: {self.support()}" + + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if isinstance(transform, ScaledTransformation | LogTransformation): + bounds = [transform.to_search(x) for x in self.support()] + return Unbounded( + initial_value=transform.to_search(self.initial_value), + lower=np.min(bounds), + upper=np.max(bounds), + ) + return None + + +class Exponential(Distribution): + """ + Represents an exponential distribution with a specified scale parameter. + + Parameters + ---------- + scale : float + The scale parameter (lambda) of the exponential distribution. + """ + + def __init__(self, scale: float, loc: float = 0): + properties = {"loc": float(loc), "scale": float(scale)} + super().__init__(stats.expon(**properties), properties=properties) + + def _dlogpdf_dx(self, x): + """Evaluates the first derivative of the log exponential distribution at x.""" + return -1 / self.properties["scale"] * np.ones_like(x) class JointDistribution(Distribution): @@ -391,23 +477,30 @@ def __init__(self, *distributions: Distribution | stats.distributions.rv_frozen) if all(distribution is None for distribution in distributions): return - if not all( - isinstance(distribution, (Distribution, stats.distributions.rv_frozen)) - for distribution in distributions - ): - raise ValueError( - "All distributions must be instances of Distribution or scipy.stats.distributions.rv_frozen" - ) + for distribution in distributions: + if not isinstance(distribution, Distribution): + raise ValueError( + "All distributions must be instances of Distribution. " + f"Received {distribution}" + ) self._n_parameters = len(distributions) - self._distributions_list: list[Distribution] = [ - distribution - if isinstance(distribution, Distribution) - else Distribution(distribution) - for distribution in distributions - ] - - def rvs(self, size=1, random_state=None): + self._distributions_list: list[Distribution] = list(distributions) + + def support(self) -> np.ndarray: + """Return the support as a numpy array of dimensions (2, n_parameters).""" + return np.asarray([d.support() for d in self._distributions_list]).T + + def mean(self): + return np.asarray([d.mean() for d in self._distributions_list]).T + + def std(self): + return np.asarray([d.std() for d in self._distributions_list]) + + def cov(self): + return np.diag(self.std()) ** 2 + + def rvs(self, size: int = 1, random_state: int | None = None): """Sample each distribution individually and then compile.""" all_samples = [] for distribution in self._distributions_list: @@ -436,8 +529,7 @@ def logpdf(self, x: float | np.ndarray) -> float: ) return sum( - distribution.logpdf(x) - for distribution, x in zip(self._distributions_list, x, strict=False) + d.logpdf(x) for d, x in zip(self._distributions_list, x, strict=False) ) def logpdfS1(self, x: float | np.ndarray) -> tuple[float, np.ndarray]: @@ -477,7 +569,22 @@ def logpdfS1(self, x: float | np.ndarray) -> tuple[float, np.ndarray]: return output, doutput.T def __repr__(self) -> str: - distributions_repr = "; ".join( - [repr(distribution) for distribution in self._distributions_list] - ) - return f"{self.__class__.__name__}(distributions: [{distributions_repr}])" + distributions_repr = "; ".join([repr(d) for d in self._distributions_list]) + return f"{self.name}(distributions: [{distributions_repr}])" + + def marginal(self, position: int): + return self._distributions_list[position] + + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + list_of_transforms = [] + for t, d in zip( + transform.transformations, self._distributions_list, strict=False + ): + transformed_dist = d.get_transformed_distribution(t) + if transformed_dist is None: + return None + else: + list_of_transforms.append(transformed_dist) + + return JointDistribution(*list_of_transforms) diff --git a/pybop/parameters/multivariate_distributions.py b/pybop/parameters/multivariate_distributions.py index c13fde329..ce69684ee 100644 --- a/pybop/parameters/multivariate_distributions.py +++ b/pybop/parameters/multivariate_distributions.py @@ -3,14 +3,13 @@ import scipy.stats as stats from scipy.linalg import sqrtm -from pybop.parameters.distributions import Distribution +from pybop.parameters.distributions import Distribution, Gaussian, LogNormal, Uniform from pybop.transformation.base_transformation import Transformation from pybop.transformation.transformations import ( ComposedTransformation, IdentityTransformation, LogTransformation, ScaledTransformation, - UnitHyperCube, ) @@ -34,67 +33,103 @@ class BaseMultivariateDistribution(Distribution): keys and their values as float values. """ + def support(self) -> np.ndarray: + """Return the support as a numpy array of dimensions (2, n_parameters).""" + return self.distribution.support(**self.properties) + def pdf(self, x): + """Calculates the probability density function (PDF) of the distribution at x.""" return self.distribution.pdf(x, **self.properties) def logpdf(self, x): + """Calculates the logarithm of the probability density function of the distribution at x.""" return self.distribution.logpdf(x, **self.properties) icdf = None """Multivariate distributions have no invertible CDF.""" - def icdf_marginal(self, q, dim): + def icdf_marginal(self, q, position: int): raise NotImplementedError def cdf(self, x): return self.distribution.cdf(x, **self.properties) - def cdf_marginal(self, x, dim): + def cdf_marginal(self, x, position: int): """ - Takes the marginal cumulative distribution function (CDF) at x - in dim. + Takes the marginal cumulative distribution function (CDF) at x. Parameters ---------- x : numpy.ndarray The point(s) at which to evaluate the CDF. - dim : int + position : int The dimension to which to reduce the CDF to. Returns ------- float - The marginal cumulative distribution function value at x in - dim. + The marginal cumulative distribution function value at x. """ - return integrate.quad(self.pdf_marginal, -np.inf, x)[0] + return integrate.quad(self.marginal(position), -np.inf, x)[0] - def rvs(self, size=1, random_state=None): + def rvs(self, size: int = 1, random_state: int | None = None): + """Generates random variates from the distribution.""" return self.distribution.rvs( **self.properties, size=size, random_state=random_state ) - def __repr__(self): - return f"{self.__class__.__name__}, properties: {self.properties}" + def cov(self): + """Get the covariance matrix.""" + return self.properties["cov"] - @property - def mean(self): - raise NotImplementedError + def std(self): + """Get the square root of the covariance matrix.""" + return sqrtm(self.cov()) + + def get_transformed_distribution( + self, transform: Transformation + ) -> Distribution | None: + """Get the transformed distribution in the search space.""" + transform = ( + transform + if isinstance(transform, ComposedTransformation) + else ComposedTransformation([transform] * self._n_parameters) + ) + self.check_consistent_transformation(transform) - @property - def sigma(self): - raise NotImplementedError + if isinstance(transform, IdentityTransformation) or all( + isinstance(t, IdentityTransformation) for t in transform.transformations + ): + return self - @property - def compatible_transformations(self): - raise NotImplementedError + transformed_distribution = self._transform(transform) - def marginal(self, position): - """Return univariate marginal distribution of parameter with index 'position'""" - raise NotImplementedError + if transformed_distribution is None: + transform_name = ( + transform.transformations[0].name + if isinstance(transform, ComposedTransformation) + else transform.name + ) + raise TypeError( + f"The transformation of a {self.name} distribution by a " + f"{transform_name} is undefined or not yet implemented." + ) + return transformed_distribution + + def check_consistent_transformation(self, transform: Transformation) -> None: + """Raise an error if a composed transformation contains incompatible types of transformation.""" + for transformation_category in [ScaledTransformation, LogTransformation]: + transformation_in_category: list[bool] = [ + isinstance(t, transformation_category) + for t in transform.transformations + ] + if any(transformation_in_category) and not all(transformation_in_category): + raise TypeError( + f"All transformations must be of a similar type. Received {[t.name for t in transform.transformations]}" + ) - def transformed_properties(self, transformation): - """Return distribution properties in search space""" + def marginal(self, position: int): + """Return univariate marginal distribution of parameter with index 'position'.""" raise NotImplementedError @@ -103,84 +138,69 @@ class MultivariateNonparametric(BaseMultivariateDistribution): Represents a "freeform" distribution, i.e., one that is defined from a random sampling and a kernel density estimate on that sampling. - This class provides methods to calculate the probability density - function (pdf), the logarithm of the pdf, and to generate random - variates (rvs) from the distribution. - Parameters ---------- samples : numpy.ndarray The random variates to base the distribution on. """ - def __init__(self, samples, random_state=None): - super().__init__(distribution=stats.gaussian_kde(samples)) - self.name = "MultivariateNonparametric" - self.properties = {} - self._n_parameters = samples.shape[0] - - @property - def compatible_transformations(self): - return Transformation + def __init__(self, samples: np.ndarray): + super().__init__( + distribution=stats.gaussian_kde(samples), n_parameters=samples.shape[0] + ) - def pdf(self, x): - return self.distribution.pdf(x) + def mean(self): + # scipy.stats.gaussian_kde does not have the method "mean" + # compute mean based on dataset + return np.mean(self.distribution.dataset) - def logpdf(self, x): - return self.distribution.logpdf(x) + def std(self): + # scipy.stats.gaussian_kde does not have the method "std" + # compute standard deviation based on dataset + return np.std(self.distribution.dataset, ddof=1) - def rvs(self, size=1, random_state=None): + def rvs(self, size: int = 1, random_state: int | None = None): return self.distribution.resample(size, random_state) - def marginal(self, position): - return self.distribution.marginal(position) - - def transformed_properties(self, transformation): - """Return distribution properties in search space""" - - # properties is an empty dictionary and hence independent of transformation - return self.properties + def marginal(self, position: int): + return MultivariateNonparametric(self.distribution.marginal(position).dataset) class MultivariateUniform(BaseMultivariateDistribution): """ Represents a multivariate uniform distribution. - This class provides methods to calculate the probability density - fuction (pdf), the logarithm of the pdf, and to generate random - variates (rvs) from the distribution. - Parameters ---------- bounds : numpy.ndarray - The lower and upper bounds for the uniform distribution. + The lower and upper bounds for the uniform distribution as an array + with dimensions (2, n_parameters). """ - def __init__(self, bounds, random_state=None): + def __init__(self, bounds: np.ndarray): + bounds = np.asarray(bounds) super().__init__( - distribution=stats.uniform( - loc=bounds[0, :], scale=bounds[1, :] - bounds[0, :] - ) + distribution=stats.uniform, + properties={"loc": bounds[0, :], "scale": bounds[1, :] - bounds[0, :]}, + n_parameters=bounds.shape[1], ) - self.name = "MultivariateUniform" - self.properties = {"bounds": bounds} - self._n_parameters = bounds.shape[1] - @property - def compatible_transformations(self): - return ( - IdentityTransformation, - ScaledTransformation, - UnitHyperCube, + def marginal(self, position: int): + # return univariate uniform distribution + return Uniform( + lower=self.properties["loc"][position], + upper=self.properties["loc"][position] + self.properties["scale"][position], ) - def marginal(self, position): - # return univariate unfiform distribution - return stats.uniform( - loc=self.properties["bounds"][0, position], - scale=self.properties["bounds"][1, position] - - self.properties["bounds"][0, position], - ) + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if all(isinstance(t, ScaledTransformation) for t in transform.transformations): + lower = transform.to_search(self.properties["loc"]) + upper = transform.to_search( + self.properties["loc"] + self.properties["scale"] + ) + return MultivariateUniform(bounds=[lower, upper]) + return None class MultivariateGaussian(BaseMultivariateDistribution): @@ -188,10 +208,6 @@ class MultivariateGaussian(BaseMultivariateDistribution): Represents a multivariate Gaussian (normal) distribution with a given mean and covariance. - This class provides methods to calculate the probability density - function (pdf), the logarithm of the pdf, and to generate random - variates (rvs) from the distribution. - Parameters ---------- mean : numpy.ndarray @@ -200,57 +216,43 @@ class MultivariateGaussian(BaseMultivariateDistribution): The covariance matrix (Σ) of the multivariate Gaussian distribution. Note that what is called σ in 1D would be the square root of Σ here. - bounds : numpy.ndarray - Lower and upper bounds (2nd dim) of the parameters (1st dim). """ - def __init__(self, mean=None, covariance=None, bounds=None, random_state=None): - super().__init__(distribution=stats.multivariate_normal) - self.name = "MultivariateGaussian" - self.properties = {"mean": np.asarray(mean), "cov": np.asarray(covariance)} - self.sigma2 = covariance - self._n_parameters = len(mean) + def __init__(self, mean: np.ndarray, covariance: np.ndarray): + super().__init__( + distribution=stats.multivariate_normal, + properties={"mean": np.asarray(mean), "cov": np.asarray(covariance)}, + n_parameters=len(mean), + ) - @property - def compatible_transformations(self): - return (IdentityTransformation, ScaledTransformation, UnitHyperCube) + def support(self) -> np.ndarray: + """Return the support as a numpy array of dimensions (2, n_parameters).""" + return np.asarray([[-np.inf, np.inf]] * self._n_parameters).T - @property def mean(self): return self.properties["mean"] - @property - def sigma(self): - return sqrtm(self.properties["cov"]) + def cov(self): + return self.properties["cov"] - def marginal(self, position): + def marginal(self, position: int): # Note that what is called σ in 1D is the square root of Σ here - return stats.norm( - loc=self.properties["mean"][position], - scale=np.sqrt(self.properties["cov"][position, position]), + return Gaussian( + mean=self.properties["mean"][position], + sigma=np.sqrt(self.properties["cov"][position, position]), ) - def transformed_properties(self, transformation): - """Return distribution properties in search space""" - if isinstance(transformation, self.compatible_transformations) or ( - isinstance(transformation, ComposedTransformation) - and all( - isinstance(trans, self.compatible_transformations) - for trans in transformation.transformations - ) - ): - mean = transformation.to_search(self.properties["mean"]) - covariance = transformation.convert_covariance_matrix( - self.properties["cov"], mean - ) - - return {"mean": mean, "cov": covariance} - - else: - raise TypeError( - "The transformation provided is not compatible with pybop.MultivariateGaussian. " - "Only pybop.IdentityTransformation, pybop.ScaledTransformation and pybop.UnitHypercube are allowed." + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if all(isinstance(t, ScaledTransformation) for t in transform.transformations): + mean = transform.to_search(self.properties["mean"]) + return MultivariateGaussian( + mean=mean, + covariance=transform.convert_covariance_matrix( + self.properties["cov"], mean + ), ) + return None class MultivariateLogNormal(BaseMultivariateDistribution): @@ -258,10 +260,6 @@ class MultivariateLogNormal(BaseMultivariateDistribution): Represents a multivariate log-normal distribution with a given mean and covariance of the normally distributed log(x). - This class provides methods to calculate the probability density - function (pdf), the logarithm of the pdf, and to generate random - variates (rvs) from the distribution. - Parameters ---------- mean_log_x : numpy.ndarray @@ -272,39 +270,44 @@ class MultivariateLogNormal(BaseMultivariateDistribution): square root of Σ here. """ - def __init__(self, mean_log_x=None, covariance_log_x=None): - self.name = "MultivariateLogNormal" + def __init__(self, mean_log_x, covariance_log_x): + super().__init__(n_parameters=len(mean_log_x)) self.distribution_log_x = stats.multivariate_normal self.properties_log_x = { "mean": np.asarray(mean_log_x), "cov": np.asarray(covariance_log_x), } - self._n_parameters = len(mean_log_x) - mean, covariance = self._get_covariance_and_mean() - self.sigma2 = covariance - self.properties = {"mean": mean, "cov": covariance} + self.properties = self._get_properties() - @property - def compatible_transformations(self): - return (LogTransformation, IdentityTransformation) - - def _get_covariance_and_mean(self): - """Method to compute the covariance and mean in the model space based on the - covariance and mean of log(x) (i.e. in the search space)""" - covariance = np.zeros((self._n_parameters, self._n_parameters)) + def _get_properties(self): + """ + Method to compute the mean and covariance of x given the mean and + covariance of log(x). + """ mean = np.zeros(self._n_parameters) for i in range(self._n_parameters): mean[i] = np.exp( self.properties_log_x["mean"][i] + 0.5 * self.properties_log_x["cov"][i, i] ) - + covariance = np.zeros((self._n_parameters, self._n_parameters)) for i in range(self._n_parameters): for j in range(self._n_parameters): covariance[i, j] = ( (np.exp(self.properties_log_x["cov"][i, j]) - 1) * mean[i] * mean[j] ) - return mean, covariance + + return {"mean": mean, "cov": covariance} + + def support(self) -> np.ndarray: + """Return the support as a numpy array of dimensions (2, n_parameters).""" + return np.asarray([[0, np.inf]] * self._n_parameters).T + + def mean(self): + return self.properties["mean"] + + def cov(self): + return self.properties["cov"] def pdf(self, x): if np.any(x <= 0): @@ -316,7 +319,7 @@ def pdf(self, x): # get pdf of normal distribution for transformed variable mvn_pdf = self.distribution_log_x.pdf(log_x, **self.properties_log_x) - # determinant of the jacobion of the transformation + # determinant of the jacobian of the transformation jacobian = np.prod(x) return mvn_pdf / jacobian @@ -344,7 +347,7 @@ def cdf(self, x): # use normal distribution to get cdf of transformed variable return self.distribution_log_x.cdf(log_x, **self.properties_log_x) - def rvs(self, size=1, random_state=None): + def rvs(self, size: int = 1, random_state: int | None = None): # use normal distribution of log(x) to compute rvs return np.exp( self.distribution_log_x.rvs( @@ -352,39 +355,21 @@ def rvs(self, size=1, random_state=None): ) ) - def marginal(self, position): + def marginal(self, position: int): # determine marginal distribution based on properties of distribution of log(x) - return stats.lognorm( - scale=np.exp(self.properties_log_x["mean"][position]), - s=np.sqrt(self.properties_log_x["cov"][position, position]), + return LogNormal( + mean_log_x=self.properties_log_x["mean"][position], + sigma=np.sqrt(self.properties_log_x["cov"][position, position]), ) - def transformed_properties(self, transformation): - """Return distribution properties in search space""" - - # only allow Indentity- and LogTransformation and use properties or properties_log_x to return correct properties - if isinstance(transformation, LogTransformation) or ( - isinstance(transformation, ComposedTransformation) - and all( - isinstance(trans, LogTransformation) - for trans in transformation.transformations - ) - ): - return self.properties_log_x - elif isinstance(transformation, IdentityTransformation) or ( - isinstance(transformation, ComposedTransformation) - and all( - isinstance(trans, IdentityTransformation) - for trans in transformation.transformations - ) - ): - return self.properties - - else: - raise TypeError( - "The transformation provided is not compatible with pybop.MultivariateLogNormal. " - "Only pybop.LogTransformation and pybop.IdentityTransformation are allowed." + def _transform(self, transform: Transformation) -> Distribution | None: + """Get the transformed distribution in the search space.""" + if all(isinstance(t, LogTransformation) for t in transform.transformations): + return MultivariateGaussian( + mean=self.properties_log_x["mean"], + covariance=self.properties_log_x["cov"], ) + return None class MarginalDistribution(Distribution): @@ -414,22 +399,18 @@ def __init__( self._position = position self.parent_distribution = parent_distribution - def mean(self): - if isinstance(self.parent_distribution, MultivariateNonparametric): - # scipy.stats.gaussian_kde does not have the method "mean" - # compute mean based on dataset - return np.mean(self.distribution.dataset) - else: - return self.distribution.mean() - - def std(self): - if isinstance(self.parent_distribution, MultivariateNonparametric): - # scipy.stats.gaussian_kde does not have the method "std" - # compute standard deviation based on dataset - return np.std(self.distribution.dataset, ddof=1) - else: - return self.distribution.std() - @property - def position(self): + def position(self) -> int: return self._position + + def get_transformed_distribution( + self, transform: Transformation + ) -> Distribution | None: + """ + Get the transformed distribution in the search space, by first transforming the + parent distribution and then fetching the marginal distribution. + """ + transformed_parent_distribution = ( + self.parent_distribution.get_transformed_distribution(transform) + ) + return MarginalDistribution(transformed_parent_distribution, self._position) diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index c3526ded3..1a2cff038 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -3,14 +3,17 @@ from collections import OrderedDict from collections.abc import Iterator, Sequence from copy import deepcopy -from dataclasses import dataclass from typing import Any import numpy as np -import scipy.stats as stats from numpy.typing import NDArray -from pybop.parameters.distributions import Distribution, JointDistribution +from pybop.parameters.distributions import ( + Distribution, + JointDistribution, + Unbounded, + Uniform, +) from pybop.parameters.multivariate_distributions import ( BaseMultivariateDistribution, MarginalDistribution, @@ -19,7 +22,6 @@ from pybop.transformation.transformations import ( ComposedTransformation, IdentityTransformation, - LogTransformation, ) # Type aliases @@ -47,61 +49,17 @@ class ParameterNotFoundError(ParameterError): pass -@dataclass(frozen=True) -class Bounds: - """ - Immutable bounds representation with validation. - - Attributes - ---------- - lower : float - Lower bound (inclusive) - upper : float - Upper bound (inclusive) - """ - - lower: float - upper: float - - def __post_init__(self) -> None: - if self.lower >= self.upper: - raise ParameterValidationError( - f"Lower bound ({self.lower}) must be less than upper bound ({self.upper})" - ) - - def contains(self, value: NumericValue) -> bool: - """Check if value is within bounds.""" - return self.lower <= value <= self.upper - - def contains_array(self, values: ArrayLike) -> bool: - """Check if all values in array are within bounds.""" - arr = np.asarray(values) - return bool(np.all((arr >= self.lower) & (arr <= self.upper))) - - def clip(self, value: NumericValue) -> float: - """Clip value to bounds.""" - return float(np.clip(value, self.lower, self.upper)) - - def clip_array(self, values: ArrayLike) -> NDArray[np.floating]: - """Clip array values to bounds.""" - return np.clip(values, self.lower, self.upper) - - def width(self) -> float: - """Return the width of the bounds.""" - return self.upper - self.lower - - class Parameter: """ Represents a parameter within the PyBOP framework. This class encapsulates the definition of a parameter, including its - initial value, bounds. + distribution, initial value and transformation. Parameters ---------- - distribution : stats.distribution.rv_frozen | Distribution, optional - Distribution of the parameter + distribution : Distribution, optional + Probability distribution for the parameter. bounds : tuple[float, float], optional Parameter bounds as (lower, upper) initial_value : NumericValue, optional @@ -112,120 +70,78 @@ class Parameter: def __init__( self, - distribution: stats.distributions.rv_frozen | Distribution | None = None, + distribution: Distribution | None = None, bounds: BoundsPair | None = None, initial_value: float = None, transformation: Transformation | None = None, - ) -> None: + ): self._distribution = distribution - self._bounds = None self._initial_value = None self._transformation = transformation or IdentityTransformation() - - # Some optimisers (EP-BOLFI) and all samplers require distribution properties in the search space - # rather than the model space. Some transformations are not suitable for some multivariate - # distributions as they are currently implemented in this context - self._check_compatible_transformation() - - if self._distribution is not None: - lower, upper = self._distribution.support() - if np.isinf(lower) and np.isinf(upper): - self._bounds = None - else: - self._bounds = Bounds(float(lower), float(upper)) + self._transformed_distribution = None if bounds is not None: - if distribution is not None: + if self._distribution is not None: raise ParameterError( "Bounds can only be set if no distribution is provided. If a bounded distribution " "is needed, please ensure the distribution itself is bounded." ) - # Set bounds with validation - self._bounds = Bounds(float(bounds[0]), float(bounds[1])) - # Add uniform distribution for finite bounds in order to sample initial values - if all(np.isfinite(np.asarray(bounds))): - self._distribution = stats.uniform( - loc=bounds[0], scale=bounds[1] - bounds[0] - ) - - # Set and validate initial value - self.update_initial_value(value=initial_value) - - def _check_compatible_transformation(self) -> None: - """Raise an error if the transformation is not compatible with the distribution.""" - if isinstance(self._distribution, MarginalDistribution): - allowed_transformations = ( - self._distribution.parent_distribution.compatible_transformations - ) - if not isinstance(self._transformation, allowed_transformations): - raise TypeError( - "The transformation provided is not compatible with " - f"pybop.{self._distribution.parent_distribution.name}. Only " - + ", ".join([trans.__name__ for trans in allowed_transformations]) - + " are allowed." + # Validate bounds + if bounds[0] >= bounds[1]: + raise ParameterValidationError( + f"Lower bound ({bounds[0]}) must be less than upper bound ({bounds[1]})" ) - def sample_from_distribution( - self, - n_samples: int = 1, - *, - random_state: int | None = None, - transformed: bool = False, - ) -> NDArray[np.floating] | None: - """ - Sample from parameter's distribution. - - Parameters - ---------- - n_samples : int - Number of samples to draw (default: 1). - random_state : int, optional - Random seed for reproducibility. - transformed : bool, optional - Whether to apply transformation to samples (default: False). + # Use a uniform or unbounded distribution to represent the bounds + if all(np.isfinite(np.asarray(bounds))): + self._distribution = Uniform(lower=bounds[0], upper=bounds[1]) + else: + self._distribution = Unbounded( + initial_value=initial_value, lower=bounds[0], upper=bounds[1] + ) - Returns - ------- - NDArray[np.floating] or None - Array of samples, or None if no distribution exists exists - """ if self._distribution is None: - return None - - samples = self._distribution.rvs(n_samples, random_state=random_state) - samples = np.atleast_1d(samples).astype(float) + if initial_value is not None: + self._distribution = Unbounded(initial_value=initial_value) + else: + self._distribution = Distribution() - if transformed: - samples = np.array([self._transformation.to_search(s)[0] for s in samples]) + # Set and validate initial value + self.update_initial_value(value=initial_value) - return samples + # Set the transformed distribution + self._transformed_distribution = ( + self._distribution.get_transformed_distribution(self._transformation) + ) def update_initial_value(self, value: NumericValue | None) -> None: """Update the initial parameter value.""" self._initial_value = float(value) if value is not None else None - self._validate_initial_value_within_bounds() + + if not self.value_within_bounds(self._initial_value): + raise ParameterValidationError( + f"Initial value {value} is outside the parameter bounds {self.bounds}." + ) def __repr__(self) -> str: """String representation of the parameter.""" return f"Parameter - Distribution: {self._distribution}, Bounds: {self.bounds}, Initial value: {self._initial_value}" - def _validate_initial_value_within_bounds(self) -> None: - """Validate that the initial value is within the bounds.""" - if self._initial_value is not None: - if self._bounds is None: - return + def value_within_bounds(self, value: float = None) -> bool: + """Check if the value is within the bounds.""" + if value is None or self.bounds is None: + return True - if not self._bounds.contains(self._initial_value): - raise ParameterValidationError( - f"Initial value {self._initial_value} is outside bounds {self.bounds}" - ) + if self.bounds[0] <= value <= self.bounds[1]: + return True + return False def get_initial_value(self, transformed: bool = False) -> NDArray | None: """Get initial value in either the model space or the transformed search space.""" if self._initial_value is None and self._distribution is not None: # Try to sample from distribution if available - self.update_initial_value(self.sample_from_distribution(1)[0]) + self.update_initial_value(self._distribution.rvs(1)[0]) if self._initial_value is None: # If still None, just return this @@ -237,35 +153,13 @@ def get_initial_value(self, transformed: bool = False) -> NDArray | None: def get_mean(self, transformed: bool = False): """Get the mean of each parameter, or its initial value.""" - if self.distribution is not None: - mean = self.distribution.mean() - elif self.bounds is not None and np.isfinite(self.bounds[1] - self.bounds[0]): - lower, upper = self.bounds - mean = (lower + upper) / 2 - else: - mean = self.get_initial_value() - - if transformed and mean is not None: - mean = self.transformation.to_search(np.asarray(mean)) - - return np.asarray(mean).item() + dist = self._transformed_distribution if transformed else self._distribution + return dist.mean() def get_std(self, transformed: bool = False): """Get the standard deviation, or an estimate of it.""" - if self.distribution is not None: - std = self.distribution.std() - elif self.bounds is not None and np.isfinite(self.bounds[1] - self.bounds[0]): - lower, upper = self.bounds - std = 0.05 * (upper - lower) - else: - std = 0.05 * self.get_initial_value() - - if transformed and std is not None: - std = self.transformation.convert_standard_deviation( - std, self.get_mean(transformed=True) - ) - - return np.asarray(std).item() + dist = self._transformed_distribution if transformed else self._distribution + return dist.std() def __call__(self, *unused_args, **unused_kwargs) -> float | None: """Return the initial value. The unused arguments are to pass pybamm.ParameterValues checks.""" @@ -276,20 +170,27 @@ def initial_value(self) -> float | None: return self._initial_value @property - def bounds(self) -> tuple | None: + def bounds(self) -> tuple[float] | None: """Parameter bounds as (lower, upper) tuple.""" - return ( - None if self._bounds is None else (self._bounds.lower, self._bounds.upper) - ) + lower, upper = self._distribution.support() + + if np.isinf(lower) and np.isinf(upper): + return None + else: + return (lower, upper) @property - def distribution(self) -> Any | None: + def distribution(self) -> Distribution: return self._distribution @property def transformation(self) -> Transformation: return self._transformation + @property + def transformed_distribution(self) -> Distribution: + return self._transformed_distribution + class Parameters: """ @@ -306,24 +207,23 @@ def __init__(self, parameters: dict | Parameters | None = None) -> None: ) self._parameters = None self._distribution = None - self._multivariate = None - self._transform = None + self._transformation = None + self._transformed_distribution = None self._collect_parameters(parameters) - self._transform = self.construct_transformation() - self.update_distribution() def _collect_parameters(self, parameters): parameters = parameters or {} self._parameters = OrderedDict() for name, param in parameters.items(): - self._add(name, param, update_transform=False, update_distribution=False) + self.add(name, param, update_attributes=False) + self._update_attributes() def __getitem__(self, name: str) -> Parameter: return self.get(name) def __setitem__(self, name: str, param: Parameter) -> None: - self.set(name, param) + self.set(name, param, update_attributes=True) def __len__(self) -> int: return len(self._parameters) @@ -339,13 +239,7 @@ def names(self) -> list[str]: def __iter__(self) -> Iterator[Parameter]: return iter(self._parameters.values()) - def add( - self, name: str, parameter: Parameter, update_distribution: bool = True - ) -> None: - """Add a parameter to the collection.""" - self._add(name, parameter, update_distribution=update_distribution) - - def update_distribution(self): + def _update_attributes(self): """ Method to determine whether to construct a JointDistribution or a MultivariateDistribution and to set up the distribution. @@ -354,14 +248,16 @@ def update_distribution(self): marginal distribution. The pybop.MarginalDistribution class retains the underlying pybop.MultivariateDistribution in the parent_distribution property. """ + self._transformation = self.construct_transformation() + # check if any distribution is a pybop.MarginalDistribution - self._multivariate = any( + multivariate = any( isinstance(param.distribution, MarginalDistribution) for param in self ) # if there is a pybop.MarginalDistribution ensure all distributions are marginal # distributions of the same parent_distribution - if self._multivariate: + if multivariate: if not all( isinstance(param.distribution, MarginalDistribution) for param in self ): @@ -384,27 +280,24 @@ def update_distribution(self): index = np.argsort([p.distribution.position for p in self.__iter__()]) parameter_order = [parameter_list[i] for i in index] self._parameters = {key: self._parameters[key] for key in parameter_order} - self._transform._transformations = [ # noqa: SLF001 - self._transform.transformations[i] for i in index - ] + self._transformation = ComposedTransformation( + [self._transformation.transformations[i] for i in index] + ) else: list_of_distributions = [ - param.distribution - for param in self._parameters.values() - if param.distribution is not None + param.distribution for param in self._parameters.values() ] - if len(list_of_distributions) == len(self): + if len(list_of_distributions) > 0: self._distribution = JointDistribution(*list_of_distributions) - else: - self._distribution = None - def _add( - self, - name: str, - parameter: Parameter, - update_transform: bool = True, - update_distribution: bool = True, + if self._distribution is not None: + self._transformed_distribution = ( + self._distribution.get_transformed_distribution(self._transformation) + ) + + def add( + self, name: str, parameter: Parameter, update_attributes: bool = True ) -> None: """ Internal method to add a parameter to the collection. @@ -415,10 +308,8 @@ def _add( Name of the parameter. parameter : pybop.Parameter The parameter to add. - update_transform : bool, optional - Whether to update the transformation after adding (default: True). - update_distribution : bool, optional - Whether to update the joint or multivariate distribution (default: True). + update_attributes : bool, optional + Whether to update the transformation and distributions after adding (default: True). """ if not isinstance(parameter, Parameter): raise TypeError("Expected Parameter instance") @@ -427,13 +318,9 @@ def _add( raise ParameterError(f"Parameter for '{name}' already exists") self._parameters[name] = parameter + self._update_attributes() - if update_transform: - self._transform = self.construct_transformation() - if update_distribution: - self.update_distribution() - - def join(self, parameters=None): + def join(self, parameters) -> None: """ Join two Parameters objects into the first by copying across each Parameter. @@ -444,12 +331,12 @@ def join(self, parameters=None): for name, param in parameters.items(): if name not in self._parameters.keys(): self.add( - name, param, update_distribution=False + name, param, update_attributes=False ) # don't update every time else: print(f"Discarding duplicate {name}.") - self.update_distribution() # update once when all parameters are added + self._update_attributes() # update once when all parameters are added def get(self, name: str) -> Parameter: """Get a parameter by name.""" @@ -457,17 +344,15 @@ def get(self, name: str) -> Parameter: raise ParameterNotFoundError(f"Parameter for '{name}' not found") return self._parameters[name] - def set( - self, name: str, param: Parameter, update_distribution: bool = True - ) -> None: + def set(self, name: str, param: Parameter, update_attributes: bool = True) -> None: """Set 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({"Parameter must be of type pybop.Parameter"}) self._parameters[name] = param - if update_distribution: - self.update_distribution() + if update_attributes: + self._update_attributes() def get_bounds(self, transformed: bool = False) -> dict: """ @@ -478,27 +363,9 @@ def get_bounds(self, transformed: bool = False) -> dict: transformed : bool If True, the transformation is applied to the output (default: False). """ - bounds = {"lower": [], "upper": []} - for param in self._parameters.values(): - lower, upper = param.bounds or (-np.inf, np.inf) - - if transformed and param.bounds is not None: - if isinstance(param.transformation, LogTransformation) and lower == 0: - bound_one = -np.inf - else: - bound_one = float(param.transformation.to_search(lower)[0]) - bound_two = float(param.transformation.to_search(upper)[0]) + bounds_array = self.get_bounds_array(transformed=transformed).T - if np.isnan(bound_one) or np.isnan(bound_two): - raise ValueError("Transformed bounds resulted in NaN values.") - - lower = np.minimum(bound_one, bound_two) - upper = np.maximum(bound_one, bound_two) - - bounds["lower"].append(lower) - bounds["upper"].append(upper) - - return bounds + return {"lower": bounds_array[0], "upper": bounds_array[1]} def get_bounds_array(self, transformed: bool = False) -> np.ndarray: """ @@ -509,12 +376,11 @@ def get_bounds_array(self, transformed: bool = False) -> np.ndarray: bounds : numpy.ndarray An array of shape (n_parameters, 2) containing the bounds for each parameter. """ - bounds = self.get_bounds(transformed=transformed) - return np.column_stack([bounds["lower"], bounds["upper"]]) + dist = self._transformed_distribution if transformed else self._distribution + return dist.support().T def update( self, - *, initial_values: ArrayLike | Inputs | None = None, **individual_updates: dict[str, Any], ) -> None: @@ -563,7 +429,6 @@ def _bulk_update_initial_values(self, values: ArrayLike | Inputs) -> None: def sample_from_distribution( self, n_samples: int = 1, - *, random_state: int | None = None, transformed: bool = False, ) -> NDArray[np.floating] | None: @@ -584,16 +449,10 @@ def sample_from_distribution( NDArray[np.floating] or None Array of shape (n_samples, n_parameters) or None if any distribution is missing """ - if self._distribution is None: - return None - - samples = self._distribution.rvs(n_samples, random_state=random_state) - samples = np.atleast_2d(samples) + dist = self._transformed_distribution if transformed else self._distribution - if transformed: - samples = np.asarray([self.transformation.to_search(s) for s in samples]) - - return samples + samples = dist.rvs(n_samples, random_state=random_state) + return np.atleast_2d(samples) def get_mean(self, transformed: bool = False): """ @@ -604,19 +463,10 @@ def get_mean(self, transformed: bool = False): transformed : bool, optional If True, the transformation is applied to the output (default: False). """ - if self._multivariate: - if transformed: - return self.transformed_distribution_properties["mean"] - else: - return self.distribution.properties["mean"] - - else: - means = [] - for param in self._parameters.values(): - means.append(param.get_mean(transformed=transformed)) - return np.asarray(means).T + dist = self._transformed_distribution if transformed else self._distribution + return dist.mean() - def get_std(self, transformed: bool = False) -> list: + def get_std(self, transformed: bool = False) -> np.ndarray: """ Get the standard deviation, or an estimate of it, for each parameter. @@ -625,11 +475,8 @@ def get_std(self, transformed: bool = False) -> list: transformed : bool, optional If True, the transformation is applied to the output (default: False). """ - standard_deviations = [] - - for param in self._parameters.values(): - standard_deviations.append(param.get_std(transformed=transformed)) - return standard_deviations + dist = self._transformed_distribution if transformed else self._distribution + return dist.std() def get_covariance(self, transformed: bool = False): """ @@ -640,15 +487,8 @@ def get_covariance(self, transformed: bool = False): transformed : bool, optional If True, the transformation is applied to the output (default: False). """ - if self._multivariate: - if transformed: - return self.transformed_distribution_properties["cov"] - else: - return self.distribution.properties["cov"] - - else: - standard_deviations = self.get_std(transformed=transformed) - return (np.eye(len(self)) * np.asarray(standard_deviations)) ** 2 + dist = self._transformed_distribution if transformed else self._distribution + return dist.cov() @property def distribution( @@ -657,7 +497,7 @@ def distribution( """Return the joint or multivariate distribution.""" return self._distribution - def get_initial_values(self, *, transformed: bool = False) -> NDArray[np.floating]: + def get_initial_values(self, transformed: bool = False) -> NDArray[np.floating]: """ Get initial values as array. @@ -684,7 +524,7 @@ def get_initial_values(self, *, transformed: bool = False) -> NDArray[np.floatin @property def transformation(self) -> Transformation: """Get the transformation for the parameters.""" - return self._transform + return self._transformation def construct_transformation(self) -> Transformation: """ @@ -700,15 +540,6 @@ def construct_transformation(self) -> Transformation: return ComposedTransformation(transformations) - @property - def transformed_distribution_properties(self): - # retrieve properties of the distribution in the search space - # needed for ep-bolfi optimiser - if self._multivariate: - return self.distribution.transformed_properties(self._transform) - else: - raise NotImplementedError - def get_bounds_for_plotly(self, transformed: bool = False) -> np.ndarray: """ Retrieve parameter bounds in the format expected by Plotly. @@ -718,13 +549,13 @@ def get_bounds_for_plotly(self, transformed: bool = False) -> np.ndarray: bounds : numpy.ndarray An array of shape (n_parameters, 2) containing the bounds for each parameter. """ - bounds = self.get_bounds(transformed=transformed) + bounds_array = self.get_bounds_array(transformed=transformed) # Validate that all parameters have bounds - if not np.isfinite(list(bounds.values())).all(): + if not np.isfinite(bounds_array).all(): raise ValueError("All parameters require bounds for plot.") - return np.asarray(list(bounds.values())).T + return bounds_array def to_dict(self, values: str | ArrayLike | None = None) -> Inputs: """ @@ -760,7 +591,7 @@ def verify_inputs(self, inputs: Inputs) -> bool: valid = True for name, param in self._parameters.items(): if param.bounds is not None: - if not param._bounds.contains(inputs[name]): # noqa: SLF001 + if not param.value_within_bounds(inputs[name]): valid = False return valid diff --git a/pybop/transformation/base_transformation.py b/pybop/transformation/base_transformation.py index f694e2a31..90792998b 100644 --- a/pybop/transformation/base_transformation.py +++ b/pybop/transformation/base_transformation.py @@ -122,3 +122,7 @@ def verify_input( raise ValueError(f"Transform must have {self._n_parameters} elements") return input_array + + @property + def name(self): + return self.__class__.__name__ diff --git a/pybop/transformation/transformations.py b/pybop/transformation/transformations.py index f9b921273..a1876cbb7 100644 --- a/pybop/transformation/transformations.py +++ b/pybop/transformation/transformations.py @@ -3,62 +3,6 @@ from pybop import Transformation -class IdentityTransformation(Transformation): - """ - This class implements a trivial transformation where the model parameter space - and the search space are identical. It extends the base Transformation class. - - The transformation is defined as: - - to_search: y = x - - to_model: x = y - - Key properties: - 1. Jacobian: Identity matrix - 2. Log determinant of Jacobian: Always 0 - 3. Elementwise: True (each output dimension depends only on the corresponding input dimension) - - Use cases: - 1. When no transformation is needed between spaces - 2. As a placeholder in composite transformations - 3. For testing and benchmarking other transformations - - Note: While this transformation doesn't change the parameters, it still provides - all the methods required by the Transformation interface, making it useful in - scenarios where a transformation object is expected but no actual transformation - is needed. - - Initially based on pints.IdentityTransformation method. - """ - - def __init__(self, n_parameters: int = 1): - self._n_parameters = n_parameters - - def is_elementwise(self) -> bool: - """See :meth:`Transformation.is_elementwise()`.""" - return True - - def jacobian(self, q: np.ndarray) -> np.ndarray: - """See :meth:`Transformation.jacobian()`.""" - return np.eye(self._n_parameters) - - def jacobian_S1(self, q: np.ndarray) -> tuple[np.ndarray, np.ndarray]: - """See :meth:`Transformation.jacobian_S1()`.""" - n = self._n_parameters - return self.jacobian(q), np.zeros((n, n, n)) - - def log_jacobian_det(self, q: np.ndarray) -> float: - """See :meth:`Transformation.log_jacobian_det()`.""" - return 0.0 - - def log_jacobian_det_S1(self, q: np.ndarray) -> tuple[float, np.ndarray]: - """See :meth:`Transformation.log_jacobian_det_S1()`.""" - return self.log_jacobian_det(q), np.zeros(self._n_parameters) - - def _transform(self, x: np.ndarray, method: str) -> np.ndarray: - """See :meth:`Transformation._transform`.""" - return self.verify_input(x) - - class ScaledTransformation(Transformation): """ This class implements a linear transformation between the model parameter space @@ -66,19 +10,24 @@ class ScaledTransformation(Transformation): It extends the base Transformation class. The transformation is defined as: - - to_search: y = coefficient * (x + intercept) - - to_model: x = y / coefficient - intercept + - to_search: q = coefficient * (x + intercept) + - to_model: x = q / coefficient - intercept - Where: + where: - x is in the model parameter space - - y is in the search space - - coefficient is the scaling factor - - intercept is the offset + - q is in the search space This transformation is useful for scaling and shifting parameters to a more suitable range for optimisation algorithms. Based on pints.ScaledTransformation class. + + Attributes + ---------- + coefficient : np.ndarray + The scaling coefficient for each parameter. + intercept : np.ndarray + The intercept to shift the input parameters. """ def __init__( @@ -124,6 +73,45 @@ def _transform(self, x: np.ndarray, method: str) -> np.ndarray: raise ValueError(f"Unknown method: {method}") +class IdentityTransformation(ScaledTransformation): + """ + This class implements a trivial transformation where the model parameter space + and the search space are identical. It extends the base Transformation class. + + Key properties: + 1. Jacobian: Identity matrix + 2. Log determinant of Jacobian: Always 0 + 3. Elementwise: True (each output dimension depends only on the corresponding input dimension) + + Use cases: + 1. When no transformation is needed between spaces + 2. As a placeholder in composite transformations + 3. For testing and benchmarking other transformations + + Note: While this transformation doesn't change the parameters, it still provides + all the methods required by the Transformation interface, making it useful in + scenarios where a transformation object is expected but no actual transformation + is needed. + + Initially based on pints.IdentityTransformation method. + """ + + def __init__(self, n_parameters: int = 1): + super().__init__( + coefficient=np.ones(n_parameters), + intercept=np.zeros(n_parameters), + n_parameters=n_parameters, + ) + + def log_jacobian_det(self, q: np.ndarray) -> float: + """See :meth:`Transformation.log_jacobian_det()`.""" + return 0.0 + + def _transform(self, x: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" + return self.verify_input(x) + + class UnitHyperCube(ScaledTransformation): """ A class that implements a linear transformation between the model parameter space @@ -136,41 +124,27 @@ class UnitHyperCube(ScaledTransformation): Parameters ---------- - lower : float or array-like of shape (n,) + lower : numeric or array-like of shape (n,) The lower bound(s) of the model parameter space. - upper : float or array-like of shape (n,) + upper : numeric or array-like of shape (n,) The upper bound(s) of the model parameter space. - - Attributes - ---------- - lower : np.ndarray - The lower bound of the input space for each parameter. - upper : np.ndarray - The upper bound of the input space for each parameter. - coeff : np.ndarray - The scaling coefficient (1 / (upper - lower)) for each parameter. - inter : np.ndarray - The intercept (-lower) to shift the input parameters. """ def __init__( self, - lower: float | list | np.ndarray, - upper: float | list | np.ndarray, + lower: int | float | list | np.ndarray, + upper: int | float | list | np.ndarray, ): - self.lower = lower - self.upper = upper - # Validate that upper > lower for all elements - if np.any(self.upper <= self.lower): + if np.any(upper <= lower): raise ValueError( "All elements of upper bounds must be greater than lower bounds." ) - - # Compute the scaling coefficient (1 / (upper - lower)) and intercept (-lower) - self.coeff = 1 / (self.upper - self.lower) - self.inter = -self.lower - super().__init__(coefficient=self.coeff, intercept=self.inter) + super().__init__( + coefficient=1 / (upper - lower), + intercept=-lower, + n_parameters=1 if isinstance(lower, int | float) else len(lower), + ) class LogTransformation(Transformation): @@ -179,12 +153,12 @@ class LogTransformation(Transformation): and a search space. It extends the base Transformation class. The transformation is defined as: - - to_search: y = log(x) - - to_model: x = exp(y) + - to_search: q = log(x) + - to_model: x = exp(q) Where: - x is in the model parameter space (strictly positive) - - y is in the search space (can be any real number) + - q is in the search space (can be any real number) This transformation is particularly useful for: 1. Parameters that are strictly positive and may span several orders of magnitude. diff --git a/tests/integration/test_eis_parameterisation.py b/tests/integration/test_eis_parameterisation.py index eca589588..3683ce362 100644 --- a/tests/integration/test_eis_parameterisation.py +++ b/tests/integration/test_eis_parameterisation.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -40,12 +39,10 @@ def parameter_values(self): def parameters(self): return { "Negative electrode active material volume fraction": pybop.Parameter( - distribution=stats.uniform(loc=0.3, scale=0.9 - 0.3), - initial_value=stats.uniform(loc=0.4, scale=0.75 - 0.4).rvs(), + distribution=pybop.Uniform(0.3, 0.9) ), "Positive electrode active material volume fraction": pybop.Parameter( - distribution=stats.uniform(loc=0.3, scale=0.9 - 0.3), - initial_value=stats.uniform(loc=0.4, scale=0.75 - 0.4).rvs(), + distribution=pybop.Uniform(0.3, 0.9) ), } diff --git a/tests/integration/test_half_cell_model.py b/tests/integration/test_half_cell_model.py index d727e731b..9ea80e87d 100644 --- a/tests/integration/test_half_cell_model.py +++ b/tests/integration/test_half_cell_model.py @@ -2,7 +2,6 @@ import pybamm import pytest from pybamm import Parameter -from scipy import stats import pybop @@ -64,8 +63,7 @@ def parameter_values(self, model): def parameters(self): return { "Positive electrode active material volume fraction": pybop.Parameter( - stats.uniform(0.4, 0.75 - 0.4), - # no bounds + distribution=pybop.Uniform(0.4, 0.75) ), } diff --git a/tests/integration/test_hessian.py b/tests/integration/test_hessian.py index cb53ac2a8..e3e91cf86 100644 --- a/tests/integration/test_hessian.py +++ b/tests/integration/test_hessian.py @@ -26,18 +26,10 @@ def parameters(self, request): self.ground_truth = request.param return { "R0 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian( - 0.05, - 0.01, - truncated_at=[0.02, 0.08], - ) + distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[0.02, 0.08]) ), "R1 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian( - 0.05, - 0.01, - truncated_at=[0.02, 0.08], - ) + distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[0.02, 0.08]) ), } @@ -53,13 +45,12 @@ def parameter_values(self, model, parameters): { "Open-circuit voltage [V]": model.default_parameter_values[ "Open-circuit voltage [V]" - ] + ], + "C1 [F]": 1000, + "R0 [Ohm]": self.ground_truth[0], + "R1 [Ohm]": self.ground_truth[1], } ) - parameter_values.update({"C1 [F]": 1000}) - parameter_values.update( - {"R0 [Ohm]": self.ground_truth[0], "R1 [Ohm]": self.ground_truth[1]} - ) return parameter_values @pytest.fixture diff --git a/tests/integration/test_model_experiment_changes.py b/tests/integration/test_model_experiment_changes.py index 5d5359971..21ae49b06 100644 --- a/tests/integration/test_model_experiment_changes.py +++ b/tests/integration/test_model_experiment_changes.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -19,9 +18,7 @@ class TestModelAndExperimentChanges: { "Negative particle radius [m]": pybop.Parameter( # geometric parameter distribution=pybop.Gaussian( - 6e-06, - 0.1e-6, - truncated_at=[1e-6, 9e-6], + 6e-06, 0.1e-6, truncated_at=[1e-6, 9e-6] ), initial_value=5.86e-6, ), @@ -32,9 +29,7 @@ class TestModelAndExperimentChanges: { "Positive particle diffusivity [m2.s-1]": pybop.Parameter( # non-geometric parameter distribution=pybop.Gaussian( - 3.43e-15, - 1e-15, - truncated_at=[1e-15, 5e-15], + 3.43e-15, 1e-15, truncated_at=[1e-15, 5e-15] ), initial_value=4e-15, ), @@ -158,7 +153,7 @@ def test_multi_fitting_problem(self, solver): parameter_values.update( { "Negative electrode active material volume fraction": pybop.Parameter( - distribution=stats.norm(loc=0.68, scale=0.05), + distribution=pybop.Gaussian(0.68, 0.05) ) } ) @@ -181,7 +176,7 @@ def test_multi_fitting_problem(self, solver): parameter_values.update( { "Negative electrode active material volume fraction": pybop.Parameter( - distribution=stats.norm(loc=0.68, scale=0.05), + distribution=pybop.Gaussian(0.68, 0.05) ) } ) diff --git a/tests/integration/test_monte_carlo.py b/tests/integration/test_monte_carlo.py index 244db6fd4..9e45a8616 100644 --- a/tests/integration/test_monte_carlo.py +++ b/tests/integration/test_monte_carlo.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop from pybop import ( @@ -23,6 +22,7 @@ class Test_Sampling_SPM: @pytest.fixture(autouse=True) def setup(self): + self.sigma = 1e-3 self.ground_truth = np.clip( pybop.add_noise(np.asarray([0.55, 0.55, 3e-3]), [5e-2, 5e-2, 1e-4]), a_min=[0.4, 0.4, 1e-5], @@ -51,13 +51,10 @@ def model_and_parameter_values(self): def parameters(self): return { "Negative electrode active material volume fraction": pybop.Parameter( - pybop.Gaussian(0.575, 0.05, truncated_at=[0.375, 0.725]), - initial_value=stats.uniform(0.4, 0.7 - 0.4).rvs(), + distribution=pybop.Gaussian(0.575, 0.05, truncated_at=[0.375, 0.725]) ), "Positive electrode active material volume fraction": pybop.Parameter( - stats.norm(loc=0.525, scale=0.05), - initial_value=stats.uniform(0.4, 0.7 - 0.4).rvs(), - # no bounds + distribution=pybop.Gaussian(0.525, 0.05) # no bounds ), } @@ -72,16 +69,13 @@ def log_pdf(self, model_and_parameter_values, parameters): simulator = pybop.pybamm.Simulator( model, parameter_values=parameter_values, protocol=dataset ) - likelihood = pybop.GaussianLogLikelihood(dataset, sigma=0.002 * 1.2) + likelihood = pybop.GaussianLogLikelihood(dataset, sigma=self.sigma * 1.2) return pybop.LogPosterior(simulator, likelihood) @pytest.fixture def map_estimate(self, log_pdf): - options = pybop.PintsOptions( - max_iterations=100, - max_unchanged_iterations=35, - ) - optim = pybop.CMAES(log_pdf, options=options) + options = pybop.SciPyMinimizeOptions(maxiter=50) + optim = pybop.SciPyMinimize(log_pdf, options=options) result = optim.run() return result.x @@ -118,12 +112,15 @@ def test_sampling_spm(self, quick_sampler, log_pdf, map_estimate): x = np.mean(result.chains, axis=1) for i in range(len(x)): np.testing.assert_allclose(x[i], self.ground_truth, atol=1.6e-2) + np.testing.assert_allclose( + result.chains[i][-1], self.ground_truth, atol=2e-2 + ) def get_data(self, model, parameter_values): experiment = pybamm.Experiment( [ - "Discharge at 0.5C for 4 minutes (12 second period)", - "Charge at 0.5C for 4 minutes (12 second period)", + "Discharge at 2C for 2 minutes (12 second period)", + "Charge at 2C for 2 minutes (12 second period)", ] ) solution = pybamm.Simulation( @@ -133,6 +130,8 @@ def get_data(self, model, parameter_values): { "Time [s]": solution["Time [s]"].data, "Current [A]": solution["Current [A]"].data, - "Voltage [V]": pybop.add_noise(solution["Voltage [V]"].data, 0.002), + "Voltage [V]": pybop.add_noise( + solution["Voltage [V]"].data, self.sigma + ), } ) diff --git a/tests/integration/test_monte_carlo_thevenin.py b/tests/integration/test_monte_carlo_thevenin.py index cbeffefc7..6c5a49dfd 100644 --- a/tests/integration/test_monte_carlo_thevenin.py +++ b/tests/integration/test_monte_carlo_thevenin.py @@ -3,12 +3,10 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop from pybop import ( MALAMCMC, - DramACMC, HamiltonianMCMC, MonomialGammaHamiltonianMCMC, RaoBlackwellACMC, @@ -31,13 +29,6 @@ def setup(self): self.ground_truth = np.clip( pybop.add_noise(np.asarray([0.05, 0.05]), 0.01), a_min=1e-4, a_max=0.1 ) - self.fast_samplers = [ - MALAMCMC, - RaoBlackwellACMC, - SliceDoublingMCMC, - SliceStepoutMCMC, - DramACMC, - ] @pytest.fixture def model(self): @@ -51,12 +42,8 @@ def parameter_values(self, model): { "Open-circuit voltage [V]": model.default_parameter_values[ "Open-circuit voltage [V]" - ] - } - ) - parameter_values.update( - { - "C1 [F]": 1000, + ], + "C1 [F]": 50 / self.ground_truth[1], "R0 [Ohm]": self.ground_truth[0], "R1 [Ohm]": self.ground_truth[1], } @@ -67,14 +54,12 @@ def parameter_values(self, model): def parameters(self): return { "R0 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(5e-2, 5e-3, truncated_at=[1e-4, 1e-1]), + distribution=pybop.LogNormal(mean_log_x=np.log(0.05), sigma=0.02), transformation=pybop.LogTransformation(), - initial_value=stats.uniform(2e-3, 8e-2 - 2e-3).rvs(), ), "R1 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(5e-2, 5e-3, truncated_at=[1e-4, 1e-1]), + distribution=pybop.LogNormal(mean_log_x=np.log(0.05), sigma=0.02), transformation=pybop.LogTransformation(), - initial_value=stats.uniform(2e-3, 8e-2 - 2e-3).rvs(), ), } @@ -93,8 +78,8 @@ def log_pdf(self, model, parameter_values, parameters): @pytest.fixture def map_estimate(self, log_pdf): - options = pybop.PintsOptions(max_iterations=80) - optim = pybop.CMAES(log_pdf, options=options) + options = pybop.SciPyMinimizeOptions(maxiter=50) + optim = pybop.SciPyMinimize(log_pdf, options=options) result = optim.run() return result.x @@ -142,7 +127,10 @@ def test_sampling_thevenin(self, sampler, log_pdf, map_estimate): def get_data(self, model, parameter_values): experiment = pybamm.Experiment( - ["Discharge at 0.5C for 3 minutes (20 second period)"] + [ + "Discharge at 2C for 2 minutes (12 seconds period)", + "Rest for 20 seconds (4 seconds period)", + ] ) solution = pybamm.Simulation( model, parameter_values=parameter_values, experiment=experiment diff --git a/tests/integration/test_optimisation_options.py b/tests/integration/test_optimisation_options.py index 4d3d252e7..1e41d2d51 100644 --- a/tests/integration/test_optimisation_options.py +++ b/tests/integration/test_optimisation_options.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -42,8 +41,7 @@ def parameters(self): distribution=pybop.Gaussian(0.55, 0.05, truncated_at=[0.375, 0.75]), ), "Positive electrode active material volume fraction": pybop.Parameter( - stats.norm(loc=0.55, scale=0.05), - # no bounds + distribution=pybop.Gaussian(0.55, 0.05) # no bounds ), } diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 9c390ac2e..8ddec01fa 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -42,25 +41,10 @@ def model_and_parameter_values(self): def parameters(self): return { "Negative electrode active material volume fraction": pybop.Parameter( - stats.uniform(0.3, 0.9 - 0.3), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(0.3, 0.9) ), "Positive electrode active material volume fraction": pybop.Parameter( - stats.uniform(0.3, 0.9 - 0.3), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), - ), - } - - @pytest.fixture - def priors(self): - return { - "Negative electrode active material volume fraction": pybop.Parameter( - pybop.Uniform(0.3, 0.9), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), - ), - "Positive electrode active material volume fraction": pybop.Parameter( - pybop.Uniform(0.3, 0.9), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(0.3, 0.9) ), } @@ -142,18 +126,12 @@ def optim(self, optimiser, model_and_parameter_values, parameters, cost_class): ) problem.parameters["Negative electrode active material volume fraction"] = ( pybop.Parameter( - stats.uniform( - bounds["lower"][0], bounds["upper"][0] - bounds["lower"][0] - ), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(bounds["lower"][0], bounds["upper"][0]) ) ) problem.parameters["Positive electrode active material volume fraction"] = ( pybop.Parameter( - stats.uniform( - bounds["lower"][1], bounds["upper"][1] - bounds["lower"][1] - ), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(bounds["lower"][1], bounds["upper"][1]) ) ) @@ -258,18 +236,12 @@ def test_multiple_signals(self, multi_optimiser, two_signal_problem): two_signal_problem.parameters[ "Negative electrode active material volume fraction" ] = pybop.Parameter( - stats.uniform( - bounds["lower"][0], bounds["upper"][0] - bounds["lower"][0] - ), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(bounds["lower"][0], bounds["upper"][0]) ) two_signal_problem.parameters[ "Positive electrode active material volume fraction" ] = pybop.Parameter( - stats.uniform( - bounds["lower"][1], bounds["upper"][1] - bounds["lower"][1] - ), - initial_value=stats.uniform(0.4, 0.75 - 0.4).rvs(), + distribution=pybop.Uniform(bounds["lower"][1], bounds["upper"][1]) ) # Test each optimiser diff --git a/tests/integration/test_thevenin_parameterisation.py b/tests/integration/test_thevenin_parameterisation.py index a626329b1..1230186c8 100644 --- a/tests/integration/test_thevenin_parameterisation.py +++ b/tests/integration/test_thevenin_parameterisation.py @@ -17,7 +17,7 @@ class TestTheveninParameterisation: @pytest.fixture(autouse=True) def setup(self): self.ground_truth = np.clip( - pybop.add_noise(np.asarray([0.05, 0.05]), 0.01), a_min=0.0, a_max=0.1 + pybop.add_noise(np.asarray([0.05, 0.05]), 0.01), a_min=1e-4, a_max=0.1 ) @pytest.fixture @@ -32,12 +32,8 @@ def parameter_values(self, model): { "Open-circuit voltage [V]": model.default_parameter_values[ "Open-circuit voltage [V]" - ] - } - ) - parameter_values.update( - { - "C1 [F]": 1000, + ], + "C1 [F]": 50 / self.ground_truth[1], "R0 [Ohm]": self.ground_truth[0], "R1 [Ohm]": self.ground_truth[1], } @@ -48,12 +44,12 @@ def parameter_values(self, model): def parameters(self): return { "R0 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[1e-6, 0.1]), - transformation=pybop.LogTransformation(), + distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[0.0, 0.1]), + transformation=pybop.UnitHyperCube(0.04, 0.05), ), "R1 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[1e-6, 0.1]), - transformation=pybop.LogTransformation(), + distribution=pybop.Gaussian(0.05, 0.01, truncated_at=[0.0, 0.1]), + transformation=pybop.UnitHyperCube(0.04, 0.05), ), } @@ -61,10 +57,6 @@ def parameters(self): def dataset(self, model, parameter_values): return self.get_data(model, parameter_values) - @pytest.mark.parametrize( - "cost_class", - [pybop.RootMeanSquaredError, pybop.SumSquaredError], - ) @pytest.mark.parametrize( "optimiser, method", [ @@ -76,30 +68,21 @@ def dataset(self, model, parameter_values): ], ) def test_optimisers_on_thevenin_model( - self, - model, - parameter_values, - parameters, - dataset, - cost_class, - optimiser, - method, + self, model, parameter_values, parameters, dataset, optimiser, method ): parameter_values.update(parameters) simulator = pybop.pybamm.Simulator( model, parameter_values=parameter_values, protocol=dataset ) # Define the cost to optimise - cost = cost_class(dataset) + cost = pybop.SumSquaredError(dataset) problem = pybop.Problem(simulator, cost) x0 = problem.parameters.get_initial_values() if optimiser is pybop.SciPyMinimize: - options = pybop.SciPyMinimizeOptions(maxiter=150) + options = pybop.SciPyMinimizeOptions(maxiter=150, method=method) else: options = pybop.PintsOptions(max_iterations=150) - if optimiser is pybop.SciPyMinimize: - options.method = method if method == "L-BFGS-B": options.jac = True optim = optimiser(problem, options=options) @@ -126,9 +109,8 @@ def test_optimisers_on_thevenin_model( def get_data(self, model, parameter_values): experiment = pybamm.Experiment( [ - "Discharge at 0.5C for 6 minutes (12 seconds period)", + "Discharge at 2C for 2 minutes (12 seconds period)", "Rest for 20 seconds (4 seconds period)", - "Charge at 0.5C for 6 minutes (12 seconds period)", ] ) solution = pybamm.Simulation( diff --git a/tests/integration/test_transformation.py b/tests/integration/test_transformation.py index 707c5e3b0..93f23277a 100644 --- a/tests/integration/test_transformation.py +++ b/tests/integration/test_transformation.py @@ -31,7 +31,7 @@ class TestTransformation: def setup(self): self.sigma = 2e-3 self.ground_truth = np.clip( - pybop.add_noise(np.asarray([0.005, 0.005]), 0.001), a_min=0.0, a_max=0.01 + pybop.add_noise(np.asarray([0.005, 0.005]), 0.001), a_min=1e-4, a_max=0.01 ) @pytest.fixture @@ -46,11 +46,7 @@ def parameter_values(self, model): { "Open-circuit voltage [V]": model.default_parameter_values[ "Open-circuit voltage [V]" - ] - } - ) - parameter_values.update( - { + ], "C1 [F]": 50 / self.ground_truth[1], "R0 [Ohm]": self.ground_truth[0], "R1 [Ohm]": self.ground_truth[1], @@ -62,11 +58,13 @@ def parameter_values(self, model): def parameters(self, transformation_r0, transformation_r1): return { "R0 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(0.005, 0.002, truncated_at=[1e-6, 0.01]), + distribution=pybop.LogUniform(1e-4, 0.01), + initial_value=0.005, transformation=transformation_r0, ), "R1 [Ohm]": pybop.Parameter( - distribution=pybop.Gaussian(0.005, 0.002, truncated_at=[1e-6, 0.01]), + distribution=pybop.LogUniform(1e-4, 0.01), + initial_value=0.005, transformation=transformation_r1, ), } @@ -106,7 +104,7 @@ def problem(self, model, parameter_values, parameters, cost_class): itertools.product( [ pybop.IdentityTransformation(), - pybop.UnitHyperCube(1e-4, 0.01), + pybop.UnitHyperCube(5e-5, 0.01), pybop.LogTransformation(), ], repeat=2, diff --git a/tests/integration/test_weighted_cost.py b/tests/integration/test_weighted_cost.py index fc9185194..37f7ffdc3 100644 --- a/tests/integration/test_weighted_cost.py +++ b/tests/integration/test_weighted_cost.py @@ -2,7 +2,6 @@ import pybamm import pytest from pybamm import Parameter -from scipy import stats import pybop @@ -33,6 +32,7 @@ def model(self): @pytest.fixture def parameter_values(self): parameter_values = pybamm.ParameterValues("Chen2020") + x = self.ground_truth parameter_values.update( { "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), @@ -50,11 +50,6 @@ def parameter_values(self): ), "Cell mass [kg]": pybop.pybamm.cell_mass(), "Cell volume [m3]": pybop.pybamm.cell_volume(), - } - ) - x = self.ground_truth - parameter_values.update( - { "Negative electrode active material volume fraction": x[0], "Positive electrode active material volume fraction": x[1], } @@ -65,11 +60,10 @@ def parameter_values(self): def parameters(self): return { "Negative electrode active material volume fraction": pybop.Parameter( - stats.uniform(0.4, 0.75 - 0.4), + distribution=pybop.Uniform(0.4, 0.75) ), "Positive electrode active material volume fraction": pybop.Parameter( - stats.uniform(0.4, 0.75 - 0.4), - # no bounds + distribution=pybop.Uniform(0.4, 0.75) # no bounds ), } diff --git a/tests/unit/test_classifier.py b/tests/unit/test_classifier.py index 401b4aade..4beeef14f 100644 --- a/tests/unit/test_classifier.py +++ b/tests/unit/test_classifier.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -29,11 +28,7 @@ def problem(self): parameter_values = model.default_parameter_values parameter_values.update( - { - "R0 [Ohm]": pybop.Parameter( - distribution=stats.uniform(loc=0.001, scale=0.1 - 0.001) - ) - } + {"R0 [Ohm]": pybop.Parameter(pybop.Uniform(lower=0.001, upper=0.1))} ) simulator = pybop.pybamm.Simulator( model, parameter_values=parameter_values, protocol=dataset diff --git a/tests/unit/test_distributions.py b/tests/unit/test_distributions.py index 3c522bf11..1fe779255 100644 --- a/tests/unit/test_distributions.py +++ b/tests/unit/test_distributions.py @@ -15,7 +15,7 @@ class TestDistributions: """ - A class to test the distribution. + A class to test the distributions. """ pytestmark = pytest.mark.unit @@ -146,12 +146,15 @@ def test_distributions( JointDistribution1.logpdfS1([0.4]) # Test properties - assert Uniform.mean() == (Uniform.upper - Uniform.lower) / 2 + assert ( + Uniform.mean() + == Uniform.properties["loc"] + Uniform.properties["scale"] / 2 + ) np.testing.assert_allclose( - Uniform.std(), (Uniform.upper - Uniform.lower) / (2 * np.sqrt(3)), atol=1e-8 + Uniform.std(), Uniform.mean() / np.sqrt(3), atol=1e-8 ) - assert Exponential.mean() == Exponential.scale - assert Exponential.std() == Exponential.scale + assert Exponential.mean() == Exponential.properties["scale"] + assert Exponential.std() == Exponential.properties["scale"] def test_gaussian_rvs(self, Gaussian): samples = Gaussian.rvs(size=500) @@ -161,7 +164,7 @@ def test_gaussian_rvs(self, Gaussian): assert abs(std - 1) < 0.2 def test_incorrect_rvs(self, Gaussian): - with pytest.raises(ValueError): + with pytest.raises(TypeError): Gaussian.rvs(size="a") with pytest.raises(ValueError): Gaussian.rvs(size=(1, 2, -1)) @@ -177,12 +180,14 @@ def test_exponential_rvs(self, Exponential): assert abs(mean - 1) < 0.2 def test_repr(self, Gaussian, Uniform, Exponential, JointDistribution1): - assert repr(Gaussian) == "Gaussian, mean: 0.5, standard deviation: 1.0" - assert repr(Uniform) == "Uniform, lower: 0, upper: 1" - assert repr(Exponential) == "Exponential, loc: 0, scale: 1" + assert repr(Gaussian) == "Gaussian, properties: {'loc': 0.5, 'scale': 1.0}" + assert repr(Uniform) == "Uniform, bounds: (0.0, 1.0)" + assert ( + repr(Exponential) == "Exponential, properties: {'loc': 0.0, 'scale': 1.0}" + ) assert ( repr(JointDistribution1) - == "JointDistribution(distributions: [Gaussian, mean: 0.5, standard deviation: 1.0; Uniform, lower: 0, upper: 1])" + == "JointDistribution(distributions: [Gaussian, properties: {'loc': 0.5, 'scale': 1.0}; Uniform, bounds: (0.0, 1.0)])" ) def test_invalid_size(self, Gaussian, Uniform, Exponential): @@ -226,16 +231,21 @@ def test_multivariate_distributions( == MultivariateGaussian.rvs(10).shape[1] ) np.testing.assert_allclose( - MultivariateUniform.properties["bounds"], np.asarray([[0, 0], [1, 2]]) + MultivariateUniform.properties["loc"], np.asarray([0, 0]) ) - np.testing.assert_allclose(MultivariateGaussian.mean, np.asarray([0, 1])) np.testing.assert_allclose( - MultivariateGaussian.sigma, sqrtm(np.asarray([[0.2, 0.0], [0.0, 2.0]])) + MultivariateUniform.properties["scale"], np.asarray([1, 2]) + ) + np.testing.assert_allclose(MultivariateGaussian.mean(), np.asarray([0, 1])) + np.testing.assert_allclose( + MultivariateGaussian.std(), sqrtm(np.asarray([[0.2, 0.0], [0.0, 2.0]])) ) # Test marginal distributions of multivariate distributions - assert isinstance(MultivariateNonparametric.marginal(0), stats.gaussian_kde) - marginal = pybop.MarginalDistribution(MultivariateNonparametric, 1) + assert isinstance( + MultivariateNonparametric.marginal(0).distribution, stats.gaussian_kde + ) + marginal = MultivariateNonparametric.marginal(1) assert pytest.approx(marginal.mean()) == np.mean(random_dataset[1, :]) assert pytest.approx(marginal.std()) == np.std(random_dataset[1, :], ddof=1) @@ -251,7 +261,7 @@ def test_multivariate_uniform(self, bounds): dist = pybop.MultivariateUniform(bounds=bounds) # check bounds set correctly - support = dist.distribution.support() + support = dist.support() np.testing.assert_array_equal(support[0], bounds[0, :]) np.testing.assert_array_equal(support[1], bounds[1, :]) @@ -323,74 +333,53 @@ def test_multivariate_log_normal(self, X): def test_transformed_properties( self, transformation, MultivariateNonparametric, MultivariateGaussian ): - # mulitvariate log-normal distribution + # multivariate log-normal distribution mean = np.asarray([0, 1]) covariance = np.asarray([[0.2, 0.0], [0.0, 2.0]]) dist_lognorm = pybop.MultivariateLogNormal( mean_log_x=mean, covariance_log_x=covariance ) - if not isinstance(transformation, dist_lognorm.compatible_transformations): - with pytest.raises( - TypeError, - match="The transformation provided is not compatible with pybop.MultivariateLogNormal. " - "Only pybop.LogTransformation and pybop.IdentityTransformation are allowed.", - ): - dist_lognorm.transformed_properties(transformation) - elif isinstance(transformation, IdentityTransformation): - properties = dist_lognorm.transformed_properties(transformation) + if isinstance(transformation, IdentityTransformation): + t_dist = dist_lognorm.get_transformed_distribution(transformation) dist_m_0 = dist_lognorm.marginal(0) dist_m_1 = dist_lognorm.marginal(1) - assert pytest.approx(properties["mean"][0]) == dist_m_0.mean() - assert pytest.approx(properties["mean"][1]) == dist_m_1.mean() - assert pytest.approx(np.sqrt(properties["cov"][0, 0])) == dist_m_0.std() - assert pytest.approx(np.sqrt(properties["cov"][1, 1])) == dist_m_1.std() + assert pytest.approx(t_dist.mean()[0]) == dist_m_0.mean() + assert pytest.approx(t_dist.mean()[1]) == dist_m_1.mean() + assert pytest.approx(np.sqrt(t_dist.cov()[0, 0])) == dist_m_0.std() + assert pytest.approx(np.sqrt(t_dist.cov()[1, 1])) == dist_m_1.std() - properties = dist_lognorm.transformed_properties( + t_dist = dist_lognorm.get_transformed_distribution( ComposedTransformation([transformation, transformation]) ) - assert pytest.approx(properties["mean"][0]) == dist_m_0.mean() - assert pytest.approx(properties["mean"][1]) == dist_m_1.mean() - assert pytest.approx(np.sqrt(properties["cov"][0, 0])) == dist_m_0.std() - assert pytest.approx(np.sqrt(properties["cov"][1, 1])) == dist_m_1.std() + assert pytest.approx(t_dist.mean()[0]) == dist_m_0.mean() + assert pytest.approx(t_dist.mean()[1]) == dist_m_1.mean() + assert pytest.approx(np.sqrt(t_dist.cov()[0, 0])) == dist_m_0.std() + assert pytest.approx(np.sqrt(t_dist.cov()[1, 1])) == dist_m_1.std() with pytest.raises( TypeError, - match="The transformation provided is not compatible with pybop.MultivariateLogNormal. " - "Only pybop.LogTransformation and pybop.IdentityTransformation are allowed.", + match="All transformations must be of a similar type. Received ", ): - dist_lognorm.transformed_properties( + dist_lognorm.get_transformed_distribution( ComposedTransformation([transformation, LogTransformation()]) ) elif isinstance(transformation, LogTransformation): - properties = dist_lognorm.transformed_properties(transformation) - np.testing.assert_allclose(properties["cov"], covariance, atol=1e-7) - properties = dist_lognorm.transformed_properties( + t_dist = dist_lognorm.get_transformed_distribution(transformation) + np.testing.assert_allclose(t_dist.cov(), covariance, atol=1e-7) + t_dist = dist_lognorm.get_transformed_distribution( ComposedTransformation([transformation, transformation]) ) - np.testing.assert_allclose(properties["mean"], mean, atol=1e-7) - with pytest.raises( - TypeError, - match="The transformation provided is not compatible with pybop.MultivariateLogNormal. " - "Only pybop.LogTransformation and pybop.IdentityTransformation are allowed.", - ): - dist_lognorm.transformed_properties( - ComposedTransformation([transformation, IdentityTransformation()]) - ) + np.testing.assert_allclose(t_dist.mean(), mean, atol=1e-7) # multivariate non-parametric distribution - assert isinstance( - transformation, MultivariateNonparametric.compatible_transformations - ) - assert MultivariateNonparametric.transformed_properties(transformation) == {} + assert MultivariateNonparametric.properties == {} - # multivariate guassian - if not isinstance( - transformation, MultivariateGaussian.compatible_transformations - ): + # multivariate gaussian + if not isinstance(transformation, ScaledTransformation): with pytest.raises( TypeError, - match="The transformation provided is not compatible with pybop.MultivariateGaussian. " - "Only pybop.IdentityTransformation, pybop.ScaledTransformation and pybop.UnitHypercube are allowed.", + match="The transformation of a MultivariateGaussian distribution by a " + "LogTransformation is undefined or not yet implemented.", ): - MultivariateGaussian.transformed_properties(transformation) + MultivariateGaussian.get_transformed_distribution(transformation) diff --git a/tests/unit/test_evaluation.py b/tests/unit/test_evaluation.py index aae829e1c..24ead7651 100644 --- a/tests/unit/test_evaluation.py +++ b/tests/unit/test_evaluation.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -36,7 +35,7 @@ def parameters(self): ), ), "Positive electrode Bruggeman coefficient (electrode)": pybop.Parameter( - distribution=stats.norm(loc=1.5, scale=0.1), + distribution=pybop.LogNormal(mean_log_x=1.5, sigma=0.1), transformation=pybop.LogTransformation(), ), } diff --git a/tests/unit/test_likelihoods.py b/tests/unit/test_likelihoods.py index 841c084af..04ab7c3ef 100644 --- a/tests/unit/test_likelihoods.py +++ b/tests/unit/test_likelihoods.py @@ -1,7 +1,6 @@ import numpy as np import pybamm import pytest -from scipy import stats import pybop @@ -108,7 +107,7 @@ def test_gaussian_log_likelihood(self, simulator, dataset): assert grad_likelihood[key[1]] <= 0 # Test construction with sigma as a Parameter - sigma = pybop.Parameter(stats.uniform(loc=0.4, scale=0.6 - 0.4)) + sigma = pybop.Parameter(pybop.Uniform(lower=0.4, upper=0.6)) likelihood = pybop.GaussianLogLikelihood(dataset, sigma=sigma) # Test invalid sigma diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index e46540132..e8574a074 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -41,18 +41,10 @@ def one_parameter(self): def two_parameters(self): return { "Negative electrode active material volume fraction": pybop.Parameter( - distribution=pybop.Gaussian( - 0.6, - 0.02, - truncated_at=[0.58, 0.62], - ) + distribution=pybop.Gaussian(0.6, 0.02, truncated_at=[0.58, 0.62]) ), "Positive electrode active material volume fraction": pybop.Parameter( - distribution=pybop.Gaussian( - 0.5, - 0.05, - truncated_at=[0.48, 0.52], - ) + distribution=pybop.Gaussian(0.5, 0.05, truncated_at=[0.48, 0.52]) ), } @@ -633,7 +625,7 @@ def test_halting(self, problem, model): "Positive electrode active material volume fraction" ]._transformation = pybop.IdentityTransformation() - # Test max evalutions + # Test max evaluations options = pybop.PintsOptions(max_evaluations=1, verbose=True) optim = pybop.GradientDescent(problem, options=options) result = optim.run() @@ -740,28 +732,6 @@ def test_optimisation_result(self, result, problem): assert result.n_evaluations in result._n_evaluations assert result.x0 in result._x0 - def test_multistart_fails_without_distribution(self, model, dataset): - # parameter with inifinite bound (no distribution) - parameter_values = model.default_parameter_values - param = pybop.Parameter(bounds=(0.5, np.inf), initial_value=0.8) - parameter_values.update( - {"Positive electrode active material volume fraction": param} - ) - simulator = pybop.pybamm.Simulator( - model, parameter_values=parameter_values, protocol=dataset - ) - cost = pybop.SumSquaredError(dataset) - problem = pybop.Problem(simulator, cost) - - # Setup optimiser - options = pybop.PintsOptions(max_iterations=1, multistart=3) - optim = pybop.XNES(problem, options=options) - - with pytest.raises( - RuntimeError, match="Distributions must be provided for multi-start" - ): - optim.run() - def compare_result_data(self, result1, result2): assert result1.method_name == result2.method_name assert result1.n_runs == result2.n_runs diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py index 07d14b667..948f537c5 100644 --- a/tests/unit/test_parameters.py +++ b/tests/unit/test_parameters.py @@ -38,23 +38,19 @@ def test_parameter_construction(self, parameter): # test error if bounds and distribution with pytest.raises( ParameterError, - match="Bounds can only be set if no distribution is provided. If a bounded distribution is needed, please ensure the distribution itself is bounded.", + match="Bounds can only be set if no distribution is provided. If a bounded " + "distribution is needed, please ensure the distribution itself is bounded.", ): - pybop.Parameter(distribution=stats.norm(0.3, 0.1), bounds=(0.4, 0.8)) + pybop.Parameter( + distribution=pybop.Distribution(stats.norm(0.3, 0.1)), bounds=(0.4, 0.8) + ) def test_parameter_repr(self, parameter): assert ( repr(parameter) - == f"Parameter - Distribution: Gaussian, mean: {parameter.distribution.mean()}, standard deviation: {parameter.distribution.std()}, Bounds: (0.375, 0.7), Initial value: 0.6" + == f"Parameter - Distribution: {repr(parameter.distribution)}, Bounds: (0.375, 0.7), Initial value: 0.6" ) - def test_parameter_sampling(self, parameter): - samples = parameter.sample_from_distribution(n_samples=500) - assert (samples >= 0.375).all() and (samples <= 0.7).all() - - parameter = pybop.Parameter(bounds=(0, np.inf)) - assert parameter.sample_from_distribution() is None - def test_parameter_update(self, parameter): # Test initial value update parameter.update_initial_value(value=0.654) @@ -76,10 +72,12 @@ def test_invalid_inputs(self, parameter): ): pybop.Parameter(bounds=[0.7, 0.3]) - # Intiial value outside bounds + # Initial value outside bounds with pytest.raises( ParameterValidationError, - match=re.escape("Initial value 0.2 is outside bounds (0.3, 0.7)"), + match=re.escape( + "Initial value 0.2 is outside the parameter bounds (0.3, 0.7)." + ), ): pybop.Parameter(bounds=[0.3, 0.7], initial_value=0.2) @@ -90,12 +88,6 @@ def test_sample_initial_values(self): sample = parameter.get_initial_value() assert (sample >= 0.375) and (sample <= 0.7) - def test_get_mean(self): - # Test mean based on bounds with no distribution - param = pybop.Parameter(bounds=[0.3, 0.7]) - param._distribution = None - assert pytest.approx(param.get_mean()) == 0.5 - class TestParameters: """ @@ -131,9 +123,7 @@ def test_parameters_construction(self, name, parameter): { name: parameter, "Positive electrode active material volume fraction": pybop.Parameter( - distribution=pybop.Gaussian( - 0.6, 0.02, truncated_at=[0.375, 0.7] - ), + pybop.Gaussian(0.6, 0.02, truncated_at=[0.375, 0.7]), initial_value=0.6, ), } @@ -167,7 +157,7 @@ def test_parameters_transformation(self, name): params = pybop.Parameters( { "LogParam": pybop.Parameter( - bounds=[0, 1], + distribution=pybop.LogUniform(1e-4, 1), transformation=pybop.LogTransformation(), ), "ScaledParam": pybop.Parameter( @@ -187,24 +177,9 @@ def test_parameters_transformation(self, name): # Test transformed bounds bounds = params.get_bounds(transformed=True) - np.testing.assert_allclose(bounds["lower"], [-np.inf, 0.5, 0, -1]) + np.testing.assert_allclose(bounds["lower"], [np.log(1e-4), 0.5, 0, -1]) np.testing.assert_allclose(bounds["upper"], [np.log(1), 1.5, 1, 0]) - # Test unbounded transformation causes ValueError due to NaN - params = pybop.Parameters( - { - name: pybop.Parameter( - distribution=pybop.Gaussian(0.01, 0.2, truncated_at=[-1, 1]), - transformation=pybop.LogTransformation(), - ) - } - ) - - with pytest.raises( - ValueError, match="Transformed bounds resulted in NaN values." - ): - params.get_bounds(transformed=True) - def test_parameters_sampling(self, name, parameter): parameter._transformation = pybop.ScaledTransformation( coefficient=0.2, intercept=-1 @@ -218,25 +193,22 @@ def test_parameters_sampling(self, name, parameter): param = pybop.Parameter(initial_value=0.5) params = pybop.Parameters({name: param}) - samples = params.sample_from_distribution(n_samples=500, transformed=True) - assert samples is None + with pytest.raises(NotImplementedError): + params.sample_from_distribution(n_samples=500, transformed=True) - def test_get_sigma(self, name): - parameter = pybop.Parameter(stats.norm(loc=0.6, scale=0.02)) + def test_get_std(self, name): + parameter = pybop.Parameter(pybop.Distribution(stats.norm(loc=0.6, scale=0.02))) params = pybop.Parameters({name: parameter}) assert params.get_std() == pytest.approx([0.02]) parameter = pybop.Parameter(bounds=(0.375, 0.7)) - parameter._distribution = None params = pybop.Parameters({name: parameter}) - assert params.get_std() == [0.05 * (parameter.bounds[1] - parameter.bounds[0])] + assert params.get_std() == [parameter.distribution.std()] def test_initial_values_without_attributes(self): # Test without initial values - parameter = pybop.Parameters( - {"Negative electrode conductivity [S.m-1]": pybop.Parameter()} - ) - with pytest.raises(ParameterError, match="has no initial value"): + parameter = pybop.Parameters({"Param": pybop.Parameter()}) + with pytest.raises(NotImplementedError): parameter.get_initial_values() def test_get_initial_values_if_none(self, name, parameter): @@ -265,7 +237,8 @@ def test_parameters_repr(self, name, parameter): params = pybop.Parameters({name: parameter}) assert ( repr(params) - == f"Parameters(1):\n Negative electrode active material volume fraction: Parameter - Distribution: Gaussian, mean: {parameter.distribution.mean()}, standard deviation: {parameter.distribution.std()}, Bounds: (0.375, 0.7), Initial value: 0.6" + == "Parameters(1):\n Negative electrode active material volume fraction: Parameter - " + f"Distribution: {repr(params[name].distribution)}, Bounds: (0.375, 0.7), Initial value: 0.6" ) @@ -309,12 +282,14 @@ def multivariate_parameters(self, distribution): def test_compatible_transformation(self, distribution1): with pytest.raises( TypeError, - match="The transformation provided is not compatible with pybop.MultivariateLogNormal. Only LogTransformation, IdentityTransformation are allowed.", + match="The transformation of a MultivariateGaussian distribution by a " + "LogTransformation is undefined or not yet implemented.", ): pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 1), - initial_value=1e-15, - transformation=pybop.ScaledTransformation(0.5, 1.0), + distribution=pybop.MarginalDistribution( + pybop.MultivariateGaussian([0.2, 0.5], [[10, 0.0], [0.0, 10]]), 1 + ), + transformation=pybop.LogTransformation(), ) def test_rvs(self, distribution2): @@ -350,19 +325,16 @@ def test_get_covariance(self, distribution1): def test_input_checks_multivariate_parameters(self, distribution1): with pytest.raises( TypeError, - match="A Parameters object with a MarginalDistribution cannot be combined with parameters with other types of distributions", + match="A Parameters object with a MarginalDistribution cannot be combined with " + "parameters with other types of distributions", ): pybop.Parameters( { "Negative particle diffusivity [m2.s-1]": pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 0), - initial_value=3.9e-14, - transformation=pybop.LogTransformation(), + distribution=pybop.MarginalDistribution(distribution1, 0) ), "Positive particle diffusivity [m2.s-1]": pybop.Parameter( - initial_value=1e-15, - bounds=[1e-16, 1e-14], - transformation=pybop.LogTransformation(), + bounds=[1e-16, 1e-14] ), }, ) @@ -370,26 +342,21 @@ def test_input_checks_multivariate_parameters(self, distribution1): 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(), + bounds=[3.9e-15, 3.9e-13] ), "Positive particle diffusivity [m2.s-1]": pybop.Parameter( - initial_value=1e-15, - bounds=[1e-16, 1e-14], - transformation=pybop.LogTransformation(), + bounds=[1e-16, 1e-14] ), }, ) with pytest.raises( TypeError, - match="A Parameters object with a MarginalDistribution cannot be combined with parameters with other types of distributions", + match="A Parameters object with a MarginalDistribution cannot be combined with " + "parameters with other types of distributions", ): params["Negative particle diffusivity [m2.s-1]"] = pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 0), - initial_value=3.9e-14, - transformation=pybop.LogTransformation(), + distribution=pybop.MarginalDistribution(distribution1, 0) ) distribution2 = pybop.MultivariateUniform(np.asarray([[0, 0], [1, 2]])) @@ -400,14 +367,10 @@ def test_input_checks_multivariate_parameters(self, distribution1): pybop.Parameters( { "Negative particle diffusivity [m2.s-1]": pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 0), - initial_value=3.9e-14, - transformation=pybop.LogTransformation(), + distribution=pybop.MarginalDistribution(distribution1, 0) ), "Positive particle diffusivity [m2.s-1]": pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution2, 1), - initial_value=1e-15, - transformation=pybop.IdentityTransformation(), + distribution=pybop.MarginalDistribution(distribution2, 1) ), }, ) @@ -417,14 +380,10 @@ def test_parameter_order(self, distribution1): params = pybop.Parameters( { "Negative particle diffusivity [m2.s-1]": pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 1), - initial_value=3.9e-14, - transformation=pybop.LogTransformation(), + distribution=pybop.MarginalDistribution(distribution1, 1) ), "Positive particle diffusivity [m2.s-1]": pybop.Parameter( - distribution=pybop.MarginalDistribution(distribution1, 0), - initial_value=1e-15, - transformation=pybop.IdentityTransformation(), + distribution=pybop.MarginalDistribution(distribution1, 0) ), }, ) diff --git a/tests/unit/test_problem.py b/tests/unit/test_problem.py index 18f75bb9a..2e15d7a1d 100644 --- a/tests/unit/test_problem.py +++ b/tests/unit/test_problem.py @@ -22,18 +22,10 @@ def model(self): def parameters(self): return { "Negative particle radius [m]": pybop.Parameter( - distribution=pybop.Gaussian( - 2e-05, - 0.1e-5, - truncated_at=[1e-6, 5e-5], - ) + distribution=pybop.Gaussian(2e-05, 0.1e-5, truncated_at=[1e-6, 5e-5]) ), "Positive particle radius [m]": pybop.Parameter( - distribution=pybop.Gaussian( - 0.5e-05, - 0.1e-5, - truncated_at=[1e-6, 5e-5], - ) + distribution=pybop.Gaussian(0.5e-05, 0.1e-5, truncated_at=[1e-6, 5e-5]) ), } diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index e5bb18f4d..e54b038c6 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -21,18 +21,22 @@ def parameters(self): return pybop.Parameters( { "Identity": pybop.Parameter( + distribution=pybop.LogUniform(1e-4, 1), transformation=pybop.IdentityTransformation(), ), "Scaled": pybop.Parameter( + distribution=pybop.LogUniform(1e-4, 1), transformation=pybop.ScaledTransformation( coefficient=2.0, intercept=1 ), ), "Log": pybop.Parameter( + distribution=pybop.LogUniform(1e-4, 1), transformation=pybop.LogTransformation(), ), "UnitHyperCube": pybop.Parameter( - transformation=pybop.UnitHyperCube(10, 100) + distribution=pybop.LogUniform(1e-4, 1), + transformation=pybop.UnitHyperCube(5e-5, 1), ), } ) @@ -89,10 +93,10 @@ def test_scaled_transformation(self, parameters): def test_hypercube_transformation(self, parameters): q = np.asarray([0.5]) - coeff = 1 / (100 - 10) + coeff = 1 / (1 - 5e-5) transformation = parameters["UnitHyperCube"].transformation p = transformation.to_model(q) - assert np.allclose(p, (q / coeff) + 10) + assert np.allclose(p, (q / coeff) + 5e-5) assert transformation.n_parameters == 1 assert transformation.is_elementwise()