Exact cover I

,

1. Introduction

Karp (1972) introduces the name ‘exact cover’ for the following decision problem:

INPUT: Family $$\{S_j\}$$ of subsets of a set $$\{u_i, i=1,2,\ldots,t\}$$

PROPERTY: There is a subfamily $$\{T_h\} ⊆ \{S_j\}$$ such that the sets $$T_h$$ are disjoint and $$∪T_h = ∪S_j = \{u_i, i=1,2,\ldots,t\}$$

(The problem being to determine, for an arbitrary input, whether it possesses the property.) Karp goes on to prove that this problem is NP-complete, using a reduction from chromatic number, which has a reduction from satisfiability with at most three literals per clause (3SAT for short), which in turn has a reduction from satisfiability, which he proves is NP-complete directly.

I’m going to call the subsets $$S_j$$ ‘choices’ and the $$u_i$$ ‘constraints’, and to say that $$S_j$$ ‘satisfies’ $$u_i$$ iff $$u_i ∈ S_j$$. The problem is then to determine if there is a set of choices that satisfies all the constraints exactly once each.1

Instances of this problem often arise in the solution of combinatorial problems, except that in this context we are interested in computing or counting the solutions, not just deciding if the problem has a solution. For example, consider the problem of tiling a region of 60 squares with the twelve pentominoes:2

In this problem the constraints are the 60 squares to be tiled (each square must be covered exactly once), and the 12 pentominoes to be placed (each must be placed exactly once). The choices are the different locations and orientations at which each pentomino can be placed.

Pentomino tiling problems were among the first combinatorial problems to be solved using a backtracking algorithm. Scott (1958) describes how he programmed the MANIAC II to find the 65 essentially different tilings of an 8×8 square (minus the central 2×2 square) by the twelve pentominoes, two of which are shown here.3

2. ‘Algorithm X’

How might this be solved? Well, Knuth (2000) gives the following algorithm, which he calls ‘Algorithm X’:4

1. If the set of unsatisfied constraints is empty, the problem is solved. Yield the solution and return.
2. Otherwise, choose an unsatisfied constraint c.
3. If there are no choices that satisfy this constraint, the problem cannot be solved from this position. Return.
4. Otherwise, for each choice i that satisfies c:
1. For each constraint j satisfied by i:
1. Remove j from the set of unsatisfied constraints.
2. For each choice k that satisfies j, remove k from the set of choices.
2. Recursively run this algorithm on the reduced set of choices and constraints.
3. For each constraint j satisfied by i:
1. Restore j to the set of unsatisfied constraints.
2. For each choice k that satisfies j, restore k to the set of choices.

Knuth comments on this:

Algorithm X is simply a statement of the obvious trial-and-error approach. (Indeed, I can’t think of any other reasonable way to do the job, in general.)

3. Python implementation

An obvious question is, why Python, since combinatorial problems are computationally intensive, and Python is pretty slow? Well, Python’s combination of simplicity, terseness, and flexibility makes it a good language for experimenting with and modifying code. And Python’s native sets can be used to implement the sets of constraints and choices, meaning that I will be able to represent the constraints and choices themselves by any Python objects, avoiding the necessity for complex encoding and decoding. The algorithm needs:

• A collection of unsatisfied constraints (see step 2). This must support efficient addition and removal of constraints (at steps 4.1.1 and 4.3.1). This is implemented as a set, unsatisfied.
• An efficient way to find the choices that satisfy a constraint (at steps 4, 4.1.2 and 4.3.2). This is implemented as a dictionary, choices, mapping each constraint to the set of choices satisfying that constraint.
• An efficient way to find the constraints satisfied by a choice (at steps 4.1 and 4.3). This is implemented as a dictionary, constraints, mapping each choice to the set of constraints satisfied by that choice.
from collections import defaultdict
from collections.abc import Iterator
from random import shuffle

class ExactCover(Iterator):
"""An iterator that yields solutions to an EXACT COVER problem.

An EXACT COVER problem consists of a set of "choices". Each choice
satisfies one or more "constraints". Choices and constraints may
be represented by any hashable Python objects, with the proviso
that all choices must be distinct, as must all constraints.

A solution is a list of choices such that each constraint is
satisfied by exactly one of the choices in the solution.

"""
# This implements Donald Knuth's "Algorithm X"
# http://lanl.arxiv.org/pdf/cs/0011047.pdf

