Module `conduction.interpolation.interpolator`

This file is part of Conduction https://git.dias.ie/itherc/conduction/

Conduction is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version.

Conduction is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with Conduction. If not, see http://www.gnu.org/licenses/.

Expand source code

"""
Copyright 2017 Ben Mather

This file is part of Conduction <https://git.dias.ie/itherc/conduction/>

Conduction is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or any later version.

Conduction is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with Conduction.  If not, see <http://www.gnu.org/licenses/>.
"""

from scipy.interpolate import RegularGridInterpolator as RGI
import itertools
import numpy as np

class RegularGridInterpolator(RGI):

    def __init__(self, points, values, method="linear", bounds_error=False, fill_value=np.nan):

        super(RegularGridInterpolator, self).__init__(points, values, method, bounds_error, fill_value)


    def adjoint(self, xi, dxi, method=None):
        """
        Interpolation adjoint using the derivatives dxi at coordinates xi
        Parameters
        ----------
        xi : ndarray of shape (..., ndim)
            The coordinates to sample the gridded data at
        dxi : ndarray of shape (..., ndim)
             The derivatives at the coordinates xi
        method : str
                The method of interplolation to perform.
                Supports either 'linear' or 'nearest'
        """
        if method is None:
            method = self.method
            
        self.derivative = dxi
        indices, norm_distances, out_of_bounds = self._find_indices(xi.T)
        if method == "linear":
            result = self._evaluate_linear_adjoint(indices, norm_distances, out_of_bounds)
        elif method == "nearest":
            result = self._evaluate_nearest_adjoint(indices, norm_distances, out_of_bounds)
        return result

    def _evaluate_linear_adjoint(self, indices, norm_distances, out_of_bounds):
        # slice for broadcasting over trailing dimensions in self.values
        vslice = (slice(None),) + (None,)*(self.values.ndim - len(indices))
        # find relevant values
        # each i and i+1 represents a edge
        edges = itertools.product(*[[i, i + 1] for i in indices])
        values = np.zeros_like(self.values)
        for edge_indices in edges:
            weight = 1.
            for ei, i, yi in zip(edge_indices, indices, norm_distances):
                weight *= np.where(ei == i, 1 - yi, yi)
            values[edge_indices] += np.asarray(self.derivative) * weight[vslice]
        return values

    def _evaluate_nearest_adjoint(self, indices, norm_distances, out_of_bounds):
        idx_res = []
        values = np.zeros_like(self.values)
        for i, yi in zip(indices, norm_distances):
            idx_res.append(np.where(yi <= 0.5, i, i + 1))
        values[idx_res] = self.derivative
        return values


