Python/Scipy: найти ограниченный минимум/максимум матрицы

Думаю, проще всего конкретизировать мою проблему, обобщенный случай объяснить сложно.

Скажем, у меня есть матрица

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]):

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 больше, чем самое новое наблюдение)

    • Нет элемента i в a, который удовлетворяет обоим

      -- 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 = {2,3}, которое удовлетворяет i > 1. - если среди некоторых (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]
        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], моя вторая попытка тоже не увенчалась успехом.


python numpy scipy maximization
person FooBar    schedule 12.11.2013    source источник
comment
Можете ли вы привести пример входных данных (в идеале в виде кода Python — жестко запрограммировать некоторые массивы) и ожидаемый результат для иллюстрации?   -  person YXD    schedule 12.11.2013
comment
Ладно, так бы и было :)   -  person FooBar    schedule 12.11.2013
comment
Хм, не могу заставить его работать полностью прямо сейчас, но попробуйте поиграть с end_slice = a[..., -1]; b = np.sort(a, axis=-1); b >= end_slice[..., None]; indices = np.argmax(b >= end_slice[..., None], axis=-1) (точки с запятой там, где у нас должны быть новые строки...)   -  person YXD    schedule 12.11.2013

Ответы (1)


arrow_upward
2
arrow_downward

Я думаю, мистер Э. на правильном пути. Вы определенно начинаете с сортировки массива без этого последнего значения времени:

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

т. е. самый большой среди меньших — это элемент, стоящий прямо перед самым маленьким среди тех, которые равны или больше. Этот второй случай делает более ясным, что этот подход дает мусорный результат, если нет элемента, удовлетворяющего условию. Во втором случае, когда это произойдет, вы получите -1 для индекса, чтобы вы могли проверить правильность результатов, выполнив np.any(index == -1).

Вы можете установить индекс равным -1, если условие не может быть выполнено для первого случая, выполнив

mask = b > a[..., -1:]
wrong = np.all(~mask, axis=-1)
index = np.argmax(mask, axis=-1)
index[wrong] = -1
person Jaime    schedule 12.11.2013
comment
Превосходно. Я думаю, что понимаю, что происходит в начале, но я просто приму всю магию индексов как должное. А если бы, напротив, у меня были самые большие числа, которые меньше, чем последнее наблюдение, я бы изменил только первые три строки следующим образом? Большое спасибо! b = sort(oldPrices[...,:-1], axis=-1) b = b[...,::-1] mask = b ‹ oldPrices[..., -1:] - person FooBar; 12.11.2013
comment
(См. также последнее добавленное мной изменение) - person FooBar; 12.11.2013
comment
Упс, я исправил ошибку: индексация должна выполняться по b, отсортированному массиву, а не по исходному массиву. Смотрите мои правки во второй части вашего вопроса, и когда этот подход терпит неудачу. - person Jaime; 12.11.2013
comment
Спасибо. Я здесь только догадываюсь, но я думаю, вы имеете в виду mask = b ‹= a[..., -1:] ? - person FooBar; 12.11.2013
comment
Могу я обратить ваше внимание на связанный с этим вопрос? stackoverflow.com/questions/24098205/ - person FooBar; 09.06.2014