import JuMP as J
import SCIP
import GLPK
import Plots
import Combinatorics
import Random

function tsp_mtz(n, dists)
    m = J.Model(GLPK.Optimizer)
    
    # variable: connection i->j is used.
    J.@variable(m, tour[i=1:n, j=1:n], Bin)
    # variable: sequence number of visited nodes
    J.@variable(m, sequence[i=1:n] >= 0, Int)
        
    # facility 1 is the start.
    J.fix(sequence[1], 1, force=true)

    # Miller/Tucker/Zemlin (1960)
    
    J.@objective(m, Min, sum(tour[i,j] * dists[i,j] for i = 1:n for j = 1:n))
    
    # no self-visiting.
    J.@constraint(m, facility_no_self[i=1:n], tour[i,i] == 0)
    # for each facility: exactly one incoming edge
    J.@constraint(m, facility_incoming[i=1:n], sum(tour[:, i]) == 1)
    # for each facility: exactly one outgoing edge
    J.@constraint(m, facility_outgoing[i=1:n], sum(tour[i, :]) == 1)

    # reformulated:
    J.@constraint(m, monoton[i=1:n, j=2:n], sequence[i] - sequence[j] + n * tour[i,j] <= n-1)

    J.optimize!(m)
    
    J.value.(tour), J.objective_value(m)
end

tsp_mtz (generic function with 1 method)

function tsp_dfj_subsets(n, dists)
    m = J.Model(GLPK.Optimizer)
    
    # variable: connection i->j is used.
    J.@variable(m, tour[i=1:n, j=1:n], Bin)
    
    J.@objective(m, Min, sum(tour[i,j] * dists[i,j] for i = 1:n for j = 1:n))
    
    # no self-visiting.
    J.@constraint(m, facility_no_self[i=1:n], tour[i,i] == 0)
    # for each facility: exactly one incoming edge
    J.@constraint(m, facility_incoming[i=1:n], sum(tour[:, i]) == 1)
    # for each facility: exactly one outgoing edge
    J.@constraint(m, facility_outgoing[i=1:n], sum(tour[i, :]) == 1)

    vertices = collect(1:n)
    for ps in Combinatorics.powerset(vertices)
        if length(ps) < n && length(ps) > 0
           J.@constraint(m, sum(tour[i,j] for i in ps for j in ps if i != j) <= length(ps) - 1) 
        end
    end

    J.optimize!(m)
    
    J.value.(tour), J.objective_value(m)
end

tsp_dfj_subsets (generic function with 1 method)

function tsp_dfj_short(n, dists)
    m = J.Model(GLPK.Optimizer)
    
    # variable: connection i->j is used.
    J.@variable(m, tour[i=1:n, j=i:n], Bin)
    
    J.@objective(m, Min, sum(tour[i,j] * dists[i,j] for i = 1:n for j = i:n))
    
    J.@constraint(m, enterexit[i=1:n], sum(tour[k, l] for k in 1:n for l in 1:n if (i == k || i == l) && k < l) == 2)
    
    for ps in Combinatorics.powerset(collect(1:n), 2)
        if length(ps) == n
            continue
        end
        J.@constraint(m, sum(tour[i,j] for i in 1:n for j in 1:n if (in(i, ps) + in(j, ps)) == 1 && i < j) >= 2)
    end
    
    J.optimize!(m)
    
    J.value.(tour), J.objective_value(m)
end

tsp_dfj_short (generic function with 1 method)

function run(n, solverfunc; asym=false)
    Random.seed!(1)
    facilities = rand(n, 2) .* 50

    dists = zeros(n, n)
    for i = 1:n
        for j = 1:n
            dists[i,j] = sqrt(sum((facilities[i, :] .- facilities[j, :]).^2))
        end
    end

    @time tour, obj = solverfunc(n, dists[1:n, 1:n])
    println("n = $n, obj = $obj")

    Plots.scatter(facilities[:,1], facilities[:,2], legend=false)
    for i in 1:n
        for j in 1:n
            if (asym || i < j) && tour[i,j] > 0
                Plots.plot!([facilities[i,1],facilities[j,1]], [facilities[i,2],facilities[j,2]],
                    color=:red, arrow=true)
            end
        end
    end
    Plots.current()
end

run (generic function with 1 method)

run(12, tsp_dfj_short, asym=false)

  0.118678 seconds (286.95 k allocations: 28.615 MiB, 24.58% gc time)
n = 12, obj = 190.56658600378748

run(12, tsp_dfj_subsets, asym=true)

  0.092210 seconds (145.42 k allocations: 31.570 MiB)
n = 12, obj = 190.56658600378773

run(12, tsp_mtz, asym=true)

  0.546728 seconds (14.85 k allocations: 942.938 KiB)
n = 12, obj = 190.5665860037875