Подтвердить что ты не робот

Быстрая реализация в Mathematica для Position2D

Я ищу быстрое выполнение для следующего, я назову его Position2D из-за отсутствия лучшего термина:

Position2D[ matrix, sub_matrix ]

который находит расположение sub_matrix внутри matrix и возвращает верхнюю левую и нижнюю правую строку/столбец совпадения.

Например, это:

Position2D[{
 {0, 1, 2, 3},
 {1, 2, 3, 4},
 {2, 3, 4, 5},
 {3, 4, 5, 6}
}, {
 {2, 3},
 {3, 4}
}]

должен вернуть это:

{
 {{1, 3}, {2, 4}},
 {{2, 2}, {3, 3}},
 {{3, 1}, {4, 2}} 
}

Он должен быть достаточно быстрым, чтобы быстро работать с матрицами 3000x2000 с субматрицами 100x100. Для простоты достаточно рассматривать только целые матрицы.

4b9b3361

Ответ 1

Алгоритм

Следующий код основан на эффективной пользовательской позиционной функции для нахождения позиций (возможно перекрывающихся) целых последовательностей в большом целочисленном списке. Основная идея заключается в том, что мы можем сначала попытаться найти позиции, в которых первая строка подматрицы находится в большой матрице Flatten -ed, а затем отфильтровывает их, извлекает полные подматрицы и сравнивает их с суб- представляющая интерес матрица. Это будет эффективно для большинства случаев, за исключением очень патологических (те, для которых эта процедура будет генерировать огромное количество потенциальных кандидатов в потенциальных кандидатах, в то время как истинное количество записей в подматрице будет намного меньше, но такие случаи кажутся маловероятными как правило, и затем можно сделать дополнительные усовершенствования этой простой схемы).

Для больших матриц предлагаемое решение будет примерно в 15-25 раз быстрее, чем решение @Szabolcs, когда используется скомпилированная версия функции позиционирования последовательностей и в 3-5 раз быстрее для реализации положений последовательности на верхнем уровне - функция поиска. Фактическое ускорение зависит от размеров матрицы, больше для больших матриц. Код и контрольные показатели приведены ниже.

код

Обычно эффективная функция для нахождения позиций суб-списка (последовательности)

Эти вспомогательные функции связаны с Норбертом Позаром и взяты из этого потока Mathgroup. Они используются для эффективного поиска стартовых позиций целочисленной последовательности в более крупном списке (подробнее см. Упомянутую запись).

