Tuesday, November 7, 2017

Suguru: GAMS/MIP vs Python/Z3

 For many logical puzzles a constraint solver like Z3 may be more suitable than a Mathematical Programming approach using a MIP solver. In [1] a MIP formulation was presented for Suguru puzzle. Now let’s see how much difference we observe when we use Z3.

There are two major implications for the model itself:

  1. We can use integers directly, i.e.  \(x_{i,j} \in \{1,\dots,u_{i,j}\}\) instead of \(x_{i,j,k} \in \{0,1\}\).
  2. We have access to constructs like Distinct and != (unequal).

However most of the code is devoted to data manipulation as we will see below. Although the approaches are different (different languages, different solver technologies), the actual code is remarkably similar.

Data entry: GAMS

First we need to get the data into GAMS. Unfortunately we need to enter two matrices: one for the blocks and one for the initial values.

set
  i
'rows' /i1*i7/
  j
'columns' /j1*j9/
  k
'values' /1*9/
  b
'blocks' /b1*b15/
;
table blockno(i,j)
     
j1 j2 j3 j4 j5 j6 j7 j8 j9
 
i1   1  1  2  2  2  3  3  4  6
 
i2   1  1  2  2  3  3  4  4  6
 
i3   1  7  8  8  3  4  4  5  6
 
i4   7  7  7  8  8  9  9  6  6
 
i5   7 10 10  8 12  9  9 14 15
 
i6  11 11 10 10 12 12  9 15 15
 
i7  11 11 10 12 12 13 13 15 15
;
table initval(i,j)
     
j1 j2 j3 j4 j5 j6 j7 j8 j9
 
i1   1        1
 
i2      2                 3
 
i3   4     4              1  4
 
i4            1
 
i5   4     5     4     5  1
 
i6                        3
 
i7               2           5
;

The set k is a bit of an over-estimate: we actually need only values 1 through 5 for this problem (the maximum size of a block is 5).

Data entry: Python

Entering the same boards in Python can look like:

from z3 import *

blockno = [ [
1, 1, 2, 2, 2, 3, 3, 4, 6],
             [
1, 1, 2, 2, 3, 3, 4, 4, 6],
             [
1, 7, 8, 8, 3, 4, 4, 5, 6],
             [
7, 7, 7, 8, 8, 9, 9, 6, 6],
             [
7,10,10, 8,12, 9, 9,14,15],
             [
11,11,10,10,12,12, 9,15,15],
             [
11,11,10,12,12,13,13,15,15] ]

initval = [ [
1, 0, 0, 1, 0, 0, 0, 0, 0],
             [
0, 2, 0, 0, 0, 0, 0, 3, 0],
             [
4, 0, 4, 0, 0, 0, 0, 1, 4],
             [
0, 0, 0, 1, 0, 0, 0, 0, 0],
             [
4, 0, 5, 0, 4, 0, 5, 1, 0],
             [
0, 0, 0, 0, 0, 0, 0, 3, 0],
             [
0, 0, 0, 0, 2, 0, 0, 0, 5] ]

m = len(blockno) 
# number of rows
n = len(blockno[
0]) # number of columns
R = range(m)
C = range(n)
bmax = max([blockno[i][j]
for i in R for j in C])
B = range(bmax)

The order is a bit reversed here: first enter the boards with the block numbers and the initial values and then calculate the dimensions.

Derived data: GAMS

We need to calculate a number of things before we can actually write the equations.

alias (i,i2),(j,j2);

sets
    bmap(b,i,j)
'maps between b <-> (i,j)'
    bk(b,k)
'allowed (b,k) combinations'
    ijk(i,j,k)
'allowed values'
    n(i,j,i2,j2)
'neighbors'
;
parameters
    len(b)
'number of cells in block'
;

bmap(b,i,j) = blockno(i,j)=
ord(b);
len(b) =
sum(bmap(b,i,j),1);
bk(b,k) =
ord(k)<=len(b);
ijk(i,j,k) =
sum((bmap(b,i,j),bk(b,k)),1);

n(i,j,i,j+1)   = blockno(i,j)<>blockno(i,j+1);
n(i,j,i+1,j-1) = blockno(i,j)<>blockno(i+1,j-1);
n(i,j,i+1,j)   = blockno(i,j)<>blockno(i+1,j);
n(i,j,i+1,j+1) = blockno(i,j)<>blockno(i+1,j+1);

Note that the set n (indicating neighbors we need to inspect for different values) only looks “forward” from cell \((i,j)\). This is to prevent double-counting: we don’t want to compare neighboring cells twice in the equations.

Derived data: Python

The Python/Z3 model needs to calculate similar derived data:

# derived data
maxb = [
0]*bmax    # number of cells in each block
bij = [[]
for b in B]  # list of cells in each block
nij = [ [ []
for j in C ] for i in R ] # neighbors
for i in R:
    
for j in C:
         bn = blockno[i][j]-
1
         maxb[bn] +=
1
         bij[bn] += [(i,j)]
        
if j+1and blockno[i][j]!=blockno[i][j+1]:
             nij[i][j] += [(i,j+
1)]
        
if i+1and j-1>=0 and blockno[i][j]!=blockno[i+1][j-1]:
             nij[i][j] += [(i+
1,j-1)]
        
