Алгоритм «танцующих ссылок» на Julia: реализация и влияние типизации на производительность

В этой статье я (опять) хочу рассмотреть алгоритм поиска решения задачи полного покрытия, теперь уже с нормальной реализацией через структуру "танцующих ссылок". Заодно на этом примере хочу показать, где и зачем указание типов в Julia критично для производительности, а где оно не обязательно.

Задача точного покрытия множества $$ некоторым набором $$ его подмножеств состоит в выборе попарно непересекающихся подмножеств $$, в сумме составляющих $$. Например, если множество $$ — это геометрическая фигура, а $$ — это кусочки "паркета", то задача точного покрытия спрашивает, можно ли замостить фигуру таким паркетом.

Как пример можно рассмотреть складывание геометрических фигур из фигурок пентамино

Некоторые решения для прямоугольников (из Википедии):

Для конечного множества $$ и конечного числа подмножеств $$ задача может быть представлена как матрица логических значений (0/1), где в $$-й строке единицы соответствуют элементам множества $$, содержащимся в $$. Задача тогда состоит в выборе таких строк, чтобы в получившейся матрице в каждом столбце была ровно одна единица.

Пример матрицы:

Строки 1, 4 и 5 составляют точное покрытие.

Задача является NP-полной, поэтому единственным гарантированным способом решения является поиск с возвратом. Особенности задачи позволяют построить специальную структуру данных для такого поиска, популяризованную Дональдом Кнутом, — систему "танцующих ссылок".

"Алгоритм X"

Алгоритм поиска точного покрытия, описанный Кнутом в работе "Танцующие ссылки", представляет собой простой поиск с возвратом. Псевдокод этой процедуры выглядит следующим образом:

Вход:     `matr` - матрица из непокрытых столбцов     `solution` - массив, содержащий строки решения     `k` - указатель на текущий индекс массива  Процедура `Algorithm_x`     Если `matr` не содержит непокрытых столбцов:         Вернуть (или вывести) `solution`.     Иначе выбрать столбец `c`.     Покрыть столбец `c`.     Для всех строк `r`, пересекающихся со стобцом `c`:         Записать `solution[k] = r`.         Для всех столбцов `j`, пересекающихся с `r`:             Покрыть столбец `j`.         Вызвать `Algorithm_x(matr, solution, k+1)`.         Если алгоритм нашёл решение, вернуть его.         (бэктрекинг)         Для всех столбцов `j`, пересекающихся с `r`:             Восстановить столбец `j`.     Восстановить столбец `c`.     Возврат из процедуры.

"Покрытие" столбца c означает вычёркивание из матрицы его и всех строк, с которыми он пересекается (поскольку если столбец покрыт какой-либо строкой r, то другими строками его покрыть уже нельзя). Восстановление означает вставку удалённого столбца и строк обратно.

Танцующие ссылки

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

Идея танцующих ссылок основана на том, что удаление узла node из двусвязного списка делается очень просто: если элемент имеет ссылку prev на предыдущий и next на следующий элемент, то удаление — это

node.next.prev = node.prev node.prev.next = node.next

При этом сам узел списка сохраняет информацию о том, какие узлы были его соседями. Благодаря этой информации узел можно восстановить:

node.next.prev = node node.prev.next = node

Если из списка удалено несколько элементов, то вставка их в порядке, обратном порядку удаления, восстанавливает исходное состояние списка.

Для булевой матрицы можно представлять единичные элементы узлами списка, связанного по двум направлениям. Для приведённой выше матрицы список будет выглядеть так (из статьи Кнута):

h — входной узел для структуры, в узлах A — G хранятся ссылки на соседей и число строк, пересекающих столбец, в остальных узлах хранятся только ссылки на ненулевые элементы в той же строке и в том же столбце. Также внутренние узлы хранят ссылки на заголовки столбцов, на рисунке не показанные.

Процедура покрытия столбца в этой структуре требует следующих действий:

• Отцепить заголовок столбца от правого и левого соседей.
• Проходя по ссылкам вниз, для каждой пересекаемой строки $r$ отцепить от соседей сверху и снизу все ссылки в той же строке, кроме собственно $r$.

Восстановление выполняется в обратном порядке — сначала восстановить для каждой строки (проходя вверх по ссылкам от заголовка столбца) соседей по вертикали, затем восстановить заголовок столбца в ряду. Чтобы структура корректно восстанавливалась, покрытие и восстановление столбцов должны вызываться в строго обратном порядке.

Рассмотрим реализацию этой структуры и алгоритма поиска точного покрытия на языке Julia.

Узел списка содержит ссылки на узлы слева, справа, сверху и снизу и на головной узел столбца, к которому он относится:

