Source code for colour.algebra.interpolation

# -*- coding: utf-8 -*-
"""
Interpolation
=============

Defines classes for interpolating variables.

-   :class:`colour.KernelInterpolator`: 1-D function generic interpolation with
    arbitrary kernel.
-   :class:`colour.LinearInterpolator`: 1-D function linear interpolation.
-   :class:`colour.SpragueInterpolator`: 1-D function fifth-order polynomial
    interpolation using *Sprague (1880)* method.
-   :class:`colour.CubicSplineInterpolator`: 1-D function cubic spline
    interpolation.
-   :class:`colour.PchipInterpolator`: 1-D function piecewise cube Hermite
    interpolation.
-   :class:`colour.NullInterpolator`: 1-D function null interpolation.
-   :func:`colour.lagrange_coefficients`: Computation of
    *Lagrange Coefficients*.

References
----------
-   :cite:`Burger2009b` : Burger, W., & Burge, M. J. (2009). Principles of
    Digital Image Processing. London: Springer London.
    doi:10.1007/978-1-84800-195-4
-   :cite:`CIETC1-382005f` : CIE TC 1-38. (2005). 9.2.4 Method of
    interpolation for uniformly spaced independent variable. In CIE 167:2005
    Recommended Practice for Tabulating Spectral Data for Use in Colour
    Computations (pp. 1-27). ISBN:978-3-901-90641-1
-   :cite:`CIETC1-382005h` : CIE TC 1-38. (2005). Table V. Values of the
    c-coefficients of Equ.s 6 and 7. In CIE 167:2005 Recommended Practice for
    Tabulating Spectral Data for Use in Colour Computations (p. 19).
    ISBN:978-3-901-90641-1
-   :cite:`Fairman1985b` : Fairman, H. S. (1985). The calculation of weight
    factors for tristimulus integration. Color Research & Application, 10(4),
    199-203. doi:10.1002/col.5080100407
-   :cite:`Westland2012h` : Westland, S., Ripamonti, C., & Cheung, V. (2012).
    Interpolation Methods. In Computational Colour Science Using MATLAB
    (2nd ed., pp. 29-37). ISBN:978-0-470-66569-5
-   :cite:`Wikipedia` : Wikipedia. (n.d.). Lanczos resampling. Retrieved
    October 14, 2017, from https://en.wikipedia.org/wiki/Lanczos_resampling
-   :cite:`Wikipediacb` : Wikipedia. (n.d.). Lagrange polynomial - Definition.
    Retrieved January 20, 2016, from https://en.wikipedia.org/wiki/\
Lagrange_polynomial#Definition
"""

from __future__ import division, unicode_literals

import numpy as np
import scipy.interpolate
from collections import OrderedDict, Mapping
from six.moves import reduce

from colour.constants import DEFAULT_FLOAT_DTYPE
from colour.utilities import (as_numeric, interval, is_integer, is_numeric,
                              closest_indexes, warning)

__author__ = 'Colour Developers'
__copyright__ = 'Copyright (C) 2013-2018 - Colour Developers'
__license__ = 'New BSD License - http://opensource.org/licenses/BSD-3-Clause'
__maintainer__ = 'Colour Developers'
__email__ = 'colour-science@googlegroups.com'
__status__ = 'Production'

__all__ = [
    'kernel_nearest_neighbour', 'kernel_linear', 'kernel_sinc',
    'kernel_lanczos', 'kernel_cardinal_spline', 'KernelInterpolator',
    'LinearInterpolator', 'SpragueInterpolator', 'CubicSplineInterpolator',
    'PchipInterpolator', 'NullInterpolator', 'lagrange_coefficients'
]


