From 403a280c210fd37006ae165b453238a257d6a804 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 30 May 2024 17:22:31 +0100 Subject: [PATCH 01/18] feat: initial transformation class --- examples/scripts/spm_CMAES.py | 4 +- examples/standalone/cost.py | 4 +- pybop/__init__.py | 6 + pybop/costs/_likelihoods.py | 8 +- pybop/costs/base_cost.py | 16 +-- pybop/costs/design_costs.py | 12 +- pybop/costs/fitting_costs.py | 11 +- pybop/optimisers/base_optimiser.py | 4 + pybop/optimisers/base_pints_optimiser.py | 55 ++++++-- pybop/transformations/__init__.py | 157 ++++++++++++++++++++++ pybop/transformations/_transformations.py | 56 ++++++++ 11 files changed, 288 insertions(+), 45 deletions(-) create mode 100644 pybop/transformations/__init__.py create mode 100644 pybop/transformations/_transformations.py diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index b60bc019..a71de0bb 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -42,7 +42,9 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) cost = pybop.SumSquaredError(problem) -optim = pybop.CMAES(cost, max_iterations=100) +optim = pybop.CMAES( + cost, max_iterations=100, transformation=pybop.IdentityTransformation(2) +) # Run the optimisation x, final_cost = optim.run() diff --git a/examples/standalone/cost.py b/examples/standalone/cost.py index c3763ff4..7befc6f7 100644 --- a/examples/standalone/cost.py +++ b/examples/standalone/cost.py @@ -27,7 +27,7 @@ class StandaloneCost(pybop.BaseCost): Methods ------- - __call__(x, grad=None) + __call__(x) Calculate the cost for a given parameter value. """ @@ -48,7 +48,7 @@ def __init__(self, problem=None): upper=[10], ) - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the cost for a given parameter value. diff --git a/pybop/__init__.py b/pybop/__init__.py index ecd42019..b020a0c5 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -73,6 +73,12 @@ GaussianLogLikelihoodKnownSigma, ) +# +# Transformation classes +# +from .transformations import Transformation +from .transformations._transformations import IdentityTransformation + # # Dataset class # diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 91374cc0..20ef32bc 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -60,7 +60,7 @@ def __init__(self, problem, sigma): self.sigma2 = self.sigma0**-2 self._dl = np.ones(self._n_parameters) - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calls the problem.evaluate method and calculates the log-likelihood @@ -86,7 +86,7 @@ def _evaluate(self, x, grad=None): else: return np.sum(e) - def _evaluateS1(self, x, grad=None): + def _evaluateS1(self, x): """ Calls the problem.evaluateS1 method and calculates the log-likelihood @@ -120,7 +120,7 @@ def __init__(self, problem): self._logpi = -0.5 * self.n_time_data * np.log(2 * np.pi) self._dl = np.ones(self._n_parameters + self.n_outputs) - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Evaluates the Gaussian log-likelihood for the given parameters. @@ -159,7 +159,7 @@ def _evaluate(self, x, grad=None): else: return np.sum(e) - def _evaluateS1(self, x, grad=None): + def _evaluateS1(self, x): """ Calls the problem.evaluateS1 method and calculates the log-likelihood diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 025cb2e4..85eb55d9 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -52,13 +52,13 @@ def __init__(self, problem=None, sigma=None): def n_parameters(self): return self._n_parameters - def __call__(self, x, grad=None): + def __call__(self, x): """ Call the evaluate function for a given set of parameters. """ - return self.evaluate(x, grad) + return self.evaluate(x) - def evaluate(self, x, grad=None): + def evaluate(self, x): """ Call the evaluate function for a given set of parameters. @@ -66,9 +66,6 @@ def evaluate(self, x, grad=None): ---------- x : array-like The parameters for which to evaluate the cost. - grad : array-like, optional - An array to store the gradient of the cost function with respect - to the parameters. Returns ------- @@ -81,7 +78,7 @@ def evaluate(self, x, grad=None): If an error occurs during the calculation of the cost. """ try: - return self._evaluate(x, grad) + return self._evaluate(x) except NotImplementedError as e: raise e @@ -89,7 +86,7 @@ def evaluate(self, x, grad=None): except Exception as e: raise ValueError(f"Error in cost calculation: {e}") - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the cost function value for a given set of parameters. @@ -99,9 +96,6 @@ def _evaluate(self, x, grad=None): ---------- x : array-like The parameters for which to evaluate the cost. - grad : array-like, optional - An array to store the gradient of the cost function with respect - to the parameters. Returns ------- diff --git a/pybop/costs/design_costs.py b/pybop/costs/design_costs.py index f6364cdc..811f0c38 100644 --- a/pybop/costs/design_costs.py +++ b/pybop/costs/design_costs.py @@ -65,7 +65,7 @@ def update_simulation_data(self, initial_conditions): self.problem._target = {key: solution[key] for key in self.problem.signal} self.dt = solution["Time [s]"][1] - solution["Time [s]"][0] - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Computes the value of the cost function. @@ -75,8 +75,6 @@ def _evaluate(self, x, grad=None): ---------- x : array The parameter set for which to compute the cost. - grad : array, optional - Gradient information, not used in this method. Raises ------ @@ -99,7 +97,7 @@ class GravimetricEnergyDensity(DesignCost): def __init__(self, problem, update_capacity=False): super(GravimetricEnergyDensity, self).__init__(problem, update_capacity) - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Computes the cost function for the energy density. @@ -107,8 +105,6 @@ def _evaluate(self, x, grad=None): ---------- x : array The parameter set for which to compute the cost. - grad : array, optional - Gradient information, not used in this method. Returns ------- @@ -158,7 +154,7 @@ class VolumetricEnergyDensity(DesignCost): def __init__(self, problem, update_capacity=False): super(VolumetricEnergyDensity, self).__init__(problem, update_capacity) - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Computes the cost function for the energy density. @@ -166,8 +162,6 @@ def _evaluate(self, x, grad=None): ---------- x : array The parameter set for which to compute the cost. - grad : array, optional - Gradient information, not used in this method. Returns ------- diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index b7b26659..d4fd2331 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -23,7 +23,7 @@ def __init__(self, problem): # Default fail gradient self._de = 1.0 - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the root mean square error for a given set of parameters. @@ -138,7 +138,7 @@ def __init__(self, problem): # Default fail gradient self._de = 1.0 - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the sum of squared errors for a given set of parameters. @@ -146,9 +146,6 @@ def _evaluate(self, x, grad=None): ---------- x : array-like The parameters for which to evaluate the cost. - grad : array-like, optional - An array to store the gradient of the cost function with respect - to the parameters. Returns ------- @@ -236,7 +233,7 @@ def __init__(self, observer: Observer): super().__init__(problem=observer) self._observer = observer - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the observer cost for a given set of parameters. @@ -314,7 +311,7 @@ def __init__(self, problem, likelihood, sigma=None): ): raise ValueError(f"{self.likelihood} must be a subclass of BaseLikelihood") - def _evaluate(self, x, grad=None): + def _evaluate(self, x): """ Calculate the maximum a posteriori cost for a given set of parameters. diff --git a/pybop/optimisers/base_optimiser.py b/pybop/optimisers/base_optimiser.py index 713df4d4..9d4cd919 100644 --- a/pybop/optimisers/base_optimiser.py +++ b/pybop/optimisers/base_optimiser.py @@ -54,6 +54,7 @@ def __init__( self.verbose = False self.log = [] self.minimising = True + self.transformation = None self.physical_viability = False self.allow_infeasible_solutions = False self.default_max_iterations = 1000 @@ -103,6 +104,9 @@ def set_base_options(self): self.sigma0 = self.unset_options.pop("sigma0", self.sigma0) self.verbose = self.unset_options.pop("verbose", self.verbose) self.minimising = self.unset_options.pop("minimising", self.minimising) + self.transformation = self.unset_options.pop( + "transformation", self.transformation + ) if "allow_infeasible_solutions" in self.unset_options.keys(): self.set_allow_infeasible_solutions( self.unset_options.pop("allow_infeasible_solutions") diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index 543d32b6..074719e3 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -1,7 +1,7 @@ import numpy as np import pints -from pybop import BaseOptimiser +from pybop import BaseLikelihood, BaseOptimiser, Transformation class BasePintsOptimiser(BaseOptimiser): @@ -38,11 +38,11 @@ def __init__(self, cost, pints_optimiser, **optimiser_kwargs): self._evaluations = None self._iterations = None - # PyBOP doesn't currently support the PINTS transformation class - self._transformation = None - self.pints_optimiser = pints_optimiser super().__init__(cost, **optimiser_kwargs) + self.f = self.cost + if self.transformation is not None: + self.set_transformation(self.transformation) def _set_up_optimiser(self): """ @@ -145,6 +145,39 @@ def _sanitise_inputs(self): self.bounds["lower"], self.bounds["upper"] ) + def set_transformation(self, transformation: Transformation): + """ + Apply the given transformation to the optimizer's settings. + + Initial credit: Pints team + + Parameters + ---------- + transformation : pybop.Transformation + The transformation object to be applied. + """ + # Convert cost or log pdf + if isinstance(self.cost, BaseLikelihood): + self.f = transformation.convert_log_pdf(self.cost) + else: + self.f = transformation.convert_cost(self.cost) + + # Convert initial position + self.x0 = transformation.to_search(self.x0) + + # Convert sigma0, if provided + if self.sigma0 is not None: + self.sigma0 = transformation.convert_standard_deviation( + self.sigma0, self.x0 + ) + + # Convert boundaries, if provided + if self._boundaries: + self._boundaries = transformation.convert_boundaries(self._boundaries) + + # Store the transformation for later detransformation + self.transformation = transformation + def name(self): """ Provides the name of the optimisation strategy. @@ -191,12 +224,12 @@ def _run(self): if self._needs_sensitivities: def f(x): - L, dl = self.cost.evaluateS1(x) + L, dl = self.f.evaluateS1(x) return (L, dl) if self.minimising else (-L, -dl) else: - def f(x, grad=None): - return self.cost(x, grad) if self.minimising else -self.cost(x, grad) + def f(x): + return self.f(x) if self.minimising else -self.f(x) # Create evaluator object if self._parallel: @@ -316,8 +349,8 @@ def f(x, grad=None): # Show current parameters x_user = self.pints_optimiser.x_guessed() - if self._transformation is not None: - x_user = self._transformation.to_model(x_user) + if self.transformation is not None: + x_user = self.transformation.to_model(x_user) for p in x_user: print(pints.strfloat(p)) print("-" * 40) @@ -339,8 +372,8 @@ def f(x, grad=None): f = self.pints_optimiser.f_best() # Inverse transform search parameters - if self._transformation is not None: - x = self._transformation.to_model(x) + if self.transformation is not None: + x = self.transformation.to_model(x) # Store result final_cost = f if self.minimising else -f diff --git a/pybop/transformations/__init__.py b/pybop/transformations/__init__.py new file mode 100644 index 00000000..074693b9 --- /dev/null +++ b/pybop/transformations/__init__.py @@ -0,0 +1,157 @@ +from abc import ABC, abstractmethod +from typing import Tuple, Union, Sequence +import numpy as np + +from pybop import BaseCost + +class Transformation(ABC): + """ + Abstract base class for transformations between two parameter spaces: the model + parameter space and a search space. + + If `trans` is an instance of a `Transformation` class, you can apply the + transformation of a parameter vector from the model space `p` to the search + space `q` using `q = trans.to_search(p)` and the inverse using `p = trans.to_model(q)`. + + Based on pints.transformation method. + + References + ---------- + .. [1] Erik Jorgensen and Asger Roer Pedersen. "How to Obtain Those Nasty Standard Errors From Transformed Data." + http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.47.9023 + .. [2] Kaare Brandt Petersen and Michael Syskind Pedersen. "The Matrix Cookbook." 2012. + """ + + def convert_log_pdf(self, log_pdf): + """Returns a transformed log-PDF class.""" + return TransformedLogPDF(log_pdf, self) + + def convert_log_prior(self, log_prior): + """Returns a transformed log-prior class.""" + return TransformedLogPrior(log_prior, self) + + def convert_cost(self, cost): + """Returns a transformed cost class.""" + return TransformedCost(cost, self) + + def convert_boundaries(self, boundaries): + """Returns a transformed boundaries class.""" + return TransformedBoundaries(boundaries, self) + + def convert_covariance_matrix(self, covariance: np.ndarray, q: np.ndarray) -> np.ndarray: + """ + Converts a covariance matrix `covariance` from the model space to the search space + around a parameter vector `q` in the search space. + """ + jac_inv = np.linalg.pinv(self.jacobian(q)) + return jac_inv @ covariance @ jac_inv.T + + def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarray) -> np.ndarray: + """ + Converts standard deviation `std`, either a scalar or a vector, from the model space + to the search space around a parameter vector `q` in the search space. + """ + jac_inv = np.linalg.pinv(self.jacobian(q)) + cov = jac_inv @ jac_inv.T + return std * np.sqrt(np.diagonal(cov)) + + @abstractmethod + def jacobian(self, q: np.ndarray) -> np.ndarray: + """Returns the Jacobian matrix of the transformation at the parameter vector `q`.""" + pass + + def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, Sequence[np.ndarray]]: + """ + Computes the Jacobian matrix and its partial derivatives at the parameter vector `q`. + + Returns a tuple `(jacobian, jacobian_derivatives)`. + """ + raise NotImplementedError("jacobian_S1 method must be implemented if used.") + + def log_jacobian_det(self, q: np.ndarray) -> float: + """ + Returns the logarithm of the absolute value of the determinant of the Jacobian matrix + at the parameter vector `q`. + """ + return np.log(np.abs(np.linalg.det(self.jacobian(q)))) + + def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + """ + Computes the logarithm of the absolute value of the determinant of the Jacobian, + and returns it along with its partial derivatives. + """ + jacobian, jacobian_derivatives = self.jacobian_S1(q) + jacobian_inv = np.linalg.pinv(jacobian) + derivatives = np.array([np.trace(jacobian_inv @ djac) for djac in jacobian_derivatives]) + return self.log_jacobian_det(q), derivatives + + @abstractmethod + def n_parameters(self) -> int: + """Returns the dimension of the parameter space this transformation is defined over.""" + pass + + @abstractmethod + def to_model(self, q: np.ndarray) -> np.ndarray: + """Transforms a parameter vector `q` from the search space to the model space.""" + pass + + @abstractmethod + def to_search(self, p: np.ndarray) -> np.ndarray: + """Transforms a parameter vector `p` from the model space to the search space.""" + pass + + def is_elementwise(self) -> bool: + """ + Returns `True` if the transformation is element-wise, meaning it can be applied + element-by-element independently. + """ + raise NotImplementedError("is_elementwise method must be implemented if used.") + +class TransformedLogPDF(BaseCost): + """Transformed log-PDF class.""" + def __init__(self, log_pdf, transformation): + self._log_pdf = log_pdf + self._transformation = transformation + + def __call__(self, q): + p = self._transformation.to_model(q) + log_pdf = self._log_pdf(p) + + # Calculate the PDF using change of variable + # Wikipedia: https://w.wiki/UsJ + log_jacobian_det = self._transformation.log_jacobian_det(q) + return log_pdf + log_jacobian_det + + def _evaluateS1(self, x): + p = self._transformation.to_model(x) + log_pdf, log_pdf_derivatives = self._log_pdf._evaluateS1(p) + log_jacobian_det, log_jacobian_det_derivatives = self._transformation.log_jacobian_det_S1(x) + return log_pdf + log_jacobian_det, log_pdf_derivatives + log_jacobian_det_derivatives + +class TransformedLogPrior: + """Transformed log-prior class.""" + def __init__(self, log_prior, transformation): + self._log_prior = log_prior + self._transformation = transformation + + def __call__(self, q): + return self + +class TransformedCost(BaseCost): + """Transformed cost class.""" + def __init__(self, cost, transformation): + self._cost = cost + self._transformation = transformation + + def __call__(self, q ): + p = self._transformation.to_model(q) + return self._cost(p) + + def _evaluateS1(self, x): + p = self._transformation.to_model(x) + return self._cost._evaluateS1(x) + +class TransformedBoundaries: + """Transformed boundaries class.""" + def __init__(self, boundaries, transformation): + self._boundaries = boundaries diff --git a/pybop/transformations/_transformations.py b/pybop/transformations/_transformations.py new file mode 100644 index 00000000..a7a1541c --- /dev/null +++ b/pybop/transformations/_transformations.py @@ -0,0 +1,56 @@ +from typing import Sequence, Tuple, Union + +import numpy as np + +from pybop import Transformation + + +class IdentityTransformation(Transformation): + """ + Identity transformation between two parameter spaces: the model parameter space and a search space. + + Based on pints.IdentityTransformation method. + """ + + def __init__(self, n_parameters: int): + self._n_parameters = n_parameters + + def elementwise(self) -> bool: + """See :meth:`Transformation.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 n_parameters(self) -> int: + """See :meth:`Transformation.n_parameters()`.""" + return self._n_parameters + + def to_model(self, q: Union[np.ndarray, Sequence, float, int]) -> np.ndarray: + """See :meth:`Transformation.to_model()`.""" + return np.asarray(q) + + def to_search(self, p: Union[np.ndarray, Sequence, float, int]) -> np.ndarray: + """See :meth:`Transformation.to_search()`.""" + return np.asarray(p) + + +# TODO: Implement the following classes: +# class LogTransformation(Transformation): +# class LogitTransformation(Transformation): +# class ComposedTransformation(Transformation): +# class ScaledTransformation(Transformation): From 8530224c83c3c8a0ca0e13e75f648402eae32c32 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Tue, 18 Jun 2024 09:08:40 +0100 Subject: [PATCH 02/18] feat: Integrate transformations into Parameters API, adds log, scaled, and composed transformations, updts for example API usage --- examples/scripts/spm_CMAES.py | 6 +- pybop/__init__.py | 17 +- pybop/costs/base_cost.py | 22 +- pybop/optimisers/base_optimiser.py | 3 + pybop/optimisers/base_pints_optimiser.py | 40 +--- pybop/parameters/parameter.py | 52 ++++- pybop/transformations/__init__.py | 112 ++++----- pybop/transformations/_transformations.py | 265 ++++++++++++++++++++-- 8 files changed, 387 insertions(+), 130 deletions(-) diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index 43d816d5..ff4e0f3a 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -13,12 +13,14 @@ prior=pybop.Gaussian(6e-06, 0.1e-6), bounds=[1e-6, 9e-6], true_value=parameter_set["Negative particle radius [m]"], + transformation=pybop.ScaledTransformation(scale=2.0), ), pybop.Parameter( "Positive particle radius [m]", prior=pybop.Gaussian(4.5e-06, 0.1e-6), bounds=[1e-6, 9e-6], true_value=parameter_set["Positive particle radius [m]"], + transformation=pybop.LogTransformation(), ), ) @@ -42,9 +44,7 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) cost = pybop.SumSquaredError(problem) -optim = pybop.CMAES( - cost, max_iterations=100, transformation=pybop.IdentityTransformation(2) -) +optim = pybop.CMAES(cost, max_iterations=100, min_iterations=40) # Run the optimisation x, final_cost = optim.run() diff --git a/pybop/__init__.py b/pybop/__init__.py index 7d86b88e..82d527b0 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -78,6 +78,17 @@ from .problems.fitting_problem import FittingProblem from .problems.design_problem import DesignProblem +# +# Transformation classes +# +from .transformations import Transformation +from .transformations._transformations import ( + IdentityTransformation, + ScaledTransformation, + LogTransformation, + ComposedTransformation, +) + # # Cost function class # @@ -99,12 +110,6 @@ GaussianLogLikelihoodKnownSigma, ) -# -# Transformation classes -# -from .transformations import Transformation -from .transformations._transformations import IdentityTransformation - # # Optimiser class # diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 7f316a72..1d45b04b 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -1,4 +1,4 @@ -from pybop import BaseProblem +from pybop import BaseProblem, ComposedTransformation, IdentityTransformation class BaseCost: @@ -25,6 +25,7 @@ class BaseCost: def __init__(self, problem=None): self.parameters = None + self.transformation = None self.problem = problem self.x0 = None if isinstance(self.problem, BaseProblem): @@ -34,6 +35,13 @@ def __init__(self, problem=None): self.n_outputs = self.problem.n_outputs self.signal = self.problem.signal + # Construct ComposedTransformation from list of transformations + self.transformations = [ + t if t is not None else IdentityTransformation() + for t in self.parameters.get_transformations() + ] + self.transformation = ComposedTransformation(self.transformations) + @property def n_parameters(self): return len(self.parameters) @@ -42,7 +50,11 @@ def __call__(self, x): """ Call the evaluate function for a given set of parameters. """ - return self.evaluate(x) + if self.transformation: + p = self.transformation.to_model(x) + return self.evaluate(p) + else: + return self.evaluate(x) def evaluate(self, x): """ @@ -116,7 +128,11 @@ def evaluateS1(self, x): If an error occurs during the calculation of the cost or gradient. """ try: - return self._evaluateS1(x) + if self.transformation: + p = self.transformation.to_model(x) + return self._evaluateS1(p) + else: + return self._evaluateS1(x) except NotImplementedError as e: raise e diff --git a/pybop/optimisers/base_optimiser.py b/pybop/optimisers/base_optimiser.py index 7fff3bc3..921bf254 100644 --- a/pybop/optimisers/base_optimiser.py +++ b/pybop/optimisers/base_optimiser.py @@ -75,6 +75,9 @@ def __init__( # Set default initial standard deviation (for all or no parameters) self.sigma0 = cost.parameters.get_sigma0() or self.sigma0 + # Set transformation + self.transformation = cost.transformation + else: try: cost_test = cost(optimiser_kwargs.get("x0", [])) diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index 94376618..84a0ba71 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -10,7 +10,10 @@ from pints import SequentialEvaluator as PintsSequentialEvaluator from pints import strfloat as PintsStrFloat -from pybop import BaseLikelihood, BaseOptimiser, Result, Transformation +from pybop import ( + BaseOptimiser, + Result, +) class BasePintsOptimiser(BaseOptimiser): @@ -51,8 +54,6 @@ def __init__(self, cost, pints_optimiser, **optimiser_kwargs): self.pints_optimiser = pints_optimiser super().__init__(cost, **optimiser_kwargs) self.f = self.cost - if self.transformation is not None: - self.set_transformation(self.transformation) def _set_up_optimiser(self): """ @@ -154,39 +155,6 @@ def _sanitise_inputs(self): self.bounds["lower"], self.bounds["upper"] ) - def set_transformation(self, transformation: Transformation): - """ - Apply the given transformation to the optimizer's settings. - - Initial credit: Pints team - - Parameters - ---------- - transformation : pybop.Transformation - The transformation object to be applied. - """ - # Convert cost or log pdf - if isinstance(self.cost, BaseLikelihood): - self.f = transformation.convert_log_pdf(self.cost) - else: - self.f = transformation.convert_cost(self.cost) - - # Convert initial position - self.x0 = transformation.to_search(self.x0) - - # Convert sigma0, if provided - if self.sigma0 is not None: - self.sigma0 = transformation.convert_standard_deviation( - self.sigma0, self.x0 - ) - - # Convert boundaries, if provided - if self._boundaries: - self._boundaries = transformation.convert_boundaries(self._boundaries) - - # Store the transformation for later detransformation - self.transformation = transformation - def name(self): """ Provides the name of the optimisation strategy. diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 76754847..94549a31 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -32,7 +32,13 @@ class Parameter: """ def __init__( - self, name, initial_value=None, true_value=None, prior=None, bounds=None + self, + name, + initial_value=None, + true_value=None, + prior=None, + bounds=None, + transformation=None, ): """ Construct the parameter class with a name, initial value, prior, and bounds. @@ -42,8 +48,9 @@ def __init__( self.true_value = true_value self.initial_value = initial_value self.value = initial_value - self.set_bounds(bounds) + self.transformation = transformation self.margin = 1e-4 + self.set_bounds(bounds) def rvs(self, n_samples, random_state=None): """ @@ -64,6 +71,9 @@ def rvs(self, n_samples, random_state=None): """ samples = self.prior.rvs(n_samples, random_state=random_state) + if self.transformation is not None: + samples = self.transformation.to_search(samples) + # Constrain samples to be within bounds if self.bounds is not None: offset = self.margin * (self.upper_bound - self.lower_bound) @@ -142,11 +152,19 @@ def set_bounds(self, bounds=None): if bounds is not None: if bounds[0] >= bounds[1]: raise ValueError("Lower bound must be less than upper bound") + elif self.transformation is not None: + self.lower_bound = np.ndarray.item( + self.transformation.to_search(bounds[0]) + ) + self.upper_bound = np.ndarray.item( + self.transformation.to_search(bounds[1]) + ) else: self.lower_bound = bounds[0] self.upper_bound = bounds[1] - - self.bounds = bounds + self.bounds = [self.lower_bound, self.upper_bound] + else: + self.bounds = bounds class Parameters: @@ -166,6 +184,7 @@ def __init__(self, *args): self.param = OrderedDict() for param in args: self.add(param) + self.initial_value() def __getitem__(self, key: str) -> Parameter: """ @@ -309,12 +328,20 @@ def get_sigma0(self) -> List: for param in self.param.values(): if hasattr(param.prior, "sigma"): - sigma0.append(param.prior.sigma) + if param.transformation is None: + sigma0.append(param.prior.sigma) + else: + sigma0.append( + np.ndarray.item( + param.transformation.convert_standard_deviation( + param.prior.sigma, param.initial_value + ) + ) + ) else: all_have_sigma = False if not all_have_sigma: sigma0 = None - return sigma0 def initial_value(self) -> List: @@ -324,7 +351,7 @@ def initial_value(self) -> List: initial_values = [] for param in self.param.values(): - if param.initial_value is None: + if param.initial_value is None and param.prior is not None: initial_value = param.rvs(1) param.update(initial_value=initial_value[0]) initial_values.append(param.initial_value) @@ -353,6 +380,17 @@ def true_value(self) -> List: return true_values + def get_transformations(self): + """ + Get the transformations for each parameter. + """ + transformations = [] + + for param in self.param.values(): + transformations.append(param.transformation) + + return transformations + def get_bounds_for_plotly(self): """ Retrieve parameter bounds in the format expected by Plotly. diff --git a/pybop/transformations/__init__.py b/pybop/transformations/__init__.py index 074693b9..f1ad8cfc 100644 --- a/pybop/transformations/__init__.py +++ b/pybop/transformations/__init__.py @@ -1,8 +1,8 @@ from abc import ABC, abstractmethod from typing import Tuple, Union, Sequence import numpy as np +from pints import TransformedBoundaries as PintsTransformedBoundaries -from pybop import BaseCost class Transformation(ABC): """ @@ -22,18 +22,10 @@ class Transformation(ABC): .. [2] Kaare Brandt Petersen and Michael Syskind Pedersen. "The Matrix Cookbook." 2012. """ - def convert_log_pdf(self, log_pdf): - """Returns a transformed log-PDF class.""" - return TransformedLogPDF(log_pdf, self) - def convert_log_prior(self, log_prior): """Returns a transformed log-prior class.""" return TransformedLogPrior(log_prior, self) - def convert_cost(self, cost): - """Returns a transformed cost class.""" - return TransformedCost(cost, self) - def convert_boundaries(self, boundaries): """Returns a transformed boundaries class.""" return TransformedBoundaries(boundaries, self) @@ -51,6 +43,8 @@ def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarra Converts standard deviation `std`, either a scalar or a vector, from the model space to the search space around a parameter vector `q` in the search space. """ + if isinstance(q, (int, float)): + q = np.asarray([q]) jac_inv = np.linalg.pinv(self.jacobian(q)) cov = jac_inv @ jac_inv.T return std * np.sqrt(np.diagonal(cov)) @@ -64,7 +58,7 @@ def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, Sequence[np.ndarray]]: """ Computes the Jacobian matrix and its partial derivatives at the parameter vector `q`. - Returns a tuple `(jacobian, jacobian_derivatives)`. + Returns a tuple `(jacobian, hessian)`. """ raise NotImplementedError("jacobian_S1 method must be implemented if used.") @@ -80,9 +74,9 @@ def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: Computes the logarithm of the absolute value of the determinant of the Jacobian, and returns it along with its partial derivatives. """ - jacobian, jacobian_derivatives = self.jacobian_S1(q) + jacobian, hessian = self.jacobian_S1(q) jacobian_inv = np.linalg.pinv(jacobian) - derivatives = np.array([np.trace(jacobian_inv @ djac) for djac in jacobian_derivatives]) + derivatives = np.array([np.trace(jacobian_inv @ djac) for djac in hessian]) return self.log_jacobian_det(q), derivatives @abstractmethod @@ -90,14 +84,20 @@ def n_parameters(self) -> int: """Returns the dimension of the parameter space this transformation is defined over.""" pass - @abstractmethod def to_model(self, q: np.ndarray) -> np.ndarray: """Transforms a parameter vector `q` from the search space to the model space.""" - pass + return self._transform(q, method="to_model") - @abstractmethod def to_search(self, p: np.ndarray) -> np.ndarray: """Transforms a parameter vector `p` from the model space to the search space.""" + return self._transform(p, method="to_search") + + @abstractmethod + def _transform(self, x: np.ndarray, method: str) -> np.ndarray: + """ + Transforms a parameter vector `x` from the search space to the model space if `method` + is "to_model", or from the model space to the search space if `method` is "to_search". + """ pass def is_elementwise(self) -> bool: @@ -107,51 +107,39 @@ def is_elementwise(self) -> bool: """ raise NotImplementedError("is_elementwise method must be implemented if used.") -class TransformedLogPDF(BaseCost): - """Transformed log-PDF class.""" - def __init__(self, log_pdf, transformation): - self._log_pdf = log_pdf - self._transformation = transformation - - def __call__(self, q): - p = self._transformation.to_model(q) - log_pdf = self._log_pdf(p) - - # Calculate the PDF using change of variable - # Wikipedia: https://w.wiki/UsJ - log_jacobian_det = self._transformation.log_jacobian_det(q) - return log_pdf + log_jacobian_det - - def _evaluateS1(self, x): - p = self._transformation.to_model(x) - log_pdf, log_pdf_derivatives = self._log_pdf._evaluateS1(p) - log_jacobian_det, log_jacobian_det_derivatives = self._transformation.log_jacobian_det_S1(x) - return log_pdf + log_jacobian_det, log_pdf_derivatives + log_jacobian_det_derivatives - -class TransformedLogPrior: - """Transformed log-prior class.""" - def __init__(self, log_prior, transformation): - self._log_prior = log_prior - self._transformation = transformation - - def __call__(self, q): - return self - -class TransformedCost(BaseCost): - """Transformed cost class.""" - def __init__(self, cost, transformation): - self._cost = cost - self._transformation = transformation - - def __call__(self, q ): - p = self._transformation.to_model(q) - return self._cost(p) - - def _evaluateS1(self, x): - p = self._transformation.to_model(x) - return self._cost._evaluateS1(x) - -class TransformedBoundaries: - """Transformed boundaries class.""" +class TransformedBoundaries(PintsTransformedBoundaries): + """Transformed boundaries class inherited from Pints TransformedBoundaries.""" def __init__(self, boundaries, transformation): - self._boundaries = boundaries + super().__init__(boundaries, transformation) + + +# ---- To be implemented with Monte Carlo PR ------ # +# class TransformedLogPDF(BaseCost): +# """Transformed log-PDF class.""" +# def __init__(self, log_pdf, transformation): +# self._log_pdf = log_pdf +# self._transformation = transformation + +# def __call__(self, q): +# p = self._transformation.to_model(q) +# log_pdf = self._log_pdf(p) + +# # Calculate the PDF using change of variable +# # Wikipedia: https://w.wiki/UsJ +# log_jacobian_det = self._transformation.log_jacobian_det(q) +# return log_pdf + log_jacobian_det + +# def _evaluateS1(self, x): +# p = self._transformation.to_model(x) +# log_pdf, log_pdf_derivatives = self._log_pdf._evaluateS1(p) +# log_jacobian_det, log_jacobian_det_derivatives = self._transformation.log_jacobian_det_S1(x) +# return log_pdf + log_jacobian_det, log_pdf_derivatives + log_jacobian_det_derivatives + +# class TransformedLogPrior: +# """Transformed log-prior class.""" +# def __init__(self, log_prior, transformation): +# self._log_prior = log_prior +# self._transformation = transformation + +# def __call__(self, q): +# return self diff --git a/pybop/transformations/_transformations.py b/pybop/transformations/_transformations.py index a7a1541c..32f38327 100644 --- a/pybop/transformations/_transformations.py +++ b/pybop/transformations/_transformations.py @@ -1,4 +1,4 @@ -from typing import Sequence, Tuple, Union +from typing import List, Tuple, Union import numpy as np @@ -12,11 +12,11 @@ class IdentityTransformation(Transformation): Based on pints.IdentityTransformation method. """ - def __init__(self, n_parameters: int): + def __init__(self, n_parameters: int = 1): self._n_parameters = n_parameters - def elementwise(self) -> bool: - """See :meth:`Transformation.elementwise()`.""" + def is_elementwise(self) -> bool: + """See :meth:`Transformation.is_elementwise()`.""" return True def jacobian(self, q: np.ndarray) -> np.ndarray: @@ -40,17 +40,256 @@ def n_parameters(self) -> int: """See :meth:`Transformation.n_parameters()`.""" return self._n_parameters - def to_model(self, q: Union[np.ndarray, Sequence, float, int]) -> np.ndarray: - """See :meth:`Transformation.to_model()`.""" - return np.asarray(q) + def _transform(self, x: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" + return np.asarray(x) - def to_search(self, p: Union[np.ndarray, Sequence, float, int]) -> np.ndarray: - """See :meth:`Transformation.to_search()`.""" - return np.asarray(p) + +class ScaledTransformation(Transformation): + """ + Scaled transformation between two parameter spaces: the model parameter space and a search space. + + Based on pints.ScaledTransformation method. + """ + + def __init__( + self, + scale: Union[list, float, np.ndarray], + translate: Union[list, float, np.ndarray] = 0, + n_parameters: int = 1, + ): + self._translate = translate + self._scale = np.asarray([scale]) + self._inverse_scale = np.reciprocal(self._scale) + 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.diag(self._inverse_scale) + + 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 np.sum(np.log(np.abs(self._scale))) + + 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 n_parameters(self) -> int: + """See :meth:`Transformation.n_parameters()`.""" + return self._n_parameters + + def _transform(self, x: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" + x = np.asarray(x) + if method == "to_model": + return x * self._inverse_scale - self._translate + elif method == "to_search": + return self._scale * (x + self._translate) + else: + raise ValueError(f"Unknown method: {method}") + + +class LogTransformation(Transformation): + """ + Log transformation between two parameter spaces: the model parameter space and a search space. + + Based on pints.LogTransformation 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.diag(1 / q) + + def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """See :meth:`Transformation.jacobian_S1()`.""" + return np.diag(1 / q), np.diag(-1 / q**2) + + def log_jacobian_det(self, q: np.ndarray) -> float: + """See :meth:`Transformation.log_jacobian_det()`.""" + return np.sum(-np.log(q)) + + 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), -1 / q + + def n_parameters(self) -> int: + """See :meth:`Transformation.n_parameters()`.""" + return self._n_parameters + + def _transform(self, x: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" + if isinstance(x, (int, float)): + x = np.asarray([x]) + if method == "to_model": + return np.exp(x) + elif method == "to_search": + return np.log(x) + else: + raise ValueError(f"Unknown method: {method}") + + +class ComposedTransformation(Transformation): + """ + N-dimensional Transformation composed of one or more other N_i-dimensional + sub-transformations, where the sum of N_i equals N. + + The dimensionality of the individual transformations does not have to be + the same, i.e., N_i != N_j is allowed. + + Extends pybop.Transformation. Based on pints.ComposedTransformation method. + """ + + def __init__(self, transformations: Transformation): + if not [*transformations]: + raise ValueError("Must have at least one sub-transformation.") + self._transformations: List[Transformation] = transformations + self._n_parameters: int = len(transformations) + self._is_elementwise: bool = all(t.is_elementwise() for t in transformations) + + if self._is_elementwise: + self._jacobian = self._elementwise_jacobian + self._log_jacobian_det = self._elementwise_log_jacobian_det + self._log_jacobian_det_S1 = self._elementwise_log_jacobian_det_S1 + else: + self._jacobian = self._general_jacobian + self._log_jacobian_det = self._general_log_jacobian_det + self._log_jacobian_det_S1 = self._general_log_jacobian_det_S1 + + def is_elementwise(self) -> bool: + return self._is_elementwise + + def jacobian(self, q: np.ndarray) -> np.ndarray: + return self._jacobian(q) + + def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + q = np.asarray(q) + output_S1 = np.zeros( + (self._n_parameters, self._n_parameters, self._n_parameters) + ) + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters() + _, jac_S1 = transformation.jacobian_S1(q[lo:hi]) + for i, jac_S1_i in enumerate(jac_S1): + output_S1[lo + i, lo:hi, lo:hi] = jac_S1_i + lo = hi + + return self.jacobian(q), output_S1 + + def log_jacobian_det(self, q: np.ndarray) -> float: + return self._log_jacobian_det(q) + + def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + return self._log_jacobian_det_S1(q) + + def _transform(self, data: np.ndarray, method: str) -> np.ndarray: + data = np.asarray(data) + output = np.zeros_like(data) + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters() + output[lo:hi] = getattr(transformation, method)(data[lo:hi]) + lo = hi + + return output + + def n_parameters(self) -> int: + return self._n_parameters + + def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: + q = np.asarray(q) + diag = np.zeros(self._n_parameters) + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters() + diag[lo:hi] = np.diagonal(transformation.jacobian(q[lo:hi])) + lo = hi + + return np.diag(diag) + + def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: + q = np.asarray(q) + return sum( + transformation.log_jacobian_det(q[lo : lo + transformation.n_parameters()]) + for lo, transformation in self._iter_transformations() + ) + + def _elementwise_log_jacobian_det_S1( + self, q: np.ndarray + ) -> Tuple[float, np.ndarray]: + q = np.asarray(q) + output = 0.0 + output_S1 = np.zeros(self._n_parameters) + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters() + j, j_S1 = transformation.log_jacobian_det_S1(q[lo:hi]) + output += j + output_S1[lo:hi] = j_S1 + lo = hi + + return output, output_S1 + + def _general_jacobian(self, q: np.ndarray) -> np.ndarray: + q = np.asarray(q) + jacobian_blocks = [] + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters() + jacobian_blocks.append(transformation.jacobian(q[lo:hi])) + lo = hi + + return np.block( + [ + [ + b if i == j else np.zeros_like(b) + for j, b in enumerate(jacobian_blocks) + ] + for i, _ in enumerate(jacobian_blocks) + ] + ) + + def _general_log_jacobian_det(self, q: np.ndarray) -> float: + return np.log(np.abs(np.linalg.det(self.jacobian(q)))) + + def _general_log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + q = np.asarray(q) + jac, jac_S1 = self.jacobian_S1(q) + out_S1 = np.zeros(self._n_parameters) + + for i, jac_S1_i in enumerate(jac_S1): + out_S1[i] = np.trace(np.matmul(np.linalg.pinv(jac), jac_S1_i)) + + return self.log_jacobian_det(q), out_S1 + + def _iter_transformations(self): + lo = 0 + for transformation in self._transformations: + yield lo, transformation + lo += transformation.n_parameters() # TODO: Implement the following classes: -# class LogTransformation(Transformation): # class LogitTransformation(Transformation): -# class ComposedTransformation(Transformation): -# class ScaledTransformation(Transformation): From 7a12d6b866ea8934fddef903aa82b95b2cb0670d Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 21 Jun 2024 17:28:54 +0100 Subject: [PATCH 03/18] tests: add tests, remove redundant methods, add catch on x0==ground_truth --- pybop/transformations/__init__.py | 20 ++-- pybop/transformations/_transformations.py | 2 +- .../integration/test_optimisation_options.py | 2 + .../integration/test_spm_parameterisations.py | 2 + .../test_thevenin_parameterisation.py | 2 + tests/integration/test_transformation.py | 102 ++++++++++++++++++ tests/unit/test_transformations.py | 71 ++++++++++++ 7 files changed, 187 insertions(+), 14 deletions(-) create mode 100644 tests/integration/test_transformation.py create mode 100644 tests/unit/test_transformations.py diff --git a/pybop/transformations/__init__.py b/pybop/transformations/__init__.py index f1ad8cfc..ab87fd93 100644 --- a/pybop/transformations/__init__.py +++ b/pybop/transformations/__init__.py @@ -9,9 +9,9 @@ class Transformation(ABC): Abstract base class for transformations between two parameter spaces: the model parameter space and a search space. - If `trans` is an instance of a `Transformation` class, you can apply the + If `transform` is an instance of a `Transformation` class, you can apply the transformation of a parameter vector from the model space `p` to the search - space `q` using `q = trans.to_search(p)` and the inverse using `p = trans.to_model(q)`. + space `q` using `q = transform.to_search(p)` and the inverse using `p = transform.to_model(q)`. Based on pints.transformation method. @@ -21,14 +21,14 @@ class Transformation(ABC): http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.47.9023 .. [2] Kaare Brandt Petersen and Michael Syskind Pedersen. "The Matrix Cookbook." 2012. """ - - def convert_log_prior(self, log_prior): - """Returns a transformed log-prior class.""" - return TransformedLogPrior(log_prior, self) + # ---- To be implemented with Monte Carlo PR ------ # + # def convert_log_prior(self, log_prior): + # """Returns a transformed log-prior class.""" + # return TransformedLogPrior(log_prior, self) def convert_boundaries(self, boundaries): """Returns a transformed boundaries class.""" - return TransformedBoundaries(boundaries, self) + return PintsTransformedBoundaries(boundaries, self) def convert_covariance_matrix(self, covariance: np.ndarray, q: np.ndarray) -> np.ndarray: """ @@ -107,12 +107,6 @@ def is_elementwise(self) -> bool: """ raise NotImplementedError("is_elementwise method must be implemented if used.") -class TransformedBoundaries(PintsTransformedBoundaries): - """Transformed boundaries class inherited from Pints TransformedBoundaries.""" - def __init__(self, boundaries, transformation): - super().__init__(boundaries, transformation) - - # ---- To be implemented with Monte Carlo PR ------ # # class TransformedLogPDF(BaseCost): # """Transformed log-PDF class.""" diff --git a/pybop/transformations/_transformations.py b/pybop/transformations/_transformations.py index 32f38327..8b9dd8b3 100644 --- a/pybop/transformations/_transformations.py +++ b/pybop/transformations/_transformations.py @@ -59,7 +59,7 @@ def __init__( n_parameters: int = 1, ): self._translate = translate - self._scale = np.asarray([scale]) + self._scale = np.asarray([float(scale)]) self._inverse_scale = np.reciprocal(self._scale) self._n_parameters = n_parameters diff --git a/tests/integration/test_optimisation_options.py b/tests/integration/test_optimisation_options.py index dcd94276..3103ef35 100644 --- a/tests/integration/test_optimisation_options.py +++ b/tests/integration/test_optimisation_options.py @@ -104,6 +104,8 @@ def test_optimisation_f_guessed(self, f_guessed, spm_costs): assert initial_cost > final_cost else: assert initial_cost < final_cost + else: + raise ValueError("Initial value is the same as the ground truth value.") np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) def get_data(self, model, parameters, x, init_soc): diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 9ae2b421..60ef4381 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -122,6 +122,8 @@ def test_spm_optimisers(self, optimiser, spm_costs): assert initial_cost > final_cost else: assert initial_cost < final_cost + else: + raise ValueError("Initial value is the same as the ground truth value.") np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) @pytest.fixture diff --git a/tests/integration/test_thevenin_parameterisation.py b/tests/integration/test_thevenin_parameterisation.py index 57bb0689..37d67b62 100644 --- a/tests/integration/test_thevenin_parameterisation.py +++ b/tests/integration/test_thevenin_parameterisation.py @@ -90,6 +90,8 @@ def test_optimisers_on_simple_model(self, optimiser, cost): assert initial_cost > final_cost else: assert initial_cost < final_cost + else: + raise ValueError("Initial value is the same as the ground truth value.") np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) def get_data(self, model, parameters, x): diff --git a/tests/integration/test_transformation.py b/tests/integration/test_transformation.py new file mode 100644 index 00000000..45feff80 --- /dev/null +++ b/tests/integration/test_transformation.py @@ -0,0 +1,102 @@ +import numpy as np +import pytest + +import pybop + + +class TestTransformation: + """ + A class for integration tests of the transformation methods. + """ + + @pytest.fixture(autouse=True) + def setup(self): + self.ground_truth = np.array([0.5, 0.2]) + np.random.normal( + loc=0.0, scale=0.05, size=2 + ) + + @pytest.fixture + def model(self): + parameter_set = pybop.ParameterSet.pybamm("Chen2020") + return pybop.lithium_ion.SPMe(parameter_set=parameter_set) + + @pytest.fixture + def parameters(self): + return pybop.Parameters( + pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Uniform(0.4, 0.7), + bounds=[0.375, 0.725], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode conductivity [S.m-1]", + prior=pybop.Uniform(0.05, 0.45), + bounds=[0.04, 0.5], + transformation=pybop.LogTransformation(), + ), + ) + + @pytest.fixture(params=[0.4, 0.7]) + def init_soc(self, request): + return request.param + + def noise(self, sigma, values): + return np.random.normal(0, sigma, values) + + @pytest.fixture + def cost(self, model, parameters, init_soc): + # Form dataset + solution = self.get_data(model, parameters, self.ground_truth, init_soc) + dataset = pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": solution["Voltage [V]"].data + + self.noise(0.002, len(solution["Time [s]"].data)), + } + ) + + # Define the cost to optimise + problem = pybop.FittingProblem(model, parameters, dataset, init_soc=init_soc) + return pybop.RootMeanSquaredError(problem) + + @pytest.mark.parametrize( + "optimiser", + [ + pybop.AdamW, + pybop.CMAES, + ], + ) + @pytest.mark.integration + def test_spm_optimisers(self, optimiser, cost): + x0 = cost.parameters.initial_value() + optim = optimiser( + cost=cost, + max_iterations=125, + ) + + initial_cost = optim.cost(x0) + x, final_cost = optim.run() + + # Assertions + if not np.allclose(x0, self.ground_truth, atol=1e-5): + assert initial_cost > final_cost + np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) + else: + raise ValueError("Initial value is the same as the ground truth value.") + + def get_data(self, model, parameters, x, init_soc): + model.parameters = parameters + experiment = pybop.Experiment( + [ + ( + "Discharge at 1C for 3 minutes (4 second period)", + "Charge at 1C for 3 minutes (4 second period)", + ), + ] + ) + sim = model.predict( + init_soc=init_soc, experiment=experiment, inputs=parameters.as_dict(x) + ) + return sim diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py new file mode 100644 index 00000000..60dd318e --- /dev/null +++ b/tests/unit/test_transformations.py @@ -0,0 +1,71 @@ +import numpy as np +import pytest + +import pybop + + +class TestTransformations: + """ + A class to test the transformations. + """ + + @pytest.fixture + def parameters(self): + return pybop.Parameters( + pybop.Parameter( + "Identity", + transformation=pybop.IdentityTransformation(), + ), + pybop.Parameter( + "Scaled", + transformation=pybop.ScaledTransformation(scale=2.0, translate=1), + ), + pybop.Parameter( + "Log", + transformation=pybop.LogTransformation(), + ), + ) + + @pytest.mark.unit + def test_identity_transformation(self, parameters): + q = np.array([5.0]) + transformation = parameters["Identity"].transformation + assert np.array_equal(transformation.to_model(q), q) + assert np.array_equal(transformation.to_search(q), q) + assert transformation.log_jacobian_det(q) == 0.0 + + jac, jac_S1 = transformation.jacobian_S1(q) + assert np.array_equal(jac, np.eye(1)) + assert np.array_equal(jac_S1, np.zeros((1, 1, 1))) + + @pytest.mark.unit + def test_scaled_transformation(self, parameters): + q = np.array([2.5]) + transformation = parameters["Scaled"].transformation + p = transformation.to_model(q) + assert np.allclose(p, (q / 2.0) - 1.0) + + q_transformed = transformation.to_search(p) + assert np.allclose(q_transformed, q) + assert np.allclose( + transformation.log_jacobian_det(q), np.sum(np.log(np.abs(2.0))) + ) + + jac, jac_S1 = transformation.jacobian_S1(q) + assert np.array_equal(jac, np.diag([0.5])) + assert np.array_equal(jac_S1, np.zeros((1, 1, 1))) + + @pytest.mark.unit + def test_log_transformation(self, parameters): + q = np.array([10]) + transformation = parameters["Log"].transformation + p = transformation.to_model(q) + assert np.allclose(p, np.exp(q)) + + q_transformed = transformation.to_search(p) + assert np.allclose(q_transformed, q) + assert np.allclose(transformation.log_jacobian_det(q), -np.sum(np.log(q))) + + jac, jac_S1 = transformation.jacobian_S1(q) + assert np.array_equal(jac, np.diag(1 / q)) + assert np.array_equal(jac_S1, np.diag(-1 / q**2)) From 9a4b34bf6d3e6666b71efa9ed31e9106b3650e23 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sat, 22 Jun 2024 11:37:35 +0100 Subject: [PATCH 04/18] fix: add _verify_inputs, increase coverage, add temp fix to GaussianLogLikelihood --- pybop/__init__.py | 4 +- pybop/costs/_likelihoods.py | 4 + .../__init__.py | 38 ++++---- .../_transformation.py} | 86 ++++++++++--------- tests/unit/test_likelihoods.py | 6 +- tests/unit/test_transformations.py | 64 +++++++++++++- 6 files changed, 139 insertions(+), 63 deletions(-) rename pybop/{transformations => transformation}/__init__.py (82%) rename pybop/{transformations/_transformations.py => transformation/_transformation.py} (82%) diff --git a/pybop/__init__.py b/pybop/__init__.py index 82d527b0..5755a4b3 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -81,8 +81,8 @@ # # Transformation classes # -from .transformations import Transformation -from .transformations._transformations import ( +from .transformation import Transformation +from .transformation._transformation import ( IdentityTransformation, ScaledTransformation, LogTransformation, diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 67ad97ec..615184bc 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -1,5 +1,6 @@ import numpy as np +from pybop import IdentityTransformation from pybop.costs.base_cost import BaseCost @@ -124,6 +125,9 @@ def __init__(self, problem): super(GaussianLogLikelihood, self).__init__(problem) self._logpi = -0.5 * self.n_time_data * np.log(2 * np.pi) self._dl = np.ones(self.n_parameters + self.n_outputs) + self.transformation.append( + IdentityTransformation() + ) # Temporary fix, ahead of #338 def _evaluate(self, x): """ diff --git a/pybop/transformations/__init__.py b/pybop/transformation/__init__.py similarity index 82% rename from pybop/transformations/__init__.py rename to pybop/transformation/__init__.py index ab87fd93..20c7819f 100644 --- a/pybop/transformations/__init__.py +++ b/pybop/transformation/__init__.py @@ -1,7 +1,6 @@ from abc import ABC, abstractmethod -from typing import Tuple, Union, Sequence +from typing import Tuple, Union, Sequence, List import numpy as np -from pints import TransformedBoundaries as PintsTransformedBoundaries class Transformation(ABC): @@ -26,17 +25,13 @@ class Transformation(ABC): # """Returns a transformed log-prior class.""" # return TransformedLogPrior(log_prior, self) - def convert_boundaries(self, boundaries): - """Returns a transformed boundaries class.""" - return PintsTransformedBoundaries(boundaries, self) - - def convert_covariance_matrix(self, covariance: np.ndarray, q: np.ndarray) -> np.ndarray: + def convert_covariance_matrix(self, cov: np.ndarray, q: np.ndarray) -> np.ndarray: """ Converts a covariance matrix `covariance` from the model space to the search space around a parameter vector `q` in the search space. """ jac_inv = np.linalg.pinv(self.jacobian(q)) - return jac_inv @ covariance @ jac_inv.T + return jac_inv @ cov @ jac_inv.T def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarray) -> np.ndarray: """ @@ -67,22 +62,18 @@ def log_jacobian_det(self, q: np.ndarray) -> float: Returns the logarithm of the absolute value of the determinant of the Jacobian matrix at the parameter vector `q`. """ - return np.log(np.abs(np.linalg.det(self.jacobian(q)))) + raise NotImplementedError("log_jacobian_det method must be implemented if used.") def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: """ Computes the logarithm of the absolute value of the determinant of the Jacobian, and returns it along with its partial derivatives. """ - jacobian, hessian = self.jacobian_S1(q) - jacobian_inv = np.linalg.pinv(jacobian) - derivatives = np.array([np.trace(jacobian_inv @ djac) for djac in hessian]) - return self.log_jacobian_det(q), derivatives + raise NotImplementedError("log_jacobian_det_S1 method must be implemented if used.") - @abstractmethod - def n_parameters(self) -> int: - """Returns the dimension of the parameter space this transformation is defined over.""" - pass + @property + def n_parameters(self): + return self._n_parameters def to_model(self, q: np.ndarray) -> np.ndarray: """Transforms a parameter vector `q` from the search space to the model space.""" @@ -107,6 +98,19 @@ def is_elementwise(self) -> bool: """ raise NotImplementedError("is_elementwise method must be implemented if used.") + def _verify_input(self, input: Union[List[float], float, np.ndarray]) -> None: + """Set and validate the translate parameter.""" + if isinstance(input, (float, int)): + input = np.full(self._n_parameters, float(input)) + elif isinstance(input, (list, np.ndarray)): + if len(input) != self._n_parameters: + raise ValueError(f"Translate must have {self._n_parameters} elements") + input = np.asarray(input, dtype=float) + else: + raise TypeError("Translate must be a float, list, or numpy array") + + return input + # ---- To be implemented with Monte Carlo PR ------ # # class TransformedLogPDF(BaseCost): # """Transformed log-PDF class.""" diff --git a/pybop/transformations/_transformations.py b/pybop/transformation/_transformation.py similarity index 82% rename from pybop/transformations/_transformations.py rename to pybop/transformation/_transformation.py index 8b9dd8b3..37db5a3e 100644 --- a/pybop/transformations/_transformations.py +++ b/pybop/transformation/_transformation.py @@ -36,10 +36,6 @@ 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 n_parameters(self) -> int: - """See :meth:`Transformation.n_parameters()`.""" - return self._n_parameters - def _transform(self, x: np.ndarray, method: str) -> np.ndarray: """See :meth:`Transformation._transform`.""" return np.asarray(x) @@ -58,10 +54,10 @@ def __init__( translate: Union[list, float, np.ndarray] = 0, n_parameters: int = 1, ): - self._translate = translate - self._scale = np.asarray([float(scale)]) - self._inverse_scale = np.reciprocal(self._scale) self._n_parameters = n_parameters + self._translate = self._verify_input(translate) + self._scale = self._verify_input(scale) + self._inverse_scale = np.reciprocal(self._scale) def is_elementwise(self) -> bool: """See :meth:`Transformation.is_elementwise()`.""" @@ -84,13 +80,9 @@ 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 n_parameters(self) -> int: - """See :meth:`Transformation.n_parameters()`.""" - return self._n_parameters - def _transform(self, x: np.ndarray, method: str) -> np.ndarray: """See :meth:`Transformation._transform`.""" - x = np.asarray(x) + x = self._verify_input(x) if method == "to_model": return x * self._inverse_scale - self._translate elif method == "to_search": @@ -129,14 +121,9 @@ 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), -1 / q - def n_parameters(self) -> int: - """See :meth:`Transformation.n_parameters()`.""" - return self._n_parameters - def _transform(self, x: np.ndarray, method: str) -> np.ndarray: """See :meth:`Transformation._transform`.""" - if isinstance(x, (int, float)): - x = np.asarray([x]) + x = self._verify_input(x) if method == "to_model": return np.exp(x) elif method == "to_search": @@ -156,13 +143,24 @@ class ComposedTransformation(Transformation): Extends pybop.Transformation. Based on pints.ComposedTransformation method. """ - def __init__(self, transformations: Transformation): - if not [*transformations]: + def __init__(self, transformations: List[Transformation]): + if not transformations: raise ValueError("Must have at least one sub-transformation.") - self._transformations: List[Transformation] = transformations - self._n_parameters: int = len(transformations) - self._is_elementwise: bool = all(t.is_elementwise() for t in transformations) - + self._transformations = [] + self._n_parameters = 0 + self._is_elementwise = True + for transformation in transformations: + self._append_transformation(transformation) + self._update_methods() + + def _append_transformation(self, transformation: Transformation): + if not isinstance(transformation, Transformation): + raise ValueError("The appended object must be a Transformation.") + self._transformations.append(transformation) + self._n_parameters += transformation.n_parameters + self._is_elementwise = self._is_elementwise and transformation.is_elementwise() + + def _update_methods(self): if self._is_elementwise: self._jacobian = self._elementwise_jacobian self._log_jacobian_det = self._elementwise_log_jacobian_det @@ -172,6 +170,16 @@ def __init__(self, transformations: Transformation): self._log_jacobian_det = self._general_log_jacobian_det self._log_jacobian_det_S1 = self._general_log_jacobian_det_S1 + def append(self, transformation: Transformation): + """ + Append a new transformation to the existing composition. + + Args: + transformation (Transformation): The transformation to append. + """ + self._append_transformation(transformation) + self._update_methods() + def is_elementwise(self) -> bool: return self._is_elementwise @@ -179,14 +187,14 @@ def jacobian(self, q: np.ndarray) -> np.ndarray: return self._jacobian(q) def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: - q = np.asarray(q) + q = self._verify_input(q) output_S1 = np.zeros( (self._n_parameters, self._n_parameters, self._n_parameters) ) lo = 0 for transformation in self._transformations: - hi = lo + transformation.n_parameters() + hi = lo + transformation.n_parameters _, jac_S1 = transformation.jacobian_S1(q[lo:hi]) for i, jac_S1_i in enumerate(jac_S1): output_S1[lo + i, lo:hi, lo:hi] = jac_S1_i @@ -195,40 +203,38 @@ def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: return self.jacobian(q), output_S1 def log_jacobian_det(self, q: np.ndarray) -> float: + """See :meth:`Transformation.log_jacobian_det()`.""" return self._log_jacobian_det(q) def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: return self._log_jacobian_det_S1(q) def _transform(self, data: np.ndarray, method: str) -> np.ndarray: - data = np.asarray(data) + data = self._verify_input(data) output = np.zeros_like(data) lo = 0 for transformation in self._transformations: - hi = lo + transformation.n_parameters() + hi = lo + transformation.n_parameters output[lo:hi] = getattr(transformation, method)(data[lo:hi]) lo = hi return output - def n_parameters(self) -> int: - return self._n_parameters - def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: - q = np.asarray(q) + q = self._verify_input(q) diag = np.zeros(self._n_parameters) lo = 0 for transformation in self._transformations: - hi = lo + transformation.n_parameters() + hi = lo + transformation.n_parameters diag[lo:hi] = np.diagonal(transformation.jacobian(q[lo:hi])) lo = hi return np.diag(diag) def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: - q = np.asarray(q) + q = self._verify_input(q) return sum( transformation.log_jacobian_det(q[lo : lo + transformation.n_parameters()]) for lo, transformation in self._iter_transformations() @@ -237,13 +243,13 @@ def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: def _elementwise_log_jacobian_det_S1( self, q: np.ndarray ) -> Tuple[float, np.ndarray]: - q = np.asarray(q) + q = self._verify_input(q) output = 0.0 output_S1 = np.zeros(self._n_parameters) lo = 0 for transformation in self._transformations: - hi = lo + transformation.n_parameters() + hi = lo + transformation.n_parameters j, j_S1 = transformation.log_jacobian_det_S1(q[lo:hi]) output += j output_S1[lo:hi] = j_S1 @@ -252,12 +258,12 @@ def _elementwise_log_jacobian_det_S1( return output, output_S1 def _general_jacobian(self, q: np.ndarray) -> np.ndarray: - q = np.asarray(q) + q = self._verify_input(q) jacobian_blocks = [] lo = 0 for transformation in self._transformations: - hi = lo + transformation.n_parameters() + hi = lo + transformation.n_parameters jacobian_blocks.append(transformation.jacobian(q[lo:hi])) lo = hi @@ -275,7 +281,7 @@ def _general_log_jacobian_det(self, q: np.ndarray) -> float: return np.log(np.abs(np.linalg.det(self.jacobian(q)))) def _general_log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: - q = np.asarray(q) + q = self._verify_input(q) jac, jac_S1 = self.jacobian_S1(q) out_S1 = np.zeros(self._n_parameters) @@ -288,7 +294,7 @@ def _iter_transformations(self): lo = 0 for transformation in self._transformations: yield lo, transformation - lo += transformation.n_parameters() + lo += transformation.n_parameters # TODO: Implement the following classes: diff --git a/tests/unit/test_likelihoods.py b/tests/unit/test_likelihoods.py index 41ee3667..cc55ca73 100644 --- a/tests/unit/test_likelihoods.py +++ b/tests/unit/test_likelihoods.py @@ -86,7 +86,7 @@ def test_base_likelihood_call_raises_not_implemented_error( ): likelihood = pybop.BaseLikelihood(one_signal_problem) with pytest.raises(NotImplementedError): - likelihood(np.array([0.5, 0.5])) + likelihood(np.array([0.5])) @pytest.mark.unit def test_set_get_sigma(self, one_signal_problem): @@ -128,8 +128,8 @@ def test_gaussian_log_likelihood_known_sigma(self, problem_name, request): @pytest.mark.unit def test_gaussian_log_likelihood(self, one_signal_problem): likelihood = pybop.GaussianLogLikelihood(one_signal_problem) - result = likelihood(np.array([0.5, 0.5])) - grad_result, grad_likelihood = likelihood.evaluateS1(np.array([0.5, 0.5])) + result = likelihood(np.array([0.8, 0.2])) + grad_result, grad_likelihood = likelihood.evaluateS1(np.array([0.8, 0.2])) assert isinstance(result, float) np.testing.assert_allclose(result, grad_result, atol=1e-5) assert np.all(grad_likelihood <= 0) diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index 60dd318e..5b917cea 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -1,10 +1,12 @@ +import inspect + import numpy as np import pytest import pybop -class TestTransformations: +class TestTransformation: """ A class to test the transformations. """ @@ -33,17 +35,25 @@ def test_identity_transformation(self, parameters): assert np.array_equal(transformation.to_model(q), q) assert np.array_equal(transformation.to_search(q), q) assert transformation.log_jacobian_det(q) == 0.0 + assert transformation.n_parameters == 1 jac, jac_S1 = transformation.jacobian_S1(q) assert np.array_equal(jac, np.eye(1)) assert np.array_equal(jac_S1, np.zeros((1, 1, 1))) + # Test covariance transformation + cov = np.array([[0.5]]) + q = np.array([5.0]) + cov_transformed = transformation.convert_covariance_matrix(cov, q) + assert np.array_equal(cov_transformed, cov) + @pytest.mark.unit def test_scaled_transformation(self, parameters): q = np.array([2.5]) transformation = parameters["Scaled"].transformation p = transformation.to_model(q) assert np.allclose(p, (q / 2.0) - 1.0) + assert transformation.n_parameters == 1 q_transformed = transformation.to_search(p) assert np.allclose(q_transformed, q) @@ -55,12 +65,18 @@ def test_scaled_transformation(self, parameters): assert np.array_equal(jac, np.diag([0.5])) assert np.array_equal(jac_S1, np.zeros((1, 1, 1))) + # Test covariance transformation + cov = np.array([[0.5]]) + cov_transformed = transformation.convert_covariance_matrix(cov, q) + assert np.array_equal(cov_transformed, cov * transformation._scale**2) + @pytest.mark.unit def test_log_transformation(self, parameters): q = np.array([10]) transformation = parameters["Log"].transformation p = transformation.to_model(q) assert np.allclose(p, np.exp(q)) + assert transformation.n_parameters == 1 q_transformed = transformation.to_search(p) assert np.allclose(q_transformed, q) @@ -69,3 +85,49 @@ def test_log_transformation(self, parameters): jac, jac_S1 = transformation.jacobian_S1(q) assert np.array_equal(jac, np.diag(1 / q)) assert np.array_equal(jac_S1, np.diag(-1 / q**2)) + + # Test covariance transformation + cov = np.array([[0.5]]) + cov_transformed = transformation.convert_covariance_matrix(cov, q) + assert np.array_equal(cov_transformed, cov * q**2) + + +class TestBaseTransformation: + """ + A class to test the abstract base transformation class. + """ + + @pytest.mark.unit + def test_abstract_base_transformation(self): + with pytest.raises(TypeError): + pybop.Transformation() + + @pytest.mark.unit + def test_abstract_methods(self): + abstract_methods = ["jacobian", "_transform"] + for method in abstract_methods: + assert hasattr(pybop.Transformation, method) + assert getattr(pybop.Transformation, method).__isabstractmethod__ + + @pytest.mark.unit + def test_concrete_methods(self): + concrete_methods = [ + "convert_covariance_matrix", + "convert_standard_deviation", + "log_jacobian_det", + "to_model", + "to_search", + ] + for method in concrete_methods: + assert hasattr(pybop.Transformation, method) + assert not inspect.isabstract(getattr(pybop.Transformation, method)) + + @pytest.mark.unit + def test_not_implemented_methods(self): + not_implemented_methods = [ + "jacobian_S1", + "log_jacobian_det_S1", + "is_elementwise", + ] + for method in not_implemented_methods: + assert hasattr(pybop.Transformation, method) From aab8995fc7bd6911af7ead9946f34f48ae3ba5ad Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sat, 22 Jun 2024 19:16:25 +0100 Subject: [PATCH 05/18] docs: update docstrings --- pybop/transformation/_transformation.py | 208 +++++++++++++++++++++--- 1 file changed, 186 insertions(+), 22 deletions(-) diff --git a/pybop/transformation/_transformation.py b/pybop/transformation/_transformation.py index 37db5a3e..7746ddc7 100644 --- a/pybop/transformation/_transformation.py +++ b/pybop/transformation/_transformation.py @@ -7,9 +7,29 @@ class IdentityTransformation(Transformation): """ - Identity transformation between two parameter spaces: the model parameter space and a search space. + This class implements a trivial transformation where the model parameter space + and the search space are identical. It extends the base Transformation class. - Based on pints.IdentityTransformation method. + 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): @@ -43,21 +63,36 @@ def _transform(self, x: np.ndarray, method: str) -> np.ndarray: class ScaledTransformation(Transformation): """ - Scaled transformation between two parameter spaces: the model parameter space and a search space. + This class implements a linear transformation between the model parameter space + and a search space, using a coefficient (scale factor) and an intercept (offset). + It extends the base Transformation class. + + The transformation is defined as: + - to_search: y = coefficient * (x + intercept) + - to_model: x = y / coefficient - intercept + + Where: + - x is in the model parameter space + - y is in the search space + - coefficient is the scaling factor + - intercept is the offset + + This transformation is useful for scaling and shifting parameters to a more + suitable range for optimisation algorithms. - Based on pints.ScaledTransformation method. + Based on pints.ScaledTransformation class. """ def __init__( self, - scale: Union[list, float, np.ndarray], - translate: Union[list, float, np.ndarray] = 0, + coefficient: Union[list, float, np.ndarray], + intercept: Union[list, float, np.ndarray] = 0, n_parameters: int = 1, ): self._n_parameters = n_parameters - self._translate = self._verify_input(translate) - self._scale = self._verify_input(scale) - self._inverse_scale = np.reciprocal(self._scale) + self._intercept = self._verify_input(intercept) + self._coefficient = self._verify_input(coefficient) + self._inverse_coeff = np.reciprocal(self._coefficient) def is_elementwise(self) -> bool: """See :meth:`Transformation.is_elementwise()`.""" @@ -65,7 +100,7 @@ def is_elementwise(self) -> bool: def jacobian(self, q: np.ndarray) -> np.ndarray: """See :meth:`Transformation.jacobian()`.""" - return np.diag(self._inverse_scale) + return np.diag(self._inverse_coeff) def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """See :meth:`Transformation.jacobian_S1()`.""" @@ -74,7 +109,7 @@ def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: def log_jacobian_det(self, q: np.ndarray) -> float: """See :meth:`Transformation.log_jacobian_det()`.""" - return np.sum(np.log(np.abs(self._scale))) + return np.sum(np.log(np.abs(self._coefficient))) def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: """See :meth:`Transformation.log_jacobian_det_S1()`.""" @@ -84,18 +119,38 @@ def _transform(self, x: np.ndarray, method: str) -> np.ndarray: """See :meth:`Transformation._transform`.""" x = self._verify_input(x) if method == "to_model": - return x * self._inverse_scale - self._translate + return x * self._inverse_coeff - self._intercept elif method == "to_search": - return self._scale * (x + self._translate) + return self._coefficient * (x + self._intercept) else: raise ValueError(f"Unknown method: {method}") class LogTransformation(Transformation): """ - Log transformation between two parameter spaces: the model parameter space and a search space. + This class implements a logarithmic transformation between the model parameter space + and a search space. It extends the base Transformation class. - Based on pints.LogTransformation method. + The transformation is defined as: + - to_search: y = log(x) + - to_model: x = exp(y) + + Where: + - x is in the model parameter space (strictly positive) + - y 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. + 2. Converting multiplicative processes to additive ones in the search space. + 3. Ensuring positivity constraints without explicit bounds in optimisation. + + Note: Care should be taken when using this transformation, as it can introduce + bias in the parameter estimates if not accounted for properly in the likelihood + or cost function. Simply, E[log(x)] <= log(E[x]) as per to Jensen's inequality. + For more information, see Jensen's inequality: + https://en.wikipedia.org/w/index.php?title=Jensen%27s_inequality&oldid=1212437916#Probabilistic_form + + Initially based on pints.LogTransformation class. """ def __init__(self, n_parameters: int = 1): @@ -137,10 +192,14 @@ class ComposedTransformation(Transformation): N-dimensional Transformation composed of one or more other N_i-dimensional sub-transformations, where the sum of N_i equals N. + This class allows for the composition of multiple transformations, each potentially + operating on a different number of dimensions. The total dimensionality of the + composed transformation is the sum of the dimensionalities of its components. + The dimensionality of the individual transformations does not have to be the same, i.e., N_i != N_j is allowed. - Extends pybop.Transformation. Based on pints.ComposedTransformation method. + Extends pybop.Transformation. Initially based on pints.ComposedTransformation class. """ def __init__(self, transformations: List[Transformation]): @@ -154,6 +213,19 @@ def __init__(self, transformations: List[Transformation]): self._update_methods() def _append_transformation(self, transformation: Transformation): + """ + Append a transformation to the internal list of transformations. + + Parameters + ---------- + transformation : Transformation + Transformation to append. + + Raises + ------ + ValueError + If the appended object is not a Transformation. + """ if not isinstance(transformation, Transformation): raise ValueError("The appended object must be a Transformation.") self._transformations.append(transformation) @@ -161,6 +233,9 @@ def _append_transformation(self, transformation: Transformation): self._is_elementwise = self._is_elementwise and transformation.is_elementwise() def _update_methods(self): + """ + Update the internal methods based on whether the transformation is elementwise. + """ if self._is_elementwise: self._jacobian = self._elementwise_jacobian self._log_jacobian_det = self._elementwise_log_jacobian_det @@ -174,19 +249,24 @@ def append(self, transformation: Transformation): """ Append a new transformation to the existing composition. - Args: - transformation (Transformation): The transformation to append. + Parameters + ---------- + transformation : Transformation + The transformation to append. """ self._append_transformation(transformation) self._update_methods() def is_elementwise(self) -> bool: + """See :meth:`Transformation.is_elementwise()`.""" return self._is_elementwise def jacobian(self, q: np.ndarray) -> np.ndarray: + """See :meth:`Transformation.jacobian()`.""" return self._jacobian(q) def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """See :meth:`Transformation.jacobian_S1()`.""" q = self._verify_input(q) output_S1 = np.zeros( (self._n_parameters, self._n_parameters, self._n_parameters) @@ -207,9 +287,11 @@ def log_jacobian_det(self, q: np.ndarray) -> float: return self._log_jacobian_det(q) 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_S1(q) def _transform(self, data: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" data = self._verify_input(data) output = np.zeros_like(data) lo = 0 @@ -222,6 +304,19 @@ def _transform(self, data: np.ndarray, method: str) -> np.ndarray: return output def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: + """ + Compute the elementwise Jacobian of the composed transformation. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + np.ndarray + Diagonal matrix representing the elementwise Jacobian. + """ q = self._verify_input(q) diag = np.zeros(self._n_parameters) lo = 0 @@ -234,6 +329,19 @@ def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: return np.diag(diag) def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: + """ + Compute the elementwise logarithm of the determinant of the Jacobian. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + float + Sum of log determinants of individual transformations. + """ q = self._verify_input(q) return sum( transformation.log_jacobian_det(q[lo : lo + transformation.n_parameters()]) @@ -243,6 +351,19 @@ def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: def _elementwise_log_jacobian_det_S1( self, q: np.ndarray ) -> Tuple[float, np.ndarray]: + """ + Compute the elementwise logarithm of the determinant of the Jacobian and its first-order sensitivities. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + Tuple[float, np.ndarray] + Tuple of sum of log determinants and concatenated first-order sensitivities. + """ q = self._verify_input(q) output = 0.0 output_S1 = np.zeros(self._n_parameters) @@ -258,6 +379,19 @@ def _elementwise_log_jacobian_det_S1( return output, output_S1 def _general_jacobian(self, q: np.ndarray) -> np.ndarray: + """ + Compute the general Jacobian of the composed transformation. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + np.ndarray + Block diagonal matrix of the Jacobian of each transformation. + """ q = self._verify_input(q) jacobian_blocks = [] lo = 0 @@ -278,9 +412,35 @@ def _general_jacobian(self, q: np.ndarray) -> np.ndarray: ) def _general_log_jacobian_det(self, q: np.ndarray) -> float: + """ + Compute the general logarithm of the determinant of the Jacobian. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + float + Logarithm of the absolute value of the determinant of the Jacobian. + """ return np.log(np.abs(np.linalg.det(self.jacobian(q)))) def _general_log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + """ + Compute the general logarithm of the determinant of the Jacobian and its first-order sensitivities. + + Parameters + ---------- + q : np.ndarray + Input array. + + Returns + ------- + Tuple[float, np.ndarray] + Tuple of log determinant and its first-order sensitivities. + """ q = self._verify_input(q) jac, jac_S1 = self.jacobian_S1(q) out_S1 = np.zeros(self._n_parameters) @@ -291,11 +451,15 @@ def _general_log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray return self.log_jacobian_det(q), out_S1 def _iter_transformations(self): + """ + Iterate over the transformations in the composition. + + Yields + ------ + Tuple[int, Transformation] + Tuple of starting index and transformation object for each sub-transformation. + """ lo = 0 for transformation in self._transformations: yield lo, transformation lo += transformation.n_parameters - - -# TODO: Implement the following classes: -# class LogitTransformation(Transformation): From feae89fd58505cd6d3f0abf26237ca9fccc1518a Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sun, 23 Jun 2024 17:52:20 +0100 Subject: [PATCH 06/18] fix: add catch for transformation == None, updt tests for arg rename --- pybop/costs/_likelihoods.py | 7 ++++--- pybop/costs/base_cost.py | 25 ++++++++++++++++--------- tests/unit/test_transformations.py | 4 ++-- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 615184bc..7aa027e2 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -125,9 +125,10 @@ def __init__(self, problem): super(GaussianLogLikelihood, self).__init__(problem) self._logpi = -0.5 * self.n_time_data * np.log(2 * np.pi) self._dl = np.ones(self.n_parameters + self.n_outputs) - self.transformation.append( - IdentityTransformation() - ) # Temporary fix, ahead of #338 + if self.transformation: + self.transformation.append( + IdentityTransformation() + ) # Temporary fix, ahead of #338 def _evaluate(self, x): """ diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 1d45b04b..525ab9f5 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -34,17 +34,20 @@ def __init__(self, problem=None): self.x0 = self.problem.x0 self.n_outputs = self.problem.n_outputs self.signal = self.problem.signal + self.transformation = self.construct_transformation() - # Construct ComposedTransformation from list of transformations - self.transformations = [ - t if t is not None else IdentityTransformation() - for t in self.parameters.get_transformations() - ] - self.transformation = ComposedTransformation(self.transformations) + def construct_transformation(self): + """ + Create a ComposedTransformation object from the individual parameters transformations. + """ + transformations = self.parameters.get_transformations() + if not transformations or all(t is None for t in transformations): + return None - @property - def n_parameters(self): - return len(self.parameters) + valid_transformations = [ + t if t is not None else IdentityTransformation() for t in transformations + ] + return ComposedTransformation(valid_transformations) def __call__(self, x): """ @@ -161,3 +164,7 @@ def _evaluateS1(self, x): If the method has not been implemented by the subclass. """ raise NotImplementedError + + @property + def n_parameters(self): + return len(self.parameters) diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index 5b917cea..fe4f8898 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -20,7 +20,7 @@ def parameters(self): ), pybop.Parameter( "Scaled", - transformation=pybop.ScaledTransformation(scale=2.0, translate=1), + transformation=pybop.ScaledTransformation(coefficient=2.0, intercept=1), ), pybop.Parameter( "Log", @@ -68,7 +68,7 @@ def test_scaled_transformation(self, parameters): # Test covariance transformation cov = np.array([[0.5]]) cov_transformed = transformation.convert_covariance_matrix(cov, q) - assert np.array_equal(cov_transformed, cov * transformation._scale**2) + assert np.array_equal(cov_transformed, cov * transformation._coefficient**2) @pytest.mark.unit def test_log_transformation(self, parameters): From 6c697efa3999b5823ebd33a2d3bffd7867cd113b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Mon, 24 Jun 2024 11:20:22 +0100 Subject: [PATCH 07/18] fix: old args for ScaledTransformation --- examples/scripts/spm_CMAES.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index ff4e0f3a..1fac5d4f 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -13,7 +13,7 @@ prior=pybop.Gaussian(6e-06, 0.1e-6), bounds=[1e-6, 9e-6], true_value=parameter_set["Negative particle radius [m]"], - transformation=pybop.ScaledTransformation(scale=2.0), + transformation=pybop.ScaledTransformation(coefficient=2.0), ), pybop.Parameter( "Positive particle radius [m]", From 4e1f8454f960b8464a4349478717b29ca6374b23 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 3 Jul 2024 09:51:32 +0100 Subject: [PATCH 08/18] tests: add ComposedTransformation unit tests, increase coverage --- tests/unit/test_transformations.py | 31 +++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index fe4f8898..953fee22 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -30,12 +30,13 @@ def parameters(self): @pytest.mark.unit def test_identity_transformation(self, parameters): - q = np.array([5.0]) + q = np.asarray([5.0]) transformation = parameters["Identity"].transformation assert np.array_equal(transformation.to_model(q), q) assert np.array_equal(transformation.to_search(q), q) assert transformation.log_jacobian_det(q) == 0.0 assert transformation.n_parameters == 1 + assert transformation.is_elementwise() jac, jac_S1 = transformation.jacobian_S1(q) assert np.array_equal(jac, np.eye(1)) @@ -49,11 +50,12 @@ def test_identity_transformation(self, parameters): @pytest.mark.unit def test_scaled_transformation(self, parameters): - q = np.array([2.5]) + q = np.asarray([2.5]) transformation = parameters["Scaled"].transformation p = transformation.to_model(q) assert np.allclose(p, (q / 2.0) - 1.0) assert transformation.n_parameters == 1 + assert transformation.is_elementwise() q_transformed = transformation.to_search(p) assert np.allclose(q_transformed, q) @@ -72,7 +74,7 @@ def test_scaled_transformation(self, parameters): @pytest.mark.unit def test_log_transformation(self, parameters): - q = np.array([10]) + q = np.asarray([10]) transformation = parameters["Log"].transformation p = transformation.to_model(q) assert np.allclose(p, np.exp(q)) @@ -91,6 +93,29 @@ def test_log_transformation(self, parameters): cov_transformed = transformation.convert_covariance_matrix(cov, q) assert np.array_equal(cov_transformed, cov * q**2) + @pytest.mark.unit + def test_composed_transformation(self, parameters): + q = np.asarray([5.0, 2.5]) + transformation = pybop.ComposedTransformation( + [parameters["Identity"].transformation, parameters["Scaled"].transformation] + ) + p = transformation.to_model(q) + np.testing.assert_allclose(p, np.asarray([5.0, ((q[1] / 2.0) - 1.0)])) + jac = transformation.jacobian(q) + jac_S1 = transformation.jacobian_S1(q) + np.testing.assert_allclose(jac, np.asarray([[1, 0], [0, 0.5]])) + np.testing.assert_allclose(jac_S1[0], jac) + np.testing.assert_allclose( + jac_S1[1], np.asarray([[[0, 0], [0, 0]], [[0, 0], [0, 0]]]) + ) + + transformation.append(parameters["Log"].transformation) + q = np.asarray([5.0, 2.5, 10]) + p = transformation.to_model(q) + np.testing.assert_allclose( + p, np.asarray([5.0, ((q[1] / 2.0) - 1.0), np.exp(q[2])]) + ) + class TestBaseTransformation: """ From d12173d2b757e3595b6b212550ca2b55b03fe49b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Tue, 9 Jul 2024 17:36:23 +0100 Subject: [PATCH 09/18] fix: leftover merge items --- pybop/costs/_likelihoods.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index ff15a79a..fce0608c 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -173,7 +173,7 @@ def dsigma_scale(self, new_value): raise ValueError("dsigma_scale must be non-negative") self._dsigma_scale = new_value - def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> float: + def _evaluate(self, inputs: Inputs) -> float: """ Evaluates the Gaussian log-likelihood for the given parameters. @@ -293,9 +293,6 @@ def _evaluate(self, inputs: Inputs) -> float: ---------- inputs : Inputs The parameters for which to evaluate the cost. - grad : array-like, optional - An array to store the gradient of the cost function with respect - to the parameters. Returns ------- From 87dae1c91965db05e05aac20ab8ddc13b07723f7 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 10 Jul 2024 11:13:28 +0100 Subject: [PATCH 10/18] tests: increase coverage, add condition on prior bounds creation, bugfix plotly bounds method --- pybop/parameters/parameter.py | 9 +- pybop/transformation/__init__.py | 31 ++++--- pybop/transformation/_transformation.py | 9 +- tests/unit/test_transformations.py | 111 ++++++++++++++++++++++-- 4 files changed, 131 insertions(+), 29 deletions(-) diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index c82b67cc..7fe793ff 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -172,16 +172,17 @@ def set_bounds(self, bounds=None, boundary_multiplier=6): else: self.lower_bound = bounds[0] self.upper_bound = bounds[1] - bounds = [self.lower_bound, self.upper_bound] elif self.prior is not None: self.applied_prior_bounds = True self.lower_bound = self.prior.mean - boundary_multiplier * self.prior.sigma self.upper_bound = self.prior.mean + boundary_multiplier * self.prior.sigma - bounds = [self.lower_bound, self.upper_bound] print("Default bounds applied based on prior distribution.") + else: + self.bounds = None + return - self.bounds = bounds + self.bounds = [self.lower_bound, self.upper_bound] def get_initial_value(self) -> float: """ @@ -467,7 +468,7 @@ def get_bounds_for_plotly(self): UserWarning, stacklevel=2, ) - elif param.bounds is not None: + if param.bounds is not None: bounds[i] = param.bounds else: raise ValueError("All parameters require bounds for plotting.") diff --git a/pybop/transformation/__init__.py b/pybop/transformation/__init__.py index e96e9fea..08ba5e17 100644 --- a/pybop/transformation/__init__.py +++ b/pybop/transformation/__init__.py @@ -47,7 +47,6 @@ def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarra @abstractmethod def jacobian(self, q: np.ndarray) -> np.ndarray: """Returns the Jacobian matrix of the transformation at the parameter vector `q`.""" - pass def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, Sequence[np.ndarray]]: """ @@ -89,7 +88,6 @@ def _transform(self, x: np.ndarray, method: str) -> np.ndarray: Transforms a parameter vector `x` from the search space to the model space if `method` is "to_model", or from the model space to the search space if `method` is "to_search". """ - pass def is_elementwise(self) -> bool: """ @@ -98,22 +96,23 @@ def is_elementwise(self) -> bool: """ raise NotImplementedError("is_elementwise method must be implemented if used.") - def _verify_input(self, input: Union[List[float], float, np.ndarray]) -> None: - """Set and validate the translate parameter.""" + def _verify_input(self, input: Union[float, int, List[float], np.ndarray, dict[str, float]]) -> np.ndarray: + """Set and validate the transformation parameter.""" if isinstance(input, (float, int)): - input = np.full(self._n_parameters, float(input)) + return np.full(self._n_parameters, float(input)) + if isinstance(input, dict): - if len(input) != self._n_parameters: - raise ValueError(f"Translate must have {self._n_parameters} elements") - input = np.asarray([k for k in input.values()], dtype=float) - elif isinstance(input, (list, np.ndarray)): - if len(input) != self._n_parameters: - raise ValueError(f"Translate must have {self._n_parameters} elements") - input = np.asarray(input, dtype=float) - else: - raise TypeError("Translate must be a float, list, or numpy array") - - return input + input = list(input.values()) + + try: + input_array = np.asarray(input, dtype=float) + except (ValueError, TypeError): + raise TypeError("Transform must be a float, int, list, numpy array, or dictionary") + + if input_array.size != self._n_parameters: + raise ValueError(f"Transform must have {self._n_parameters} elements") + + return input_array # ---- To be implemented with Monte Carlo PR ------ # # class TransformedLogPDF(BaseCost): diff --git a/pybop/transformation/_transformation.py b/pybop/transformation/_transformation.py index 7746ddc7..6af57260 100644 --- a/pybop/transformation/_transformation.py +++ b/pybop/transformation/_transformation.py @@ -227,15 +227,18 @@ def _append_transformation(self, transformation: Transformation): If the appended object is not a Transformation. """ if not isinstance(transformation, Transformation): - raise ValueError("The appended object must be a Transformation.") + raise TypeError("The appended object must be a Transformation.") self._transformations.append(transformation) self._n_parameters += transformation.n_parameters self._is_elementwise = self._is_elementwise and transformation.is_elementwise() - def _update_methods(self): + def _update_methods(self, elementwise: bool = None): """ Update the internal methods based on whether the transformation is elementwise. """ + if elementwise is not None: + self._is_elementwise = elementwise + if self._is_elementwise: self._jacobian = self._elementwise_jacobian self._log_jacobian_det = self._elementwise_log_jacobian_det @@ -344,7 +347,7 @@ def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: """ q = self._verify_input(q) return sum( - transformation.log_jacobian_det(q[lo : lo + transformation.n_parameters()]) + transformation.log_jacobian_det(q[lo : lo + transformation.n_parameters]) for lo, transformation in self._iter_transformations() ) diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index 953fee22..9131e622 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -35,6 +35,7 @@ def test_identity_transformation(self, parameters): assert np.array_equal(transformation.to_model(q), q) assert np.array_equal(transformation.to_search(q), q) assert transformation.log_jacobian_det(q) == 0.0 + assert transformation.log_jacobian_det_S1(q) == (0.0, np.zeros(1)) assert transformation.n_parameters == 1 assert transformation.is_elementwise() @@ -62,6 +63,9 @@ def test_scaled_transformation(self, parameters): assert np.allclose( transformation.log_jacobian_det(q), np.sum(np.log(np.abs(2.0))) ) + log_jac_det_S1 = transformation.log_jacobian_det_S1(q) + assert log_jac_det_S1[0] == np.sum(np.log(np.abs(2.0))) + assert log_jac_det_S1[1] == np.zeros(1) jac, jac_S1 = transformation.jacobian_S1(q) assert np.array_equal(jac, np.diag([0.5])) @@ -72,6 +76,10 @@ def test_scaled_transformation(self, parameters): cov_transformed = transformation.convert_covariance_matrix(cov, q) assert np.array_equal(cov_transformed, cov * transformation._coefficient**2) + # Test incorrect transform + with pytest.raises(ValueError, match="Unknown method:"): + transformation._transform(q, "bad-string") + @pytest.mark.unit def test_log_transformation(self, parameters): q = np.asarray([10]) @@ -84,6 +92,10 @@ def test_log_transformation(self, parameters): assert np.allclose(q_transformed, q) assert np.allclose(transformation.log_jacobian_det(q), -np.sum(np.log(q))) + log_jac_det_S1 = transformation.log_jacobian_det_S1(q) + assert log_jac_det_S1[0] == -np.sum(np.log(q)) + assert log_jac_det_S1[1] == -1 / q + jac, jac_S1 = transformation.jacobian_S1(q) assert np.array_equal(jac, np.diag(1 / q)) assert np.array_equal(jac_S1, np.diag(-1 / q**2)) @@ -93,12 +105,38 @@ def test_log_transformation(self, parameters): cov_transformed = transformation.convert_covariance_matrix(cov, q) assert np.array_equal(cov_transformed, cov * q**2) + # Test incorrect transform + with pytest.raises(ValueError, match="Unknown method:"): + transformation._transform(q, "bad-string") + @pytest.mark.unit def test_composed_transformation(self, parameters): - q = np.asarray([5.0, 2.5]) + # Test elementwise transformations transformation = pybop.ComposedTransformation( [parameters["Identity"].transformation, parameters["Scaled"].transformation] ) + self._test_elementwise_transformations(transformation) + + # Test appending a non-elementwise transformation + transformation.append(parameters["Log"].transformation) + self._test_non_elementwise_transformations(transformation) + + # Test composed with no transformations + with pytest.raises( + ValueError, match="Must have at least one sub-transformation." + ): + pybop.ComposedTransformation([]) + + # Test incorrect append object + with pytest.raises( + TypeError, match="The appended object must be a Transformation." + ): + pybop.ComposedTransformation( + [parameters["Identity"].transformation, "string"] + ) + + def _test_elementwise_transformations(self, transformation): + q = np.asarray([5.0, 2.5]) p = transformation.to_model(q) np.testing.assert_allclose(p, np.asarray([5.0, ((q[1] / 2.0) - 1.0)])) jac = transformation.jacobian(q) @@ -108,20 +146,77 @@ def test_composed_transformation(self, parameters): np.testing.assert_allclose( jac_S1[1], np.asarray([[[0, 0], [0, 0]], [[0, 0], [0, 0]]]) ) + assert transformation.is_elementwise() is True - transformation.append(parameters["Log"].transformation) + def _test_non_elementwise_transformations(self, transformation): q = np.asarray([5.0, 2.5, 10]) p = transformation.to_model(q) np.testing.assert_allclose( p, np.asarray([5.0, ((q[1] / 2.0) - 1.0), np.exp(q[2])]) ) + assert_log_jac_det = np.sum(np.log(np.abs(2.0))) + np.sum(-np.log(10)) + + log_jac_det = transformation.log_jacobian_det(q) + assert log_jac_det == assert_log_jac_det + + log_jac_det_S1 = transformation.log_jacobian_det_S1(q) + assert log_jac_det_S1[0] == assert_log_jac_det + np.testing.assert_allclose(log_jac_det_S1[1], np.asarray([0.0, 0.0, -0.1])) + + # Test general transformation + transformation._update_methods(elementwise=False) + assert transformation.is_elementwise() is False + + jac_general = transformation.jacobian(q) + assert np.array_equal(jac_general, np.diag([1, 0.5, 1 / 10])) + + assert_log_jac_det_general = -np.sum(np.log(np.abs(2.0))) + np.sum(-np.log(10)) + log_jac_det_general = transformation.log_jacobian_det(q) + np.testing.assert_almost_equal(log_jac_det_general, assert_log_jac_det_general) + + log_jac_det_S1_general = transformation.log_jacobian_det_S1(q) + np.testing.assert_almost_equal( + log_jac_det_S1_general[0], assert_log_jac_det_general + ) + np.testing.assert_allclose( + log_jac_det_S1_general[1], np.asarray([0.0, 0.0, -0.1]) + ) + + @pytest.mark.unit + def test_verify_input(self, parameters): + q = np.asarray([5.0]) + q_dict = {"Identity": q[0]} + transformation = parameters["Scaled"].transformation + assert transformation._verify_input(q) == q + assert transformation._verify_input(q.tolist()) == q + assert transformation._verify_input(q_dict) == q + + with pytest.raises( + TypeError, match="Transform must be a float, int, list, numpy array," + ): + transformation._verify_input("string") + + with pytest.raises(ValueError, match="Transform must have"): + transformation._verify_input(np.array([1.0, 2.0, 3.0])) + class TestBaseTransformation: """ A class to test the abstract base transformation class. """ + @pytest.fixture + def ConcreteTransformation(self): + class ConcreteTransformation(pybop.Transformation): + def jacobian(self, q): + pass + + def _transform(self, q, method): + pass + + return ConcreteTransformation() + @pytest.mark.unit def test_abstract_base_transformation(self): with pytest.raises(TypeError): @@ -148,11 +243,15 @@ def test_concrete_methods(self): assert not inspect.isabstract(getattr(pybop.Transformation, method)) @pytest.mark.unit - def test_not_implemented_methods(self): + def test_not_implemented_methods(self, ConcreteTransformation): not_implemented_methods = [ "jacobian_S1", "log_jacobian_det_S1", - "is_elementwise", ] - for method in not_implemented_methods: - assert hasattr(pybop.Transformation, method) + for method_name in not_implemented_methods: + with pytest.raises(NotImplementedError): + method = getattr(ConcreteTransformation, method_name) + method(np.array([1.0])) + + with pytest.raises(NotImplementedError): + ConcreteTransformation.is_elementwise() From 6a72568e7fb4e7237903ede011d851852c8af3cd Mon Sep 17 00:00:00 2001 From: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> Date: Wed, 10 Jul 2024 11:35:27 +0100 Subject: [PATCH 11/18] Apply suggestions from code review --- examples/scripts/spm_CMAES.py | 2 +- pybop/optimisers/base_pints_optimiser.py | 10 +++------- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index 7416f924..91c4c3c3 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -44,7 +44,7 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) cost = pybop.SumSquaredError(problem) -optim = pybop.CMAES(cost, max_iterations=100, min_iterations=40) +optim = pybop.CMAES(cost, max_iterations=100) # Run the optimisation x, final_cost = optim.run() diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index cb8abd91..8bf856c4 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -10,10 +10,7 @@ from pints import SequentialEvaluator as PintsSequentialEvaluator from pints import strfloat as PintsStrFloat -from pybop import ( - BaseOptimiser, - Result, -) +from pybop import BaseOptimiser, Result class BasePintsOptimiser(BaseOptimiser): @@ -53,7 +50,6 @@ def __init__(self, cost, pints_optimiser, **optimiser_kwargs): self.pints_optimiser = pints_optimiser super().__init__(cost, **optimiser_kwargs) - self.f = self.cost def _set_up_optimiser(self): """ @@ -197,12 +193,12 @@ def _run(self): if self._needs_sensitivities: def f(x): - L, dl = self.f.evaluateS1(x) + L, dl = self.cost.evaluateS1(x) return (L, dl) if self.minimising else (-L, -dl) else: def f(x): - return self.f(x) if self.minimising else -self.f(x) + return self.cost(x) if self.minimising else -self.cost(x) # Create evaluator object if self._parallel: From 4cec9fa3a94b8e97d6ab1d9ebfddbf127cdfe4f1 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 10 Jul 2024 12:38:09 +0100 Subject: [PATCH 12/18] refacotr: remove general transformation methods for ComposedTransformation class --- pybop/transformation/_transformation.py | 154 ++++-------------------- tests/unit/test_transformations.py | 78 ++++-------- 2 files changed, 51 insertions(+), 181 deletions(-) diff --git a/pybop/transformation/_transformation.py b/pybop/transformation/_transformation.py index 6af57260..29eb4493 100644 --- a/pybop/transformation/_transformation.py +++ b/pybop/transformation/_transformation.py @@ -210,7 +210,6 @@ def __init__(self, transformations: List[Transformation]): self._is_elementwise = True for transformation in transformations: self._append_transformation(transformation) - self._update_methods() def _append_transformation(self, transformation: Transformation): """ @@ -232,22 +231,6 @@ def _append_transformation(self, transformation: Transformation): self._n_parameters += transformation.n_parameters self._is_elementwise = self._is_elementwise and transformation.is_elementwise() - def _update_methods(self, elementwise: bool = None): - """ - Update the internal methods based on whether the transformation is elementwise. - """ - if elementwise is not None: - self._is_elementwise = elementwise - - if self._is_elementwise: - self._jacobian = self._elementwise_jacobian - self._log_jacobian_det = self._elementwise_log_jacobian_det - self._log_jacobian_det_S1 = self._elementwise_log_jacobian_det_S1 - else: - self._jacobian = self._general_jacobian - self._log_jacobian_det = self._general_log_jacobian_det - self._log_jacobian_det_S1 = self._general_log_jacobian_det_S1 - def append(self, transformation: Transformation): """ Append a new transformation to the existing composition. @@ -258,55 +241,12 @@ def append(self, transformation: Transformation): The transformation to append. """ self._append_transformation(transformation) - self._update_methods() def is_elementwise(self) -> bool: """See :meth:`Transformation.is_elementwise()`.""" return self._is_elementwise def jacobian(self, q: np.ndarray) -> np.ndarray: - """See :meth:`Transformation.jacobian()`.""" - return self._jacobian(q) - - def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: - """See :meth:`Transformation.jacobian_S1()`.""" - q = self._verify_input(q) - output_S1 = np.zeros( - (self._n_parameters, self._n_parameters, self._n_parameters) - ) - lo = 0 - - for transformation in self._transformations: - hi = lo + transformation.n_parameters - _, jac_S1 = transformation.jacobian_S1(q[lo:hi]) - for i, jac_S1_i in enumerate(jac_S1): - output_S1[lo + i, lo:hi, lo:hi] = jac_S1_i - lo = hi - - return self.jacobian(q), output_S1 - - def log_jacobian_det(self, q: np.ndarray) -> float: - """See :meth:`Transformation.log_jacobian_det()`.""" - return self._log_jacobian_det(q) - - 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_S1(q) - - def _transform(self, data: np.ndarray, method: str) -> np.ndarray: - """See :meth:`Transformation._transform`.""" - data = self._verify_input(data) - output = np.zeros_like(data) - lo = 0 - - for transformation in self._transformations: - hi = lo + transformation.n_parameters - output[lo:hi] = getattr(transformation, method)(data[lo:hi]) - lo = hi - - return output - - def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: """ Compute the elementwise Jacobian of the composed transformation. @@ -331,7 +271,24 @@ def _elementwise_jacobian(self, q: np.ndarray) -> np.ndarray: return np.diag(diag) - def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: + def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """See :meth:`Transformation.jacobian_S1()`.""" + q = self._verify_input(q) + output_S1 = np.zeros( + (self._n_parameters, self._n_parameters, self._n_parameters) + ) + lo = 0 + + for transformation in self._transformations: + hi = lo + transformation.n_parameters + _, jac_S1 = transformation.jacobian_S1(q[lo:hi]) + for i, jac_S1_i in enumerate(jac_S1): + output_S1[lo + i, lo:hi, lo:hi] = jac_S1_i + lo = hi + + return self.jacobian(q), output_S1 + + def log_jacobian_det(self, q: np.ndarray) -> float: """ Compute the elementwise logarithm of the determinant of the Jacobian. @@ -351,9 +308,7 @@ def _elementwise_log_jacobian_det(self, q: np.ndarray) -> float: for lo, transformation in self._iter_transformations() ) - def _elementwise_log_jacobian_det_S1( - self, q: np.ndarray - ) -> Tuple[float, np.ndarray]: + def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: """ Compute the elementwise logarithm of the determinant of the Jacobian and its first-order sensitivities. @@ -381,77 +336,18 @@ def _elementwise_log_jacobian_det_S1( return output, output_S1 - def _general_jacobian(self, q: np.ndarray) -> np.ndarray: - """ - Compute the general Jacobian of the composed transformation. - - Parameters - ---------- - q : np.ndarray - Input array. - - Returns - ------- - np.ndarray - Block diagonal matrix of the Jacobian of each transformation. - """ - q = self._verify_input(q) - jacobian_blocks = [] + def _transform(self, data: np.ndarray, method: str) -> np.ndarray: + """See :meth:`Transformation._transform`.""" + data = self._verify_input(data) + output = np.zeros_like(data) lo = 0 for transformation in self._transformations: hi = lo + transformation.n_parameters - jacobian_blocks.append(transformation.jacobian(q[lo:hi])) + output[lo:hi] = getattr(transformation, method)(data[lo:hi]) lo = hi - return np.block( - [ - [ - b if i == j else np.zeros_like(b) - for j, b in enumerate(jacobian_blocks) - ] - for i, _ in enumerate(jacobian_blocks) - ] - ) - - def _general_log_jacobian_det(self, q: np.ndarray) -> float: - """ - Compute the general logarithm of the determinant of the Jacobian. - - Parameters - ---------- - q : np.ndarray - Input array. - - Returns - ------- - float - Logarithm of the absolute value of the determinant of the Jacobian. - """ - return np.log(np.abs(np.linalg.det(self.jacobian(q)))) - - def _general_log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: - """ - Compute the general logarithm of the determinant of the Jacobian and its first-order sensitivities. - - Parameters - ---------- - q : np.ndarray - Input array. - - Returns - ------- - Tuple[float, np.ndarray] - Tuple of log determinant and its first-order sensitivities. - """ - q = self._verify_input(q) - jac, jac_S1 = self.jacobian_S1(q) - out_S1 = np.zeros(self._n_parameters) - - for i, jac_S1_i in enumerate(jac_S1): - out_S1[i] = np.trace(np.matmul(np.linalg.pinv(jac), jac_S1_i)) - - return self.log_jacobian_det(q), out_S1 + return output def _iter_transformations(self): """ diff --git a/tests/unit/test_transformations.py b/tests/unit/test_transformations.py index 9131e622..9a86d1f0 100644 --- a/tests/unit/test_transformations.py +++ b/tests/unit/test_transformations.py @@ -115,73 +115,47 @@ def test_composed_transformation(self, parameters): transformation = pybop.ComposedTransformation( [parameters["Identity"].transformation, parameters["Scaled"].transformation] ) - self._test_elementwise_transformations(transformation) - - # Test appending a non-elementwise transformation + # Test appending transformations transformation.append(parameters["Log"].transformation) - self._test_non_elementwise_transformations(transformation) - - # Test composed with no transformations - with pytest.raises( - ValueError, match="Must have at least one sub-transformation." - ): - pybop.ComposedTransformation([]) - - # Test incorrect append object - with pytest.raises( - TypeError, match="The appended object must be a Transformation." - ): - pybop.ComposedTransformation( - [parameters["Identity"].transformation, "string"] - ) - - def _test_elementwise_transformations(self, transformation): - q = np.asarray([5.0, 2.5]) - p = transformation.to_model(q) - np.testing.assert_allclose(p, np.asarray([5.0, ((q[1] / 2.0) - 1.0)])) - jac = transformation.jacobian(q) - jac_S1 = transformation.jacobian_S1(q) - np.testing.assert_allclose(jac, np.asarray([[1, 0], [0, 0.5]])) - np.testing.assert_allclose(jac_S1[0], jac) - np.testing.assert_allclose( - jac_S1[1], np.asarray([[[0, 0], [0, 0]], [[0, 0], [0, 0]]]) - ) assert transformation.is_elementwise() is True - def _test_non_elementwise_transformations(self, transformation): q = np.asarray([5.0, 2.5, 10]) p = transformation.to_model(q) np.testing.assert_allclose( p, np.asarray([5.0, ((q[1] / 2.0) - 1.0), np.exp(q[2])]) ) + jac = transformation.jacobian(q) + jac_S1 = transformation.jacobian_S1(q) + np.testing.assert_allclose( + jac, np.asarray([[1, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.1]]) + ) + np.testing.assert_allclose(jac_S1[0], jac) + assert jac_S1[1].shape == (3, 3, 3) + assert jac_S1[1][2, 2, 2] == -0.01 + np.testing.assert_allclose(jac_S1[1][0, :, :], np.zeros((3, 3))) + np.testing.assert_allclose(jac_S1[1][1, :, :], np.zeros((3, 3))) - assert_log_jac_det = np.sum(np.log(np.abs(2.0))) + np.sum(-np.log(10)) - + correct_output = np.sum(np.log(np.abs(2.0))) + np.sum(-np.log(10)) log_jac_det = transformation.log_jacobian_det(q) - assert log_jac_det == assert_log_jac_det + assert log_jac_det == correct_output log_jac_det_S1 = transformation.log_jacobian_det_S1(q) - assert log_jac_det_S1[0] == assert_log_jac_det + assert log_jac_det_S1[0] == correct_output np.testing.assert_allclose(log_jac_det_S1[1], np.asarray([0.0, 0.0, -0.1])) - # Test general transformation - transformation._update_methods(elementwise=False) - assert transformation.is_elementwise() is False - - jac_general = transformation.jacobian(q) - assert np.array_equal(jac_general, np.diag([1, 0.5, 1 / 10])) - - assert_log_jac_det_general = -np.sum(np.log(np.abs(2.0))) + np.sum(-np.log(10)) - log_jac_det_general = transformation.log_jacobian_det(q) - np.testing.assert_almost_equal(log_jac_det_general, assert_log_jac_det_general) + # Test composed with no transformations + with pytest.raises( + ValueError, match="Must have at least one sub-transformation." + ): + pybop.ComposedTransformation([]) - log_jac_det_S1_general = transformation.log_jacobian_det_S1(q) - np.testing.assert_almost_equal( - log_jac_det_S1_general[0], assert_log_jac_det_general - ) - np.testing.assert_allclose( - log_jac_det_S1_general[1], np.asarray([0.0, 0.0, -0.1]) - ) + # Test incorrect append object + with pytest.raises( + TypeError, match="The appended object must be a Transformation." + ): + pybop.ComposedTransformation( + [parameters["Identity"].transformation, "string"] + ) @pytest.mark.unit def test_verify_input(self, parameters): From 77a9f413ea9574285132f845f8085c10af99d1aa Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 11 Jul 2024 09:15:51 +0100 Subject: [PATCH 13/18] fix: apply ruff linting --- pybop/transformation/_transformation.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pybop/transformation/_transformation.py b/pybop/transformation/_transformation.py index 29eb4493..782fe46b 100644 --- a/pybop/transformation/_transformation.py +++ b/pybop/transformation/_transformation.py @@ -1,4 +1,4 @@ -from typing import List, Tuple, Union +from typing import Union import numpy as np @@ -43,7 +43,7 @@ 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]: + 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)) @@ -52,7 +52,7 @@ 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]: + 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) @@ -102,7 +102,7 @@ def jacobian(self, q: np.ndarray) -> np.ndarray: """See :meth:`Transformation.jacobian()`.""" return np.diag(self._inverse_coeff) - def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + 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)) @@ -111,7 +111,7 @@ def log_jacobian_det(self, q: np.ndarray) -> float: """See :meth:`Transformation.log_jacobian_det()`.""" return np.sum(np.log(np.abs(self._coefficient))) - def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + 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) @@ -164,7 +164,7 @@ def jacobian(self, q: np.ndarray) -> np.ndarray: """See :meth:`Transformation.jacobian()`.""" return np.diag(1 / q) - def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + def jacobian_S1(self, q: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """See :meth:`Transformation.jacobian_S1()`.""" return np.diag(1 / q), np.diag(-1 / q**2) @@ -172,7 +172,7 @@ def log_jacobian_det(self, q: np.ndarray) -> float: """See :meth:`Transformation.log_jacobian_det()`.""" return np.sum(-np.log(q)) - def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + 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), -1 / q @@ -202,7 +202,7 @@ class ComposedTransformation(Transformation): Extends pybop.Transformation. Initially based on pints.ComposedTransformation class. """ - def __init__(self, transformations: List[Transformation]): + def __init__(self, transformations: list[Transformation]): if not transformations: raise ValueError("Must have at least one sub-transformation.") self._transformations = [] @@ -271,7 +271,7 @@ def jacobian(self, q: np.ndarray) -> np.ndarray: return np.diag(diag) - def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + def jacobian_S1(self, q: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """See :meth:`Transformation.jacobian_S1()`.""" q = self._verify_input(q) output_S1 = np.zeros( @@ -308,7 +308,7 @@ def log_jacobian_det(self, q: np.ndarray) -> float: for lo, transformation in self._iter_transformations() ) - def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + def log_jacobian_det_S1(self, q: np.ndarray) -> tuple[float, np.ndarray]: """ Compute the elementwise logarithm of the determinant of the Jacobian and its first-order sensitivities. From 105bda25c56d97a98d4ec9fef67d97bfb38d7923 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 31 Jul 2024 17:04:37 +0100 Subject: [PATCH 14/18] fix: review suggestions. --- examples/standalone/cost.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/standalone/cost.py b/examples/standalone/cost.py index bf5d600f..1103ea89 100644 --- a/examples/standalone/cost.py +++ b/examples/standalone/cost.py @@ -54,9 +54,6 @@ def _evaluate(self, inputs): ---------- inputs : Dict The parameters for which to evaluate the cost. - grad : array-like, optional - Unused parameter, present for compatibility with gradient-based - optimizers. Returns ------- From 53d544cc08be3d3b036c2b3cf4a8602f5e0523f3 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 31 Jul 2024 17:26:39 +0100 Subject: [PATCH 15/18] fix: remainder suggestions. --- examples/scripts/spm_CMAES.py | 2 +- pybop/__init__.py | 22 ++++++++++----------- pybop/costs/base_cost.py | 37 ++++++++++------------------------- pybop/parameters/parameter.py | 15 +++++++++++++- 4 files changed, 36 insertions(+), 40 deletions(-) diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index 91c4c3c3..4376fb46 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -44,7 +44,7 @@ # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) cost = pybop.SumSquaredError(problem) -optim = pybop.CMAES(cost, max_iterations=100) +optim = pybop.CMAES(cost, sigma0=0.25, max_unchanged_iterations=20, max_iterations=50) # Run the optimisation x, final_cost = optim.run() diff --git a/pybop/__init__.py b/pybop/__init__.py index f3a7bb57..24c392c1 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -55,6 +55,17 @@ # from ._dataset import Dataset +# +# Transformation classes +# +from .transformation import Transformation +from .transformation._transformation import ( + IdentityTransformation, + ScaledTransformation, + LogTransformation, + ComposedTransformation, +) + # # Parameter classes # @@ -78,17 +89,6 @@ from .problems.fitting_problem import FittingProblem from .problems.design_problem import DesignProblem -# -# Transformation classes -# -from .transformation import Transformation -from .transformation._transformation import ( - IdentityTransformation, - ScaledTransformation, - LogTransformation, - ComposedTransformation, -) - # # Cost function class # diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 733ef202..f8630678 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -1,4 +1,4 @@ -from pybop import BaseProblem, ComposedTransformation, IdentityTransformation +from pybop import BaseProblem from pybop.parameters.parameter import Inputs, Parameters @@ -31,30 +31,13 @@ def __init__(self, problem=None): self.parameters.join(self.problem.parameters) self.n_outputs = self.problem.n_outputs self.signal = self.problem.signal - self.transformation = self.construct_transformation() - - def construct_transformation(self): - """ - Create a ComposedTransformation object from the individual parameters transformations. - """ - transformations = self.parameters.get_transformations() - if not transformations or all(t is None for t in transformations): - return None - - valid_transformations = [ - t if t is not None else IdentityTransformation() for t in transformations - ] - return ComposedTransformation(valid_transformations) + self.transformation = self.parameters.construct_transformation() def __call__(self, x): """ Call the evaluate function for a given set of parameters. """ - if self.transformation: - p = self.transformation.to_model(x) - return self.evaluate(p) - else: - return self.evaluate(x) + return self.evaluate(x) def evaluate(self, x): """ @@ -75,7 +58,9 @@ def evaluate(self, x): ValueError If an error occurs during the calculation of the cost. """ - inputs = self.parameters.verify(x) + if self.transformation: + p = self.transformation.to_model(x) + inputs = self.parameters.verify(p if self.transformation else x) try: return self._evaluate(inputs) @@ -129,14 +114,12 @@ def evaluateS1(self, x): ValueError If an error occurs during the calculation of the cost or gradient. """ - inputs = self.parameters.verify(x) + if self.transformation: + p = self.transformation.to_model(x) + inputs = self.parameters.verify(p if self.transformation else x) try: - if self.transformation: - p = self.transformation.to_model(inputs) - return self._evaluateS1(p) - else: - return self._evaluateS1(inputs) + return self._evaluateS1(inputs) except NotImplementedError as e: raise e diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 0ac47243..d219df93 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -4,6 +4,7 @@ import numpy as np +from pybop import ComposedTransformation, IdentityTransformation from pybop._utils import is_numeric Inputs = dict[str, float] @@ -57,7 +58,6 @@ def __init__( self.applied_prior_bounds = False self.set_bounds(bounds) self.margin = 1e-4 - self.set_bounds(bounds) def rvs(self, n_samples, random_state=None): """ @@ -449,6 +449,19 @@ def get_transformations(self): return transformations + def construct_transformation(self): + """ + Create a ComposedTransformation object from the individual parameter transformations. + """ + transformations = self.get_transformations() + if not transformations or all(t is None for t in transformations): + return None + + valid_transformations = [ + t if t is not None else IdentityTransformation() for t in transformations + ] + return ComposedTransformation(valid_transformations) + def get_bounds_for_plotly(self): """ Retrieve parameter bounds in the format expected by Plotly. From 0a5216c9ef537068f4ddc147eb344be514467d80 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 31 Jul 2024 21:20:57 +0100 Subject: [PATCH 16/18] tests: update transformation integration parameters, remove sigma0 value for multi_signal w/ gaussian likelihood --- pybop/costs/_weighted_cost.py | 7 ++----- tests/integration/test_spm_parameterisations.py | 1 - tests/integration/test_transformation.py | 1 + 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/pybop/costs/_weighted_cost.py b/pybop/costs/_weighted_cost.py index effa5a51..216e0303 100644 --- a/pybop/costs/_weighted_cost.py +++ b/pybop/costs/_weighted_cost.py @@ -60,7 +60,7 @@ def __init__(self, *args, weights: Optional[list[float]] = None): for cost in self.costs: self.parameters.join(cost.parameters) - def _evaluate(self, inputs: Inputs, grad=None): + def _evaluate(self, inputs: Inputs): """ Calculate the weighted cost for a given set of parameters. @@ -68,9 +68,6 @@ def _evaluate(self, inputs: Inputs, grad=None): ---------- inputs : Inputs The parameters for which to compute the cost. - grad : array-like, optional - An array to store the gradient of the cost function with respect - to the parameters. Returns ------- @@ -90,7 +87,7 @@ def _evaluate(self, inputs: Inputs, grad=None): cost._current_prediction = cost.problem.evaluate(inputs) else: cost._current_prediction = self._current_prediction - e[i] = cost._evaluate(inputs, grad) + e[i] = cost._evaluate(inputs) return np.dot(e, self.weights) diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 5e9d6b00..b3188e19 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -196,7 +196,6 @@ def test_multiple_signals(self, multi_optimiser, spm_two_signal_cost): # Test each optimiser optim = multi_optimiser( cost=spm_two_signal_cost, - sigma0=0.03, max_iterations=250, ) diff --git a/tests/integration/test_transformation.py b/tests/integration/test_transformation.py index 45feff80..4ca5a0c1 100644 --- a/tests/integration/test_transformation.py +++ b/tests/integration/test_transformation.py @@ -73,6 +73,7 @@ def test_spm_optimisers(self, optimiser, cost): x0 = cost.parameters.initial_value() optim = optimiser( cost=cost, + max_unchanged_iterations=35, max_iterations=125, ) From 54c0afe55cefd05e84b21da183d5e48ae112850f Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 1 Aug 2024 10:21:48 +0100 Subject: [PATCH 17/18] Apply suggestions from review --- pybop/parameters/parameter.py | 13 ++++++++++--- tests/unit/test_parameters.py | 17 +++++++++++++++++ 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 6976f31e..cb1a1936 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -409,9 +409,16 @@ def initial_value(self) -> np.ndarray: initial_values = [] for param in self.param.values(): - if param.initial_value is None and param.prior is not None: - initial_value = param.rvs(1)[0] - param.update(initial_value=initial_value) + if param.initial_value is None: + if param.prior is not None: + initial_value = param.rvs(1)[0] + param.update(initial_value=initial_value) + else: + warnings.warn( + "Initial value or Prior are None, proceeding without initial value.", + UserWarning, + stacklevel=2, + ) initial_values.append(param.initial_value) return np.asarray(initial_values) diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py index 90c43622..e48cda8e 100644 --- a/tests/unit/test_parameters.py +++ b/tests/unit/test_parameters.py @@ -1,3 +1,4 @@ +import numpy as np import pytest import pybop @@ -212,3 +213,19 @@ def test_get_sigma(self, parameter): parameter.prior = None assert params.get_sigma0() is None + + @pytest.mark.unit + def test_initial_values_without_attributes(self): + # Test without initial values + parameter = pybop.Parameters( + pybop.Parameter( + "Negative electrode conductivity [S.m-1]", + ) + ) + with pytest.warns( + UserWarning, + match="Initial value or Prior are None, proceeding without initial value.", + ): + sample = parameter.initial_value() + + np.testing.assert_equal(sample, np.array([None])) From 8dca6d1818192f7ec0e1b756a88a521d95f46453 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 1 Aug 2024 14:45:49 +0100 Subject: [PATCH 18/18] Update file names, add changelog entry --- CHANGELOG.md | 1 + pybop/__init__.py | 4 +- .../{__init__.py => base_transformation.py} | 42 ++++++++++++------- ...{_transformation.py => transformations.py} | 0 4 files changed, 31 insertions(+), 16 deletions(-) rename pybop/transformation/{__init__.py => base_transformation.py} (82%) rename pybop/transformation/{_transformation.py => transformations.py} (100%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 80ed158f..567f88c8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#357](https://github.com/pybop-team/PyBOP/pull/357) - Adds `Transformation()` class with `LogTransformation()`, `IdentityTransformation()`, and `ScaledTransformation()`, `ComposedTransformation()` implementations with corresponding examples and tests. - [#427](https://github.com/pybop-team/PyBOP/issues/427) - Adds the nbstripout pre-commit hook to remove unnecessary metadata from notebooks. - [#327](https://github.com/pybop-team/PyBOP/issues/327) - Adds the `WeightedCost` subclass, defines when to evaluate a problem and adds the `spm_weighted_cost` example script. - [#393](https://github.com/pybop-team/PyBOP/pull/383) - Adds Minkowski and SumofPower cost classes, with an example and corresponding tests. diff --git a/pybop/__init__.py b/pybop/__init__.py index 66c93e33..ff8c08a8 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -58,8 +58,8 @@ # # Transformation classes # -from .transformation import Transformation -from .transformation._transformation import ( +from .transformation.base_transformation import Transformation +from .transformation.transformations import ( IdentityTransformation, ScaledTransformation, LogTransformation, diff --git a/pybop/transformation/__init__.py b/pybop/transformation/base_transformation.py similarity index 82% rename from pybop/transformation/__init__.py rename to pybop/transformation/base_transformation.py index 08ba5e17..c9905775 100644 --- a/pybop/transformation/__init__.py +++ b/pybop/transformation/base_transformation.py @@ -1,5 +1,7 @@ from abc import ABC, abstractmethod -from typing import Tuple, Union, Sequence, List +from collections.abc import Sequence +from typing import Union + import numpy as np @@ -20,6 +22,7 @@ class Transformation(ABC): http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.47.9023 .. [2] Kaare Brandt Petersen and Michael Syskind Pedersen. "The Matrix Cookbook." 2012. """ + # ---- To be implemented with Monte Carlo PR ------ # # def convert_log_prior(self, log_prior): # """Returns a transformed log-prior class.""" @@ -33,7 +36,9 @@ def convert_covariance_matrix(self, cov: np.ndarray, q: np.ndarray) -> np.ndarra jac_inv = np.linalg.pinv(self.jacobian(q)) return jac_inv @ cov @ jac_inv.T - def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarray) -> np.ndarray: + def convert_standard_deviation( + self, std: Union[float, np.ndarray], q: np.ndarray + ) -> np.ndarray: """ Converts standard deviation `std`, either a scalar or a vector, from the model space to the search space around a parameter vector `q` in the search space. @@ -48,7 +53,7 @@ def convert_standard_deviation(self, std: Union[float, np.ndarray], q: np.ndarra def jacobian(self, q: np.ndarray) -> np.ndarray: """Returns the Jacobian matrix of the transformation at the parameter vector `q`.""" - def jacobian_S1(self, q: np.ndarray) -> Tuple[np.ndarray, Sequence[np.ndarray]]: + def jacobian_S1(self, q: np.ndarray) -> tuple[np.ndarray, Sequence[np.ndarray]]: """ Computes the Jacobian matrix and its partial derivatives at the parameter vector `q`. @@ -61,14 +66,18 @@ def log_jacobian_det(self, q: np.ndarray) -> float: Returns the logarithm of the absolute value of the determinant of the Jacobian matrix at the parameter vector `q`. """ - raise NotImplementedError("log_jacobian_det method must be implemented if used.") + raise NotImplementedError( + "log_jacobian_det method must be implemented if used." + ) - def log_jacobian_det_S1(self, q: np.ndarray) -> Tuple[float, np.ndarray]: + def log_jacobian_det_S1(self, q: np.ndarray) -> tuple[float, np.ndarray]: """ Computes the logarithm of the absolute value of the determinant of the Jacobian, and returns it along with its partial derivatives. """ - raise NotImplementedError("log_jacobian_det_S1 method must be implemented if used.") + raise NotImplementedError( + "log_jacobian_det_S1 method must be implemented if used." + ) @property def n_parameters(self): @@ -96,24 +105,29 @@ def is_elementwise(self) -> bool: """ raise NotImplementedError("is_elementwise method must be implemented if used.") - def _verify_input(self, input: Union[float, int, List[float], np.ndarray, dict[str, float]]) -> np.ndarray: + def _verify_input( + self, inputs: Union[float, int, list[float], np.ndarray, dict[str, float]] + ) -> np.ndarray: """Set and validate the transformation parameter.""" - if isinstance(input, (float, int)): - return np.full(self._n_parameters, float(input)) + if isinstance(inputs, (float, int)): + return np.full(self._n_parameters, float(inputs)) - if isinstance(input, dict): - input = list(input.values()) + if isinstance(inputs, dict): + inputs = list(inputs.values()) try: - input_array = np.asarray(input, dtype=float) - except (ValueError, TypeError): - raise TypeError("Transform must be a float, int, list, numpy array, or dictionary") + input_array = np.asarray(inputs, dtype=float) + except (ValueError, TypeError) as e: + raise TypeError( + "Transform must be a float, int, list, numpy array, or dictionary" + ) from e if input_array.size != self._n_parameters: raise ValueError(f"Transform must have {self._n_parameters} elements") return input_array + # ---- To be implemented with Monte Carlo PR ------ # # class TransformedLogPDF(BaseCost): # """Transformed log-PDF class.""" diff --git a/pybop/transformation/_transformation.py b/pybop/transformation/transformations.py similarity index 100% rename from pybop/transformation/_transformation.py rename to pybop/transformation/transformations.py