基于最小值/最大值的 Numpy 动态数组切片

Question

我有一个 hape (365, x, y) 的 3 维数组，其中 36 对应于 =daily 数据。在某些情况下，时间轴上的所有元素axis=0 都是np.nan。

axis=0 沿线每个点的时间序列如下所示：

我需要找到最大值（峰值数据）出现的索引，然后是峰值两侧的两个最小值。

import numpy as np

a = np.random.random(365, 3, 3) * 10
a[:, 0, 0] = np.nan

peak_mask = np.ma.masked_array(a, np.isnan(a))
peak_indexes = np.nanargmax(peak_mask, axis=0)

我可以使用以下方法找到峰值之前的最小值：

early_minimum_indexes = np.full_like(peak_indexes, fill_value=0)

for i in range(peak_indexes.shape[0]):
    for j in range(peak_indexes.shape[1]):
        if peak_indexes[i, j] == 0:
            early_minimum_indexes[i, j] = 0
        else:
            early_mask = np.ma.masked_array(a, np.isnan(a))
            early_loc = np.nanargmin(early_mask[:peak_indexes[i, j], i, j], axis=0)   
            early_minimum_indexes[i, j] = early_loc

得到的波峰和波谷如​​下所示：

对于大型数组（1m+ 个元素），这种方法在时间上是非常不合理的。有没有更好的方法来使用 numpy 做到这一点？

Answer 1

虽然在这种情况下使用掩码数组可能不是最最有效的解决方案，但它可以让您在特定轴上执行掩码操作，同时或多或少地保留形状，这非常方便.请记住，在许多情况下，被屏蔽的函数最终仍会复制被屏蔽的数据。

您当前的代码中的想法基本正确，但您错过了一些技巧，例如能够否定和组合掩码。此外，预先将掩码分配为布尔值更有效，并且像np.full(..., 0) -> np.zeros(..., dtype=bool) 这样的小挑剔。

让我们反过来解决这个问题。假设你有一个表现良好的一维数组，它有一个峰值，比如a1。您可以使用掩码轻松找到最大值和最小值（或索引），如下所示：

peak_index = np.nanargmax(a1)
mask = np.zeros(a1.size, dtype=np.bool)
mask[peak:] = True
trough_plus = np.nanargmin(np.ma.array(a1, mask=~mask))
trough_minus = np.nanargmin(np.ma.array(a1, mask=mask))

这尊重掩码数组相对于普通 numpy 布尔索引翻转掩码的意义这一事实。最大值出现在trough_plus的计算中也是可以的，因为它保证不是最小值（除非你有全南的情况）。

现在如果a1 已经是一个掩码数组（但仍然是一维数组），您可以做同样的事情，但暂时合并掩码。例如：

a1 = np.ma.array(a1, mask=np.isnan(a1))
peak_index = a1.argmax()
mask = np.zeros(a1.size, dtype=np.bool)
mask[peak:] = True
trough_plus = np.ma.masked_array(a1, mask=a.mask | ~mask).argmin()
trough_minus  (np.ma.masked_array(a1, mask=a.mask | mask).argmin()

同样，由于掩码数组具有反转掩码，因此使用 | 而不是 & 组合掩码非常重要，就像使用普通的 numpy 布尔掩码一样。在这种情况下，不需要调用 nan 版本的 argmax 和 argmin，因为所有的 nan 都已经被屏蔽掉了。

鉴于axis 关键字在 numpy 函数中的普遍存在，希望从这里对多维度的泛化变得清晰：

a = np.ma.array(a, mask=np.isnan(a))
peak_indices = a.argmax(axis=0).reshape(1, *a.shape[1:])
mask = np.arange(a.shape[0]).reshape(-1, *(1,) * (a.ndim - 1)) >= peak_indices

trough_plus = np.ma.masked_array(a, mask=~mask | a.mask).argmin(axis=0)
trough_minus = np.ma.masked_array(a, mask=mask | a.mask).argmin(axis=0)

N 维掩蔽技术来自Fill mask efficiently based on start indices，专门为此目的而询问。

Answer 2

这是一个方法

1. 复制数据
2. 保存所有 nan 位置并将所有 nan 替换为全局 min-1
3. 找到按行的 argmax
4. 从整行中减去它的值
   • 请注意，现在每一行只有非正值，最大值现在为零
5. 将所有 nan 位置归零
6. 翻转最大值右侧所有值的符号
   • 这是主要思想；它在右手最小值之前的位置创建一个新的行全局最大值；同时它确保左侧的最小值现在是行全局的
7. 检索按行的 argmin 和 argmax，这些是原始数组中左右 mins 的位置
8. 查找全南行并用 INVALINT 覆盖这些位置的最大和最小索引

代码：

INVALINT = -9999
t,x,y = a.shape
t,x,y = np.ogrid[:t,:x,:y]
inval = np.isnan(a)
b = np.where(inval,np.nanmin(a)-1,a)
pk = b.argmax(axis=0)
pkval = b[pk,x,y]
b -= pkval
b[inval] = 0
b[t>pk[None]] *= -1
ltr = b.argmin(axis=0)
rtr = b.argmax(axis=0)
del b
inval = inval.all(axis=0)
pk[inval] = INVALINT
ltr[inval] = INVALINT
rtr[inval] = INVALINT

# result is now in ltr ("left trough"), pk ("peak") and rtr