ENH: Expand Polation Options for ND Functions. #691
Conversation
phmbressan
added
the
Enhancement
|
I will update the Documentation regarding the changes shortly, but feature wise the PR should be ready to review. |
| @@ -32,6 +37,7 @@ | |||
| "akima": 2, | |||
| "spline": 3, | |||
| "shepard": 4, | |||
| "rbf": 5, | |||
There was a problem hiding this comment.
don't you have to also include the nearest option here?
There was a problem hiding this comment.
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.
| self._domain = source[:, :-1] | ||
| self._image = source[:, -1] |
There was a problem hiding this comment.
isn't these lines equivalent to x_array and y_array?
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.
rocketpy/mathutils/function.py
Outdated
| 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".
rocketpy/mathutils/function.py
Outdated
| 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 |
There was a problem hiding this comment.
Missing a check for "rbf" in this metod. Rigth now setting extrapolation to rbf with interpolation linear raises a warning
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: ... 'rbf' will be used to extrapolate.". If rbf is not an extrapolation method, this warning is wrong.
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:
"Extrapolation 'natural' for multidimensional linear "
"interpolation uses a 'rbf' algorithm for smoother results."
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.
There was a problem hiding this comment.
There seems to be a minimum amount of points for the ND linear interpolation to be used. Could this be checked for when initializing the Function object?
This will work:
alpha,mach,cL
0.10471975511965975,0.0,0.2094395102393195
0.10471975511965975,0.5,0.2094395102393195
0.20943951023931956,0.0,0.4188790204786391
0.20943951023931956,0.5,0.4188790204786391
But this won't:
alpha,mach,cL
0.10471975511965975,0.0,0.2094395102393195
0.20943951023931956,0.0,0.4188790204786391
Here is the error raised:
{
"name": "QhullError",
"message": "QH6214 qhull input error: not enough points(2) to construct initial simplex (need 4)
While executing: | qhull d Qbb Qc Qz Qt Q12
Options selected for Qhull 2019.1.r 2019/06/21:
run-id 213006174 delaunay Qbbound-last Qcoplanar-keep Qz-infinity-point
Qtriangulate Q12-allow-wide _pre-merge _zero-centrum Qinterior-keep
_maxoutside 0
",
"stack": "---------------------------------------------------------------------------
QhullError Traceback (most recent call last)
Cell In[36], line 1
----> 1 gennose = GenericSurface(
2 reference_area=np.pi * calisto.radius**2,
3 reference_length=2 * calisto.radius,
4 coefficients={
5 \"cL\": \"nose_cL.csv\",
6 \"cQ\": \"nose_cQ.csv\",
7 },
8 center_of_pressure=(0, 0, 0),
9 name=\"nose\",
10 )
File c:\\mateus\\github\\rocketpy\\rocketpy\\rocket\\aero_surface\\generic_surface.py:91, in GenericSurface.__init__(self, reference_area, reference_length, coefficients, center_of_pressure, name)
89 coefficients = self._complete_coefficients(coefficients, default_coefficients)
90 for coeff, coeff_value in coefficients.items():
---> 91 value = self._process_input(coeff_value, coeff)
92 setattr(self, coeff, value)
File c:\\mateus\\github\\rocketpy\\rocketpy\\rocket\\aero_surface\\generic_surface.py:330, in GenericSurface._process_input(self, input_data, coeff_name)
313 \"\"\"Process the input data, either as a CSV file or a callable function.
314
315 Parameters
(...)
326 pitch_rate, yaw_rate, roll_rate).
327 \"\"\"
328 if isinstance(input_data, str):
329 # Input is assumed to be a file path to a CSV
--> 330 return self.__load_csv(input_data, coeff_name)
331 elif isinstance(input_data, Function):
332 if input_data.__dom_dim__ != 7:
File c:\\mateus\\github\\rocketpy\\rocketpy\\rocket\\aero_surface\\generic_surface.py:420, in GenericSurface.__load_csv(self, file_path, coeff_name)
417 raise ValueError(f\"No independent variables found in {coeff_name} CSV.\")
419 # Initialize the CSV-based function
--> 420 csv_func = Function(
421 file_path,
422 interpolation='linear',
423 extrapolation='rbf',
424 )
426 # Create a mask for the presence of each independent variable
427 # save on self to avoid loss of scope
428 _mask = [1 if col in present_columns else 0 for col in independent_vars]
File c:\\mateus\\github\\rocketpy\\rocketpy\\mathutils\\function.py:137, in Function.__init__(self, source, inputs, outputs, interpolation, extrapolation, title)
134 self.__img_dim__ = 1 # always 1, here for backwards compatibility
136 # args must be passed from self.
--> 137 self.set_source(self.source)
138 self.set_inputs(self.__inputs__)
139 self.set_outputs(self.__outputs__)
File c:\\mateus\\github\\rocketpy\\rocketpy\\mathutils\\function.py:254, in Function.set_source(self, source)
251 self.get_value_opt = self.__get_value_opt_nd
253 self.source = source
--> 254 self.set_interpolation(self.__interpolation__)
255 self.set_extrapolation(self.__extrapolation__)
256 return self
File c:\\mateus\\github\\rocketpy\\rocketpy\\mathutils\\function.py:298, in Function.set_interpolation(self, method)
296 self.__interpolation__ = self.__validate_interpolation(method)
297 self.__update_interpolation_coefficients(self.__interpolation__)
--> 298 self.__set_interpolation_func()
299 return self
File c:\\mateus\\github\\rocketpy\\rocketpy\\mathutils\\function.py:358, in Function.__set_interpolation_func(self)
355 return (x - x_left) * (dy / dx) + y_left
357 else:
--> 358 interpolator = LinearNDInterpolator(self._domain, self._image)
360 def linear_interpolation(
361 x, x_min, x_max, x_data, y_data, coeffs
362 ): # pylint: disable=unused-argument
363 return interpolator(x)
File interpnd.pyx:301, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__()
File interpnd.pyx:92, in scipy.interpolate.interpnd.NDInterpolatorBase.__init__()
File interpnd.pyx:305, in scipy.interpolate.interpnd.LinearNDInterpolator._calculate_triangulation()
File _qhull.pyx:1885, in scipy.spatial._qhull.Delaunay.__init__()
File _qhull.pyx:352, in scipy.spatial._qhull._Qhull.__init__()
QhullError: QH6214 qhull input error: not enough points(2) to construct initial simplex (need 4)
While executing: | qhull d Qbb Qc Qz Qt Q12
Options selected for Qhull 2019.1.r 2019/06/21:
run-id 213006174 delaunay Qbbound-last Qcoplanar-keep Qz-infinity-point
Qtriangulate Q12-allow-wide _pre-merge _zero-centrum Qinterior-keep
_maxoutside 0
"
}
I have introduced a check into commit c04d401. Let me know if this solution is enough. |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #691 +/- ##
===========================================
+ Coverage 75.42% 75.45% +0.03%
===========================================
Files 96 96
Lines 10841 10887 +46
===========================================
+ Hits 8177 8215 +38
- Misses 2664 2672 +8 ☔ View full report in Codecov by Sentry. |
Seems good enough to me |
Pull request type
Checklist
black rocketpy/ tests/) has passed locallypytest tests -m slow --runslow) have passed locallyCHANGELOG.mdhas been updated (if relevant)Current behavior
Options of ND interpolation and extrapolation of ND Functions are rather limited. Only "shepard" is supported, which might not be adequate for some datasets.
New behavior
This PR introduces the following polation methods for ND Functions:
Interpolation:
Extrapolation:
Note: I could not find an easy / reasonable approach for strictly linear extrapolation, since the dataset dealt by the Function class is generally "unstructed" (check out scipy tools) which prevent the usage of utilities that make use of "structed" data (like interpn). It should be possible to build a Delaunay triangulation with scipy tools to find the nearest points and manually build the linear extrapolation, but the accuracy and applicability of the results will likely vary a lot.
Therefore, due to having more smooth and predictable results, I have decided to apply "rbf" extrapolation to the "linear" option (with a explicit warning). This might be worth of an improvement in the future.
Result Examples
Sample plots comparing interpolation results for the paraboloid$f(x,y) = x^2 + y^2$ :
Lambda Function
Rbf Interpolation
Linear Interpolation
Shepard Interpolation
Breaking change