class KDTreeInterpolator(object):
    """
    This object implements scipy.spatial.cKDTree for fast nearest-neighbour lookup
    and interpolates the coordinate based on the next closes value

    Arguments
    ---------
     points        : ndarray(n,dim) arbitrary set of points to interpolate over
     values        : ndarray(n,) values at the points
     bounds_error  : bool, optional
     fill_value    : number, optional

    Methods
    -------
     __call__
     adjoint

    """
    def __init__(self, points, values, bounds_error=False, fill_value=np.nan):
        from scipy.spatial import cKDTree
        self.tree = cKDTree(points)
        data = self.tree.data
        npoints, ndim = data.shape
        
        self.fill_value = fill_value
        self.bounds_error = bounds_error
        
        bbox = []
        for i in range(0, ndim):
            bbox.append((data[:,i].min(), data[:,i].max()))
        
        self.npoints = npoints
        self.ndim = ndim
        self.bbox = bbox
        
        self._values = np.ravel(values)
        
    def __call__(self, xi, *args, **kwargs):
        idx, d, bmask = self._find_indices(xi)
        
        if self.bounds_error and bmask.any():
            bidx = np.nonzero(bmask)[0]
            raise ValueError("Coordinates in xi are out of bounds in:\n {}".format(bidx))
        
        vi = self.values[idx]
        if not self.bounds_error and self.fill_value is not None:
            vi[bmask] = self.fill_value
        
        return vi
    
    def adjoint(self, xi, dxi, *args, **kwargs):
        idx, d, bmask = self._find_indices(xi)
        
        if self.bounds_error and bmask.any():
            bidx = np.nonzero(bmask)[0]
            raise ValueError("Coordinates in xi are out of bounds in:\n {}".format(bidx))
        
        dv = np.zeros_like(self.values)

        # remove indices that are out of bounds
        idx_inbounds = idx[~bmask]
        
        ux = np.unique(idx_inbounds)
        for u in ux:
            dv[u] = dxi[idx==u].sum()

        if not self.bounds_error and self.fill_value is not None:
            dv[idx][bmask] = self.fill_value
            
        return dv
    
    def _find_indices(self, xi):
        d, idx = self.tree.query(xi)
        bbox = self.bbox
        ndim = self.ndim
        
        bounds = np.zeros(idx.size, dtype=bool)
        for i in range(0, ndim):
            bmin, bmax = bbox[i]
            mask = np.logical_or(xi[:,i] < bmin, xi[:,i] > bmax)
            bounds[mask] = True
            
        return idx, d, bounds


    @property
    def values(self):
        return self._values
    @values.setter
    def values(self, value):
        self._values = value.ravel()
    @values.getter
    def values(self):
        return self._values

Classes

class KDTreeInterpolator (points, values, bounds_error=False, fill_value=nan)

This object implements scipy.spatial.cKDTree for fast nearest-neighbour lookup and interpolates the coordinate based on the next closes value

Arguments

points : ndarray(n,dim) arbitrary set of points to interpolate over values : ndarray(n,) values at the points bounds_error : bool, optional fill_value : number, optional

Methods

call adjoint

Expand source code

class KDTreeInterpolator(object):
    """
    This object implements scipy.spatial.cKDTree for fast nearest-neighbour lookup
    and interpolates the coordinate based on the next closes value

    Arguments
    ---------
     points        : ndarray(n,dim) arbitrary set of points to interpolate over
     values        : ndarray(n,) values at the points
     bounds_error  : bool, optional
     fill_value    : number, optional

    Methods
    -------
     __call__
     adjoint

    """
    def __init__(self, points, values, bounds_error=False, fill_value=np.nan):
        from scipy.spatial import cKDTree
        self.tree = cKDTree(points)
        data = self.tree.data
        npoints, ndim = data.shape
        
        self.fill_value = fill_value
        self.bounds_error = bounds_error
        
        bbox = []
        for i in range(0, ndim):
            bbox.append((data[:,i].min(), data[:,i].max()))
        
        self.npoints = npoints
        self.ndim = ndim
        self.bbox = bbox
        
        self._values = np.ravel(values)
        
    def __call__(self, xi, *args, **kwargs):
        idx, d, bmask = self._find_indices(xi)
        
        if self.bounds_error and bmask.any():
            bidx = np.nonzero(bmask)[0]
            raise ValueError("Coordinates in xi are out of bounds in:\n {}".format(bidx))
        
        vi = self.values[idx]
        if not self.bounds_error and self.fill_value is not None:
            vi[bmask] = self.fill_value
        
        return vi
    
    def adjoint(self, xi, dxi, *args, **kwargs):
        idx, d, bmask = self._find_indices(xi)
        
        if self.bounds_error and bmask.any():
            bidx = np.nonzero(bmask)[0]
            raise ValueError("Coordinates in xi are out of bounds in:\n {}".format(bidx))
        
        dv = np.zeros_like(self.values)

        # remove indices that are out of bounds
        idx_inbounds = idx[~bmask]
        
        ux = np.unique(idx_inbounds)
        for u in ux:
            dv[u] = dxi[idx==u].sum()

        if not self.bounds_error and self.fill_value is not None:
            dv[idx][bmask] = self.fill_value
            
        return dv
    
    def _find_indices(self, xi):
        d, idx = self.tree.query(xi)
        bbox = self.bbox
        ndim = self.ndim
        
        bounds = np.zeros(idx.size, dtype=bool)
        for i in range(0, ndim):
            bmin, bmax = bbox[i]
            mask = np.logical_or(xi[:,i] < bmin, xi[:,i] > bmax)
            bounds[mask] = True
            
        return idx, d, bounds


    @property
    def values(self):
        return self._values
    @values.setter
    def values(self, value):
        self._values = value.ravel()
    @values.getter
    def values(self):
        return self._values

