sysid blog

Numberlink Puzzle

May 30, 2020 ☕️ 4 min read



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



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.


k:number kp,q:cellsk: \text{number k}\\ p,q: \text{cells}\\


xp,k={1 if cell p has value k0 otherwisex_{p,k} = \begin{cases}1 \text{ if cell $p$ has value $k$}\\ 0 \text{ otherwise}\end{cases}


Npq={1,  true if p,q are neighbors 0,  else cpk={1,  1 if cell p is an end-point with value k 0,  else 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}


qN(p,q)xq,k=1,  if cell p is an end-point with value k\sum_{q|N(p,q)} x_{q,k}=1, \ \text{ if cell $p$ is an end-point with value $k$}\\ 2xp,kqNp,qxq,k2xp,k+Mp(1xp,k),  if cell p is not an end-point2x_{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} \\Mp=N(p,q)=qN(p,q)1,  number of neighbors of cell p, at most 4M_p=|N(p,q)|=\sum_{q|N(p,q)} 1, \ \text{ number of neighbors of cell p, at most 4}\\xp,k=cp,k,  if cell p is an end-pointxp,k{0,1}x_{p,k} = c_{p,k}, \ \text{ if cell p is an end-point}\\ x_{p,k} \in \{0,1\}\\


We need the set of neighbors, which we use as bespoke index NpgN_{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 pp 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
        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])
        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)


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



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

always be building, by sysid.

© 2022