Быстрое внедрение в Mathematica для Position2D - PullRequest
20 голосов
/ 03 декабря 2011

Я ищу быструю реализацию для следующего, я назову это 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 подматрицами.Для простоты достаточно рассмотреть только целочисленные матрицы.

Ответы [ 4 ]

29 голосов
/ 03 декабря 2011

Алгоритм

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

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

код

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

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

Clear[seqPos];
fdz[v_] := Rest@DeleteDuplicates@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,
     Rest@seq
  ] - Length[seq];

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

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

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

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

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 = Reverse@Dimensions[what]},
         posFlat = sp[Flatten@m, First@what];
         pos2D = 
            Pick[Transpose[#], Total[Clip[Reverse@dm - # - dwr + 2, {0, 1}]],2] &@
                {Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm];
         Transpose[{#, Transpose[Transpose[#] + dwr - 1]}] &@
            Select[pos2D,
                m[[Last@# ;; Last@# + Last@dwr - 1, 
                   First@# ;; First@# + First@dwr - 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]
dwr=Reverse@Dimensions[what]

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

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

In[77]:= posFlat = seqPos[Flatten@m, First@what]
Out[77]= {3, 6, 9}

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

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

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

In[84]:= 
Select[posInterm,
   m[[Last@#;;Last@#+Last@dwr-1,First@#;;First@#+First@dwr-1]]==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[Reverse@dm - # - dwr + 2, {0, 1}]], 2] &@. В частности, рассмотрим

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

Координаты левого верхнего угла должны оставаться в пределах разницы размеров матрицы и подматрицы. Вышеуказанные подсписки были составлены из x и y координат верхнего левого угла. Я добавил 2, чтобы сделать все действительные результаты строго положительными. Мы должны выбрать только координаты в тех позициях в Transpose@{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[Reverse@dm-#-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[[Last@#;;Last@#+Last@dwr-1,First@#;;First@#+First@dwr-1]]==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,Last@Dimensions[largeTestMat] - Last@Dimensions[what] + 1}, nmat];
rposy = RandomInteger[{1,First@Dimensions[largeTestMat] - First@Dimensions[what] + 1},nmat];
(* insert the submatrix nmat times *)
With[{dwr = Reverse@Dimensions[what]},
    Do[largeTestMat[[Last@p ;; Last@p + Last@dwr - 1, 
                     First@p ;; First@p + First@dwr - 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]:= Sort@ps1===Sort[Reverse[ps2,{3}]]
Out[360]= True

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

13 голосов
/ 03 декабря 2011

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

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, чтобы получить список потенциальных позиций, а затем фильтрует те, которые действительно соответствуют.Вероятно, быстрее на упакованных реальных матрицах.

8 голосов
/ 07 декабря 2011

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

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

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

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

6 голосов
/ 03 декабря 2011

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

Position2D[bigMat_?MatrixQ, smallMat_?MatrixQ] := 
 Module[{pos, sdim = Dimensions[smallMat] - 1}, 
  pos = Position[bigMat, smallMat[[1, 1]]];
  Quiet[Select[pos, (MatchQ[
      bigMat[[Sequence@@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}}} *)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...