Python - Intersection of 2D Numpy Arrays

Question

I'm desperately searching for an efficient way to check if two 2D numpy Arrays intersect.

So what I have is two arrays with an arbitrary amount of 2D arrays like:

A=np.array([[2,3,4],[5,6,7],[8,9,10]])
B=np.array([[5,6,7],[1,3,4]])
C=np.array([[1,2,3],[6,6,7],[10,8,9]])

All I need is a True if there is at least one vector intersecting with another one of the other array, otherwise a false. So it should give results like this:

f(A,B)  -> True
f(A,C)  -> False

I'm kind of new to python and at first I wrote my program with Python lists, which works but of course is very inefficient. The Program takes days to finish so I am working on a numpy.array solution now, but these arrays really are not so easy to handle.

Here's Some Context about my Program and the Python List Solution:

What i'm doing is something like a self-avoiding random walk in 3 Dimensions. http://en.wikipedia.org/wiki/Self-avoiding_walk. But instead of doing a Random walk and hoping that it will reach a desirable length (e.g. i want chains build up of 1000 beads) without reaching a dead end i do the following:

I create a "flat" Chain with the desired length N:

X=[]
for i in range(0,N+1):
    X.append((i,0,0))

Now i fold this flat chain:

1. randomly choose one of the elements ("pivotelement")
2. randomly choose one direction ( either all elements to the left or to the right of the pivotelment)
3. randomly choose one out of 9 possible rotations in space (3 axes * 3 possible rotations 90°,180°,270°)
4. rotate all the elements of the chosen direction with the chosen rotation
5. Check if the new elements of the chosen direction intersect with the other direction
6. No intersection -> accept the new configuration, else -> keep the old chain.

Steps 1.-6. have to be done a large amount of times (e.g. for a chain of length 1000, ~5000 Times) so these steps have to be done efficiently. My List-based solution for this is the following:

def PivotFold(chain):
randPiv=random.randint(1,N)  #Chooses a random pivotelement, N is the Chainlength
Pivot=chain[randPiv]  #get that pivotelement
C=[]  #C is going to be a shifted copy of the chain
intersect=False
for j in range (0,N+1):   # Here i shift the hole chain to get the pivotelement to the origin, so i can use simple rotations around the origin
    C.append((chain[j][0]-Pivot[0],chain[j][1]-Pivot[1],chain[j][2]-Pivot[2]))
rotRand=random.randint(1,18)  # rotRand is used to choose a direction and a Rotation (2 possible direction * 9 rotations = 18 possibilitys)
#Rotations around Z-Axis
if rotRand==1:
    for j in range (randPiv,N+1):
        C[j]=(-C[j][1],C[j][0],C[j][2])
        if C[0:randPiv].__contains__(C[j])==True:
            intersect=True
            break
elif rotRand==2:
    for j in range (randPiv,N+1):
        C[j]=(C[j][1],-C[j][0],C[j][2])
        if C[0:randPiv].__contains__(C[j])==True:
            intersect=True
            break
...etc
if intersect==False: # return C if there was no intersection in C
    Shizz=C
else:
    Shizz=chain
return Shizz

The Function PivotFold(chain) will be used on the initially flat chain X a large amount of times. it's pretty naivly written so maybe you have some protips to improve this ^^ I thought that numpyarrays would be good because i can efficiently shift and rotate entire chains without looping over all the elements ...

Solution 1

This should do it:

In [11]:

def f(arrA, arrB):
    return not set(map(tuple, arrA)).isdisjoint(map(tuple, arrB))
In [12]:

f(A, B)
Out[12]:
True
In [13]:

f(A, C)
Out[13]:
False
In [14]:

f(B, C)
Out[14]:
False

To find intersection? OK, set sounds like a logical choice. But numpy.array or list are not hashable? OK, convert them to tuple. That is the idea.

A numpy way of doing involves very unreadable boardcasting:

In [34]:

(A[...,np.newaxis]==B[...,np.newaxis].T).all(1)
Out[34]:
array([[False, False],
       [ True, False],
       [False, False]], dtype=bool)