mutable struct LinkNode     left     right     above     below     col     # конструктор с привязкой к столбцу     function LinkNode(col)         self = new()         self.above = self.below = self.right = self.left = self         self.col = col         return self     end     # конструктор с пустым столбцом (для входного узла)     function LinkNode()         self = new()         self.above = self.below = self.right = self.left = self         return self     end end

Головной узел столбца содержит ссылки на соседние элементы, размер и идентификатор.

mutable struct LinkColumn     links     size     id     function LinkColumn(id)         self = new()         links = LinkNode(self)         self.links = links         self.size = 0         self.id = id         return self     end end

Сама матрица представляется входным узлом, от которого можно перейти к списку матриц.

struct LinkMatrix     root_node     #=     конструктор создает список столбцов из итерируемого объекта,     содержащего идентификаторы столбцов     =#     function LinkMatrix(ids)         root = LinkNode()         self = new(root)         prev = root         for id in ids             col = LinkColumn(id)             collinks = col.links             collinks.left = prev             collinks.right = root             prev.right = collinks             root.left = collinks             prev = collinks         end         return self     end end

Алгоритм поиска точного покрытия записывается в соответствии с псевдокодом:

function algorithm_x!(matr, solution=FullRow[])     if iscovered(matr)         return solution     end     col = choose_col(matr)     cover!(col)     for node in col         push!(solution, FullRow(node))         for j in RestRow(node)             cover!(j.col)         end         if !isnothing(algorithm_x!(matr, solution))             return solution         end         node = pop!(solution).node         # важен порядок восстановления         for j in Iterators.reverse(RestRow(node))             uncover!(j.col)         end     end     uncover!(col)     return end

Вспомогательные структуры данных и процедуры:

# проверка, покрыта ли матрица полностью iscovered(matr) = matr.root === matr.root.right  # итератор по столбцам матрицы function Base.iterate(matr::LinkMatrix, next = matr.root.right)     next === matr.root && return     return (next.col, next.right) end  #= структура данных, представляющая все узлы в строке, кроме заданного =# struct RestRow     node end  # итерация по всем узлам в строке, кроме заданного function Base.iterate(row::RestRow, next = row.node.right)     next === row.node && return     return (next, next.right) end  # итерация справа налево function Base.iterate(row::Iterators.Reverse{<:RestRow}, next = row.itr.node.left)     next === row.itr.node && return     return (next, next.left) end  # строка, содержащая заданный узел struct FullRow     node end  # итерация по строке, начиная с заданного узла function Base.iterate(row::FullRow)     node = row.node     return (node, node.right) end  function Base.iterate(row::FullRow, next)     next === row.node && return     return (next, next.right) end

Выбор столбца, который пытаемся накрыть следующим, делается на основе "S-эвристики" — выбирается столбец, для которого меньше всего возможных накрывающих его строк:

function choose_col(matr)     bestcol = matr.root.right.col     s = bestcol.size     for col in matr         l = col.size         l < 2 && return col # если возможен только 1 вариант или ни одного, выбираем сразу         if l < s             s, bestcol = l, col         end     end     return bestcol end

Процедура cover!(col) "отцепляет" столбец от правого и левого соседей (что важно — оставляя связи по вертикали без изменений, чтобы можно было потом восстановить), затем проходится по каждой строке и отцепляет узлы от соседей сверху и снизу (кроме тех ссылок, которые по цепочке ведут к столбцу col). Когда узел отцепляется от соседей по вертикали, размер столбца, где он находился, уменьшаем на единицу.

function cover!(col)     detach_rl!(col.links)     for node in col         for elt in RestRow(node)             detach_ab!(elt)         end     end end  # отсоединяет от соседей слева и справа function detach_rl!(node)     node.right.left = node.left     node.left.right = node.right     return node end  # отсоединяет от соседей сверху и снизу function detach_ab!(node)     node.above.below = node.below     node.below.above = node.above     node.col.size -= 1     return node end

Процедура восстановления, uncover!(col) возвращает всё назад:

function uncover!(col)     # важен порядок прохода по столбцу     for node in Iterators.reverse(col)         # порядок прохода по строке не важен         for elt in RestRow(node)             restore_ab!(elt)         end     end     restore_rl!(col.links) end  function restore_rl!(node)     node.right.left = node     node.left.right = node     return node end  function restore_ab!(node)     node.above.below = node     node.below.above = node     node.col.size += 1     return node end

Всё, что осталось для удобной жизни — это добавить процедуру добавления строки, которая накрывает столбцы с заданными идентификаторами:

#= matr - матрица, куда вставляется строка col_ids - коллекция с идентификаторами =# function insert_row!(matr, col_ids)     isempty(col_ids) && return     num_inserted = 0     first_in_row = nothing     prevnode = nothing     # идём по столбцам, если id подходит - создаём новый узел, прицепляем к строке     for col in matr         if col.id in col_ids             newnode = LinkNode(col)             if isnothing(first_in_row)                 first_in_row = newnode             end             first_in_row.left = newnode             newnode.right = first_in_row             if !isnothing(prevnode)                 prevnode.right = newnode                 newnode.left = prevnode             end             prevnode = newnode             num_inserted += 1             # если все идентификаторы найдены - подцепляем строку к матрице             if num_inserted == length(col_ids)                 for node in FullRow(first_in_row)                     col = node.col                     collinks = col.links                     lastnode = collinks.above                     col.size += 1                     node.above = lastnode                     node.below = collinks                     lastnode.below = node                     collinks.above = node                 end                 return first_in_row             end         end     end     #=     сюда попадаем, если не нашли все идентификаторы     в этом случае строка содержит идентификатор, которого нет в столбцах     игнорируем эту строку и выходим     =#     return end

Скорость работы

Проверим алгоритм на решении судоку и замерим скорость.

Судоку преобразуется в следующую задачу полного покрытия:

• (:row, row, num) — в строке row есть число num
• (:col, col, num) — в столбце col есть число num
• (:block, col, num) — в столбце block есть число num
• (самое хитрое условие) (:fill, row, col) — на пересечении столбца col и строки row есть число

Наличие уже заполненных клеток приводит к тому, что часть столбцов уже "накрыта", и их в матрицу задачи вносить не нужно (или можно внести и сразу же "закрыть").

function sudoku2xcover(fill)     ncells = size(fill)     side = ncells[1]     side == ncells[2] || throw(ArgumentError("Invalid matrix size $ncells"))     blockside = isqrt(side)     blockside^2 == side || throw(ArgumentError("Side $side is not a full square"))     row_constr = trues(ncells)     col_constr = trues(ncells)     blk_constr = trues(ncells)     for col in 1:side, row in 1:side         num = fill[row,col]         if num in 1:side             block = col - (col - 1) % blockside + (row - 1) ÷ blockside             row_constr[row, num] = false             col_constr[col, num] = false             blk_constr[block, num] = false         end     end     constr = collect((:row, row, num) for row in 1:side, num in 1:side if row_constr[row,num])     append!(constr, (:col, col, num) for col in 1:side, num in 1:side if col_constr[col,num])     append!(constr, (:block, block, num) for block in 1:side, num in 1:side if blk_constr[block,num])     append!(constr, (:fill, row, col) for row in 1:side, col in 1:side if !(fill[row,col] in 1:side))     matr = LinkMatrix(constr)     for row in 1:side, col in 1:side, num in 1:side         block = col - (col - 1) % blockside + (row - 1) ÷ blockside         if row_constr[row,num] && col_constr[col,num] && blk_constr[block,num] && !(fill[row,col] in 1:side)             col_ids = (:row, row, num), (:col, col, num), (:block, block, num), (:fill, row, col)             insert_row!(matr, col_ids)         end     end     return matr end  function sol2matr(soln, fieldsize)     ans = zeros(Int, fieldsize)     for node in soln         indval = row2indval(node)         ans[indval.ind] = indval.val     end     return ans end  function row2indval(dlrow)     row, col, n = 0,0,0     for j in dlrow         constr = j.col.id         if constr[1] === :fill             _, row, col = constr         else             _, _, n = constr         end         all(>(0), (row, col, n)) && break     end     return (ind = CartesianIndex(row, col), val = n) end  function sudoku(fill)     cover = sudoku2xcover(fill)     soln = sol2matr(algorithm_x!(cover), size(fill))     soln .+= fill     return soln end

Посмотрим на время решения "самого сложного судоку" под названием "Золотой самородок"

const testpuzzle = [         0 0 0 0 0 0 0 3 9;         0 0 0 0 1 0 0 0 5;         0 0 3 0 0 5 8 0 0;         0 0 8 0 0 9 0 0 6;         0 7 0 0 2 0 0 0 0;         1 0 0 4 0 0 0 0 0;         0 0 9 0 0 8 0 5 0;         0 2 0 0 0 0 6 0 0;         4 0 0 7 0 0 0 0 0     ] julia> using BenchmarkTools  julia> @btime sudoku(testpuzzle);   33.012 ms (115869 allocations: 3.40 MiB)

(ответ не показан, чтобы не спойлить)
Решение требует порядка 0.1 секунды. Быстрее, чем вручную, но как-то не впечатляет. Попробуем "переписать всё в сишном стиле", чтобы ублажить компилятор и ускорить выполнение.

Ускоряем работу, проставив типы

