Есть ли лучший способ реализации гистограммы?

Question

У меня есть 2D массив uint16 numpy, я хочу вычислить гистограмму для этого массива. Я использую функцию:

def calc_hist(source):
    hist = np.zeros(2**16, dtype='uint16')
    for i in range(source.shape[0]):
       for j in range(source.shape[1]):
           hist[source[i, j]] = hist[source[i, j]] + 1
    return hist 

Эта функция выполняется слишком долго. Как я понимаю, в модуле numpy есть функция гистограммы, но я не могу понять, как ее использовать. Я пробовал:

hist,_ = np.histogram(source.flatten(), bins=range(2**16))

Но я получаю другие результаты, чем моя собственная функция. Как я могу вызвать numpy .histogram для достижения того же результата? или есть еще варианты?

Answer 1

Для ввода с типом данных uint16, numpy.bincount должно работать хорошо:

hist = np.bincount(source.ravel(), minlength=2**16)

Ваша функция делает почти то же, что bincount, но bincount реализовано в C.

Например, следующие проверки проверяют, что использование bincount дает тот же результат, что и ваша calc_hist функция:

In [159]: rng = np.random.default_rng()

In [160]: x = rng.integers(0, 2**16, size=(1000, 1000))                                      

In [161]: h1 = calc_hist(x)

In [162]: h2 = np.bincount(x.ravel(), minlength=2**16)

In [163]: (h1 == h2).all()  # Verify that they are the same.
Out[163]: True

Проверьте производительность с помощью команды i python %timeit. Вы можете видеть, что использование bincount намного быстрее .

In [164]: %timeit calc_hist(x)
2.66 s ± 21.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [165]: %timeit np.bincount(x.ravel(), minlength=2**16)
3.13 ms ± 100 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Answer 2

Как указал Корли Бригман, пропуск bins=range(x) определяет края корзины [1]. Таким образом, вы получите x-1 с соответствующими ребрами [0, 1), [1, 2), ..., [x-1, x].

В вашем случае вы будете есть 2 ^ 16-1 бункеров. Чтобы исправить это, просто используйте range(2**16+1).

[1] https://numpy.org/doc/stable/reference/generated/numpy.histogram.html?highlight=histogram#numpy .histogram