[docs]def kernel_nearest_neighbour(x): """ Returns the *nearest-neighbour* kernel evaluated at given samples. Parameters ---------- x : array_like Samples at which to evaluate the *nearest-neighbour* kernel. Returns ------- ndarray The *nearest-neighbour* kernel evaluated at given samples. References ---------- - :cite:`Burger2009b` Examples -------- >>> kernel_nearest_neighbour(np.linspace(0, 1, 10)) array([1, 1, 1, 1, 1, 0, 0, 0, 0, 0]) """ return np.where(np.abs(x) < 0.5, 1, 0)
[docs]def kernel_linear(x): """ Returns the *linear* kernel evaluated at given samples. Parameters ---------- x : array_like Samples at which to evaluate the *linear* kernel. Returns ------- ndarray The *linear* kernel evaluated at given samples. References ---------- - :cite:`Burger2009b` Examples -------- >>> kernel_linear(np.linspace(0, 1, 10)) # doctest: +ELLIPSIS array([ 1. , 0.8888888..., 0.7777777..., \ 0.6666666..., 0.5555555..., 0.4444444..., 0.3333333..., 0.2222222..., \ 0.1111111..., 0. ]) """ return np.where(np.abs(x) < 1, 1 - np.abs(x), 0)
[docs]def kernel_sinc(x, a=3): """ Returns the *sinc* kernel evaluated at given samples. Parameters ---------- x : array_like Samples at which to evaluate the *sinc* kernel. a : int, optional Size of the *sinc* kernel. Returns ------- ndarray The *sinc* kernel evaluated at given samples. References ---------- - :cite:`Burger2009b` Examples -------- >>> kernel_sinc(np.linspace(0, 1, 10)) # doctest: +ELLIPSIS array([ 1.0000000...e+00, 9.7981553...e-01, 9.2072542...e-01, 8.2699334...e-01, 7.0531659...e-01, 5.6425327...e-01, 4.1349667...e-01, 2.6306440...e-01, 1.2247694...e-01, 3.8981718...e-17]) """ assert a >= 1, '"a" must be equal or superior to 1!' return np.where(np.abs(x) < a, np.sinc(x), 0)
[docs]def kernel_lanczos(x, a=3): """ Returns the *lanczos* kernel evaluated at given samples. Parameters ---------- x : array_like Samples at which to evaluate the *lanczos* kernel. a : int, optional Size of the *lanczos* kernel. Returns ------- ndarray The *lanczos* kernel evaluated at given samples. References ---------- - :cite:`Wikipedia` Examples -------- >>> kernel_lanczos(np.linspace(0, 1, 10)) # doctest: +ELLIPSIS array([ 1.0000000...e+00, 9.7760615...e-01, 9.1243770...e-01, 8.1030092...e-01, 6.8012706...e-01, 5.3295773...e-01, 3.8071690...e-01, 2.3492839...e-01, 1.0554054...e-01, 3.2237621...e-17]) """ assert a >= 1, '"a" must be equal or superior to 1!' return np.where(np.abs(x) < a, np.sinc(x) * np.sinc(x / a), 0)
[docs]def kernel_cardinal_spline(x, a=0.5, b=0.0): """ Returns the *cardinal spline* kernel evaluated at given samples. Notable *cardinal spline* :math:`a` and :math:`b` parameterizations: - *Catmull-Rom*: :math:`(a=0.5, b=0)` - *Cubic B-Spline*: :math:`(a=0, b=1)` - *Mitchell-Netravalli*: :math:`(a=\cfrac{1}{3}, b=\cfrac{1}{3})` Parameters ---------- x : array_like Samples at which to evaluate the *cardinal spline* kernel. a : int, optional :math:`a` control parameter. b : int, optional :math:`b` control parameter. Returns ------- ndarray The *cardinal spline* kernel evaluated at given samples. References ---------- - :cite:`Burger2009b` Examples -------- >>> kernel_cardinal_spline(np.linspace(0, 1, 10)) # doctest: +ELLIPSIS array([ 1. , 0.9711934..., 0.8930041..., \ 0.7777777..., 0.6378600..., 0.4855967..., 0.3333333..., 0.1934156..., \ 0.0781893..., 0. ]) """ x_abs = np.abs(x) y = np.where(x_abs < 1, (-6 * a - 9 * b + 12) * x_abs ** 3 + (6 * a + 12 * b - 18) * x_abs ** 2 - 2 * b + 6, (-6 * a - b) * x_abs ** 3 + (30 * a + 6 * b) * x_abs ** 2 + (-48 * a - 12 * b) * x_abs + 24 * a + 8 * b) y[x_abs >= 2] = 0 return 1 / 6 * y
[docs]class KernelInterpolator(object): """ Kernel based interpolation of a 1-D function. The reconstruction of a continuous signal can be described as a linear convolution operation. Interpolation can be expressed as a convolution of the given discrete function :math:`g(x)` with some continuous interpolation kernel :math:`k(w)`: :math:`\hat{g}(w_0) = [k * g](w_0) = \ \sum_{x=-\infty}^{\infty}k(w_0 - x)\cdot g(x)` Parameters ---------- x : array_like Independent :math:`x` variable values corresponding with :math:`y` variable. y : array_like Dependent and already known :math:`y` variable values to interpolate. window : int, optional Width of the window in samples on each side. kernel : callable, optional Kernel to use for interpolation. kernel_args : dict, optional Arguments to use when calling the kernel. padding_args : dict, optional Arguments to use when padding :math:`y` variable values with the :func:`np.pad` definition. dtype : type Data type used for internal conversions. Attributes ---------- x y window kernel kernel_args padding_args Methods ------- __call__ References ---------- - :cite:`Burger2009b` - :cite:`Wikipedia` Examples -------- Interpolating a single numeric variable: >>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100, ... 69.5900, 27.8007, 86.0500]) >>> x = np.arange(len(y)) >>> f = KernelInterpolator(x, y) >>> f(0.5) # doctest: +ELLIPSIS 6.9411400... Interpolating an *array_like* variable: >>> f([0.25, 0.75]) # doctest: +ELLIPSIS array([ 6.1806208..., 8.0823848...]) Using a different *lanczos* kernel: >>> f = KernelInterpolator(x, y, kernel=kernel_sinc) >>> f([0.25, 0.75]) # doctest: +ELLIPSIS array([ 6.5147317..., 8.3965466...]) Using a different window size: >>> f = KernelInterpolator( ... x, ... y, ... window=16, ... kernel=kernel_lanczos, ... kernel_args={'a': 16}) >>> f([0.25, 0.75]) # doctest: +ELLIPSIS array([ 5.396179..., 5.652109...]) """
[docs] def __init__(self, x, y, window=3, kernel=kernel_lanczos, kernel_args=None, padding_args=None, dtype=DEFAULT_FLOAT_DTYPE): self._x_p = None self._y_p = None self._x = None self._y = None self._window = None self._padding_args = {'pad_width': (window, window), 'mode': 'reflect'} self._dtype = dtype self.x = x self.y = y self.window = window self.padding_args = padding_args self._kernel = None self.kernel = kernel self._kernel_args = {} self.kernel_args = kernel_args self._validate_dimensions()
@property def x(self): """ Getter and setter property for the independent :math:`x` variable. Parameters ---------- value : array_like Value to set the independent :math:`x` variable with. Returns ------- array_like Independent :math:`x` variable. """ return self._x @x.setter def x(self, value): """ Setter for the **self.x** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"x" independent variable must have exactly one dimension!') value_interval = interval(value) if value_interval.size != 1: warning(('"x" independent variable is not uniform, ' 'unpredictable results may occur!')) self._x = value if self._window is not None: self._x_p = np.pad( self._x, (self._window, self._window), 'linear_ramp', end_values=( np.min(self._x) - self._window * value_interval[0], np.max(self._x) + self._window * value_interval[0])) @property def y(self): """ Getter and setter property for the dependent and already known :math:`y` variable. Parameters ---------- value : array_like Value to set the dependent and already known :math:`y` variable with. Returns ------- array_like Dependent and already known :math:`y` variable. """ return self._y @y.setter def y(self, value): """ Setter for the **self.y** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"y" dependent variable must have exactly one dimension!') self._y = value if self._window is not None: self._y_p = np.pad(self._y, **self._padding_args) @property def window(self): """ Getter and setter property for the window. Parameters ---------- value : int Value to set the window with. Returns ------- int Window. """ return self._window @window.setter def window(self, value): """ Setter for the **self.window** property. """ if value is not None: assert is_integer(value), '"window" must be an integer!' assert value >= 1, '"window" must be equal or superior to 1!' self._window = value # Triggering "self._x_p" update. if self._x is not None: self.x = self._x # Triggering "self._y_p" update. if self._y is not None: self.y = self._y @property def kernel(self): """ Getter and setter property for the kernel callable. Parameters ---------- value : callable Value to set the kernel callable. Returns ------- callable Kernel callable. """ return self._kernel @kernel.setter def kernel(self, value): """ Setter for the **self.kernel** property. """ if value is not None: assert hasattr( value, '__call__'), ('"{0}" attribute: "{1}" is not callable!'.format( 'kernel', value)) self._kernel = value @property def kernel_args(self): """ Getter and setter property for the kernel call time arguments. Parameters ---------- value : dict Value to call the interpolation kernel with. Returns ------- dict Kernel call time arguments. """ return self._kernel_args @kernel_args.setter def kernel_args(self, value): """ Setter for the **self.kernel_args** property. """ if value is not None: assert isinstance(value, (dict, OrderedDict)), ( '"{0}" attribute: "{1}" type is not "dict" or "OrderedDict"!' ).format('kernel_args', value) self._kernel_args = value @property def padding_args(self): """ Getter and setter property for the kernel call time arguments. Parameters ---------- value : dict Value to call the interpolation kernel with. Returns ------- dict Kernel call time arguments. """ return self._padding_args @padding_args.setter def padding_args(self, value): """ Setter for the **self.padding_args** property. """ if value is not None: assert isinstance(value, Mapping), ( '"{0}" attribute: "{1}" type is not a "Mapping" instance!' ).format('padding_args', value) self._padding_args = value # Triggering "self._y_p" update. if self._y is not None: self.y = self._y
[docs] def __call__(self, x): """ Evaluates the interpolator at given point(s). Parameters ---------- x : numeric or array_like Point(s) to evaluate the interpolant at. Returns ------- float or ndarray Interpolated value(s). """ x = np.atleast_1d(x).astype(self._dtype) xi = as_numeric(self._evaluate(x)) return xi
def _evaluate(self, x): """ Performs the interpolator evaluation at given points. Parameters ---------- x : ndarray Points to evaluate the interpolant at. Returns ------- ndarray Interpolated points values. """ self._validate_dimensions() self._validate_interpolation_range(x) x_interval = interval(self._x)[0] x_f = np.floor(x / x_interval) windows = (x_f[:, np.newaxis] + np.arange(-self._window + 1, self._window + 1)) clip_l = min(self._x_p) / x_interval clip_h = max(self._x_p) / x_interval windows = np.clip(windows, clip_l, clip_h) - clip_l windows = np.around(windows).astype(np.int_) return np.sum( self._y_p[windows] * self._kernel(x[:, np.newaxis] / x_interval - windows - min(self._x_p) / x_interval, **self._kernel_args), axis=-1) def _validate_dimensions(self): """ Validates variables dimensions to be the same. """ if len(self._x) != len(self._y): raise ValueError( ('"x" independent and "y" dependent variables have different ' 'dimensions: "{0}", "{1}"').format( len(self._x), len(self._y))) def _validate_interpolation_range(self, x): """ Validates given point to be in interpolation range. """ below_interpolation_range = x < self._x[0] above_interpolation_range = x > self._x[-1] if below_interpolation_range.any(): raise ValueError('"{0}" is below interpolation range.'.format(x)) if above_interpolation_range.any(): raise ValueError('"{0}" is above interpolation range.'.format(x))
[docs]class LinearInterpolator(object): """ Linearly interpolates a 1-D function. Parameters ---------- x : array_like Independent :math:`x` variable values corresponding with :math:`y` variable. y : array_like Dependent and already known :math:`y` variable values to interpolate. dtype : type Data type used for internal conversions. Attributes ---------- x y Methods ------- __call__ Notes ----- This class is a wrapper around *numpy.interp* definition. Examples -------- Interpolating a single numeric variable: >>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100, ... 69.5900, 27.8007, 86.0500]) >>> x = np.arange(len(y)) >>> f = LinearInterpolator(x, y) >>> # Doctests ellipsis for Python 2.x compatibility. >>> f(0.5) # doctest: +ELLIPSIS 7.64... Interpolating an *array_like* variable: >>> f([0.25, 0.75]) array([ 6.7825, 8.5075]) """
[docs] def __init__(self, x, y, dtype=DEFAULT_FLOAT_DTYPE): self._x = None self._y = None self._dtype = dtype self.x = x self.y = y self._validate_dimensions()
@property def x(self): """ Getter and setter property for the independent :math:`x` variable. Parameters ---------- value : array_like Value to set the independent :math:`x` variable with. Returns ------- array_like Independent :math:`x` variable. """ return self._x @x.setter def x(self, value): """ Setter for the **self.x** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"x" independent variable must have exactly one dimension!') self._x = value @property def y(self): """ Getter and setter property for the dependent and already known :math:`y` variable. Parameters ---------- value : array_like Value to set the dependent and already known :math:`y` variable with. Returns ------- array_like Dependent and already known :math:`y` variable. """ return self._y @y.setter def y(self, value): """ Setter for the **self.y** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"y" dependent variable must have exactly one dimension!') self._y = value
[docs] def __call__(self, x): """ Evaluates the interpolating polynomial at given point(s). Parameters ---------- x : numeric or array_like Point(s) to evaluate the interpolant at. Returns ------- float or ndarray Interpolated value(s). """ x = np.atleast_1d(x).astype(self._dtype) xi = as_numeric(self._evaluate(x)) return xi
def _evaluate(self, x): """ Performs the interpolating polynomial evaluation at given points. Parameters ---------- x : ndarray Points to evaluate the interpolant at. Returns ------- ndarray Interpolated points values. """ self._validate_dimensions() self._validate_interpolation_range(x) return np.interp(x, self._x, self._y) def _validate_dimensions(self): """ Validates variables dimensions to be the same. """ if len(self._x) != len(self._y): raise ValueError( ('"x" independent and "y" dependent variables have different ' 'dimensions: "{0}", "{1}"').format( len(self._x), len(self._y))) def _validate_interpolation_range(self, x): """ Validates given point to be in interpolation range. """ below_interpolation_range = x < self._x[0] above_interpolation_range = x > self._x[-1] if below_interpolation_range.any(): raise ValueError('"{0}" is below interpolation range.'.format(x)) if above_interpolation_range.any(): raise ValueError('"{0}" is above interpolation range.'.format(x))
[docs]class SpragueInterpolator(object): """ Constructs a fifth-order polynomial that passes through :math:`y` dependent variable. *Sprague (1880)* method is recommended by the *CIE* for interpolating functions having a uniformly spaced independent variable. Parameters ---------- x : array_like Independent :math:`x` variable values corresponding with :math:`y` variable. y : array_like Dependent and already known :math:`y` variable values to interpolate. dtype : type Data type used for internal conversions. Attributes ---------- x y Methods ------- __call__ Notes ----- The minimum number :math:`k` of data points required along the interpolation axis is :math:`k=6`. References ---------- - :cite:`CIETC1-382005f` - :cite:`Westland2012h` Examples -------- Interpolating a single numeric variable: >>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100, ... 69.5900, 27.8007, 86.0500]) >>> x = np.arange(len(y)) >>> f = SpragueInterpolator(x, y) >>> f(0.5) # doctest: +ELLIPSIS 7.2185025... Interpolating an *array_like* variable: >>> f([0.25, 0.75]) # doctest: +ELLIPSIS array([ 6.7295161..., 7.8140625...]) """ SPRAGUE_C_COEFFICIENTS = np.array([ [884, -1960, 3033, -2648, 1080, -180], [508, -540, 488, -367, 144, -24], [-24, 144, -367, 488, -540, 508], [-180, 1080, -2648, 3033, -1960, 884], ]) """ Defines the coefficients used to generate extra points for boundaries interpolation. SPRAGUE_C_COEFFICIENTS : array_like, (4, 6) References ---------- - :cite:`CIETC1-382005h` """
[docs] def __init__(self, x, y, dtype=DEFAULT_FLOAT_DTYPE): self._xp = None self._yp = None self._x = None self._y = None self._dtype = dtype self.x = x self.y = y self._validate_dimensions()
@property def x(self): """ Getter and setter property for the independent :math:`x` variable. Parameters ---------- value : array_like Value to set the independent :math:`x` variable with. Returns ------- array_like Independent :math:`x` variable. """ return self._x @x.setter def x(self, value): """ Setter for the **self.x** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"x" independent variable must have exactly one dimension!') value_interval = interval(value)[0] xp1 = value[0] - value_interval * 2 xp2 = value[0] - value_interval xp3 = value[-1] + value_interval xp4 = value[-1] + value_interval * 2 self._xp = np.concatenate(((xp1, xp2), value, (xp3, xp4))) self._x = value @property def y(self): """ Getter and setter property for the dependent and already known :math:`y` variable. Parameters ---------- value : array_like Value to set the dependent and already known :math:`y` variable with. Returns ------- array_like Dependent and already known :math:`y` variable. """ return self._y @y.setter def y(self, value): """ Setter for the **self.y** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"y" dependent variable must have exactly one dimension!') assert len(value) >= 6, ( '"y" dependent variable values count must be in domain [6:]!') yp1 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[0], np.array(value[0:6]).reshape( (6, 1)))) / 209)[0] yp2 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[1], np.array(value[0:6]).reshape( (6, 1)))) / 209)[0] yp3 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[2], np.array(value[-6:]).reshape( (6, 1)))) / 209)[0] yp4 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[3], np.array(value[-6:]).reshape( (6, 1)))) / 209)[0] self._yp = np.concatenate(((yp1, yp2), value, (yp3, yp4))) self._y = value
[docs] def __call__(self, x): """ Evaluates the interpolating polynomial at given point(s). Parameters ---------- x : numeric or array_like Point(s) to evaluate the interpolant at. Returns ------- numeric or ndarray Interpolated value(s). """ return self._evaluate(x)
def _evaluate(self, x): """ Performs the interpolating polynomial evaluation at given point. Parameters ---------- x : numeric Point to evaluate the interpolant at. Returns ------- float Interpolated point values. """ x = np.asarray(x) self._validate_dimensions() self._validate_interpolation_range(x) i = np.searchsorted(self._xp, x) - 1 X = (x - self._xp[i]) / (self._xp[i + 1] - self._xp[i]) r = self._yp a0p = r[i] a1p = ((2 * r[i - 2] - 16 * r[i - 1] + 16 * r[i + 1] - 2 * r[i + 2]) / 24) # yapf: disable a2p = ((-r[i - 2] + 16 * r[i - 1] - 30 * r[i] + 16 * r[i + 1] - r[i + 2]) / 24) # yapf: disable a3p = ((-9 * r[i - 2] + 39 * r[i - 1] - 70 * r[i] + 66 * r[i + 1] - 33 * r[i + 2] + 7 * r[i + 3]) / 24) a4p = ((13 * r[i - 2] - 64 * r[i - 1] + 126 * r[i] - 124 * r[i + 1] + 61 * r[i + 2] - 12 * r[i + 3]) / 24) a5p = ((-5 * r[i - 2] + 25 * r[i - 1] - 50 * r[i] + 50 * r[i + 1] - 25 * r[i + 2] + 5 * r[i + 3]) / 24) y = (a0p + a1p * X + a2p * X ** 2 + a3p * X ** 3 + a4p * X ** 4 + a5p * X ** 5) return y def _validate_dimensions(self): """ Validates variables dimensions to be the same. """ if len(self._x) != len(self._y): raise ValueError( ('"x" independent and "y" dependent variables have different ' 'dimensions: "{0}", "{1}"').format( len(self._x), len(self._y))) def _validate_interpolation_range(self, x): """ Validates given point to be in interpolation range. """ below_interpolation_range = x < self._x[0] above_interpolation_range = x > self._x[-1] if below_interpolation_range.any(): raise ValueError('"{0}" is below interpolation range.'.format(x)) if above_interpolation_range.any(): raise ValueError('"{0}" is above interpolation range.'.format(x))
class CubicSplineInterpolator(scipy.interpolate.interp1d): """ Interpolates a 1-D function using cubic spline interpolation. Notes ----- This class is a wrapper around *scipy.interpolate.interp1d* class. """ def __init__(self, *args, **kwargs): super(CubicSplineInterpolator, self).__init__( kind='cubic', *args, **kwargs)
[docs]class PchipInterpolator(scipy.interpolate.PchipInterpolator): """ Interpolates a 1-D function using Piecewise Cubic Hermite Interpolating Polynomial interpolation. Attributes ---------- y Notes ----- - This class is a wrapper around *scipy.interpolate.PchipInterpolator* class. """
[docs] def __init__(self, x, y, *args, **kwargs): super(PchipInterpolator, self).__init__(x, y, *args, **kwargs) self._y = y
@property def y(self): """ Getter and setter property for the dependent and already known :math:`y` variable. Parameters ---------- value : array_like Value to set the dependent and already known :math:`y` variable with. Returns ------- array_like Dependent and already known :math:`y` variable. """ return self._y
[docs]class NullInterpolator(object): """ Performs 1-D function null interpolation, i.e. a call within given tolerances will return existing :math:`y` variable values and ``default`` if outside tolerances. Parameters ---------- x : ndarray Independent :math:`x` variable values corresponding with :math:`y` variable. y : ndarray Dependent and already known :math:`y` variable values to interpolate. absolute_tolerance : numeric, optional Absolute tolerance. relative_tolerance : numeric, optional Relative tolerance. default : numeric, optional Default value for interpolation outside tolerances. dtype : type Data type used for internal conversions. Attributes ---------- x y relative_tolerance absolute_tolerance default Methods ------- __call__ Examples -------- >>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100, ... 69.5900, 27.8007, 86.0500]) >>> x = np.arange(len(y)) >>> f = NullInterpolator(x, y) >>> f(0.5) nan >>> f(1.0) # doctest: +ELLIPSIS 9.3699999... >>> f = NullInterpolator(x, y, absolute_tolerance=0.01) >>> f(1.01) # doctest: +ELLIPSIS 9.3699999... """
[docs] def __init__(self, x, y, absolute_tolerance=10e-7, relative_tolerance=10e-7, default=np.nan, dtype=DEFAULT_FLOAT_DTYPE): self._x = None self._y = None self._absolute_tolerance = None self._relative_tolerance = None self._default = None self._dtype = dtype self.x = x self.y = y self.absolute_tolerance = absolute_tolerance self.relative_tolerance = relative_tolerance self.default = default self._validate_dimensions()
@property def x(self): """ Getter and setter property for the independent :math:`x` variable. Parameters ---------- value : array_like Value to set the independent :math:`x` variable with. Returns ------- array_like Independent :math:`x` variable. """ return self._x @x.setter def x(self, value): """ Setter for the **self.x** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"x" independent variable must have exactly one dimension!') self._x = value @property def y(self): """ Getter and setter property for the dependent and already known :math:`y` variable. Parameters ---------- value : array_like Value to set the dependent and already known :math:`y` variable with. Returns ------- array_like Dependent and already known :math:`y` variable. """ return self._y @y.setter def y(self, value): """ Setter for the **self.y** property. """ if value is not None: value = np.atleast_1d(value).astype(self._dtype) assert value.ndim == 1, ( '"y" dependent variable must have exactly one dimension!') self._y = value @property def relative_tolerance(self): """ Getter and setter property for the relative tolerance. Parameters ---------- value : numeric Value to set the relative tolerance with. Returns ------- numeric Relative tolerance. """ return self._relative_tolerance @relative_tolerance.setter def relative_tolerance(self, value): """ Setter for the **self.relative_tolerance** property. """ if value is not None: assert is_numeric(value), ( '"relative_tolerance" variable must be a "numeric"!') self._relative_tolerance = value @property def absolute_tolerance(self): """ Getter and setter property for the absolute tolerance. Parameters ---------- value : numeric Value to set the absolute tolerance with. Returns ------- numeric Absolute tolerance. """ return self._absolute_tolerance @absolute_tolerance.setter def absolute_tolerance(self, value): """ Setter for the **self.absolute_tolerance** property. """ if value is not None: assert is_numeric(value), ( '"absolute_tolerance" variable must be a "numeric"!') self._absolute_tolerance = value @property def default(self): """ Getter and setter property for the default value for call outside tolerances. Parameters ---------- value : numeric Value to set the default value with. Returns ------- numeric Default value. """ return self._default @default.setter def default(self, value): """ Setter for the **self.default** property. """ if value is not None: assert is_numeric(value), ( '"default" variable must be a "numeric"!') self._default = value
[docs] def __call__(self, x): """ Evaluates the interpolator at given point(s). Parameters ---------- x : numeric or array_like Point(s) to evaluate the interpolant at. Returns ------- float or ndarray Interpolated value(s). """ x = np.atleast_1d(x).astype(self._dtype) xi = as_numeric(self._evaluate(x)) return xi
def _evaluate(self, x): """ Performs the interpolator evaluation at given points. Parameters ---------- x : ndarray Points to evaluate the interpolant at. Returns ------- ndarray Interpolated points values. """ self._validate_dimensions() self._validate_interpolation_range(x) indexes = closest_indexes(self._x, x) values = self._y[indexes] values[~np.isclose( self._x[indexes], x, rtol=self._absolute_tolerance, atol=self._relative_tolerance)] = self._default return values def _validate_dimensions(self): """ Validates variables dimensions to be the same. """ if len(self._x) != len(self._y): raise ValueError( ('"x" independent and "y" dependent variables have different ' 'dimensions: "{0}", "{1}"').format( len(self._x), len(self._y))) def _validate_interpolation_range(self, x): """ Validates given point to be in interpolation range. """ below_interpolation_range = x < self._x[0] above_interpolation_range = x > self._x[-1] if below_interpolation_range.any(): raise ValueError('"{0}" is below interpolation range.'.format(x)) if above_interpolation_range.any(): raise ValueError('"{0}" is above interpolation range.'.format(x))
[docs]def lagrange_coefficients(r, n=4): """ Computes the *Lagrange Coefficients* at given point :math:`r` for degree :math:`n`. Parameters ---------- r : numeric Point to get the *Lagrange Coefficients* at. n : int, optional Degree of the *Lagrange Coefficients* being calculated. Returns ------- ndarray References ---------- - :cite:`Fairman1985b` - :cite:`Wikipediacb` Examples -------- >>> lagrange_coefficients(0.1) array([ 0.8265, 0.2755, -0.1305, 0.0285]) """ r_i = np.arange(n) L_n = [] for j in range(len(r_i)): basis = [(r - r_i[i]) / (r_i[j] - r_i[i]) for i in range(len(r_i)) if i != j] L_n.append(reduce(lambda x, y: x * y, basis)) # noqa return np.array(L_n)