ENH: integrate.trapezoid: add array API standard support

[Armavica]

Reference issue

Towards gh-20930.

What does this implement/fix?

Implement the Array API for integrate.trapezoid.

Additional information

Hello! I met @lucascolley at the NumFocus summit this week, and he advised me to start with integrate.trapezoid to help with the Array API implementation effort. I am submitting this as a draft PR because I don't understand yet everything that is going on here, and I am looking for feedback and advice.

I have several questions so far:

• There is a masked array test, which does not work with Array API. I understand that supporting masked arrays is out of scope of Array API. How should we then proceed with this test?
• I remember that Lucas told me that np.diff would probably need to be reimplemented, so I don't really understand why my initial strategy of blindly replacing np with xp seems to pass some tests.

[lucascolley]

thanks for working on this @Armavica !

There is a masked array test, which does not work with Array API. I understand that supporting masked arrays is out of scope of Array API. How should we then proceed with this test?

You can decorate the test with skip_xp_invalid_arg from scipy.conftest - see https://scipy.github.io/devdocs/dev/api-dev/array_api.html#adding-tests.

I remember that Lucas told me that np.diff would probably need to be reimplemented, so I don't really understand why my initial strategy of blindly replacing np with xp seems to pass some tests.

Let’s see if we get any failures when you push the converted tests :)

[Armavica]

Thank you for your comments and for the pointers to documentation! I addressed your comments, it looks like tests are passing :)

[lucascolley]

We also need to convert the tests here as described on https://scipy.github.io/devdocs/dev/api-dev/array_api.html#adding-tests. At the moment, all of the tests still only input NumPy arrays, which is why xp.diff and xp.asanyarray aren't failing.

[Armavica]

Ah, that's why! I had the feeling something was still wrong despite tests passing. Thank you, I get it now.

[Armavica]

I fixed the array API commit, and added a commit to fix a bug which was already present in the original version, in cases where 1 < x.ndim ≤ axis < y.ndim.

[lucascolley]

thanks, getting close!

[lucascolley]

Could you add this test file to the XP_TESTS list in https://github.com/scipy/scipy/blob/main/.github/workflows/array_api.yml?

[Armavica]

I added the line to the file as requested.

It looks like this change automatically triggered review requests on my behalf?

[lucascolley]

yep, no worries, that is https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-code-owners

[lucascolley]

thanks @Armavica ! I pushed a small update (explained inline) and I'm happy with this now. Will leave open for a while in case anyone has comments about the bug-fix!

[Armavica]

Thank you for your help on this @lucascolley!

[Armavica]

Note that another approach for what I fixed would be to do like in simpson:

scipy/scipy/integrate/_quadrature.py

Lines 459 to 461 in 8e4d71a

	elif len(x.shape) != len(y.shape):
	raise ValueError("If given, shape of x must be 1-D or the "
	"same as y.")

[lucascolley]

Actually, one more thing! Could you add an array-like test, along the lines of the tests added in gh-21209? It doesn't have to be complicated, just a sanity check that we are still accepting Python lists.

[Armavica]

I added the test but I am probably missing something because I don't think the test runs. At least on my machine, it never does (even with -b numpy or -b all).

[lucascolley]

I think it is just over-indented! It's currently sitting inside test_masked.

[Armavica]

Gosh that's embarrassing… 😅 It's fixed and the test seems to run now.

[lucascolley]

LGTM, thanks again!

[lucascolley]

Note that another approach for what I fixed would be to do like in simpson

What was the behaviour before your fix, by the way? Just nonsense returned?

[mdhaber]

@Armavica thanks for working on this.

@lucascolley I saw that you requested my review. I took a close look and didn't find much wrong with the Array API translation, but mostly because I was distracted by the function's nonstandard broadcasting rules. They don't follow either those of NumPy linalg or scipy.stats (which differ only for positive axis when the arrays don't have the same number of dimensions), and they aren't documented. I'd suggest that as a follow-up we emit a FutureWarning for input arrays that are not of the same shape. Afterwards, this function could be greatly simplified by following standard broadcasting rules.

The one thing I did notice is that this preserved dtype before but doesn't now. This is because it uses sum, and in the 2022.12 of the array API standard (which is still in effect in array_api_compat), sum doesn't maintain dtype. Until 2023.12 takes effect, I'd suggest taking the result_type and passing that explicitly to sum.

[Armavica]

What was the behaviour before your fix, by the way? Just nonsense returned?

No, an indexing error.

[Armavica]

The one thing I did notice is that this preserved dtype before but doesn't now. This is because it uses sum, and in the 2022.12 of the array API standard (which is still in effect in array_api_compat), sum doesn't maintain dtype. Until 2023.12 takes effect, I'd suggest taking the result_type and passing that explicitly to sum.

@mdhaber I added this. It turned out to break the previous version of test_array_like (and the masked array test too, apparently), because the input had an integer dtype, and it looks like xp.sum makes rounding errors in this case (returning 20 instead of the correct answer 22). So I also changed the test from y = [t*t for t in x] to y = [float(t*t) for t in x].
Perhaps another option would be to let the sum return a float, that we could only then convert to the original dtype, resulting in possibly fewer rounding errors even on integer inputs?

[mdhaber]

Thanks; I wasn't explicit enough. Please pass the arrays into _array_api.xp_broadcast_promote (gh-21528) with force_floating=True to determine the result type. (Unfortunately, you cannot just use the broadcasted arrays that are returned because trapezoid does not follow normal broadcasting rules; you'll need to apply their dtype to the input arrays with xp.astype(..., copy=False).)

[lucascolley]

CI failures look unrelated - thanks again @Armavica !

[lucascolley]

looks like this x.ndim == 1 guard is what has gone missing for gh-21908.