稀疏数据差异的有效直方图

Question

我想计算一个数组 A 中的所有元素与另一个数组 B 中的所有元素之间的差异的直方图。

所以我想要以下数据的直方图：

Delta1 = A1-B1
Delta2 = A1-B2
Delta3 = A1-B3
...
DeltaN = A2-B1
DeltaN+1 = A2-B2
DeltaN+2 = A2-B3 
...

这个计算的重点是表明这些数据具有相关性，即使不是每个数据点在另一个数组中都有一个“伙伴”，并且在实践中相关性相当嘈杂。

问题在于这些文件实际上非常大，几 GB，并且向量的所有条目都是 64 位整数，差异非常大。 对我来说，将这些数据转换为二进制数组以便能够使用相关函数和傅立叶变换来计算它似乎是不可行的。

这是一个小示例，可以让您更好地了解我正在查看的内容。 这种在 for 循环中进行 numpy 搜索排序的实现相当慢。

import numpy as np
import matplotlib.pyplot as plt

timetagsA = [668656283,974986989,1294941174,1364697327,\
1478796061,1525549542,1715828978,2080480431,2175456303,2921498771,3671218524,\
4186901001,4444689281,5087334517,5467644990,5836391057,6249837363,6368090967,8344821453,\
8933832044,9731229532]


timetagsB = [13455,1294941188,1715828990,2921498781,5087334530,5087334733,6368090978,9731229545,9731229800,9731249954]

max_delta_t = 500
nbins = 10000 
histo=np.zeros((nbins,2), dtype = float) 
histo[:,0]=np.arange(0,nbins)   

for i in range(0,int(len(timetagsA))):
    delta_t = 0
    j = np.searchsorted(timetagsB,timetagsA[i]) 
    while (np.round(delta_t) < max_delta_t and j<len(timetagsB)):
        delta_t = timetagsB[j] - timetagsA[i] 

        if(delta_t<max_delta_t):
            histo[int(delta_t),1]+=1 

        j = j+1

plt.plot(histo[0:50,1])
plt.show()

如果有人可以帮助我找到一种更快的计算方法，那就太好了。提前致谢！

Answer 1

编辑

下面的解决方案是假设你的数据太大以至于你不能使用np.substract和np.outer然后切片你想要保留的值：

arr_diff = np.subtract.outer(arrB, arrA)
print (arr_diff[(0<arr_diff ) &(arr_diff <max_delta_t)])
# array([ 14,  12,  10,  13, 216,  11,  13, 268], dtype=int64)

使用您的示例数据它可以工作，但不适用于太大的数据集

原始解决方案

我们首先假设您的max_delta_t 小于timetagsB 中两个连续值之间的差值，这是一种简单的方法（然后我们可以尝试推广它）。

#create the array instead of list
arrA = np.array(timetagsA)
arrB = np.array(timetagsB)
max_delta_t = np.diff(arrB).min() - 1 #here it's 202 just for the explanation

您可以以矢量化方式使用np.searchsorted：

# create the array of search 
arr_search = np.searchsorted(arrB, arrA) # the position of each element of arrA in arrB
print (arr_search)
# array([1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 6, 6, 6, 6, 7, 7, 7],dtype=int64)

您可以通过将arrB与arr_search切片来计算arrB的元素与arrA的每个元素对应的差异

# calculate the difference
arr_diff = arrB[arr_search] - arrA
print (arr_diff[arr_diff<max_delta_t]) # finc the one smaller than max_delta_t
# array([14, 12, 10, 13, 11, 13], dtype=int64)

所以你要找的东西然后由np.bincount计算出来

arr_bins = np.bincount(arr_diff[arr_diff<max_delta_t])
#to make it look like histo but not especially necessary
histo = np.array([range(len(arr_bins)),arr_bins]).T

现在的问题是，当max_delta_t 大于arrB 中的两个连续值时，arrA 和arrB 之间存在一些无法使用此方法获得的差异值。这是一种方法，可能不是最有效的，具体取决于数据的值。对于max_delta_t 的任何值

#need an array with the number of elements in arrB for each element of arrA 
# within a max_delta_t range
arr_diff_search = np.searchsorted(arrB, arrA + max_delta_t)- np.searchsorted(arrB, arrA)
#do a loop to calculate all the values you are interested in
list_arr = []
for i in range(arr_diff_search.max()+1):
    arr_diff = arrB[(arr_search+i)%len(arrB)][(arr_diff_search>=i)] - arrA[(arr_diff_search>=i)]
    list_arr.append(arr_diff[(0<arr_diff)&(arr_diff<max_delta_t)])

现在你可以np.concatenate list_arr 并使用np.bincount 如：

arr_bins = np.bincount(np.concatenate(list_arr))
histo = np.array([range(len(arr_bins)),arr_bins]).T