From b01db5508e9285b2799db76185d945ce5d67b2c1 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Fri, 16 Jun 2023 16:28:38 +0200 Subject: [PATCH 01/25] added tests for regular grid interpolation --- tests/test_regular_grid.py | 75 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 tests/test_regular_grid.py diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py new file mode 100644 index 0000000..6904770 --- /dev/null +++ b/tests/test_regular_grid.py @@ -0,0 +1,75 @@ +import numpy as np +import warnings + +from unittest import TestCase +from ezyrb import Linear, Database, POD, ReducedOrderModel, RegularGrid + + +class TestRegularGrid(TestCase): + def test_params(self): + reg = RegularGrid(fill_value=0) + assert reg.fill_value == 0 + + def test_default_params(self): + reg = RegularGrid() + assert np.isnan(reg.fill_value) + + # def test_predict(self): + # reg = Linear() + # reg.fit([[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[1, 0], [20, 5], [8, 6]]) + # result = reg.predict([[1,2], [8,9], [6,7]]) + # assert (result[0] == [1,0]).all() + # assert (result[1] == [8,6]).all() + # assert (result[2] == [20, 5]).all() + + def test_predict1d(self): + reg = RegularGrid() + points = np.array([[1], [6], [8]]) + V = np.array([[1, 0], [20, 5], [8, 6]]) # r, n = 2, 3 + reg.fit(points, 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() + reg.fit([[1], [2]], [[1, 1], [2, 2]]) + result = reg.predict([[1.5]]) + assert (result[0] == [[1.5, 1.5]]).all() + + def test_predict1d_3(self): + reg = RegularGrid() + reg.fit([1, 2], [[1, 1], [2, 2]]) + result = reg.predict([[1.5]]) + assert (result[0] == [[1.5, 1.5]]).all() + + 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_wrong1(self): + # wrong number of params + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", category=np.VisibleDeprecationWarning) + with self.assertRaises(Exception): + reg = Linear() + reg.fit([[1, 2], [6, ], [8, 9]], [[1, 0], [20, 5], [8, 6]]) + + def test_wrong2(self): + # wrong number of values + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", category=np.VisibleDeprecationWarning) + with self.assertRaises(Exception): + reg = Linear() + reg.fit([[1, 2], [6, ], [8, 9]], [[20, 5], [8, 6]]) From b3d300195ead6002a5724ead316f2245554fe0ed Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Fri, 16 Jun 2023 16:29:38 +0200 Subject: [PATCH 02/25] added regular grid interpolation module --- ezyrb/__init__.py | 3 ++- ezyrb/regular_grid.py | 59 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 ezyrb/regular_grid.py diff --git a/ezyrb/__init__.py b/ezyrb/__init__.py index 17d87f7..1c8a017 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 0000000..6584f86 --- /dev/null +++ b/ezyrb/regular_grid.py @@ -0,0 +1,59 @@ +""" +Module for higher order interpolation on regular grids +""" + +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +from .approximation import Approximation + + +class RegularGrid(Approximation): + """ + Multidimensional interpolator on regular grids. + + :param float fill_value: value used to fill in for requested points outside + of the convex hull of the input points. If not provided, then the + default is numpy.nan. + """ + + def __init__(self, fill_value=np.nan): + self.fill_value = fill_value + self.interpolator = None + + def fit(self, grid, values): + """ + Construct the interpolator given `grid` and `values`. + + :param array_like points: the coordinates of the points. + :param array_like values: the values in the points. + """ + # the first dimension is the list of parameters, the second one is + # the dimensionality of each tuple of parameters (we look for + # parameters of dimensionality one) + + # as_np_array = np.array(points) + # if not np.issubdtype(as_np_array.dtype, np.number): + # raise ValueError('Invalid format or dimension for the argument' + # '`points`.') + + # if as_np_array.shape[-1] == 1: + # as_np_array = np.squeeze(as_np_array, axis=-1) + + # if as_np_array.ndim == 1 or (as_np_array.ndim == 2 + # and as_np_array.shape[1] == 1): + # self.interpolator = interp1d(as_np_array, values, axis=0) + # else: + # self.interpolator = LinearNDInterp(points, + # values, + # fill_value=self.fill_value) + + def predict(self, new_point): + """ + Evaluate interpolator at given `new_points`. + + :param array_like new_points: the coordinates of the given points. + :return: the interpolated values. + :rtype: numpy.ndarray + """ + return self.interpolator(new_point) From fe253045ebd7d4f6e20f30a7f1374589bee6471f Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Sat, 17 Jun 2023 12:03:12 +0200 Subject: [PATCH 03/25] added some tests --- tests/test_regular_grid.py | 80 ++++++++++++++++++++++---------------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 6904770..5de150c 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -24,9 +24,10 @@ def test_default_params(self): def test_predict1d(self): reg = RegularGrid() - points = np.array([[1], [6], [8]]) - V = np.array([[1, 0], [20, 5], [8, 6]]) # r, n = 2, 3 - reg.fit(points, V) + x1 = np.array([1, 6, 8]) + V = np.array([[1, 0], [20, 5], [8, 6]]) # n, r = 3, 2 + grid = [x1, ] + reg.fit(grid, V) result = reg.predict([[1], [8], [6]]) assert (result[0] == [1, 0]).all() assert (result[1] == [8, 6]).all() @@ -34,42 +35,53 @@ def test_predict1d(self): def test_predict1d_2(self): reg = RegularGrid() - reg.fit([[1], [2]], [[1, 1], [2, 2]]) + 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_predict1d_3(self): + def test_predict2d(self): reg = RegularGrid() - reg.fit([1, 2], [[1, 1], [2, 2]]) - result = reg.predict([[1.5]]) - assert (result[0] == [[1.5, 1.5]]).all() + x1 = [1, 2, 3] + x2 = [4, 5, 6, 7] + 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([x1, x2], V, method="linear") + result = reg.predict([[1, 4], [1, 5]]) + assert (result[0] == [1, 21]).all() + assert (result[1] == [2, 22]).all() - 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) + # TODO: test kvargs? depend on scipy version.... + # TODO: rom.fit() does not work, use reduction.fit() and approximation.fit() instead. + + # def test_with_db_predict(self): + # reg = RegularGrid() + # pod = POD() + # x1 = np.array([1, 2, 3]) + # xx, _ = np.meshgrid(x1, 1, indexing="ij") + # db = Database(xx, + # 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 + # rom.fit() + # assert rom.predict([1]) == 1 + # assert rom.predict([2]) == 5 + # assert rom.predict([3]) == 3 - def test_wrong1(self): - # wrong number of params - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", category=np.VisibleDeprecationWarning) - with self.assertRaises(Exception): - reg = Linear() - reg.fit([[1, 2], [6, ], [8, 9]], [[1, 0], [20, 5], [8, 6]]) + # def test_wrong1(self): + # # wrong number of params + # with warnings.catch_warnings(): + # warnings.filterwarnings( + # "ignore", category=np.VisibleDeprecationWarning) + # with self.assertRaises(Exception): + # reg = Linear() + # reg.fit([[1, 2], [6, ], [8, 9]], [[1, 0], [20, 5], [8, 6]]) - def test_wrong2(self): - # wrong number of values - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", category=np.VisibleDeprecationWarning) - with self.assertRaises(Exception): - reg = Linear() - reg.fit([[1, 2], [6, ], [8, 9]], [[20, 5], [8, 6]]) + # def test_wrong2(self): + # # wrong number of values + # with warnings.catch_warnings(): + # warnings.filterwarnings( + # "ignore", category=np.VisibleDeprecationWarning) + # with self.assertRaises(Exception): + # reg = Linear() + # reg.fit([[1, 2], [6, ], [8, 9]], [[20, 5], [8, 6]]) From 880962b85b3623bbeda923f27a020d9f862c012f Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Sat, 17 Jun 2023 12:04:15 +0200 Subject: [PATCH 04/25] wrapper for scipy's regular grid interpolation --- ezyrb/regular_grid.py | 46 +++++++++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 17 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 6584f86..bf344b2 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -21,7 +21,7 @@ def __init__(self, fill_value=np.nan): self.fill_value = fill_value self.interpolator = None - def fit(self, grid, values): + def fit(self, grid, values, **kvargs): """ Construct the interpolator given `grid` and `values`. @@ -32,23 +32,30 @@ def fit(self, grid, values): # the dimensionality of each tuple of parameters (we look for # parameters of dimensionality one) - # as_np_array = np.array(points) - # if not np.issubdtype(as_np_array.dtype, np.number): - # raise ValueError('Invalid format or dimension for the argument' - # '`points`.') + # we have two options + # 1.: we could make an interpolator for every mode and its coefficients + # or 2.: we "interpolate" the mode number + # option 1 is cleaner, but option 2 performs better + # X = U S VT, X being shaped m, n - # if as_np_array.shape[-1] == 1: - # as_np_array = np.squeeze(as_np_array, axis=-1) + self.dim = len(grid) + vals = np.asarray(values) + r, n = vals.T.shape # vals = (S*VT).T + self.mode_nr = np.arange(r) + extended_grid = [self.mode_nr, *grid] + shape = [r, ] + for i in range(self.dim): + shape.append(len(grid[i])) + assert np.prod(shape) == vals.size, "Values don't match grid. "\ + "Make sure to pass a grid, not a list of points!\n"\ + "HINT: did you use rom.fit()? This method does not work with a "\ + "grid. Use reduction.fit(...) and approximation.fit(...) instead." + self.interpolator = RegularGridInterpolator(extended_grid, + vals.T.reshape(shape), + fill_value=self.fill_value, + **kvargs) - # if as_np_array.ndim == 1 or (as_np_array.ndim == 2 - # and as_np_array.shape[1] == 1): - # self.interpolator = interp1d(as_np_array, values, axis=0) - # else: - # self.interpolator = LinearNDInterp(points, - # values, - # fill_value=self.fill_value) - - def predict(self, new_point): + def predict(self, new_points): """ Evaluate interpolator at given `new_points`. @@ -56,4 +63,9 @@ def predict(self, new_point): :return: the interpolated values. :rtype: numpy.ndarray """ - return self.interpolator(new_point) + dim = self.dim + xi_extended = np.zeros((len(self.mode_nr), len(new_points), dim+1)) + xi_extended[:, :, 0] = self.mode_nr[:, None] + for i in range(dim): + xi_extended[:, :, i+1] = np.array(new_points)[:, i] + return self.interpolator(xi_extended).T From ce6b8e4b512ae9e187af7bdb0a1ed51360ff0c77 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Sat, 17 Jun 2023 12:39:09 +0200 Subject: [PATCH 05/25] removed some comments --- ezyrb/regular_grid.py | 8 ++----- tests/test_regular_grid.py | 44 +------------------------------------- 2 files changed, 3 insertions(+), 49 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index bf344b2..d90d194 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -1,10 +1,8 @@ """ Module for higher order interpolation on regular grids """ - import numpy as np from scipy.interpolate import RegularGridInterpolator - from .approximation import Approximation @@ -20,6 +18,7 @@ class RegularGrid(Approximation): def __init__(self, fill_value=np.nan): self.fill_value = fill_value self.interpolator = None + return def fit(self, grid, values, **kvargs): """ @@ -28,10 +27,6 @@ def fit(self, grid, values, **kvargs): :param array_like points: the coordinates of the points. :param array_like values: the values in the points. """ - # the first dimension is the list of parameters, the second one is - # the dimensionality of each tuple of parameters (we look for - # parameters of dimensionality one) - # we have two options # 1.: we could make an interpolator for every mode and its coefficients # or 2.: we "interpolate" the mode number @@ -54,6 +49,7 @@ def fit(self, grid, values, **kvargs): vals.T.reshape(shape), fill_value=self.fill_value, **kvargs) + return def predict(self, new_points): """ diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 5de150c..8686d98 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -1,8 +1,6 @@ import numpy as np -import warnings - from unittest import TestCase -from ezyrb import Linear, Database, POD, ReducedOrderModel, RegularGrid +from ezyrb import RegularGrid # Database, POD, ReducedOrderModel class TestRegularGrid(TestCase): @@ -14,14 +12,6 @@ def test_default_params(self): reg = RegularGrid() assert np.isnan(reg.fill_value) - # def test_predict(self): - # reg = Linear() - # reg.fit([[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[1, 0], [20, 5], [8, 6]]) - # result = reg.predict([[1,2], [8,9], [6,7]]) - # assert (result[0] == [1,0]).all() - # assert (result[1] == [8,6]).all() - # assert (result[2] == [20, 5]).all() - def test_predict1d(self): reg = RegularGrid() x1 = np.array([1, 6, 8]) @@ -53,35 +43,3 @@ def test_predict2d(self): # TODO: test kvargs? depend on scipy version.... # TODO: rom.fit() does not work, use reduction.fit() and approximation.fit() instead. - - # def test_with_db_predict(self): - # reg = RegularGrid() - # pod = POD() - # x1 = np.array([1, 2, 3]) - # xx, _ = np.meshgrid(x1, 1, indexing="ij") - # db = Database(xx, - # 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_wrong1(self): - # # wrong number of params - # with warnings.catch_warnings(): - # warnings.filterwarnings( - # "ignore", category=np.VisibleDeprecationWarning) - # with self.assertRaises(Exception): - # reg = Linear() - # reg.fit([[1, 2], [6, ], [8, 9]], [[1, 0], [20, 5], [8, 6]]) - - # def test_wrong2(self): - # # wrong number of values - # with warnings.catch_warnings(): - # warnings.filterwarnings( - # "ignore", category=np.VisibleDeprecationWarning) - # with self.assertRaises(Exception): - # reg = Linear() - # reg.fit([[1, 2], [6, ], [8, 9]], [[20, 5], [8, 6]]) From 8df746976d5ed47fb23616e1ee4857b8e648a4ee Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Sat, 17 Jun 2023 17:27:08 +0200 Subject: [PATCH 06/25] minor code style improvements --- ezyrb/regular_grid.py | 49 +++++++++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index d90d194..4b6c77b 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -9,6 +9,7 @@ class RegularGrid(Approximation): """ Multidimensional interpolator on regular grids. + See scipy's RegularGridInterpolator :param float fill_value: value used to fill in for requested points outside of the convex hull of the input points. If not provided, then the @@ -18,14 +19,30 @@ class RegularGrid(Approximation): def __init__(self, fill_value=np.nan): self.fill_value = fill_value self.interpolator = None - return + self.dim = None + self.mode_nr = None - def fit(self, grid, values, **kvargs): + def fit(self, points, values, **kvargs): """ - Construct the interpolator given `grid` and `values`. + Construct the interpolator given `points` and `values`. + + Parameters + ---------- + points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, ) + The points defining the regular grid in n dimensions. The points in + each dimension (i.e. every elements of the points tuple) must be + strictly ascending or descending. + + values : array_like, shape (m1, ..., mn, ...) + The data on the regular grid in n dimensions. Complex data can be + acceptable. + **kvargs : TYPE + see scipy.interpolate.RegularGridInterpolator + + Returns + ------- + None. - :param array_like points: the coordinates of the points. - :param array_like values: the values in the points. """ # we have two options # 1.: we could make an interpolator for every mode and its coefficients @@ -33,14 +50,14 @@ def fit(self, grid, values, **kvargs): # option 1 is cleaner, but option 2 performs better # X = U S VT, X being shaped m, n - self.dim = len(grid) + self.dim = len(points) vals = np.asarray(values) - r, n = vals.T.shape # vals = (S*VT).T + r = vals.T.shape[0] # vals = (S*VT).T self.mode_nr = np.arange(r) - extended_grid = [self.mode_nr, *grid] + extended_grid = [self.mode_nr, *points] shape = [r, ] for i in range(self.dim): - shape.append(len(grid[i])) + shape.append(len(points[i])) assert np.prod(shape) == vals.size, "Values don't match grid. "\ "Make sure to pass a grid, not a list of points!\n"\ "HINT: did you use rom.fit()? This method does not work with a "\ @@ -49,19 +66,11 @@ def fit(self, grid, values, **kvargs): vals.T.reshape(shape), fill_value=self.fill_value, **kvargs) - return - def predict(self, new_points): - """ - Evaluate interpolator at given `new_points`. - - :param array_like new_points: the coordinates of the given points. - :return: the interpolated values. - :rtype: numpy.ndarray - """ + def predict(self, new_point): dim = self.dim - xi_extended = np.zeros((len(self.mode_nr), len(new_points), dim+1)) + xi_extended = np.zeros((len(self.mode_nr), len(new_point), dim+1)) xi_extended[:, :, 0] = self.mode_nr[:, None] for i in range(dim): - xi_extended[:, :, i+1] = np.array(new_points)[:, i] + xi_extended[:, :, i+1] = np.array(new_point)[:, i] return self.interpolator(xi_extended).T From 6977a9e289af941b63419a6c8b49d3bf173aecb4 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Sun, 18 Jun 2023 18:17:11 +0200 Subject: [PATCH 07/25] adjusted docsting --- ezyrb/regular_grid.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 4b6c77b..a0448b5 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -25,30 +25,21 @@ def __init__(self, fill_value=np.nan): def fit(self, points, values, **kvargs): """ Construct the interpolator given `points` and `values`. + see scipy.interpolate.RegularGridInterpolator - Parameters - ---------- - points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, ) + :param array_like points: with shapes (m1, ), ..., (mn, ) The points defining the regular grid in n dimensions. The points in each dimension (i.e. every elements of the points tuple) must be strictly ascending or descending. - values : array_like, shape (m1, ..., mn, ...) - The data on the regular grid in n dimensions. Complex data can be - acceptable. - **kvargs : TYPE - see scipy.interpolate.RegularGridInterpolator - - Returns - ------- - None. - + :param array_like values: shape (m1, ..., mn, ...) + The data on the regular grid in n dimensions. """ # we have two options # 1.: we could make an interpolator for every mode and its coefficients # or 2.: we "interpolate" the mode number # option 1 is cleaner, but option 2 performs better - # X = U S VT, X being shaped m, n + # X = U S VT, X being shaped (m, n) self.dim = len(points) vals = np.asarray(values) From 24bee31614cb3932f09baf3548cddfeab190c234 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 19 Jun 2023 10:55:43 +0200 Subject: [PATCH 08/25] added pre-processing to map points and values onto regular grid in 2D --- ezyrb/regular_grid.py | 10 ++++++++++ tests/test_regular_grid.py | 18 ++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index a0448b5..8d48891 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -22,6 +22,16 @@ def __init__(self, fill_value=np.nan): self.dim = None self.mode_nr = None + def get_grid_axes(self, pts_scrmbld, vals_scrmbld): + # rounding errors! + # works in 2D + x1, unique_inverse1 = np.unique(pts_scrmbld[:, 0], return_inverse=True) + x2, unique_inverse2 = np.unique(pts_scrmbld[:, 1], return_inverse=True) + new_row_index = unique_inverse1*len(x2)+unique_inverse2 + reverse_scrambling = np.argsort(new_row_index) + vals_on_regular_grid = vals_scrmbld[reverse_scrambling, :] + return (x1, x2), vals_on_regular_grid + def fit(self, points, values, **kvargs): """ Construct the interpolator given `points` and `values`. diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 8686d98..84bb55d 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -41,5 +41,23 @@ def test_predict2d(self): assert (result[0] == [1, 21]).all() assert (result[1] == [2, 22]).all() + def test_points2grid_2D(): + 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() + + (x1, x2), V_unscrambled = reg.get_grid_axes(pts_scrmbld, V_scrmbld) + + assert np.allclose(V_unscrambled, V) + # TODO: test kvargs? depend on scipy version.... # TODO: rom.fit() does not work, use reduction.fit() and approximation.fit() instead. From b5835b8119754df3848aa6bd607c63384aa02d9b Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 19 Jun 2023 11:46:22 +0200 Subject: [PATCH 09/25] modified pre-processing to map points and values onto regular grid in nD --- ezyrb/regular_grid.py | 46 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 6 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 8d48891..55205e6 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -24,13 +24,23 @@ def __init__(self, fill_value=np.nan): def get_grid_axes(self, pts_scrmbld, vals_scrmbld): # rounding errors! - # works in 2D - x1, unique_inverse1 = np.unique(pts_scrmbld[:, 0], return_inverse=True) - x2, unique_inverse2 = np.unique(pts_scrmbld[:, 1], return_inverse=True) - new_row_index = unique_inverse1*len(x2)+unique_inverse2 + 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): + raise ValueError("points and values are not on a regular grid") + new_row_index = calculate_flat_index(iN, nN) reverse_scrambling = np.argsort(new_row_index) vals_on_regular_grid = vals_scrmbld[reverse_scrambling, :] - return (x1, x2), vals_on_regular_grid + return grid_axes, vals_on_regular_grid def fit(self, points, values, **kvargs): """ @@ -41,7 +51,6 @@ def fit(self, points, values, **kvargs): The points defining the regular grid in n dimensions. The points in each dimension (i.e. every elements of the points tuple) must be strictly ascending or descending. - :param array_like values: shape (m1, ..., mn, ...) The data on the regular grid in n dimensions. """ @@ -75,3 +84,28 @@ def predict(self, new_point): for i in range(dim): xi_extended[:, :, i+1] = np.array(new_point)[:, i] return self.interpolator(xi_extended).T + + +def calculate_flat_index(iN, nN): + """ + Calculates the flat index for a multidimensional array given the indices + and dimensions. + + :param list iN: indices representing the position of the element(s) + in each dimension. + :param list nN: size of the array in each dimension. + + :rtype: numpy.ndarray + """ + # index = i1 + n1 * (i2 + n2 * (... (iN-1 + nN-1 * iN) ...)) + if len(iN) != len(nN): + raise ValueError("The lengths of iN and nN should be the same.") + + if any((i < 0).any() or (i >= n).any() for i, n in zip(iN, nN)): + raise ValueError("The indices are out of bounds.") + + index = 0 + for i, n in zip(iN, nN): + index = i + n * index + + return index From a9deaa0228e146b5d173ffe78d02d181975c4617 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 19 Jun 2023 12:53:10 +0200 Subject: [PATCH 10/25] modified the regular grid interpolation so it works with the grid points, not the grid axes. --- ezyrb/regular_grid.py | 63 ++++++++++++++++++++++++++------------ tests/test_regular_grid.py | 57 +++++++++++++++++++++++++++++----- 2 files changed, 93 insertions(+), 27 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 55205e6..fe1b7c6 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -20,10 +20,18 @@ def __init__(self, fill_value=np.nan): self.fill_value = fill_value self.interpolator = None self.dim = None + self.n_modes = 0 self.mode_nr = None def get_grid_axes(self, pts_scrmbld, vals_scrmbld): - # rounding errors! + """ + calculates the grid axes from a meshed grid. The grid axes are given as + a tuple of ndarray of float, with shapes (m1, ), …, (mn, ) + The values are ordered so they fit on a mesh generated with + numpy.meshgrid(ax1, ax2, ..., axn, indexing="ij") + + """ + # be aware of floating point precision in points! grid_axes = [] iN = [] # index in dimension N nN = [] # size of dimension N @@ -45,44 +53,61 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): def fit(self, points, values, **kvargs): """ 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: with shapes (m1, ), ..., (mn, ) - The points defining the regular grid in n dimensions. The points in - each dimension (i.e. every elements of the points tuple) must be - strictly ascending or descending. - :param array_like values: shape (m1, ..., mn, ...) - The data on the regular grid in n dimensions. + :param array_like points: the coordinates of the points. + :param array_like values: the values in the points. """ # we have two options # 1.: we could make an interpolator for every mode and its coefficients # or 2.: we "interpolate" the mode number # option 1 is cleaner, but option 2 performs better # X = U S VT, X being shaped (m, n) + points = np.array(points) + if not np.issubdtype(points.dtype, np.number): + raise ValueError('Invalid format or dimension for the argument' + '`points`.') - self.dim = len(points) + if len(points.shape) == 1: + points.shape = (-1, 1) + + self.dim = len(points[0]) vals = np.asarray(values) - r = vals.T.shape[0] # vals = (S*VT).T - self.mode_nr = np.arange(r) - extended_grid = [self.mode_nr, *points] - shape = [r, ] + grid_axes, values_grd = self.get_grid_axes(points, vals) + self.n_modes = vals.T.shape[0] # vals = (S@VT).T = S@V + if self.n_modes > 1: + self.mode_nr = np.arange(self.n_modes) + extended_grid = [self.mode_nr, *grid_axes] + shape = [self.n_modes, ] + else: + extended_grid = grid_axes + shape = [] for i in range(self.dim): - shape.append(len(points[i])) - assert np.prod(shape) == vals.size, "Values don't match grid. "\ + shape.append(len(grid_axes[i])) + assert np.prod(shape) == values_grd.size, "Values don't match grid. "\ "Make sure to pass a grid, not a list of points!\n"\ "HINT: did you use rom.fit()? This method does not work with a "\ "grid. Use reduction.fit(...) and approximation.fit(...) instead." self.interpolator = RegularGridInterpolator(extended_grid, - vals.T.reshape(shape), + values_grd.T.reshape( + shape), fill_value=self.fill_value, **kvargs) def predict(self, new_point): + new_point = np.array(new_point) + if len(new_point.shape) == 1: + new_point.shape = (-1, 1) + dim = self.dim - xi_extended = np.zeros((len(self.mode_nr), len(new_point), dim+1)) - xi_extended[:, :, 0] = self.mode_nr[:, None] - for i in range(dim): - xi_extended[:, :, i+1] = np.array(new_point)[:, i] + if self.n_modes > 1: + xi_extended = np.zeros((len(self.mode_nr), len(new_point), dim+1)) + xi_extended[:, :, 0] = self.mode_nr[:, None] + for i in range(dim): + xi_extended[:, :, i+1] = np.array(new_point)[:, i] + else: + xi_extended = new_point return self.interpolator(xi_extended).T diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 84bb55d..9ddb98d 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -1,6 +1,6 @@ import numpy as np from unittest import TestCase -from ezyrb import RegularGrid # Database, POD, ReducedOrderModel +from ezyrb import RegularGrid, Database, POD, ReducedOrderModel class TestRegularGrid(TestCase): @@ -12,12 +12,19 @@ def test_default_params(self): reg = RegularGrid() assert np.isnan(reg.fill_value) + 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]) + x1 = np.array([[1], [6], [8]]) V = np.array([[1, 0], [20, 5], [8, 6]]) # n, r = 3, 2 - grid = [x1, ] - reg.fit(grid, V) + reg.fit(x1, V) result = reg.predict([[1], [8], [6]]) assert (result[0] == [1, 0]).all() assert (result[1] == [8, 6]).all() @@ -26,7 +33,7 @@ def test_predict1d(self): def test_predict1d_2(self): reg = RegularGrid() x1 = [1, 2] - reg.fit([x1, ], [[1, 1], [2, 2]]) + reg.fit(x1, [[1, 1], [2, 2]]) result = reg.predict([[1.5]]) assert (result[0] == [[1.5, 1.5]]).all() @@ -34,14 +41,16 @@ 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([x1, x2], V, method="linear") + 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_points2grid_2D(): + def test_get_grid_axes_2D(self): x1 = [.1, .2, .3] x2 = [4, 5, 6, 7] xx1, xx2 = np.meshgrid(x1, x2, indexing="ij") @@ -54,10 +63,42 @@ def test_points2grid_2D(): 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) - (x1, x2), V_unscrambled = reg.get_grid_axes(pts_scrmbld, V_scrmbld) + 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 # TODO: test kvargs? depend on scipy version.... # TODO: rom.fit() does not work, use reduction.fit() and approximation.fit() instead. From 98d1425b097465b7d2bcc2effbf0dafe85b07779 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 19 Jun 2023 12:58:34 +0200 Subject: [PATCH 11/25] removed useless fill_value argument --- ezyrb/regular_grid.py | 7 +------ tests/test_regular_grid.py | 8 -------- 2 files changed, 1 insertion(+), 14 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index fe1b7c6..55e0e25 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -11,13 +11,9 @@ class RegularGrid(Approximation): Multidimensional interpolator on regular grids. See scipy's RegularGridInterpolator - :param float fill_value: value used to fill in for requested points outside - of the convex hull of the input points. If not provided, then the - default is numpy.nan. """ - def __init__(self, fill_value=np.nan): - self.fill_value = fill_value + def __init__(self): self.interpolator = None self.dim = None self.n_modes = 0 @@ -92,7 +88,6 @@ def fit(self, points, values, **kvargs): self.interpolator = RegularGridInterpolator(extended_grid, values_grd.T.reshape( shape), - fill_value=self.fill_value, **kvargs) def predict(self, new_point): diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 9ddb98d..8070766 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -4,13 +4,6 @@ class TestRegularGrid(TestCase): - def test_params(self): - reg = RegularGrid(fill_value=0) - assert reg.fill_value == 0 - - def test_default_params(self): - reg = RegularGrid() - assert np.isnan(reg.fill_value) def test_1D_1mode(self): reg = RegularGrid() @@ -101,4 +94,3 @@ def test_with_db_predict(self): assert rom.predict([2]) == 5 assert rom.predict([3]) == 3 # TODO: test kvargs? depend on scipy version.... - # TODO: rom.fit() does not work, use reduction.fit() and approximation.fit() instead. From 761fcde1cc72c7886707e4c912d03218f369850c Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 10:21:16 +0200 Subject: [PATCH 12/25] added a doc file for the automatic doc generation of the new module --- docs/source/regular_grid.rst | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 docs/source/regular_grid.rst diff --git a/docs/source/regular_grid.rst b/docs/source/regular_grid.rst new file mode 100644 index 0000000..8c627ac --- /dev/null +++ b/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: From 9326752da3068fa2cd9f2d625d95b222469233df Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:15:30 +0200 Subject: [PATCH 13/25] improved docstring of get_grid_axes --- ezyrb/regular_grid.py | 54 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 55e0e25..d95ce95 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -4,13 +4,47 @@ import numpy as np from scipy.interpolate import RegularGridInterpolator from .approximation import Approximation +from datetime import datetime as dt class RegularGrid(Approximation): """ Multidimensional interpolator on regular grids. - See scipy's RegularGridInterpolator - + 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): @@ -21,11 +55,16 @@ def __init__(self): def get_grid_axes(self, pts_scrmbld, vals_scrmbld): """ - calculates the grid axes from a meshed grid. The grid axes are given as - a tuple of ndarray of float, with shapes (m1, ), …, (mn, ) - The values are ordered so they fit on a mesh generated with + 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) + """ # be aware of floating point precision in points! grid_axes = [] @@ -38,7 +77,6 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): grid_axes.append(xn) nN.append(len(xn)) iN.append(unique_inverse_n) - if np.prod(nN) != len(vals_scrmbld): raise ValueError("points and values are not on a regular grid") new_row_index = calculate_flat_index(iN, nN) @@ -46,7 +84,7 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): vals_on_regular_grid = vals_scrmbld[reverse_scrambling, :] return grid_axes, vals_on_regular_grid - def fit(self, points, values, **kvargs): + 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. @@ -88,7 +126,7 @@ def fit(self, points, values, **kvargs): self.interpolator = RegularGridInterpolator(extended_grid, values_grd.T.reshape( shape), - **kvargs) + **kwargs) def predict(self, new_point): new_point = np.array(new_point) From c2692fb44f5290bd3f245620d6444d40ec748af5 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:28:19 +0200 Subject: [PATCH 14/25] moved ValueError to the right place and adjusted the error message. --- ezyrb/regular_grid.py | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index d95ce95..9c4d03c 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -2,9 +2,8 @@ Module for higher order interpolation on regular grids """ import numpy as np -from scipy.interpolate import RegularGridInterpolator +from scipy.interpolate import RegularGridInterpolator as RGI from .approximation import Approximation -from datetime import datetime as dt class RegularGrid(Approximation): @@ -66,7 +65,6 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): :rtype: (list, numpy.ndarray) """ - # be aware of floating point precision in points! grid_axes = [] iN = [] # index in dimension N nN = [] # size of dimension N @@ -78,7 +76,10 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): nN.append(len(xn)) iN.append(unique_inverse_n) if np.prod(nN) != len(vals_scrmbld): - raise ValueError("points and values are not on a regular grid") + 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 = calculate_flat_index(iN, nN) reverse_scrambling = np.argsort(new_row_index) vals_on_regular_grid = vals_scrmbld[reverse_scrambling, :] @@ -93,11 +94,6 @@ def fit(self, points, values, **kwargs): :param array_like points: the coordinates of the points. :param array_like values: the values in the points. """ - # we have two options - # 1.: we could make an interpolator for every mode and its coefficients - # or 2.: we "interpolate" the mode number - # option 1 is cleaner, but option 2 performs better - # X = U S VT, X being shaped (m, n) points = np.array(points) if not np.issubdtype(points.dtype, np.number): raise ValueError('Invalid format or dimension for the argument' @@ -109,7 +105,7 @@ def fit(self, points, values, **kwargs): self.dim = len(points[0]) vals = np.asarray(values) grid_axes, values_grd = self.get_grid_axes(points, vals) - self.n_modes = vals.T.shape[0] # vals = (S@VT).T = S@V + self.n_modes = vals.T.shape[0] if self.n_modes > 1: self.mode_nr = np.arange(self.n_modes) extended_grid = [self.mode_nr, *grid_axes] @@ -119,14 +115,8 @@ def fit(self, points, values, **kwargs): shape = [] for i in range(self.dim): shape.append(len(grid_axes[i])) - assert np.prod(shape) == values_grd.size, "Values don't match grid. "\ - "Make sure to pass a grid, not a list of points!\n"\ - "HINT: did you use rom.fit()? This method does not work with a "\ - "grid. Use reduction.fit(...) and approximation.fit(...) instead." - self.interpolator = RegularGridInterpolator(extended_grid, - values_grd.T.reshape( - shape), - **kwargs) + self.interpolator = RGI(extended_grid, values_grd.T.reshape(shape), + **kwargs) def predict(self, new_point): new_point = np.array(new_point) From 3f790630603dc4a71e11a4af05eeaf918de07bcf Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:32:15 +0200 Subject: [PATCH 15/25] added docstring for predict method --- ezyrb/regular_grid.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 9c4d03c..b2398c2 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -119,6 +119,13 @@ def fit(self, points, values, **kwargs): **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.array(new_point) if len(new_point.shape) == 1: new_point.shape = (-1, 1) From eb7f1a30e99a5f0751d0d84e9d50130b6b5db906 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:33:27 +0200 Subject: [PATCH 16/25] removed 'dim' variable --- ezyrb/regular_grid.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index b2398c2..f795246 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -130,11 +130,11 @@ def predict(self, new_point): if len(new_point.shape) == 1: new_point.shape = (-1, 1) - dim = self.dim if self.n_modes > 1: - xi_extended = np.zeros((len(self.mode_nr), len(new_point), dim+1)) + shape = (len(self.mode_nr), len(new_point), self.dim+1) + xi_extended = np.zeros(shape) xi_extended[:, :, 0] = self.mode_nr[:, None] - for i in range(dim): + for i in range(self.dim): xi_extended[:, :, i+1] = np.array(new_point)[:, i] else: xi_extended = new_point From 706bae70524c917d60469cad2fc454a45092452a Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:37:25 +0200 Subject: [PATCH 17/25] removed calculate_flat_index(...) entirely --- ezyrb/regular_grid.py | 27 +-------------------------- 1 file changed, 1 insertion(+), 26 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index f795246..3400e42 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -80,7 +80,7 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): "points that are on a regular grid! Be aware of floating "\ "point precision." raise ValueError(msg) - new_row_index = calculate_flat_index(iN, nN) + 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 @@ -139,28 +139,3 @@ def predict(self, new_point): else: xi_extended = new_point return self.interpolator(xi_extended).T - - -def calculate_flat_index(iN, nN): - """ - Calculates the flat index for a multidimensional array given the indices - and dimensions. - - :param list iN: indices representing the position of the element(s) - in each dimension. - :param list nN: size of the array in each dimension. - - :rtype: numpy.ndarray - """ - # index = i1 + n1 * (i2 + n2 * (... (iN-1 + nN-1 * iN) ...)) - if len(iN) != len(nN): - raise ValueError("The lengths of iN and nN should be the same.") - - if any((i < 0).any() or (i >= n).any() for i, n in zip(iN, nN)): - raise ValueError("The indices are out of bounds.") - - index = 0 - for i, n in zip(iN, nN): - index = i + n * index - - return index From 4148b4f49919206088b0789ccf5fbb0eb8f57aa6 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:40:54 +0200 Subject: [PATCH 18/25] removed TODO comment --- tests/test_regular_grid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 8070766..487f895 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -93,4 +93,3 @@ def test_with_db_predict(self): assert rom.predict([1]) == 1 assert rom.predict([2]) == 5 assert rom.predict([3]) == 3 - # TODO: test kvargs? depend on scipy version.... From 4a770305d2967a99e9625aed628b10fa871a6bc0 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 11:47:57 +0200 Subject: [PATCH 19/25] other cosmetic changes --- ezyrb/regular_grid.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 3400e42..9d3e5b7 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -82,7 +82,7 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): 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, :] + vals_on_regular_grid = vals_scrmbld[reverse_scrambling] return grid_axes, vals_on_regular_grid def fit(self, points, values, **kwargs): @@ -94,18 +94,16 @@ def fit(self, points, values, **kwargs): :param array_like points: the coordinates of the points. :param array_like values: the values in the points. """ - points = np.array(points) + points = np.asarray(points) if not np.issubdtype(points.dtype, np.number): raise ValueError('Invalid format or dimension for the argument' '`points`.') + points = np.atleast_2d(points) - if len(points.shape) == 1: - points.shape = (-1, 1) - - self.dim = len(points[0]) + self.dim = points.shape[1] vals = np.asarray(values) grid_axes, values_grd = self.get_grid_axes(points, vals) - self.n_modes = vals.T.shape[0] + self.n_modes = vals.shape[-1] if self.n_modes > 1: self.mode_nr = np.arange(self.n_modes) extended_grid = [self.mode_nr, *grid_axes] @@ -126,16 +124,15 @@ def predict(self, new_point): :return: the interpolated values. :rtype: numpy.ndarray """ - new_point = np.array(new_point) + new_point = np.asarray(new_point) if len(new_point.shape) == 1: new_point.shape = (-1, 1) if self.n_modes > 1: shape = (len(self.mode_nr), len(new_point), self.dim+1) xi_extended = np.zeros(shape) - xi_extended[:, :, 0] = self.mode_nr[:, None] - for i in range(self.dim): - xi_extended[:, :, i+1] = np.array(new_point)[:, i] + xi_extended[..., 0] = self.mode_nr[:, None] + xi_extended[..., 1:] = new_point else: xi_extended = new_point return self.interpolator(xi_extended).T From 638af8d8c7995f1d7066b9a4df2c89ab33bcc76c Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 20:03:59 +0200 Subject: [PATCH 20/25] some simplifications, no need to operate with the transposed of the values. --- ezyrb/regular_grid.py | 31 +++++-------------------------- 1 file changed, 5 insertions(+), 26 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 9d3e5b7..f45468f 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -48,9 +48,6 @@ class RegularGrid(Approximation): def __init__(self): self.interpolator = None - self.dim = None - self.n_modes = 0 - self.mode_nr = None def get_grid_axes(self, pts_scrmbld, vals_scrmbld): """ @@ -76,6 +73,7 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): nN.append(len(xn)) iN.append(unique_inverse_n) if np.prod(nN) != len(vals_scrmbld): + print(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." @@ -99,22 +97,11 @@ def fit(self, points, values, **kwargs): raise ValueError('Invalid format or dimension for the argument' '`points`.') points = np.atleast_2d(points) - - self.dim = points.shape[1] vals = np.asarray(values) grid_axes, values_grd = self.get_grid_axes(points, vals) - self.n_modes = vals.shape[-1] - if self.n_modes > 1: - self.mode_nr = np.arange(self.n_modes) - extended_grid = [self.mode_nr, *grid_axes] - shape = [self.n_modes, ] - else: - extended_grid = grid_axes - shape = [] - for i in range(self.dim): - shape.append(len(grid_axes[i])) - self.interpolator = RGI(extended_grid, values_grd.T.reshape(shape), - **kwargs) + 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): """ @@ -127,12 +114,4 @@ def predict(self, new_point): new_point = np.asarray(new_point) if len(new_point.shape) == 1: new_point.shape = (-1, 1) - - if self.n_modes > 1: - shape = (len(self.mode_nr), len(new_point), self.dim+1) - xi_extended = np.zeros(shape) - xi_extended[..., 0] = self.mode_nr[:, None] - xi_extended[..., 1:] = new_point - else: - xi_extended = new_point - return self.interpolator(xi_extended).T + return self.interpolator(new_point) From e49ad1883cb870d757bb469cd8a92722ac23cd97 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Mon, 26 Jun 2023 20:08:51 +0200 Subject: [PATCH 21/25] reverted suggested change, as it changed the behaviour --- ezyrb/regular_grid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index f45468f..2a8a7f2 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -73,7 +73,6 @@ def get_grid_axes(self, pts_scrmbld, vals_scrmbld): nN.append(len(xn)) iN.append(unique_inverse_n) if np.prod(nN) != len(vals_scrmbld): - print(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." @@ -96,7 +95,8 @@ def fit(self, points, values, **kwargs): if not np.issubdtype(points.dtype, np.number): raise ValueError('Invalid format or dimension for the argument' '`points`.') - points = np.atleast_2d(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] From 3f77b6a1ee8fd8828624e0d91f707a554e381c64 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Tue, 27 Jun 2023 11:56:52 +0200 Subject: [PATCH 22/25] improved code style --- ezyrb/regular_grid.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ezyrb/regular_grid.py b/ezyrb/regular_grid.py index 2a8a7f2..01f2870 100644 --- a/ezyrb/regular_grid.py +++ b/ezyrb/regular_grid.py @@ -95,8 +95,8 @@ def fit(self, points, values, **kwargs): 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) + 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] @@ -112,6 +112,6 @@ def predict(self, new_point): :rtype: numpy.ndarray """ new_point = np.asarray(new_point) - if len(new_point.shape) == 1: - new_point.shape = (-1, 1) + if new_point.ndim == 1: + new_point = new_point[:, None] return self.interpolator(new_point) From 48f0c424cf4644964c4de31bdc80bb6608dc00c9 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Tue, 27 Jun 2023 13:07:37 +0200 Subject: [PATCH 23/25] updated autogenerated rst and added it to the code documentation --- docs/source/code.rst | 1 + docs/source/regular_grid.rst | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/source/code.rst b/docs/source/code.rst index 635e01f..844f9fd 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 index 8c627ac..bb7cdf9 100644 --- a/docs/source/regular_grid.rst +++ b/docs/source/regular_grid.rst @@ -10,10 +10,10 @@ RegularGrid :nosignatures: RegularGrid + RegularGrid.def RegularGrid.get_grid_axes RegularGrid.fit RegularGrid.predict - RegularGrid.calculate_flat_index .. autoclass:: RegularGrid :members: From 1d8b6902b3cf45ded076cee88ad420a50c136f49 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Wed, 28 Jun 2023 11:24:14 +0200 Subject: [PATCH 24/25] added test to make sure exception is raised. --- tests/test_regular_grid.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 487f895..a5540d0 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -1,5 +1,5 @@ import numpy as np -from unittest import TestCase +from unittest import TestCase, main from ezyrb import RegularGrid, Database, POD, ReducedOrderModel @@ -93,3 +93,15 @@ def test_with_db_predict(self): 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() From 2ad23aee374656fa0a9afb0b788aa5e04f5819c8 Mon Sep 17 00:00:00 2001 From: Florian Arbes Date: Wed, 28 Jun 2023 12:01:57 +0200 Subject: [PATCH 25/25] removed 'def' member manually --- docs/source/regular_grid.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/source/regular_grid.rst b/docs/source/regular_grid.rst index bb7cdf9..2b1901e 100644 --- a/docs/source/regular_grid.rst +++ b/docs/source/regular_grid.rst @@ -10,7 +10,6 @@ RegularGrid :nosignatures: RegularGrid - RegularGrid.def RegularGrid.get_grid_axes RegularGrid.fit RegularGrid.predict