In [1]:
import itertools
import math
from scipy.spatial.distance import euclidean
%pylab inline
In [2]:
randst = np.random.RandomState(123456789)
S = randst.rand(20).reshape(-1, 2)
S = np.vstack((S,np.array([0.65, 0.68])))
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')
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)
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
In [7]:
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')
Each core knows all horizontals and it's bounding verticals. Therefore, compute $I_{ij}$ using the available subsets of information
In [8]:
Iij = np.array(list(itertools.product(vert, horiz)))
print Iij
for i in range(len(Iij)):
lookup[tuple(Iij[i])] = -i
Compute $C_{ij}$ as the set of all points in $V_{ij}$ which are closer to a member of $I_{ij}$ than to a member of $V_{j}$ an is not a member of $H_{i}$.
In [9]:
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)
In [12]:
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(r'$V_{j}$', fontsize=18)
ylabel(r'$H_{i}$', fontsize=18)
xticks([])
yticks([])
Cij
Out[12]:
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.
$H_{1} = {2, 0, 1}\\ H_{2} = {3, 4, 5}\\ H_{3} = {10, 9, 7, 6, 8}$
$V_{1} = {2, 3, 10, 9}\\ V_{2} = {4, 6, 7}\\ V_{3} = {0, 1, 5, 8}$
$ C_{3,1} = \{ \} C_{3,2} = \{\} C_{3,3} = \{5\}\\ C_{2,1} = \{2\} C_{2,2} = \{6\} C_{2,3} = \{8\}\\ C_{1,1} = \{3\} C_{1,2} = \{2\} C_{1,3} = \{5,8\}\\$
By way of example, consider the following set computations for the $H_{1}$ row:
$C_{1} = \cup_{j=1}^p C_{1j} = \{3, 2, 5, 8\}\\ H_{1} \cup C{1} = \{2, 0, 1, 3, 5, 8\}$
Nearest Neighbor 2, 0, 1 have already been computer, so a reduced number of computations is required.
How is it that the maximum cardinality of $C_{ij}$ is 8?
In [17]:
In [ ]: