From c57c20072a6000605e549964becee03457541f26 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:09:26 -0300 Subject: [PATCH 1/8] ENH: add multiple argument support for shepard opt. - Minor bug in computation if distance zero solved --- rocketpy/mathutils/function.py | 80 +++++++++++++++++++--------------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 251d58199..8d97565ef 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -482,25 +482,33 @@ def get_value_opt(x): return y elif self.__interpolation__ == "shepard": - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - len_y_data = len(y_data) # A little speed up - # change the function's name to avoid mypy's error def get_value_opt_multiple(*args): - x = np.array([[float(val) for val in args]]) + x_data = self.source[:, 0:-1] # Support for N-Dimensions + y_data = self.source[:, -1] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=1) + distances_squared = np.sum(sub_matrix**2, axis=2) - zero_distance_index = np.where(distances_squared == 0)[0] - if len(zero_distance_index) > 0: - return y_data[zero_distance_index[0]] + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False - weights = distances_squared ** (-1.5) - numerator_sum = np.sum(y_data * weights) - denominator_sum = np.sum(weights) + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] - return numerator_sum / denominator_sum + return result get_value_opt = get_value_opt_multiple @@ -872,28 +880,32 @@ def get_value(self, *args): # Returns value for shepard interpolation elif self.__interpolation__ == "shepard": - if all(isinstance(arg, Iterable) for arg in args): - x = list(np.column_stack(args)) - else: - x = [[float(x) for x in list(args)]] - ans = x - x_data = self.source[:, 0:-1] + x_data = self.source[:, 0:-1] # Support for N-Dimensions y_data = self.source[:, -1] - for i, _ in enumerate(x): - numerator_sum = 0 - denominator_sum = 0 - for o, _ in enumerate(y_data): - sub = x_data[o] - x[i] - distance = (sub.dot(sub)) ** (0.5) - if distance == 0: - numerator_sum = y_data[o] - denominator_sum = 1 - break - weight = distance ** (-3) - numerator_sum = numerator_sum + y_data[o] * weight - denominator_sum = denominator_sum + weight - ans[i] = numerator_sum / denominator_sum - return ans if len(ans) > 1 else ans[0] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + + sub_matrix = x_data - x + distances_squared = np.sum(sub_matrix**2, axis=2) + + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False + + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] + + return result if len(result) > 1 else result[0] + # Returns value for polynomial interpolation function type elif self.__interpolation__ == "polynomial": if isinstance(args[0], (int, float)): From 3cd456f945c993afde829ff6c7d3802a6c2d8d08 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:10:01 -0300 Subject: [PATCH 2/8] TST: improve testing of higher dimensional functions. --- tests/test_function.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/tests/test_function.py b/tests/test_function.py index 183888b96..5b7daab5a 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -262,6 +262,7 @@ def test_multivariable_function_plot(mock_show): @pytest.mark.parametrize( "x,y,z_expected", [ + (0, 0, 1), (1, 0, 0), (0, 1, 0), (0, 0, 1), @@ -270,15 +271,40 @@ def test_multivariable_function_plot(mock_show): ([0, 0.5], [0, 0.5], [1, 1 / 3]), ], ) -def test_shepard_interpolation(x, y, z_expected): +def test_2d_shepard_interpolation(x, y, z_expected): """Test the shepard interpolation method of the Function class.""" # Test plane x + y + z = 1 source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] func = Function(source=source, inputs=["x", "y"], outputs=["z"]) z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() assert np.isclose(z, z_expected, atol=1e-8).all() +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize("other", [1, 0.1, np.int_(1), np.float_(0.1), np.array([1])]) def test_sum_arithmetic_priority(other): """Test the arithmetic priority of the add operation of the Function class, From 954be1dad6c1dba92fa53ddf1fd8844cd750f549 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:22:51 -0300 Subject: [PATCH 3/8] FIX: interpolation return type if single value array. --- rocketpy/mathutils/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 8d97565ef..0e8213ac1 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -508,7 +508,7 @@ def get_value_opt_multiple(*args): result[valid_indexes] = numerator_sum / denominator_sum result[~valid_indexes] = y_data[zero_distances[1]] - return result + return result if len(result) > 1 else result[0] get_value_opt = get_value_opt_multiple From 7eab8ed751f1443070fbb410c99c38b75e7be90b Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:24:56 -0300 Subject: [PATCH 4/8] TST: move testing to unit module. --- tests/test_function.py | 46 ------------------------------------- tests/unit/test_function.py | 46 +++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 46 deletions(-) diff --git a/tests/test_function.py b/tests/test_function.py index 5b7daab5a..0b31d21ca 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -259,52 +259,6 @@ def test_multivariable_function_plot(mock_show): assert func.plot() == None -@pytest.mark.parametrize( - "x,y,z_expected", - [ - (0, 0, 1), - (1, 0, 0), - (0, 1, 0), - (0, 0, 1), - (0.5, 0.5, 1 / 3), - (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), - ([0, 0.5], [0, 0.5], [1, 1 / 3]), - ], -) -def test_2d_shepard_interpolation(x, y, z_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z = 1 - source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] - func = Function(source=source, inputs=["x", "y"], outputs=["z"]) - z = func(x, y) - z_opt = func.get_value_opt(x, y) - assert np.isclose(z, z_opt, atol=1e-8).all() - assert np.isclose(z, z_expected, atol=1e-8).all() - - -@pytest.mark.parametrize( - "x,y,z,w_expected", - [ - (0, 0, 0, 1), - (1, 0, 0, 0), - (0, 1, 0, 0), - (0, 0, 1, 0), - (0.5, 0.5, 0.5, 1 / 4), - (0.25, 0.25, 0.25, 0.700632626832), - ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), - ], -) -def test_3d_shepard_interpolation(x, y, z, w_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z + w = 1 - source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] - func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) - w = func(x, y, z) - w_opt = func.get_value_opt(x, y, z) - assert np.isclose(w, w_opt, atol=1e-8).all() - assert np.isclose(w_expected, w, atol=1e-8).all() - - @pytest.mark.parametrize("other", [1, 0.1, np.int_(1), np.float_(0.1), np.array([1])]) def test_sum_arithmetic_priority(other): """Test the arithmetic priority of the add operation of the Function class, diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index e41d30f04..ea52a88e6 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -48,6 +48,52 @@ def test_integral_spline_interpolation(request, func, a, b): ) +@pytest.mark.parametrize( + "x,y,z_expected", + [ + (0, 0, 1), + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 1 / 3), + (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), + ([0, 0.5], [0, 0.5], [1, 1 / 3]), + ], +) +def test_2d_shepard_interpolation(x, y, z_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z = 1 + source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] + func = Function(source=source, inputs=["x", "y"], outputs=["z"]) + z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() + assert np.isclose(z, z_expected, atol=1e-8).all() + + +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize( "func_input, derivative_input, expected_derivative", [ From 8fdfa555430795895be5eedd6710f70468928ba1 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:29:06 -0300 Subject: [PATCH 5/8] MNT: update CHANGELOG with changes. --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a8cbb9ed6..33926e5ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,7 +32,7 @@ straightforward as possible. ### Added -- +- ENH: Shepard Optimized Interpolation - Multiple Inputs Support [#515](https://github.com/RocketPy-Team/RocketPy/pull/515) ### Changed From 34ed384b50aff46a02173d98d8f3d34024d20769 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Tue, 19 Dec 2023 14:59:53 -0300 Subject: [PATCH 6/8] MNT: solve minor test placement issue. --- tests/test_function.py | 50 +++++++++++++++++++++++++++++++++++++ tests/unit/test_function.py | 46 ---------------------------------- 2 files changed, 50 insertions(+), 46 deletions(-) diff --git a/tests/test_function.py b/tests/test_function.py index 0b31d21ca..9cc2e68f0 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -215,6 +215,56 @@ def test_multivariable_dataset(a, b): assert np.isclose(func(a, b), a + b, atol=1e-6) +@pytest.mark.parametrize( + "x,y,z_expected", + [ + (0, 0, 1), + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 1 / 3), + (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), + ([0, 0.5], [0, 0.5], [1, 1 / 3]), + ], +) +def test_2d_shepard_interpolation(x, y, z_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z = 1 + source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] + func = Function( + source=source, inputs=["x", "y"], outputs=["z"], interpolation="shepard" + ) + z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() + assert np.isclose(z_expected, z, atol=1e-8).all() + + +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function( + source=source, inputs=["x", "y", "z"], outputs=["w"], interpolation="shepard" + ) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize("a", [-1, -0.5, 0, 0.5, 1]) @pytest.mark.parametrize("b", [-1, -0.5, 0, 0.5, 1]) def test_multivariable_function(a, b): diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index ea52a88e6..e41d30f04 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -48,52 +48,6 @@ def test_integral_spline_interpolation(request, func, a, b): ) -@pytest.mark.parametrize( - "x,y,z_expected", - [ - (0, 0, 1), - (1, 0, 0), - (0, 1, 0), - (0, 0, 1), - (0.5, 0.5, 1 / 3), - (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), - ([0, 0.5], [0, 0.5], [1, 1 / 3]), - ], -) -def test_2d_shepard_interpolation(x, y, z_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z = 1 - source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] - func = Function(source=source, inputs=["x", "y"], outputs=["z"]) - z = func(x, y) - z_opt = func.get_value_opt(x, y) - assert np.isclose(z, z_opt, atol=1e-8).all() - assert np.isclose(z, z_expected, atol=1e-8).all() - - -@pytest.mark.parametrize( - "x,y,z,w_expected", - [ - (0, 0, 0, 1), - (1, 0, 0, 0), - (0, 1, 0, 0), - (0, 0, 1, 0), - (0.5, 0.5, 0.5, 1 / 4), - (0.25, 0.25, 0.25, 0.700632626832), - ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), - ], -) -def test_3d_shepard_interpolation(x, y, z, w_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z + w = 1 - source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] - func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) - w = func(x, y, z) - w_opt = func.get_value_opt(x, y, z) - assert np.isclose(w, w_opt, atol=1e-8).all() - assert np.isclose(w_expected, w, atol=1e-8).all() - - @pytest.mark.parametrize( "func_input, derivative_input, expected_derivative", [ From 4ed16bee1d88fe781b6fad5bff73aa9d49720898 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Tue, 19 Dec 2023 18:41:00 -0300 Subject: [PATCH 7/8] MNT: avoid code interpolation code repetition. --- rocketpy/mathutils/function.py | 93 ++++++++++++++++------------------ 1 file changed, 43 insertions(+), 50 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 0e8213ac1..ceeb6a2e2 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -484,31 +484,7 @@ def get_value_opt(x): elif self.__interpolation__ == "shepard": # change the function's name to avoid mypy's error def get_value_opt_multiple(*args): - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - - arg_stack = np.column_stack(args) - arg_qty, arg_dim = arg_stack.shape - result = np.zeros(arg_qty) - - # Reshape to vectorize calculations - x = arg_stack.reshape(arg_qty, 1, arg_dim) - - sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=2) - - # Remove zero distances from further calculations - zero_distances = np.where(distances_squared == 0) - valid_indexes = np.ones(arg_qty, dtype=bool) - valid_indexes[zero_distances[0]] = False - - weights = distances_squared[valid_indexes] ** (-1.5) - numerator_sum = np.sum(y_data * weights, axis=1) - denominator_sum = np.sum(weights, axis=1) - result[valid_indexes] = numerator_sum / denominator_sum - result[~valid_indexes] = y_data[zero_distances[1]] - - return result if len(result) > 1 else result[0] + return self.__interpolate_shepard__(args) get_value_opt = get_value_opt_multiple @@ -880,31 +856,7 @@ def get_value(self, *args): # Returns value for shepard interpolation elif self.__interpolation__ == "shepard": - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - - arg_stack = np.column_stack(args) - arg_qty, arg_dim = arg_stack.shape - result = np.zeros(arg_qty) - - # Reshape to vectorize calculations - x = arg_stack.reshape(arg_qty, 1, arg_dim) - - sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=2) - - # Remove zero distances from further calculations - zero_distances = np.where(distances_squared == 0) - valid_indexes = np.ones(arg_qty, dtype=bool) - valid_indexes[zero_distances[0]] = False - - weights = distances_squared[valid_indexes] ** (-1.5) - numerator_sum = np.sum(y_data * weights, axis=1) - denominator_sum = np.sum(weights, axis=1) - result[valid_indexes] = numerator_sum / denominator_sum - result[~valid_indexes] = y_data[zero_distances[1]] - - return result if len(result) > 1 else result[0] + return self.__interpolate_shepard__(args) # Returns value for polynomial interpolation function type elif self.__interpolation__ == "polynomial": @@ -1656,6 +1608,47 @@ def __interpolate_akima__(self): coeffs[4 * i : 4 * i + 4] = np.linalg.solve(matrix, result) self.__akima_coefficients__ = coeffs + def __interpolate_shepard__(self, args): + """Calculates the shepard interpolation from the given arguments. + The shepard interpolation is computed by a inverse distance weighting + in a vectorized manner. + + Parameters + ---------- + args : scalar, list + Values where the Function is to be evaluated. + + Returns + ------- + result : scalar, list + The result of the interpolation. + """ + x_data = self.source[:, 0:-1] # Support for N-Dimensions + y_data = self.source[:, -1] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + + sub_matrix = x_data - x + distances_squared = np.sum(sub_matrix**2, axis=2) + + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False + + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] + + return result if len(result) > 1 else result[0] + def __neg__(self): """Negates the Function object. The result has the same effect as multiplying the Function by -1. From acbc14dae924d874e6a2803e2f494be92eec7f09 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 18 Jan 2024 20:56:04 -0300 Subject: [PATCH 8/8] MNT: remove repeated test case. --- tests/test_function.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_function.py b/tests/test_function.py index 9cc2e68f0..352feb7b3 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -218,7 +218,6 @@ def test_multivariable_dataset(a, b): @pytest.mark.parametrize( "x,y,z_expected", [ - (0, 0, 1), (1, 0, 0), (0, 1, 0), (0, 0, 1),