diff --git a/docs/source/code.rst b/docs/source/code.rst index 635e01fe..844f9fd2 100644 --- a/docs/source/code.rst +++ b/docs/source/code.rst @@ -17,3 +17,4 @@ Code Documentation ae podae reducedordermodel + regular_grid diff --git a/docs/source/regular_grid.rst b/docs/source/regular_grid.rst new file mode 100644 index 00000000..2b1901ec --- /dev/null +++ b/docs/source/regular_grid.rst @@ -0,0 +1,22 @@ +RegularGrid +===================== + +.. currentmodule:: ezyrb.regular_grid + +.. automodule:: ezyrb.regular_grid + +.. autosummary:: + :toctree: _summaries + :nosignatures: + + RegularGrid + RegularGrid.get_grid_axes + RegularGrid.fit + RegularGrid.predict + +.. autoclass:: RegularGrid + :members: + :private-members: + :undoc-members: + :show-inheritance: + :noindex: diff --git a/ezyrb/__init__.py b/ezyrb/__init__.py index 17d87f7a..1c8a017e 100644 --- a/ezyrb/__init__.py +++ b/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 diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py new file mode 100644 index 00000000..01f28700 --- /dev/null +++ b/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 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) + + 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) diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py new file mode 100644 index 00000000..a5540d07 --- /dev/null +++ b/tests/test_regular_grid.py @@ -0,0 +1,107 @@ +import numpy as np +from unittest import TestCase, main +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 + + def test_fails(self): + reg = RegularGrid() + reg = RegularGrid() + p = [[1, 2]] + V = [[1, 1], [2, 2]] + with self.assertRaises(ValueError): + reg.fit(p, V) + + +if __name__ == "__main__": + main()