ENH: Add derivatives for BarycentricInterpolator #18197

jcwhitehead · 2023-03-25T14:25:04Z

What does this implement/fix?

Implements derivatives for BarycentricInterpolator.

I've extended the existing tests for the KroghInterpolator derivatives to also test the new methods, and added a few additional tests.

Reference issue

This should close #2725.

Additional information

I haven't contributed to scipy before – let me know if I've missed anything! Happy to implement changes before it can be merged.

…ut and Trefethen (cf scipy#2725)

j-bowhay

Thanks @jcwhitehead, didn't have time to check the maths in any detail but a few suggestions to start

scipy/interpolate/_polyint.py

scipy/interpolate/tests/test_polyint.py

scipy/interpolate/_polyint.py

tupui

Thank you @jcwhitehead. Some additional suggestions to help with the review.

(Thanks @j-bowhay! A few things for you to note 😉)

scipy/interpolate/_polyint.py

jcwhitehead · 2023-03-28T17:58:25Z

I just removed @ev-br as a requested reviewer by accident - apologies (I don't seem to have a way to undo this).

I think I've addressed all the points raised but if there's any more feedback I can keep iterating further. Please do just 'unresolve' a conversation I've marked as resolved if you disagree and think there's more to be done.

tupui · 2023-03-28T18:03:25Z

I just removed @ev-br as a requested reviewer by accident - apologies (I don't seem to have a way to undo this).

No worries, this has no impact for us and how we approach PR reviews. (We don't use these features that much.)

I think I've addressed all the points raised but if there's any more feedback I can keep iterating further. Please do just 'unresolve' a conversation I've marked as resolved if you disagree and think there's more to be done.

In general please don't resolve discussions at all. This is important for maintainers to know where they are at in their review. You can add a quick comment or some thumbs up if you want on specific discussions to indicate you took care of it. Or make a general message listing what you did and what's left to be done.

scipy/interpolate/_polyint.py

tupui · 2023-03-28T18:09:54Z

scipy/interpolate/_polyint.py

+        x : array_like
+            1D array of points at which to evaluate the derivatives
+        der : integer, optional
+            The number of derivatives to evaluate, from 'order 0' (der=1) to order der-1.


Our linter is a bit lax at the moment, but we still would like to have everything under 80 chars (79 max).

tupui · 2023-03-28T18:12:42Z

scipy/interpolate/_polyint.py

@@ -665,8 +738,8 @@ def _evaluate(self, x):
            c = x[..., np.newaxis] - self.xi
            z = c == 0
            c[z] = 1
-            c = self.wi/c
-            with np.errstate(divide='ignore'):
+            c = self.wi / c


Not 100% sure about this one, but noting this here. Please refrain from updating anything else than what your work requires. Unless of course the linter is unhappy and ask you to do it.

tupui · 2023-03-28T18:13:20Z

scipy/interpolate/_polyint.py

-            c = self.wi/c
-            with np.errstate(divide='ignore'):
+            c = self.wi / c
+            with np.errstate(divide='ignore', invalid='ignore'):


Where is the warning coming from? We try to avoid things like that and if we do, it's best to document

Thanks for spotting this - I think this is probably now unnecessary after I added the ValueError in __init__ to catch the special (nonsensical) case in which nodes coincide, so I will remove it.

tupui · 2023-03-28T18:16:59Z

scipy/interpolate/_polyint.py

-    def __init__(self, xi, yi=None, axis=0):
-        _Interpolator1D.__init__(self, xi, yi, axis)
+    def __init__(self, xi, yi=None, axis=0, *, wi=None):
+        _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)


Suggested change

_Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)

super().__init__(self, xi, yi, axis)

Should the corresponding line in KroghInterpolator also be updated to use super() or is this kind of style change only for new code?

scipy/interpolate/_polyint.py

…more specific with match='...'

jcwhitehead · 2023-04-01T08:45:15Z

Those latest three CI failures don't look like they're related to anything that's changed in the commit:

FAILED ../../stats/tests/test_continuous_basic.py::test_cont_basic[500-200-geninvgauss-arg34]
FAILED ../../stats/tests/test_continuous_basic.py::test_methods_with_lists[geninvgauss-args34-ppf]
FAILED ../../stats/tests/test_distributions.py::TestGenInvGauss::test_rvs_p_zero

From my side I think everything raised so far has now been addressed - if not or if there's more feedback I'm happy to keep iterating to get it approved and merged.

j-bowhay · 2023-04-01T17:48:51Z

Thanks @jcwhitehead this is on my list to have another look at but it might be a couple of weeks before I get enough time

j-bowhay

A few minor changes but otherwise I think this is looking good. See https://numpydoc.readthedocs.io/en/latest/format.html#common-rest-concepts for when to use single vs double ticks.

Finally, as this is a new feature please could you write a post to the mailing list inviting others to take a look.

scipy/interpolate/_polyint.py

j-bowhay · 2023-04-10T15:57:41Z

scipy/interpolate/_polyint.py

+            # fill in correct diagonal entries: each column sums to 0
+            np.fill_diagonal(c, 0)
+            # c[j,j] = -sum_{i != j} c[i,j]
+            d = -c.sum(axis=1)


Equation numbers of the original paper would be great to add actually.

scipy/interpolate/_polyint.py

j-bowhay · 2023-06-15T05:46:08Z

Friendly ping @jcwhitehead, are you interested in finishing this up? This is looking in good shape with only some minor changes required so it would be great to see this in the next SciPy release

[skip ci]

j-bowhay · 2023-07-21T08:47:46Z

@ev-br @Kai-Striega @AtsushiSakai I have pushed some changes to hopefully get this over the line. Since I have made some edits would one of you be able to have a look over this

ev-br

LGTM modulo minor nits.

scipy/interpolate/_polyint.py

ev-br · 2023-07-21T14:05:23Z

scipy/interpolate/_polyint.py

+            c = self.xi[..., np.newaxis] - self.xi
+            # c[i,j] = (w[j] / w[i]) / (xi[i] - xi[j])
+            c = np.divide(self.wi, c * self.wi[..., np.newaxis],
+                          out=np.zeros_like(c), where=c != 0)


Not sure what is going on here. out=np.zeros_like(c) signals "throw away the result", but then it is returned, ouch. It is probably deliberately ingenious (which it is!), but can we please write it in a more boring way?

This was one of the changes I made to avoid rewriting the diagonal back a forth (as c is 0 along the diagonal). Here the division is performed normally apart from where c==0 where the result is instead left unchanged in the out array (hence is zero as desired).

Would you be happy with a comment explaining what is happening?

Unfortunately I haven't had any time to keep up with this lately so I'm out of the loop, but I'd actually be a bit surprised if the new version with its quadratic comparisons was faster than the previous version with two linear fills.

If it's a meaningful difference it's probably worth propagating to the older barycentric interpolation code which uses the same two diagonal fills.

@j-bowhay I guess you've checked the new version is faster - do you have some numbers? I'd be especially interested in the high derivatives, huge xi case

I have to admit I thought I benchmarked this but clearly not. There is no applicable difference so I will revert the change:

import numpy as np import perfplot def setup(n): xi = np.random.random(n) c = xi - xi[:, np.newaxis] np.fill_diagonal(c, 0) return c def old(x): c = np.copy(x) np.fill_diagonal(c, 1) c = 1/c np.fill_diagonal(c,0) return c def new(x): c = np.copy(x) return np.divide(1, c, out=np.zeros_like(c), where=c!=0) perfplot.show(setup=setup, kernels=[old, new], n_range=[2**k for k in range(15)])

I'm fine with just adding a comment, FTR.

ev-br · 2023-07-21T14:09:28Z

scipy/interpolate/_polyint.py

+
+        if self._diff_cij is None:
+            # c[i,j] = xi[i] - xi[j]
+            c = self.xi[..., np.newaxis] - self.xi


Nit: The comment hints that self.xi is 1D, is this always correct? If it is, I'd ask to consider writing a more explicit xi[:, np.newaxis].

Rationale: an Ellipsis signals an unspecified number of multiple dimensions. So if it really possibly multidim, could you please make it more explicit in the comment line above.

Ah yes good spot self.xi is the x coordinates of the points the polynomial should pass through so is 1d

ev-br · 2023-07-21T14:12:06Z

scipy/interpolate/_polyint.py

+            # fill in correct diagonal entries: each column sums to 0
+            np.fill_diagonal(c, 0)
+            # c[j,j] = -sum_{i != j} c[i,j]
+            d = -c.sum(axis=1)


Equation numbers of the original paper would be great to add actually.

j-bowhay · 2023-08-01T05:16:58Z

Thanks @ev-br @jcwhitehead, I have hopefully addressed all comments including reverting the use of np.divide and adding a random_state argument

ev-br · 2023-08-11T10:36:56Z

This now LGTM. Anything you want to add/tweak @jcwhitehead @j-bowhay ?

One minor hiccup is 11 commits looking a bit too much here, would be nice to squash some, while preserving authorship.

scipy/interpolate/_polyint.py

[skip actions] [skip cirrus]

j-bowhay · 2023-08-11T18:39:06Z

I had one final look and all happy from my end.

I am happy just to squash merge if that works for everyone (my contributions were mostly just housekeeping)

ev-br · 2023-08-13T15:46:41Z

OK, let's merge this then if everyone's happy. Thank you @jcwhitehead , @j-bowhay !
Fun fact: the original issue this PR closed was opened exactly one day short of ten years ago :-). Great to see it fixed!

P.S. Turns out GH adds a co-authored-by stanza to commit messages on squash-merges.

j-bowhay · 2023-08-13T15:48:27Z

Thanks @jcwhitehead for getting the ball rolling on this, it was a great first contribution. I hope we can tempt you back for more when you have the time!

Implemented derivatives for BarycentricInterpolator according to Berr…

c562194

…ut and Trefethen (cf scipy#2725)

jcwhitehead requested a review from ev-br as a code owner March 25, 2023 14:25

j-bowhay added enhancement A new feature or improvement scipy.interpolate labels Mar 25, 2023

jcwhitehead added 3 commits March 25, 2023 17:18

Make docstring a raw string to fix docs error

277d8c9

Fix docstring import error

2b45707

Remove duplicated line in docstring

3303239

j-bowhay requested changes Mar 25, 2023

View reviewed changes

tupui reviewed Mar 26, 2023

View reviewed changes

jcwhitehead added 2 commits March 27, 2023 22:42

Updated tests, docstrings, signatures etc in line with PR feedback

5315f78

Fix docstring indentation

9758ec8

jcwhitehead requested review from j-bowhay and removed request for ev-br March 28, 2023 17:53

tupui reviewed Mar 28, 2023

View reviewed changes

tupui requested a review from ev-br March 28, 2023 18:20

Split long lines, remove errstate='invalid', make test_repeated_node …

d501320

…more specific with match='...'

j-bowhay added this to In review in scipy.interpolate Apr 6, 2023

j-bowhay added this to the 1.11.0 milestone Apr 6, 2023

j-bowhay requested changes Apr 10, 2023

View reviewed changes

j-bowhay reviewed Apr 20, 2023

View reviewed changes

scipy/interpolate/_polyint.py Outdated Show resolved Hide resolved

ev-br added the needs-work Items that are pending response from the author label May 5, 2023

h-vetinari modified the milestones: 1.11.0, 1.12.0 May 20, 2023

j-bowhay and others added 2 commits July 21, 2023 09:03

DOC: fix formatting

a933122

[skip ci]

Merge remote-tracking branch 'upstream/main' into jcwhitehead/main

cb5773b

use np.divide to avoid rewriting diagonal

64de284

ev-br reviewed Jul 21, 2023

View reviewed changes

j-bowhay added 6 commits July 31, 2023 08:52

MAINT: make explicit calculation involves 1d arrays

59af078

MAINT: remove unseeded randomness

c30ecce

DOC: add equation numbers

1073789

revert c calculations to fill diagonal

336a3f5

DOC: use seed in example

38cc844

fix doc build

14feb82

j-bowhay reviewed Aug 11, 2023

View reviewed changes

scipy/interpolate/_polyint.py Outdated Show resolved Hide resolved

DOC: add missing backtics

bfe4d54

[skip actions] [skip cirrus]

ev-br merged commit 5ebc7af into scipy:main Aug 13, 2023
6 checks passed

scipy.interpolate automation moved this from In review to Done Aug 13, 2023

ev-br removed the needs-work Items that are pending response from the author label Aug 13, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: Add derivatives for BarycentricInterpolator #18197

ENH: Add derivatives for BarycentricInterpolator #18197

jcwhitehead commented Mar 25, 2023 •

edited

j-bowhay left a comment

tupui left a comment

jcwhitehead commented Mar 28, 2023 •

edited

tupui commented Mar 28, 2023

tupui Mar 28, 2023

tupui Mar 28, 2023

tupui Mar 28, 2023

jcwhitehead Mar 28, 2023

tupui Mar 28, 2023

jcwhitehead Mar 28, 2023

jcwhitehead commented Apr 1, 2023 •

edited

j-bowhay commented Apr 1, 2023

j-bowhay left a comment

j-bowhay Apr 10, 2023

ev-br Jul 21, 2023

j-bowhay commented Jun 15, 2023

j-bowhay commented Jul 21, 2023

ev-br left a comment

ev-br Jul 21, 2023

j-bowhay Jul 31, 2023

jcwhitehead Jul 31, 2023 •

edited

j-bowhay Jul 31, 2023 •

edited

ev-br Jul 31, 2023

ev-br Jul 21, 2023

j-bowhay Jul 31, 2023

ev-br Jul 21, 2023

j-bowhay commented Aug 1, 2023

ev-br commented Aug 11, 2023

j-bowhay commented Aug 11, 2023

ev-br commented Aug 13, 2023

j-bowhay commented Aug 13, 2023

	_Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)
	super().__init__(self, xi, yi, axis)

ENH: Add derivatives for BarycentricInterpolator #18197

ENH: Add derivatives for BarycentricInterpolator #18197

Conversation

jcwhitehead commented Mar 25, 2023 • edited

What does this implement/fix?

Reference issue

Additional information

j-bowhay left a comment

Choose a reason for hiding this comment

tupui left a comment

Choose a reason for hiding this comment

jcwhitehead commented Mar 28, 2023 • edited

tupui commented Mar 28, 2023

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

jcwhitehead commented Apr 1, 2023 • edited

j-bowhay commented Apr 1, 2023

j-bowhay left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

j-bowhay commented Jun 15, 2023

j-bowhay commented Jul 21, 2023

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

jcwhitehead Jul 31, 2023 • edited

Choose a reason for hiding this comment

j-bowhay Jul 31, 2023 • edited

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

j-bowhay commented Aug 1, 2023

ev-br commented Aug 11, 2023

j-bowhay commented Aug 11, 2023

ev-br commented Aug 13, 2023

j-bowhay commented Aug 13, 2023

jcwhitehead commented Mar 25, 2023 •

edited

jcwhitehead commented Mar 28, 2023 •

edited

jcwhitehead commented Apr 1, 2023 •

edited

jcwhitehead Jul 31, 2023 •

edited

j-bowhay Jul 31, 2023 •

edited