ENH: add more interpolation methods to RegularGridInterpolator. #7903

thearn · 2017-09-20T23:38:37Z

This includes the addition of 5 3 new spline-based interpolators for
scipy.interpolate.RegularGridInterpolator. These allow for smooth interpolation
on a regular grid in an arbitrary number of dimensions. These also provide
gradient calculations of the interpolated result with respect to an
evaluation point or set of points.

This includes the addition of 5 new spline-based interpolators for scipy.interpolate.RegularGridInterpolator. These allow for smooth interpolation on a regular grid in an arbitrary number of dimensions. These also provide gradient calculations of the interpolated result with respect to an evaluation point or set of points.

thearn · 2017-09-20T23:38:59Z

[Updated] A bit more detail: This includes the addition of 5 3 new spline-based interpolators for
scipy.interpolate.RegularGridInterpolator:

slinear
~~quadratic~~
cubic
~~quartic~~
quintic

These each allow for smooth interpolation
on a regular grid in an arbitrary number of dimensions. These methods also allow for
gradient calculations of the interpolated result with respect to an
evaluation point or set of points, via call to a new method called gradient:

interp = RegularGridInterpolator(points, values, method='cubic')
result = interp(x)
deriv = interp.gradient(x)

Note that gradients cannot be computed for the two currently existing interpolation methods, linear and nearest.

These additional 5 methods are examples of tensor spline interpolation. They work by applying instances of ~~UnivariateSpline~~ BSpline to each 1D sequence present in the last dimensions of an array of data at the point of interest, "collapsing" these values into the next dimension, and continues
until a final value (and its gradient) is reached. The interpolation can be understood as being conducted by functions in a basis constructed from the tensor product of the 1D splines.

This is useful to my work for various surrogate modeling purposes, but is
general enough that it could be useful to others. I found no other existing
n-dimensional interpolation library in Python for structured data that includes gradients.

Though I've implemented and tested the functionality and API as I'm currently using it in my own projects, I am certainly open to comments and changes for consideration for inclusion into scipy.interpolate.

Proposed updated docs: https://thearn.github.io/docs/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator

pv · 2017-09-21T08:25:01Z

Quick comment: please don't do unrelated whitespace cleanup in a PR that does something else. I'd also prefer not using UnivariateSpline in the implementation. That is a semi-deprecated interface to FITPACK. If possible, consider using make_interp_spline so that FITPACK is not used at all (if you want to use FITPACK, please use splrep+BSpline or splrep+splev). The issue with FITPACK is that some people may expect that the knots of the spline coincide with grid points --- IIRC this is not guaranteed with FITPACK with its automatic grid selection.

thearn · 2017-09-21T10:50:59Z

Sorry, auto-pep8 plugin issue with the whitespace. I can revert that.

I think make_interp_spline should work fine as well, let me test that on my end. I agree that user expectations would likely be that the knots would match the given data.

UPDATE make_interp_spline works, and its axis kwarg gave vectorization potential for the algorithm. Only odd-order splines are supported by BSpline, so the added list of methods is now slinear, cubic, and quintic. Interpolation tests pass and the graphical results visually match.

