Tuesday, September 28, 2021

Another Rook Tour on the Chessboard



In [1] the following problem is described:

  1. The chessboard shown above has \(8\times 8 = 64\) cells.
  2. Place numbers 1 through 64 in the cells such that the next number is either directly left, right, below, or above the current number. 
  3. The dark cells must be occupied by prime numbers. Note there are 18 dark cells, and 18 primes in the sequence \(1,\dots,64\). These prime numbers are: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61.

This problem looks very much like the Numbrix puzzle [3]. We don't have the given cells that are provided in Numbrix, but instead, we have some cells that need to be covered by prime numbers.

The path formed by the consecutive values looks like how a rook moves on a chessboard. Which explains the title of this post (and of the original post). If we require a proper tour we need values 1 and 64 to be neighbors. Otherwise, if this is not required, I think a better term is: Hamiltonian path. I'll discuss both cases, starting with the Hamiltonian path problem.

Mathematical Model


As usual for this type of problems, we use the following decision variables: \[\color{darkred}x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has a value $k$} \\ 0 & \text{otherwise}\end{cases}\] With this, our model can look like: 


Hamiltonian Path Model 
\[\begin{align} & \sum_k \color{darkred}x_{i,j,k} = 1 && \forall i,j &&\text{One value in a cell}\\ & \sum_{i,j} \color{darkred}x_{i,j,k} = 1 && \forall k &&\text{Unique values} \\ & \color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} && \forall i,j,k\lt 64 && \text{Next value is a neighbor} \\ & \sum_p \color{darkred}x_{i,j,p} = 1 && \text{for dark cells $(i,j)$} &&\text{Prime values} \\ & \color{darkred}x_{i,j,k} \in \{0,1\} \end{align} \]

Here \(p\) is a subset of \(k\): all the prime values. The next-value constraint is explained in more detail in [3]. This is a fascinating constraint. The only thing we are saying here is: the next value is found among the neighbors of \((i,j)\). Or to be more precise: \[\color{darkred}x_{i,j,k}=1\Rightarrow \color{darkred}x_{i-1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i+1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j-1,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j+1,k+1}=1\]There is no concept of a path here. Rather, the formation of a path is a natural consequence. I assume in this constraint that when we address outside the board, the corresponding variable is zero. I think this model is quite elegant. It precisely states the conditions we must obey for a valid solution. 

Notes:

  • This model has no objective.
  • We don't require a tour here. There is no restriction that says that values 1 and 64 must be neighbors. That version of the problem is looked at further down.
  • This is not a very small model. The number of \(\color{darkred}x_{i,j,k}\) variables is \(8 \times 8 \times 64 = 4096\). The model size (after adding a dummy objective) is:

MODEL STATISTICS

BLOCKS OF EQUATIONS           6     SINGLE EQUATIONS        4,179
BLOCKS OF VARIABLES           2     SINGLE VARIABLES        4,097
NON ZERO ELEMENTS        26,661     DISCRETE VARIABLES      4,096


A final question is: is the solution unique or are there multiple solutions? To forbid a previously found solution, we can add the cut: \[\sum_{i,j,k} \color{darkblue}x^*_{i,j,k} \cdot \color{darkred}x_{i,j,k} \le  {\bf{card}}(k)-1\] where \(\color{darkblue}x^*_{i,j,k}\) is such a previously found solution. This constraint says: at least one cell should change. Note that we never can only change a 0 value to a 1: the total number of ones is given by \({\bf{card}}(k)=64\). So we only need to keep track of the ones. And that is what this constraint is doing: force to flip at least one 1 to a 0. We repeat solving and adding a cut to the model until the model becomes infeasible. An alternative approach would be to use the solution pool to multiple optimal solutions. The solution pool is a technology found in some advanced solvers. As the number of solutions is small, I use here the DIY approach using the "no good" cuts.

The feasible paths we find are:


----    106 PARAMETER result  negative numbers indicate prime cells

INDEX 1 = solution1

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          51          50         -47          -2          -5           8           9          16
i6          52          49          48           1           6          -7          10          15
i7         -53          56          57          58         -59          60         -11          14
i8          54          55          64          63          62         -61          12         -13

INDEX 1 = solution2

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          51          50         -47          -2          -5           8           9          16
i6          52          49          48           1           6          -7          10          15
i7         -53          64          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13

INDEX 1 = solution3

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          49          48         -47          -2          -5           8           9          16
i6          50          51          64           1           6          -7          10          15
i7         -53          52          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13

