# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Interpolation utilities"""
import warnings
import numpy as np
import scipy.interpolate
from astropy import units as u
__all__ = [
"ScaledRegularGridInterpolator",
"interpolation_scale",
"interpolate_likelihood_profile",
]
[docs]class ScaledRegularGridInterpolator:
"""Thin wrapper around `scipy.interpolate.RegularGridInterpolator`.
The values are scaled before the interpolation and back-scaled after the
interpolation.
Parameters
----------
points : tuple of `~numpy.ndarray` or `~astropy.units.Quantity`
Tuple of points passed to `RegularGridInterpolator`.
values : `~numpy.ndarray`
Values passed to `RegularGridInterpolator`.
points_scale : tuple of str
Interpolation scale used for the points.
values_scale : {'lin', 'log', 'sqrt'}
Interpolation scaling applied to values. If the values vary over many magnitudes
a 'log' scaling is recommended.
axis : int or None
Axis along which to interpolate.
**kwargs : dict
Keyword arguments passed to `RegularGridInterpolator`.
"""
def __init__(
self,
points,
values,
points_scale=None,
values_scale="lin",
extrapolate=True,
axis=None,
**kwargs,
):
if points_scale is None:
points_scale = ["lin"] * len(points)
self.scale_points = [interpolation_scale(scale) for scale in points_scale]
self.scale = interpolation_scale(values_scale)
points_scaled = tuple([scale(p) for p, scale in zip(points, self.scale_points)])
values_scaled = self.scale(values)
self.axis = axis
if extrapolate:
kwargs.setdefault("bounds_error", False)
kwargs.setdefault("fill_value", None)
if axis is None:
self._interpolate = scipy.interpolate.RegularGridInterpolator(
points=points_scaled, values=values_scaled, **kwargs
)
else:
self._interpolate = scipy.interpolate.interp1d(
points_scaled[0], values_scaled, axis=axis
)
[docs] def __call__(self, points, method="linear", clip=True, **kwargs):
"""Interpolate data points.
Parameters
----------
points : tuple of `~numpy.ndarray` or `~astropy.units.Quantity`
Tuple of coordinate arrays of the form (x_1, x_2, x_3, ...). Arrays are
broadcasted internally.
method : {"linear", "nearest"}
Linear or nearest neighbour interpolation.
clip : bool
Clip values at zero after interpolation.
"""
points = tuple([scale(p) for scale, p in zip(self.scale_points, points)])
if self.axis is None:
points = np.broadcast_arrays(*points)
points_interp = np.stack([_.flat for _ in points]).T
values = self._interpolate(points_interp, method, **kwargs)
values = self.scale.inverse(values.reshape(points[0].shape))
else:
values = self._interpolate(points[0])
values = self.scale.inverse(values)
if clip:
values = np.clip(values, 0, np.inf)
return values
[docs]def interpolation_scale(scale="lin"):
"""Interpolation scaling.
Parameters
----------
scale : {"lin", "log", "sqrt"}
Choose interpolation scaling.
"""
if scale in ["lin", "linear"]:
return LinearScale()
elif scale == "log":
return LogScale()
elif scale == "sqrt":
return SqrtScale()
else:
raise ValueError(f"Not a valid value scaling mode: '{scale}'.")
class InterpolationScale:
"""Interpolation scale base class."""
def __call__(self, values):
if hasattr(self, "_unit"):
values = values.to_value(self._unit)
else:
if isinstance(values, u.Quantity):
self._unit = values.unit
values = values.value
return self._scale(values)
def inverse(self, values):
values = self._inverse(self, values)
if hasattr(self, "_unit"):
return u.Quantity(values, self._unit, copy=False)
else:
return values
class LogScale(InterpolationScale):
"""Logarithmic scaling"""
tiny = np.finfo(np.float32).tiny
def _scale(self, values):
values = np.clip(values, self.tiny, np.inf)
return np.log(values)
@staticmethod
def _inverse(self, values):
output = np.exp(values)
is_tiny = abs(output) - self.tiny <= self.tiny
if np.any(is_tiny):
try:
output[is_tiny] = 0.0
except (TypeError):
output = 0.0
warnings.warn(
"Interpolated values reached float32 precision limit", Warning
)
# for example TemplateSpectralModel used to define diffuse models
# could require large precision so users may want to redefine unit scaling.
return output
class SqrtScale(InterpolationScale):
"""Sqrt scaling"""
@staticmethod
def _scale(values):
sign = np.sign(values)
return sign * np.sqrt(sign * values)
@staticmethod
def _inverse(self, values):
return np.power(values, 2)
class LinearScale(InterpolationScale):
"""Linear scaling"""
@staticmethod
def _scale(values):
return values
@staticmethod
def _inverse(self, values):
return values
[docs]def interpolate_likelihood_profile(value_scan, dloglike_scan, interp_scale="sqrt"):
"""Helper function to interpolate likelihood profiles.
Parameters
----------
value_scan : `~numpy.ndarray`
Array of parameter values.
dloglike_scan : `~numpy.ndarray`
Array of delta log-likelihood values, with respect to the minimum.
interp_scale : {"sqrt", "lin"}
Interpolation scale applied to the likelihood profile. If the profile is
of parabolic shape, a "sqrt" scaling is recommended. In other cases or
for fine sampled profiles a "lin" can also be used.
Returns
-------
interp : `ScaledRegularGridInterpolator`
Interpolator instance.
"""
# likelihood profiles are typically of parabolic shape, so we use a
# sqrt scaling of the values and perform linear interpolation on the scaled
# values
sign = np.sign(np.gradient(dloglike_scan))
return ScaledRegularGridInterpolator(
points=(value_scan,), values=sign * dloglike_scan, values_scale=interp_scale
)