I'll update the PR shortly, but I'm suddenly getting a failure for scipy/special/tests/test_basic.py::TestBessel::()::test_yv_cephes_vs_amos that only shows up on Travis-CI in the 2nd test job (the Python 2 config run with TESTMODE=fast COVERAGE= NUMPYSPEC="--pre --upgrade --timeout=60 -f $PRE_WHEELS numpy".

I'm not sure how that's coming about, my changes are only within scipy.interpolate. I'll see if I can reproduce that locally.

…rpolation (#1) * FIX: switched from UnivariateSpline to make_interp_spline for 1D interpolation * FIX: can vectorize inner loop with _bsplines calls * FIX: revert auto-pep8 whitespace edits * FIX: removed FITPACK specific wording and transforms

* FIX: updated docstrings to reflect method changes * FIX: removed unnecessary reshaping when using `make_interp_spline` * FIX: grid is now an unneeded paramater to _evaluate_seperable()

MAINT: refactored spline calls

rgommers · 2017-09-22T09:22:02Z

@thearn that test failure in scipy.special is unrelated, due to a recent change in numpy that wasn't quite correct.

* FIX: Vectorized entire first dimension gradient call

ev-br

Had a cursory look, left a few comments, mostly on cosmetics. Am intending to give it a proper code review, but definitely not in time for 1.2.0, sadly.

Overall I'm a big +1 for the functionality.
Not sure how to deal with gradients: having them for some interpolation kinds and not others feels clunky.
A bigger picture question is whether to separate evaluations into an NdBSpline class (or TensorProductSpline or some such) --- I don't have a firm opinion though.

ev-br · 2018-11-05T22:30:44Z

scipy/interpolate/interpolate.py

+        is used, then if any dimension has fewer points than `k` + 1, an error
+        will be raised. If spline_dim_error=False, then the spline interpolant
+        order will be reduced as needed on a per-dimension basis. Default
+        is True (raise an exception).


This option feels borderline extraneous. Maybe better raise unconditionally if there are not enough points?

ev-br · 2018-11-05T22:32:37Z

scipy/interpolate/interpolate.py


    Methods
    -------
    __call__
+    gradient
+    methods


This method should not be public (API surface).

ev-br · 2018-11-05T22:41:34Z

scipy/interpolate/interpolate.py

+    strategy, and are best used when the fitted data (the points and values
+    arrays) are large relative to the number of points to be interpolated.
+    As such, they are not an efficient choice for some tasks, such as
+    re-sampling of image data.


Not sure I understand the comment about the values of the fitted data being larger than the number of data points. Rephrase?

Also for resampling, the See Also section could list relevant signal/ndimage routines.

ev-br · 2018-11-05T22:42:35Z

scipy/interpolate/interpolate.py

@@ -2319,12 +2335,24 @@ class RegularGridInterpolator(object):
    return an array of `nan` values. Nearest-neighbor interpolation will work
    as usual in this case.

+    The 'slinear', 'cubic', 'quintic' methods are spline-based interpolators.
+    Use of the spline interpolations allows for getting gradient
+    values via the ``gradient`` method.


These could use the .. versionadded:: directive (otherwise it looks like they are available from scipy 0.14.0)

scipy/interpolate/interpolate.py

ev-br · 2018-11-05T22:56:44Z

scipy/interpolate/tests/test_interpolate.py

+        value4 = interp(x)
+
+        # values from different methods should be different
+        assert_raises(AssertionError, assert_equal, value1, value2)


numpy.testing.assert_(value1 != value2)

ev-br · 2018-11-05T22:57:33Z

scipy/interpolate/tests/test_interpolate.py

+        values = values - 2j * values
+        sample = np.asarray([[0.1, 0.1, 1., .9]])
+
+        # spline methods dont support complex values


Typo: either do not or with an apostrophe, don't.

ev-br · 2018-11-05T23:00:33Z

scipy/interpolate/tests/test_interpolate.py

+        result1 = interp(x)
+        gradient1 = interp.gradient(x)
+        result_actual_1 = 0.2630309995970872
+        result_gradient_1 = np.array([0.22505535, -0.46465198, 0.02523666])


What are these values and why are they correct? Maybe this test is trying to do two things at once: check against an independently calculated values (this is good! only please add at least a comment on how to reproduce these values), and somehow validate a random input. I'd say these should be two separate tests.

ev-br · 2018-11-05T23:01:39Z

scipy/interpolate/tests/test_interpolate.py

+        result_gradient_1 = np.array([0.22505535, -0.46465198, 0.02523666])
+
+        assert_almost_equal(result1, result_actual_1)
+        assert_almost_equal(gradient1, result_gradient_1)


Bikeshedding: maybe use assert_allclose with appropriate atol and rtol?

ev-br · 2018-11-05T23:04:06Z

scipy/interpolate/tests/test_interpolate.py

-                       [1, 3.3, 1.2, 4.0, 5.0, 1.0, 3]]).T
+        points, z, f, df = self._sample_2d_data_smooth()
+        x, y = points
+        np.random.seed(3)


A very minor nitpick, optional: a recommended idiom is to create a RandomState instance instead of seeding the global generator. We do not really follow this, and there is no need for a code churn on this account alone, but in new code maybe better follow the recommendation.

ev-br · 2019-05-20T21:18:01Z

Needs a rebase

caniko · 2020-11-13T14:37:45Z

This PR has been silent for some time now, any ETA on the new methods?

rgommers · 2020-11-14T18:54:50Z

This PR has been silent for some time now, any ETA on the new methods?

No, it has stalled - if you'd be interested in helping move it forward, that would be useful.

caniko · 2020-11-19T11:04:15Z

This PR has been silent for some time now, any ETA on the new methods?

No, it has stalled - if you'd be interested in helping move it forward, that would be useful.

I'll take a look, if help still is needed, some time in the future after I am done with my current projects.

ev-br · 2022-05-01T14:18:00Z

Now that gh-15357 is in, we can carry on looking at further enhancements from this PR. As I see, there are two major things, which ~~could~~should probably be done in two separate PRs:

what to do when the number of points along a dimension is smaller then the interpolation degree plus one. ISTM a binary switch kwarg is suboptimal, so the options are to either always raise a hard error (this is what ENH: interpolate: add new methods for RegularGridInterpolator. #15357 does), or adjust the spline degree along the dimension. The latter is similar to what e.g. pchip does in 1D for just two points for a cubic polynomial.
Gradients. This is huge! Maybe the right API is to mirror what minimize(..., jac=True) does: if a compute_gradients switch is True, calling the RGI instance returns a tuple of arrays, for the values and gradients. This has a downside of course that the return type of RGI.__call__ depends on the arguments. But this is conceptually simpler and way better performance-wise then a requiring separate gradients call?

Also, from a consistency POV it seems better to have it compute gradients for all methods, with method=linear being piecewise-constant gradients and method=nearest is all zeros. And ISTM gradients only need to be implemented in RGI itself, not interpn.

Thoughts?

ev-br · 2022-05-01T14:25:57Z

@AtsushiSakai @thearn would you be interested in continuing this? I'll gladly review --- or I can take over if you prefer.

AtsushiSakai · 2022-05-02T02:18:45Z

Thank you. I can continue these, but I need time to implement big feature like gradients because I don’t have enough spare time these days. So, It is ok for me that you can take over these. If you have something which has higher priority, you can leave these. I will try these when I have time.

ev-br

This pseudo-review is a way to record of my tracing through this code (which is quite ingenious!), trying to understand the details of what is what.

To avoid rebasing this branch, I simply dropped the full RegularGridInterpolator class from this PR into _rgi.py, renamed it RGI2, and dropped a breakpoint into the __call__ method.

Comments below use the following toy example:

>>> X, Y = np.arange(5), np.arange(6)
>>> values = (X**3 + 2*X) * (Y**3)[:, None]
>>> rgi2 = RGI2((X, Y), values.T, method='cubic')

>>> x, y = 1.5, 0.5     # the evaluation point
>>> rgi2((x, y))

Then the result is (x**3 + 2*x) * y**3 and the gradient is [ (3*x**2 + 2) * y**3, (x**3 + 2*x) * 3*y**2 ]

ev-br · 2022-11-04T05:32:47Z

scipy/interpolate/interpolate.py

+                                                         data_values,
+                                                         xi[:, i],
+                                                         ki[i],
+                                                         compute_gradients)


This line produces first_values = (X**3 + 2*X) * y**3 and first_derivs = (X**3 + 2* X) * 3*y**2 --- the 0th dimension is still a vector, and the second dimension (last in the array, first in the iteration) has y instead of Y.

ev-br · 2022-11-04T05:39:10Z

scipy/interpolate/interpolate.py

+                                              ki,
+                                              compute_gradients=False,
+                                              first_dim_gradient=True)
+            all_gradients[:, -1] = top_grad[-1]


This if clause seems unreachable: the default value of first_dim_gradient is False, and the only way to call it with first_dim_gradient=True is guarded by if first_dim_gradient. Removing this clause completely makes things run fine.

ev-br · 2022-11-04T05:45:25Z

scipy/interpolate/interpolate.py

+                                                               x[i],
+                                                               ki[i],
+                                                               compute_gradients)
+


at this point, values = (X**3 + 2*X) * 0.5**3 and local_derivs = (X**3 + 2*X) * 3*(0.5)**2 : these are 1D arrays where the first dimension is an array, and the second dimension is a scalar, the 2nd element of the input xi array.

ev-br · 2022-11-04T05:49:34Z

scipy/interpolate/interpolate.py

+                                                         ki,
+                                                         compute_gradients=False)
+                else:
+                    gradient[i] = local_derivs


This incantation interpolates the local_derivs --- which is a 1D array of derivatives at x[-1] (this is already computed!) and the X array with respect to the first dimension. The result is evaluated at x[0], which gives the derivative w.r.t. the first dimension at x[0].

ev-br · 2022-11-04T06:26:32Z

So, the current status of porting the gradient computations from this PR: the approach in this PR is quite ingenious indeed, and it can be used in user code today by just copy-pasting the RGI class from this PR, which is pure python. This is probably very slow since

the approach is recursive in the number of dimensions, and
each 1D interpolation constructs a throw-away BSpline object.

And the code get significantly more complicated as compared to the current RGI class which does evaluations but not gradients.

The long-term plan for scipy.interpolate is to finally implement an NdBSpline class to do evaluations of tensor product splines. This would allow to compute derivatives much more easily and efficiently.
Therefore, I think we should hold on porting this, and re-evaluate once NdBSpline is implemented.
(Help with this is appreciated BTW).

I'm going to leave this PR open though, to make it more discoverable and to boost the "we want this" signal.

ev-br · 2024-02-10T16:28:04Z

A belated thank you @thearn for starting this project. While it took us an embarrassing amount of time to act, and not all of your code made it, the idea did strike a chord; the lengths you had to go to work around limitations of then-current scipy.interpolate served as a helpful guidance for improving the library. Thank you!

thearn · 2024-02-11T02:18:46Z

Thank you, @ev-br !

thearn changed the title ~~ENH: add more methods to RegularGridInterpolator.<interpolate>.~~ ENH: add more methods to RegularGridInterpolator. Sep 20, 2017

