Найдите, является ли матрица 2d подмножеством другой 2d-матрицы

Ответ 1

Хорошо, вот тогда 2D Рабин-Карп.

В следующем обсуждении предположим, что мы хотим найти подматрицу (m, m) внутри a ( n, n). (Эта концепция работает и для прямоугольных матриц, но у меня кончились индексы.)

Идея заключается в том, что для каждой возможной подматрицы мы вычисляем хэш. Только если этот хеш совпадает с хешем матрицы, которую мы хотим найти, мы сравним элемент-мудрый.

Чтобы сделать это эффективным, мы должны избегать повторного вычисления всего хэша субматрицы каждый раз. Поскольку сегодня я мало спал, единственная функция хэша, для которой я мог бы легко понять, как это сделать, - это сумма 1 с в соответствующей подматрице. Я оставляю это как упражнение кому-то умнее, чем мне, чтобы понять лучшую функцию хэширования.

Теперь, если мы только что проверили подматрицу от (i, j) до ( i + m - 1, j + m - 1) и знаем, что внутри него находится x 1s, мы можем вычислить число 1s в подматрица одна вправо, т.е. из ( i, j + 1) в ( i + m - 1, j + m) - путем вычитания числа 1s в под-векторе из ( i, j) до ( i + m - 1, j) и добавив число 1s в подвектор из ( i, j + m) до ( i + m - 1, j + m).

Если мы ударим по правому краю большой матрицы, мы сдвинем окно вниз на один, а затем обратно на левое поле, а затем снова вниз на один, а затем снова вправо и т.д.

Обратите внимание, что для этого требуется выполнение операций O (m), а не O ( m ²) для каждого кандидата. Если мы сделаем это для каждой пары индексов, получим O ( m n ²). Таким образом, умновав окно размера потенциальной подматрицы через большую матрицу, мы можем уменьшить объем работы на множитель m. То есть, если мы не получаем слишком много хэш-коллизий.

Вот изображение:

Когда мы сдвигаем текущее окно вправо, мы вычитаем число 1s в красном столбце с левой стороны и добавим число 1s в зеленый вектор столбца с правой стороны, чтобы получить число 1s в новом окне.

Я реализовал быструю демонстрацию этой идеи, используя большую библиотеку шаблонов Eigen. В этом примере также используются некоторые вещи из 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;
}

Пока он компилируется и запускается, рассмотрите вышеприведенное как псевдокод. Я почти не пытался его оптимизировать. Это было всего лишь доказательством концепции для меня.

Ответ 2

Я собираюсь представить алгоритм, который занимает O(n*n) время в худшем случае, когда k = O(sqrt(n)) и O(n*n + n*k*k) вообще. Это расширение Ахо-Корасика до 2D. Напомним, что Aho-Corasick находит все вхождения набора шаблонов в целевой строке T, и он делает это во времени, линейном по длине шаблона, длине T и количеству вхождений.

Введем некоторую терминологию. Стог сена - это большая матрица, в которой мы ищем, и игла является матрицей шаблонов. Стол сена является матрицей nxn, а игла представляет собой матрицу kxk. Набор узоров, которые мы собираемся использовать в Ахо-Корасике, представляет собой набор рядов иглы. Этот набор содержит не более k строк и будет иметь меньше, если есть повторяющиеся строки.

Мы собираемся построить автомат Aho-Corasick (который является дополнением Trie с ссылками на отказ), а затем запустить алгоритм поиска на каждой строке стога сена. Поэтому мы берем каждый ряд иглы и ищем ее в каждом ряду стога сена. Для этого мы можем использовать алгоритм линейного времени 1D, но это все равно будет неэффективным. Преимущество Aho-Corasick заключается в том, что он ищет все шаблоны сразу.

Во время поиска мы собираемся заполнить матрицу A, которую мы будем использовать позже. Когда мы ищем в первом ряду стога сена, первая строка A заполняется вхождениями строк иглы в первом ряду стога сена. Таким образом, мы получим первую строку A, которая выглядит, например, как 2 - 0 - - 1. Это означает, что номер строки 0 иглы появляется в позиции 2 в первом ряду стога сена; номер строки 1 появляется в позиции 5; номер строки 2 появляется в позиции 0. Записи - - это позиции, которые не совпадали. Продолжайте делать это для каждой строки.

