Skip to content

Implemented np.allclose in numba/np/arraymath.py #7940

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sklam merged 16 commits into from
Apr 21, 2022
Merged

Implemented np.allclose in numba/np/arraymath.py #7940

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
e209bc1
Added np_allclose in numba/np/arraymath.py
czgdp1807 Mar 26, 2022
d16cd63
Added test for np_allclose
czgdp1807 Mar 26, 2022
eee11d0
Added broadcasting support and stronger tests
czgdp1807 Apr 1, 2022
2dc5bb7
Handeled scalars and improved tests
czgdp1807 Apr 2, 2022
ca98b63
Added test data from np.allclose tests
czgdp1807 Apr 2, 2022
6542301
Applied linting
czgdp1807 Apr 2, 2022
f09d8ff
Call a separate function to ensure same rank of arrays
czgdp1807 Apr 2, 2022
50367ef
Applied linting
czgdp1807 Apr 2, 2022
0cc9234
Merge branch 'main' into np_allclose
czgdp1807 Apr 2, 2022
9ee9842
Use np.broadcast_arrays and remove utility functions
czgdp1807 Apr 2, 2022
4d90833
Separate numpy tests for allclose
czgdp1807 Apr 7, 2022
ead6817
Improved tests of np_allclose
czgdp1807 Apr 8, 2022
9c1d44f
Corrected error messages in np_allclose
czgdp1807 Apr 8, 2022
1b24788
Added numpy.allclose to numpysupported.rst
czgdp1807 Apr 9, 2022
a93b0ef
Merge branch 'main' into np_allclose
czgdp1807 Apr 20, 2022
66969ce
Corrected argument number in error message
czgdp1807 Apr 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions docs/source/reference/numpysupported.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@ The following methods of Numpy arrays are supported:

* :meth:`~numpy.ndarray.transpose`
* :meth:`~numpy.ndarray.view` (only the 1-argument form)
* :meth:`~numpy.ndarray.__contains__`
* :meth:`~numpy.ndarray.__contains__`

Where applicable, the corresponding top-level NumPy functions (such as
:func:`numpy.argmax`) are similarly supported.
Expand Down Expand Up @@ -416,6 +416,7 @@ Other functions

The following top-level functions are supported:

* :func:`numpy.allclose`
* :func:`numpy.append`
* :func:`numpy.arange`
* :func:`numpy.argsort` (``kind`` key word argument supported for values
Expand Down Expand Up @@ -655,7 +656,7 @@ Permutations
Distributions
'''''''''''''

The following functions support all arguments.
The following functions support all arguments.

* :func:`numpy.random.beta`
* :func:`numpy.random.binomial`
Expand Down
94 changes: 94 additions & 0 deletions numba/np/arraymath.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -864,6 +864,100 @@ def flat_all(a):
return flat_all


@register_jitable
def _allclose_scalars(a_v, b_v, rtol=1e-05, atol=1e-08, equal_nan=False):
a_v_isnan = np.isnan(a_v)
b_v_isnan = np.isnan(b_v)

# only one of the values is NaN and the
# other is not.
if ( (not a_v_isnan and b_v_isnan) or
(a_v_isnan and not b_v_isnan) ):
return False

# either both of the values are NaN
# or both are numbers
if a_v_isnan and b_v_isnan:
if not equal_nan:
return False
else:
if np.isinf(a_v) or np.isinf(b_v):
return a_v == b_v

if np.abs(a_v - b_v) > atol + rtol * np.abs(b_v * 1.0):
return False

return True


@overload(np.allclose)
@overload_method(types.Array, "allclose")
def np_allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False):

if not type_can_asarray(a):
raise TypeError('The first argument "a" must be array-like')

if not type_can_asarray(b):
raise TypeError('The first argument "b" must be array-like')

@sklam sklam Apr 14, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
raise TypeError('The first argument "b" must be array-like')
raise TypeError('The second argument "b" must be array-like')

@czgdp1807 czgdp1807 Apr 20, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have addressed the review and the tests are passing. Please let me know if anything else is to be done here. Thanks.


if not isinstance(rtol, types.Float):
raise TypeError('The third argument "rtol" must be a '
'floating point')

if not isinstance(atol, types.Float):
raise TypingError('The fourth argument "atol" must be a '
'floating point')

if not isinstance(equal_nan, types.Boolean):
raise TypeError('The fifth argument "equal_nan" must be a '
'boolean')

is_a_scalar = isinstance(a, types.Number)
is_b_scalar = isinstance(b, types.Number)

