Python:scipy/numpy 两个一维向量之间的所有对计算

Question

我有两个列表 l1 和 l2 ，其中包含可能具有不同长度的整数，并且我想在这两个向量之间的每个可能配对之间执行计算。

具体来说，我正在检查每对之间的汉明距离，如果距离足够小，我想“计算”它。

天真的，这可以实现

def hamming_distance(n1: int, n2: int) -> float:
    return bin(n1 ^ n2).count('1')/32.0

matches = 0

for n1 in l1:
    for n2 in l2:
        sim = 1 - hamming_distance(n1, n2)

        if sim >= threshold:
            matches += 1

但这并不是很快。

我没有成功地尝试利用scipy.spatial.distance.cdist，我认为我首先要计算所有对之间的汉明距离，如 scipy.spatial.cdist documentation states它将

Compute distance between each pair of the two collections of inputs.

然后计算满足谓词1 - d >=阈值的元素数量，其中d是汉明距离，即

from scipy.spatial.distance import cdist

l1 = l1.reshape(-1, 2)  # After np.array
l2 = l2.reshape(-1, 2)
r = cdist(l1, l2, 'hamming')
matches = np.count_nonzero(1 - r >= threshold)

但是各个解找到的匹配数是不同的。我注意到可以使用函数 cdist(XA, XB, f) 调用 cdist 但我还没有成功编写 的实现hamming_distance 以便其正确广播。

我看过this question/answer但它假设两个列表的长度相同，但这里的情况并非如此。

Answer 1

以下是使用的三种方法

1. scipy.spatial.KDTree
2. scipy.spatial.distance.cdist
3. 查找表

在一对长度为 100 和 200 的 32 位 int 向量上，它们都给出相同的结果；他们的速度比较如下:

count_sim_kd 16.408800622448325 ms
count_sim_cd 12.41896384395659 ms
count_sim_lu 0.8755046688020229 ms

因此，在这个问题规模下，查找以巨大优势获胜。

代码:

import numpy as np
from scipy.spatial import cKDTree as KDTree
from scipy.spatial.distance import cdist

l1 = np.random.randint(0,2**32,100)
l2 = np.random.randint(0,2**32,200)
threshold = 10/32

def hamming_distance(n1: int, n2: int) -> float:
    return bin(n1 ^ n2).count('1')/32.0

matches = 0

for n1 in l1:
    for n2 in l2:
        sim = 1 - hamming_distance(n1, n2)

        if sim >= threshold:
            matches += 1

def count_sim_kd(a,b,th):
    A,B = (KDTree(np.unpackbits(x[:,None].view(np.uint8),axis=1))
           for x in (a,b))
    return A.sparse_distance_matrix(B,max_distance=32-int(32*th),p=1).nnz

def count_sim_cd(a,b,th):
    A,B = (np.unpackbits(x[:,None].view(np.uint8),axis=1) for x in (a,b))
    return np.count_nonzero(cdist(A,B,"minkowski",p=1)<=32-int(32*th))


lu = sum(np.unravel_index(np.arange(256),8*(2,)))
def count_sim_lu(a,b,th):
    return np.count_nonzero(lu[(a[:,None,None]^b[None,:,None])
                               .view(np.uint8)].sum(2)<=32-int(32*th))

from timeit import timeit
for f in (count_sim_kd,count_sim_cd,count_sim_lu):
    assert f(l1,l2,threshold)==matches
    print(f.__name__,timeit(lambda:f(l1,l2,threshold),number=100)*10,'ms')