Louis Abraham's Home Page

Exact cover and applications

16 Feb 2023

A few months ago, I bought a set of pentominos because it reminded me of a book I had read as a child: The Wright 3. I solved a few configurations but soon my computer scientist education and my laziness took over and whispered: exact cover. If you are looking for a set cover solver, simply use from malib import exact_cover.

What is exact cover?

Exact cover is one of the most classical problems in computer science. It appears very naturally in, like, many areas – ok, naturally as in logical games.

Exact cover problems always have that same vibe of fulfilling existence and uniqueness constraints. For example, in Sudoku, you are required to fill squares such that every column (resp. line and area) contains (existence) exactly one (uniqueness) digit of each type.

More formally, you are required to cover some constraints by choosing pieces. A piece is a choice that will fulfill some constraints. An additional property of exact cover is that in the solution, each constraint must be fulfilled by exactly one piece.

It is usual to associate to a piece the set of constraints it fulfills, as we will do in the code. We can also associate a constraint with the piece that can fulfill it, and we will also do it in the code. Actually, I think the best way to visualize a set cover instance is as a bipartite graph:

Here, the problem with pieces $A=\{1\}, B=\{2, 4\}, C=\{3, 5\}, D=\{5\}$ is represented by the graph. The left part is the set of pieces, the right part is the set of constraints. A piece is connected to a constraint if it fulfills it. The solution is obviously $A, B, C$ since it covers $\{1, 2, 3, 4, 5\}$ without any overlap.

Encoding Sudoku

In Sudoku, the pieces are meant to represent putting some digit in a square. Filling a square with a digit impacts 4 constraints:

Hence the pieces are sets of 4 elements.

from collections import defaultdict

