在一个二维的Numpy数组中计算数字对的频率的最有效方法

Question

假设我有以下2D数组：

import numpy as np

np.random.seed(123)
a = np.random.randint(1, 6, size=(5, 3))

生成：

In [371]: a
Out[371]:
array([[3, 5, 3],
       [2, 4, 3],
       [4, 2, 2],
       [1, 2, 2],
       [1, 1, 2]])

有没有更高效的方法（如Numpy，Pandas等）来计算所有数字对的频率比以下解决方案？

from collections import Counter
from itertools import combinations

def pair_freq(a, sort=False, sort_axis=-1):
    a = np.asarray(a)
    if sort:
        a = np.sort(a, axis=sort_axis)
    res = Counter()
    for row in a:
        res.update(combinations(row, 2))
    return res

res = pair_freq(a)

生产类似于那样的东西：

In [38]: res
Out[38]:
Counter({(3, 5): 1,
         (3, 3): 1,
         (5, 3): 1,
         (2, 4): 1,
         (2, 3): 1,
         (4, 3): 1,
         (4, 2): 2,
         (2, 2): 2,
         (1, 2): 4,
         (1, 1): 1})

或者：

In [39]: res.most_common()
Out[39]:
[((1, 2), 4),
 ((4, 2), 2),
 ((2, 2), 2),
 ((3, 5), 1),
 ((3, 3), 1),
 ((5, 3), 1),
 ((2, 4), 1),
 ((2, 3), 1),
 ((4, 3), 1),
 ((1, 1), 1)]

PS：结果数据集可能会看起来不同 - 例如像一个多级Pandas DataFrame或其他内容。

我试图增加a数组的维度，并使用np.isin()和一个由所有成对组合的列表一起，但我仍然无法摆脱循环。

更新：

(a) 你只关心2个数字的组合的频率（而不关心3个数字的组合频率）吗？

是的，我只关心成对（2个数字）的组合。

(b) 你是否认为(3,5)和(5,3)是不同的，还是把它们看作是相同的两次出现？

实际上这两种方法都可以 - 如果需要，我可以事先对数组进行排序：

a = np.sort(a, axis=1)

更新2：

您是否希望仅由a和b的源列引起(a,b)和(b,a)之间的区别，还是其他方式也可以？ 请考虑三行[[1,2,1], [3,1,2], [1,2,5]]以理解此问题。在这里，您认为输出应该是什么？哪些应该是不同的2元组以及它们的频率应该是什么？

In [40]: a = np.array([[1,2,1],[3,1,2],[1,2,5]])

In [41]: a
Out[41]:
array([[1, 2, 1],
       [3, 1, 2],
       [1, 2, 5]])

我期望得到以下结果：

In [42]: pair_freq(a).most_common()
Out[42]:
[((1, 2), 3),
 ((1, 1), 1),
 ((2, 1), 1),
 ((3, 1), 1),
 ((3, 2), 1),
 ((1, 5), 1),
 ((2, 5), 1)]

因为更加灵活，所以我想将 (a, b) 和 (b,a) 视为同一对元素。 我可以这样做：

In [43]: pair_freq(a, sort=True).most_common()
Out[43]: [((1, 2), 4), ((1, 1), 1), ((1, 3), 1), ((2, 3), 1), ((1, 5), 1), ((2, 5), 1)]

- MaxU - stand with Ukraine

---

(a) 你只对两个数字的组合频率感兴趣（不关心三个数字的组合频率）吗？ (b) 你想把(3,5)和(5,3)视为不同的事物，还是将它们视为相同的两个出现次数？ - fountainhead

相关问题。行[1,2,2]生成了两个2元组(1,2)（因为数字2出现在不同的两列中），而你目前的解决方案将这两个2元组视为不同。这种区别重要吗？ - fountainhead

@fountainhead，是的，我认为将(a, b)和(b, a)视为不同的对会给我们更多的灵活性。 - MaxU - stand with Ukraine

@fountainhead。行[1,2,2]生成了2元组(1,2)两次(因为2出现在两列中)，而您当前的解决方案将这两个2元组视为不同。这种区别很重要吗？- 是的，在这种情况下，我想把(1,2)计算两次。 - MaxU - stand with Ukraine

你想让(a,b)和(b,a)之间的区别仅由a和b的源列引起，还是其他情况也可以？为了理解这个问题，请考虑三行[1,2,1]、[3,1,2]和[1,2,5]。你认为这里的输出应该是什么？应该是哪些不同的二元组以及它们的频率？ - fountainhead

Answer 1

如果您的元素不是太大的非负整数，则bincount速度快：最初的回答。

from collections import Counter
from itertools import combinations
import numpy as np

def pairs(a):
    M = a.max() + 1
    a = a.T
    return sum(np.bincount((M * a[j] + a[j+1:]).ravel(), None, M*M)
               for j in range(len(a) - 1)).reshape(M, M)

def pairs_F_3(a):
    M = a.max() + 1
    return (np.bincount(a[1:].ravel() + M*a[:2].ravel(), None, M*M) +
            np.bincount(a[2].ravel() + M*a[0].ravel(), None, M*M))