Предположим теперь, что в игле нет повторяющихся строк. Присвойте 0 первому ряду иглы, 1 ко второму и т.д. Теперь мы будем искать шаблон [0 1 2 ... k-1] в каждом столбце матрицы A с использованием линейного алгоритма поиска 1D (например, KMP). Напомним, что каждая строка A хранит позиции, в которых появляются строки иглы. Поэтому, если столбец содержит шаблон [0 1 2 ... k-1], это означает, что номер строки 0 иглы появляется в некотором ряду стога сена, номер строки 1 иглы находится чуть ниже нее и так далее. Это именно то, что мы хотим. Если есть повторяющиеся строки, просто назначьте уникальный номер каждой уникальной строке.

Поиск в столбце принимает O(n) с использованием линейного алгоритма времени. Поэтому для поиска всех столбцов требуется O(n*n). Мы заполняем матрицу во время поиска, ищем каждую строку стога сена (есть строки n), а поиск в строке занимает O(n+k*k). Итак, O(n(n+k*k)) в целом.

Итак, идея заключалась в том, чтобы найти эту матрицу, а затем свести проблему к сопоставлению 1D-шаблонов. Ахо-Корасик просто эффективен, я не знаю, есть ли еще один эффективный способ найти матрицу.

EDIT: добавлена ​​реализация.

Вот моя реализация на С++. Максимальное значение n установлено равным 100, но вы можете его изменить.

Программа начинается с чтения двух целых чисел n k (размеры матриц). Затем он читает строки n, каждая из которых содержит строку из 0 и 1 длины n. Затем он читает строки k, каждая из которых содержит строку из 0 и 1 длины k. Вывод - верхняя левая координата всех совпадений. Например, для следующего ввода.

12 2
101110111011
111010111011
110110111011
101110111010
101110111010
101110111010
101110111010
111010111011
111010111011
111010111011
111010111011
111010111011
11
10

Программа выведет:

match at (2,0)
match at (1,1)
match at (0,2)
match at (6,2)
match at (2,10)

#include <cstdio>
#include <cstring>
#include <string>
#include <queue>
#include <iostream>

using namespace std;

const int N = 100;
const int M = N;
int n, m;

string haystack[N], needle[M];
int A[N][N]; /* filled by successive calls to match */
int p[N]; /* pattern to search for in columns of A */

struct Node
{   Node *a[2]; /* alphabet is binary */
    Node *suff; /* pointer to node whose prefix = longest proper suffix of this node */
    int flag;

    Node()
    {   a[0] = a[1] = 0;
        suff = 0;
        flag = -1;
    }
};

void insert(Node *x, string s)
{   static int id = 0;
    static int p_size = 0;

    for(int i = 0; i < s.size(); i++)
    {   char c = s[i];
        if(x->a[c - '0'] == 0)
            x->a[c - '0'] = new Node;
        x = x->a[c - '0'];
    }
    if(x->flag == -1)
        x->flag = id++;

    /* update pattern */
    p[p_size++] = x->flag;
}

Node *longest_suffix(Node *x, int c)
{   while(x->a[c] == 0)
        x = x->suff;
    return x->a[c];
}

Node *mk_automaton(void)
{   Node *trie = new Node;
    for(int i = 0; i < m; i++)
    {   insert(trie, needle[i]);
    }

    queue<Node*> q;

    /* level 1 */
    for(int i = 0; i < 2; i++)
    {   if(trie->a[i])
        {   trie->a[i]->suff = trie;
            q.push(trie->a[i]);
        }
        else trie->a[i] = trie;
    }

    /* level > 1 */
    while(q.empty() == false)
    {   Node *x = q.front(); q.pop();
        for(int i = 0; i < 2; i++)
        {   if(x->a[i] == 0) continue;
            x->a[i]->suff = longest_suffix(x->suff, i);
            q.push(x->a[i]);
        }
    }

    return trie;
}

/* search for patterns in haystack[j] */
void match(Node *x, int j)
{   for(int i = 0; i < n; i++)
    {   x = longest_suffix(x, haystack[j][i] - '0');
        if(x->flag != -1)
        {   A[j][i-m+1] = x->flag;
        }
    }
}

