The next three commits refactor a bit further but show how we could do the dispatching with a decorator. Even though we'd have to add decorated functions to the tests manually, it's a lot simpler overall because this approach wouldn't require changes to __init__.py. You wouldn't even need a _support_alternative_backends.py file for each subpackage, and generic implementations can be placed right next to the original function.

The latest commit resolves #21012 (comment) for the case in which array arguments are specified as consecutive positional or keyword arguments, which is the most common case. data-apis/array-api-compat#153 is an issue if the arguments are specified in an unusual way.

At least one of the property-based array API tests will fail because det returns a float64 given float32 NumPy arrays. Looks like this is documented, so an exception can be made in the test.

Thanks @lucascolley. Do you want all functions to support check_finite, or should I just work on emitting a nicer error if the user passes arguments not supported by the array API extension?

Let's start with just emitting nice errors. We should figure out how to start supporting other keywords in an agnostic way, though. Can we somehow only support check_finite if it is in the function's signature to start with?

The existing SciPy test suite serves as a contract to what the the SciPy functionality actually is. If we are to dispatch to other backends, we should IMO keep the contract and be explicit of what we xfail/skip.

I think that it would be a good idea to convert the tests eventually¹, but while the behaviour is experimental and opt-in via an env variable, surely we are fine to punt on it for now. Any users who are going out of their way to help us test this experimental dispatching should be smart enough to send bug reports to the appropriate library, rather than to us.

For the dispatch itself, I've two high-level comments first.

1. The current approach seems to involve populating a submodule-size dict per submodule, is this correct?

How about using the fact that a module is a namespace itself, and do something along the lines of

- submodule/__init__.py:

from _support_alternative_backends import ... # explicit names, not import *
__all__ = _submodule_api.__all__


- submodule/_support_alternative_backends.py:

from ._submodule_api import * 

# wrap/dispatch on all imported names here


- submodule/_submodule_api.py:

from _private_modules import *    # or explicit names

# populate the __all__ list

2. The dispatch approach still seems to rely on either a number of leading arguments (or kwargs?). I don't believe it scales to all scipy APIs. How does, e.g. https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.convolve1d.html#scipy.ndimage.convolve1d fits?
   Note the output optional which comes after axis.

I think that it would be a good idea to convert the tests eventually1, but while the behaviour is experimental and opt-in via an env variable, surely we are fine to punt on it for now.

I disagree.

I disagree.

Okay, let's split this discussion out to gh-21031.

1. The current approach seems to involve populating a submodule-size dict per submodule, is this correct?

IIUC you're referring to the array_special_func_map and similar, and you're proposing that we don't need to define such a map, we can just use __all__? I think there are a few reasons to define the map: we don't necessarily want to dispatch everything (some specific attention is desirable), and we need to provide a little information about the function's signature (unless we want to do some fancy introspection).

However, this map is just what I did for special. I also show in this PR how we could use a decorator instead of a such a map. If we go with the decorator approach, we can eliminate _support_alternative_backend.py from each submodule.

2. The dispatch approach still seems to rely on either a number of leading arguments (or kwargs?). I don't believe it scales to all scipy APIs.

No, if we're not allowed to make any modifications, it will not work for all APIs. But "I think this will work for all linalg functions. I found a handful of signal functions that don't fit this mold, and that's fine - we just need to look closer to know whether it's worth generalizing or whether these can be treated ad hoc." (#21012 (comment).). We can move quickly through the low-hanging fruit, identify the functions that don't fit the mold, and make the generalizations as needed. If we find that we need to add support for functions like your convolve_dispatcher, etc., to handle the special cases, we can certainly do that. This just proposes that we don't need all the boilerplate for simple functions.

As for convolve1d in particular, I'm not sure. It will not pass that argument to array_namespace, so perhaps the user would get a different error message if, say, the output is not a CuPy array but the first two arguments are. But I think otherwise it will work for CuPy, which seems to be the only other backend with a direct convolve1d equivalent. For other backends, one would have to write custom code anyway if we want to support the output argument.

So for convolve1d it works in the sense that:

import cupy as xp
from scipy.ndimage import convolve1d
input = [2, 8, 0, 4, 1, 9, 9, 0]
weights = [1, 3]
output = xp.empty(8, dtype=int)
convolve1d(xp.asarray(input), weights=xp.asarray(weights), output=output)
print(output, type(output))
# [14 24  4 13 12 36 27  0] <class 'cupy.ndarray'>

However, if the user were to pass an invalid argument for output, they would get whatever error CuPy provides, not e.g. "TypeError: Multiple namespaces for array inputs:".

An example of a function that wouldn't work for all arguments as written is sosfilt, where axis separates two array arguments. Currently, array_namespace produces an error if one of the arguments is an int. This already causes problems when we're translating functions and arguments can be either Python scalars or arrays, but if it were not changed, then passing the axis argument as an int would make array_namespace raise an error, and passing axis as a 0d CuPy array makes CuPy raise an error since it doesn't interpret this as an int. One can imagine variations to the code that would get arround this, e.g. having array_namespace ignore Python scalars as it ignores None, or having a way of skipping arguments to be passed to array_namespace.

For the sake of argument, I can also imagine examples of signatures in stats in which weights is accepted as a keyword-only argument. (We will actually be translating these functions manually, not dispatching, but there might be similar cases in other modules that we'd want to dispatch.)

• If the alternative backend has the desired function, the only problem is that a weights array would not be passed into array_namespace for validation.
• If the alternative backend does not have the desired function and we need to fall back to the SciPy implementation with a NumPy array, the problem is that the current code would not convert weights to a NumPy array before passing it into the SciPy function. So if the SciPy function itself does not call np.asarray on weights and the conversion doesn't happen implicitly, there could be a problem.

Again, if we find a large class of functions like this that we want to dispatch, we either ensure that the SciPy versions of the functions convert to NumPy arrays, or can generalize the code to handle keyword-only array arguments.

The current approach seems to involve populating a submodule-size dict per submodule, is this correct?

IIUC you're referring to the array_special_func_map and similar, and you're proposing that we don't need to define such a map, we can just use all? I think there are a few reasons to define the map: we don't necessarily want to dispatch everything (some specific attention is desirable), and we need to provide a little information about the function's signature (unless we want to do some fancy introspection).

Here's what I meant: #21091
This can be probably simplified quite a bit, but the main idea is to have a namespace instead of a dict.
Maybe we can avoid messing with sys.modules and use vanilla imports even?

We can move quickly through the low-hanging fruit, identify the functions that don't fit the mold, and make the generalizations as needed.

gh-21091 has the whole of ndimage, and it's explicit which ones need what. Not sure how to simplify it further.

Here's what I meant

OK, so it looks like you're using _ndimage_api.__all__ plus scipy/ndimage/_dispatchers.py instead of a dictionary defined in _support_alternative_backends. Sure, that's fine. The other approach shown here was to apply a decorator directly to the function, no dictionary or namespace file necessary. There are pros/cons of each.

...has the whole of ndimage

OK, since the first argument of all of those functions is an array, I think that functions wrapped per this PR would probably pass all the same tests with CuPy as those wrapped by gh-21091. The difference is just that if the user were to pass arrays with different backends, this PR would not raise an error immediately in array_namespace.

IIUC, gh-21091 does not convert other CPU arrays (e.g. torch CPU) to NumPy, run the SciPy function, and convert back. Is that right?