can I do fast set difference with floats using numpy if elements are equal up to some tolerance

Question

I have two lists of float numbers, and I want to calculate the set difference between them.

With numpy I originally wrote the following code:

aprows = allpoints.view([('',allpoints.dtype)]*allpoints.shape[1])
rprows = toberemovedpoints.view([('',toberemovedpoints.dtype)]*toberemovedpoints.shape[1])
diff = setdiff1d(aprows, rprows).view(allpoints.dtype).reshape(-1, 2)

This works well for things like integers. In case of 2d points with float coordinates that are the result of some geometrical calculations, there's a problem of finite precision and rounding errors causing the set difference to miss some equalities. For now I resorted to the much, much slower:

diff = []
for a in allpoints: 
    remove = False
    for p in toberemovedpoints:
        if norm(p-a) < 0.1:
            remove = True
    if not remove:
        diff.append(a)
return array(diff)

But is there a way to write this with numpy and gain back the speed?

Note that I want the remaining points to still have their full precision, so first rounding the numbers and then do a set difference probably is not the way forward (or is it? :) )

Edited to add an solution based on scipy.KDTree that seems to work:

def remove_points_fast(allpoints, toberemovedpoints):
    diff = []
    removed = 0
    # prepare a KDTree
    from scipy.spatial import KDTree
    tree = KDTree(toberemovedpoints,  leafsize=allpoints.shape[0]+1)
    for p in allpoints:
        distance, ndx = tree.query([p], k=1)
        if distance < 0.1:
            removed += 1
        else:
            diff.append(p)
    return array(diff), removed

Answer 1

If you want to do this with the matrix form, you have a lot of memory consumption with larger arrays. If that does not matter, then you get the difference matrix by:

diff_array = allpoints[:,None] - toberemovedpoints[None,:]

The resulting array has as many rows as there are points in allpoints, and as many columns as there are points in toberemovedpoints. Then you can manipulate this any way you want (e.g. calculate the absolute value), which gives you a boolean array. To find which rows have any hits (absolute difference < .1), use numpy.any:

hits = numpy.any(numpy.abs(diff_array) < .1, axis=1)

Now you have a vector which has the same number of items as there were rows in the difference array. You can use that vector to index all points (negation because we wanted the non-matching points):

return allpoints[-hits]

This is a numpyish way of doing this. But, as I said above, it takes a lot of memory.

---

If you have larger data, then you are better off doing it point by point. Something like this:

return allpoints[-numpy.array([numpy.any(numpy.abs(a-toberemoved) < .1) for a in allpoints ])]

This should perform well in most cases, and the memory use is much lower than with the matrix solution. (For stylistic reasons you may want to use numpy.all instead of numpy.any and turn the comparison around to get rid of the negation.)

(Beware, there may be pritning mistakes in the code.)