Недавно вышла статья о реализации задачи по поиску максимального по площади квадрата внутри массива из единиц и нулей. Автор предлагал её решить через префиксные суммы.
Я какое-то время обдумывал её, а потом полез читать комментарии. Там показали способ решения с использованием динамического программирования (DP). Но меня всё еще не устраивало, что нам приходится обрабатывать полностью все поля (точнее итеративно N*M и при найденной единице "количество единиц * 4").
Сначала я придумал что-то такое:
Бежим построчно (на рисунке выше строки выделены разным цветом), считаем последовательности единиц, потом склеиваем построчно, в местах склеек будут образовываться квадраты. Но это усложняло поиск, пришлось бы хранить отдельно информацию об интервалах и никуда не убирало необходимость полного прохода N*M.
Тогда мои мысли навели меня на другое решение — мы будем искать сразу максимальный квадрат, но не перебирая все ячейки, а создавая списки из квадратов, вложенных друг в друга и если в каком‑то квадрате будет обнаружен ноль, то он удалит из списка все квадраты, которые содержат в себе эту ячейку. Вероятно хватило бы максимальной кучи, после прохода по массиву в вершине оставался бы ответ.
Посмотрите на картинку ниже, очевидно что квадрат 9×9 существовать не может, так как при чтении левого верхнего и правого нижнего углов мы получаем нули. Исследуемая зона сужается (на картинке ниже видно, что слева сверху и справа снизу мы установили ограничение)
Таким образом в работе мы сможем продолжить поиски только в зоне (2, 1), (9,8)
и (1,2), (8, 9)
Так как мы движемся слева направо и сверху вниз, и проверяем только левый верхний и правый нижний элементы, то у нас на (2, 1) единица, а вот на (9, 8) уже 0.
Дальше этих размышлений я не зашёл и конечно до практической реализации не добрался, так как мне показалось, что нужно сделать что-то попроще и я опять подумал о выделении зон строками, а чтобы не выискивать какой-то новый алгоритм склеивания строк я подумал об уже готовой реализации из комментария под статьёй с тем самым динамическим программированием. Огромный плюс DP в том, что мы используем ранее посчитанные данные. Но как нам уйти от необходимости полного обхода массива? Для DP узким местом всегда будет массив в большей степени заполненный единицами, а не нулями.
И тут я подумал, что для начала нужно отфильтровать массив очень грубо и так как полный перебор нас не устраивает, то можно просматривать таблицу через сито 2х2, то есть каждый четвёртый элемент. Что нам это даёт?
Не забываем что мы движемся слева направо сверху вниз с шагом 2. Если мы находим 0, то в этом блоке точно нет квадрата 2х2 и тем более 3х3 или 4х4. Как видно на рисунке ниже при проверке (2, 2) мы получаем ноль, поэтому ячейки (1, 1), (2, 1) и (1, 2) мы более не рассматриваем.
В каждой итерации мы фиксируем левый верхний угол и если вдруг нашли 1, то мы точно знаем, что этот блок нужно исследовать подробнее, но не сразу. Мы уменьшаем шаг до 1 и сдвигаемся вправо, пока не встретим 0 или край массива справа. Так как для работы реализации через динамическое программирование важно чтобы в трёх ячейках слева, сверху и слева-сверху была актуальная информация, то мы прогоняем подсчёт для сектора в 3 строки, с ограничением по X. Давайте поэтапно посмотрим на картинках
Так выглядит наш массив с обозначением тех мест, куда мы попадаем, если будем просто двигаться с шагом по 2 ячейки, а ниже то, как работает алгоритм
На первом шаге мы проверяем желтую ячейку фиксируя левую верхнюю координату возможной зоны, где потребуется расчет, обозначено зеленым. Заметьте, если желтую ячейку мы проверяем на 0, то о зеленой мы ничего не знаем.
На втором шаге мы находим 1 на желтом поле и опять зеленое поле - мы фиксируем координаты. Так как мы нашли 1, то будем двигаться вправо c шагом 1, а не 2.
На рисунке выше мы двигались по 1 вправо и попали на 0 в ячейке (7,2). Пунктиром я выделил зону, по которой мы проходим и подсчитываем maxS (сохраню название переменной из прошлой статьи), то есть максимальный квадрат из единиц. Обратите внимание, что считаем мы на одну строку ниже, это необходимо, чтобы этот алгоритм работал и maxS читался верно, только так мы можем передать информацию о найденных квадратах на строки ниже.
Условно все эти действия назовём тюнингом по отношению к оригинальному алгоритму. И это — тюнинг 2. Первым считаем забег по ячейкам с шагом 2.
Так как все строки, в которых мы производим расчёт имеют высоту 3, то в таких полосках массива, как может показаться, мы найдём только квадраты размером 3. Но помните я в самом начале рассказывал о полосках и их склеивании. Тут мы этот механизм получаем из коробки. Если блоки ниже будут иметь общую площадь, то расчеты с верхних ячеек перейдут вниз.
Тюнинг 3 — если мы уже нашли какое‑то значение maxS (напомню, это имя переменной, обозначающее размер стороны максимального квадрата в массиве в оригинальной статьи), то нет смысла обсчитывать зоны, которые по ширине меньше или равны maxS. Например мы нашли уже квадрат 5×5 и если нам попадаются зоны 3×3, 4×3, или 5×3 (вы же помните что 3 здесь это высота зон, которые мы обрабатываем), то обсчитывать их нет смысла, так как они не могут быть частью квадрата 6×6, который нам мог бы быть нужен.
На картинке выше уже обнаруженная зона A, мы находим желтую зону B, она шириной 3, что меньше уже найденной стороны maxS == 5. А вот для зоны C мы хоть и нашли «дурацкий» огрызок 2×2 в желтой зоне, но сама зона шириной 6 и ниже может скрываться квадрат 6×6, если мы его тут не посчитаем, то мы его пропустим.
Тюнинг 4 — для обсчёта прямоугольников мы пишем функцию:
int scanRectangle(int startX, int startY, int endX)
{
int maxS = 0;
for (int x = startX; x <= endX; x++)
{
square[x, startY] += Math.Min(
Math.Min(square[x, startY - 1], square[x - 1, startY]),
square[x - 1, startY - 1]);
maxS = Math.Max(maxS, square[x, startY]);
}
startY++;
if (startY < square.GetLength(0))
for (int x = startX; x <= endX; x++)
if (square[x, startY] == 1)
{
square[x, startY] += Math.Min(
Math.Min(square[x, startY - 1], square[x - 1, startY]),
square[x - 1, startY - 1]);
maxS = Math.Max(maxS, square[x, startY]);
}
return maxS;
}Зная высоту такого блока нам не требуется вторая координата Y, но так же мы точно знаем, что во второй строке все ячейки заполнены единицами:
Так происходит потому что наш алгоритм фиксирует первый ноль, отмечен желтым на (2, 2), потом на следующей итерации ловит 1 на (4, 2), в этом месте происходит фиксация координат (4, 1) как левый верхний угол для зоны обсчёта (на рисунке зона обозначена пунктиром). Далее мы меняем шаг на 1 и фиксируем все единицы в 5 и 6 столбцах, пока не ловим снова ноль в 7 столбце. Вот и получается что в такой зоне всегда гарантированно во второй строке будут только единицы. А то что на рисунке обозначено светло-коричневым - это те ячейки, значения которых мы будем читать и использовать для заполнения зоны алгоритмом DP. Заранее знать о том нули там или единицы нам не нужно, DP сам разберётся. После обработки массив будет выглядеть так:
И maxS будет равен 2.
Ниже код на C# реализации этого алгоритма:
const int scale = 1;
int[,] square = new int[100 * scale, 100 * scale];
public int Zara6502()
{
int maxS = 0;
int step = 2;
int currX = 0;
int Ysize = square.GetLength(0);
int Xsize = square.GetLength(1);
for (int currY = 1; currY < Ysize; currY += 2)
{
int tempX0 = 0;
step = 2;
for (currX = 1; currX < Xsize; currX += step)
{
if (step == 2) tempX0 = currX;
if (square[currX, currY] == 0)
{
if (step == 1 && (currX - tempX0) > maxS)
maxS = Math.Max(maxS, scanRectangle(tempX0, currY, currX - 1));
step = 2;
}
else step = 1;
}
if (step == 1 && (Xsize - tempX0) > maxS)
maxS = Math.Max(maxS, scanRectangle(tempX0, currY, Xsize - 1));
}
if (maxS == 0)
{
for (int currY = 0; currY < Ysize; currY += 2)
for (currX = 0; currX < Xsize; currX++)
if (square[currX, currY] == 1)
{
maxS = 1;
return maxS;
}
for (int currY = 1; currY < Ysize; currY += 2)
for (currX = 0; currX < Xsize; currX += 2)
if (square[currX, currY] == 1)
{
maxS = 1;
return maxS;
}
}
return maxS;
}
int scanRectangle(int startX, int startY, int endX)
{
int maxS = 0;
for (int x = startX; x <= endX; x++)
{
square[x, startY] += Math.Min(
Math.Min(square[x, startY - 1], square[x - 1, startY]),
square[x - 1, startY - 1]);
maxS = Math.Max(maxS, square[x, startY]);
}
startY++;
if (startY < square.GetLength(0))
for (int x = startX; x <= endX; x++)
if (square[x, startY] == 1)
{
square[x, startY] += Math.Min(
Math.Min(square[x, startY - 1], square[x - 1, startY]),
square[x - 1, startY - 1]);
maxS = Math.Max(maxS, square[x, startY]);
}
return maxS;
}Ну а теперь тестируем решения: динамическое программирование и авторский. Я выбрал такие шаблоны:
random 50% — заполнение массива единицами в случайные позиции, соотношение 0 и 1 всегда строго 50%
random squares with noises — генерируем квадраты размером от 2 до 10 в количестве (width * 0.1), для scale=1 это 10 квадратов на квадрат 100×100. Но перед размещением квадратов генерируется 30% шум из единиц
empty — пустой массив с одной единицей в правом нижнем углу
fullfill — полностью заполненный единицами массив
random squares no noises — такой же как п.2, но без шума
Ну и результаты:
random 50%:
===========
| Method | scale | Mean | Error | StdDev | Ratio | RatioSD | Allocated | Alloc Ratio |
|-------------- |------ |----------------:|--------------:|--------------:|------:|--------:|----------:|------------:|
| Zara6502 | 1 | 19.85 us | 0.210 us | 0.196 us | 1.00 | 0.01 | - | NA |
| petropavelMod | 1 | 94.99 us | 0.585 us | 0.518 us | 4.79 | 0.05 | - | NA |
| | | | | | | | | |
| Zara6502 | 10 | 7,690.34 us | 151.352 us | 155.427 us | 1.00 | 0.03 | 6 B | 1.00 |
| petropavelMod | 10 | 26,028.98 us | 132.892 us | 124.307 us | 3.39 | 0.07 | 4 B | 0.67 |
| | | | | | | | | |
| Zara6502 | 100 | 1,097,519.83 us | 10,810.220 us | 10,111.886 us | 1.00 | 0.01 | 400 B | 1.00 |
| petropavelMod | 100 | 3,357,881.64 us | 18,745.504 us | 16,617.401 us | 3.06 | 0.03 | 400 B | 1.00 |random squares with noises:
===========================
| Method | scale | Mean | Error | StdDev | Ratio | RatioSD | Allocated | Alloc Ratio |
|-------------- |------ |-----------------:|---------------:|---------------:|------:|--------:|----------:|------------:|
| Zara6502 | 1 | 9.584 us | 0.1051 us | 0.0983 us | 1.00 | 0.01 | - | NA |
| petropavelMod | 1 | 49.209 us | 0.1180 us | 0.0985 us | 5.13 | 0.05 | - | NA |
| | | | | | | | | |
| Zara6502 | 10 | 4,552.811 us | 27.2865 us | 24.1888 us | 1.00 | 0.01 | 3 B | 1.00 |
| petropavelMod | 10 | 19,623.032 us | 12.4785 us | 9.7424 us | 4.31 | 0.02 | 12 B | 4.00 |
| | | | | | | | | |
| Zara6502 | 100 | 778,358.150 us | 5,154.4114 us | 4,569.2514 us | 1.00 | 0.01 | 400 B | 1.00 |
| petropavelMod | 100 | 2,662,157.893 us | 21,424.1769 us | 20,040.1884 us | 3.42 | 0.03 | 400 B | 1.00 |
empty:
======
| Method | scale | Mean | Error | StdDev | Ratio | RatioSD | Allocated | Alloc Ratio |
|-------------- |------ |-----------------:|---------------:|---------------:|------:|--------:|----------:|------------:|
| Zara6502 | 1 | 5.831 us | 0.0177 us | 0.0165 us | 1.00 | 0.00 | - | NA |
| petropavelMod | 1 | 19.090 us | 0.0546 us | 0.0511 us | 3.27 | 0.01 | - | NA |
| | | | | | | | | |
| Zara6502 | 10 | 1,217.029 us | 12.5058 us | 11.0861 us | 1.00 | 0.01 | 1 B | 1.00 |
| petropavelMod | 10 | 15,574.495 us | 114.6756 us | 107.2676 us | 12.80 | 0.14 | 6 B | 6.00 |
| | | | | | | | | |
| Zara6502 | 100 | 492,762.913 us | 7,201.3221 us | 6,736.1212 us | 1.00 | 0.02 | 400 B | 1.00 |
| petropavelMod | 100 | 1,914,521.580 us | 29,803.5592 us | 27,878.2678 us | 3.89 | 0.07 | 400 B | 1.00 |
fullfill:
=========
| Method | scale | Mean | Error | StdDev | Ratio | RatioSD | Allocated | Alloc Ratio |
|-------------- |------ |----------------:|--------------:|--------------:|------:|--------:|----------:|------------:|
| Zara6502 | 1 | 11.37 us | 0.026 us | 0.023 us | 1.00 | 0.00 | - | NA |
| petropavelMod | 1 | 19.09 us | 0.036 us | 0.032 us | 1.68 | 0.00 | - | NA |
| | | | | | | | | |
| Zara6502 | 10 | 7,560.39 us | 35.708 us | 33.401 us | 1.00 | 0.01 | 3 B | 1.00 |
| petropavelMod | 10 | 15,161.26 us | 32.087 us | 28.444 us | 2.01 | 0.01 | 6 B | 2.00 |
| | | | | | | | | |
| Zara6502 | 100 | 1,152,131.31 us | 14,010.416 us | 13,105.352 us | 1.00 | 0.02 | 400 B | 1.00 |
| petropavelMod | 100 | 2,157,440.79 us | 3,971.820 us | 3,100.936 us | 1.87 | 0.02 | 112 B | 0.28 |
random squares no noises:
=========================
| Method | scale | Mean | Error | StdDev | Ratio | RatioSD | Allocated | Alloc Ratio |
|-------------- |------ |-----------------:|---------------:|---------------:|------:|--------:|----------:|------------:|
| Zara6502 | 1 | 6.608 us | 0.0354 us | 0.0331 us | 1.00 | 0.01 | - | NA |
| petropavelMod | 1 | 20.060 us | 0.0230 us | 0.0204 us | 3.04 | 0.01 | - | NA |
| | | | | | | | | |
| Zara6502 | 10 | 1,502.952 us | 11.9526 us | 10.5957 us | 1.00 | 0.01 | 1 B | 1.00 |
| petropavelMod | 10 | 15,548.515 us | 96.1894 us | 89.9756 us | 10.35 | 0.09 | 6 B | 6.00 |
| | | | | | | | | |
| Zara6502 | 100 | 604,071.771 us | 5,431.8570 us | 4,815.1998 us | 1.00 | 0.01 | 400 B | 1.00 |
| petropavelMod | 100 | 2,194,348.127 us | 13,435.8735 us | 12,567.9244 us | 3.63 | 0.03 | 112 B | 0.28 |
Спасибо что дочитали до конца. Готов отбивать летящие помидоры :-) Всем бобра!
P.S.: может кто подскажет какой массив для авторского алгоритма станет бутылочным горлышком? и почему для массива 1000х1000 авторский алгоритм иногда бустился до х10 по сравнению с DP?