import warnings
import numpy as np
from scipy.optimize import minimize
from pybop.costs.base_cost import BaseCost
from pybop.costs.evaluation import Evaluation
from pybop.parameters.parameter import Inputs
from pybop.processing.dataset import Dataset
from pybop.simulators.solution import Solution
[docs]
def indices_of(values, target):
roots = (np.sign(values[1:] - target) - np.sign(values[:-1] - target)).nonzero()[0]
nearest_roots = np.unique(
np.where(
np.abs(values[roots] - target) < np.abs(values[roots + 1]), roots, roots + 1
)
)
return nearest_roots
[docs]
class FeatureDistance(BaseCost):
"""
Base for defining cost functions based on comparing fit functions.
Parameters
----------
dataset : pybop.Dataset
Dataset object containing the target data.
target : str, optional
The name of the target variable.
feature : str, optional
The name of the parameter to use for fitting, must be one of the supported features.
time_start : float, optional
Set the time (in seconds) from which onwards the data shall be
fitted, counted from the start of the data. Default is the start.
time_end : float, optional
Set the time (in seconds) until which the data shall be fitted,
counted from the start of the data. Default is the end.
"""
_supported_features = []
def __init__(
self,
dataset: Dataset,
target: str = None,
feature: str = None,
time_start: float = None,
time_end: float = None,
):
super().__init__()
if feature not in self._supported_features:
raise ValueError(
f"Feature '{feature}' not supported. Options: "
+ str(self._supported_features)
)
self.set_target(target, dataset)
if len(self._target) != 1:
raise ValueError("Feature distances require exactly one target variable.")
self.feature = feature
self.time_start = time_start
self.time_end = time_end
self.start_index = (
indices_of(self.domain_data, self.time_start)[0] if self.time_start else 0
)
self.end_index = (
indices_of(self.domain_data, self.time_end)[0] if self.time_end else None
)
with warnings.catch_warnings():
# Suppress SciPy's UserWarning about delta_grad == 0.
warnings.simplefilter("ignore")
self.data_fit = self._fit(
self.domain_data[self.start_index : self.end_index],
self.target_data[self._target[0]][self.start_index : self.end_index],
)
[docs]
def evaluate(
self,
solution: Solution,
inputs: Inputs | None = None,
calculate_sensitivities: bool = False,
) -> Evaluation:
"""Evaluate the feature distance for the given solution."""
return Evaluation(self.__call__(y=solution[self.target[0]].data))
[docs]
def _inverse_fit_function(self, y, *args):
return NotImplementedError
[docs]
def _fit_guess(self, t, y):
return NotImplementedError
[docs]
def _feature_selection(self, fit):
return NotImplementedError
[docs]
def _fit(self, t, y):
"""
Uses SciPy to fit data. For numerical reasons, the fitting
involves applying the fit function to data and comparing to identity.
"""
t = t - t[0]
fit_guess = self._fit_guess(t, y)
return self._feature_selection(
minimize(
lambda x: np.sum((t - self._inverse_fit_function(y, *x)) ** 2) ** 0.5,
x0=fit_guess,
method="trust-constr",
).x
)
[docs]
def __call__(
self,
y: np.ndarray,
dy: np.ndarray | None = None,
) -> float | tuple[float, np.ndarray]:
with warnings.catch_warnings():
# Suppress SciPy's UserWarning about delta_grad == 0.
warnings.simplefilter("ignore")
# Handle FailedSolution states
if len(y[self.start_index : self.end_index]) == 0:
return self.failure(
self.parameters.names, calculate_sensitivities=False
)
error = np.abs(
np.asarray(
[
(
self._fit(
self.domain_data[self.start_index : self.end_index],
y[self.start_index : self.end_index],
)
- self.data_fit
)
/ self.data_fit
]
)
)
return error.item()
[docs]
class SquareRootFeatureDistance(FeatureDistance):
"""
Square-root fit cost function.
Fits a square-root fit function and compares either its offset or
its slope between model predictions and target data.
Supported features:
- "offset": The value of the square-root fit at the start.
- "slope": The prefactor of the square-root over time.
- "inverse_slope": 1 over "slope"; may perform better.
"""
_supported_features = ["offset", "slope", "inverse_slope"]
def __init__(
self,
dataset: Dataset,
target: str = None,
feature: str = "inverse_slope",
time_start: float = None,
time_end: float = None,
):
super().__init__(dataset, target, feature, time_start, time_end)
[docs]
def _inverse_fit_function(self, y, b, c):
"""Square function to transform data for a linear fit."""
return ((y - b) / c) ** 2
[docs]
def _fit_guess(self, t, y):
return [y[0], (y[-1] - y[0]) / (t[-1] - t[0]) ** 0.5]
[docs]
def _feature_selection(self, fit):
if self.feature == "offset":
return fit[0]
elif self.feature == "slope":
return fit[1]
elif self.feature == "inverse_slope":
return 1 / fit[1]
[docs]
class ExponentialFeatureDistance(FeatureDistance):
"""
Exponential fit cost function.
Fits an exponential and compares either its asymptote, its magnitude,
or its timescale between model predictions and target data.
Supported features:
- "asymptote": The exponential fit value at infinite time.
- "magnitude": The prefactor of the exponential term.
- "timescale": The denominator in the exponential argument.
- "inverse_timescale": 1 over "timescale"; may perform better.
"""
_supported_features = ["asymptote", "magnitude", "timescale", "inverse_timescale"]
def __init__(
self,
dataset: Dataset,
target: str = None,
feature: str = "inverse_timescale",
time_start: float = None,
time_end: float = None,
):
super().__init__(dataset, target, feature, time_start, time_end)
[docs]
def _inverse_fit_function(self, y, b, c, d):
"""Logarithm function to transform data for a linear fit."""
log_arg = (y - b) / c
log_arg[log_arg <= 0] = 0.1**d
return -d * np.log(log_arg)
[docs]
def _fit_guess(self, t, y):
constant = y[-1]
difference = (y[-1] - y[0]) / (t[-1] - t[0])
second_difference = (y[-1] - 2 * y[len(y) // 2] + y[0]) / (
t[len(t) // 2] - t[0]
) ** 2
return [
constant,
difference**2 / second_difference,
-difference / second_difference,
]
[docs]
def _feature_selection(self, fit):
if self.feature == "asymptote":
return fit[0]
elif self.feature == "magnitude":
return fit[1]
elif self.feature == "timescale":
return fit[2]
elif self.feature == "inverse_timescale":
return 1 / fit[2]