## Friday, August 3, 2018

### Sudoku++, MIP vs Constraint Solvers

In  a variant of the Sudoku puzzle is presented:

Sudoku Variant
Place numbers in the grid so that each row and each column contain the numbers 1 through 5 and the sums of numbers in the outlined regions are all different.

When we want to model this problem, we basically have three sets of constraints:

1. Each row $$i$$ has unique values $$1,\dots,5$$
2. Each column $$j$$ has unique values $$1,\dots,5$$
3. Each region $$r$$ has unique sums
This leads to repeated use of the all-different constraint.

Below I will try to model this problem in different ways:
1. A MIP formulation using GAMS. When using integer programming we don't have an all-different construct. However, with a little bit of effort, we can implement this using linear constraints.
2. A Python implementation using the Z3 SMT solver from Microsoft.
3. A Minizinc implementation using the Gecode constraint programming solver.
The last two approaches allow us to use all-different directly.

In addition, I will discuss how we can verify that there there is only one unique solution to this problem.

### MIP Model

As usual we will use a binary variable : $x_{i,j,k} = \begin{cases} 1 & \text{if cell (i,j) has value k}\\ 0 & \text{otherwise}\end{cases}$

Using this definition, the row and column uniqueness constraints are easy using assignment constraints: \begin{align} & \sum_i x_{i,j,k} = 1 & \forall j,k \\ & \sum_j x_{i,j,k} = 1 & \forall i,k \\ &\sum_k x_{i,j,k} = 1 &\forall i,j\end{align} These assignment constraint effectively implement what is known as all-different constraints in constraint programming. See  for different MIP formulations for these all-different constraints.

Let's indicate a region by $$r$$.  We can calculate the value of a region by: $v_r = \sum_{(i,j)|R_{i,j,r}} \sum_k k\> x_{i,j,k}$ We sum over cells in region $$r$$. Now we have to add an all-different constraint: $\text{all-different}(v_k)$ This is a slightly different all-different constraint from the ones we discussed before. Here we don't enumerate all possible values.

