From 2b4472307d5678c4d86b466fc0110b18417bb19c Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Sat, 14 Sep 2024 15:49:21 -0300 Subject: [PATCH 1/8] ENH: expand polation options for ND Functions. --- rocketpy/mathutils/function.py | 267 ++++++++++++++++++++++----------- 1 file changed, 183 insertions(+), 84 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 095b82268..e3f904a9b 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -16,6 +16,11 @@ import matplotlib.pyplot as plt import numpy as np from scipy import integrate, linalg, optimize +from scipy.interpolate import ( + LinearNDInterpolator, + NearestNDInterpolator, + RBFInterpolator, +) # Numpy 1.x compatibility, # TODO: remove these lines when all dependencies support numpy>=2.0.0 @@ -32,6 +37,7 @@ "akima": 2, "spline": 3, "shepard": 4, + "rbf": 5, } EXTRAPOLATION_TYPES = {"zero": 0, "natural": 1, "constant": 2} @@ -224,6 +230,8 @@ def set_source(self, source): # pylint: disable=too-many-statements else: # Evaluate dimension self.__dom_dim__ = source.shape[1] - 1 + self._domain = source[:, :-1] + self._image = source[:, -1] # set x and y. If Function is 2D, also set z if self.__dom_dim__ == 1: @@ -330,22 +338,29 @@ def set_extrapolation(self, method="constant"): def __set_interpolation_func(self): # pylint: disable=too-many-statements """Defines interpolation function used by the Function. Each - interpolation method has its own function with exception of shepard, - which has its interpolation/extrapolation function defined in - ``Function.__interpolate_shepard__``. The function is stored in - the attribute _interpolation_func.""" + interpolation method has its own function`. + The function is stored in the attribute _interpolation_func.""" interpolation = INTERPOLATION_TYPES[self.__interpolation__] if interpolation == 0: # linear + if self.__dom_dim__ == 1: - def linear_interpolation( - x, x_min, x_max, x_data, y_data, coeffs - ): # pylint: disable=unused-argument - x_interval = bisect_left(x_data, x) - x_left = x_data[x_interval - 1] - y_left = y_data[x_interval - 1] - dx = float(x_data[x_interval] - x_left) - dy = float(y_data[x_interval] - y_left) - return (x - x_left) * (dy / dx) + y_left + def linear_interpolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + x_interval = bisect_left(x_data, x) + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + return (x - x_left) * (dy / dx) + y_left + + else: + interpolator = LinearNDInterpolator(self._domain, self._image) + + def linear_interpolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + return interpolator(x) self._interpolation_func = linear_interpolation @@ -383,8 +398,41 @@ def spline_interpolation( self._interpolation_func = spline_interpolation - elif interpolation == 4: # shepard does not use interpolation function - self._interpolation_func = None + elif interpolation == 4: # shepard + + # pylint: disable=unused-argument + def shepard_interpolation(x, x_min, x_max, x_data, y_data, _): + arg_qty, arg_dim = x.shape + result = np.empty(arg_qty) + x = x.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 + + self._interpolation_func = shepard_interpolation + + elif interpolation == 5: # RBF + + interpolator = RBFInterpolator(self._domain, self._image, neighbors=100) + + def rbf_interpolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + return interpolator(x) + + self._interpolation_func = rbf_interpolation def __set_extrapolation_func(self): # pylint: disable=too-many-statements """Defines extrapolation function used by the Function. Each @@ -393,10 +441,7 @@ def __set_extrapolation_func(self): # pylint: disable=too-many-statements interpolation = INTERPOLATION_TYPES[self.__interpolation__] extrapolation = EXTRAPOLATION_TYPES[self.__extrapolation__] - if interpolation == 4: # shepard does not use extrapolation function - self._extrapolation_func = None - - elif extrapolation == 0: # zero + if extrapolation == 0: # zero def zero_extrapolation( x, x_min, x_max, x_data, y_data, coeffs @@ -407,15 +452,27 @@ def zero_extrapolation( elif extrapolation == 1: # natural if interpolation == 0: # linear - def natural_extrapolation( - x, x_min, x_max, x_data, y_data, coeffs - ): # pylint: disable=unused-argument - x_interval = 1 if x < x_min else -1 - x_left = x_data[x_interval - 1] - y_left = y_data[x_interval - 1] - dx = float(x_data[x_interval] - x_left) - dy = float(y_data[x_interval] - y_left) - return (x - x_left) * (dy / dx) + y_left + if self.__dom_dim__ == 1: + + def natural_extrapolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + x_interval = 1 if x < x_min else -1 + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + return (x - x_left) * (dy / dx) + y_left + + else: + interpolator = RBFInterpolator( + self._domain, self._image, neighbors=100 + ) + + def natural_extrapolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + return interpolator(x) elif interpolation == 1: # polynomial @@ -445,13 +502,55 @@ def natural_extrapolation( x = x - x_data[-2] return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + elif interpolation == 4: # shepard + + # pylint: disable=unused-argument + def natural_extrapolation(x, x_min, x_max, x_data, y_data, _): + arg_qty, arg_dim = x.shape + result = np.empty(arg_qty) + x = x.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 + + elif interpolation == 5: # RBF + + interpolator = RBFInterpolator(self._domain, self._image, neighbors=100) + + def natural_extrapolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + return interpolator(x) + self._extrapolation_func = natural_extrapolation elif extrapolation == 2: # constant - def constant_extrapolation( - x, x_min, x_max, x_data, y_data, coeffs - ): # pylint: disable=unused-argument - return y_data[0] if x < x_min else y_data[-1] + if self.__dom_dim__ == 1: + + def constant_extrapolation( + x, x_min, x_max, x_data, y_data, coeffs + ): # pylint: disable=unused-argument + return y_data[0] if x < x_min else y_data[-1] + + else: + + extrapolator = NearestNDInterpolator(self._domain, self._image) + + def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + # pylint: disable=unused-argument + return extrapolator(x) self._extrapolation_func = constant_extrapolation @@ -496,10 +595,41 @@ def __get_value_opt_1d(self, x): return y def __get_value_opt_nd(self, *args): - """Evaluate the Function at a single point (x, y, z). This method is - used when the Function is N-D.""" - # always use shepard for N-D functions - return self.__interpolate_shepard__(args) + """Evaluate the Function in a vectorized fashion for ND domains. + + Parameters + ---------- + args : tuple + Values where the Function is to be evaluated. + + Returns + ------- + result : scalar, ndarray + Value of the Function at the specified points. + """ + args = np.column_stack(args) + arg_qty = len(args) + result = np.empty(arg_qty) + + min_domain = self._domain.T.min(axis=1) + max_domain = self._domain.T.max(axis=1) + + lower, upper = args < min_domain, args > max_domain + extrap = np.logical_or(lower.any(axis=1), upper.any(axis=1)) + + if extrap.any(): + result[extrap] = self._extrapolation_func( + args[extrap], min_domain, max_domain, self._domain, self._image, None + ) + if (~extrap).any(): + result[~extrap] = self._interpolation_func( + args[~extrap], min_domain, max_domain, self._domain, self._image, None + ) + + if arg_qty == 1: + return float(result[0]) + + return result def set_discrete( self, @@ -898,7 +1028,7 @@ def get_value(self, *args): if all(isinstance(arg, Iterable) for arg in args): return [self.source(*arg) for arg in zip(*args)] - elif self.__dom_dim__ > 1: # deals with nd functions and shepard interp + elif self.__dom_dim__ > 1: # deals with nd functions return self.get_value_opt(*args) # Returns value for other interpolation type @@ -1765,47 +1895,6 @@ 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. @@ -3238,14 +3327,17 @@ def __validate_interpolation(self, interpolation): interpolation = "spline" ## multiple dimensions elif self.__dom_dim__ > 1: - if interpolation not in [None, "shepard"]: + if interpolation is None: + interpolation = "shepard" + if interpolation not in ["shepard", "linear", "rbf"]: warnings.warn( ( - "Interpolation method set to 'shepard'. Only 'shepard' " - "interpolation is supported for multiple dimensions." + "Interpolation method set to 'shepard'. The methods " + "'linear' and 'shepard' are supported for multiple " + "dimensions." ), ) - interpolation = "shepard" + interpolation = "shepard" return interpolation def __validate_extrapolation(self, extrapolation): @@ -3261,12 +3353,19 @@ def __validate_extrapolation(self, extrapolation): ## multiple dimensions elif self.__dom_dim__ > 1: - if extrapolation not in [None, "natural"]: + if extrapolation is None: + extrapolation = "natural" + if extrapolation == "natural" and self.__interpolation__ == "linear": + warnings.warn( + "Extrapolation 'natural' is not supported for multidimensional " + "linear interpolation. 'rbf' will be used to extrapolate." + ) + elif extrapolation not in ["constant", "natural", "zero"]: warnings.warn( "Extrapolation method set to 'natural'. Other methods " "are not supported yet." ) - extrapolation = "natural" + extrapolation = "natural" return extrapolation From 66e9b8ec9110c0a1ca7041d5254408444376fde2 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Sat, 14 Sep 2024 15:49:48 -0300 Subject: [PATCH 2/8] TST: introduce testing for ND linear and rbf interpolations. --- tests/unit/test_function.py | 124 ++++++++++++++++++++++++++++++------ 1 file changed, 104 insertions(+), 20 deletions(-) diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 9efb64c0c..df540c1ae 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -502,6 +502,110 @@ def test_3d_shepard_interpolation(x, y, z, w_expected): assert np.isclose(w_expected, w, atol=1e-8).all() +@pytest.mark.parametrize( + "x,y,z_expected", + [ + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 0), + (0.25, 0.25, 0.5), + ([0, 0.5], [0, 0.5], [1, 0]), + ], +) +def test_2d_rbf_interpolation(x, y, z_expected): + """Test the rbf 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="rbf" + ) + 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, -0.5), + (0.25, 0.25, 0.25, 0.25), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, -0.5]), + ], +) +def test_3d_rbf_interpolation(x, y, z, w_expected): + """Test the rbf 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="rbf" + ) + 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( + "x,y,z_expected", + [ + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 0), + (0.25, 0.25, 0.5), + ([0, 0.5], [0, 0.5], [1, 0]), + ], +) +def test_2d_linear_interpolation(x, y, z_expected): + """Test the linear 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="linear" + ) + 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", + [ + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0, 0, 0, 1), + (0.5, 0.5, 0.5, -0.5), + (0.25, 0.25, 0.25, 0.25), + ([0, 0.25], [0, 0.25], [0, 0.25], [1, 0.25]), + ], +) +def test_3d_linear_interpolation(x, y, z, w_expected): + """Test the linear 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), + (0.5, 0.5, 0.5, -0.5), + ] + func = Function( + source=source, inputs=["x", "y", "z"], outputs=["w"], interpolation="linear" + ) + 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_multivariate_function(a, b): @@ -565,26 +669,6 @@ def test_set_discrete_based_on_2d_model(func_2d_from_csv): assert discretized_func.__extrapolation__ == func_2d_from_csv.__extrapolation__ -@pytest.mark.parametrize( - "x,y,z_expected", - [ - (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_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) - assert np.isclose(z, z_expected, atol=1e-8).all() - - @pytest.mark.parametrize("other", [1, 0.1, np.int_(1), np.float64(0.1), np.array([1])]) def test_sum_arithmetic_priority(other): """Test the arithmetic priority of the add operation of the Function class, From d61c5c2761a10342ce6fe8eb3a62efb36c743a49 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Sat, 14 Sep 2024 16:17:00 -0300 Subject: [PATCH 3/8] MNT: add Function interpolation options to CHANGELOG. --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f54edd029..6f6247ee3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ Attention: The newest changes should be on top --> ### Added +- ENH: Expand Polation Options for ND Functions. [#691](https://github.com/RocketPy-Team/RocketPy/pull/691) - DOC : Cavour Flight Example [#682](https://github.com/RocketPy-Team/RocketPy/pull/682) - DOC: Halcyon Flight Example [#681](https://github.com/RocketPy-Team/RocketPy/pull/681) - ENH: Adds GenericMotor.load_from_eng_file() method [#676](https://github.com/RocketPy-Team/RocketPy/pull/676) From 86458aab47331eaeb59f9c61795f84548c5eb90c Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 19 Sep 2024 18:54:50 -0300 Subject: [PATCH 4/8] DOC: improve docstring regarding N-D polations. --- rocketpy/mathutils/function.py | 60 ++++++++++++++++++++-------------- 1 file changed, 36 insertions(+), 24 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index e3f904a9b..97b5eded6 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -92,14 +92,17 @@ def __init__( interpolation : string, optional Interpolation method to be used if source type is ndarray. For 1-D functions, linear, polynomial, akima and spline are - supported. For N-D functions, only shepard is supported. - Default for 1-D functions is spline. + supported. For N-D functions, linear, shepard and rbf are + supported. + Default for 1-D functions is spline and for N-D functions is + shepard. extrapolation : string, optional Extrapolation method to be used if source type is ndarray. Options are 'natural', which keeps interpolation, 'constant', - which returns the value of the function at the edge of the interval, - and 'zero', which returns zero for all points outside of source - range. Default for 1-D functions is constant. + which returns the value of the function at the nearest edge of + the domain, and 'zero', which returns zero for all points outside + of source range. Default for 1-D functions is constant and for + N-D functions is natural. title : string, optional Title to be displayed in the plots' figures. If none, the title will be constructed using the inputs and outputs arguments in the form @@ -285,8 +288,10 @@ def set_interpolation(self, method="spline"): method : string, optional Interpolation method to be used if source type is ndarray. For 1-D functions, linear, polynomial, akima and spline is - supported. For N-D functions, only shepard is supported. - Default is 'spline'. + supported. For N-D functions, linear, shepard and rbf are + supported. + Default for 1-D functions is spline and for N-D functions is + shepard. Returns ------- @@ -322,9 +327,10 @@ def set_extrapolation(self, method="constant"): extrapolation : string, optional Extrapolation method to be used if source type is ndarray. Options are 'natural', which keeps interpolation, 'constant', - which returns the value of the function at the edge of the interval, - and 'zero', which returns zero for all points outside of source - range. Default is 'constant'. + which returns the value of the function at the nearest edge of + the domain, and 'zero', which returns zero for all points outside + of source range. Default for 1-D functions is constant and for + N-D functions is natural. Returns ------- @@ -595,7 +601,7 @@ def __get_value_opt_1d(self, x): return y def __get_value_opt_nd(self, *args): - """Evaluate the Function in a vectorized fashion for ND domains. + """Evaluate the Function in a vectorized fashion for N-D domains. Parameters ---------- @@ -660,15 +666,18 @@ def set_discrete( Number of samples to be taken from inside range. Default is 200. interpolation : string Interpolation method to be used if source type is ndarray. - For 1-D functions, linear, polynomial, akima and spline is - supported. For N-D functions, only shepard is supported. - Default is 'spline'. + For 1-D functions, linear, polynomial, akima and spline are + supported. For N-D functions, linear, shepard and rbf are + supported. + Default for 1-D functions is spline and for N-D functions is + shepard. extrapolation : string, optional Extrapolation method to be used if source type is ndarray. Options are 'natural', which keeps interpolation, 'constant', - which returns the value of the function at the edge of the interval, - and 'zero', which returns zero for all points outside of source - range. Default is 'constant'. + which returns the value of the function at the nearest edge of + the domain, and 'zero', which returns zero for all points outside + of source range. Default for 1-D functions is constant and for + N-D functions is natural. one_by_one : boolean, optional If True, evaluate Function in each sample point separately. If False, evaluates Function in vectorized form. Default is True. @@ -3329,12 +3338,12 @@ def __validate_interpolation(self, interpolation): elif self.__dom_dim__ > 1: if interpolation is None: interpolation = "shepard" - if interpolation not in ["shepard", "linear", "rbf"]: + if interpolation.lower() not in ["shepard", "linear", "rbf"]: warnings.warn( ( "Interpolation method set to 'shepard'. The methods " - "'linear' and 'shepard' are supported for multiple " - "dimensions." + "'linear', 'shepard' and 'rbf' are supported for " + "multiple dimensions." ), ) interpolation = "shepard" @@ -3355,15 +3364,18 @@ def __validate_extrapolation(self, extrapolation): elif self.__dom_dim__ > 1: if extrapolation is None: extrapolation = "natural" - if extrapolation == "natural" and self.__interpolation__ == "linear": + if ( + extrapolation.lower() == "natural" + and self.__interpolation__ == "linear" + ): warnings.warn( "Extrapolation 'natural' is not supported for multidimensional " "linear interpolation. 'rbf' will be used to extrapolate." ) - elif extrapolation not in ["constant", "natural", "zero"]: + elif extrapolation.lower() not in ["constant", "natural", "zero"]: warnings.warn( - "Extrapolation method set to 'natural'. Other methods " - "are not supported yet." + "Extrapolation method set to 'natural' because the " + f"{extrapolation} method is not supported." ) extrapolation = "natural" return extrapolation From c04d401947d716f40f594318c482774603921bd2 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 19 Sep 2024 19:26:16 -0300 Subject: [PATCH 5/8] MNT: add check for degenerate domains. --- rocketpy/mathutils/function.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 97b5eded6..1ec1fb3c7 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -3239,6 +3239,15 @@ def __validate_source(self, source): # pylint: disable=too-many-statements raise ValueError( "Source must be a 2D array in the form [[x1, x2 ..., xn, y], ...]." ) + + source_len, source_dim = source.shape + + if source_len < source_dim: + raise ValueError( + "Too few data points to define a domain. The number of rows " + "must be greater than or equal to the number of columns." + ) + return source if isinstance(source, NUMERICAL_TYPES): From 12ac468ce8824ade4fd231563be7cee5061881c4 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 19 Sep 2024 19:38:09 -0300 Subject: [PATCH 6/8] STY: fix trailing whitespaces in docstrings. --- rocketpy/mathutils/function.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 1ec1fb3c7..7e0acf8da 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -94,14 +94,14 @@ def __init__( For 1-D functions, linear, polynomial, akima and spline are supported. For N-D functions, linear, shepard and rbf are supported. - Default for 1-D functions is spline and for N-D functions is + Default for 1-D functions is spline and for N-D functions is shepard. extrapolation : string, optional Extrapolation method to be used if source type is ndarray. Options are 'natural', which keeps interpolation, 'constant', which returns the value of the function at the nearest edge of - the domain, and 'zero', which returns zero for all points outside - of source range. Default for 1-D functions is constant and for + the domain, and 'zero', which returns zero for all points outside + of source range. Default for 1-D functions is constant and for N-D functions is natural. title : string, optional Title to be displayed in the plots' figures. If none, the title will From 339d7f02d0457fcc9a1fd78570852dab2d952ac4 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 19 Sep 2024 19:53:40 -0300 Subject: [PATCH 7/8] MNT: fix doctest typing casting. --- rocketpy/mathutils/function.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 7e0acf8da..b75a00347 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -1004,13 +1004,13 @@ def get_value(self, *args): ... [(0, 0, 0), (1, 1, 1), (1, 2, 2), (2, 4, 8), (3, 9, 27)] ... ) >>> f4.get_value(1, 1) - np.float64(1.0) + 1.0 >>> f4.get_value(2, 4) - np.float64(8.0) + 8.0 >>> abs(f4.get_value(1, 1.5) - 1.5) < 1e-2 # the interpolation is not perfect - np.True_ + True >>> f4.get_value(3, 9) - np.float64(27.0) + 27.0 """ if len(args) != self.__dom_dim__: raise ValueError( From bae7caef22ac94a103c1ed345a5654e4072b24d5 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Fri, 20 Sep 2024 14:59:10 -0300 Subject: [PATCH 8/8] MNT: improve N-D extrapolation docs. --- rocketpy/mathutils/function.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index b75a00347..9e6e07b9d 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -101,8 +101,11 @@ def __init__( Options are 'natural', which keeps interpolation, 'constant', which returns the value of the function at the nearest edge of the domain, and 'zero', which returns zero for all points outside - of source range. Default for 1-D functions is constant and for - N-D functions is natural. + of source range. + Multidimensional 'natural' extrapolation for 'linear' interpolation + use a 'rbf' algorithm for smoother results. + Default for 1-D functions is constant and for N-D functions + is natural. title : string, optional Title to be displayed in the plots' figures. If none, the title will be constructed using the inputs and outputs arguments in the form @@ -329,8 +332,11 @@ def set_extrapolation(self, method="constant"): Options are 'natural', which keeps interpolation, 'constant', which returns the value of the function at the nearest edge of the domain, and 'zero', which returns zero for all points outside - of source range. Default for 1-D functions is constant and for - N-D functions is natural. + of source range. + Multidimensional 'natural' extrapolation for 'linear' interpolation + use a 'rbf' algorithm for smoother results. + Default for 1-D functions is constant and for N-D functions + is natural. Returns ------- @@ -3373,15 +3379,7 @@ def __validate_extrapolation(self, extrapolation): elif self.__dom_dim__ > 1: if extrapolation is None: extrapolation = "natural" - if ( - extrapolation.lower() == "natural" - and self.__interpolation__ == "linear" - ): - warnings.warn( - "Extrapolation 'natural' is not supported for multidimensional " - "linear interpolation. 'rbf' will be used to extrapolate." - ) - elif extrapolation.lower() not in ["constant", "natural", "zero"]: + if extrapolation.lower() not in ["constant", "natural", "zero"]: warnings.warn( "Extrapolation method set to 'natural' because the " f"{extrapolation} method is not supported."