INDEX 1 = solution4

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -43          44          45           6          -7           8           9          10
i2          42         -47          46          -5          64          25          24         -11
i3         -41          48          49           4          63          26         -23          12
i4          40          51          50          -3          62          27          22         -13
i5          39          52         -53          -2         -61          28          21          14
i6          38          55          54           1          60         -29          20          15
i7         -37          56          57          58         -59          30         -19          16
i8          36          35          34          33          32         -31          18         -17


Or plotted using R/ggplot:


All feasible paths



Finding feasible tours


Looking at these solutions, we see that just one solution forms a proper tour. The obvious question is: can we add a constraint to the model that just delivers this tour. Actually, this is rather simple. We need to change the equation next. The mathematical formulation looks like: \[\begin{cases}\color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} & \forall i,j,k\lt 64 \\ \color{darkred}x_{i-1,j,1}+\color{darkred}x_{i+1,j,1}+\color{darkred}x_{i,j-1,1}+\color{darkred}x_{i,j+1,1} \ge \color{darkred}x_{i,j,64} & \forall i,j\end{cases}\]

Here is how we can do this a bit more easily in GAMS. The original version had (see the appendix below for the complete model):

next(i,j,k+1)..   
   x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k);

This version has k+1 as equation index. That means we generate this equation for all but the last k. If we want to require that value 1 and 64 are neighbors, we can simply write:

next(i,j,k)..   
   x(i-1,j,k++1)+x(i+1,j,k++1)+x(i,j-1,k++1)+x(i,j+1,k++1) =g= x(i,j,k);


Here the ++ operator is used to implement a circular lead. Also, note that the equation next is now indexed by all k. The net effect is that we add 64 rows to the model (one for each cell). If we run the model with this equation we get a single solution:


----    106 PARAMETER result  negative numbers indicate prime cells

INDEX 1 = solution1

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          49          48         -47          -2          -5           8           9          16
i6          50          51          64           1           6          -7          10          15
i7         -53          52          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13


This is the third path from our earlier picture.

Conclusions


The problem of finding a Hamiltonian path while obeying the prime-value-only cells, is a fairly close variant of the Numbrix model discussed in [3]. Finding a tour can be handled by changing a single equation (but in an interesting way).

This was an intriguing modeling exercise. Both the underlying mathematical model and the implementation are not very complicated but somewhat surprising. Some puzzles can be elegantly stated as mathematical models and simply solved with off-the-shelf solvers. Probably this was not intended, but I often find doing these modeling exercises just as rewarding as solving these puzzles by hand. It forces you to think about the problem in a different way, totally different from solving it by hand.

References



Appendix: GAMS model

$ontext

  
A Numbrix-like puzzle.

  
Place the numbers 1..64 on an 8x8 board such that:

     
1. next values are neighbors. A neighbor is defined as
        
below, left, right, or above.
     
2. certain cells need to contain a prime number

  
https://puzzling.stackexchange.com/questions/111818/another-rooks-tour-of-the-chessboard

$offtext

*----------------------------------------------------------------------
* data
*----------------------------------------------------------------------

sets
  i
'rows'    /i1*i8/
  j
'columns' /j1*j8/
  k
'values'  /1*64/
  kprime(k)
'primes' /
            
2,3,5,7,11,13,17,19,23,29,
            
31,37,41,43,47,53,59,61
   
/
;

acronym p;
table board(i,j) 'locations with a p should get a prime value'
       
j1 j2 j3 j4 j5 j6 j7 j8
 
i1     p           p
 
i2        p     p           p
 
i3     p                 p
 
i4              p           p
 
i5           p  p  p
 
i6                    p
 
i7     p           p     p
 
i8                    p     p
;

*----------------------------------------------------------------------
* model
*----------------------------------------------------------------------

set sol /solution1*solution10/;
set solfound(sol) 'solutions found';
parameter solution(sol,i,j,k) 'solutions found so far';
solfound(sol) =
no;
solution(sol,i,j,k) = 0;


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

equations
   onevalue(i,j)
'each cell has one value (between 1 and 81)'
   unique(k)    
'each value must be in exactly one cell'
   next(i,j,k)  
'x(i,j,k)=1 => a neighbor is k+1'
   isprime(i,j) 
'we need a prime on (i,j)'
   cut(sol)     
'cut off previous solutions'
   obj          
'dummy'
;

obj.. dummy =e= 0;

onevalue(i,j)..
sum(k, x(i,j,k)) =e= 1;
unique(k)..    
sum((i,j),x(i,j,k)) =e= 1;

next(i,j,k+1)..
    x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k);

isprime(i,j)$(board(i,j)=p)..
   
sum(kprime, x(i,j,kprime)) =e= 1;

cut(solfound)..
   
sum((i,j,k),solution(solfound,i,j,k)*x(i,j,k)) =l= card(k)-1;