def __init__(self, constraints, initial=(), random=False):
"""Construct an ExactCover solver,

Required argument:
constraints -- A map from each choice to an iterable of
constraints satisfied by that choice.

Optional arguments:
initial -- An iterable of choices that must appear in every
solution. (Default: no choices.)
random -- Generate solutions in random order? (Default: False.)

For example:

>>> next(ExactCover(dict(A = [1, 4, 7],
...                      B = [1, 4],
...                      C = [4, 5, 7],
...                      D = [3, 5, 6],
...                      E = [2, 3, 6, 7],
...                      F = [2, 7])))
['B', 'D', 'F']

"""
self.random = random
self.constraints = constraints

# A map from constraint to the set of choices that satisfy
# that constraint.
self.choices = defaultdict(set)
for i in self.constraints:
for j in self.constraints[i]:

# The set of constraints which are currently unsatisfied.
self.unsatisfied = set(self.choices)

# The partial solution currently under consideration,
# implemented as a stack of choices.
self.solution = []

# Make all the initial choices.
try:
for i in initial:
self._choose(i)
self.iter = self._solve()
except KeyError:
# Initial choices were contradictory, so there are no solutions.
self.iter = iter(())

def __next__(self):
return next(self.iter)

def _solve(self):
if not self.unsatisfied:
# No remaining unsatisfied constraints.
yield list(self.solution)
return

# Pick the constraint with the fewest remaining choices
# (Knuth's "S heuristic").
best = min(self.unsatisfied, key=lambda j:len(self.choices[j]))
choices = list(self.choices[best])
if self.random:
shuffle(choices)

# Try each choice in turn and recurse.
for i in choices:
self._choose(i)
yield from self._solve()
self._unchoose(i)

def _choose(self, i):
"""Make choice i; mark constraints satisfied; and remove any
choices that clash with it.

"""
self.solution.append(i)
for j in self.constraints[i]:
self.unsatisfied.remove(j)
for k in self.choices[j]:
for l in self.constraints[k]:
if l != j:
self.choices[l].remove(k)

def _unchoose(self, i):
"""Unmake choice i; restore constraints and choices."""
last = self.solution.pop()
assert i == last
for j in self.constraints[i]:
for k in self.choices[j]:
for l in self.constraints[k]:
if l != j:


This code is pretty compact, but it has a small flaw, namely the selection of the constraint with the fewest remaining choices (Knuth’s ‘S heuristic’). This needs to take time $$O(\log n)$$ but here it takes $$O(n)$$ and so will become a bottleneck as the number of constraints increases. However, it’s not easy to make a data structure that supports the necessary operations efficiently.5 I’ll try to revisit this in a later article in this series.

4. Pentomino tiling

First, input and output. It’s a good idea to arrange for these to be easily readable and checkable. I don’t want to have to type in lists of coordinates of tiles in pentominoes; I’d much rather draw a picture:

TETROMINOES = '''
IIII  LLL  OO  TTT  ZZ
L    OO   T    ZZ
'''

PENTOMINOES = '''
I
I     L      N                      Y
I  FF L  PP  N TTT       V   W  X  YY ZZ
I FF  L  PP NN  T  U U   V  WW XXX  Y  Z
I  F  LL P  N   T  UUU VVV WW   X   Y  ZZ
'''


If some other form is more convenient for processing, let the computer do the work:

from collections import defaultdict

_EMPTY = ' '

def encode_polyominoes(picture):
"""Given a "picture" of some polyominoes (in ASCII-art form), return a
list of tuples (name, tiles) for each polyomino, where tiles are
the coordinates (i, j) of its tiles.

>>> sorted(encode_polyominoes(TETROMINOES))
... # doctest: +NORMALIZE_WHITESPACE
[('I', ((1, 0), (1, 1), (1, 2), (1, 3))),
('L', ((1, 6), (1, 7), (1, 8), (2, 6))),
('O', ((1, 11), (1, 12), (2, 11), (2, 12))),
('T', ((1, 15), (1, 16), (1, 17), (2, 16))),
('Z', ((1, 20), (1, 21), (2, 21), (2, 22)))]

"""
pieces = defaultdict(list)
for i, row in enumerate(picture.split('\n')):
for j, c in enumerate(row):
if c != _EMPTY:
pieces[c].append((i, j))
return [(c, tuple(tiles)) for c, tiles in pieces.items()]

def solution_tiles(solution):
"""Generate all the tiles in a polyomino solution."""
return (t for _, tiles in solution for t in tiles)

def solution_bounds(solution):
"""Return the bounds of the region covered by a polyomino solution."""
h, w = map(max, zip(*solution_tiles(solution)))
return h + 1, w + 1

