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

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

by: Laurent Lessard

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 = 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
;

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