model m /all/;
option mip=cplex,threads=8;

*----------------------------------------------------------------------
* solve loop
*----------------------------------------------------------------------


loop(sol,
   
solve m minimizing dummy using mip;
   
break$(m.modelstat <> %modelStat.optimal%
          
and m.modelstat <> %modelStat.integerSolution%);
    solfound(sol) =
yes;
    solution(sol,i,j,k) = round(x.l(i,j,k));
   
display x.l;
);


*----------------------------------------------------------------------
* reporting
*----------------------------------------------------------------------

parameter result(sol,i,j) 'negative numbers indicate prime cells';
result(sol,i,j) =
sum(k, k.val*solution(sol,i,j,k));
result(sol,i,j)$(board(i,j)=p) = -result(sol,i,j);
option result:0:1:1;
display result;


Notes:
  1. An acronym allows a string to be used in a numerical table (like inf, eps etc.).
  2. The equation next is indexed by \(k+1\). This means the equation is generated for all but the last k.
  3. If \(i-1\),\(i+1\),\(j-1\),\(j+1\) is outside the domain (i.e. outside the board), GAMS will return zero for the variable (or in other words: ignore that entry). That is exactly what we need here. So we don't need to add additional protection to remain inside the board.
  4. This version of the model generates all feasible Hamiltonian paths.
  5. We could have reduced the number of variables a bit by only generating valid \(\color{darkred}x_{i,j,k}\): for prime cells, non-prime \(k\)'s are not allowed. That would also eliminate the isprime equation.

4 comments:

  1. https://gist.github.com/zayenz/46fb486a4454aac72127376436b823ef is a MiniZinc model of the proper tour version of the problem that uses the global circuit constraint to model the tour-part.

    Here is the output from finding all the solutions using MiniZinc 2.5.4 and Gecode 6.3.0 as the solver
    ```
    Running prime-position-tour.mzn
    Board:
    [| 41, 40, 39, 30, 29, 28, 21, 20 |
    42, 37, 38, 31, 32, 27, 22, 19 |
    43, 36, 35, 34, 33, 26, 23, 18 |
    44, 45, 46, 3, 4, 25, 24, 17 |
    49, 48, 47, 2, 5, 8, 9, 16 |
    50, 51, 64, 1, 6, 7, 10, 15 |
    53, 52, 63, 62, 61, 60, 11, 14 |
    54, 55, 56, 57, 58, 59, 12, 13 |]
    Next:
    [| 9, 1, 2, 12, 4, 5, 15, 7 |
    17, 11, 3, 13, 21, 6, 23, 8 |
    25, 10, 18, 19, 20, 14, 31, 16 |
    26, 27, 35, 29, 37, 22, 30, 24 |
    41, 33, 34, 28, 45, 39, 47, 32 |
    42, 50, 44, 36, 46, 38, 55, 40 |
    57, 49, 43, 51, 52, 53, 63, 48 |
    58, 59, 60, 61, 62, 54, 64, 56 |]
    ----------
    ==========
    %%%mzn-stat: initTime=0.003434
    %%%mzn-stat: solveTime=0.29611
    %%%mzn-stat: solutions=1
    %%%mzn-stat: variables=256
    %%%mzn-stat: propagators=236
    %%%mzn-stat: propagations=1076445
    %%%mzn-stat: nodes=3713
    %%%mzn-stat: failures=1856
    %%%mzn-stat: restarts=0
    %%%mzn-stat: peakDepth=22
    %%%mzn-stat-end
    Finished in 656msec
    ```

    ReplyDelete
    Replies
    1. That is very fast. Of course the model is smaller (say in terms of number of variables) and you are using a special purpose constraint (circuit). It demonstrates that a good CP model is quite different than a MIP model.

      Delete
    2. How long did it take for your model to solve the problem? This is probably a case that really benefits from having the right global constraint available.

      The model uses a mod-constraint, which is not the best in terms of propagation. It might be beneficial to replace it with an extensional constraint using a table, just listing the valid pairs of values ([<1,2>,<2,3>,...,<64,1>]) instead. It was fast enough that it didn't matter here.

      Note that it is very important to run the search on the next-variables in the model. If the board-variables (with the numbers for each position) are used for searhc, the problem is much harder. I would guess that an even faster solution would be to use Warnsdorffs heuristic for the search (follow the current path in the next-variables, and go to the next possibility with the fewest remaining possibilities), since that works well for Kinghts tour problems. For that one would have to implement it in a solver like Gecode. Might not work for this case, but could be interesting to try some time.

      Delete
    3. Cplex took 6 seconds to find the tour (default settings on a slow laptop).

      Delete