Instance variables

var values

Expand source code

@values.getter
def values(self):
    return self._values

Methods

def adjoint(self, xi, dxi, *args, **kwargs)

Expand source code

def adjoint(self, xi, dxi, *args, **kwargs):
    idx, d, bmask = self._find_indices(xi)
    
    if self.bounds_error and bmask.any():
        bidx = np.nonzero(bmask)[0]
        raise ValueError("Coordinates in xi are out of bounds in:\n {}".format(bidx))
    
    dv = np.zeros_like(self.values)

    # remove indices that are out of bounds
    idx_inbounds = idx[~bmask]
    
    ux = np.unique(idx_inbounds)
    for u in ux:
        dv[u] = dxi[idx==u].sum()

    if not self.bounds_error and self.fill_value is not None:
        dv[idx][bmask] = self.fill_value
        
    return dv

class RegularGridInterpolator (points, values, method='linear', bounds_error=False, fill_value=nan)

Interpolation on a regular grid in arbitrary dimensions

The data must be defined on a regular grid; the grid spacing however may be uneven. Linear and nearest-neighbor interpolation are supported. After setting up the interpolator object, the interpolation method (linear or nearest) may be chosen at each evaluation.

Parameters

points : tuple of ndarray of float, with shapes (m1, ), …, (mn, ): The points defining the regular grid in n dimensions.
values : array_like, shape (m1, …, mn, …): The data on the regular grid in n dimensions.
method : str, optional: The method of interpolation to perform. Supported are "linear" and "nearest". This parameter will become the default for the object's __call__ method. 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 : number, optional: If provided, the value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated.

Methods

call

Notes

Contrary to LinearNDInterpolator and NearestNDInterpolator, this class avoids expensive triangulation of the input data by taking advantage of the regular grid structure.

If any of points have a dimension of size 1, linear interpolation will return an array of nan values. Nearest-neighbor interpolation will work as usual in this case.

Added in version: 0.14

Examples

Evaluate a simple example function on the points of a 3-D grid:

>>> from scipy.interpolate import RegularGridInterpolator
>>> def f(x, y, z):
...     return 2 * x**3 + 3 * y**2 - z
>>> x = np.linspace(1, 4, 11)
>>> y = np.linspace(4, 7, 22)
>>> z = np.linspace(7, 9, 33)
>>> data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))

data is now a 3-D array with data[i,j,k] = f(x[i], y[j], z[k]). Next, define an interpolating function from this data:

>>> my_interpolating_function = RegularGridInterpolator((x, y, z), data)

Evaluate the interpolating function at the two points (x,y,z) = (2.1, 6.2, 8.3) and (3.3, 5.2, 7.1):

>>> pts = np.array([[2.1, 6.2, 8.3], [3.3, 5.2, 7.1]])
>>> my_interpolating_function(pts)
array([ 125.80469388,  146.30069388])

which is indeed a close approximation to [f(2.1, 6.2, 8.3), f(3.3, 5.2, 7.1)].

References

