提高重复检查的效率

Question

我是一名研究人员，使用Python处理气候模型输出，以发现某些类型的风暴。我有8个大的numpy阵列（尺寸是109574 x 52 x 57）。这些数组用1表示那天有风暴（第一维度是时间），0表示没有风暴。另外两个维度是纬度和经度。你知道吗

我必须从这些阵列中消除连续的几天。例如，如果第1天和第2天有风暴，我只想数一次风暴。如果第1天、第2天和第3天有风暴，我只想数1和3，总共有两个风暴，第1-4天有两个风暴，以此类推。我发现风暴的最后np.总和沿时间轴计算数组中的1。你知道吗

我运行以下代码来实现这一点，但我面临的问题是，它是非常缓慢的。因为我将不得不对其他数据集重复这个过程，所以我想知道是否有一种方法可以加快这个过程以提高效率。我下面有我的代码，我非常乐意澄清任何事情。你知道吗

# If there is a storm that overlaps two two-day periods, only count it once
print("Eliminating doubles...")
for i in range(52):
    for j in range(57):
        print(i,j)
        for k in range(109573):
            if((storms1[k,i,j]) == 1 and (storms1[k+1,i,j] == 1)):
                storms1[k,i,j] = 0
            if((storms2[k,i,j]) == 1 and (storms2[k+1,i,j] == 1)):
                storms2[k,i,j] = 0
            if((storms3[k,i,j]) == 1 and (storms3[k+1,i,j] == 1)):
                storms3[k,i,j] = 0
            if((storms4[k,i,j]) == 1 and (storms4[k+1,i,j] == 1)):
                storms4[k,i,j] = 0
            if((storms5[k,i,j]) == 1 and (storms5[k+1,i,j] == 1)):
                storms5[k,i,j] = 0
            if((storms6[k,i,j]) == 1 and (storms6[k+1,i,j] == 1)):
                storms6[k,i,j] = 0
            if((storms7[k,i,j]) == 1 and (storms7[k+1,i,j] == 1)):
                storms7[k,i,j] = 0
            if((storms8[k,i,j]) == 1 and (storms8[k+1,i,j] == 1)):
                storms8[k,i,j] = 0

在有人建议使用循环遍历数组之前，为了提出这个问题，我更改了变量名以简化它们。你知道吗

谢谢你的帮助。你知道吗

Answer 1

使用模拟第一个轴的一维数组的示例。首先，找出1的组从哪里开始。接下来，找出每组的长度。最后，根据您的逻辑计算事件数：

import numpy

a = numpy.random.randint(0,2,20)

# Add an initial 0
a1 = numpy.r_[0, a]

# Mark the start of each group of 1's
d1 = numpy.diff(a1) > 0

# Indices of the start of groups of 1's
w1 = numpy.arange(len(d1))[d1]

# Length of each group
cs = numpy.cumsum(a)
c = numpy.diff(numpy.r_[cs[w1], cs[-1]+1])

# Apply the counting logic
storms = c - c//2

print(a)
>>> array([0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1])
print(c)
>>> array([1, 2, 4, 1, 3])
print(storms)
>>> array([1, 1, 2, 1, 2])

通过在不再需要变量名之后重用它们，可以比我在这里展示的节省更多的内存，等等

Answer 2

以下是一个矢量化函数，可以替换最内部的循环：

def do(KK):
    # find stretches of ones
    switch_points = np.where(np.diff(np.r_[0, KK, 0]))[0]
    switch_points.shape = -1, 2
    # isolate stretches starting on odd days and create mask
    odd_starters = switch_points[switch_points[:, 0] % 2 == 1, :]
    odd_mask = np.zeros((KK.shape[0] + 1,), dtype=KK.dtype)
    odd_mask[odd_starters] = 1, -1
    odd_mask = np.add.accumulate(odd_mask[:-1])
    # apply global 1,0,1,0,1,0,... mask
    KK[1::2] = 0
    # invert stretches starting on odd days
    KK ^= odd_mask

从外部一对循环（i和j）中调用它：

do(storms1[:, i, j])
do(storms2[:, i, j])
etc.

它将改变现有的阵列。你知道吗

这应该比循环快得多（两个外部循环没有区别）。你知道吗

工作原理：

它找到一个块的起点和终点。我们知道，在每一个这样的块中，每一个都必须归零。 使用全局1,0,1,0,1,0，。。。屏蔽算法每隔一天归零一次。你知道吗

产生

• 从偶数天开始的块中的正确结果
• 街区外无变化
• 以及从奇数天开始的正确模式的补码

算法的最后一步是反转这些奇数起始块。你知道吗

Answer 3

所以我想你想要：

storms_in[:,i,j] = [0,0,1,1,0,1,1,1,0,1,0,1,1,1,1,0]
storms_out[:,i,j]= [0,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0]

这不是您的代码示例所做的，而是您在第二段中所说的要做的。你知道吗

要做到这一点，你需要两个步骤

def storms_disc(storms):  # put the whole array here, boolean-safe
    z = np.zeros((1,) + storms.shape[1:]) # zero-pads for the ends
    changes = np.r_[storms.astype('int8') ,z] - np.r_[z, storms.astype('int8')]  #find where the weather changes
    changes=((changes[:-1] == 1) | (changes[1:] == -1)).astype('int8') # reduce dimension
    return ((np.r_[changes, z] - np.r_[z, changes])[:-1] == 1).astype(storms.dtype) #find the first of successive changes

它将整个过程矢量化，您只需调用它8次。astype调用是因为减去布尔值会导致错误，即使它们的值是1和0

测试：

storms=np.random.randint(0,2,90).reshape(10,3,3)
storms.T

array([[[1, 0, 0, 1, 1, 1, 1, 1, 1, 0],
        [0, 0, 1, 1, 0, 1, 1, 0, 0, 1],
        [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]],

       [[0, 0, 0, 1, 0, 1, 0, 0, 0, 0],
        [0, 1, 0, 0, 1, 1, 1, 0, 0, 0],
        [0, 1, 0, 0, 1, 0, 1, 0, 1, 1]],

       [[0, 1, 0, 1, 0, 1, 1, 0, 0, 0],
        [0, 1, 0, 1, 0, 1, 0, 0, 1, 1],
        [0, 0, 0, 1, 1, 1, 0, 0, 1, 0]]], dtype=int8)

storms_disc(storms).T

array([[[1, 0, 0, 1, 0, 0, 0, 0, 1, 0],
        [0, 0, 1, 0, 0, 1, 0, 0, 0, 1],
        [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]],

       [[0, 0, 0, 1, 0, 1, 0, 0, 0, 0],
        [0, 1, 0, 0, 1, 0, 1, 0, 0, 0],
        [0, 1, 0, 0, 1, 0, 1, 0, 1, 0]],

       [[0, 1, 0, 1, 0, 1, 0, 0, 0, 0],
        [0, 1, 0, 1, 0, 1, 0, 0, 1, 0],
        [0, 0, 0, 1, 0, 1, 0, 0, 1, 0]]], dtype=int8)