Re-posted from: https://bkamins.github.io/julialang/2022/09/30/knights.html
Introduction
Recently David Amos started posting interesting puzzles on Twitter. One of them
can be found here and asks to find the smallest number of knights that
you need to place on a chess board so that every square is either occupied by a
knight or attacked by a knight.
It looks as an interesting combinatorial problem. However, if one is in a rush
(e.g. I am usually quite busy with fixing bugs DataFrames.jl :)),
one wants to find a solution without having to write a lot of code and
implementing complex algorithms.
We notice that the posed problem is from optimization domain, so using
JuMP.jl to solve it immediately comes to ones mind.
Let us give it a shot!
In the post I use Julia 1.8.1, Cbc.jl 1.0.1, and JuMP.jl v1.3.1.
Solution
In the solution I will use 1 to 8 coordinates for both dimensions of the board
(the original Tweet asked for standard chess coordinates, like e.g. b3 or e5,
but they will not be needed in the solution so I skip this requirement).
First let us load the packages we need:
julia> using Cbc
julia> using JuMP
Next, as required in the Tweet, we define the attacked_by function. It takes
horizontal and vertical coordinate of a square and returns a vector of cells by
which a given cell is attacked via a knight move:
julia> attacked_by(i, j) =
[(p, q) for p in 1:8, q in 1:8
if abs((p - i) * (q - j)) == 2]
attacked_by (generic function with 1 method)
In the definition we use the fact that the only possibility that the absolute
value of the product of two integers is equal to 2 is that one of them must have
absolute value equal to 1 and other absolute value equal to 2.
(If you think that the function does not have an optimal performance I agree
with you, but performance does not matter as the problem we have is small.)
Let us check if it works correctly:
julia> attacked_by(1, 1)
2-element Vector{Tuple{Int64, Int64}}:
(3, 2)
(2, 3)
julia> attacked_by(2, 2)
4-element Vector{Tuple{Int64, Int64}}:
(4, 1)
(4, 3)
(1, 4)
(3, 4)
julia> attacked_by(4, 4)
8-element Vector{Tuple{Int64, Int64}}:
(3, 2)
(5, 2)
(2, 3)
(6, 3)
(2, 5)
(6, 5)
(3, 6)
(5, 6)
It seems all is good.
Now the fun stuff begins. First define our optimization problem:
julia> model = Model(Cbc.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: COIN Branch-and-Cut (Cbc)
julia> @variable(model, x[1:8, 1:8], Bin)
8×8 Matrix{VariableRef}:
x[1,1] x[1,2] x[1,3] x[1,4] x[1,5] x[1,6] x[1,7] x[1,8]
x[2,1] x[2,2] x[2,3] x[2,4] x[2,5] x[2,6] x[2,7] x[2,8]
x[3,1] x[3,2] x[3,3] x[3,4] x[3,5] x[3,6] x[3,7] x[3,8]
x[4,1] x[4,2] x[4,3] x[4,4] x[4,5] x[4,6] x[4,7] x[4,8]
x[5,1] x[5,2] x[5,3] x[5,4] x[5,5] x[5,6] x[5,7] x[5,8]
x[6,1] x[6,2] x[6,3] x[6,4] x[6,5] x[6,6] x[6,7] x[6,8]
x[7,1] x[7,2] x[7,3] x[7,4] x[7,5] x[7,6] x[7,7] x[7,8]
x[8,1] x[8,2] x[8,3] x[8,4] x[8,5] x[8,6] x[8,7] x[8,8]
julia> @objective(model, Min, sum(x))
x[1,1] + x[2,1] + x[3,1] + x[4,1] + x[5,1] + x[6,1] + x[7,1] + x[8,1] +
x[1,2] + x[2,2] + x[3,2] + x[4,2] + x[5,2] + x[6,2] + x[7,2] + x[8,2] +
x[1,3] + x[2,3] + x[3,3] + x[4,3] + x[5,3] + x[6,3] + x[7,3] + x[8,3] +
x[1,4] + x[2,4] + x[3,4] + x[4,4] + x[5,4] + x[6,4] + x[7,4] + x[8,4] +
x[1,5] + x[2,5] + x[3,5] + x[4,5] + x[5,5] + x[6,5] + x[7,5] + x[8,5] +
x[1,6] + x[2,6] + x[3,6] + x[4,6] + x[5,6] + x[6,6] + x[7,6] + x[8,6] +
x[1,7] + x[2,7] + x[3,7] + x[4,7] + x[5,7] + x[6,7] + x[7,7] + x[8,7] +
x[1,8] + x[2,8] + x[3,8] + x[4,8] + x[5,8] + x[6,8] + x[7,8] + x[8,8]
julia> @constraint(model, c[i=1:8, j=1:8],
x[i, j] + sum(x[p...] for p in attacked_by(i, j)) >= 1)
8×8 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape}}:
c[1,1] : x[1,1] + x[3,2] + x[2,3] >= 1.0 … c[1,8] : x[2,6] + x[3,7] + x[1,8] >= 1.0
c[2,1] : x[2,1] + x[4,2] + x[1,3] + x[3,3] >= 1.0 c[2,8] : x[1,6] + x[3,6] + x[4,7] + x[2,8] >= 1.0
c[3,1] : x[3,1] + x[1,2] + x[5,2] + x[2,3] + x[4,3] >= 1.0 c[3,8] : x[2,6] + x[4,6] + x[1,7] + x[5,7] + x[3,8] >= 1.0
c[4,1] : x[4,1] + x[2,2] + x[6,2] + x[3,3] + x[5,3] >= 1.0 c[4,8] : x[3,6] + x[5,6] + x[2,7] + x[6,7] + x[4,8] >= 1.0
c[5,1] : x[5,1] + x[3,2] + x[7,2] + x[4,3] + x[6,3] >= 1.0 c[5,8] : x[4,6] + x[6,6] + x[3,7] + x[7,7] + x[5,8] >= 1.0
c[6,1] : x[6,1] + x[4,2] + x[8,2] + x[5,3] + x[7,3] >= 1.0 … c[6,8] : x[5,6] + x[7,6] + x[4,7] + x[8,7] + x[6,8] >= 1.0
c[7,1] : x[7,1] + x[5,2] + x[6,3] + x[8,3] >= 1.0 c[7,8] : x[6,6] + x[8,6] + x[5,7] + x[7,8] >= 1.0
c[8,1] : x[8,1] + x[6,2] + x[7,3] >= 1.0 c[8,8] : x[7,6] + x[6,7] + x[8,8] >= 1.0
julia> model
A JuMP Model
Minimization problem with:
Variables: 64
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.GreaterThan{Float64}`: 64 constraints
`VariableRef`-in-`MathOptInterface.ZeroOne`: 64 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: COIN Branch-and-Cut (Cbc)
Names registered in the model: c, x
As the variable in the model we set a matrix called x, which will hold 1
in cells where a knight is placed and 0 in cells where no knight is placed.
As an objective we set that we want to minimize the number of used knights.
Finally we add 64 constraints saying that all cells must be attacked or covered
at least once.
We are ready for the final blow:
julia> optimize!(model)
Welcome to the CBC MILP Solver
Version: 2.10.5
Build Date: Jan 1 1970
command line - Cbc_C_Interface -solve -quit (default strategy 1)
Continuous objective value is 12 - 0.00 seconds
Cgl0004I processed model has 64 rows, 64 columns (64 integer (64 of which binary)) and 400 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 12
Cbc0038I Before mini branch and bound, 64 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)
Cbc0038I After 0.02 seconds - Feasibility pump exiting with objective of 12 - took 0.01 seconds
Cbc0012I Integer solution of 12 found by feasibility pump after 0 iterations and 0 nodes (0.03 seconds)
Cbc0001I Search completed - best objective 12, took 0 iterations and 0 nodes (0.03 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 12 to 12
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 12.00000000
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.10
Time (Wallclock seconds): 0.10
Total time (CPU seconds): 0.11 (Wallclock seconds): 0.11
Bang! We are done. We get a lot of output. It looks like the solver found
the optimal solution to the problem that has value 12, so we need at least 12
knights. The report that we see indicates that the problem was relatively easy
to solve.
We know that we need 12 knights, but where to put them? Let us check:
julia> Bool.(value.(x))
8×8 BitMatrix:
0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0
0 1 1 0 1 1 0 0
0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0
0 0 1 1 0 1 1 0
0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0
The optimal solution forms a nice symmetric pattern.
Conclusions
Removing line breaks, that I added to make the code better fit the screen, our
solution has 6 lines (plus 3 lines for lading packages and printing the result).
JuMP.jl is amazingly expressive and easy to use. At the same time notice how
nicely it prints the elements of the optimization model (this is best
appreciated in the terminal in our case as the generated expressions are long).
If you never used it I recommend you check it out.
I can’t wait for next challenges from David Amos.