在Mathematica中快速实现Position2D

Question

我正在寻找以下快速实现,我会称之为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子矩阵的矩阵.为简单起见,仅考虑整数矩阵就足够了.

Answer 1

算法

以下代码基于有效的自定义位置函数,以在大整数列表中查找(可能重叠)整数序列的位置.主要思想是我们可以首先尝试有效地找到子矩阵的第一行在Flatten-ed大矩阵中的位置,然后对这些位置进行滤波,提取完整的子矩阵并与感兴趣的子矩阵进行比较.这对于大多数情况都是有效的,除了非常病态的情况(这些程序将产生大量潜在的候选位置,而子矩阵的真实条目数将小得多.但这种情况似乎不太可能通常,然后,可以进一步改进这种简单的方案).

对于大型矩阵,当使用序列位置函数的编译版本时,所提出的解决方案将比@Szabolcs解决方案快约15-25倍,并且对于序列位置的顶级实现快3-5倍 - 查找函数.实际加速取决于矩阵大小,对于更大的矩阵更多.代码和基准如下.

码

用于查找子列表(序列)位置的通常有效的函数

这些辅助函数归功于Norbert Pozar并取自此 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版本,它为我的代码提供了另外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 = 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}

这将找到ed主矩阵中第一行({2,3}此处)的起始位置列表Flatten.这些位置同时是子矩阵左上角的"平坦"候选位置:

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,并将Total2个子列表中的数字相加.获得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}}}

请注意,我的索引约定与@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}}}}

(子矩阵的实际数量小于我们试图生成的数量,因为它们中的许多重叠并"破坏"先前插入的数据 - 这是因为子矩阵大小是矩阵大小的一个相当大的一部分.我们的基准).

为了比较,我们应该在其中一个解决方案(级别3)中反转xy索引,并对两个列表进行排序,因为位置可能以不同的顺序获得:

In[360]:= Sort@ps1===Sort[Reverse[ps2,{3}]]
Out[360]= True

我不排除可能进一步优化的可能性.

Answer 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获取潜在位置列表,然后过滤那些实际匹配的位置.它在打包的真实矩阵上可能更快.

Answer 3

根据列昂尼德的建议,这是我的解决方案.我知道它效率不高(当我计时时它比Leonid慢大约600倍)但是它非常短暂,可记忆,并且很好地说明了很少使用的功能,PartitionMap.它来自Developer包,因此需要Needs["Developer`"]先调用.

鉴于此,Position2D可以定义为:

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

这只给出了左上角的坐标.我觉得右下坐标实际上是多余的,因为子矩阵的尺寸是已知的,但如果需要,可以通过前面{#, Dimensions[k] + # - {1, 1}} & /@的定义添加到输出中.

Answer 4

怎么样的

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}}} *)