thearn changed the title ~~ENH: add more methods to RegularGridInterpolator.~~ ENH: add more interpolation methods to RegularGridInterpolator. Sep 20, 2017

FIX: Grid points do not need to be copied

4e64658

thearn added 3 commits September 21, 2017 09:15

FIX: Vectorization and documentation updates

679e78f

* FIX: updated docstrings to reflect method changes * FIX: removed unnecessary reshaping when using `make_interp_spline` * FIX: grid is now an unneeded paramater to _evaluate_seperable()

Feature (#3)

10dade4

MAINT: refactored spline calls

rgommers added enhancement A new feature or improvement scipy.interpolate labels Sep 22, 2017

thearn added 4 commits September 22, 2017 18:44

FIX: Vectorized entire first dimension gradient call

94510c4

* FIX: Vectorized entire first dimension gradient call

simplified

9f1bde8

FIX: Removed unused arguments

319e011

FIX: Vectorized entire first dimension gradient call (#5)

665960d

* FIX: Vectorized entire first dimension gradient call

oesteban mentioned this pull request Aug 9, 2018

Regular Grid, ND, multivariate data B-Spline interpolation #9123

Open

ev-br reviewed Nov 5, 2018

View reviewed changes

pv added the needs-work Items that are pending response from the author label Aug 7, 2019

AtsushiSakai mentioned this pull request Jan 4, 2022

ENH: interpolate: add new methods for RegularGridInterpolator. #15357

Merged

ev-br added this to N-D regular grid in scipy.interpolate Apr 29, 2022

ev-br mentioned this pull request Apr 29, 2022

(2D) Interpolation functions should work with complex numbers #2650

Closed

ev-br reviewed Nov 4, 2022

View reviewed changes

ev-br moved this from N-D regular grid to triaged, waiting-for-contributor in scipy.interpolate Nov 6, 2022

ev-br mentioned this pull request May 20, 2023

Add NdBSpline: n-dim tensor product b-spline object #18492

Merged

ev-br mentioned this pull request Dec 4, 2023

ENH: use NdBSpline in RegularGridInterpolator to speed up evaluations #19633

Merged

tylerjereddy closed this in #19633 Jan 29, 2024

scipy.interpolate automation moved this from triaged, waiting-for-contributor to Abandoned / rejected / supeseded Jan 29, 2024

ev-br removed the needs-work Items that are pending response from the author label Jan 29, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: add more interpolation methods to RegularGridInterpolator. #7903

ENH: add more interpolation methods to RegularGridInterpolator. #7903

thearn commented Sep 20, 2017 •

edited

Loading

thearn commented Sep 20, 2017 •

edited

Loading

pv commented Sep 21, 2017 via email

thearn commented Sep 21, 2017 •

edited

Loading

rgommers commented Sep 22, 2017

ev-br left a comment

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br Nov 5, 2018 •

edited

Loading

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br Nov 5, 2018

ev-br commented May 20, 2019

caniko commented Nov 13, 2020

rgommers commented Nov 14, 2020

caniko commented Nov 19, 2020

ev-br commented May 1, 2022 •

edited

Loading

ev-br commented May 1, 2022

AtsushiSakai commented May 2, 2022

ev-br left a comment •

edited

Loading

ev-br Nov 4, 2022

ev-br Nov 4, 2022

ev-br Nov 4, 2022

ev-br Nov 4, 2022

ev-br commented Nov 4, 2022 •

edited

Loading

ev-br commented Feb 10, 2024

thearn commented Feb 11, 2024

ENH: add more interpolation methods to RegularGridInterpolator. #7903

ENH: add more interpolation methods to RegularGridInterpolator. #7903

Conversation

thearn commented Sep 20, 2017 • edited Loading

thearn commented Sep 20, 2017 • edited Loading

pv commented Sep 21, 2017 via email

thearn commented Sep 21, 2017 • edited Loading

rgommers commented Sep 22, 2017

ev-br left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br Nov 5, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br commented May 20, 2019

caniko commented Nov 13, 2020

rgommers commented Nov 14, 2020

caniko commented Nov 19, 2020

ev-br commented May 1, 2022 • edited Loading

ev-br commented May 1, 2022

AtsushiSakai commented May 2, 2022

ev-br left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br commented Nov 4, 2022 • edited Loading

ev-br commented Feb 10, 2024

thearn commented Feb 11, 2024

thearn commented Sep 20, 2017 •

edited

Loading

thearn commented Sep 20, 2017 •

edited

Loading

thearn commented Sep 21, 2017 •

edited

Loading

ev-br Nov 5, 2018 •

edited

Loading

ev-br commented May 1, 2022 •

edited

Loading

ev-br left a comment •

edited

Loading

ev-br commented Nov 4, 2022 •

edited

Loading