if i+1and blockno[i][j]!=blockno[i+1][j]:
             nij[i][j] += [(i+
1,j)]
        
if i+1and j+1and blockno[i][j]!=blockno[i+1][j+1]:
             nij[i][j] += [(i+
1,j+1)]

# upperbound on X[i,j]
maxij = [ [ maxb[blockno[i][j]-
1] for j in C ] for i in R ]

I use here nested lists. bij has a list of cells for each block. nij contains a list of neighboring cells for each cell.

Model constraints: GAMS

The decision variables and linear constraints are described in [1]. The GAMS representation is very close the mathematical formulation.

variable dummy 'objective variable';
binary variable x(i,j,k);

* fix initial values
x.fx(i,j,k)$(initval(i,j)=
ord(k)) = 1;


equation
    oneval(i,j)  
'pick one value per cell'
    unique(b,k)  
'unique in a block'
    neighbors(i,j,i2,j2,k)
'neighboring cells must be different'
    edummy       
'dummy objective'
;

oneval(i,j)..
sum(ijk(i,j,k), x(i,j,k)) =e= 1;

unique(bk(b,k))..
sum(bmap(b,i,j), x(i,j,k)) =e= 1;

neighbors(i,j,i2,j2,k)$(n(i,j,i2,j2)
and ijk(i,j,k) and ijk(i2,j2,k))..
    x(i,j,k)+x(i2,j2,k) =l= 1;

edummy.. dummy =e= 0;

Notice we added a dummy objective to the model.

Model constraints: Z3
# Model
X = [ [ Int(
"x%d.%d" % (i+1, j+1)) for j in C ] for i in R ]

Xlo =
And([X[i][j] >= 1 for i in R for j in C])
Xup =
And([X[i][j] <= maxij[i][j] for i in R for j in C])

Uniq =
And([ Distinct([X[i][j] for (i,j) in bij[b]]) for b in B])
Nbor =
And([ And([X[i][j] != X[i2][j2] for (i2,j2) in nij[i][j] ]) for i in R for j in C if len(nij[i][j])>0 ])

Xfix =
And([X[i][j] == initval[i][j] for i in R for j in C if initval[i][j]>0])

Cons = [Xlo,Xup,Uniq,Nbor,Xfix]

Most of this is boilerplate. We use Distinct to make the contents of the cells unique within a block. The != operator is used to make sure neighbors of different blocks do not have the same value. Both these constructs are not so easy to implement in a MIP model, but there we were rescued by using binary variables instead.

Solving and printing results: GAMS

We only look for a feasible solution (any feasible integer solution will be globally optimal immediately). This means we don’t have to worry about setting the OPTCR option in GAMS. The default gap of 10% is irrelevant here as the solver will stop after finding the first integer solution.

model m /all/;
solve m minimizing dummy using mip;

parameter v(i,j) 'value of each cell';
v(i,j) =
sum(ijk(i,j,k), x.l(i,j,k)*ord(k));
option v:0;
display v;

To report a solution we recover the integer cell values. The results look like:

----     81 PARAMETER v  value of each cell

            j1          j2          j3          j4          j5          j6          j7          j8          j9

i1           1           5           4           1           5           1           5           2           1
i2           3           2           3           2           3           2           4           3           5
i3           4           1           4           5           4           1           5           1           4
i4           5           2           3           1           3           2           3           2           3
i5           4           1           5           2           4           1           5           1           4
i6           2           3           4           3           5           3           4           3           2
i7           4           1           2           1           2           1           2           1           5

Solving and printing results: Python/Z3
# solve and print
s = Solver()
s.add(Cons)
if s.check() == sat:
     m = s.model()
     sol = [ [ m.evaluate(X[i][j])
for j in C] for i in R]
    
print(sol)

The solution looks like:

[[1, 5, 4, 1, 5, 1, 5, 2, 1], [3, 2, 3, 2, 3, 2, 4, 3, 5], [4, 1, 4, 5, 4, 1, 5, 1, 4], 
[5, 2, 3, 1, 3, 2, 3, 2, 3], [4, 1, 5, 2, 4, 1, 5, 1, 4], [2, 3, 4, 3, 5, 3, 4, 3, 2],
[4, 1, 2, 1, 2, 1, 2, 1, 5]]
Check for uniqueness: GAMS

To check if the solution is unique we add a cut and resolve. If this is infeasible we know there was only one solution.


*
* test if solution is unique
*

equation cut 'forbid current solution';
cut..
sum(ijk(i,j,k), x.l(i,j,k)*x(i,j,k)) =l= card(i)*card(j)-1;

model m2 /all/;
solve m2 minimizing dummy using mip;
abort$(m2.modelstat=1 or m2.modelstat=8) "Unexpected solution found",x.l;

We abort if the model status is not "optimal" (1) or "integer solution found" (8). (It should be enough to check on optimal only, we err on the safe side here).

Check for uniqueness: Z3

# Check for uniqueness
Cut =
Or([X[i][j] != sol[i][j] for i in R for j in C])
s.add(Cut)
if s.check() == sat:
    
print("There is an alternative solution.")
else:
    
print("Solution is unique.")
References
  1. Suguru, http://yetanothermathprogrammingconsultant.blogspot.com/2017/11/suguru.html