int match2d(Node *x)
{   int matches = 0;
    static int z[M+N];
    static int z_str[M+N+1];

    /* init */
    memset(A, -1, sizeof(A));

    /* fill the A matrix */
    for(int i = 0; i < n; i++)
    {   match(x, i);
    }

    /* build string for z algorithm */
    z_str[n+m] = -2; /* acts like `\0` for strings */
    for(int i = 0; i < m; i++)
    {   z_str[i] = p[i];
    }

    for(int i = 0; i < n; i++)
    {   /* search for pattern in column i */
        for(int j = 0; j < n; j++)
        {   z_str[j + m] = A[j][i];
        }

        /* run z algorithm */
        int l, r;
        l = r = 0;
        z[0] = n + m;
        for(int j = 1; j < n + m; j++)
        {   if(j > r)
            {   l = r = j;
                while(z_str[r] == z_str[r - l]) r++;
                z[j] = r - l;
                r--;
            }
            else
            {   if(z[j - l] < r - j + 1)
                {   z[j] = z[j - l];
                }
                else
                {   l = j;
                    while(z_str[r] == z_str[r - l]) r++;
                    z[j] = r - l;
                    r--;
                }
            }
        }

        /* locate matches */
        for(int j = m; j < n + m; j++)
        {   if(z[j] >= m)
            {   printf("match at (%d,%d)\n", j - m, i);
                matches++;
            }
        }
    }

    return matches;
}

int main(void)
{   cin >> n >> m;
    for(int i = 0; i < n; i++)
    {   cin >> haystack[i];
    }
    for(int i = 0; i < m; i++)
    {   cin >> needle[i];
    }

    Node *trie = mk_automaton();
    match2d(trie);

    return 0;
}

Ответ 3

Начнем с решения O(N * N * K). Я буду использовать следующие обозначения: A является матрицей шаблонов, B является большой матрицей (той, которую мы будем искать вхождения шаблона в).

• Мы можем исправить верхнюю строку матрицы B (т.е. будем искать все вхождения, которые начинаются в позиции (эта строка, любой столбец). Позвоните в эту строку a topRow. Теперь мы можем взять срез этой матрицы, содержащий строки [topRow; topRow + K) и все столбцы.

• Создайте новую матрицу в результате конкатенации A + column + the slice, где a column - столбец с элементами K, которые не присутствуют в A или B (если A и B состоят из 0 и 1, мы можем использовать, например, -1). Теперь мы можем рассматривать столбцы этой новой матрицы как буквы и запускать алгоритм Кнута-Морриса-Пратта. Для сравнения двух букв требуется O(K) время, поэтому временная сложность этого шага O(N * K).

Существует несколько способов исправления верхней строки O(N), поэтому общая временная сложность O(N * N * K). Это уже лучше, чем решение грубой силы, но мы еще не закончили. Теоретическая нижняя граница O(N * N) (я предполагаю, что N >= K), и я хочу ее достичь.

Посмотрим, что можно улучшить здесь. Если бы мы могли сравнить два столбца матрицы в O(1) времени вместо O(K), мы бы достигли желаемой временной сложности. Пусть объединяются все столбцы как A, так и B, вставляя некоторый разделитель после каждого столбца. Теперь у нас есть строка, и нам нужно сравнить ее подстроки (потому что столбцы и их части теперь подстроки). Пусть построено дерево суффикса в линейном времени (с использованием алгоритма Укконена). Теперь сравнение двух подстрок - это поиск высоты наименьшего общего предка (LCA) двух узлов в этом дереве. Существует алгоритм который позволяет нам делать это с использованием линейного времени предварительной обработки и O(1) времени на запрос LCA. Это означает, что мы можем сравнивать две подстроки (или столбцы) в постоянное время! Таким образом, общая временная сложность O(N * N). Существует еще один способ достижения этой временной сложности: мы можем построить массив суффикса в линейном времени и ответить на самые длинные общие запросы префикса в постоянное время (с линейной предварительной обработкой). Тем не менее, обе эти решения O(N * N) выглядят довольно сложно, и они будут иметь большую константу.

P.S Если у нас есть полиномиальная хеш-функция, которую мы можем полностью доверять (или у нас все хорошо с несколькими ложными срабатываниями), мы можем получить гораздо более простое решение O(N * N) с использованием двухмерных полиномиальных хэшей.