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

Python: scipy/numpy all pairs computation between two 1-D vectors

我有两个列表 l1l2 包含可能具有不同长度的整数,我想在这两个向量之间的每个可能配对之间执行计算。

具体来说,我正在检查每对之间的汉明距离,如果距离足够小,我想 "count" 它。

天真地,这可以实现

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 >= threshold谓词的元素个数,其中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 的实现以使其正确广播。

我查看了 ,但它假定两个列表的长度相同,但此处并非如此。

您可以将 np.bitwise_xor.outernp.binary_reprnp.char.count 一起使用:

import numpy as np

a = np.random.randint(0, 10, size=5)
b = np.random.randint(0, 10, size=5)

binary_repr = np.vectorize(np.binary_repr)
distance = np.char.count(binary_repr(np.bitwise_xor.outer(a, b)), '1') / 32

然后获取匹配项:

matches = np.sum(distance >= threshold)

以下是使用

的三种方法
  1. scipy.spatial.KDTree
  2. scipy.spatial.distance.cdist
  3. 查找table

在一对长度为 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')