4

我认为最容易说明我的问题，笼统的情况很难解释。

说我有一个矩阵

a with dimensions NxMxT,

可以将 T 视为时间维度（使问题更容易）。令 (n,m) 为通过 NxM 的索引。我可以称 (n,m) 为状态空间标识符。然后我需要找到 python/scipy 等价物

for each (n,m):
     find a*(n,m) = min(a(n,m,:) s.t. a*(n,m) > a(n,m,T)

也就是说，对于整个状态空间，找到仍然高于最后（时间维度中）观察值的最小状态空间值。

我的第一次尝试是首先解决内部问题（找到高于 a[...,-1] 的 a）：

aHigherThanLast = a[ a > a[...,-1][...,newaxis] ]

然后我想为每个 (n,m) 找到所有这些中最小的。不幸的是，aHigherThanLast 现在包含所有这些值的一维数组，所以我不再有 (n,m) 对应关系。有什么更好的方法来解决这个问题？

作为一个额外的问题：状态空间是可变的，它也可能是 3 维或更多维（NxMxKx ...），我无法对此进行硬编码。所以任何一种

for (n,m,t) in nditer(a):

是不可行的。

非常感谢！

/编辑：

a = array([[[[[[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]],



          [[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]]],




         [[[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]],



          [[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]]]]]])
# a.shape = (1L, 1L, 2L, 2L, 1L, 1L, 10L, 3L). so in this case, T = 3.
# expected output would be the sort of
# b.shape = (1L, 1L, 2L, 2L, 1L, 1L, 10L), which solves

  • b[a,b,c,d,e,f,g] > a[a,b,c,d,e,f,g,-1] （b 高于最新观察值）

    • a 中的 i 中没有任何元素同时满足

      -- a[a,b,c,d,e,f,g,t] > a[a,b,c,d,e,f,g,-1]

      -- a[a,b,c,d,e,f,g,t] < b[a,b,c,d,e,f,g] （b 是高于最新观测值的最小元素)

所以，假设前面的数组是一个简单的堆栈，如果 [0,2,1] 沿着最后一次观察，我希望

b = ones((1,1,2,2,1,1,10))*2

但是，-如果在一些 (a,b,c,d,e,f,g) 中，不仅有 {0,1,2} 的值，还有 {3}，那么我仍然想要2（因为它是满足 i > 1 的 i = {2,3} 中的较小者。 - 如果在一些 (a,b,c,d,e,f,g) 中只有值 {0,1 ,3}，我想要 3，因为 i = 3 将是满足 i > 1 的最小数字。

希望能澄清一点？

/编辑2：

非常感谢答案，它有效。如果我想要相反的情况，即较小的那些中最大的，我将如何调整它？我没有尝试通过那个复杂的索引逻辑，所以我（弱）尝试只更改前三行没有成功：

        b = sort(a[...,:-1], axis=-1)
        b = b[...,::-1]
        mask = b < a[..., -1:]
        index = argmax(mask, axis=-1)
        indices = tuple([arange(j) for j in a.shape[:-1]])
        indices = meshgrid(*indices, indexing='ij', sparse=True)
        indices.append(index)
        indices = tuple(indices)
        a[indices]

此外，我的第二次尝试 a[...,::-1][indices] 也没有结果。

4

1 回答 1

2

我认为 E 先生是在正确的轨道上。您肯定首先对没有最后一次值的数组进行排序：

b = np.sort(a[..., :-1], axis=-1)

您现在最好使用`np.searchsorted查找大于最终值的第一项在哪里，但不幸的是np.searchsorted仅适用于扁平数组，因此我们必须做更多的工作，例如创建一个布尔掩码，然后找到第一个True使用np.argmax：

mask = b > a[..., -1:]
index = np.argmax(mask, axis=-1)

您现在有了索引，要提取实际值，您需要做一些索引魔术：

indices = tuple([np.arange(j) for j in b.shape[:-1]])
indices = np.meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)

你现在终于可以做到了：

>>> b[indices]
array([[[[[[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]],


          [[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]]],



         [[[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]],


          [[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]]]]]])
>>> b[indices].shape
(1L, 1L, 2L, 2L, 1L, 1L, 10L)

要在较小的那些中获得最大的，您可以执行以下操作：

mask = b >= a[..., -1:]
index = np.argmax(mask, axis=-1) - 1

即那些较小的项目中最大的是，在相等或更大的那些中最小的项目之前。第二种情况更清楚地表明，如果没有满足条件的项目，这种方法会给出垃圾结果。在第二种情况下，当发生这种情况时，您将获得一个-1for 索引，因此您可以检查结果是否有效np.any(index == -1)。

如果第一种情况不能满足条件，您可以将索引设置为 -1

mask = b > a[..., -1:]
wrong = np.all(~mask, axis=-1)
index = np.argmax(mask, axis=-1)
index[wrong] = -1
于 2013-11-12T14:43:29.553 回答