Regular grid interpolation

docs/source/regular_grid.rst

@@ -0,0 +1,23 @@
+RegularGrid
+=====================
+
+.. currentmodule:: ezyrb.regular_grid
+
+.. automodule:: ezyrb.regular_grid
+
+.. autosummary::
+    :toctree: _summaries
+    :nosignatures:
+
+    RegularGrid
+    RegularGrid.get_grid_axes
+    RegularGrid.fit
+    RegularGrid.predict
+    RegularGrid.calculate_flat_index
+
+.. autoclass:: RegularGrid
+    :members:
+    :private-members:
+    :undoc-members:
+    :show-inheritance:
+    :noindex:

ezyrb/__init__.py

@@ -3,7 +3,7 @@
 __all__ = [
     'Database', 'Reduction', 'POD', 'Approximation', 'RBF', 'Linear', 'GPR',
     'ANN', 'KNeighborsRegressor', 'RadiusNeighborsRegressor', 'AE',
-    'ReducedOrderModel', 'PODAE'
+    'ReducedOrderModel', 'PODAE', 'RegularGrid'
 ]
 
 from .meta import *
@@ -15,6 +15,7 @@
 from .approximation import Approximation
 from .rbf import RBF
 from .linear import Linear
+from .regular_grid import RegularGrid
 from .gpr import GPR
 from .reducedordermodel import ReducedOrderModel
 from .ann import ANN

ezyrb/regular_grid.py

@@ -0,0 +1,117 @@
+"""
+Module for higher order interpolation on regular grids
+"""
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator as RGI
+from .approximation import Approximation
+
+
+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
+
+    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
+
+    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 len(points.shape) == 1:
+            points.shape = (-1, 1)
+        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)
+
+    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 len(new_point.shape) == 1:
+            new_point.shape = (-1, 1)
+        return self.interpolator(new_point)

tests/test_regular_grid.py

@@ -0,0 +1,95 @@
+import numpy as np
+from unittest import TestCase
+from ezyrb import RegularGrid, Database, POD, ReducedOrderModel
+
+
+class TestRegularGrid(TestCase):
+
+    def test_1D_1mode(self):
+        reg = RegularGrid()
+        reg.fit(np.array([1, 2, 3])[:, None],
+                np.array([1, 5, 3])[:, None])
+        assert reg.predict([1]) == 1
+        assert reg.predict([2]) == 5
+        assert reg.predict([3]) == 3
+
+    def test_predict1d(self):
+        reg = RegularGrid()
+        x1 = np.array([[1], [6], [8]])
+        V = np.array([[1, 0], [20, 5], [8, 6]])  # n, r = 3, 2
+        reg.fit(x1, V)
+        result = reg.predict([[1], [8], [6]])
+        assert (result[0] == [1, 0]).all()
+        assert (result[1] == [8, 6]).all()
+        assert (result[2] == [20, 5]).all()
+
+    def test_predict1d_2(self):
+        reg = RegularGrid()
+        x1 = [1, 2]
+        reg.fit(x1, [[1, 1], [2, 2]])
+        result = reg.predict([[1.5]])
+        assert (result[0] == [[1.5, 1.5]]).all()
+
+    def test_predict2d(self):
+        reg = RegularGrid()
+        x1 = [1, 2, 3]
+        x2 = [4, 5, 6, 7]
+        xx1, xx2 = np.meshgrid(x1, x2, indexing="ij")
+        points = np.c_[xx1.ravel(), xx2.ravel()]
+        V = [[1, 21], [2, 22], [3, 23], [4, 24], [5, 25], [6, 26],
+             [7, 27], [8, 28], [9, 29], [10, 30], [11, 31], [12, 32]]
+        reg.fit(points, V, method="linear")
+        result = reg.predict([[1, 4], [1, 5]])
+        assert (result[0] == [1, 21]).all()
+        assert (result[1] == [2, 22]).all()
+
+    def test_get_grid_axes_2D(self):
+        x1 = [.1, .2, .3]
+        x2 = [4, 5, 6, 7]
+        xx1, xx2 = np.meshgrid(x1, x2, indexing="ij")
+        V = np.arange(len(x1)*len(x2)*2).reshape(-1, 2)
+        pts = np.c_[xx1.ravel(), xx2.ravel()]
+
+        random_order = np.arange(len(pts))
+        np.random.shuffle(random_order)
+        pts_scrmbld = pts[random_order, :]
+        V_scrmbld = V[random_order, :]
+
+        reg = RegularGrid()
+        grid_axes, V_unscrambled = reg.get_grid_axes(pts_scrmbld, V_scrmbld)
+
+        assert np.allclose(V_unscrambled, V)
+        assert np.allclose(grid_axes[0], x1)
+        assert np.allclose(grid_axes[1], x2)
+
+    def test_get_grid_axes_3D(self):
+        x1 = [.1, .2, .3]
+        x2 = [4, 5, 6, 7]
+        x3 = [.12, .34, .56, .78, .90]
+        xx1, xx2, xx3 = np.meshgrid(x1, x2, x3, indexing="ij")
+        V = np.arange(len(x1)*len(x2)*len(x3)*2).reshape(-1, 2)
+        pts = np.c_[xx1.ravel(), xx2.ravel(), xx3.ravel()]
+
+        random_order = np.arange(len(pts))
+        np.random.shuffle(random_order)
+        pts_scrmbld = pts[random_order, :]
+        V_scrmbld = V[random_order, :]
+
+        reg = RegularGrid()
+        grid_axes, V_unscrambled = reg.get_grid_axes(pts_scrmbld, V_scrmbld)
+
+        assert np.allclose(V_unscrambled, V)
+        assert np.allclose(grid_axes[0], x1)
+        assert np.allclose(grid_axes[1], x2)
+        assert np.allclose(grid_axes[2], x3)
+
+    def test_with_db_predict(self):
+        reg = RegularGrid()
+        pod = POD()
+        db = Database(np.array([1, 2, 3])[:, None],
+                      np.array([1, 5, 3])[:, None])
+        rom = ReducedOrderModel(db, pod, reg)
+        rom.fit()
+        assert rom.predict([1]) == 1
+        assert rom.predict([2]) == 5
+        assert rom.predict([3]) == 3