There are two major implications for the model itself:

- 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\}\).
- 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.

seti '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 rowsn = len(blockno[ 0]) # number of
columnsR = 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);setsbmap(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'; parameterslen(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 datamaxb = [ 0]*bmax # number of cells in each blockbij = [[] for b in B] # list of cells in each blocknij = [ [ [] for j in C ] for i in R ] # neighborsfor i in R:for j in C:bn = blockno[i][j]- 1maxb[bn] += 1bij[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 valuesx.fx(i,j,k)$(initval(i,j)= ord(k)) = 1;equationoneval(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

# ModelX = [ [ 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 |

##### Solving and printing results: Python/Z3

# solve and
prints = 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], |

##### 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
uniquenessCut = 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.") |