c++ - 查找二维矩阵是否是另一个二维矩阵的子集

Question

最近我参加了一个黑客马拉松，我开始了解一个问题，该问题试图在二维矩阵中找到网格形式的模式。模式可以是 U、H 和 T，将由 3* 表示3矩阵 假设如果我想呈现 H 和 U

+--+--+--+            +--+--+--+
|1 |0 |1 |            |1 |0 |1 |
+--+--+--+            +--+--+--+
|1 |1 |1 |  --> H     |1 |0 |1 |    -> U
+--+--+--+            +--+--+--+
|1 |0 |1 |            |1 |1 |1 |
+--+--+--+            +--+--+--+

现在我需要将其搜索到包含 0 和 1 的 10*10 矩阵。最接近且唯一的解决方案我可以获得 O(n^4) 的蛮力算法。在 MATLAB 和R 有非常微妙的方法可以做到这一点，但在 C、C++ 中却没有。我尝试了很多在谷歌和 SO 上搜索这个解决方案。但我能得到的最接近的是这个 SO POST其中讨论了关于实现 Rabin-Karp 字符串搜索算法 的问题。但是没有伪代码或任何对此进行解释的帖子。有人可以帮助或提供任何链接、pdf 或一些逻辑来简化这个吗？

编辑

作为 Eugene Sh.评论说如果 N 是大矩阵 (NxN) 和 k - 小矩阵 (kxk) 的大小，buteforce 算法应该采用 O((Nk)^2)。由于 k 是固定的，因此它减少到 O(N^2)。是的，绝对正确。 但是，如果 N 和 K 很大，是否有任何通用的方法？

Answer 1

好吧，这就是 2D Rabin-Karp方法。

对于下面的讨论，假设我们想在一个(n, >n) 矩阵。 (这个概念同样适用于矩形矩阵，但我用完了索引。)

我们的想法是，对于每个可能的子矩阵，我们计算一个散列。仅当该散列与我们要查找的矩阵的散列相匹配时，我们才会按元素进行比较。

为了提高效率，我们必须避免每次都重新计算子矩阵的整个散列。因为今晚我睡得很少，所以我唯一可以弄清楚如何轻松做到这一点的哈希函数是各个子矩阵中 1 的总和。我把它作为练习留给比我聪明的人找出更好的滚动哈希函数。

现在，如果我们刚刚检查了从 (i, j) 到 (i + m – 1, j + m – 1) 并且知道里面有 x 个 1，我们可以计算其中 1 的个数子矩阵向右 - 即从 (i, j + 1) 到 (i + m – 1, j + m) – 通过从 (i, j ) 到 (i + m – 1, j) 并从 (< i>i, j + m) 到 (i + m – 1, j + m).

如果我们碰到大矩阵的右边距，我们将窗口向下移动一位，然后回到左边距，然后再次向下移动一位，然后再次向右移动，依此类推。

请注意，这需要 O(m) 次操作，而不是每个候选人的 O(m²) 次。如果我们对每一对索引都这样做，我们得到 O(mn²) 的工作。因此，通过巧妙地通过大矩阵移动一个潜在子矩阵大小的窗口，我们可以将工作量减少 m 倍。也就是说，如果我们没有遇到太多哈希冲突。

这是一张图片:

随着我们将当前窗口右移一位，减去左边红色列 vector 中1的个数，加上右边绿色列 vector 中1的个数，得到1的个数在新窗口中。

我已经使用很棒的 Eigen 实现了这个想法的快速演示C++ 模板库。该示例还使用了 Boost 中的一些内容，但仅用于参数解析和输出格式化，因此如果您没有 Boost 但想尝试代码，则可以轻松摆脱它。索引摆弄有点困惑，但我将在此处不作进一步解释。上面的散文应该足以涵盖它。

#include <cassert>
#include <cstddef>
#include <cstdlib>
#include <iostream>
#include <random>
#include <type_traits>
#include <utility>

#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

#include <Eigen/Dense>

#define PROGRAM "submatrix"
#define SEED_CSTDLIB_RAND 1

using BitMatrix = Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>;
using Index1D = BitMatrix::Index;
using Index2D = std::pair<Index1D, Index1D>;

std::ostream&
operator<<(std::ostream& out, const Index2D& idx)
{
  out << "(" << idx.first << ", " << idx.second << ")";
  return out;
}