def pairs_F(a):
    M = a.max() + 1
    a = np.ascontiguousarray(a.T) # contiguous columns (rows after .T)
                                  # appear to be typically perform better
                                  # thanks @ning chen
    return sum(np.bincount((M * a[j] + a[j+1:]).ravel(), None, M*M)
               for j in range(len(a) - 1)).reshape(M, M)

def pairs_dict(a):
    p = pairs_F(a)
    # p is a 2D table with the frequency of (y, x) at position y, x
    y, x = np.where(p)
    c = p[y, x]
    return {(yi, xi): ci for yi, xi, ci in zip(y, x, c)}

def pair_freq(a, sort=False, sort_axis=-1):
    a = np.asarray(a)
    if sort:
        a = np.sort(a, axis=sort_axis)
    res = Counter()
    for row in a:
        res.update(combinations(row, 2))
    return res


from timeit import timeit
A = [np.random.randint(0, 1000, (1000, 120)),
     np.random.randint(0, 100, (100000, 12))]
for a in A:
    print('shape:', a.shape, 'range:', a.max() + 1)
    res2 = pairs_dict(a)
    res = pair_freq(a)
    print(f'results equal: {res==res2}')
    print('bincount', timeit(lambda:pairs(a), number=10)*100, 'ms')
    print('bc(F)   ', timeit(lambda:pairs_F(a), number=10)*100, 'ms')
    print('bc->dict', timeit(lambda:pairs_dict(a), number=10)*100, 'ms')
    print('Counter ', timeit(lambda:pair_freq(a), number=4)*250,'ms')

示例运行：

shape: (1000, 120) range: 1000
results equal: True
bincount 461.14772390574217 ms
bc(F)    435.3669326752424 ms
bc->dict 932.1215840056539 ms
Counter  3473.3258984051645 ms
shape: (100000, 12) range: 100
results equal: True
bincount 89.80463854968548 ms
bc(F)    43.449611216783524 ms
bc->dict 46.470773220062256 ms
Counter  1987.6734036952257 ms

- Paul Panzer

---

这个解决方案可能很快，但解决方案提供者更快！ - fountainhead

假设 a 的形状为 m,n。当 m >> n 时，使用 'F' 排序的 a 代码运行更快。 - Happy Boy

@ningchen 实际上，改为列主序似乎大多数情况下都是值得的；对我尝试的所有参数都有效！ - Paul Panzer

Answer 2

我有一个想法，代码如下。 我的代码最大的缺点是随着列数的增加而运行非常缓慢，并且比@Paul Panzer的代码慢。我向Paul Panzer道歉。如果你想更快地运行，只需忽略num_to_items函数。因为(1, 1)等于1*2 ** 20 + 1。

import numpy as np
from random import choice
from itertools import izip
from scipy.sparse import csr_matrix, csc_matrix
from scipy import sparse as sp


c_10 = np.array([[choice(range(1, 10)) for _ in range(3)] for _ in range(1000)])
c_1000 = np.array([[choice(range(1, 1000)) for _ in range(3)] for _ in range(1000)])

def _bit_to_items(num):
    return (num >> 20, num & 0b1111111111111111111)


def unique_bit_shit(c):
    cc = c << 20 # suppose that: 2**20 > max(c)

    dialog_mtx_1 = np.array([[1, 0, 0],
                             [1, 0, 0],
                             [0, 1, 0]])

    dialog_mtx_2 = np.array([[0, 1, 0],
                             [0, 0, 1],
                             [0, 0, 1]])

    dialog_mtx_1 = dialog_mtx_1.T
    dialog_mtx_2 = dialog_mtx_2.T

    pairs = cc.dot(dialog_mtx_1) + c.dot(dialog_mtx_2)
    pairs_num, count = np.unique(pairs, return_counts=True)
    return [(_bit_to_items(num), v) for num, v in izip(pairs_num, count)]



def _dot_to_items(num):
    # 2**20 is 1048576
    return (num / 1048576, num % 1048576) 


def unique_dot(c):
    dialog_mtx_3 = np.array([[2**20, 2**20, 0],
                             [1,     0,     2**20],
                             [0,     1,     1]])

    pairs = c.dot(dialog_mtx_3)

    pairs_num, count = np.unique(pairs, return_counts=True)
    return [(_dot_to_items(num), v) for num, v in izip(pairs_num, count)]

- Happy Boy

---

当a中元素的范围非常大时，实际上正确的条件应该是“如果范围与数组大小相比较大”。例如，如果您将数组的行数增加到10000，则bincount会更快。此外，您的代码是否适用于除3列之外的其他内容？ - Paul Panzer

正如你所说，我的测试并没有覆盖所有的情况。“我的代码更快”并不准确。 - Happy Boy

@Paul Panzer。随着列数的增加，np.dot和np.unique的计算时间也会增加。因此，当列数超过100时，unique_dot将变得非常缓慢。使用csc_matrix.dot可以提高速度，但仍然比您的代码慢。 - Happy Boy