Add np.allclose and np.isclose support ref. issue #4074

[jeertmans]

Hi everyone,

As proposed in issue #4074, I added support for the np.allclose function.

I also wrote unit test. One of them could be enhanced (about checking data type) but support to numpy.issubdtype(a, b) should be first brought. Anyway, I think that the function works quite as expected.

If I made any mistake or forgot anything, don't hesitate :)

If you think the implementation is correct, I can easily implement the numpy.isclose function.

I timed the function and it seems to work pretty well:

a = np.random.rand(100, 100)
np.allclose(a, a)
>>> True
allclose(a, a)
>>> True
%timeit allclose(a, a)
>>> 10.9 µs ± 31.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit np.allclose(a, a)
>>> 59.4 µs ± 221 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

[esc]

@jeertmans thanks for submitting this! I have scheduled it for review.

[esc]

@jeertmans so, we discussed this in the core-dev meeting yesterday and I think the implementation may be overly complex. I don't think that much of the error checking will be needed as much of the functionality may already be present? Probably the core of the algorithm can be done using the vectorized:

return False if any(abs(a - b) > (atol + rtol * abs(b))) else return True

Where both a and b are arrays.

With perhaps a special case for Nan? Or did I miss something?

[esc]

I'm just looking into this now so I'll leave messages here as I come across things.

For example, I am not not convinced we need an extra _broadcastable_, because using subtraction on arrays will yield it:

In [13]: from numba import njit

In [14]: @njit
    ...: def sub(a, b):
    ...:     return (a-b)
    ...:

In [15]: sub(a, b)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-0a247c27eacb> in <module>
----> 1 sub(a, b)

ValueError: unable to broadcast argument 1 to output array
File "/Users/vhaenel/git/numba/numba/np/npyimpl.py", line 228,

Granted, the error message isn't the same as numpy:

In [16]: a - b
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-16-09bd029d0285> in <module>
----> 1 a - b

ValueError: operands could not be broadcast together with shapes (4,3) (3,5)

But that is probably something that should be fixed elsewhere in Numba.

[jeertmans]

Thank you @esc !
I will take a deeper look at it and I keep you up to date.
About vectorizing, I did not know the any() function was implemented, my bad :)

For the special case of NaN, I still have to think about how to do this.

[esc]

Another thought is, to look at the relevant Numpy functions:

• https://github.com/numpy/numpy/blob/master/numpy/core/numeric.py#L2256
• https://github.com/numpy/numpy/blob/master/numpy/core/numeric.py#L2344

As far as I can tell allclose is implemented using isclose and isclose seems a little more involved. 😁

[jeertmans]

So @esc , I looked a bit into it and it seems that using vectorized function will just be slower than the current implementation.
Even without taking into account the special case of NaNs, using np.all(abs(a - b) > (atol + rtol * abs(b))) takes much more time. Same happens if a swap np.all to np.any.

The differences in time are shown:

a = np.random.rand(1000, 1000)
b = a  # All values are the same
...
(numbaenv) jerome@jerome-N501VW:~/Documents/programmation/numba$ python3 ../allclose_perf.py
>>> 0.23322463035583496 (using vectorised op.)
>>> 0.15993094444274902 (my implementation)
>>> 0.9217619895935059 (np.allclose)
b = -a  # All values are different
...
(numbaenv) jerome@jerome-N501VW:~/Documents/programmation/numba$ python3 ../allclose_perf.py
>>> 0.2327868938446045 (using vectorised op.)
>>> 3.337860107421875e-06 (my implementation)
>>> 0.9791421890258789 (np.allclose)

From my point of view, the problem is that the np.any (resp. the np.all) function is not implemented "smartly".
The np.allclose implementation uses np.iclose function. But we actually don't need this O(n) memory allocated by this function. It is actually worse that O(n) because, in order to handle NaN, they allocate several arrays for this purpose.

So, if you want to purely replicate Numpy's behavior, I can easily do this, but I think it will come with a slowdown in performances.

About the _broadcastable_ function, I agree that it is maybe not the best solution, but that is the only one I saw that was without any useless memory allocation, which Numpy does.

[esc]

@jeertmans excellent, thank you for looking into this with such detail. I'll give it some thought and get back to you!

[esc]

@jeertmans I gave this some thought and discussed with some of the other developers. We believe that correctness is more important than performance, initially. That is to say, performance does of course matter greatly as Numba is all about performance, but the implementation needs to be 100% correct (i.e. matching Numpy behaviour) first. If you would like to add code that diverges from the Numpy approach, that is of course encouraged, especially if it improves performance, however the tests must be very rigorous and thorough in such cases and demonstrate that correctness was not sacrificed. As a start, the _broadcastable_ function will need to be extensively tested and covered by test cases. Also, it is probably worth looking into why np.all(abs(a - b) > (atol + rtol * abs(b))) is slower and by how much. Having a vectorized form is significantly beneficial in many cases as it helps to avoid temporaries and lends itself to be being parallelized with Numba's parallel accelerator (aka parfors).

