This program computes the longest path that doesn't intersect itself.
by: Laurent Lessard
using PyPlot
using PyCall
@pyimport matplotlib.patches as patch
# declare some global variables (these will be changed later)
global N1,N2,N,M,A,B,C,hx,hy
# returns 1 if (x,y) is a valid cell, and 0 otherwise
function iscell(x,y)
1 <= x && x <= N1 && 1 <= y && y <= N2
end
# linear index corresponding to the cell with coordinates (p,q)
function coord2cell(p,q)
sub2ind((N1,N2),p,q)
end
# coorindates of the cell with linear index c
function cell2coord(c)
ind2sub((N1,N2),c)
end
;
# convert a move to a pair of coordinates (uses incidence matrix)
function move2coords(m)
cells = find(B[:,m])
p,q = cell2coord(cells[1])
r,s = cell2coord(cells[2])
p,q,r,s
end
# generate adjacency, incidence, and conflict matrix
function generateABC()
global N1,N2,N,M,A,B,C,hx,hy
# create an adjacency matrix A for admissible knight moves (cells x cells)
# A[c1,c2] is equal to 1 if a knight can move from cell c1 to cell c2
N = N1*N2
A = spzeros(N,N)
for i = 1:N1
for j = 1:N2 # loop over (i,j) starting coordinate
c1 = coord2cell(i,j)
for moves in [(hx,hy),(hy,hx),(-hx,hy),(-hy,hx),(-hx,-hy),(-hy,-hx),(hx,-hy),(hy,-hx)]
p,q = i+moves[1],j+moves[2]
if iscell(p,q)
c2 = coord2cell(p,q)
if c1 <= c2
A[c1,c2] = 1
end
end
end
end
end
# create an incidence matrix B for admissible knights moves (moves x cells)
# B[c,m] is equal to 1 if move m involves cell c. (exactly two ones per column)
moves = find(A)
M = length(moves)
B = spzeros(N,M)
for (i,m) in enumerate(moves)
p,q = ind2sub((N,N),m)
B[p,i] = 1
B[q,i] = 1
end
# create a conflict matrix C for non-intersecting paths (moves x moves)
# C[m1,m2] is equal to 1 if moves m1 and m2 intersect.
# this matrix does not include moves that share the same start/end cells
C = spzeros(M,M)
for m = 1:M
p,q,r,s = move2coords(m)
c1,c2 = coord2cell(p,q),coord2cell(r,s) # cell identifiers
mx,my = (p+r)/2,(q+s)/2 # midpoints
s1 = (q-s)/(p-r) # slope
for mm = 1:M
u,v,x,y = move2coords(mm)
d1,d2 = coord2cell(u,v),coord2cell(x,y) # cell identifiers
mmx,mmy = (u+x)/2,(v+y)/2 # midpoints
s2 = (v-y)/(u-x) # slope
# slope must be different
if s1 != s2
α = [p-r x-u; q-s y-v]\[p-u; q-v]
if (0 ≤ α[1] ≤ 1) & (0 ≤ α[2] ≤ 1) # they intersect!
# shouldn't share any cells in common
if (c1 != d1) & (c1 != d2) & (c2 != d1) & (c2 != d2)
C[m,mm] = 1
end
end
end
end
end
;
end
function initializeEverything(n1,n2,h1,h2)
global N1,N2,N,M,A,B,C,hx,hy
N1,N2 = n1,n2
hx,hy = h1,h2
generateABC()
end
;
# HELPER FUNCTION: identify all cycles in the solution
function getCycles(xsol)
used_edges = find(xsol)
Bsol = B[:,used_edges]
# identify the open path (won't have any cycles)
start,ending = find(getvalue(w))
open_path = []
first_node = start
first_edge = find(Bsol[first_node,:])[1]
push!(open_path,used_edges[first_edge])
while true
next_node = setdiff(find(Bsol[:,first_edge]),first_node)[1]
next_edges = find(Bsol[next_node,:])
if length(next_edges) == 1
break
else
next_edge = setdiff(next_edges,first_edge)[1]
push!(open_path,used_edges[next_edge])
first_node,first_edge = next_node,next_edge
end
end
# remove open path, identify cycles
used_edges = setdiff(used_edges,open_path)
cycles = []
while !isempty(used_edges)
Bsol = B[:,used_edges]
start,ending = find(B[:,used_edges[1]])
cycle = []
first_node = start
initial_edge = find(Bsol[first_node,:])[1]
first_edge = initial_edge
push!(cycle,used_edges[first_edge])
while true
next_node = setdiff(find(Bsol[:,first_edge]),first_node)[1]
next_edge = setdiff(find(Bsol[next_node,:]),first_edge)[1]
push!(cycle,used_edges[next_edge])
if next_edge == initial_edge
break
else
first_node,first_edge = next_node,next_edge
end
end
push!(cycles,cycle)
used_edges = setdiff(used_edges,cycle)
end
cycles
end
;
# HELPER FUNCTION: generate a list of coordinates from a solution
function sol2coordlist(xsol)
used_edges = find(xsol)
Bsol = B[:,used_edges]
start,ending = find(getvalue(w))
open_path = []
first_node = start
first_edge = find(Bsol[first_node,:])[1]
push!(open_path,cell2coord(first_node))
while true
next_node = setdiff(find(Bsol[:,first_edge]),first_node)[1]
next_edges = find(Bsol[next_node,:])
if length(next_edges) == 1
break
else
next_edge = setdiff(next_edges,first_edge)[1]
push!(open_path,cell2coord(next_node))
first_node,first_edge = next_node,next_edge
end
end
open_path
end
;
# function that draws the board
function drawboard(cfig=[],a=1,b=1,c=1)
if cfig == []
cfig = figure(figsize=(4,4))
end
ax = cfig[:add_subplot](a,b,c)
ax[:set_aspect]("equal")
#ax[:axis]("off")
for i = 1:N1
for j = 1:N2
col = mod(i+j,2)==0 ? .7*[1,1,1] : .9*[1,1,1]
c = patch.Rectangle([i-.5,j-.5],1,1,fc=col)
ax[:add_artist](c)
end
end
dx = 0.45
axis([dx,N1+1-dx,dx,N2+1-dx])
tight_layout()
return
end
# plot a binary vector of length M as a sequence of moves
function plotmoves(moves,cfig=[],a=1,b=1,c=1,ms=10)
drawboard(cfig,a,b,c)
for m in find(moves)
p,q,r,s = move2coords(m)
plot([p,r],[q,s],"r.-",markersize=ms)
end
return
end
;
# simple example: 6x6 board, using knight moves (2,1)
initializeEverything(6,6,2,1);
using JuMP, Gurobi
m = Model(solver=GurobiSolver(OutputFlag=0))
@variable(m, x[1:M], Bin)
@variable(m, z[1:N], Bin)
@variable(m, w[1:N], Bin)
## zero out all odd coordinates (use for parity-preserving move types only)
#for i = 1:N1
# for j =1:N2
# if mod(i+j,2) == 0
# c = coord2cell(i,j)
# @constraint(m, z[c] == 0)
# @constraint(m, w[c] == 0)
# end
# end
#end
H = maximum(sum(C,1)) # M-trick upper bound
@constraint(m, sum(w) ≤ 2) # at most two special nodes (start and end)
@constraint(m, B*x .== 2*z - w) # at most two paths incident on each cell (except start and end)
@constraint(m, C*x .≤ H*(1-x)) # no intersecting paths (replace 37 by max())
@objective(m, Max, sum(x))
@time solve(m)
xsol = getvalue(x)
cfig = figure(figsize=(4,4))
plotmoves(xsol,cfig)
getobjectivevalue(m)
2.862246 seconds (170 allocations: 132.156 KB)
17.0
i = 0
while true
i = i+1
cycles = getCycles(xsol)
println(cycles)
println(getobjectivevalue(m))
if isempty(cycles)
cfig = figure(figsize=(4,4))
plotmoves(xsol,cfig)
break
end
for c in cycles
L = length(c)-1
@constraint(m, sum(x[c[1:L]]) ≤ L-1)
end
println(" ")
println("ITERATION ", i)
@time solve(m)
xsol = getvalue(x)
end
# print out details of solution (coordinate list and edge list)
println(" ")
sol_coords = sol2coordlist(xsol)
show(sol_coords')
println(" ")
show(find(xsol))
Any[Any[5,6,23,55,57,77,78,64,63,30,10,9,5],Any[17,48,49,37,17]] 17.0 ITERATION 1 3.389343 seconds (58 allocations: 18.750 KB) Any[] 17.0
Any[(4,4) (2,5) (3,3) (1,4) (2,6) (4,5) (5,3) (3,4) (4,2) (2,3) (1,1) (3,2) (5,1) (4,3) (6,2) (5,4) (4,6)] [3,4,11,13,19,21,28,33,35,41,48,49,55,65,67,73,75]
"""
CONFIRMED SOLUTIONS for camels (1,3)
4x4 --> (3, optimal)
Any[(1,1) (4,2) (1,3)]
[2,3,12]
5x5 --> (5, optimal)
Any[(1,4) (4,3) (1,2) (4,1) (5,4)]
[1,7,10,19,25]
6x6 --> (7, optimal)
Any[(1,1) (2,4) (3,1) (6,2) (5,5) (4,2) (3,5)]
[6,15,16,35,40,41,60]
7x7 --> (13, optimal)
Any[(5,4) (2,5) (5,6) (2,7) (1,4) (2,1) (3,4) (6,3) (3,2) (6,1) (7,4) (4,5) (7,6)]
[3,15,17,22,24,35,41,48,71,76,79,81,88]
8x8 --> (17, optimal)
Any[(6,3) (7,6) (4,5) (5,8) (8,7) (7,4) (8,1) (5,2) (6,5) (3,4) (4,7) (1,8) (2,5) (1,2) (4,1) (5,4) (2,3)]
[1,7,33,35,41,47,61,63,74,88,90,101,115,117,118,129,132]
"""
"""
CONFIRMED SOLUTIONS for zebra (2,3)
4x4 --> (2, optimal)
Any[(4,3) (1,1)]
[2,6]
5x5 --> (4, optimal)
Any[(5,3) (2,1) (4,4) (1,2)]
[4,11,12,19]
6x6 --> (7, optimal)
Any[(3,6) (1,3) (4,5) (2,2) (5,4) (3,1) (6,3)]
[6,17,18,28,30,39,43]
7x7 --> (11, optimal)
Any[(3,1) (1,4) (4,6) (2,3) (5,1) (3,4) (6,2) (4,5) (7,7) (5,4) (7,1)]
[2,5,9,14,15,21,35,52,54,79,80]
8x8 --> (17, optimal)
Any[(6,2) (3,4) (1,1) (4,3) (7,1) (5,4) (8,2) (6,5) (8,8) (5,6) (7,3) (4,5) (1,3) (3,6) (6,4) (4,7) (1,5)]
[4,5,15,17,23,25,42,43,49,59,61,67,85,86,103,119,120]
"""
"""
CONFIRMED SOLUTIONS for giraffe (1,4)
5x5 --> (4, optimal)
Any[(5,1) (1,2) (5,3) (1,4)]
[1,4,5,16]
6x6 --> (7, optimal)
Any[(6,1) (2,2) (6,3) (2,4) (6,5) (2,6) (1,2)]
[2,7,8,10,26,29,31]
7x7 --> (10, optimal)
Any[(3,3) (7,2) (6,6) (2,5) (6,4) (2,3) (6,2) (2,1) (1,5) (5,6)]
[5,8,9,17,19,23,49,51,52,56]
8x8 --> (15, optimal)
Any[(4,1) (8,2) (4,3) (8,4) (7,8) (3,7) (7,6) (3,5) (7,4) (3,3) (7,2) (3,1) (2,5) (6,6) (2,7)]
[7,8,11,12,23,24,28,32,63,66,71,73,76,109,110]
"""
"""
CONFIRMED SOLUTIONS for knights (1,2)
3x3 --> (2, optimal)
(2,1) (3,3)
[7,8]
4x4 --> (5, optimal)
(1,2) (3,1) (4,3) (2,4) (3,2)
[1,6,13,18,19]
5x5 --> (10, optimal) -- w[1] = 1
Any[(1,1) (2,3) (1,5) (3,4) (5,5) (4,3) (2,4) (3,2) (5,3) (4,1)]
[2,9,19,20,24,25,35,36,47,48]
6x6 --> (17, optimal) -- w[15] = 1 {corner not optimal!}
(2,1) (1,3) (3,2) (5,1) (4,3) (6,2) (5,4) (6,6) (4,5) (5,3) (3,4) (4,2) (2,3) (1,5) (3,6) (4,4) (2,5)
[4,9,10,13,19,21,33,35,41,45,48,49,55,69,70,79,80]
7x7 --> (24, optimal) -- required 9 subtour eliminations
(2,1) (1,3) (3,4) (2,2) (4,3) (3,1) (5,2) (7,1) (6,3) (4,2) (5,4) (3,3) (4,5) (2,4) (1,6) (3,7) (2,5) (4,6) (6,7) (7,5) (5,6) (4,4) (6,5) (5,3)
[7,8,11,20,22,29,30,38,40,46,48,54,64,66,72,74,77,88,90,93,104,106,117,118]
8x8 --> (35, optimal) -- required 7 subtour eliminations
(2,3) (1,1) (3,2) (5,3) (4,1) (6,2) (8,1) (7,3) (5,2) (6,4) (4,3) (2,2) (3,4) (1,3) (2,5) (3,7) (1,6) (2,8) (4,7) (5,5) (3,6) (2,4) (4,5) (3,3) (5,4) (7,5) (6,3) (8,4) (7,6) (8,8) (6,7) (4,8) (5,6) (7,7) (6,5)
[3,9,10,15,24,26,28,35,36,44,46,54,56,58,64,67,74,76,84,86,88,96,99,113,122,124,127,138,140,145,147,153,155,167,168]
8x5 --> (19, optimal)
Any[(2,3) (3,1) (1,2) (2,4) (3,2) (5,1) (6,3) (4,2) (5,4) (3,3) (2,5) (4,4) (6,5) (5,3) (7,4) (6,2) (8,1) (7,3) (8,5)]
[1,4,10,16,30,32,35,41,42,52,54,60,62,68,69,82,84,89,90]
10x4 --> (17, optimal)
Any[(2,4) (1,2) (3,1) (2,3) (4,4) (3,2) (5,1) (4,3) (6,2) (5,4) (7,3) (6,1) (8,2) (7,4) (9,3) (8,1) (10,2)]
[1,4,13,16,20,27,29,38,46,53,60,62,65,67,73,75,81]
12x4 --> (21, optimal)
Any[(2,1) (3,3) (1,2) (2,4) (4,3) (3,1) (5,2) (4,4) (6,3) (5,1) (7,2) (6,4) (8,3) (7,1) (9,2) (8,4) (10,3) (9,1) (11,2) (10,4) (12,3)]
[7,11,15,19,26,28,30,38,46,54,61,65,67,73,75,81,83,89,91,97,99]
""";
cfig = figure(figsize=(8,8))
sizes = [5,6,7,8]
edgelists = (
[2,9,19,20,24,25,35,36,47,48],
[4,9,10,13,19,21,33,35,41,45,48,49,55,69,70,79,80],
[7,8,11,20,22,29,30,38,40,46,48,54,64,66,72,74,77,88,90,93,104,106,117,118],
[3,9,10,15,24,26,28,35,36,44,46,54,56,58,64,67,74,76,84,86,88,96,99,113,122,124,127,138,140,145,147,153,155,167,168]
)
for i = 1:4
initializeEverything(sizes[i],sizes[i],2,1)
edgelist = edgelists[i]
ysol = spzeros(M,1)
ysol[edgelist]=1
plotmoves(ysol,cfig,2,2,i,7)
axis("off")
title(string(sizes[i], "x", sizes[i], " longest knight path = ", length(edgelist)))
end
tight_layout()
savefig("nonintersecting_knight_paths.png")
cfig = figure(figsize=(8,3))
edgelist = [7,11,15,19,26,28,30,38,46,54,61,65,67,73,75,81,83,89,91,97,99]
initializeEverything(12,4,2,1)
ysol = spzeros(M,1)
ysol[edgelist]=1
plotmoves(ysol,cfig,1,1,1,7)
axis("off")
title(string("12x4 longest knight path = ", length(edgelist)))
tight_layout()
savefig("nonintersecting_rect_path.png")
cfig = figure(figsize=(8,8))
types = ["knight", "camel", "zebra", "giraffe"]
movestyle = [(2,1), (3,1), (3,2), (4,1)]
edgelists = (
[3,9,10,15,24,26,28,35,36,44,46,54,56,58,64,67,74,76,84,86,88,96,99,113,122,124,127,138,140,145,147,153,155,167,168],
[1,7,33,35,41,47,61,63,74,88,90,101,115,117,118,129,132],
[4,5,15,17,23,25,42,43,49,59,61,67,85,86,103,119,120],
[7,8,11,12,23,24,28,32,63,66,71,73,76,109,110]
)
for i = 1:4
initializeEverything(8,8,movestyle[i][1],movestyle[i][2])
edgelist = edgelists[i]
ysol = spzeros(M,1)
ysol[edgelist]=1
plotmoves(ysol,cfig,2,2,i,7)
axis("off")
title(string("8x8 longest ", types[i], " path = ", length(edgelist)))
end
tight_layout()
savefig("nonintersecting_fairy_paths.png")