Modeling Sudoku Puzzles with Python

The popular Sudoku puzzles which appear daily in newspapers the world over have, lately, attracted the attention of mathematicians and computer scientists. There are many, difficult, unsolved problems about Sudoku puzzles and their generalizations which make them especially interesting to mathematicians. Also, as is well-known, the generalization of the Sudoku puzzle to larger dimension is an NP-complete problem and therefore of substantial interest to computer scientists. In this article we discuss the modeling of Sudoku puzzles in a variety of different mathematical domains. We show how to use existing third-party Python libraries to implement these models. Those implementations, which include translations into the domains of constraint satisfaction, integer programming, polynomial calculus and graph theory, are available in an opensource Python library sudoku.py developed by the authors and available at http://bitbucket.org/matthew/scipy2010


Sudoku puzzles
A Sudoku puzzle is shown near the top of the second column on this page.
To complete this puzzle requires the puzzler to fill every empty cell with an integer between 1 and 9 in such a way that every number from 1 up to 9 appears once in every row, every column and every one of the small 3 by 3 boxes highlighted with thick borders.
Sudoku puzzles vary widely in difficulty. Determining the hardness of Sudoku puzzles is a challenging research problem for computational scientists. Harder puzzles typically have fewer prescribed symbols. However, the number of prescribed cells is not alone responsible for the difficulty of a puzzle and it is not wellunderstood what makes a particular Sudoku puzzle hard, either for a human or for an algorithm to solve.
The Sudoku puzzles which are published for entertainment invariably have unique solutions. A Sudoku puzzle is said to be well-formed if it has a unique solution. Another challenging research problem is to determine how few cells need to be filled for a Sudoku puzzle to be well-formed. Well-formed Sudoku with 17 symbols exist. It is unknown whether or not there exists a wellformed puzzle with only 16 clues. In this paper we consider all Sudoku puzzles, as defined in the next paragraph, not only the well-formed ones.

