5

我有一个数组，我想查看该数组中的任何元素是否大于或等于该数组中的任何其他元素。我可以做两个 for 循环，但我的数组长度为 10,000 或更大，因此创建了一个非常慢的程序。无论如何，我可以更快地做到这一点吗？

[编辑] 我只需要它来查看它是否大于或等于我正在查看的元素之后的元素，如果是，我需要知道它的索引。

[编辑] 我将更清楚地重新解释我的问题，因为当前的解决方案不能满足我的需要。首先，这里有一些代码

x=linspace(-10, 10, 10000)
t=linspace(0,5,10000)

u=np.exp(-x**2)

k=u*t+x

所以我取一个 x 数组，通过将它放入高斯得到它的高度，然后基于该高度，这就是 x 值在空间中传播的速度，我用 k 找到。我的问题是，我需要找出高斯何时变为双值函数（或者换句话说，何时发生冲击）。如果我执行 argmax 解决方案，我将始终获得 k 中的最后一个值，因为它非常接近于零，我需要元素之后的第一个值，这将在我的函数中给我一个 double 值。

[编辑] 小例子

x=[0,1,2,3,4,5,6,7,8,9,10] #Input 
k=[0,1,2,3,4,5,6,5,4,10] #adjusted for speed

output I want
in this case, 5 is the first number that goes above a number that comes after it.
So I need to know the index of where 5 is located and possibly the index 
of the number that it is greater than
4

3 回答 3

5

大于后一个值的第一个值必然对应于局部最小值中的最小值：

k = np.array([0,1,2,3,4,5,6,5,4,10])
lm_i = np.where(np.diff(np.sign(np.diff(k))) > 0)[0] + 1
mlm = np.min(k[lm_i])
mlm_i = lm_i[np.argmin(k[lm_i])]

大于后一个值的第一个值的索引是大于该最小局部最小值的第一个索引：

i = np.where(k > mlm)[0][0]

解图

（忽略图形似乎与切线处的水平线相交；这只是一个显示伪影。）

作为一个单行：

np.where(k > np.min(k[np.where(np.diff(np.sign(np.diff(k))) > 0)[0] + 1]))[0][0]

请注意，这是大约。比 root 的解决方案快 1000 倍，因为它是完全矢量化的：

%timeit np.where(k > np.min(k[np.where(np.diff(np.sign(np.diff(k))) > 0)[0] + 1]))[0][0]
1000 loops, best of 3: 228 us per loop
于 2013-03-11T16:35:46.740 回答
3

矢量化解决方案，比 ecatmur 的解决方案快 25%：

np.where(k > np.min(k[np.where(np.diff(k) < 0)[0][0]:]))[0][0]

一种天真的方法：

next(i for i in np.arange(len(arr)) if arr[i:].argmin() != 0)
于 2013-03-11T15:09:31.833 回答
2

编辑 实际上，拥有 10,000 项 python for 循环比在 100,000,000 项数组上操作更便宜::

In [14]: np.where(np.array([True if np.all(k[:j] <= k[j]) else
                            False for j in xrange(len(k))]) == 0)
Out[14]: (array([5129, 5130, 5131, ..., 6324, 6325, 6326]),)

In [15]: %timeit np.where(np.array([True if np.all(k[:j] <= k[j]) else
                                    False for j in xrange(len(k))]) == 0)
1 loops, best of 3: 201 ms per loop

就内存而言，这将是昂贵的，但您可以使用广播对搜索进行矢量化。如果你这样做：

>>> k <= k[:, None]
array([[ True, False, False, ..., False, False, False],
       [ True,  True, False, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       ..., 
       [ True,  True,  True, ...,  True, False, False],
       [ True,  True,  True, ...,  True,  True, False],
       [ True,  True,  True, ...,  True,  True,  True]], dtype=bool)

返回值是一个布尔数组，其中位置项[i, j]告诉您是否k[j]小于或等于k[i]。什么时候可以使用np.cumprod如下：

>>> np.cumprod(k <= k[:, None], axis=1)
array([[1, 0, 0, ..., 0, 0, 0],
       [1, 1, 0, ..., 0, 0, 0],
       [1, 1, 1, ..., 0, 0, 0],
       ..., 
       [1, 1, 1, ..., 1, 0, 0],
       [1, 1, 1, ..., 1, 1, 0],
       [1, 1, 1, ..., 1, 1, 1]])

位置中的项目[i, j]告诉您是否k[j]小于或等于 中的所有项目k[:i]。如果您采用该矩阵的对角线：

>>> np.cumprod(k <= k[:, None], axis=1)[np.diag_indices(k.shape[0])]
array([1, 1, 1, ..., 1, 1, 1])

位置的项目[i]告诉你是否k[i]小于或等于它之前的所有项目。查找该数组为零的位置：

>>> np.where(np.cumprod(k <= k[:, None],
...                     axis=1)[np.diag_indices(k.shape[0])] == 0)
(array([5129, 5130, 5131, ..., 6324, 6325, 6326]),)

并且您将拥有满足您所需条件的所有值的索引。

如果您只对第一个感兴趣：

>>> np.argmax(np.cumprod(k <= k[:, None],
...                      axis=1)[np.diag_indices(k.shape[0])] == 0)
5129

这不是一个轻量级的操作，但如果你有足够的内存来容纳所有的布尔数组，它不会让你等待太久：

In [3]: %timeit np.argmax(np.cumprod(k <= k[:, None],
                                     axis=1)[np.diag_indices(k.shape[0])] == 0)
1 loops, best of 3: 948 ms per loop
于 2013-03-11T16:23:39.040 回答