A Dancing Links implementation for set multicover
In my previous posts on Miracle Sudoku and on my variant sudoku solver package, I mentioned implementing Knuth's dancing links algorithm for set multicover problems.
At its heart, the algorithm is a straightforward depth-first search, but which exploits the structure of a set multicover problem to efficiently eliminate branches of the search tree.
Since we're leveraging the structure of the problem, it helps to put that structure into a class.
A covering problem instance consists of a set of constraints and a set of choices. The constraints simply form an abstract set, with each choice being a subset of the constraints. The intended interpretation is that a given choice satisfies its constraints.
The goal of an exact cover problem is to pick a minimal set of choices to satisfy every constraint.
To put things more concretely, we can take the same approach as Knuth and view the covering problem instance as a 2D matrix, where columns correspond to constraints and rows correspond to choices.
Enumerating the columns \(c\) and rows \(r\), we get a matrix representation \(A\) of the covering problem by placing a \(1\) in \(A(r, c)\) if choice \(r\) belongs to constraint \(c\).
The multiset multicover problem generalizes this in two ways:
- The constraints are a multiset, with the intended interpretation being that a constraint which appears with multiplicity should be satisfied with the same multiplicity.
- The choices are multisets, with the intended interpretation being that a constraint appearing in a choice with multiplicity is satisfied by that choice with that multiplicity.
Including just the first generalization gives a set multicover problem, including the second gives a multiset multicover problem.
A third generalization, which isn't required for traditional Sudoku but which is required for combinatorial puzzles like the \(n\) queens problem and which we'll require for our variant Sudoku module, is to allow for some constraints to be optional: meaning that they can fail to be satisfied, but must not be satisfied more than the specified limit of times.
Our SetMultiCover class will take as its argument a dictionary of choices, where the keys are the choices and the corresponding values are sets of constraints satisfied by the choice.
class SetMultiCover:
def __init__(
self, choices, initial_choices=(), multiplicities=None, must_satisfy=None
):
For the dancing links aglorithm, we also need to operate along the rows: given a constraint, we need to know all the choices which satisfy that constraint. Dancing links uses this to efficiently eliminate other choices from the search space once a constraint is satisfied. So we'll need to recover this row-to-column, or constraint-to-choices, representation:
self.choice_to_constraints = choices
self.constraint_to_choices = defaultdict(set)
for choice in self.choice_to_constraints:
for constraint in self.choice_to_constraints[choice]:
self.constraint_to_choices[constraint].add(choice)
To accomodate the multicover generalization, we need to keep counters of the number of times a constraint still needs to be satisfied (for each constraint, we'll decrement this value as we satisfy it). The optional constraint generalization is easy, we simply keep separate counters for the required and the optional constraints.
To make things simpler, when a list of required constraints isn't provided, we revert to the regular set multicover problem where all constraints are required. This might seem a bit backwards, but it's the correct approach to take instead of requiring that a user who wnats a regular set multicover problem provide a redundant list of all the constraints.
# Dump constraints into required/optional Counters
self.required_constraints = Counter()
self.optional_constraints = Counter()
if must_satisfy:
self.required_constraints.update(must_satisfy)
self.optional_constraints.update(
constraint
for constraint in self.constraint_to_choices.keys()
if constraint not in self.required_constraints
)
if multiplicities is not None:
for key, value in multiplicities.items():
if key in self.required_constraints:
self.required_constraints[key] = value
else:
self.optional_constraints[key] = value
else:
# Default to all constraints must be satisfied.
self.required_constraints.update(self.constraint_to_choices.keys())
if multiplicities is not None:
for key, value in multiplicities.items():
self.required_constraints[key] = value
# Dictionary mapping constraint to which counter it belongs.
self.which_counter = dict()
self.which_counter.update(
(constraint, self.required_constraints)
for constraint in self.required_constraints
)
self.which_counter.update(
(constraint, self.optional_constraints)
for constraint in self.optional_constraints
)
As part of the implementation, we'll go ahead and build the solution by calling in to the private _dancing_links solver we're about to write. To make things a bit more Pythonic, we'll be making the objects of this class into iterators, so we provide the standard iterator methods.
# Initialize partial solution being constructed.
self.solution = []
# Dump in initial solutions and solve.
try:
for choice in initial_choices:
self._make_choice(choice)
self.it = self._dancing_links()
except KeyError:
# If we can't make all the initial choices, puzzle is unsolveable.
self.it = iter(())
def __iter__(self):
return self
def __next__(self):
return next(self.it)
Alright, now time to write the _dancing_links method.
Knuth's algorithm proceeds as follows:
1. If the matrix A has no columns, the current partial
solution is a valid solution; terminate successfully.
2. Otherwise choose a column c (deterministically).
3. Choose a row r such that A(r, c) = 1 (nondeterministically).
4. Include row r in the partial solution.
5. For each column j such that A(r, j) = 1,
for each row i such that A(i, j) = 1,
delete row i from matrix A.
delete column j from matrix A.
We execute a depth-first traversal on the search tree to find all the solutions.
The only important optimization is Knuth's "S heuristic":
When choosing a new unsatisfied column, we always pick
the column with the fewest number of 1's (constraint
with the fewest number of choices).
Once we've selected a column, we consider each choice of row in order and recurse down each of those branches of the search tree.
Because we are considering a multicover problem with multiplicity, this column may remain selected, meaning we'd have to make multiple choices of rows from this column. To eliminate branching over different orders of those choices we enforce an order on rows by Python's built in id() serialization.
def _dancing_links(self, selected_column=None, last_choice=None):
if not self.required_constraints:
# Matrix A has no columns, can yield the current partial
# solution and terminate.
yield list(self.solution)
return
if selected_column in self.required_constraints:
# Still working on the same column, so only consider row
# choices in order.
choices = [
choice
for choice in self.constraint_to_choices[selected_column]
if id(choice) > id(last_choice)
]
else:
# Need to select a new column. Pick the column according
# to S heuristic.
constrained_column = min(
self.required_constraints,
key=lambda col: len(self.constraint_to_choices[col]),
)
choices = list(self.constraint_to_choices[constrained_column])
# Recurse on each of the possile choices.
for choice in choices:
self._make_choice(choice)
yield from self._dancing_links(
selected_column=constrained_column, last_choice=choice
)
self._backtrack()
The _make_choice method needs to do a lot of book-keeping: once we've decided to choose a column \(j\), we need to go through the column and look at every corresponding row \(i\) where \(A(i,j)=1\). For each such row \(i\), if that constraint has now been already satisfied its limited number of times, we need to go through it and remove every other column (other_choice) satisfying that constraint.
def _make_choice(self, choice):
self.solution.append(choice)
for constraint in self.choice_to_constraints[choice]:
counter = self.which_counter[constraint]
counter[constraint] -= 1
if not counter[constraint]:
del counter[constraint]
for other_choice in self.constraint_to_choices[constraint]:
for other_constraint in self.choice_to_constraints[other_choice]:
if other_constraint != constraint:
# WARNING: note we introduce an asymmetry
# between constraint_to_chocies and
# choices_to_constraint
self.constraint_to_choices[other_constraint].remove(
other_choice
)
One perhaps ugly feature we've commented on above is that as we apply this algorithm, we're only actually mutating the constraint_to_choices dictionaries, so if you want to modify this algorithm yourself, keep in mind that these dictionaries don't remain symmetric once the algorithm starts running.
However, this symmetry is essential to make it easy to backtrack during our depth-first-search, since we have all the information we need to recover in the choice_to_constraints dictionary we haven't mutated:
def _backtrack(self):
choice = self.solution.pop()
for constraint in self.choice_to_constraints[choice]:
counter = self.which_counter[constraint]
if constraint not in counter:
# restore what we deleted from constraint_to_choices
for other_choice in self.constraint_to_choices[constraint]:
for other_constraint in self.choice_to_constraints[other_choice]:
if other_constraint != constraint:
self.constraint_to_choices[other_constraint].add(
other_choice
)
counter[constraint] += 1
And that's the whole implementation. I'll show how to put it into action to solve variant sudokus and other puzzles in a future post.