wolfram-mathematica – 在Mathematica中对Position2D的快速实现

我正在寻找一个快速实现以下，我称之为Position2D缺乏一个更好的术语：

Position2D[ matrix,sub_matrix ]

它找到矩阵内的子矩阵的位置，并返回匹配的左上和右下行/列。

例如：

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

应该返回：

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

它应该足够快，可以在具有100×100子矩阵的3000×2000矩阵上快速工作。为了简单起见，只考虑整数矩阵就足够了。

解决方法

算法

以下代码基于有效的定制位置函数来查找大整数列表中(可能重叠)整数序列的位置。主要思想是我们可以先尝试找到扁平化大矩阵中子矩阵的第一行的位置，然后对其进行滤波，提取完整的子矩阵并与子矩阵的子矩阵进行比较利益。大多数情况下，除非非常病态的病例(对于这个程序会产生大量潜在的职位候选人)，而且子矩阵的真实数量将小得多，这样做将是有效的，但是这种情况似乎不大可能一般而言，可以进一步改进这种简单的方案)。

对于大矩阵，当使用序列位置函数的编译版本时，所提出的解决方案将比@Szabolcs的解决方案快15-25倍，而顶层实现序列位置的时间要快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}]
Out[71]= {2,4}

用于整数的更快的位置查找功能

然而，快速seqPos可能是，它仍然是我的解决方案的主要瓶颈。这是一个编译到C的版本，这给我的代码提供了另外5倍的性能提升：

seqposC = 
  Compile[{{list,_Integer,{seq,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}]
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,5}};
what = {{2,4}};

这将计算上述的尺寸(使用子矩阵的反向尺寸更方便)：

In[78]:= 
dm=Dimensions[m]
dwr=Reverse@Dimensions[what]

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

在Flattened主矩阵中找到第一行的开始位置列表({2,3})。这些位置同时位于子矩阵左上角的“平”候选位置：

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

这将使用主矩阵的维度，以全矩阵重构子矩阵左上角的二维“候选”位置：

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

然后我们可以尝试使用Select来过滤掉它们，提取完整的子矩阵并与之进行比较，但是我们在这里遇到一个问题：

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,5}}. >>

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

除了错误信息，结果是正确的。错误消息本身是由于对于列表中的最后一个位置({1,3})，子矩阵的右下角将在主矩阵之外。我们当然可以使用安静来忽略错误消息，但这是一个坏的风格。所以我们首先将这些情况进行过滤，这就是Pick [Transpose [#]，Total [Clip [Reverse @ dm – # – dwr 2，{0，1}]]，2]& @对于。具体来说，请考虑

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

左上角的坐标应保持在矩阵和子矩阵的维度之间。上述子列表由左上角的x和y坐标组成。我添加了2，使所有有效结果严格正。我们必须在Transpose @ {Mod [posFlat，#，1]，IntegerPart [posFlat /#] 1}& @Last [dm](这是posInterm)中的这些位置上只选择协调者，在这两个子列表中有正数。我使用Total [Clip […，{0,1}]]将其重新选择为仅在这个第二个列表有2的位置(Clip将所有正整数转换为1，以及2个子列表中的总和数)。获得2的唯一方法是两个子列表中的数字都是正数)。

所以，我们有：

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

在2D位置列表被过滤之后，不存在结构上无效的位置，我们可以使用Select来提取完整​​的子矩阵并针对感兴趣的子矩阵进行测试：

In[93]:= 
finalPos = 
  Select[pos2D,First@#;;First@#+First@dwr-1]]==what&]
Out[93]= {{3,2}}

在这种情况下，这两个职位都是真实的。最后要做的是重建子矩阵的右下角的位置，并将它们添加到左上角的位置。这是由这一行完成的：

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

我可以使用Map，但是对于大量的位置，上面的代码将会更有效率。

示例和基准

原来的例子：

In[216]:= Position2D[{{0,6}},4}}]
Out[216]= {{{3,{{1,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},2560}},<<151>>,{{1900,42},{1999,141}}}}

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

(实际的子矩阵数小于我们尝试生成的数量，因为它们中的许多重叠并且“破坏”先前插入的子矩阵 – 这是因为子矩阵大小是矩阵大小的相当大的一部分我们的基准)。

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

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

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