Clear[seqPos];
fdz[v_] := [email protected]@Prepend[v, 0];
seqPos[list_List, seq_List] :=
  Fold[
     fdz[#1 (1 - Unitize[list[[#1]] - #2])] + 1 &, 
     fdz[Range[Length[list] - Length[seq] + 1] *
       (1 - Unitize[list[[;; -Length[seq]]] - seq[[1]]])] + 1,
     [email protected]
  ] - Length[seq];

Пример использования:

In[71]:= seqPos[{1,2,3,2,3,2,3,4},{2,3,2}]
Out[71]= {2,4}

Более быстрая функция определения местоположения для целых чисел

Однако быстрая seqPos может быть, это все еще является основным узким местом в моем решении. Вот скомпилированная версия этого кода, которая дает еще 5-кратное повышение производительности моего кода:

seqposC = 
  Compile[{{list, _Integer, 1}, {seq, _Integer, 1}},
    Module[{i = 1, j = 1, res = Table[0, {Length[list]}], ctr = 0},
       For[i = 1, i <= Length[list], i++,
          If[list[[i]] == seq[[1]],
             While[j < Length[seq] && i + j <= Length[list] && 
                     list[[i + j]] == seq[[j + 1]], 
                j++
             ];
             If[j == Length[seq], res[[++ctr]] = i];
             j = 1;
          ]
       ];
       Take[res, ctr]
    ], CompilationTarget -> "C", RuntimeOptions -> "Speed"]

Пример использования:

In[72]:= seqposC[{1, 2, 3, 2, 3, 2, 3, 4}, {2, 3, 2}]
Out[72]= {2, 4}

Приведенные ниже тесты были переделаны с помощью этой функции (также код для основной функции слегка изменен)

Основная функция

Это основная функция. Он находит позиции первой строки в матрице и затем фильтрует их, извлекает подматрицы в этих положениях и тестирует на полную интересующую субматрицу:

Clear[Position2D];
Position2D[m_, what_,seqposF_:Automatic] :=
  Module[{posFlat, pos2D,sp = If[seqposF === Automatic,seqposC,seqposF]},
     With[{dm = Dimensions[m], dwr = [email protected][what]},
         posFlat = sp[[email protected], [email protected]];
         pos2D = 
            Pick[Transpose[#], Total[Clip[[email protected] - # - dwr + 2, {0, 1}]],2] &@
                {Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm];
         Transpose[{#, Transpose[Transpose[#] + dwr - 1]}] &@
            Select[pos2D,
                m[[[email protected]# ;; [email protected]# + [email protected] - 1, 
                   [email protected]# ;; [email protected]# + [email protected] - 1]] == what &
            ]
     ]
  ];

Для целых списков можно использовать более быструю скомпилированную подпоследовательную функцию определения местоположения seqposC (это значение по умолчанию). Для общих списков можно указать, например. seqPos, в качестве третьего аргумента.

Как это работает

Мы будем использовать простой пример для анализа кода и объяснения его внутренних действий. Это определяет нашу тестовую матрицу и подматрицу:

m  = {{0, 1, 2, 3}, {1, 2, 3, 4}, {2, 3, 4, 5}};
what = {{2, 3}, {3, 4}}; 

Это вычисляет размеры вышеизложенного (удобнее работать с измененными размерами для подматрицы):

In[78]:= 
dm=Dimensions[m]
[email protected][what]

Out[78]= {3,4}
Out[79]= {2,2}

Здесь найден список начальных позиций первой строки ({2,3} здесь) в основной матрице Flatten ed. Эти позиции в то же время являются "плоскими" позициями кандидата в верхнем левом углу подматрицы:

In[77]:= posFlat = seqPos[[email protected], [email protected]]
Out[77]= {3, 6, 9}

Это восстановит 2D-кандидатские позиции верхнего левого угла подматрицы в полной матрице, используя размеры основной матрицы:

In[83]:= posInterm = [email protected]{Mod[posFlat,#,1],IntegerPart[posFlat/#]+1}&@Last[dm]
Out[83]= {{3,1},{2,2},{1,3}}

Затем мы попробуем использовать Select, чтобы отфильтровать их, извлекая полную подматрицу и сравнив с what, но мы столкнемся с проблемой здесь:

In[84]:= 
Select[posInterm,
   m[[[email protected]#;;[email protected]#[email protected],[email protected]#;;[email protected]#[email protected]]]==what&]

During evaluation of In[84]:= Part::take: Cannot take positions 3 through 4 
in {{0,1,2,3},{1,2,3,4},{2,3,4,5}}. >>

Out[84]= {{3,1},{2,2}}

Помимо сообщения об ошибке, результат будет правильным. Само сообщение об ошибке связано с тем, что для последней позиции ({1,3}) в списке правый нижний угол подматрицы будет находиться за пределами основной матрицы. Разумеется, мы могли бы использовать Quiet просто игнорировать сообщения об ошибках, но это плохой стиль. Итак, мы сначала отфильтровываем эти случаи, и для этого предназначена строка Pick[Transpose[#], Total[Clip[[email protected] - # - dwr + 2, {0, 1}]], 2] &@. В частности, рассмотрим

In[90]:= 
[email protected] - # - dwr + 2 &@{Mod[posFlat, #, 1],IntegerPart[posFlat/#] + 1} &@Last[dm]
Out[90]= {{1,2,3},{2,1,0}}

Координаты верхних левых углов должны оставаться в пределах разности размеров матрицы и подматрицы. Вышеуказанные подписи были сделаны из координаты x и y верхних левых углов. Я добавил 2, чтобы все достоверные результаты были строго положительными. Мы должны выбрать только координаты в тех положениях в [email protected]{Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm] (что posInterm), при которых оба подписок выше имеют строго положительные числа. Я использовал Total[Clip[...,{0,1}]], чтобы переделать его в подборку только в тех положениях, в которых этот второй список имеет 2 (Clip преобразует все положительные целые числа в 1, а Total суммирует числа в 2 подсписках. get 2 - когда числа в обоих подсписках положительны).

Итак, мы имеем:

In[92]:= 
pos2D=Pick[Transpose[#],Total[Clip[[email protected]#-dwr+2,{0,1}]],2]&@
           {Mod[posFlat,#,1],IntegerPart[posFlat/#]+1}&@Last[dm]
Out[92]= {{3,1},{2,2}} 

После того, как список 2D-позиций был отфильтрован, поэтому нет структурно-недопустимых позиций, мы можем использовать Select для извлечения полных подматриц и теста против интересующей подматрицы:

In[93]:= 
finalPos = 
  Select[pos2D,m[[[email protected]#;;[email protected]#[email protected],[email protected]#;;[email protected]#[email protected]]]==what&]
Out[93]= {{3,1},{2,2}}

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

In[94]:= Transpose[{#,Transpose[Transpose[#]+dwr-1]}]&@finalPos 
Out[94]= {{{3,1},{4,2}},{{2,2},{3,3}}}

Я мог бы использовать Map, но для большого списка позиций вышеуказанный код был бы более эффективным.

Пример и тесты

Исходный пример:

In[216]:= Position2D[{{0,1,2,3},{1,2,3,4},{2,3,4,5},{3,4,5,6}},{{2,3},{3,4}}]
Out[216]= {{{3,1},{4,2}},{{2,2},{3,3}},{{1,3},{2,4}}}

Обратите внимание, что мои условные обозначения обратимы w.r.t. Решение @Szabolcs.

Контрольные показатели для больших матриц и подматриц

Вот тест мощности:

nmat = 1000;
(* generate a large random matrix and a sub-matrix *)
largeTestMat = RandomInteger[100, {2000, 3000}];
what = RandomInteger[10, {100, 100}];
(* generate upper left random positions where to insert the submatrix *)    
rposx = RandomInteger[{1,[email protected][largeTestMat] - [email protected][what] + 1}, nmat];
rposy = RandomInteger[{1,[email protected][largeTestMat] - [email protected][what] + 1},nmat];
(* insert the submatrix nmat times *)
With[{dwr = [email protected][what]},
    Do[largeTestMat[[[email protected] ;; [email protected] + [email protected] - 1, 
                     [email protected] ;; [email protected] + [email protected] - 1]] = what, 
       {p,Transpose[{rposx, rposy}]}]]

Теперь мы проверяем:

In[358]:= (ps1 = position2D[largeTestMat,what])//Short//Timing
Out[358]= {1.39,{{{1,2461},{100,2560}},<<151>>,{{1900,42},{1999,141}}}}

In[359]:= (ps2 = Position2D[largeTestMat,what])//Short//Timing
Out[359]= {0.062,{{{2461,1},{2560,100}},<<151>>,{{42,1900},{141,1999}}}}

(фактическое число подматриц меньше числа, которое мы пытаемся сгенерировать, поскольку многие из них перекрываются и "уничтожают" ранее вставленные - это происходит потому, что размер подматрицы является значительной частью размер матрицы в нашем тесте).

Чтобы сравнить, мы должны отменить индексы x-y в одном из решений (уровень 3) и отсортировать оба списка, так как позиции могут быть получены в другом порядке:

In[360]:= [email protected]1===Sort[Reverse[ps2,{3}]]
Out[360]= True

Я не исключаю, что возможны дальнейшие оптимизации.

Ответ 2

Это моя реализация:

position2D[m_, k_] :=
 Module[{di, dj, extractSubmatrix, pos},
  {di, dj} = Dimensions[k] - 1;
  extractSubmatrix[{i_, j_}] := m[[i ;; i + di, j ;; j + dj]];
  pos = Position[ListCorrelate[k, m], ListCorrelate[k, k][[1, 1]]];
  pos = Select[pos, extractSubmatrix[#] == k &];
  {#, # + {di, dj}} & /@ pos
 ]

Он использует ListCorrelate для получения списка потенциальных позиций, а затем фильтрует те, которые на самом деле соответствуют. Это, вероятно, быстрее на упакованных реальных матрицах.

Ответ 3

В соответствии с предложением Леонида здесь мое решение. Я знаю, что он не очень эффективен (он примерно в 600 раз медленнее, чем Леонид, когда я его приурочил), но он очень короткий, запоминающийся и хорошая иллюстрация редко используемой функции PartitionMap. Это из пакета Developer, поэтому сначала требуется вызов Needs["Developer`"].

Учитывая, что Position2D можно определить как:

Position2D[m_, k_] :=  Position[PartitionMap[k == # &, m, Dimensions[k], {1, 1}], True]

Это дает только верхние левые координаты. Я чувствую, что нижние правые координаты действительно избыточны, поскольку размеры подматрицы известны, но если возникает необходимость, можно добавить их к выводу путем добавления {#, Dimensions[k] + # - {1, 1}} & /@ к указанному выше определению.

Ответ 4

Как насчет чего-то вроде

Position2D[bigMat_?MatrixQ, smallMat_?MatrixQ] := 
 Module[{pos, sdim = Dimensions[smallMat] - 1}, 
  pos = Position[bigMat, smallMat[[1, 1]]];
  Quiet[Select[pos, (MatchQ[
      bigMat[[[email protected]@Thread[Span[#, # + sdim]]]], smallMat] &)], 
    Part::take]]

который вернет верхние левые позиции подматриц. Пример:

Position2D[{{0, 1, 2, 3}, {1, 2, 3, 4}, {2, 3, 4, 5}, {3, 5, 5, 6}}, 
   {{2, 3}, {3, _}}]
(* Returns: {{1, 3}, {2, 2}, {3, 1}} *) 

И для поиска матрицы 1000x1000 она занимает около 2 секунд на моей старой машине

SeedRandom[1]
big = RandomInteger[{0, 10}, {1000, 1000}];

Position2D[big, {{1, 1, _}, {1, 1, 1}}] // Timing
(* {1.88012, {{155, 91}, {295, 709}, {685, 661}, 
              {818, 568}, {924, 45}, {981, 613}}} *)