.. [1] Python package regulargrid by Johannes Buchner, see https://pypi.python.org/pypi/regulargrid/ .. [2] Wikipedia, "Trilinear interpolation", https://en.wikipedia.org/wiki/Trilinear_interpolation .. [3] Weiser, Alan, and Sergio E. Zarantonello. "A note on piecewise linear and multilinear table interpolation in many dimensions." MATH. COMPUT. 50.181 (1988): 189-196. https://www.ams.org/journals/mcom/1988-50-181/S0025-5718-1988-0917826-0/S0025-5718-1988-0917826-0.pdf

Expand source code

class RegularGridInterpolator(RGI):

    def __init__(self, points, values, method="linear", bounds_error=False, fill_value=np.nan):

        super(RegularGridInterpolator, self).__init__(points, values, method, bounds_error, fill_value)


    def adjoint(self, xi, dxi, method=None):
        """
        Interpolation adjoint using the derivatives dxi at coordinates xi
        Parameters
        ----------
        xi : ndarray of shape (..., ndim)
            The coordinates to sample the gridded data at
        dxi : ndarray of shape (..., ndim)
             The derivatives at the coordinates xi
        method : str
                The method of interplolation to perform.
                Supports either 'linear' or 'nearest'
        """
        if method is None:
            method = self.method
            
        self.derivative = dxi
        indices, norm_distances, out_of_bounds = self._find_indices(xi.T)
        if method == "linear":
            result = self._evaluate_linear_adjoint(indices, norm_distances, out_of_bounds)
        elif method == "nearest":
            result = self._evaluate_nearest_adjoint(indices, norm_distances, out_of_bounds)
        return result

    def _evaluate_linear_adjoint(self, indices, norm_distances, out_of_bounds):
        # slice for broadcasting over trailing dimensions in self.values
        vslice = (slice(None),) + (None,)*(self.values.ndim - len(indices))
        # find relevant values
        # each i and i+1 represents a edge
        edges = itertools.product(*[[i, i + 1] for i in indices])
        values = np.zeros_like(self.values)
        for edge_indices in edges:
            weight = 1.
            for ei, i, yi in zip(edge_indices, indices, norm_distances):
                weight *= np.where(ei == i, 1 - yi, yi)
            values[edge_indices] += np.asarray(self.derivative) * weight[vslice]
        return values

    def _evaluate_nearest_adjoint(self, indices, norm_distances, out_of_bounds):
        idx_res = []
        values = np.zeros_like(self.values)
        for i, yi in zip(indices, norm_distances):
            idx_res.append(np.where(yi <= 0.5, i, i + 1))
        values[idx_res] = self.derivative
        return values

Ancestors

scipy.interpolate.interpolate.RegularGridInterpolator

Methods

def adjoint(self, xi, dxi, method=None)

Interpolation adjoint using the derivatives dxi at coordinates xi Parameters

xi : ndarray of shape (…, ndim): The coordinates to sample the gridded data at
dxi : ndarray of shape (…, ndim): The derivatives at the coordinates xi
method : str: The method of interplolation to perform. Supports either 'linear' or 'nearest'

Expand source code

def adjoint(self, xi, dxi, method=None):
    """
    Interpolation adjoint using the derivatives dxi at coordinates xi
    Parameters
    ----------
    xi : ndarray of shape (..., ndim)
        The coordinates to sample the gridded data at
    dxi : ndarray of shape (..., ndim)
         The derivatives at the coordinates xi
    method : str
            The method of interplolation to perform.
            Supports either 'linear' or 'nearest'
    """
    if method is None:
        method = self.method
        
    self.derivative = dxi
    indices, norm_distances, out_of_bounds = self._find_indices(xi.T)
    if method == "linear":
        result = self._evaluate_linear_adjoint(indices, norm_distances, out_of_bounds)
    elif method == "nearest":
        result = self._evaluate_nearest_adjoint(indices, norm_distances, out_of_bounds)
    return result

Classes

Arguments

Methods

Instance variables

Methods

Parameters

Methods

Notes

Examples

See Also

References

Ancestors

Methods