However, my recommendation in this case would actually be to take step back, look at the is_close operation of Numpy and then use that as a basis for allclose. It seems to have a significant portion of code dedicated to dealing with specific edge cases, such handling of masked arrays, Nan and inf values which I don't see in this PR at present. In any case, thank you very much for looking into this, we do appreciate your efforts! This is a long standing issue and it would be awesome to see it resolved!

[jeertmans]

@esc Thank you for you complete answer ! I will give it a moment and will reach you back soon =)

[jeertmans]

So @esc , I have implemented the function according to Numpy's source code.

It works but I have some small problems:

1. numpy.asanyarray and other functions are not implemented in Numba, so that the implementation is not a perfect mirror of Numpy's. I don't see when it should cause problem but I guess a good idea should be to implement those functions, but it is maybe to much for a single PR.
2. Numpy's way to handle 0d arrays is not compliant with Numba's way: most Numpy functions accept scalar elements and return scalar. But Numba doesn't seem to like working with both scalars and arrays:

>>> Rejected as the implementation raised a specific error:
>>>     TypingError: Failed in nopython mode pipeline (step: nopython frontend)
>>> Can't unify return type from the following types: array(bool, 1d, C), bool

For the moment, I decided to used numpy.atleast_1d function to overcome this problem. The problem with this is that the output will not have the same format as numpy.isclose function.

Thus, I also wrote similar tests for isclose (and updated PR title).

Now, performances:

a = np.random.rand(1000, 1000)
%timeit allclose(a,a)
>>> 7.49 ms ± 54.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit np.allclose(a,a)
>>> 9.27 ms ± 90.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
b = - a
%timeit allclose(a,b)
>>> 7.23 ms ± 469 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit np.allclose(a,b)
>>> 10.2 ms ± 572 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

I am of course open to any critic or suggestion :)

[esc]

@jeertmans great! Excited to see this coming along. I have marked it as "waiting on reviewer" for now. Until we get around to taking a closer look, perhaps you could fix the remaining flake8 errors in the meantime:

numba/np/arraymath.py:4495:81: E501 line too long (83 > 80 characters)
numba/np/arraymath.py:4516:81: E501 line too long (82 > 80 characters)

[jeertmans]

@esc Thanks ! I have fixed it (I guess, waiting now for the results from Flake8). I also made a slight modification but nothing drastic :)

[jeertmans]

Well there seems to be another error in the checks but I really don't think my code does not pass checks, see the errors logs:

# Install latest llvmlite build
$CONDA_INSTALL -c numba llvmlite
\dirname "${CONDA_EXE}"
\dirname "${SYSP}"
Collecting package metadata (current_repodata.json): ...working... failed

CondaHTTPError: HTTP 504 GATEWAY TIME-OUT for url <https://conda.anaconda.org/numba/linux-64/current_repodata.json>
Elapsed: 01:00.130693
CF-RAY: 5db791664bde9fe8-IAD

A remote server error occurred when trying to retrieve this URL.

A 500-type error (e.g. 500, 501, 502, 503, etc.) indicates the server failed to
fulfill a valid request.  The problem may be spurious, and will resolve itself if you
try your request again.  If the problem persists, consider notifying the maintainer
of the remote server.

(https://dev.azure.com/numba/numba/_build/results?buildId=6740&view=logs&j=f2076665-76cb-5d7c-2ee9-55e9ea449f59&t=3c2cd5ad-7adf-505c-38f1-d4511566d1f6&l=3865)

If I'm wrong, can you help me finding what is the problem ? :)

[esc]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 1 pipeline(s).

[esc]

@jeertmans that is an error/glitch from anaconda.org -- I've re-queued the tests now 🤞

[jeertmans]

Thanks @esc !

[jeertmans]

Hello @stuartarchibald , you added this PR to the "PR Backlog" milestone, but I was wondering what it was and it could not find any information about it on the internet :/
Can you help me ? :)

[stuartarchibald]

Hello @stuartarchibald , you added this PR to the "PR Backlog" milestone, but I was wondering what it was and it could not find any information about it on the internet :/
Can you help me ? :)

Hi @jeertmans, I can try and explain. "PR Backlog" is a Numba specific milestone we use to help manage the project. By default all pull requests which are not yet targeted at a specific release but are in progress are added to this milestone. If a PR is critical for a particular release it will get added to a specific milestone, for all other PRs, they are added to a release milestone when they get near completion and the release milestone to which they are added depends on where in the release cycle they happen to land, how risky the code change is, if there are any dependencies etc. Hope this explanation helps?

[jeertmans]

@stuartarchibald Yes it helped ! Thank you :)

[esc]

@jeertmans just as a heads up, we are now in the burndown phase for the next release candidate, 0.52.0RC1, which essentially means, we will have to put reviewing this on hold, until after the release. Hopefully we can resume reviewing this soon.

[jeertmans]

@esc thank you for keeping me up to date ! Go luck with the new release ;)
This PR can wait for sure

[gmarkall]

@esc What hall we do about re-reviewing this? I guess it's a bit close to 0.54 to look at, but maybe could go in the 0.55 milestone to remind us to take a look at it before then?