# Knights, zebras, camels, etc. Non-intersecting paths!¶

This program computes the longest path that doesn't intersect itself.

In [2]:
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


#### Simple functions for index manipulation¶

In [3]:
# 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
;


#### Matrices that encode constraints¶

In [17]:
# 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
;


#### Cycle detection¶

In [18]:
# 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 functions for displaying solution data¶

In [19]:
# 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
;


#### Plotting functions¶

In [20]:
# function that draws the board
function drawboard(cfig=[],a=1,b=1,c=1)
if cfig == []
cfig = figure(figsize=(4,4))
end
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)
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
;


# Initialize board and solve as MIP (solution might have subtours)¶

In [24]:
# 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)

Out[24]:
17.0

## Perform adaptive subtour elimination until we achieve optimality¶

In [25]:
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]
In [13]:
"""
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]
""";


# Plot comparing different length knight tours¶

In [26]:
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")


## Plot for a non-square board¶

In [27]:
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")


## Plot comparing different chess pieces¶

In [28]:
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")

In [ ]: