在 numpy 数组中查找相同值序列的长度（运行长度编码）

Question

在 pylab 程序（也可能是 matlab 程序）中，我有一个代表距离的 numpy 数字数组： d[t] 是时间的 距离 t （以及我的时间跨度数据是 len(d) 时间单位）。

我感兴趣的事件是距离低于某个阈值时，我想计算这些事件的持续时间。使用 b = d<threshold 很容易得到一个布尔数组，问题归结为计算 b 中只为真字的长度序列。但我不知道如何有效地做到这一点（即使用 numpy 基元），我求助于遍历数组并进行手动更改检测（即当值从 False 变为 True 时初始化计数器，只要值为 True 就增加计数器，并在值返回到 False 时将计数器输出到序列）。但这非常慢。

如何有效地检测 numpy 数组中的那种序列？

下面是一些 python 代码来说明我的问题：第四个点需要很长时间才能出现（如果没有，请增加数组的大小）

 from pylab import *

threshold = 7

print '.'
d = 10*rand(10000000)

print '.'

b = d<threshold

print '.'

durations=[]
for i in xrange(len(b)):
    if b[i] and (i==0 or not b[i-1]):
        counter=1
    if  i>0 and b[i-1] and b[i]:
        counter+=1
    if (b[i-1] and not b[i]) or i==len(b)-1:
        durations.append(counter)

print '.'

原文由 Gyom 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

虽然不是 numpy 原语， itertools 函数通常非常快，所以一定要试试这个（当然还要测量包括这个在内的各种解决方案的时间）：

 def runs_of_ones(bits):
  for bit, group in itertools.groupby(bits):
    if bit: yield sum(group)

如果您确实需要列表中的值，当然可以使用 list(runs_of_ones(bits))；但也许列表理解可能会稍微快一些：

 def runs_of_ones_list(bits):
  return [sum(g) for b, g in itertools.groupby(bits) if b]

转向“numpy-native”的可能性，怎么样：

 def runs_of_ones_array(bits):
  # make sure all runs of ones are well-bounded
  bounded = numpy.hstack(([0], bits, [0]))
  # get 1 at run starts and -1 at run ends
  difs = numpy.diff(bounded)
  run_starts, = numpy.where(difs > 0)
  run_ends, = numpy.where(difs < 0)
  return run_ends - run_starts

再次强调：一定要在适合您的实际示例中相互比较解决方案！

原文由 Alex Martelli 发布，翻译遵循 CC BY-SA 2.5 许可协议

Answer 2

适用于任何数组的完全 numpy 向量化和通用 RLE（也适用于字符串、布尔值等）。

输出运行长度、起始位置和值的元组。

 import numpy as np

def rle(inarray):
        """ run length encoding. Partial credit to R rle function.
            Multi datatype arrays catered for including non Numpy
            returns: tuple (runlengths, startpositions, values) """
        ia = np.asarray(inarray)                # force numpy
        n = len(ia)
        if n == 0:
            return (None, None, None)
        else:
            y = ia[1:] != ia[:-1]               # pairwise unequal (string safe)
            i = np.append(np.where(y), n - 1)   # must include last element posi
            z = np.diff(np.append(-1, i))       # run lengths
            p = np.cumsum(np.append(0, z))[:-1] # positions
            return(z, p, ia[i])

相当快（i7）：

 xx = np.random.randint(0, 5, 1000000)
%timeit yy = rle(xx)
100 loops, best of 3: 18.6 ms per loop

多种数据类型：

 rle([True, True, True, False, True, False, False])
Out[8]:
(array([3, 1, 1, 2]),
 array([0, 3, 4, 5]),
 array([ True, False,  True, False], dtype=bool))

rle(np.array([5, 4, 4, 4, 4, 0, 0]))
Out[9]: (array([1, 4, 2]), array([0, 1, 5]), array([5, 4, 0]))

rle(["hello", "hello", "my", "friend", "okay", "okay", "bye"])
Out[10]:
(array([2, 1, 1, 2, 1]),
 array([0, 2, 3, 4, 6]),
 array(['hello', 'my', 'friend', 'okay', 'bye'],
       dtype='|S6'))

与上面的 Alex Martelli 相同的结果：

 xx = np.random.randint(0, 2, 20)

xx
Out[60]: array([1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1])

am = runs_of_ones_array(xx)

tb = rle(xx)

am
Out[63]: array([4, 5, 2, 5])

tb[0][tb[2] == 1]
Out[64]: array([4, 5, 2, 5])

%timeit runs_of_ones_array(xx)
10000 loops, best of 3: 28.5 µs per loop

%timeit rle(xx)
10000 loops, best of 3: 38.2 µs per loop

比 Alex 稍慢（但仍然非常快），而且更灵活。

原文由 Thomas Browne 发布，翻译遵循 CC BY-SA 4.0 许可协议