Back to the simple point example of comparing two shapes for common points where there might be floating point issues. Note that in the b1 array, I introduces a really tiny difference in the second point so that it differs from its previous incarnation.... b
a = np.array([[0., 0.], [1., 4.], [4., 3.], [5., 0.],[0., 0.]])
b = np.array([[1., 0.], [1., 4.], [5., 0.], [1., 0.]])
b1 = np.array([[1., 0.], [1., 4.000001], [5., 0.], [1., 0.]])
So if we compare a to b for equality, then we would have an issue because of the small floating point difference. To account for this, we need to step up the check to account for this.
Here is the code. Read the header.
def a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False):
r"""Return indices of, or points where, two point arrays are equal.
Parameters
----------
a, b : 2d arrays
No error checking so ensure that a and b are 2d, but their shape need
not be the same.
Notes
-----
Modified from `np.isclose`, stripping out the nan and ma checks.
Adjust atol and rtol appropriately.
>>> np.nonzero(a_eq_b(a, b))[0]
One could use::
>>> np.equal(A, B).all(1).nonzero()[0]
>>> np.where((a[:, None] == b).all(-1).any(1))[0]
If you aren't worried about floating-point issues in equality checks.
"""
a = np.atleast_2d(a)
b = np.atleast_2d(b)
if b.size > 2:
b = b[:, None]
w = np.less_equal(np.abs(a - b), atol + rtol * abs(b)).all(-1)
if w.ndim > 1:
w = w.any(-1).squeeze()
if return_pnts:
return b[w].squeeze()
return w.squeeze()
The magic is to account for tolerances (covered in the NumPy documentation). This allows one to identify points that are ridiculously close that they should be considered the same. If the difference is due to errors, then they can be fixed.
So compare a to b and a to b1 with and without tolerance checks.
# -- a to b comparison, both are correct
a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False)
array([0, 1, 1, 0])
a_eq_b(a, b, atol=0., rtol=0., return_pnts=False)
array([0, 1, 1, 0])
# -- a to b1 comparison, tolerance check is correct
a_eq_b(a, b1, atol=1.0e-8, rtol=1.0e-5, return_pnts=False)
Out[160]: array([0, 1, 1, 0])
# -- exactitude check only identifies one 'equality'
a_eq_b(a, b1, atol=0., rtol=0., return_pnts=False)
Out[161]: array([0, 0, 1, 0])
So when you things are all good, make sure that you account for the small differences.
If coordinates/values are the result of a calculation, it is really important to round/truncate your results to the expected precision that is associated with the rest of your data.
Enough for now.
Addendum
see Common points for a fuller discussion
You must be a registered user to add a comment. If you've already registered, sign in. Otherwise, register and sign in.