有效地计算3D numpy阵列沿具有不同面元边缘的轴的直方图 [英] Efficiently calculate histogram of a 3D numpy array along an axis with different bin edges

问题描述

我有一个3D numpy数组，表示为 data ，形状为N x R x C，即N个样本，R行和C列.我想获取样本和行的每种组合的沿列的直方图.但是bin边缘(请参阅bins 长度固定为S的> numpy.histogram )在不同的行上会有所不同，但会在样本之间共享.以这个示例为例，对于第一个样本( data [0] )，其第一行的bin边缘序列与第二行的bin边缘序列不同，但与第一行的bin边缘序列相同来自第二个样本( data [1] ).因此，所有面元边缘序列都存储在形状为R x S的2D numpy数组中，表示为 bin_edges .

I have a 3D numpy array, denoted as data, of shape N x R x C, i.e. N samples, R rows and C columns. I would like to obtain histograms along column for each combination of sample and row. However bin edges (see argument bins in numpy.histogram), of fixed length S, will be different at different rows but are shared across samples. Consider this example for illustration, for the 1st sample (data[0]), bin edge sequence for its 1st row is different from that for its 2nd row, but is the same as that for the 1st row from the 2nd sample (data[1]). Thus all the bin edge sequences are stored in a 2D numpy array of shape R x S, denoted as bin_edges.

我的问题是如何有效地计算直方图?

My question is how to efficiently calculate the histograms?

使用 numpy.histogram ，我能够提出一个可行但相当缓慢的解决方案，如以下代码片段所示

Using numpy.histogram, I was able to come up with a working but fairly slow solution as shown in the below code snippet

```
Get dummy data

    N: number of samples
    R: number of rows (or kernels)
    C: number of columns (or pixels)
    S: number of bins
```
import numpy as np

N, R, C, S = 100, 50, 1000, 10
data = np.random.randn(N, R, C)

# for each row/kernel, pool pixels of all samples
poolsamples = np.swapaxes(data, 0, 1).reshape(R, -1)
# use quantiles as bin edges
percentiles = np.linspace(0, 100, num=(S + 1))
bin_edges = np.transpose(np.percentile(poolsamples, percentiles, axis=1))


```
A working but slow solution of getting histograms along column
```
hist = np.empty((N, R, S))
for idx in np.arange(R):
    bin_edges_i = bin_edges[idx, :]
    counts = np.apply_along_axis(
        lambda a: np.histogram(a, bins=bin_edges_i)[0],
        1, data[:, idx, :])
    hist[:, idx, :] = counts

可能的方向

• 花哨的numpy重塑形式，完全避免使用for循环
• 此问题源于为经过训练的神经网络转发的每个图像提取低端特征而引起的问题.因此，如果直方图提取可以嵌入到TensorFlow图中并最终在GPU上进行，那将是理想的选择！
• 我注意到一个python软件包 fast-histogram 声称比python快7-15倍 numpy.histogram .但是，一维直方图功能只能获取bin的数量，而不是实际bin的位置
• numexpr ?

Possible directions

• Fancy numpy reshape to avoid using for loop at all
• This problem arises from extracting low-end characteristics for each image forwarded through a trained neural network. Therefore, if the histogram extraction can be embedded in TensorFlow graph and ultimately be carried out on GPU, that would be ideal!
• I noticed a python package fast-histogram which claims to be 7-15x faster than numpy.histogram. However 1d histogram function can only takes number of bins instead of actual bin positions
• numexpr?

我很想听听任何输入！预先感谢！

I would love to hear any inputs! Thanks in advance!

推荐答案

使用 np.searchsorted 的 2D 版本:

Making use of 2D version of np.searchsorted : searchsorted2d -

def vectorized_app(data, bin_edges):
    N, R, C = data.shape
    a = np.sort(data.reshape(-1,C),1)
    b = np.repeat(bin_edges[None],N,axis=0).reshape(-1,bin_edges.shape[-1])

    idx = searchsorted2d(a,b)
    idx[:,0] = 0
    idx[:,-1] = a.shape[1]
    out = (idx[:,1:] - idx[:,:-1]).reshape(N,R,-1)
    return out

运行时测试-

In [591]: N, R, C, S = 100, 50, 1000, 10
     ...: data = np.random.randn(N, R, C)
     ...: 
     ...: # for each row/kernel, pool pixels of all samples
     ...: poolsamples = np.swapaxes(data, 0, 1).reshape(R, -1)
     ...: # use quantiles as bin edges
     ...: percentiles = np.linspace(0, 100, num=(S + 1))
     ...: bin_edges = np.transpose(np.percentile(poolsamples, percentiles, axis=1))
     ...: 

In [592]: %timeit org_app(data, bin_edges)
1 loop, best of 3: 481 ms per loop

In [593]: %timeit vectorized_app(data, bin_edges)
1 loop, best of 3: 224 ms per loop

In [595]: np.allclose(org_app(data, bin_edges), vectorized_app(data, bin_edges))
Out[595]: True

在那里的速度超过了 2x .

仔细观察后发现，提出的矢量化方法的瓶颈在于排序本身-

Closer look reveals that the bottleneck with the proposed vectorized method is the sorting itself -

In [594]: %timeit np.sort(data.reshape(-1,C),1)
1 loop, best of 3: 194 ms per loop

我们需要这种排序才能使用 searchsorted .

We need this sorting to use searchsorted.

这篇关于有效地计算3D numpy阵列沿具有不同面元边缘的轴的直方图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！