One way to model this, is the following. Interpret the problem as: $v_r \le v_{r'} - 1 \textbf{ or } v_r \ge v_{r'} + 1 \>\> \forall r\ne r'$ We can model this with binary variables $$\delta$$: \begin{align} & v_r \le v_{r'} - 1 + M\delta_{r,r'}\\ & v_r \ge v_{r'} + 1 -M (1-\delta_{r,r'})\\ & \delta_{r,r'} \in \{0,1\} \end{align} Instead of $$\forall r\ne r'$$, we can save some equations and variables by only generating these constraints $$\forall i\lt j$$. We need a good value for $$M$$. If the largest possible value of a region is $$L$$, then $$M=L+1$$. The largest region has 6 cells, so one easy bound is $$L=6 \times 5=30$$. With some more analysis we can reduce this even further, but I am happy with this bound.

Interestingly, we used in this model two different ways to implement all-different constraints: one through assignment constraints and the other by an or constraint.

#### GAMS Model: Mixed Integer Programming

The GAMS model can look like:

 sets   i 'rows' /r1*r5/   j 'columns' /c1*c5/   k 'values' /1*5/   r 'regions' /r1*r11/ ; table grid(i,j) 'encoding of regions'    c1  c2  c3  c4  c5 r1  1   1   2   2   3 r2  4   5   5   5   6 r3  4   7   7   5   6 r4  8   7   5   5  10 r5  8   9   9  10  11 ; set region(i,j,r) '(i,j) is in r'; region(i,j,r) = grid(i,j) = ord(r); alias(r,r1,r2); set rr(r1,r2) 'compare r1,r2'; rr(r1,r2) = ord(r1)

The results look like:

----     88 PARAMETER grid  encoding of regions

c1          c2          c3          c4          c5

r1           1           1           2           2           3
r2           4           5           5           5           6
r3           4           7           7           5           6
r4           8           7           5           5          10
r5           8           9           9          10          11

----     88 PARAMETER xv  cell values

c1          c2          c3          c4          c5

r1           1           3           5           4           2
r2           5           1           4           2           3
r3           3           5           2           1           4
r4           2           4           1           3           5
r5           4           2           3           5           1

----     88 VARIABLE v.L  value of a region

r1   4,    r2   9,    r3   2,    r4   8,    r5  12,    r6   7,    r7  11
r8   6,    r9   5,    r10 10,    r11  1


After playing with my crayons, we have:

This is not a super easy model to solve. Standard Sudoku models are typically solved in the presolve phase, so the MIP solver does not have to any branching. This is not the case with this model:

        Nodes                                         Cuts/
Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

0     0        0.0000    61                      0.0000       99
0     0        0.0000     2                     Cuts: 2      101
0     0        0.0000     3                     Cuts: 6      106
0     0        0.0000     3                 ZeroHalf: 2      108
0     2        0.0000     3                      0.0000      108
Elapsed time = 0.25 sec. (25.37 ticks, tree = 0.01 MB, solutions = 0)
1338   159        0.0000     3                      0.0000    13389
Cuts: 8
2012   123    infeasible                            0.0000    24769
2822    96    infeasible                            0.0000    38178
3072    93        0.0000    67                      0.0000    42948
3486   107        0.0000    38                      0.0000    50692
Cuts: 10
3909   123        0.0000    55                      0.0000    59355
4040   130        0.0000    45                      0.0000    62204
4150    39        0.0000    68                      0.0000    64390
*  4234     0      integral     0        0.0000        0.0000    66184    0.00%
Found incumbent of value 0.000000 after 9.97 sec. (2491.49 ticks)


The solver has to do some real work here: 4234 nodes and 66184 simplex iterations. As this is a feasibility problem, we can stop as soon as we find a feasible integer solution.

#### Uniqueness of the solution

It is easy to add a constraint that checks the uniqueness of this solution. Let's first record the current solution: $$a_{i,j,k} = x^*_{i,j,k}$$. Now add the constraint:$\sum_{i,j,k} a_{i,j,k} x_{i,j,k} \le 5^2-1$ Solve this model and check the status. If this new problem is infeasible, we know our original solution was unique.

When I did this the solver showed:

MIP status(119): integer infeasible or unbounded

In general, I find this always a bit a of disappointing message. The solver should be able to make up its mind (even if it takes a few micro-seconds: the user's time is really more important than a little bit of computation time for the solver).

In this case it is even more troubling. The model has a dummy objective $$z=0$$. So if the solver says, well, may be the model is unbounded, it clearly was not paying attention.

### Constraint solver: Z3/Python

A constraint solver typically has built-in support for the all-different constraint. This allows us to use $$x_{i,j} \in 1,\dots,5$$ directly as main variable. A Z3 Python model can look like:

from z3 import *

grid = [[ 1,   1,   2,   2,   3],
[ 4,   5,   5,   5,   6],
[ 4,   7,   7,   5,   6],
[ 8,   7,   5,   5,  10],
[ 8,   9,   9,  10,  11]]

ROWS = range(len(grid))     # 0..rows-1
COLS = range(len(grid))  # 0..cols-1
MAX = 5                     # 1 <= x[i,j] <= MAX

REGIONS = set([grid[i][j] for i in ROWS for j in COLS])  # region numbers

# for each region form a list of (i,j) tuples
region = [[] for r in REGIONS]
for i in ROWS:
for j in COLS:
r = grid[i][j]
region[r-1] += [(i,j)]

#
# decision variables
#
X = [ [ Int("x%d.%d" % (i+1, j+1)) for j in COLS ] for i in ROWS ]
V = [ Int("v%d" % r) for r in REGIONS]

#
# constraints
#
Bounds = And([And(X[i][j] >= 1,X[i][j] <= MAX) for i in ROWS for j in COLS ])
UniqRows = And([ Distinct([X[i][j] for j in COLS]) for i in ROWS])
UniqCols = And([ Distinct([X[i][j] for i in ROWS]) for j in COLS])
UniqRegionValues = Distinct([V[r-1] for r in REGIONS])
CalcV = And([V[r-1] == Sum([X[i][j] for (i,j) in region[r-1]]) for r in REGIONS])

Cons = [Bounds, UniqRows, UniqCols, UniqRegionValues, CalcV]

#
# solve and print solution
#
s = Solver()
if s.check() == sat:
m = s.model()
sol = [ [ m.evaluate(X[i][j]) for j in COLS] for i in ROWS]
print(sol)

#
# Check for uniqueness
#
Forbid = Or([X[i][j] != sol[i][j] for i in ROWS for j in COLS])
if s.check() == sat:
print("There is an alternative solution.")
else:
print("Solution is unique.")


The output looks like:

[[1, 3, 5, 4, 2], [5, 1, 4, 2, 3], [3, 5, 2, 1, 4], [2, 4, 1, 3, 5], [4, 2, 3, 5, 1]]
Solution is unique.


In this model I tried to be efficient w.r.t constraint Calcv. Basically I tried to make only one pass over the data grid[i][j]. As this model is small, I probably could have dropped the calculation of region and immediate form the equation:

CalcV = And([V[r-1] == Sum([X[i][j] for i in COLS for j in ROWS if grid[i][j] == r]) for r in REGIONS])

This will make multiple passes over grid[i][j], but that is not a problem for a model with this size.

### Constraint solver: Minizinc/Gecode

This model can be coded easily in Minizinc:

include "globals.mzn";

% number of rows, columns, regions
int: m = 5;
int: n = 5;
int: nr = 11;

% sets
set of int: I = 1..m;
set of int: J = 1..n;
set of int: R = 1..nr;

% grid
array[I,J] of int: grid;
grid = [| 1,   1,   2,   2,   3,
| 4,   5,   5,   5,   6,
| 4,   7,   7,   5,   6,
| 8,   7,   5,   5,  10,
| 8,   9,   9,  10,  11  |];

array[I,J] of var int: x;
array[R] of var int: v;

constraint forall (i in I, j in J) ( x[i,j] >= 1 );

constraint forall (i in I, j in J) ( x[i,j] <= 5 );

constraint forall (i in I) ( alldifferent([ x[i,j] | j in J]) );

constraint forall (j in J) ( alldifferent([ x[i,j] | i in I]) );

constraint forall (r in R) ( v[r] = sum( [ x[i,j] | i in I, j in J where grid[i,j] = r] ) );

constraint alldifferent( [v[r] | r in R] );

solve satisfy;


I did not optimize the constraint that calculates v[r]. I assume this version will make multiple passes over grid[i,j].

Gecode will try to find multiple solutions automatically. We don't need to add code for this like we did in the previous examples. Note that there is an option to limit the number of solutions (we should set this higher than one). In this case Gecode finds just one unique solution:

Compiling sudokux.mzn
Running sudokux.mzn
x = array2d(1..5 ,1..5 ,[1, 3, 5, 4, 2, 5, 1, 4, 2, 3, 3, 5, 2, 1, 4, 2, 4, 1, 3, 5, 4, 2, 3, 5, 1]);
v = array1d(1..11 ,[4, 9, 2, 8, 12, 7, 11, 6, 5, 10, 1]);
----------
==========
Finished in 333msec


Minizinc does not know how to print tables so unfortunately we get the results as a long list. It is possible to do your own formatting, but it strange that the default behavior is so user unfriendly. A better layout would really help.

#### Constraint solver: or-tools SAT/CP

Google's or-tools has a number of constraint solvers. Here we use the newer SAT/CP solver. The code is similar to the Z3 code. A little bit more code is needed to handle printing of multiple solution.

from ortools.sat.python import cp_model

#
# used in callback
#
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
"""Print intermediate solutions."""

def __init__(self,X):
self.__solution_count = 0
self.__x = X

def NewSolution(self):
self.__solution_count += 1
R = range(len(self.__x))
C = range(len(self.__x))
for i in R:
for j in C:
print(self.Value(self.__x[i][j]),end=" ")
print()
print("---------------")

def SolutionCount(self):
return self.__solution_count

grid = [[ 1,   1,   2,   2,   3],
[ 4,   5,   5,   5,   6],
[ 4,   7,   7,   5,   6],
[ 8,   7,   5,   5,  10],
[ 8,   9,   9,  10,  11]]

ROWS = range(len(grid))     # 0..rows-1
COLS = range(len(grid))  # 0..cols-1
MAX = 5                     # 1 <= x[i,j] <= MAX

REGIONS = set([grid[i][j] for i in ROWS for j in COLS])  # region numbers, i.e. 1..11

# Create the SAT/CP model.
model = cp_model.CpModel()

#
# decision variables
#
X = [ [ model.NewIntVar(1,MAX,"x%d.%d" % (i+1, j+1)) for j in COLS ] for i in ROWS ]
V = [ model.NewIntVar(1,MAX*MAX,"v%d" % r) for r in REGIONS]

#
# constraints
#
for i in ROWS:

for j in COLS:

for r in REGIONS:
model.Add(V[r-1] == sum([X[i][j] for i in COLS for j in ROWS if grid[i][j] == r ]))

#
# solve
#
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(X)
status = solver.SearchForAllSolutions(model, solution_printer)
print("Number of solutions found: %i" % solution_printer.SolutionCount())
print("Time:  %.2g seconds" % solver.WallTime())


The output looks like:
1 3 5 4 2
5 1 4 2 3
3 5 2 1 4
2 4 1 3 5
4 2 3 5 1
---------------
Number of solutions found: 1
Time:  1 seconds


### Performance

This is not a very easy model. Here are some timings:

ModelSolverSolution TimeNodesIterationsNotes
MIPCBC2761584053192290
ConstraintZ337.5
ConstraintGecode0.333
ConstraintOr-tools1

We see quite a variation in solution times. This is not unusual for combinatorial problems like this.

### References

1. https://blogs.wsj.com/puzzle/2018/05/24/varsity-math-week-141/ and https://momath.org/home/varsity-math/varsity-math-week-141/
2. https://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html
3. Williams, H. Paul and Yan, Hong (2001), Representations of the 'all_different' predicate of constraint satisfaction in integer programming, Informs Journal on Computing, 13 (2). 96-103.
4. Z3, https://github.com/Z3Prover/z3
5. Minizinc, http://www.minizinc.org/
6. Gecode, http://www.gecode.org/