In [36]:

(A[...,np.newaxis]==B[...,np.newaxis].T).all(1).any()
Out[36]:
True

Some timeit result:

In [38]:
#Dan's method
%timeit set_comp(A,B)
10000 loops, best of 3: 34.1 µs per loop
In [39]:
#Avoiding lambda will speed things up
%timeit f(A,B)
10000 loops, best of 3: 23.8 µs per loop
In [40]:
#numpy way probably will be slow, unless the size of the array is very big (my guess)
%timeit (A[...,np.newaxis]==B[...,np.newaxis].T).all(1).any()
10000 loops, best of 3: 49.8 µs per loop

Also the numpy method will be RAM hungry, as A[...,np.newaxis]==B[...,np.newaxis].T step creates a 3D array.

Solution 2

Using the same idea outlined here, you could do the following:

def make_1d_view(a):
    a = np.ascontiguousarray(a)
    dt = np.dtype((np.void, a.dtype.itemsize * a.shape[1]))
    return a.view(dt).ravel()

def f(a, b):
    return len(np.intersect1d(make_1d_view(A), make_1d_view(b))) != 0

>>> f(A, B)
True
>>> f(A, C)
False

This doesn't work for floating point types (it will not consider +0.0 and -0.0 the same value), and np.intersect1d uses sorting, so it is has linearithmic, not linear, performance. You may be able to squeeze some performance by replicating the source of np.intersect1d in your code, and instead of checking the length of the return array, calling np.any on the boolean indexing array.

Solution 3

You can also get the job done with some np.tile and np.swapaxes business!

def intersect2d(X, Y):
        """
        Function to find intersection of two 2D arrays.
        Returns index of rows in X that are common to Y.
        """
        X = np.tile(X[:,:,None], (1, 1, Y.shape[0]) )
        Y = np.swapaxes(Y[:,:,None], 0, 2)
        Y = np.tile(Y, (X.shape[0], 1, 1))
        eq = np.all(np.equal(X, Y), axis = 1)
        eq = np.any(eq, axis = 1)
        return np.nonzero(eq)[0]

To answer the question more specifically, you'd only need to check if the returned array is empty.

Solution 4

This should be much faster it is not O(n^2) like the for-loop solution, but it isn't fully numpythonic. Not sure how better to leverage numpy here

def set_comp(a, b):
   sets_a = set(map(lambda x: frozenset(tuple(x)), a))
   sets_b = set(map(lambda x: frozenset(tuple(x)), b))
   return not sets_a.isdisjoint(sets_b)

View more solutions

Comments

[Admin]

Thanks for your answer but i am looking for a solution based on numpy functions without any loops.

[Mazdak]

oh yes what about this ? : numpy.array_equal(a1, a2)[source]¶ True if two arrays have the same shape and elements, False otherwise.

[Mazdak]

this compare is for arrays and in this case u mast use for subarrays if you dont want to use loop i havnt any idea !

[Admin]

yes this one comes realy close to what i want. but array_equal works with 1d arrays. so i could only use it to check single elements like np.array_equal([1,2,3],[1,2,3]) gives true for example ... but i would have to loop over all the elements of my 2 arrays to check for an intersection which again is not what i want :/

[Mazdak]

yes and i wrote for you you must use this in loop bu if you dont want loop for now i havnt any idea ! but i think for this case we cant dont use loop if there was an function in numpy thay do this it actually use loop because the aim of loop is this !that reduce our calculate !!!

[daniel]

@bluebird7, you'll want to use (i == j).all() in this case. Just FYI

[daniel]

Ah, sorry @CTZhum just saw your response as I was writing mine. It is more succinct, nice :)

[CT Zhu]

Hey @daniel, the numpy solution is even more concise: (A[...,np.newaxis]==B[...,np.newaxis].T).all(1).any(). Cheers.

[Admin]

Hey guys, thanks for your answers! that looks good, i'll try to implement it tomorrow :)

[Admin]

Was able to implement your set method :) seems to work fine for large arrays. thanks a lot!