wolfram-mathematica - 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编码的大矩阵中的位置，然后对其进行过滤，提取出完整的子矩阵，然后与的子矩阵进行比较利益。对于大多数情况，这将是非常有效的，但非常病理的情况除外(对于这种情况，此过程将生成大量潜在的职位候选者，而子矩阵的实际条目数会少得多。但是，这种情况似乎不太可能通常，然后可以对该简单方案进行进一步的改进)。

对于大型矩阵，当使用序列位置函数的编译版本时，建议的解决方案将比@Szabolcs的解决方案快约15-25倍，而对于序列位置的顶层实现-查找功能，建议的解决方案将比@Szabolcs的解决方案快约3-5倍。 。实际的加速取决于矩阵大小，对于较大的矩阵，则更大。代码和基准如下。

码

用于查找子列表(序列)位置的一般有效功能

这些辅助函数归功于Norbert Pozar，并取自this 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(这是默认设置)。对于一般 list ，可以提供例如 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} ed主矩阵中找到第一行的起始位置(此处为 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 coordiantes组成。我加了2，使所有有效结果严格为正。我们只需要在 Transpose@{Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm](即 posInterm)中的那些位置上选择coordiantes，上面的两个子列表都具有严格的正数。我使用 Total[Clip[...,{0,1}]]将其重铸为仅在第二个列表具有 2的位置上进行选择( Clip将所有正整数转换为 1，并且 Total将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}}}}

(子矩阵的实际数量小于我们尝试生成的数量，因为许多子矩阵重叠并“破坏”了先前插入的子矩阵-之所以如此，是因为子矩阵的大小是矩阵大小的相当大的一部分)我们的基准)。

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

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

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