if is_a_scalar and is_b_scalar:
def np_allclose_impl_scalar_scalar(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
return _allclose_scalars(a, b, rtol=rtol, atol=atol,
equal_nan=equal_nan)
return np_allclose_impl_scalar_scalar
elif is_a_scalar and not is_b_scalar:
def np_allclose_impl_scalar_array(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
b = np.asarray(b)
for bv in np.nditer(b):
if not _allclose_scalars(a, bv.item(), rtol=rtol, atol=atol,
equal_nan=equal_nan):
return False
return True
return np_allclose_impl_scalar_array
Comment on lines +924 to +933

@guilhermeleobas guilhermeleobas Apr 6, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can replace this block of code by:

Suggested change
elif is_a_scalar and not is_b_scalar:
def np_allclose_impl_scalar_array(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
b = np.asarray(b)
for bv in np.nditer(b):
if not _allclose_scalars(a, bv.item(), rtol=rtol, atol=atol,
equal_nan=equal_nan):
return False
return True
return np_allclose_impl_scalar_array
elif is_a_scalar and not is_b_scalar:
def np_allclose_impl_scalar_array(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
return np.allclose(b, a, rtol=rtol, atol=atol, equal_nan=equal_nan)
return np_allclose_impl_scalar_array

@czgdp1807 czgdp1807 Apr 7, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this will give right results all the time, because the check depends on the ordering of a and b. For (a, b) the check is abs(a - b) > atol + rtol * abs(b), where as for (b, a) the check is abs(b - a) > atol + rtol * abs(a). See the RHS in both the inequalities, one uses b and the other uses a. So, np.allclose(a, b) might give different result from np.allclose(b, a) for edge cases.

Quoting from "Notes" section of https://numpy.org/doc/stable/reference/generated/numpy.allclose.html

The above equation is not symmetric in a and b, so that allclose(a, b) might be different from allclose(b, a) in some rare cases.

Please let me know if should still apply the above code suggestion.

@guilhermeleobas guilhermeleobas Apr 7, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fine. I wasn't aware of this info.

elif not is_a_scalar and is_b_scalar:
def np_allclose_impl_array_scalar(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
a = np.asarray(a)
for av in np.nditer(a):
if not _allclose_scalars(av.item(), b, rtol=rtol, atol=atol,
equal_nan=equal_nan):
return False
return True
return np_allclose_impl_array_scalar
elif not is_a_scalar and not is_b_scalar:
def np_allclose_impl_array_array(a, b, rtol=1e-05, atol=1e-08,
equal_nan=False):
a = np.asarray(a)
b = np.asarray(b)
a_a, b_b = np.broadcast_arrays(a, b)

for av, bv in np.nditer((a_a, b_b)):
if not _allclose_scalars(av.item(), bv.item(), rtol=rtol,
atol=atol, equal_nan=equal_nan):
return False

return True

return np_allclose_impl_array_array


@overload(np.any)
@overload_method(types.Array, "any")
def np_any(a):
Expand Down
182 changes: 182 additions & 0 deletions numba/tests/test_np_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -329,6 +329,10 @@ def np_trapz_x_dx(y, x, dx):
return np.trapz(y, x, dx)


def np_allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False):
return np.allclose(a, b, rtol, atol, equal_nan)


def np_average(a, axis=None, weights=None):
return np.average(a, axis=axis, weights=weights)

Expand Down Expand Up @@ -3172,6 +3176,184 @@ def test_1D_weights_axis(data, axis, weights):
#test with axis argument
test_1D_weights_axis(data, axis=1, weights=w)

def test_allclose(self):

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

min_int = np.iinfo(np.int_).min
a = np.array([min_int], dtype=np.int_)

simple_data = [
(np.asarray([1e10, 1e-7]), np.asarray([1.00001e10, 1e-8])),
(np.asarray([1e10, 1e-8]), np.asarray([1.00001e10, 1e-9])),
(np.asarray([1e10, 1e-8]), np.asarray([1.0001e10, 1e-9])),
(np.asarray([1e10]), np.asarray([1.0001e10, 1e-9])),
(1.0, 1.0),
(np.array([np.inf, 1]), np.array([0, np.inf])),
(a, a)
]

for a, b in simple_data:
py_result = pyfunc(a, b)
c_result = cfunc(a, b)
self.assertEqual(py_result, c_result)

a = np.asarray([1.0, np.nan])
b = np.asarray([1.0, np.nan])
self.assertFalse(cfunc(a, b))
self.assertEquals(pyfunc(a, b, equal_nan=True),
cfunc(a, b, equal_nan=True))

b = np.asarray([np.nan, 1.0])
self.assertEquals(pyfunc(a, b), cfunc(a, b))

noise_levels = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 0.0]
zero_array = np.zeros((25, 4))
a = np.random.ranf((25, 4))
for noise in noise_levels:
for rtol in noise_levels:
for atol in noise_levels:
py_result = pyfunc(zero_array, noise,
atol=atol, rtol=rtol)
c_result = cfunc(zero_array, noise,
atol=atol, rtol=rtol)
self.assertEqual(py_result, c_result)

py_result = pyfunc(noise, zero_array,
atol=atol, rtol=rtol)
c_result = cfunc(noise, zero_array,
atol=atol, rtol=rtol)
self.assertEqual(py_result, c_result)

py_result = pyfunc(np.asarray([noise]), zero_array,
atol=atol, rtol=rtol)
c_result = cfunc(np.asarray([noise]), zero_array,
atol=atol, rtol=rtol)
self.assertEqual(py_result, c_result)

py_result = pyfunc(a, a + noise, atol=atol, rtol=rtol)
c_result = cfunc(a, a + noise, atol=atol, rtol=rtol)
self.assertEqual(py_result, c_result)

py_result = pyfunc(a + noise, a, atol=atol, rtol=rtol)
c_result = cfunc(a + noise, a, atol=atol, rtol=rtol)
self.assertEqual(py_result, c_result)

def test_ip_allclose_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2402-L2420 # noqa: E501
pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

arr = np.array([100.0, 1000.0])
aran = np.arange(125).astype(dtype=np.float64).reshape((5, 5, 5))

atol = 1e-8
rtol = 1e-5

numpy_data = [
(np.asarray([1, 0]), np.asarray([1, 0])),
(np.asarray([atol]), np.asarray([0.0])),
(np.asarray([1.0]), np.asarray([1 + rtol + atol])),
(arr, arr + arr * rtol),
(arr, arr + arr * rtol + atol * 2),
(aran, aran + aran * rtol),
(np.inf, np.inf),
(np.inf, np.asarray([np.inf]))
]

for (x, y) in numpy_data:
self.assertEquals(pyfunc(x, y), cfunc(x, y))

def test_ip_not_allclose_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2422-L2441 # noqa: E501

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

aran = np.arange(125).astype(dtype=np.float64).reshape((5, 5, 5))

atol = 1e-8
rtol = 1e-5

numpy_data = [
(np.asarray([np.inf, 0]), np.asarray([1.0, np.inf])),
(np.asarray([np.inf, 0]), np.asarray([1.0, 0])),
(np.asarray([np.inf, np.inf]), np.asarray([1.0, np.inf])),
(np.asarray([np.inf, np.inf]), np.asarray([1.0, 0.0])),
(np.asarray([-np.inf, 0.0]), np.asarray([np.inf, 0.0])),
(np.asarray([np.nan, 0.0]), np.asarray([np.nan, 0.0])),
(np.asarray([atol * 2]), np.asarray([0.0])),
(np.asarray([1.0]), np.asarray([1 + rtol + atol * 2])),
(aran, aran + aran * atol + atol * 2),
(np.array([np.inf, 1.0]), np.array([0.0, np.inf]))
]

for (x, y) in numpy_data:
self.assertEquals(pyfunc(x, y), cfunc(x, y))

def test_return_class_is_ndarray_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2460-L2468 # noqa: E501

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

class Foo(np.ndarray):
def __new__(cls, *args, **kwargs):
return np.array(*args, **kwargs).view(cls)

a = Foo([1])
self.assertTrue(type(cfunc(a, a)) is bool)

def test_equalnan_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2456-L2458 # noqa: E501
pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

x = np.array([1.0, np.nan])

self.assertEquals(pyfunc(x, x, equal_nan=True),
cfunc(x, x, equal_nan=True))

def test_no_parameter_modification_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2443-L2448 # noqa: E501

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

x = np.array([np.inf, 1])
y = np.array([0, np.inf])

cfunc(x, y)
np.testing.assert_array_equal(x, np.array([np.inf, 1]))
np.testing.assert_array_equal(y, np.array([0, np.inf]))

def test_min_int_numpy(self):
# https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/core/tests/test_numeric.py#L2450-L2454 # noqa: E501

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

min_int = np.iinfo(np.int_).min
a = np.array([min_int], dtype=np.int_)

self.assertEquals(pyfunc(a, a), cfunc(a, a))

def test_allclose_exception(self):
self.disable_leak_check()

pyfunc = np_allclose
cfunc = jit(nopython=True)(pyfunc)

a = np.asarray([1e10, 1e-9, np.nan])
b = np.asarray([1.0001e10, 1e-9])

with self.assertRaises(ValueError) as e:
cfunc(a, b)

self.assertIn(("shape mismatch: objects cannot be broadcast to "
"a single shape"), str(e.exception))

def test_interp_basic(self):
pyfunc = interp
cfunc = jit(nopython=True)(pyfunc)
Expand Down