Source code for ezyrb.regular_grid

"""
Module for higher order interpolation on regular grids
"""
import numpy as np
from scipy.interpolate import RegularGridInterpolator as RGI
from .approximation import Approximation


[docs]class RegularGrid(Approximation): """ Multidimensional interpolator on regular grids. See also scipy's RegularGridInterpolator for information on kwargs. :param array_like points: the coordinates of the points on a regular grid. :param array_like values: The (vector-) data on the regular grid in n dimensions. :Example: >>> import ezyrb >>> import numpy as np >>> 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) >>> xg, yg, zg = np.meshgrid(x, y, z, indexing='ij') >>> points = np.c_[xg.ravel(), yg.ravel(), zg.ravel()] >>> data_mode_x = f(xg, yg, zg).reshape(-1, 1) # lets assume we have 2 modes, i.e. a rank 2 model >>> data = np.concatenate((data_mode_x, data_mode_x/10), axis=1) >>> rgi = ezyrb.RegularGrid() >>> rgi.fit(points, data, method="linear") >>> pts = np.array([[2.1, 6.2, 8.3], ... [3.3, 5.2, 7.1], ... [1., 4., 7.], ... [4., 7., 9.]]) >>> rgi.predict(pts) array([[125.80469388, 12.58046939], [146.30069388, 14.63006939], [ 43. , 4.3 ], [266. , 26.6 ]]) >>> f(pts[:, 0], pts[:, 1], pts[:, 2]) array([125.542, 145.894, 43. , 266. ]) >>> f(pts[:, 0], pts[:, 1], pts[:, 2])/10 array([12.5542, 14.5894, 4.3 , 26.6 ]) """ def __init__(self): self.interpolator = None
[docs] def get_grid_axes(self, pts_scrmbld, vals_scrmbld): """ Calculates the grid axes from a meshed grid and re-orders the values so they fit on a mesh generated with numpy.meshgrid(ax1, ax2, ..., axn, indexing="ij") :param array_like pts_scrmbld: the coordinates of the points. :param array_like vals_scrmbld: the (vector-)values in the points. :return: The grid axes given as a tuple of ndarray, with shapes (m1, ), …, (mn, ) and values mapped on the ordered grid. :rtype: (list, numpy.ndarray) """ grid_axes = [] iN = [] # index in dimension N nN = [] # size of dimension N dim = pts_scrmbld.shape[1] for i in range(dim): xn, unique_inverse_n = np.unique(pts_scrmbld[:, i], return_inverse=True) grid_axes.append(xn) nN.append(len(xn)) iN.append(unique_inverse_n) if np.prod(nN) != len(vals_scrmbld): msg = "Values don't match grid. Make sure to pass a list of "\ "points that are on a regular grid! Be aware of floating "\ "point precision." raise ValueError(msg) new_row_index = np.ravel_multi_index(iN, nN) reverse_scrambling = np.argsort(new_row_index) vals_on_regular_grid = vals_scrmbld[reverse_scrambling] return grid_axes, vals_on_regular_grid
[docs] def fit(self, points, values, **kwargs): """ Construct the interpolator given `points` and `values`. Assumes that the points are on a regular grid, fails when not. see scipy.interpolate.RegularGridInterpolator :param array_like points: the coordinates of the points. :param array_like values: the values in the points. """ points = np.asarray(points) if not np.issubdtype(points.dtype, np.number): raise ValueError('Invalid format or dimension for the argument' '`points`.') if points.ndim == 1: points = points[:, None] vals = np.asarray(values) grid_axes, values_grd = self.get_grid_axes(points, vals) shape = [len(ax) for ax in grid_axes] shape.append(-1) self.interpolator = RGI(grid_axes, values_grd.reshape(shape), **kwargs)
[docs] def predict(self, new_point): """ Evaluate interpolator at given `new_point`, can be multiple points. :param array_like new_point: the coordinates of the given point(s). :return: the interpolated values. :rtype: numpy.ndarray """ new_point = np.asarray(new_point) if new_point.ndim == 1: new_point = new_point[:, None] return self.interpolator(new_point)