def greedy(V, d):
cover = set()
uncovered = set(V)
# invariant: all(distance(x,y) > d for x in cover for y in uncovered)
while uncovered:
x = uncovered.pop()
cover.add(x)
uncovered = {y for y in uncovered if distance(x,y) > d}
uncovered.pop() 调用来使这个贪心算法更好，例如选择一个能够“覆盖”最多剩余点的向量。v中的每个点都可以被至少一个距离 <r 的 v1 中的点到达，而且我们还希望v1中的所有点相互之间的距离都大于 r。 - Pierre D首先，我们可以推导出一个无向图，其中顶点是向量v的元素，边缘是不超过r的点对。这不需要距离成为一个适当的度量：它不必满足三角不等式，只需对称(dist(a,b) == dist(b,a))并且为0，如果两个点相等。
在完成这一步之后，我们可以完全忽略距离，并专注于对图进行分区。正如其他人所指出的那样，该分区是一个Vertex Cover问题。但是，它有一个扭曲的地方：我们要求覆盖中的所有顶点都是不相交的（即：v1中的向量不能彼此在距离r以内）。
为了获得图形，我们可以使用高效的KDTree仅计算一次最近邻结构。特别地，我们将使用kd.sparse_distance_matrix(kd, r)来获取所有距离小于r的点对：
from scipy.spatial import KDTree
def get_graph(v, r):
kd = KDTree(v)
sm = kd.sparse_distance_matrix(kd, r)
graph = {i: {i} for i, _ in sm.keys()}
for i, j in sm.keys():
graph[i] |= {j}
return graph
（注意：如果我们想要使用其他距离，.sparse_distance_matrix()有一个参数p）。
图的分区可以通过多种方式完成。@DanR展示了一种贪婪的方法，非常快速，但通常在找到v1的大小时不够优秀。
下面展示了一种暴力方法，如果存在最小解，则保证找到最小解。对于相对较少的向量来说是足够的，并且可以为其他启发式算法的最优解提供基准。
为了加速组合搜索，我们首先观察到通常一些点与任何其他点的距离不超过r（即上面找到的图并未表示v的所有点）。我们可以将这些点放在一边，因为它们必然是v1的一部分，但不会干扰搜索的其余部分（下面称之为“单例”）。其次，如果它们不完全不相交，我们会削减无望的“线索”（组合扩展中的前缀）。这通过切断整个搜索空间的大片区域来显著加速完整搜索。如果没有削减前缀，则完整搜索的时间呈指数增长。
这个的速度是多少？它取决于v向量的分布，以及与距离r（强烈）相关。实际上，它取决于找到的聚类数。为了给出一个想法，在2D中均匀选择v 20个向量，我观察到大约30毫秒。对于40个向量，通常在100毫秒左右，但对于某些r值，它可能会跳到超过2秒。
我们实现了combinations的变体，它具有一个check函数来削减不太有希望的前缀：
from itertools import combinations
def _okall(tup, *args, **kwargs):
return True
def combinations_all(iterable, n0=1, n1=None, check=None, *args, **kwargs):
pool = tuple(iterable)
n = len(pool)
if n0 > n:
return
n1 = n if n1 is None else n1
check = _okall if check is None else check
if n0 < 1:
yield ()
n0 = 1
seed = list(combinations(range(n), n0-1))
for r in range(n0, n1+1):
prev_seed = seed
seed = []
for prefix in prev_seed:
for j in range(max(prefix, default=-1)+1, n):
indices = prefix + (j,)
tup = tuple(pool[i] for i in indices)
if check(tup, *args, **kwargs):
seed.append(indices)
yield tup
例子：
# list all combinations of [0,1,2,3] (of all sizes), excluding those starting
# by (0,1)
>>> list(combinations_all(range(4), check=lambda tup: tup[:2] != (0,1)))
[(0,),
(1,),
(2,),
(3,),
(0, 2),
(0, 3),
(1, 2),
(1, 3),
(2, 3),
(0, 2, 3),
(1, 2, 3)]
v1的索引：def check(tup, graph):
# refuses tup if any one is within reach of another
# optimization: tup[:-1] has already passed this test, so just check the last one
if len(tup) < 2:
return True
reach_of_last = graph[tup[-1]]
prefix = set(tup[:-1])
return prefix.isdisjoint(reach_of_last)
def brute_vq(v, r):
n = v.shape[0]
graph = get_graph(v, r)
to_part = set(graph)
singletons = set(range(n)).difference(to_part)
if not to_part:
return sorted(singletons)
for tup in combinations_all(to_part, check=check, graph=graph):
# here tup are indices that are far apart enough (guaranteed disjoint in reach)
# check if they fully cover to_part
cover = {j for i in tup for j in graph[i]}
if cover == to_part:
v1_idx = sorted(singletons.union(tup))
return v1_idx
示例:
v = np.random.uniform(size=(20, 2))
r_s = [[0.2, 0.3], [0.4, 0.5]]
fig, axes = plt.subplots(nrows=len(r_s), ncols=len(r_s[0]),
figsize=(12,12), sharex=True, sharey=True)
for r, ax in zip(np.ravel(r_s), np.ravel(axes)):
idx = brute_vq(v, r)
ax.set_aspect('equal')
ax.scatter(v[:, 0], v[:, 1])
for i, p in enumerate(v):
ax.text(*p, str(i))
for i in idx:
circ = plt.Circle(v[i], radius=r, color='g', fill=False)
ax.add_patch(circ)
plt.show();
idx 是 v 向量列表的索引，用于生成 v1。在向量量化（VQ）和聚类术语中，这样的解决方案 v1 被称为码本或原型。每张图片显示了给定 r 值的解决方案，它绘制了每个原型的绿色圆圈（以其为中心，半径为 r），以及 v 向量作为蓝点，以便通过视觉验证解决方案。 - Pierre D我的看法：
A B C D | min,max
A 0 3 2 1 | 1,3
B 3 0 5 4 | 3,5
C 2 5 0 6 | 2,6
D 1 4 6 0 | 1,6 # note: D-C violates the triangle inequality
距离矩阵，重新标记：
A D B C | min,max
A 0 1 3 2 | 1,3
D 1 0 4 6 | 1,6
B 3 4 0 5 | 3,5
C 2 6 5 0 | 2,6
现在，对于每一行（或列），取最小值的最大值（或最大值的最小值...） 距离为3时，西北部分完全连接，其余部分仍然部分连接 （对于d=4，B也将进入子集）
A D | B C | min,max
A 0 1 | 3 2 | 1,3
D 1 0 | 4 6 | 1,6
----+ ----------
B 3 4 | 0 5 | 3,5
C 2 6 | 5 0 | 2,6
v1的向量是什么？ - Pierre Dkd = KDTree（v），然后kd.sparse_distance_matrix（kd，r）是一种非常有效的方法，用于查找彼此之间距离更近的所有点对。然后这只是一个带有扭曲的图形覆盖问题：覆盖中相互排斥。 - Pierre D
v1成为一组新向量（不一定是v的成员）。我所知道的一个例外是k-medoids。但这并不完全符合您的问题，只是探索的一个领域。其他感兴趣的领域包括分层聚类。 - Pierre D