NEP-18 not working

[leofang]

I am testing the compatibility with CuPy and I notice this error:

import cupy as cp
import numpy as np
import autocorr


N = 10240
np.random.seed(0)
t = np.arange(N)
A = np.exp(-0.05 * t)[:, np.newaxis] + np.random.rand(N, 24) * 0.1

A_cp = cp.asarray(A)
g1, tau = autocorr.multitau(A_cp.T, 16)

Traceback (most recent call last):
  File "test_autocorr_cupy.py", line 13, in <module>
    g1, tau = autocorr.multitau(A_cp.T, 16)
  File "/home/leofang/autocorr/autocorr/multitau.py", line 57, in multitau
    g2[:, i] = np.mean(a[:, :N - i] * a[:, i:], axis=1) / t1 / t2
ValueError: object __array__ method not producing an array

The reasons is that internally multitau allocates two arrays g2 and tau using NumPy API, which is not (yet) smart enough to dispatch to CuPy.

I tried to refresh my memory on NEP-18, but I don't think this issue was addressed there. Looks like to accommodate different NEP-18 compliant arrays, we need to detect the array source and call the corresponding array-creation API?

if isinstance(input, numpy.ndarray):
    g2 = numpy.empty(...)
elif isinstance(input, cupy.ndarray):
    g2 = cupy.empty(...)
elif  # dask? sparse?

I may have missed something obvious, though.

[tacaswell]

can we walk back up the class / module tree of the input to sort out what version of empty to call?

[leofang]

Looks like we need to rely on __module__? Not really a fan of it...

>>> import inspect
>>>
>>> import cupy as cp
>>> a = cp.empty(10)
>>> type(a).__module__
'cupy.core.core'
>>> inspect.getmodule(a) is None  # useless
True
>>>
>>> import numpy as np
>>> b = np.empty(10)
>>> type(b).__module__
'numpy'
>>> inspect.getmodule(b) is None  # useless
True

[leofang]

For CuPy's uninformative __module__, just raised an issue: cupy/cupy#3266. Not sure whether other libraries are compliant on this aspect or not.

[danielballan]

I'm late to the party here, just noticed this. The way to handle array creation is use np.*_like functions, as in

numpy.empty_like(some_input_whose_type_will_be_imitated, dtype=..., shape=...)

In the past, the *_like functions did not allow users to customize the shape because the whole point was to return array of the same shape. Now, they may be used to return an array of the same type but, optionally, a different shape or internal dtype. This was a added in a recent release of numpy; the PR to add it to dask is not yet merged as of this writing: dask/dask#6064