sudoku.py
With sudoku.py, the process of building models of Sudoku puzzles, which can then be solved using algorithms for computing solutions of the models, is a simple matter. In order to understand how to build the models, first it is necessary to explain the two different representations of Sudoku puzzles in sudoku.py.
The dictionary representation of a puzzle is a mapping between cell labels and cell values. Cell values are integers in the range {1, . . . , n 2 } and cell labels are integers in the range {1, . . . , n 4 }. The labeling of a Sudoku puzzle of boxsize n starts with 1 in the top-left corner and moves along rows, continuing to the next row when a row is finished. So, the cell in row i and column j is labeled (i − 1)n 2 + j.
For example, the puzzle from the introduction can be represented by the dictionary >>> d = {1: 2, 2: 5, 5: 3, 7: 9, 9: 1, . In practice, however, the user mainly interacts with sudoku.py either by creating specific puzzles instances through input of puzzle strings, directly or from a text file, or by using generator functions.
The string representation of a Sudoku puzzle of boxsize n is a string of ascii characters of length n 4 . In such a string a period character represents an empty cell and other ascii characters are used to specify assigned values. Whitespace characters and newlines are ignored when Puzzle objects are built from strings.
A possible string representation of the puzzle from the introduction is: The first argument to random_puzzle is the number of prescribed cells in the puzzle. Solving of puzzles in sudoku.py is handled by the solve function. This function can use a variety of different algorithms, specified by an optional model keyword argument, to solve the puzzle. Possible values are CP for constraint propagation, lp for linear programming, graph to use a node coloring algorithm on a graph puzzle model and groebner to solve a polynomial system model via a Groebner basis algorithm. The default behavior is to use constraint propagation.

Models
In this section we introduce several models of Sudoku and show how to use existing Python components to implement these models. The models introduced here are all implemented in sudoku.py. Implementation details are discussed in this section and demonstrations of the components of sudoku.py corresponding to each of the different models are given.

Constraint models
Constraint models for Sudoku puzzles are discussed in [Sim05]. A simple model uses the AllDifferent constraint. A constraint program is a collection of constraints. A constraint restricts the values which can be assigned to certain variables in a solution of the constraint problem. The AllDifferent constraint restricts variables to having mutually different values.
Modeling Sudoku puzzles is easy with the AllDifferent constraint. To model the empty Sudoku puzzle (i.e. the puzzle with no clues) a constraint program having an AllDifferent constraint for every row, column and box is sufficient.
For example, if we let x i ∈ {1, . . . , n 2 } for 1 ≤ i ≤ n 4 , where x i = j means that cell i gets value j then the constraint model for a Sudoku puzzle of boxsize n = 3 would include constraints: These constraints ensure that, respectively, the variables in the first row, column and box get different values. The Sudoku constraint model in sudoku.py is implemented using python-constraint v1.1 by Gustavo Niemeyer. This open-source library is available at http://labix.org/pythonconstraint.
With python-constraint a Problem having variables for every cell {1, . . . , n 4 } of the Sudoku puzzle is required. The list of cell labels is given by the function cells in sudoku.py. Every variable has the same domain {1, . . . , n 2 } of symbols. The list of symbols in sudoku.py is given by the symbols function.
The Problem member function addVariables provides a convenient method for adding variables to a constraint problem object.
>>> from constraint import Problem >>> from sudoku import cells, symbols >>> cp = Problem() >>> cp.addVariables(cells(n), symbols(n)) The AllDifferent constraint in python-constraint is implemented as AllDifferentConstraint(). The addConstraint(constraint, variables) member function is used to add a constraint on variables to a constraint Problem object. So, to build an empty Sudoku puzzle constraint model we can do the following. Here the functions cells_by_row, cells_by_col and cells_by_box give the cell labels of a Sudoku puzzle ordered, respectively, by row, column and box. These three loops, respectively, add to the constraint problem object the necessary constraints on row, column and box variables. To extend this model to a Sudoku puzzle with clues requires additional constraints to ensure that the values assigned to clue variables are fixed. One possibility is to use an ExactSum constraint for each clue.
The ExactSum constraint restricts the sum of a set of variables to a precise given value. We can slightly abuse the ExactSum constraint to specify that certain individual variables are given certain specific values. In particular, if the puzzle clues are given by a dictionary d then we can complete our model by adding the following constraints.
>>> from constraint import ExactSumConstraint as Exact >>> for cell in d: To solve the Sudoku puzzle now can be done by solving the constraint model cp. The constraint propogation algorithm of python-constraint can be invoked by the getSolution member function.

Graph models
A graph model for Sudoku is presented in [Var05]. In this model, every cell of the Sudoku grid is represented by a node of the graph. The edges of the graph are given by the dependency relationships between cells. In other words, if two cells lie in the same row, column or box, then their nodes are joined by an edge in the graph. In the graph model, a Sudoku puzzle is given by a partial assignment of colors to the nodes of the graph. The color assigned to a node corresponds to a value assigned to the corresponding cell. A solution of the puzzle is given by a coloring of the nodes with colors {1, . . . , n 2 } which extends the original partial coloring. A node coloring of the Sudoku graph which corresponds to a completed puzzle has the property that adjacent vertices are colored differently. Such a node coloring is called proper.
The Sudoku graph model in sudoku.py is implemented using networkx v1.1. This open-source Python graph library is available at http://networkx.lanl.gov/ Modeling an empty Sudoku puzzle as a networkx.Graph object requires nodes for every cell and edges for every pair of dependent cells. To add nodes (respectively, edges) to a graph, networkx provides member functions add_nodes_from (respectively, add_edges_from). Cell labels are obtained from sudoku.py's cells function.
>>> import networkx >>> g = networkx.Graph() >>> g.add_nodes_from(cells(n)) Dependent cells are computed using the dependent_cells function. This function returns the list of all pairs (x, y) with x < y such that x and y either lie in the same row, same column or same box.

>>> s = solve(p, model = 'graph')
Polynomial system models The graph model above is introduced in [Var05] as a prelude to modeling Sudoku puzzles as systems of polynomial equations. The polynomial system model in [Var05] involves variables x i for i ∈ {1, . . . , n 4 } where x i = j is interpreted as the cell with label i being assigned the value j.
The Sudoku polynomial-system model in sudoku.py is implemented using sympy v0.6.7. This open-source symbolic algebra Python library is available at http://code.google.com/p/sympy/ Variables in sympy are Symbol objects. A sympy.Symbol object has a name. So, to construct the variables for our model, first we map symbol names onto each cell label.
>>> from sudoku import cell_symbol_name >>> def cell_symbol_names(n): ... return map(cell_symbol_name, cells(n)) Now, with these names for the symbols which represent the cells of our Sudoku puzzle, we can construct the cell variable symbols themselves.
>>> from sympy import Symbol >>> def cell_symbols(n): ... return map(Symbol, cell_symbol_names(n)) Finally, with these variables, we can build a Sudoku polynomial system model. This model is based on the graph model of the previous section. There are polynomials in the system for every node in the graph model and polynomials for every edge. The role of node polynomial F(x i ) is to ensure that every cell i is assigned a number from {1, . . . , n 2 } : Node polynomials, for a sympy.Symbol object x are built as follows.

return reduce(mul,[(x-s) for s in symbols(n)])
The edge polynomial G(x i , x j ) for dependent cells i and j, ensures that cells i and j are assigned different values. These polynomials have the form. : In sympy, we build edge polynomials from the node polynomial function F. >>> from sympy import cancel, expand >>> def G(x,y,n): ... return expand(cancel((F(x,n)-F(y,n))/(x-y))) The polynomial model for the empty Sudoku puzzle consists of the collection of all node polynomials for nodes in the Sudoku graph and all edge polynomials for pairs (x,y) in dependent_symbols(n). The dependent_symbols function is simply a mapping of the sympy.Symbol constructor onto the list of dependent cells. Specifying a Sudoku puzzle requires extending this model by adding polynomials to represent clues. According to the model from [Var05], if D is the set of fixed cells (i.e. cell label, value pairs) then to the polynomial system we need to add polynomials The sympy implementation of a Groebner basis algorithm can be used to find solutions of this polynomial system. The Groebner basis depends upon a variable ordering, here specified as lexicographic. Other orderings, such as degree-lexicographic, are possible.
>>> from sympy import groebner >>> h = groebner(g, cell_symbols(n), order = 'lex') The solution of the polynomial system g is a system of linear equations in the symbols x i which can be solved by the linear solver from sympy.

>>> from sympy import solve as lsolve >>> s = lsolve(h, cell_symbols(n))
To use the polynomial-system model to find a solution to Puzzle instance p call the solve function with the keyword argument model = groebner. >>> s = solve(p, model = 'groebner')

Integer programming models
In [Bar08] a model of Sudoku as an integer programming problem is presented. In this model, the variables are all binary.
x i jk ∈ {0, 1} Variable x i jk represents the assignment of symbol k to cell (i, j) in the Sudoku puzzle.
The integer programming (IP) model has a set of equations which force the assignment of a symbol to every cell.
Other equations in the IP model represent the unique occurence of every symbol in every column: x i jk = 1, 1 ≤ j ≤ n, 1 ≤ k ≤ n every symbol in every row: n ∑ j=1 x i jk = 1, 1 ≤ i ≤ n, 1 ≤ k ≤ n and every symbol in every box: The Sudoku IP model is implemented in sudoku.py using pyglpk v0.3 by Thomas Finley. This open-source mixed integer/linear programming Python library is available at http: //tfinley.net/software/pyglpk/ In pyglpk, an integer program is represented by the matrix of coefficients of a system of linear equations. Two functions of sudoku.py provide the correct dimensions of the coefficient matrix.

Experimentation
In this section we demonstrate the use of sudoku.py for creating Python scripts for experimentation with Sudoku puzzles. For the purposes of demonstration, we discuss, briefly, enumeration of Shidoku puzzles, coloring the Sudoku graph and the hardness of random puzzles.

Enumerating Shidoku
Enumeration of Sudoku puzzles is a very difficult computational problem, which has been solved by Felgenhauer and Jarvis in [Fel06]. The enumeration of Shidoku, however, is easy. To solve the enumeration problem for Shidoku, using the constraint model implemented in sudoku.py, takes only a few lines of code and a fraction of a second of computation.

Coloring the Sudoku graph
As discussed above in the section on "Graph models", a completed Sudoku puzzle is equivalent to a minimal proper node coloring of the Sudoku graph. We have experimented with several different node coloring algorithms to see which are more effective, with respect to minimizing the number of colors, at coloring the Sudoku graph. Initially, we used Joseph Culberson's graph coloring programs (http://webdocs.cs.ualberta.ca/~joe/Coloring/index.html) by writing Sudoku puzzle graphs to a file in Dimacs format (via the dimacs_string function of sudoku.py).
Of those programs we experimented with, the program implementing the saturation degree algorithm (DSatur) of Brelaz from [Bre79] seemed most effective at minimizing the number of colors.
Motivated to investigate further, with sudoku.py we implemented a general node coloring algorithm directly in Python which can reproduce the DSatur algorithm as well as several other node coloring algorithms.
Our node coloring function allows for customization of a quite general scheme. The behavior of the algorithm is specialized by two parameters. The nodes parameter is an iterable object giving a node ordering. The choose_color parameter is a visitor object which is called every time a node is visited by the algorithm.
Several node orderings and color choice selection schemes have been implemented. The simplest sequential node coloring algorithm can be reproduced, for example, by assigning nodes = InOrder and choose_color = first_available_color. A random ordering on nodes can be acheived instead by assigning nodes = RandomOrder. Importantly for our investigations, the node ordering is given by an iterable object and so, in general, can reflect upon to current graph state. This mean that online algorithms like the DSatur algorithm can be realized by our general node coloring scheme. The DSatur algorithm is obtained by assigning nodes = DSATOrder and choose_color = first_available_color.

Hardness of random puzzles
We introduced the random_puzzle function in the introduction. The method by which this function produces a random puzzle is fairly simple. A completed Sudoku puzzle is first generated by solving the empty puzzle via constraint propagation and then from this completed puzzle the appropriate number of clues is removed.
An interesting problem is to investigate the behavior of different models on random puzzles. A simple script, available in the investigations folder of the source code, has been written to time the solution of models of random puzzles and plot the timings via matplotlib.
Two plots produced by this script highlight the different behavior of the constraint model and the integer programming model.
The first plot has time on the vertical axis and the number of clues on the horizontal axis. From this plot it seems that the constraint propogation algorithm finds puzzles with many or few clues easy. The difficult problems for the constraint solver appear to be clustered in the range of 20 to 35 clues.
A different picture emerges with the linear programming model. With the same set of randomly generated puzzles it appears that the more clues the faster the solver finds a solution.

Conclusions and future work
In this article we introduced sudoku.py, an open-source Python library for modeling Sudoku puzzles. We discussed several models of Sudoku puzzles and demonstrated how to implement these models using existing Python libraries. A few simple experiments involving Sudoku puzzles were presented.
Future plans for sudoku.py are to increase the variety of models. Both by allowing for greater customization of currently implemented models and by implementing new models. For example, we can imagine several different Sudoku models as constraint programs beyond the model presented here. Another approach is to model Sudoku puzzles as exact cover problems and investigate the effectiveness of Knuth's dancing links algorithm. Also important to us is to compare all our models with models [Lyn06] from satisfiability theory. In [Kul10] a general scheme is presented which is highly effective for modeling Sudoku.
There are great many interesting, unsolved scientific problems involing Sudoku puzzles. Our hope is that sudoku.py can become a useful tool for scientists who work on these problems.