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