高效地在Python中检测符号变化

Question

我希望你能像这位用户一样做到以下内容：Python - count sign changes。然而，我需要对其进行优化以实现超快速运行。简单来说，我想要获取时间序列并告诉每次穿过零点（改变符号）的时间。我想要记录两个零点之间的时间。由于这是实际数据 (32 位浮点数)，我怀疑我将不会得到一个完全为零的数字，所以这并不重要。我目前已经有了定时程序，所以我会测试您的结果，看看谁赢了。我的解决方案提供的时间单位是微秒：

open data       8384
sign data       8123
zcd data        415466

正如您所看到的，零点检测器是较慢的部分。以下是我的代码。

import numpy, datetime

class timer():
    def __init__(self):
        self.t0 = datetime.datetime.now()
        self.t = datetime.datetime.now()
    def __call__(self,text='unknown'):
        print text,'\t',(datetime.datetime.now()-self.t).microseconds
        self.t=datetime.datetime.now()

def zcd(data,t):
    sign_array=numpy.sign(data)
    t('sign data')
    out=[]
    current = sign_array[0]
    count=0
    for i in sign_array[1:]:
        if i!=current:
            out.append(count)
            current=i
            count=0
        else: count+=1
    t('zcd data')
    return out

def main():
    t = timer()
    data = numpy.fromfile('deci.dat',dtype=numpy.float32)
    t('open data')
    zcd(data,t)

if __name__=='__main__':
    main()

- chriscauley

---

你知道有一个叫做"timeit"的模块吗？ :) - Radomir Dopieralski

有趣...我更喜欢我的方法，因为它可以在整个函数中使用。你可以每隔几行放一个t()，快速找到瓶颈。如果我只想计时函数，我会使用Linux的 $ time python zcd.py。 - chriscauley

我猜测 time('sign data') 这一行应该是想写成 t('sign data')，是吗？ - Muhammad Alkarouri

@Muhammad Alkarouri - 好的，谢谢。我会修复它。 - chriscauley

可能是Python-计算符号变化的重复问题。 - Serge Stroobandt

Answer 1

怎么样：

import numpy
a = [1, 2, 1, 1, -3, -4, 7, 8, 9, 10, -2, 1, -3, 5, 6, 7, -10]
zero_crossings = numpy.where(numpy.diff(numpy.sign(a)))[0]

输出：

> zero_crossings
array([ 3,  5,  9, 10, 11, 12, 15])

即，zero_crossings将包含零穿越发生之前元素的索引。如果想要零穿越发生后的元素，只需在该数组中加1。

- Jim Brissom

---

我认为你把它搞反了；zero_crossings包含零交叉发生的元素之前的索引，如果你想要之后的元素，需要将数组加1。否则，回答非常好，简洁明了！ - staticfloat

当数组中存在零时，此方法无法正常工作。它会将零检测两次！例如：a = [2,1,0,-1,2] 将得到 array([1, 2, 3]) 的结果。 - YuppieNetworking

如果您只对计数（而不是索引）感兴趣，那么删除所有的0就可以了。 np.where(np.diff(np.sign([i for i in a if i])))[0].shape[0] - Vishal Gupta

这个函数几乎可以实现 numpy.sign(0) = 0 和 numpy.sign(2) = 1 以及 numpy.sign(-2) = -1。因此，您可能需要使用 numpy.where(numpy.diff(numpy.sign(a) >= 0))[0]。 - HaskellElephant

这种逻辑如何应用以仅查找负到正的符号变化？例如，从序列[-2，-4，-2，1，2，8，-1，-1，0] 中，我需要输出[0，0，0，1，0，0，0，0，1]。 - user6400946

检查 np.diff 是否返回 +2 或 -2 来处理 0 怎么样？像这样 diff = np.diff(np.sign(array)), np.nonzero((diff == -2) | (diff == 2)) - dc_Bita98

Answer 2

正如Jay Borseth所指出的那样，被接受的答案无法正确处理包含0的数组。

我建议使用以下方法：

import numpy as np
a = np.array([-2, -1, 0, 1, 2])
zero_crossings = np.where(np.diff(np.signbit(a)))[0]
print(zero_crossings)
# output: [1]

因为a)使用numpy.signbit()比numpy.sign()稍微快一些，因为它的实现更简单，所以我猜想；并且b) 它正确地处理了输入数组中的零。

然而，可能存在一个缺点：如果你的输入数组以零开始和结束，它将在开头找到一个零交叉点，但不会在结尾找到。

import numpy as np
a = np.array([0, -2, -1, 0, 1, 2, 0])
zero_crossings = np.where(np.diff(np.signbit(a)))[0]
print(zero_crossings)
# output: [0 2]

- Dominik Neise

---

嗯，那么$[-2,-1,0,-1,-2,0]$怎么办……没有交叉只是接触，但确实有一个答案。把零算作正数也不是最终的解决方案，我想。 - mikuszefski

@mikuszefski 你是对的！[1, 2, 0, -1, 0, 0, -1, 2] 应该产生 2 个零点交叉，但实际上没有。 - Serge Stroobandt

Answer 3

另一种计算零交叉并从代码中挤出更多毫秒的方法是使用 `nonzero` 并直接计算符号。假设您有一个一维数组 `data`：

def crossings_nonzero_all(data):
    pos = data > 0
    npos = ~pos
    return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]

或者，如果您只想计算通过零点的特定方向（例如从正到负）的零交叉点数量，则这种方法速度更快：

def crossings_nonzero_pos2neg(data):
    pos = data > 0
    return (pos[:-1] & ~pos[1:]).nonzero()[0]

在我的计算机上，这种方法比where(diff(sign))方法快一些(对于包含20个周期、总共40个交叉点的正弦样本数组进行计时，该数组包含10000个样本)：

$ python -mtimeit 'crossings_where(data)'
10000 loops, best of 3: 119 usec per loop

$ python -mtimeit 'crossings_nonzero_all(data)'
10000 loops, best of 3: 61.7 usec per loop

$ python -mtimeit 'crossings_nonzero_pos2neg(data)'
10000 loops, best of 3: 55.5 usec per loop

- lmjohns3

---

你可以将 (pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:]) 缩短为 pos[:-1] ^ npos[1:]，其中 ^ 是异或运算符。 - Bas Swinckels

crossings_nonzero_pos2neg([1,2,-1,1,2]) 回溯（最近的调用最先）： File "<ipython-input-3-21f24f68064f>", line 1, in <module> crossings_nonzero_pos2neg([1,2,-1,1,2]) File "<ipython-input-2-80149113a324>", line 2, in crossings_nonzero_pos2neg pos = data > 0TypeError: '>' not supported between instances of 'list' and 'int' - Mainland

请使用 "numpy.asarray()" 将列表转换后再传入。 - nvd

Answer 4

如果a的值为0，Jim Brissom的答案将失败：

import numpy  
a2 = [1, 2, 1, 1, 0, -3, -4, 7, 8, 9, 10, -2, 1, -3, 5, 6, 7, -10]  
zero_crossings2 = numpy.where(numpy.diff(numpy.sign(a2)))[0]  
print zero_crossings2  
print len(zero_crossings2)  # should be 7

输出：

[ 3  4  6 10 11 12 13 16]  
8  

零交叉的次数应该是7次，但因为sign()函数在参数为0时返回0，在参数为正值时返回1，在参数为负值时返回-1，diff()函数会将包含0的转换计算两次。

另一种替代方法可能是：

a3 = [1, 2, 1, 1, 0, -3, -4, 7, 8, 9, 10, 0, -2, 0, 0, 1, 0, -3, 0, 5, 6, 7, -10]  
s3= numpy.sign(a3)  
s3[s3==0] = -1     # replace zeros with -1  
zero_crossings3 = numpy.where(numpy.diff(s3))[0]  
print s3  
print zero_crossings3  
print len(zero_crossings3)   # should be 7

给出正确答案的是：

[ 3  6 10 14 15 18 21]
7

- Jay Borseth

---

谢谢 - 我刚看到这个答案。我想知道是否有一种简单的方法来知道零点穿越（超过0或低于0）的符号？斜率可能会有所帮助。 - Amelio Vazquez-Reina

这并没有考虑到0前后的元素符号相同的情况。 - skyork

不要使用numpy.sign，它返回负数、零或正数的-1、0或1。你应该使用numpy.where(numpy.diff(a2 > 0))[0]。或者使用Dominik Neise的答案，np.signbit。 - IceArdor

不幸的是，这个解决方案不能与其他Python容器类型一起使用，例如dique。然而，另一个解决方案可以。 - Serge Stroobandt

Answer 5

另一种适用于某些应用程序的方法是扩展表达式 np.diff(np.sign(a)) 的评估。如果我们比较这个表达式对于某些情况的反应：1. 上升穿越而没有零点：np.diff(np.sign([-10, 10])) 返回 array([2])；2. 上升穿越并带有零点：np.diff(np.sign([-10, 0, 10])) 返回 array([1, 1])；3. 下降穿越而没有零点：np.diff(np.sign([10, -10])) 返回 array([-2])；4. 下降穿越并带有零点：np.diff(np.sign([10, 0, -10])) 返回 array([-1, -1])。因此，我们必须对 1 和 2 中返回的模式评估 np.diff(...)。

sdiff = np.diff(np.sign(a))
rising_1 = (sdiff == 2)
rising_2 = (sdiff[:-1] == 1) & (sdiff[1:] == 1)
rising_all = rising_1
rising_all[1:] = rising_all[1:] | rising_2

对于第3和第4个案例：

falling_1 = (sdiff == -2) #the signs need to be the opposite
falling_2 = (sdiff[:-1] == -1) & (sdiff[1:] == -1)
falling_all = falling_1
falling_all[1:] = falling_all[1:] | falling_2

在此之后，我们可以轻松地找到具有以下指标的元素

indices_rising = np.where(rising_all)[0]
indices_falling = np.where(falling_all)[0]
indices_both = np.where(rising_all | falling_all)[0]

这种方法应该相当快，因为它可以不使用“慢”循环来管理。

这结合了其他答案的方法。

- lagru

Answer 6

我看到很多人在解决问题时经常使用diff，但是异或运算似乎更快，并且对于布尔值的结果相同（一个好的指针也可能是使用diff会产生过时警告.... :)） 以下是一个例子：

positive = a2 > 0
np.where(np.bitwise_xor(positive[1:], positive[:-1]))[0]

timeit 测量显示 diff 对我来说快了一倍半 :)

如果您不关心边缘情况，最好使用。

positive = np.signbit(a2)

但是正数 = a2 >0 似乎比使用 signbit 和检查 0（例如，positive = np.bitwise_or(np.signbit(a2),np.logical_not(a2))）更快（也更简洁）。

- ntg

Answer 7

使用带有移位数组的逐元素乘法应该是最快的：

X = np.array([ -7,   5,  -9,   4, -10,   6,   3,   3,  -5,   5])
sign_changes = np.signbit(X[1:]*X[:-1]) 
#Prepend 0 to get array of the same size
sign_changes = np.insert(sign_changes, 0, 0)

- Spas