def decode_solution(solution):
"""Format a polyomino solution as a string."""
h, w = solution_bounds(solution)
grid = [[_EMPTY] * w for _ in range(h)]
for c, tiles in solution:
for i, j in tiles:
grid[i][j] = c
return '\n'.join(''.join(row) for row in grid)


Now for the encoding of a tiling problem as an exact cover instance:

ConstraintPython representation
The square with coordinates $$(i, j)$$ must be covered.The tuple (i, j).
The piece named ‘X’ must be placed somewhere.The string 'X'.
ChoicePython representation
Place the piece named ‘X’ so that it covers the squares $$(i_0, j_0)$$, $$(i_1, j_1)$$, …The tuple ('X', ((i0, j0), (i1, j1), ...)).

This encoding has the following useful properties:

1. A solution (a list of choices) is completely self-descriptive.
2. A solution is in the same format as the output of encode_polyominoes, which means I can pass the latter to any function taking a solution. This will come in handy later on when I have code for drawing solutions.
from exactcover import ExactCover
from itertools import product

def polyomino(pieces, region, random=False):
"""Return an iterator that yields the solutions to a polyomino tiling
problem.

Required arguments:
pieces -- an iterable of pieces, each being a tuple whose first
element is the name of the piece and whose second element is a
sequence of coordinates of the tiles in the piece (in an
arbitrary coordinate system).
region -- the region to fill (an iterable of coordinates x, y).

Optional argument:
random -- Generate solutions in random order? (Default: False.)

>>> region = set(product(range(3), range(7))) - {(1, 2)}
>>> sol = min(polyomino(encode_polyominoes(TETROMINOES), region))
>>> print(decode_solution(sol))
OOTTTZZ
OO TZZL
IIIILLL

"""
region = list(region)
region_set = set(region)
constraints = {}
for c, tiles in pieces:
# Origin of piece at the first tile.
oi, oj = tiles[0]
# Adjust tiles so that they are relative to the origin of the piece.
tiles = [(ti - oi, tj - oj) for ti, tj in tiles]
for (oi, oj), ri, rj, rk in product(region, (-1, +1), (-1, +1), (0, 1)):
# Place origin of piece at oi, oj, reflecting vertically
# if ri == -1, horizontally if rj == -1, and diagonally if
# rk == 1.
placing = []
for t in tiles:
p = ri * t[rk] + oi, rj * t[1-rk] + oj
if p not in region_set:
break
placing.append(p)
else:
placing = tuple(sorted(placing))
choice = c, placing
if choice not in constraints:
constraints[choice] = (c,) + placing
return ExactCover(constraints=constraints, random=random)


Now we’re all set to find tiling solutions:

>>> region = set(product(range(8), repeat=2)) - set(product((3, 4), repeat=2))
>>> solution = next(polyomino(encode_polyominoes(PENTOMINOES), region))
>>> print(decode_solution(solution))
YYYYWZZP
LNYWWZPP
LNWWZZPP
LNN  TTT
LLN  FTV
UUXFFFTV
UXXXFVVV
UUXIIIII


5. Drawing solutions

Let’s draw that solution in Scalable Vector Graphics format, using the svgwrite Python package:

from svgwrite import Drawing

def solutions_svg(solutions, filename, columns=1, size=25, padding=5,
colour=lambda _: "white", stroke_colour="black",
stroke_width=3):
"""Format polyomino tilings as an SVG image.

Required arguments:
solutions -- iterable of solutions to the tiling problem, each of
which is a sequence of piece placements, each of which is a
tuple whose first element is the name of the piece, and whose
second element is a sequence of pairs (i, j) giving the
locations of the tiles in the piece.
filename -- where to save the SVG drawing.

Optional arguments:
columns -- number of solutions per row (default: 1).
size -- width and height of each tile (default: 25).
padding -- padding around the image (default: 5)
colour -- function taking a piece name and returning its colour
(default: a function returning white for each piece).
stroke -- stroke colour (default: black).
stroke_width -- width of strokes between pieces (default: 3).

"""
solutions = list(solutions)
h, w = solution_bounds(solutions[0])
rows = (len(solutions) + columns - 1) // columns
drawing_size = (2 * padding + (columns * (w + 1) - 1) * size,
2 * padding + (rows    * (h + 1) - 1) * size)
drawing = Drawing(debug=False, filename=filename, size=drawing_size)
for i, solution in enumerate(solutions):
y, x = divmod(i, columns)
oj = padding + (x * (w + 1)) * size
oi = padding + (y * (h + 1)) * size
group = drawing.g(stroke=stroke_colour, stroke_linecap="round",
stroke_width=0.25)
grid = [[_EMPTY] * w for _ in range(h)]
for c, placing in solution:
piece = drawing.g(fill=colour(c))
for i, j in placing:
grid[i][j] = c
piece.add(drawing.rect((j * size + oj, i * size + oi),
(size, size)))
edges = drawing.path(stroke_width=stroke_width)
for i, j in product(range(h + 1), range(w)):
if ((_EMPTY if i == 0 else grid[i-1][j])
!= (_EMPTY if i == h else grid[i][j])):
edges.push(['M', j * size + oj, i * size + oi, 'l', size, 0])
for i, j in product(range(h), range(w + 1)):
if ((_EMPTY if j == 0 else grid[i][j-1])
!= (_EMPTY if j == w else grid[i][j])):
edges.push(['M', j * size + oj, i * size + oi, 'l', 0, size])
drawing.save()

