# Numberlink Puzzle

*
*

## Problem

- Connect each two cells with the same numbers by a line.
- The sections of a line run horizontally or vertically.
- Each cell must be visited exactly once by a line.

## Solution

The trick is to identify endpoints of numberlinks. They do have only one neighbor with the same value. Cells wich are part of a path must have two neighbors with the same value.

The challenge is now to encode this insight into a Mixed Integer Programming model and fire up a solver.

There are different approaches to this problem, not all linear. We are going to implement one from Rob Pratt based
on a linearization of a mixed integer quadratically constrained model.
If you are interested in the background of the model (and you should!) then please refer to the inspiration
for this article: Numberlink Models ^{1}.

## Model

$$ k: \text{number k}\\ p,q: \text{cells}\\ $$

### Variables

$$ x_{p,k} = \begin{cases}1 \text{ if cell \(p\) has value \(k\)}\\ 0 \text{ otherwise}\end{cases} $$

### Parameters

$$ N_{pq} = \begin{cases} 1, \ \text{ true if p,q are neighbors }\\ 0, \ \text{ else }\\ \end{cases}\\ c_{pk} = \begin{cases} 1, \ \text{ 1 if cell p is an end-point with value k }\\ 0, \ \text{ else }\\ \end{cases} $$

### Constraints

$$ \sum_{q|N(p,q)} x_{q,k}=1, \ \text{ if cell \(p\) is an end-point with value \(k\)}\\ $$

$$ 2x_{p,k}\le \sum_{q|N_{p,q}} x_{q,k} \le 2x_{p,k}+M_p(1-x_{p,k}), \ \text{ if cell \(p\) is not an end-point} \\ $$

$$ M_p=|N(p,q)|=\sum_{q|N(p,q)} 1, \ \text{ number of neighbors of cell p, at most 4}\\ $$

$$ x_{p,k} = c_{p,k}, \ \text{ if cell p is an end-point}\\ x_{p,k} \in {0,1}\\ $$

## Implementation

We need the set of neighbors, which we use as bespoke index \(N_{pg}\) in order to formulate the model constraints.

```
def create_neighbors(self) -> None:
N = self.dimension
M = self.dimension
for i in range(1, N + 1):
for j in range(1, M + 1):
if i > 1:
self.neighbors[i, j, i - 1, j] = 1
if i < N:
self.neighbors[i, j, i + 1, j] = 1
if j > 1:
self.neighbors[i, j, i, j - 1] = 1
if j < M:
self.neighbors[i, j, i, j + 1] = 1
model.N = Set(initialize=[n for n in self.neighbors if self.neighbors[n] == 1])
```

No we can write the main inequality constraint, i.e. cell \(p\) is not an end-point:

```
def x_lb_c(model, i, j, k):
if (i, j, k) not in model.L:
return 2 * model.x[i, j, k] <= sum(
model.x[iii, jjj, k] for (ii, jj, iii, jjj) in model.N if ii == i and jj == j
)
else:
return Constraint.Skip
model.x_lb_c = Constraint(model.I, model.J, model.K, rule=x_lb_c)
def x_ub_c(model, i, j, k):
if (i, j, k) not in model.L:
return sum(model.x[iii, jjj, k] for (ii, jj, iii, jjj) in model.N if ii == i and jj == j) <= \
2 * model.x[i, j, k] + M * (1 - model.x[i, j, k])
else:
return Constraint.Skip
model.x_ub_c = Constraint(model.I, model.J, model.K, rule=x_ub_c)
```

Additionally, we create a helper index for the given path start- and end-points, which allows us to formulate the corresponding constraint for numberlink edges concise:

```
model.L = Set(initialize=[(i, j, k) for (i, j), k in self.links.items()])
def x_fixed_c(model, i, j, k):
if (i, j, k) in model.L:
return model.x[i, j, k] == 1
return Constraint.Skip
model.x_fixed_c = Constraint(model.I, model.J, model.K, rule=x_fixed_c)
```

## Result

### 11 x 11

Running the model against this 11x11 puzzle ^{2} was easy for CBC on 6 cores on my virtual machine. The solution
took less then 1s.

### 15 x 15

A solution for 15x15 is possible:

- Number of constraints : 5625
- Number of variables : 2700
- Duration: 00:00:40

### 20 x 20

Going for a bigger prize ^{3} however did not result in a solution with CBC. After 12h I stopped the calculation.

Here the difference between Open Source and commercial solver becomes obvious. With a commercial solver (e.g. CPLEX, Gurobi) the problem can be solved within 43min:

- Number of constraints : 20400
- Number of variables : 10000
- Duration: 00:42:47

## Summary

For larger problems CBC MIP solver is not the best tool. For this kind of puzzle SAT solvers might be a better choice ^{1}.