ENH: Expand Polation Options for ND Functions. #691
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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] | ||
|
Comment on lines
+239
to
+240
There was a problem hiding this comment. isn't these lines equivalent to There was a problem hiding this comment. No, the y_array is the second input column. Therefore, in 1-D cases y_array is the output. For N-D it is only a particular column of the input domain. |
||
|
|
||
| # 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`. | ||
| 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( | ||
phmbressan marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| 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 | ||
|
|
||
| # 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 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 | ||
|
|
||
| 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 | ||
|
|
||
| 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 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 | ||
| 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 __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 is None: | ||
| interpolation = "shepard" | ||
| if interpolation not in ["shepard", "linear", "rbf"]: | ||
| warnings.warn( | ||
| ( | ||
| "Interpolation method set to 'shepard'. The methods " | ||
| "'linear' and 'shepard' are supported for multiple " | ||
| "dimensions." | ||
| ), | ||
| ) | ||
| interpolation = "shepard" | ||
phmbressan marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| 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 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." | ||
|
There was a problem hiding this comment. it says that rbf will be used, but I don't see where you are setting extrapolation to rbf There was a problem hiding this comment. No set is needed, since there is not "rbf" extrapolation option, only "constant", "natural" and "zero". The warning is there to point out that the "natural" extrapolation of the "linear" interpolation is an rbf algorithm, as I pointed out why is the PR description. I am open to ideas on improving this, but I found the way it is implemented the most consistent regarding the current behavior of interpolation and extrapolation settings of 1-D Functions. Remember, there is no "spline" or "linear" extrapolation option for 1-D Functions, those are only supported through "natural". |
||
| ) | ||
| elif extrapolation not in ["constant", "natural", "zero"]: | ||
| warnings.warn( | ||
| "Extrapolation method set to 'natural'. Other methods " | ||
phmbressan marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| "are not supported yet." | ||
| ) | ||
| extrapolation = "natural" | ||
| return extrapolation | ||
|
There was a problem hiding this comment. Missing a check for There was a problem hiding this comment. Yes, indeed, this happens because "rbf" is not a valid extrapolation input: the same way one cannot input "linear" or "spline" for 1-D Functions, a warning will also occur. As it is for 1-D Functions, the only extrapolation inputs are "constant", "natural" and "zero". There was a problem hiding this comment. The warning message says: This warning is also raised if interpolation="linear" and extrapolation="natural" There was a problem hiding this comment. The warning should only be raised if interpolation is "linear" and extrapolation "natural". The warning is there to show that a 'rbf' algorithm is used to extrapolate in this case. The reasoning behind it I have explained in the PR description and is detailed in scipy issue 6396. I can change the warning to be more precise and better phrased, like: There was a problem hiding this comment. I don't think there should be warning at all in this case. Just mention this behaviour in the docstrings. With the way it is, there will be warnings whenever we define a ND function inside the code. So there will be a lot of warnings for a desired behaviour There was a problem hiding this comment. Makes total sense, I believe the warning was a bit overzealous in this case. I introduced this behavior in commit bae7cae. |
||
|
|
||
|
|
||
|
|
||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
don't you have to also include the
nearestoption here?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As explained a bit better in the docs commit 86458aa, constant and nearest are basically two sides of the same coin, i.e., they are the same (the constant for 1-D is already the nearest in action). Therefore, I preferred to keep the same naming for the extrapolation types, since I believe adding a "nearest" would be redundant.