3

I tried to find a float number in ndarray. Due to the software package I am using (Abaqus), the precision it outputs is a little bit low. For example, 10 is something like 10.00003. Therefore, I was wondering whether there is a "correct" way to do it, that is neater than my code.

Example code:

import numpy as np

array = np.arange(10)
number = 5.00001

if I do this:

idx = np.where(number==array)[0][0]

Then the result is empty because 5.00001 does not equal to 5.

Now I am doing:

atol = 1e-3 # Absolute tolerance
idx = np.where(abs(number-array) < atol)[0][0]

which works, and is not too messy... Yet I was wondering there would be a neater way to do it. Thanks!

PS: numpy.allclose() is another way to do it, but I need to use number * np.ones([array.shape[0], array.shape[1]]) and it still seems verbose to me...


Edit: Thank you all so much for the fantastic answers! np.isclose() is the exact function that I am looking for, and I missed it since it is not in the doc... I wouldn't have realized this until they update the doc, if it weren't you guys. Thank you again!

Yuxiang Wang
  • 8,095
  • 13
  • 64
  • 95

3 Answers3

5

PS: numpy.allclose() is another way to do it, but I need to use number * np.ones([array.shape[0], array.shape[1]]) and it still seems verbose to me...

You almost never need to do anything like number * np.ones([array.shape[0], array.shape[1]]). Just as you can multiply that scalar number by that ones array to multiply all of its 1 values by number, you can pass that scalar number to allclose to compare all of the original array's values to number. For example:

>>> a = np.array([[2.000000000001, 2.0000000002], [2.000000000001, 1.999999999]])
>>> np.allclose(a, 2)
True

As a side note, if you really do need an array of all 2s, there's an easier way to do it than multiplying 2 by ones:

>>> np.tile(2, array.shape)
array([[2, 2], [2, 2]])

For that matter, I don't know why you need to do [array.shape[0], array.shape[1]]. If the array is 2D, that's exactly the same thing as array.shape. If the array might be larger, it's exactly the same as array.shape[:2].


I'm not sure this solves your actual problem, because it seems like you want to know which ones are close and not close, rather than just whether or not they all are. But the fact that you said you could use allclose if not for the fact that it's too verbose to create the array to compare with.

So, if you need whereclose rather than allclose… well, there's no such function. But it's pretty easy to build yourself, and you can always wrap it up if you're doing it repeatedly.

If you had an isclose method—like allclose, but returning a bool array instead of a single bool—you could just write:

idx = np.where(isclose(a, b, 0, atol))[0][0]

… or, if you're doing it over and over:

def whereclose(a, b, rtol=1e-05, atol=1e-08):
    return np.where(isclose(a, b, rtol, atol))

idx = whereclose(a, b, 0, atol)[0][0]

As it turns out, version 1.7 of numpy does have exactly that function (see also here), but it doesn't appear to be in the docs. If you don't want to rely on a possibly-undocumented function, or need to work with numpy 1.6, you can write it yourself trivially:

def isclose(a, b, rtol=1e-05, atol=1e-08):
    return np.abs(a-b) <= (atol + rtol * np.abs(b))
abarnert
  • 354,177
  • 51
  • 601
  • 671
  • 1
    You probably should drop the `np.` from `np.whereclose`. By the way, the `close` function you've defined is almost identical to the existing `np.isclose` except that the `np.isclose` will automatically broadcast the array shapes. – askewchan Sep 09 '13 at 21:22
  • @askewchan: Thanks on the first part. As for the second: Yes, that was the whole point. Well, also, the fact that `isclose` isn't in numpy 1.6, and isn't in the docs for 1.7 despite apparently being in the code. Although, now that you make me think about it… `close` is a really, really misleading name to use, and I'm not sure what's better… – abarnert Sep 09 '13 at 21:37
  • Yeah I was confused by your comment at first then I tried googling :P For others' reference: [`np.isclose`](https://github.com/numpy/numpy/blob/6c729b4423857850e6553cf6c2d0fc8b026036dd/numpy/core/numeric.py#L2133) – askewchan Sep 09 '13 at 21:52
  • @askewchan: Edited the answer to show how to use `isclose` if present, or how to replace it if not. – abarnert Sep 09 '13 at 23:54
  • Thank you so much @abarnet! I really appreciate your help. np.close() is exactly what I am looking for; and also, I didn't realize np.tile before, and nor did I realize how dummy I was in writing ndarray.shape[0] and ndarray.shape[1]... Thanks so much! – Yuxiang Wang Sep 10 '13 at 20:12
3

If you have up-to-date numpy (1.7), then the best way is to use np.isclose which will broadcast the shapes together automatically:

import numpy as np
a = np.arange(10)
n = 5.000001
np.isclose(a, n).nonzero()
#(array([5]),)

or, if you expect only one match:

np.isclose(a, n).nonzero()[0][0]
#5

(np.nonzero is basically the same thing as np.where except that it doesn't have the if condition then/else capability)

askewchan
  • 45,161
  • 17
  • 118
  • 134
0

The method you use above, specifically abs(A - B) < atol, is standard for doing floating point comparisons across many languages. Obviously when using numpy A and/or B can be arrays or numbers.

Here is another approach that might be useful to look at. I'm not sure it applies to your case, but it could be very helpful if you're looking for more than one number in the array (which is a common use case). It's inspired by this question which is kind of similar.

import numpy as np

def find_close(a, b, rtol=1e-05, atol=1e-08):
    tol = atol + abs(b) * rtol
    lo = b - tol
    hi = b + tol
    order = a.argsort()
    a_sorted = a[order]
    left = a_sorted.searchsorted(lo)
    right = a_sorted.searchsorted(hi, 'right')
    return [order[L:R] for L, R in zip(left, right)]

a = np.array([2., 3., 3., 4., 0., 1.])
b = np.array([1.01, 3.01, 100.01])
print find_close(a, b, atol=.1)
# [array([5]), array([1, 2]), array([], dtype=int64)]
Community
  • 1
  • 1
Bi Rico
  • 25,283
  • 3
  • 52
  • 75