

In [1]:

    
import itertools
import math
from scipy.spatial.distance import euclidean
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
randst = np.random.RandomState(123456789)
S = randst.rand(20).reshape(-1, 2)
S = np.vstack((S,np.array([0.65, 0.68])))

1.



In [3]:

    
#Globally Sort by y-coordinate
S = S[S[:, 1].argsort()]



In [4]:

    
#Create horizontal (and vertical bounding lines - used later)
horiz = vert = np.linspace(0,1.0,4)
print horiz
print vert

x, y = np.meshgrid(horiz, vert)
plot(x, y)
plot(y, x)
plot(S[:,0],S[:,1], 'o')

labels = np.arange(len(S))
for label, x, y in zip(labels, S[:, 0], S[:, 1]):
    annotate(label, 
             xy = (x, y), xytext = (6, 3),
             textcoords = 'offset points', ha = 'right', va = 'bottom')









    



[ 0.          0.33333333  0.66666667  1.        ]
[ 0.          0.33333333  0.66666667  1.        ]



In [5]:

    
#Using the global sort, processor p stores the bounding horizontals and the pts contained within
lookup = {}
for i in range(len(S)):
    lookup[tuple(S[i])] = i

Hij = {}
for i in range(len(S)):
    Hij[i] = (np.inf, -1)
    
for i in range(1,4):
    p1 = S[np.where((S[:,1] > horiz[i-1]) & (S[:,1] < horiz[i]))]

    #This is replicated for each core in the set P
    for p in itertools.combinations(p1, 2):  
        a, b = p
        ka = lookup[tuple(a)]
        kb = lookup[tuple(b)]
        dist = euclidean(a,b)
        if dist < Hij[ka][0]:
            Hij[ka] = (dist, kb)
        if dist < Hij[kb][0]:
            Hij[kb] = (dist, ka)

2



In [6]:

    
#Process all vertical columns as above
Vij = {}
for i in range(len(S)):
    Vij[i] = (np.inf, -1)

for i in range(1,4):
    p1 = S[np.where((S[:,0] > vert[i-1]) & (S[:,0] < vert[i]))]

    #This is replicated for each core in the set P
    for p in itertools.combinations(p1, 2):  
        a, b = p
        ka = lookup[tuple(a)]
        kb = lookup[tuple(b)]
        dist = euclidean(a,b)
        if dist < Vij[ka][0]:
            Vij[ka] = (dist, kb)
        if dist < Vij[kb][0]:
            Vij[kb] = (dist, ka)

print Vij









    



{0: (0.08669582752535794, 1), 1: (0.08669582752535794, 0), 2: (0.12304015044166737, 3), 3: (0.12304015044166737, 2), 4: (0.18093137062065548, 7), 5: (0.11790441517871283, 8), 6: (0.14440184923573685, 7), 7: (0.14440184923573685, 6), 8: (0.11790441517871283, 5), 9: (0.08215464548813665, 10), 10: (0.08215464548813665, 9)}



In [10]:

    
for k, v in Hij.iteritems():
    print k, "Horizontal:",v[1], "Vertical:",Vij[k][1]

x, y = np.meshgrid(horiz, vert)
plot(x, y)
plot(y, x)
plot(S[:,0],S[:,1], 'o')

labels = np.arange(len(S))
for label, x, y in zip(labels, S[:, 0], S[:, 1]):
    annotate(label, 
             xy = (x, y), xytext = (6, 3),
             textcoords = 'offset points', ha = 'right', va = 'bottom')









    



0 Horizontal: 1 Vertical: 1
1 Horizontal: 0 Vertical: 0
2 Horizontal: 0 Vertical: 3
3 Horizontal: 4 Vertical: 2
4 Horizontal: 5 Vertical: 7
5 Horizontal: 4 Vertical: 8
6 Horizontal: 8 Vertical: 7
7 Horizontal: 6 Vertical: 6
8 Horizontal: 6 Vertical: 5
9 Horizontal: 10 Vertical: 10
10 Horizontal: 9 Vertical: 9

Each core knows all horizontals and it's bounding verticals. Therefore, compute $I_{ij}$ using the available subsets of information



In [11]:

    
Iij = np.array(list(itertools.product(vert, horiz)))
print Iij
for i in range(len(Iij)):
    lookup[tuple(Iij[i])] = -i









    



[[ 0.          0.        ]
 [ 0.          0.33333333]
 [ 0.          0.66666667]
 [ 0.          1.        ]
 [ 0.33333333  0.        ]
 [ 0.33333333  0.33333333]
 [ 0.33333333  0.66666667]
 [ 0.33333333  1.        ]
 [ 0.66666667  0.        ]
 [ 0.66666667  0.33333333]
 [ 0.66666667  0.66666667]
 [ 0.66666667  1.        ]
 [ 1.          0.        ]
 [ 1.          0.33333333]
 [ 1.          0.66666667]
 [ 1.          1.        ]]

Compute $C_{ij}$ as the set of all points which are closer to a member of $I_{ij}$ than to a member of $V_{j}$



In [12]:

    
Cij = {}
for i in range(len(S)):
    Cij[i] = (np.inf, 0)

for i in range(1,4):
    p1 = S[np.where((S[:,0] > vert[i-1]) & (S[:,0] < vert[i]))]
    bounds = Iij[np.where((Iij[:,0] >= vert[i-1]) & (Iij[:,0] <= vert[i]))]

    #This is replicated for each core in the set P
    for p in itertools.product(p1, bounds):
        a, b = p
        ka = lookup[tuple(a)]
        kb = lookup[tuple(b)]
        dist = euclidean(a,b)
        if ka in Vij.keys():
            if dist < Vij[ka][0]:
                Cij[ka] = (dist, kb)
        if kb in Vij.keys():
            if dist < Vij[kb][0]:
                Cij[kb] = (dist, ka)
remove = []
for k, v in Cij.iteritems():
    if v[0] == np.inf:
        remove.append(k)
for k in remove:
    Cij.pop(k, None)

$C_{ij}$ is now a dict of those points which are closer to an



In [18]:

    
x, y = np.meshgrid(horiz, vert)
plot(x, y)
plot(y, x)
plot(x,y, 'ro')
plot(S[:,0],S[:,1], 'o')

labels = np.arange(len(S))
Iijlabel = np.arange(len(Iij))
Iijlabel *= -1

#Label the pts
for label, x, y in zip(labels, S[:, 0], S[:, 1]):
    annotate(label, 
             xy = (x, y), xytext = (6, 3),
             textcoords = 'offset points', ha = 'right', va = 'bottom')
    
for label, x, y in zip(Iijlabel, Iij[:, 0], Iij[:, 1]):
    annotate(label, 
             xy = (x, y), xytext = (-2, 2),
             textcoords = 'offset points', ha = 'right', va = 'bottom')  
xlabel('J')
ylabel('I')
xticks([])
yticks([])

#2, 3 are in the set C_{1,1}
#6 is in the set C_{3,2}
#5, 8 are in the set C_
Cij









    Out[18]:





{2: (0.10258529611325132, -1),
 3: (0.023634018340308822, -1),
 5: (0.023543386877573954, -10),
 6: (0.021343747458109505, -10),
 8: (0.09989044874270481, -10)}

Checks for Iij at the boundry corners are redundant, in that, we know that the nearest neighbor is not off the edge of the ROI.

3.



In [ ]:

    
#Compute the nearest for the union of Hi with Ci, e.g. all Hi in this row and all Ci in this row