Working with floating point NumPy arrays for comparison and related operations
<p>I have an array of random floats and I need to compare it to another one that has the same values in a different order. For that matter I use the sum, product (and other combinations depending on the dimension of the table hence the number of equations needed).</p>

<p>Nevertheless, I encountered a precision issue when I perform the sum (or product) on the array depending on the order of the values. </p>

<p>Here is a simple standalone example to illustrate this issue :</p>

<pre><code>import numpy as np

n = 10
m = 4

tag = np.random.rand(n, m)

s1 = np.sum(tag, axis=1)
s2 = np.sum(tag[:, ::-1], axis=1)

# print the number of times s1 is not equal to s2 (should be 0)
print np.nonzero(s1 != s2)[0].shape[0]
</code></pre>

<p>If you execute this code it sometimes tells you that <code>s1</code> and <code>s2</code> are not equal and the differents is of magnitude of the computer precision.</p>

<p>The problem is I need to use those in functions like <code>np.in1d</code> where I can't really give a tolerance...</p>

<p>Is there a way to avoid this issue?</p>

<p>For the listed code, you can use <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.isclose.html" rel="nofollow"><code>np.isclose</code></a> and with it tolerance values could be specified too.</p>

<p>Using the provided sample, let's see how it could be used -</p>

<pre><code>In [201]: n = 10
     ...: m = 4
     ...: 
     ...: tag = np.random.rand(n, m)
     ...: 
     ...: s1 = np.sum(tag, axis=1)
     ...: s2 = np.sum(tag[:, ::-1], axis=1)
     ...: 

In [202]: np.nonzero(s1 != s2)[0].shape[0]
Out[202]: 4

In [203]: (~np.isclose(s1,s2)).sum() # So, all matches!
Out[203]: 0
</code></pre>

<p>To make use of tolerance values in other scenarios, we need to work on a case-by-case basis. So, let's say for an implementation that involve elementwise comparison like in <code>np.in1d</code>, we can bring in <a href="http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html" rel="nofollow"><code>broadcasting</code></a> to do those elementwise equality checks for all elems in first input against all elems in the second one. Then, we use <code>np.abs</code> to get the "closeness factor" and finally compare against the input tolerance to decide the matches. As needed to simulate <code>np.in1d</code>, we do ANY operation along one of the axis. Thus, <code>np.in1d</code> with tolerance using <code>broadcasting</code> could be implemented like so -</p>

<pre><code>def in1d_with_tolerance(A,B,tol=1e-05):
    return (np.abs(A[:,None] - B) &lt; tol).any(1)
</code></pre>

<p>As suggested in the comments by OP, we can also round floating-pt numbers after scaling them up and this should be memory efficient, as being needed for working with large arrays. So, a modified version would be like so -</p>

<pre><code>def in1d_with_tolerance_v2(A,B,tol=1e-05):
    S = round(1/tol)
    return np.in1d(np.around(A*S).astype(int),np.around(B*S).astype(int))
</code></pre>

<p>Sample run -</p>

<pre><code>In [372]: A = np.random.rand(5)
     ...: B = np.random.rand(7)
     ...: B[3] = A[1] + 0.0000008
     ...: B[6] = A[4] - 0.0000007
     ...: 

In [373]: np.in1d(A,B) # Not the result we want!
Out[373]: array([False, False, False, False, False], dtype=bool)

In [374]: in1d_with_tolerance(A,B)
Out[374]: array([False,  True, False, False,  True], dtype=bool)

In [375]: in1d_with_tolerance_v2(A,B)
Out[375]: array([False,  True, False, False,  True], dtype=bool)
</code></pre>

<p>Finally, on how to make it work for other implementations and use cases - It would depend on the implementation itself. But for most cases, <code>np.isclose</code> and <code>broadcasting</code> should help.</p>

python arrays numpy random precision