COLOURS = dict(I="#EEAAAA", F="#DDBB99", L="#CCCC88",
P="#BBDD99", N="#AAEEAA", T="#99DDBB",
U="#88CCCC", V="#99BBDD", W="#AAAAEE",
X="#BB99DD", Y="#CC88CC", Z="#DD99BB")


Then, using:

>>> solutions_svg([solution], 'solution-8x8.svg', colour=COLOURS.get)

produces the image:

How many solutions are there in all?

>>> all_solutions = list(polyomino(encode_polyominoes(PENTOMINOES), region))
>>> len(all_solutions)
520


But these include solutions that are reflections and rotations of each other. To reduce this collection to the essentially different solutions, one approach is to map each solution to a canonical form (by sorting the pieces and the tiles in each piece) and then discard duplicates:

from itertools import starmap

def canonical(solution):
"""Return a canonical version of a polyomino tiling solution."""
h, w = solution_bounds(solution)
all_tiles = sorted(solution_tiles(solution))
orientations = [(False, True)] * 3
if sorted((h - i - 1, j) for i, j in all_tiles) != all_tiles:
# Can't reflect vertically.
orientations[0] = (False,)
if sorted((i, w - j - 1) for i, j in all_tiles) != all_tiles:
# Can't reflect horizontally.
orientations[1] = (False,)
if sorted((j, i) for i, j in all_tiles) != all_tiles:
# Can't reflect diagonally.
orientations[2] = (False,)
solution = sorted(solution)
def oriented(ri, rj, rk):
oriented_solution = []
for c, tiles in solution:
oriented_tiles = []
for i, j in tiles:
if ri: i = h - i - 1 # reflect vertically
if rj: j = w - j - 1 # reflect horizontally
if rk: i, j = j, i   # reflect diagonally
oriented_tiles.append((i, j))
oriented_solution.append((c, tuple(sorted(oriented_tiles))))
return tuple(oriented_solution)
return min(starmap(oriented, product(*orientations)))


This reduces the 520 solutions to the 65 essentially different solutions:

>>> unique_solutions = set(map(canonical, all_solutions))
>>> len(unique_solutions)
65
>>> solutions_svg(unique_solutions, 'all-8x8.svg', columns=13, colour=COLOURS.get)


A. References

• Richard M. Karp (1972). ‘Reducibility Among Combinatorial Problems’. In R. E. Miller and J. W. Thatcher (editors). Complexity of Computer Computations. New York: Plenum. pp. 85–103.
• Donald E. Knuth (2000). ‘Dancing Links’. In J. Davies, B. Roscoe, and J. Woodcock (editors). Millennial Perspectives in Computer Science. Basingstoke: Palgrave Macmillan. pp. 187–214.
• Dana S. Scott (1958). ‘Programming a Combinatorial Puzzle’. Technical Report 1. Princeton University Department of Electrical Engineering.

1.  I find these terms help me understand the problem and the algorithm. In Karp’s setting there are ‘sets’ and ‘subsets’, and in Knuth’s there are ‘rows’ and ‘columns’, but I find it hard to remember what these mean, especially in the context of a puzzle which may itself have rows and columns.

2.  The colours are from this image by R. A. Nonenmacher, licensed under CC-BY-SA.

3.  I haven’t been able to find a copy of this technical report online; all I know is from the summary in Knuth (2000). If you have a copy, or know where I can find one, please let me know.

4.  Knuth presents this as a nondeterministic algorithm, but I think it’s slightly clearer to present it as a deterministic backtracking algorithm that yields all the solutions.

5.  If you’re saying, “it’s easy, dummy, just drop in the SortedSet class from the sortedcontainers package,” it’s because you haven’t tried it yet.