Python/Scipy:查找“有界"矩阵的最小值/最大值 [英] Python/Scipy: Find "bounded" min/max of a matrix
问题描述
我认为最容易说明我的问题，一般情况很难解释.
假设我有一个矩阵
a 尺寸为 NxMxT，在哪里可以将 T 视为时间维度(使问题更容易).令 (n,m) 为通过 NxM 的索引.我可能将 (n,m) 称为状态空间标识符.然后我需要找到相当于
的python/scipy对于每个 (n,m):找到 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).所以在这种情况下，T = 3.# 预期的输出将是那种# b.shape = (1L, 1L, 2L, 2L, 1L, 1L, 10L)，这就解决了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]
因此，鉴于前一个数组是一个简单的堆栈 if [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.
希望能稍微澄清一下吗?
/edit2:
非常感谢答案，它有效.如果我想要相反的东西，即较小的那些中最大的，我将如何调整它?我没有尝试通过那个复杂的索引逻辑，所以我只改变前三行的(弱)尝试没有成功:
b = sort(a[...,:-1], axis=-1)b = b[...,::-1]掩码 = b <[..., -1:]索引 = argmax(掩码，轴=-1)索引 = 元组([arange(j) for j in a.shape[:-1]])指数 = 网格(*指数，索引='ij'，稀疏=真)索引.附加(索引)指数 = 元组(指数)a[指数]此外，我的第二次尝试 [...,::-1][indices] 也没有结果.
我认为 E 先生走在正确的轨道上.您肯定会首先对没有最后时间值的数组进行排序:
b = np.sort(a[..., :-1], axis=-1)您现在最好使用 `np.searchsorted 来查找大于最终值的第一项的位置，但不幸的是 np.searchsorted 仅适用于扁平数组，所以我们必须做更多的工作，比如创建一个布尔掩码，然后使用 np.argmax 找到第一个 True:
mask = b >[..., -1:]index = np.argmax(掩码，轴=-1)您现在有了索引，要提取实际值，您需要做一些索引魔术:
indices = tuple([np.arange(j) for j in b.shape[:-1]])指数 = np.meshgrid(*indices, indexing='ij', sparse=True)索引.附加(索引)指数 = 元组(指数)你现在终于可以做到了:
<预><代码>>>>b[指数]数组([[[[[[[ 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[指数].形状(1L, 1L, 2L, 2L, 1L, 1L, 10L)<小时>
要在较小的那些中获得最大的，您可以执行以下操作:
mask = b >= a[..., -1:]index = np.argmax(mask,axis=-1) - 1即较小的项中最大的是相等或更大的项中最小项之前的项.第二种情况更清楚地表明，如果没有满足条件的项目，这种方法会给出垃圾结果.在第二种情况下，当发生这种情况时，您将获得索引的 -1，因此您可以检查结果是否有效 np.any(index == -1).
如果第一种情况不能满足条件，您可以将索引设置为-1
mask = b >[..., -1:]错误 = np.all(~mask,axis=-1)index = np.argmax(掩码，轴=-1)索引[错误] = -1I think it is easiest to specify my problem, the generalized case is difficult to explain.
Say I have a matrix
a with dimensions NxMxT,
where one can think about T as a time-dimension (to make the question easier). Let (n,m) be the indices through NxM. I might call (n,m) the state-space identifier. Then I need to find the python/scipy equivalent of
for each (n,m):
find a*(n,m) = min(a(n,m,:) s.t. a*(n,m) > a(n,m,T)
That is, find the smallest state-space value that is still higher than the last (among time dimension) observation - for the whole state-space.
My first attempt was to first solve the inner problem (find a that is higher than a[...,-1]):
aHigherThanLast = a[ a > a[...,-1][...,newaxis] ]
And then I wanted to find the smallest among all of these for each (n,m). Unfortunately, aHigherThanLast now contains a 1-D array of all these values, so I don't have the (n,m) correspondence anymore. What would be a better approach to this?
As an additional problem: The state-space is variable, it could also be 3 or more dimensions (NxMxKx...), and I cannot hard-code this. So any kind of
for (n,m,t) in nditer(a):
is not feasible.
Many thanks!
/edit:
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 is higher than the newest observation)
There is no element in i in a that satisfies both
-- 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 is the smallest element that is higher than the newest observation)
So, given that the previous array is a simple stack if [0,2,1] along the last observation, I would expect
b = ones((1,1,2,2,1,1,10))*2
however, - if, among some (a,b,c,d,e,f,g), there was not only the value of either {0,1,2}, but also {3}, then I would still want the 2 (as it is the smaller of i = {2,3} that satisfies i > 1. - if among some (a,b,c,d,e,f,g) there was only the value {0,1,3}, I would want the 3, as i = 3 would be the smallest number that satisfies i > 1.
Hope that cleared it up a bit?
/edit2:
Appreciate the answer a lot, it works. How would I adjust it if I wanted the opposite, i.e. the largest among those that are smaller? I didn't try to get through that complicated indexing logic, so my (weak) attempt of only changing the first three lines did not succeed:
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]
Also, a[...,::-1][indices], my second attempt, was not fruitful either.
I think Mr. E is on the right track. You definitely start by sorting the array without that last time value:
b = np.sort(a[..., :-1], axis=-1)
You would now ideally use `np.searchsorted to find where the first item larger than the end value is, but unfortunately np.searchsorted only works on flattened arrays, so we have to do some more work, like creating a boolean mask, then finding the first True using np.argmax:
mask = b > a[..., -1:]
index = np.argmax(mask, axis=-1)
You now have the indices, to extract the actual values, you need to do some indexing magic:
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)
And you can now finally do:
>>> 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)
To get the largest among those that are smaller, you could do something like:
mask = b >= a[..., -1:]
index = np.argmax(mask, axis=-1) - 1
i.e. the largest among those that are smaller is, the item right before the smallest among those that are equal or larger. This second case makes it more clear that this approach gives garbage result if there is no item that fulfills the condition. In this second case, when that happens, you will get a -1 for an index, so you could check that the results are valid doing np.any(index == -1).
You can set the index to -1 if the condition cannot be satisfied for the first case by doing
mask = b > a[..., -1:]
wrong = np.all(~mask, axis=-1)
index = np.argmax(mask, axis=-1)
index[wrong] = -1
这篇关于Python/Scipy:查找“有界"矩阵的最小值/最大值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