Главная проблема описанной реализации с точки зрения производительности — структуры данных имеют поля с абстрактным типом (когда тип данных поля не указан, он ставится как Any). Поэтому при обращении к полям постоянно требуется боксинг и анбоксинг данных.

Конкретизируем типы полей данных, по возможности. И сразу сталкиваемся с весёлой проблемой. Определение

# так не выйдет mutable struct LinkNode     left::LinkNode     right::LinkNode     above::LinkNode     below::LinkNode     col::LinkColumn end  mutable struct LinkColumn     links::LinkNode     size::Int     id::Any end

не проходит, поскольку LinkColumn — неизвестный тип на этапе определения типа LinkNode. По той же самой причине не проходит определение сперва LinkColumn, а потом LinkNode.

Решением в данном случае будут параметрические типы. Поставим тип поля col в структуре LinkNode как параметр, а поле links в LinkColumn будет иметь тип LinkNode{LinkColumn}. Чтобы избавиться от абстрактного типа для id, параметризуем и LinkColumn типом этого поля.

mutable struct LinkNode{C}     left::LinkNode{C}     right::LinkNode{C}     above::LinkNode{C}     below::LinkNode{C}     col::C     function LinkNode(col::C) where {C}         self = new{C}()         self.above = self.below = self.right = self.left = self         self.col = col         return self     end     function LinkNode{C}() where {C}         self = new{C}()         self.above = self.below = self.right = self.left = self         return self     end end  mutable struct LinkColumn{T}     links::LinkNode{LinkColumn{T}}     size::Int     id::T     function LinkColumn{T}(id) where {T}         self = new{T}()         links = LinkNode(self)         self.links = links         self.size = 0         self.id = id         return self     end end  LinkColumn(id::T) where {T} = LinkColumn{T}(id)

С LinkMatrix проделываем ту же процедуру.

struct LinkMatrix{T}     root::LinkNode{LinkColumn{T}}     function LinkMatrix{T}(ids) where {T}         root = LinkNode{LinkColumn{T}}()         self = new{T}(root)         prev = root         for id in ids             col = LinkColumn{T}(id)             collinks = col.links             collinks.left = prev             collinks.right = root             prev.right = collinks             root.left = collinks             prev = collinks         end         return self     end end  LinkMatrix(ids) = LinkMatrix{eltype(ids)}(ids)

Нужно не забыть указать типы полей во вспомогательных структурах:

struct RestRow{L<:LinkNode}     node::L end  struct FullRow{L<:LinkNode}     node::L end

Проверяем:

julia> @btime sudoku(testpuzzle);   777.426 μs (6803 allocations: 209.91 KiB)

Время выполнения улучшилось с 33 мс до 0.8 мс при том, что типы переменных не указаны нигде, кроме полей структур, конструкторов и перегруженных методов для iterate, где аннотации типов необходимы для правильной диспетчеризации.

Однако стоит дописать аннотацию типа ещё в одном месте — для аргументов процедуры algorithm_x!, где в качестве аргумента по умолчанию указан контейнер с абстрактным типом элементов. Благодаря множественной диспетчеризации можно напрямую перегрузить метод, вызываемый с одним аргументом.

function algorithm_x!(matr::LinkMatrix{T}) where {T}     RowType = FullRow{LinkNode{LinkColumn{T}}}     solution=RowType[]     return algorithm_x!(matr, solution) end

Проверяем:

julia> @btime sudoku(testpuzzle);   547.812 μs (3603 allocations: 129.16 KiB)

Итого получили ускорение больше чем в полсотни раз. При этом каких-то нечеловеческих усилий для расстановки типов прилагать вовсе не приходится, достаточно это сделать только в ключевых местах, в первую очередь — в полях структур (где, в большинстве случаев, о предполагаемых типах полей так и так нужно задумываться). Сильно в этом помогают параметрические типы — если заранее неизвестно, какого типа может быть поле структуры, чаще всего его стоит не оставлять без типа (т.е. с типом Any), а делать его тип параметром:

struct SomeStruct{T}     data::T end

Так поле data может быть чем угодно, но не придётся жертвовать скоростью выполнения.

Заключение

Дизайн языка Julia позволяет использовать его не только как "язык для математиков". Разумная система типов и активное использование аннотаций типов при JIT компиляции позволяют эффективно реализовывать любые классические Computer Science алгоритмы. Особо заморачиваться с типизацией для получения быстрого кода при этом не нужно, что показано примером реализации алгоритма DLX Кнута. Но выставление аннотаций типов в стратегически выверенных точках вполне может дать прирост скорости на 1-2 порядка, поэтому тем, кто начинает программировать на Julia (и тем, кто пробует язык после других языков с динамической типизацией), стоит потратить некоторое время на понимание концепций абстрактных и конкретных типов, параметрических типов, боксинга и устойчивости типов.