# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tabular models.
Tabular models of any dimension can be created using `tabular_model`.
For convenience `Tabular1D` and `Tabular2D` are provided.
Examples
--------
>>> table = np.array([[ 3., 0., 0.],
... [ 0., 2., 0.],
... [ 0., 0., 0.]])
>>> points = ([1, 2, 3], [1, 2, 3])
>>> t2 = Tabular2D(points, lookup_table=table, bounds_error=False,
... fill_value=None, method='nearest')
"""
# pylint: disable=invalid-name
import numpy as np
from astropy import units as u
from .core import Model
try:
from scipy.interpolate import interpn
has_scipy = True
except ImportError:
has_scipy = False
__all__ = ["tabular_model", "Tabular1D", "Tabular2D"]
__doctest_requires__ = {"tabular_model": ["scipy"]}
class _Tabular(Model):
"""
Returns an interpolated lookup table value.
Parameters
----------
points : tuple of ndarray of float, optional
The points defining the regular grid in n dimensions.
ndarray must have shapes (m1, ), ..., (mn, ),
lookup_table : array-like
The data on a regular grid in n dimensions.
Must have shapes (m1, ..., mn, ...)
method : str, optional
The method of interpolation to perform. Supported are "linear" and
"nearest", and "splinef2d". "splinef2d" is only supported for
2-dimensional data. Default is "linear".
bounds_error : bool, optional
If True, when interpolated values are requested outside of the
domain of the input data, a ValueError is raised.
If False, then ``fill_value`` is used.
fill_value : float or `~astropy.units.Quantity`, optional
If provided, the value to use for points outside of the
interpolation domain. If None, values outside
the domain are extrapolated. Extrapolation is not supported by method
"splinef2d". If Quantity is given, it will be converted to the unit of
``lookup_table``, if applicable.
Returns
-------
value : ndarray
Interpolated values at input coordinates.
Raises
------
ImportError
Scipy is not installed.
Notes
-----
Uses `scipy.interpolate.interpn`.
"""
linear = False
fittable = False
standard_broadcasting = False
_is_dynamic = True
_id = 0
def __init__(
self,
points=None,
lookup_table=None,
method="linear",
bounds_error=True,
fill_value=np.nan,
**kwargs,
):
n_models = kwargs.get("n_models", 1)
if n_models > 1:
raise NotImplementedError("Only n_models=1 is supported.")
super().__init__(**kwargs)
self.outputs = ("y",)
if lookup_table is None:
raise ValueError("Must provide a lookup table.")
if not isinstance(lookup_table, u.Quantity):
lookup_table = np.asarray(lookup_table)
if self.lookup_table.ndim != lookup_table.ndim:
raise ValueError(
"lookup_table should be an array with "
f"{self.lookup_table.ndim} dimensions."
)
if points is None:
points = tuple(np.arange(x, dtype=float) for x in lookup_table.shape)
else:
if lookup_table.ndim == 1 and not isinstance(points, tuple):
points = (points,)
npts = len(points)
if npts != lookup_table.ndim:
raise ValueError(
"Expected grid points in "
f"{lookup_table.ndim} directions, got {npts}."
)
if (
npts > 1
and isinstance(points[0], u.Quantity)
and len({getattr(p, "unit", None) for p in points}) > 1
):
raise ValueError("points must all have the same unit.")
if isinstance(fill_value, u.Quantity):
if not isinstance(lookup_table, u.Quantity):
raise ValueError(
f"fill value is in {fill_value.unit} but expected to be unitless."
)
fill_value = fill_value.to(lookup_table.unit).value
self.points = points
self.lookup_table = lookup_table
self.bounds_error = bounds_error
self.method = method
self.fill_value = fill_value
def __repr__(self):
return (
f"<{self.__class__.__name__}(points={self.points}, "
f"lookup_table={self.lookup_table})>"
)
def __str__(self):
default_keywords = [
("Model", self.__class__.__name__),
("Name", self.name),
("N_inputs", self.n_inputs),
("N_outputs", self.n_outputs),
("Parameters", ""),
(" points", self.points),
(" lookup_table", self.lookup_table),
(" method", self.method),
(" fill_value", self.fill_value),
(" bounds_error", self.bounds_error),
]
parts = [
f"{keyword}: {value}"
for keyword, value in default_keywords
if value is not None
]
return "\n".join(parts)
@property
def input_units(self):
pts = self.points[0]
if not isinstance(pts, u.Quantity):
return None
return {x: pts.unit for x in self.inputs}
@property
def return_units(self):
if not isinstance(self.lookup_table, u.Quantity):
return None
return {self.outputs[0]: self.lookup_table.unit}
@property
def bounding_box(self):
"""
Tuple defining the default ``bounding_box`` limits,
``(points_low, points_high)``.
Examples
--------
>>> from astropy.modeling.models import Tabular1D, Tabular2D
>>> t1 = Tabular1D(points=[1, 2, 3], lookup_table=[10, 20, 30])
>>> t1.bounding_box
ModelBoundingBox(
intervals={
x: Interval(lower=1, upper=3)
}
model=Tabular1D(inputs=('x',))
order='C'
)
>>> t2 = Tabular2D(points=[[1, 2, 3], [2, 3, 4]],
... lookup_table=[[10, 20, 30], [20, 30, 40]])
>>> t2.bounding_box
ModelBoundingBox(
intervals={
x: Interval(lower=1, upper=3)
y: Interval(lower=2, upper=4)
}
model=Tabular2D(inputs=('x', 'y'))
order='C'
)
"""
bbox = [(min(p), max(p)) for p in self.points][::-1]
if len(bbox) == 1:
bbox = bbox[0]
return bbox
def evaluate(self, *inputs):
"""
Return the interpolated values at the input coordinates.
Parameters
----------
inputs : list of scalar or list of ndarray
Input coordinates. The number of inputs must be equal
to the dimensions of the lookup table.
"""
inputs = np.broadcast_arrays(*inputs)
shape = inputs[0].shape
inputs = [inp.flatten() for inp in inputs[: self.n_inputs]]
inputs = np.array(inputs).T
if not has_scipy: # pragma: no cover
raise ImportError("Tabular model requires scipy.")
result = interpn(
self.points,
self.lookup_table,
inputs,
method=self.method,
bounds_error=self.bounds_error,
fill_value=self.fill_value,
)
# return_units not respected when points has no units
if isinstance(self.lookup_table, u.Quantity) and not isinstance(
self.points[0], u.Quantity
):
result = result * self.lookup_table.unit
if self.n_outputs == 1:
result = result.reshape(shape)
else:
result = [r.reshape(shape) for r in result]
return result
@property
def inverse(self):
if self.n_inputs == 1:
# If the wavelength array is descending instead of ascending, both
# points and lookup_table need to be reversed in the inverse transform
# for scipy.interpolate to work properly
if np.all(np.diff(self.lookup_table) > 0):
# ascending case
points = self.lookup_table
lookup_table = self.points[0]
elif np.all(np.diff(self.lookup_table) < 0):
# descending case, reverse order
points = self.lookup_table[::-1]
lookup_table = self.points[0][::-1]
else:
# equal-valued or double-valued lookup_table
raise NotImplementedError
return Tabular1D(
points=points,
lookup_table=lookup_table,
method=self.method,
bounds_error=self.bounds_error,
fill_value=self.fill_value,
)
raise NotImplementedError(
"An analytical inverse transform has not been implemented for this model."
)
[docs]def tabular_model(dim, name=None):
"""
Make a ``Tabular`` model where ``n_inputs`` is
based on the dimension of the lookup_table.
This model has to be further initialized and when evaluated
returns the interpolated values.
Parameters
----------
dim : int
Dimensions of the lookup table.
name : str
Name for the class.
Examples
--------
>>> table = np.array([[3., 0., 0.],
... [0., 2., 0.],
... [0., 0., 0.]])
>>> tab = tabular_model(2, name='Tabular2D')
>>> print(tab)
<class 'astropy.modeling.tabular.Tabular2D'>
Name: Tabular2D
N_inputs: 2
N_outputs: 1
>>> points = ([1, 2, 3], [1, 2, 3])
Setting fill_value to None, allows extrapolation.
>>> m = tab(points, lookup_table=table, name='my_table',
... bounds_error=False, fill_value=None, method='nearest')
>>> xinterp = [0, 1, 1.5, 2.72, 3.14]
>>> m(xinterp, xinterp) # doctest: +FLOAT_CMP
array([3., 3., 3., 0., 0.])
"""
if dim < 1:
raise ValueError("Lookup table must have at least one dimension.")
table = np.zeros([2] * dim)
members = {"lookup_table": table, "n_inputs": dim, "n_outputs": 1}
if dim == 1:
members["_separable"] = True
else:
members["_separable"] = False
if name is None:
model_id = _Tabular._id
_Tabular._id += 1
name = f"Tabular{model_id}"
model_class = type(str(name), (_Tabular,), members)
model_class.__module__ = "astropy.modeling.tabular"
return model_class
Tabular1D = tabular_model(1, name="Tabular1D")
Tabular2D = tabular_model(2, name="Tabular2D")
_tab_docs = """
method : str, optional
The method of interpolation to perform. Supported are "linear" and
"nearest", and "splinef2d". "splinef2d" is only supported for
2-dimensional data. Default is "linear".
bounds_error : bool, optional
If True, when interpolated values are requested outside of the
domain of the input data, a ValueError is raised.
If False, then ``fill_value`` is used.
fill_value : float, optional
If provided, the value to use for points outside of the
interpolation domain. If None, values outside
the domain are extrapolated. Extrapolation is not supported by method
"splinef2d".
Returns
-------
value : ndarray
Interpolated values at input coordinates.
Raises
------
ImportError
Scipy is not installed.
Notes
-----
Uses `scipy.interpolate.interpn`.
"""
Tabular1D.__doc__ = (
"""
Tabular model in 1D.
Returns an interpolated lookup table value.
Parameters
----------
points : array-like of float of ndim=1.
The points defining the regular grid in n dimensions.
lookup_table : array-like, of ndim=1.
The data in one dimensions.
"""
+ _tab_docs
)
Tabular2D.__doc__ = (
"""
Tabular model in 2D.
Returns an interpolated lookup table value.
Parameters
----------
points : tuple of ndarray of float, optional
The points defining the regular grid in n dimensions.
ndarray with shapes (m1, m2).
lookup_table : array-like
The data on a regular grid in 2 dimensions.
Shape (m1, m2).
"""
+ _tab_docs
)