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)

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))

"""
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")