We shall solve Sudoku puzzles using two approaches:
In [1]:
# Read the sudoku
using CSV
In [2]:
s=CSV.read("files/sudoku.csv",header=0)
Out[2]:
In [3]:
s=convert(Matrix,s)
Out[3]:
In [4]:
# Function to print sudoku from sudoku.jl
function print_sudoku(sudoku::Matrix)
# println("Sudoku:")
println("[-----------------------]")
for row in 1:9
print("[ ")
for col in 1:9
print(sudoku[row, col], " ")
if col % 3 == 0 && col < 9
print("| ")
end
end
println("]")
if row % 3 == 0
println("[-----------------------]")
end
end
end
Out[4]:
In [5]:
print_sudoku(s)
In [6]:
using JuMP
using GLPK
We modify the function example_sudoku() from the file sudoku.jl so that we can use it directly on the array.
In [7]:
function solve_sudoku(sudoku::Matrix)
initial_grid=sudoku
model = Model(with_optimizer(GLPK.Optimizer)) #GLPK
@variable(model, x[1:9, 1:9, 1:9], Bin)
@constraints(model, begin
# Constraint 1 - Only one value appears in each cell
cell[i in 1:9, j in 1:9], sum(x[i, j, :]) == 1
# Constraint 2 - Each value appears in each row once only
row[i in 1:9, k in 1:9], sum(x[i, :, k]) == 1
# Constraint 3 - Each value appears in each column once only
col[j in 1:9, k in 1:9], sum(x[:, j, k]) == 1
# Constraint 4 - Each value appears in each 3x3 subgrid once only
subgrid[i=1:3:7, j=1:3:7, val=1:9], sum(x[i:i + 2, j:j + 2, val]) == 1
end)
# Initial solution
for row in 1:9, col in 1:9
if initial_grid[row, col] != 0
@constraint(model, x[row, col, initial_grid[row, col]] == 1)
end
end
# Solve it
JuMP.optimize!(model)
term_status = JuMP.termination_status(model)
primal_status = JuMP.primal_status(model)
is_optimal = term_status == MOI.OPTIMAL
# Check solution
if is_optimal
mip_solution = JuMP.value.(x)
sol = zeros(Int, 9, 9)
for row in 1:9, col in 1:9, val in 1:9
if mip_solution[row, col, val] >= 0.9
sol[row, col] = val
end
end
return sol
else
error("The solver did not find an optimal solution.")
end
end
Out[7]:
In [8]:
# Solve the sudoku with the method from the file
# MIP (Mixed Integer Programing) from GLPK
sol=solve_sudoku(s)
Out[8]:
In [9]:
# Nicer print
print_sudoku(sol)
The formulation from P. Babu, K. Pelckmans, P. Stoica and J. Li, Linear Systems, Sparse Solutions, and Sudoku is essentially the same. We first ceate the constraints matrix for rows, columns, boxes and cells (all except clues).
In [10]:
# Convenience function
unblock(A) = mapreduce(identity, hcat, [mapreduce(identity, vcat, A[:,i])
for i = 1:size(A,2)])
using Plots
using LinearAlgebra
In [11]:
n=9
n2=81
n3=729
# Row constraints
I9=Array{Int}(I,9,9)
Rows=[repeat(I9,1,9) zeros(Int,n,n3-n2)]
for i=1:8
Rows=vcat(Rows,[zeros(Int,n,i*n2) repeat(I9,1,9) zeros(Int,n,(n-i-1)*n2)])
end
# Column constraints
Cols=repeat([I9 zeros(Int,9,72)],1,9)
for i=1:8
Cols=vcat(Cols,[zeros(Int,9,9*i) repeat([I9 zeros(Int,9,72)],1,8) I9 zeros(Int,9,72-9*i)])
end
# Box constraints
z9=zeros(Int,9,9)
Box=[z9 for i=1:n, j=1:n2]
for i=1:3, j=0:2
Box[i,9*j+1+3*(i-1):9*j+3+3*(i-1)]=[I9 for k=1:3]
end
for i=1:3, j=0:2
Box[i+3,27+9*j+1+3*(i-1):27+9*j+3+3*(i-1)]=[I9 for k=1:3]
end
for i=1:3, j=0:2
Box[i+6,54+9*j+1+3*(i-1):54+9*j+3+3*(i-1)]=[I9 for k=1:3]
end
Box=unblock(Box)
# Cell constraints
Cell=zeros(Int,n2,n3)
for i=1:n2
Cell[i,9*(i-1)+1:9*i]=ones(Int,1,9)
end
C=vcat(Rows,Cols,Box,Cell)
Out[11]:
In [12]:
using SparseArrays
spy(sparse(C),color=:pu_or,clim=(0.0,1))
Out[12]:
In [13]:
using MathProgBase
using GLPKMathProgInterface
In [14]:
function solve_sudoku(sudoku::Matrix, C::Matrix)
# Clue constraints
clue=count(!iszero,sudoku)
Clue=zeros(Int,clue,n3)
ind=findall(!iszero,sudoku'[:])
for i=1:clue
Clue[i,9*(ind[i]-1)+sudoku'[ind[i]]]=1
end
C=vcat(C,Clue)
# Solve it
l1=linprog(ones(n3),C,'=',ones(size(C,1)),0,Inf,GLPKSolverLP())
# Reconstruct the solution
if l1.status == :Infeasible
error("No solution found!")
else
return Matrix(reshape(map(Int,round.(collect(1:9)'*reshape(l1.sol,9,81))),9,9)'), l1.status
end
end
Out[14]:
In [15]:
@time sol₁,status=solve_sudoku(s,C)
Out[15]:
In [16]:
print_sudoku(sol₁)
In [17]:
sol-sol₁
Out[17]:
Let us try two harder puzzles.
In [18]:
s=CSV.read("files/sudokuhard.csv",header=0)
s=convert(Matrix,s)
print_sudoku(s)
In [19]:
@time sol=solve_sudoku(s)
print_sudoku(sol)
In [20]:
@time sol₁,status=solve_sudoku(s,C)
print_sudoku(sol₁), status
Out[20]:
In [21]:
# New puzzle
s=[ 0 0 0 6 0 1 4 0 0;
9 0 0 2 4 0 7 0 0;
0 0 3 0 0 0 0 0 9;
0 6 0 0 0 5 0 0 0;
3 1 0 0 0 0 0 0 4;
0 0 0 7 0 0 0 6 0;
2 0 0 0 0 0 9 0 0;
4 0 0 0 6 2 0 0 1;
0 0 1 9 0 8 0 0 0]
print_sudoku(s)
In [22]:
# Mixed integer programming is fine,
@time sol=solve_sudoku(s)
print_sudoku(sol)
In [23]:
# But!
@time sol₁,status=solve_sudoku(s,C)
print_sudoku(sol₁), status
Out[23]:
We see that the $l_1$ sudoku solution is not always possible.
In [ ]: