Numpy 分块缩减运算

9

我是一名经验丰富的Numpy用户,但我无法找到以下问题的解决方案。假设有以下数组:

# sorted array of times
t = numpy.cumsum(numpy.random.random(size = 100))
#  some values associated with the times
x = numpy.random.random(size=100)
# some indices into the time/data array
indices = numpy.cumsum(numpy.random.randint(low = 1, high=10,size = 20)) 
indices = indices[indices <90] # respect size of 100
if len(indices) % 2: # make number of indices even
    indices = indices[:-1]

# select some starting and end indices
istart = indices[0::2]
iend   = indices[1::2]

我现在想要的是根据istartiend所表示的间隔来缩小值数组x。即:
# e.g. use max reduce, I'll probably also need mean and stdv
what_i_want = numpy.array([numpy.max(x[is:ie]) for is,ie in zip(istart,iend)])

我已经谷歌了很多,但我只能找到通过stride_tricks进行分块操作的方法,这仅允许使用常规块。我无法找到一种不需要执行Python循环的解决方案:-(

在我的实际应用中,数组要大得多,并且性能确实很重要,所以我目前使用numba.jit

是否有任何我遗漏的numpy函数可以做到这一点?


x 的值是否总是在 [0,1) 范围内的浮点数? - Divakar
不管怎样,x通常都是一个更复杂的数组结构。 - Marti Nito
3个回答

7

您是否看过ufunc.reduceat?使用np.maximum,您可以进行如下操作:

>>> np.maximum.reduceat(x, indices)

该函数返回沿着切片 x[indices[i]:indices[i+1]] 的最大值。如果你想获取 (x[indices[2i]:indices[2i+1]),可以执行以下操作:

>>> np.maximum.reduceat(x, indices)[::2]

如果您不介意进行额外的计算,可以使用x[inidices[2i-1]:indices[2i]]。这将产生以下结果:

>>> numpy.array([numpy.max(x[ib:ie]) for ib,ie in zip(istart,iend)])
array([ 0.60265618,  0.97866485,  0.78869449,  0.79371198,  0.15463711,
        0.72413702,  0.97669218,  0.86605981])

>>> np.maximum.reduceat(x, indices)[::2]
array([ 0.60265618,  0.97866485,  0.78869449,  0.79371198,  0.15463711,
        0.72413702,  0.97669218,  0.86605981])

太棒了,这正是我在寻找的。我可以将所有索引保存在一个数组中,这样就不需要进行任何额外的计算了 :) 也许我应该提高我的搜索技巧... - Marti Nito

0
你可以这样使用 numpy.r_
what_i_want = np.array([np.max(x[np.r_[ib:ie]]) for ib,ie in zip(istart,iend)])

0

(一种非numpy解决方案,使用astropy

import numpy as np
from astropy.nddata.utils import block_reduce
data = np.arange(16).reshape(4, 4)
block_reduce(data, 2)    

将会转换:

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

至:

array([[10, 18],
       [42, 50]])

请查看this以获取更多示例。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接