BitMatrix
get_random_bit_matrix(const Index1D rows, const Index1D cols)
{
  auto matrix = BitMatrix {rows, cols};
  matrix.setRandom();
  return matrix;
}

Index2D
findSubMatrix(const BitMatrix& haystack,
              const BitMatrix& needle,
              Index1D *const collisions_ptr = nullptr) noexcept
{
  static_assert(std::is_signed<Index1D>::value, "unsigned index type");
  const auto end = Index2D {haystack.rows(), haystack.cols()};
  const auto hr = haystack.rows();
  const auto hc = haystack.cols();
  const auto nr = needle.rows();
  const auto nc = needle.cols();
  if (nr > hr || nr > hc)
    return end;
  const auto target = needle.count();
  auto current = haystack.block(0, 0, nr - 1, nc).count();
  auto j = Index1D {0};
  for (auto i = Index1D {0}; i <= hr - nr; ++i)
    {
      if (j == 0)  // at left margin
        current += haystack.block(i + nr - 1, 0, 1, nc).count();
      else if (j == hc - nc)  // at right margin
        current += haystack.block(i + nr - 1, hc - nc, 1, nc).count();
      else
        assert(!"this should never happen");
      while (true)
        {
          if (i % 2 == 0)  // moving right
            {
              if (j > 0)
                current += haystack.block(i, j + nc - 1, nr, 1).count();
            }
          else  // moving left
            {
              if (j < hc - nc)
                current += haystack.block(i, j, nr, 1).count();
            }
          assert(haystack.block(i, j, nr, nc).count() == current);
          if (current == target)
            {
              // TODO: There must be a better way than using cwiseEqual().
              if (haystack.block(i, j, nr, nc).cwiseEqual(needle).all())
                return Index2D {i, j};
              else if (collisions_ptr)
                *collisions_ptr += 1;
            }
          if (i % 2 == 0)  // moving right
            {
              if (j < hc - nc)
                {
                  current -= haystack.block(i, j, nr, 1).count();
                  ++j;
                }
              else break;
            }
          else  // moving left
            {
              if (j > 0)
                {
                  current -= haystack.block(i, j + nc - 1, nr, 1).count();
                  --j;
                }
              else break;
            }
        }
      if (i % 2 == 0)  // at right margin
        current -= haystack.block(i, hc - nc, 1, nc).count();
      else  // at left margin
        current -= haystack.block(i, 0, 1, nc).count();
    }
  return end;
}

int
main(int argc, char * * argv)
{
  if (SEED_CSTDLIB_RAND)
    {
      std::random_device rnddev {};
      srand(rnddev());
    }
  if (argc != 5)
    {
      std::cerr << "usage: " << PROGRAM
                << " ROWS_HAYSTACK COLUMNS_HAYSTACK"
                << " ROWS_NEEDLE COLUMNS_NEEDLE"
                << std::endl;
      return EXIT_FAILURE;
    }
  auto hr = boost::lexical_cast<Index1D>(argv[1]);
  auto hc = boost::lexical_cast<Index1D>(argv[2]);
  auto nr = boost::lexical_cast<Index1D>(argv[3]);
  auto nc = boost::lexical_cast<Index1D>(argv[4]);
  const auto haystack = get_random_bit_matrix(hr, hc);
  const auto needle = get_random_bit_matrix(nr, nc);
  auto collisions = Index1D {};
  const auto idx = findSubMatrix(haystack, needle, &collisions);
  const auto end = Index2D {haystack.rows(), haystack.cols()};
  std::cout << "This is the haystack:\n\n" << haystack << "\n\n";
  std::cout << "This is the needle:\n\n" << needle << "\n\n";
  if (idx != end)
    std::cout << "Found as sub-matrix at " << idx << ".\n";
  else
    std::cout << "Not found as sub-matrix.\n";
  std::cout << boost::format("There were %d (%.2f %%) hash collisions.\n")
    % collisions
    % (100.0 * collisions / ((hr - nr) * (hc - nc)));
  return (idx != end) ? EXIT_SUCCESS : EXIT_FAILURE;
}

在编译和运行时，请将以上内容视为伪代码。我几乎没有尝试优化它。这只是我自己的一个概念验证。