piece_to_constraints = defaultdict(set)
for i in range(9):
    for j in range(9):
        for digit in range(1, 10):
            piece_to_constraints[digit, i, j].add(("square", i, j))
            piece_to_constraints[digit, i, j].add(("line", i, digit))
            piece_to_constraints[digit, i, j].add(("column", j, digit))
            piece_to_constraints[digit, i, j].add(("area", i//3, j//3, digit))

We just defined the constraints of any Sudoku grid. Solving a Sudoku means completing a set of already chosen pieces. This leads us to the next section: how can we select a piece? How can we represent the state so that this operation is efficient as well as the backtracking?

Solving exact cover

There is only one way to solve exact cover, despite what you can hear elsewhere: backtracking. Backtracking simply means selecting pieces while we can. If we cannot anymore, either we have fulfilled all constraints or there are constraints that cannot be fulfilled anymore, and we change the last choice we made.

bipartite graph

However, there are a lot of optimizations that can improve the performance of backtracking, mostly by making choices in a smarter way or making the data structures more efficient.

To optimize the choices, the most standard heuristic is to first try to fulfill the constraints that are the most difficult to fulfill. This is called minimum remaining values as the most difficult constraints to fulfill are the ones that have the least number of possible pieces.

As for the data structures, the most famous algorithm is Dancing Links. It is a very efficient algorithm that uses a doubly linked list to represent the constraints and the pieces. However, I find it hard to implement correctly. The implementation I present is inspired from Ali Assaf’s post with a few tweaks to make it simpler by using only one data structure instead of two.

Instead, we are going to use a constraint_to_pieces dictionary. It maps any constraint key to the set of pieces that fulfill this constraint. Pieces are represented with both their key (only used to print the solution) and the set (represented with a tuple) of constraints they fulfill. This data structure allows us to:

Let us see what the data structure looks like on a simple example:

The constraints that are tracked are shown on the right. They link back to pieces that fulfill them. The pieces are shown on the left as illustration but they are only stored through their constraints.

Here is how to convert our piece_to_constraints dictionary to a constraint_to_pieces dictionary while using the same tuple for every instance of a piece:

from collections import defaultdict, namedtuple

Piece = namedtuple("Piece", "name constraints")

def convert(piece_to_constraints):
    """Take a dictionnary with iterable keys
    Return a dictionnary constraint_to_pieces containing sets of Piece
    """
    constraint_to_pieces = defaultdict(set)
    for piece_name, constraints in piece_to_constraints.items():
        piece = Piece(piece_name, tuple(constraints))
        for constraint in constraints:
            constraint_to_pieces[constraint].add(piece)
    return constraint_to_pieces

When you think about it, Piece tuples are just like pointers that allow to retrieve the constraints of a piece. This is why we can use the same tuple for every instance of a piece. We could have used a dictionnary as in Ali Assaf’s version but this solution makes it clear that we do not need to store the dict index.

We can now write a select function that selects a piece to add to the solution. This function takes the constraints that this piece fulfills as an argument.

What the function does is remove each contraint that the piece fulfills and adapt the problem in consequence. We are going to remove constraint 3 in the example above to illustrate the inner loop with constraint = 3, even though it is not possible to only remove $3$.

If constraint $3$ is fulfilled, then pieces $C$ and $D$ cannot be used anymore (in practice, one of them has been selected but we do not really care for this example). Hence, we have to first remove those from the constraints $2$ and $5$ in purple (constraint_to_pieces[other_constraint].remove(piece)). Finally we delete constraint 3 in green with constraint_to_pieces.pop(constraint). However, we store pieces $C$ and $D$ in the other_pieces list so that we can restore the state of the problem when we backtrack.

def select(constraint_to_pieces, constraints):
    """Suppose a list of constraints is fulfilled (for example because we added a piece).
    This implies that all pieces that have these constraints cannnot be used anymore so we remove them.
    Finally we remove the constraints from the dictionnary as well and return the pieces that were removed to be able to backtrack.
    """
    other_pieces = []
    for constraint in constraints:
        # this constraint is now fulfilled:
        # all pieces that have this constraint can be removed from the other constraints
        for piece in constraint_to_pieces[constraint]:
            for other_constraint in piece.constraints:
                if other_constraint != constraint:
                    constraint_to_pieces[other_constraint].remove(piece)
        # remove the constraint and store it for backtracking
        other_pieces.append(constraint_to_pieces.pop(constraint))
    return other_pieces

The other_pieces list that we return is useful to restore the state of the problem when we backtrack. We do this with a deselect function that does exactly the same operations in reverse:

def deselect(constraint_to_pieces, constraints, other_pieces):
    for constraint in reversed(constraints):
        constraint_to_pieces[constraint] = other_pieces.pop()
        for other_piece in constraint_to_pieces[constraint]:
            for other_constraint in other_piece.constraints:
                if other_constraint != constraint:
                    constraint_to_pieces[other_constraint].add(other_piece)

Finally, we can write a solve function that uses these two functions to solve the problem:

def solve(constraint_to_pieces, solution=None):
    if solution is None:
        solution = []
    if not constraint_to_pieces:
        # make the solution a tuple so that it cannot be modified anymore
        yield tuple(solution)
        return

    # heuristic to minimize the branching factor
    constraint = min(constraint_to_pieces, key=lambda c: len(constraint_to_pieces[c]))
    for piece in list(constraint_to_pieces[constraint]):
        solution.append(piece.name)
        other_pieces = select(constraint_to_pieces, piece.constraints)
        for s in solve(constraint_to_pieces, solution):
            yield s
        deselect(constraint_to_pieces, piece.constraints, other_pieces)
        solution.pop()

solve is a generator so that we can find all solutions.

piece_to_constraints = {"A": {1}, "B": {2, 4}, "C": {2,3,5}, "D": {3, 5}}
list(solve(convert(piece_to_constraints)))
# [('A', 'B', 'D')]

Solving Sudoku

We can now solve Sudoku with this algorithm. We just need to convert the pieces to constraints dictionary to a constraint to pieces dictionary:

def encode(grid):
    # same code as above
    piece_to_constraints = defaultdict(set)
    for i in range(9):
        for j in range(9):
            for digit in range(1, 10):
                piece_to_constraints[digit, i, j].add(("square", i, j))
                piece_to_constraints[digit, i, j].add(("line", i, digit))
                piece_to_constraints[digit, i, j].add(("column", j, digit))
                piece_to_constraints[digit, i, j].add(("area", i // 3, j // 3, digit))

    constraint_to_pieces = convert(piece_to_constraints)

    # add the pieces that are already in the grid
    solution = []
    for i in range(9):
        for j in range(9):
            if grid[i][j] != "X":
                digit = int(grid[i][j])
                select(constraint_to_pieces, piece_to_constraints[digit, i, j])
                solution.append((digit, i, j))
    
    return constraint_to_pieces, solution

# hardest sudoku ever
grid = """8XXXXXXXX
XX36XXXXX
X7XX9X2XX
X5XXX7XXX
XXXX457XX
XXX1XXX3X
XX1XXXX68
XX85XXXX1
X9XXXX4XX
""".split()

We can now solve the problem:

constraint_to_pieces, solution = encode(grid)
next(solve(constraint_to_pieces, solution))
solution_grid = [["X" for _ in range(9)] for _ in range(9)]
for digit, i, j in solution:
    solution_grid[i][j] = digit
for i in range(9):
    print("".join(str(solution_grid[i][j]) for j in range(9)))

819273645
423651879
675498213
254387196
136945782
987162534
541729368
768534921
392816457

We can even count the number of solutions (you have to re-run encode because the state was modified by solve):

constraint_to_pieces, solution = encode(grid)
print(sum(1 for _ in solve(constraint_to_pieces)))
# returns 80

Solving Kakuro

Kakuro is another japanese game. The goal is to fill a grid with digits so that each digit is used at most once in each line and columns and such that the sum of the numbers in each row and column fulfills some constraints.

You can find a whole solver I wrote on github. I will not write the whole solver here but only discuss how we can formulate it as an exact cover problem.

It is interesting that there is a unicity constraint, but no existence constraint. Actually you cannot use all digits except if the line has 9 slots. The solution is to add some dummy slots so that the new problem only has existence and unicity constraints.

Furthermore, we need to ensure that the sum of a line is equal to some predefined value. An elegant solution is to create some pieces that will both fill the dummy slots and ensure this sum constraint.

As an example, suppose we need to fill a line with 7 slots and total 40. We can add 2 dummy slots and their sum will be $45 - 40 = 5$. There are 2 possibilities: $\{1,4\}$ and $\{2, 3\}$. Hence we get again a problem with existence and unicity constraints by creating the dummy pieces.

To enumerate all parititions of a number in distinct integers, a simple dynamic programming approach works well.

Solving pentominos and pentacubes

Finally, we see how to solve pentominos. I kept this example for the end because the word “piece” would be confusing.

The pieces used by our algorithm are choices, hence they correspond to a pentominos in some orientation and at a definite position.

The constraints are simple:

Here, pieces can be represented by their constraints. For example, an example piece can be: ((0, 0, 1), (0, 1, 0), (0, 1, 1), (0, 2, 1), (0, 2, 2), 'F'). The first 5 elements are the coordinates of the pentacube and the last element is the letter of the pentacube. They all correspond to constraints since all coordinated must be filled and all pentacubes must be used once.

All the work is to actually compute orientations and positions of the different pentominos, which becomes a bit harder for pentacubes.

The whole code is available on github and can find solutions for both 2D and 3D shapes:

WVVVFFYYYY
WWXVFFUYUI
WWXVTFUUUI
PXXXTNNNNI
PPLTTTZNZI
PPLLLLZZZI