julia系列17: tsp问题代码整理

1. 常用库和基础函数

这里是优化版的函数：

using TSPLIB,LKH,Distances,PyPlot
MaxNum = 10000
tsp=readTSPLIB(:att48)
dist = [round.(Int,euclidean(tsp.nodes[i,:],tsp.nodes[j,:])) for i in 1:tsp.dimension,j in 1:tsp.dimension];
pos(tsp::TSP,t::Vector{Int},i::Int)=tsp.nodes[t[i]==n+1 ? 1 : t[i],:]
total_length(tsp::TSP,t::Vector{Int64})=sum([round.(Int,euclidean(pos(tsp,t,i),pos(tsp,t,i+1))) for i in 1:length(t)-1])
function minimum2(vec)
    i2 = (-1,-1);v2=(MaxNum,MaxNum);for (i,v) in enumerate(vec);if v2[1]>v;v2=(v,v2[1]);i2=(i,i2[1]);elseif v2[2]>v;v2=(v2[1],v);i2=(i2[1],i);end;end
    return i2,v2
end
function greedyConstruct(dist)
    path = [1];for _ in 1:size(dist)[1];c_dist = dist[path[end],:];c_dist[path].=MaxNum;next_id = findmin(c_dist)[2];push!(path,next_id);end
    return path
end
function plot_path(tsp,path)
    t = tsp.nodes[path,:];for p in path;text(tsp.nodes[p,1],tsp.nodes[p,2],p);end;scatter(t[:,1],t[:,2]);plot(t[:,1],t[:,2])
end
function plot_tree(tsp,mst)
    s = tsp.nodes;for p in 1:tsp.dimension;text(s[p,1],s[p,2],p);end;scatter(s[:,1],s[:,2]);for m in mst;plot(s[[m...],1],s[[m...],2],c="b");end
end
function minspantree(dm::AbstractMatrix{T}) where {T<:Real} # accepts views
    mst_edges = Vector{Tuple{Int, Int}}();mst_cost = zero(T);n = size(dm, 1)
    tree_dists = dm[1,:];closest_tree_verts = ones(Int, n);tree_dists[1] = MaxNum
    for _ in 1:(n-1) # need to add n - 1 other verts to tree
        cost, newvert = findmin(tree_dists)
        treevert = closest_tree_verts[newvert];mst_cost += cost
        push!(mst_edges, tuple(sort([treevert, newvert])...))
        tree_dists[newvert] = MaxNum
        for i in 1:n
            if tree_dists[i] >= MaxNum;continue;end
            if tree_dists[i] > dm[i, newvert];tree_dists[i] = dm[i, newvert];closest_tree_verts[i] = newvert;end
        end
    end
    return mst_edges, mst_cost
end

这里是老版的函数

using TSPLIB
using LKH
using Distances
using Random

function pos(tsp::TSP,t::Vector{Int},i::Int)
    n = length(t)
    i = mod(i, n)
    if i == 0
        return tsp.nodes[t[n],1:2]
    else
        return tsp.nodes[t[i],1:2]
    end
end

function pos(b::Vector{Int},i::Int)
    n_b = length(b)
    i = mod(i, n_b)
    if i == 0
        return b[n_b]
    else
        return b[i]
    end
end

function tour_length(tsp::TSP,t::Vector{Int64})
    l = 0
    for i in 1:length(t)
        l += Int.(round.(euclidean(pos(tsp,t,i),pos(tsp,t,i+1))))
    end
    l
end

function greedyConstruct(tsp)
    n = tsp.dimension
    uv = Set{Int}(2:n)
    t = Array{Int64}(undef, n)
    t[1] = 1
    for i in 2:n
        cur = t[i-1]
        mindist = Inf
        ind = 1
        for j in collect(uv)
            dis = round.(Int, euclidean(tsp.nodes[cur,1:2],tsp.nodes[j,1:2]))
            if dis < mindist
                mindist = dis
                ind = j
            end
        end
        t[i] = ind
        pop!(uv,ind)
    end
    return t
end

num = 30
f = "bbz25234"
tsp = readTSP("/Users/chen/Downloads/vlsi/"*f*".tsp")
n = tsp.dimension

# 解析问题文件
t = parse.(Int, readlines("xrb14233.tour")[6:end-3])
totaldist(t)

用LKH求解器先给出一些初始解：

println("start...")
for i in 1:num
    @time t,v = solve_tsp("/Users/chen/Downloads/vlsi/"*f*".tsp";RUNS=1,MAX_TRIALS =150, MAX_SWAPS=100,SEED=i,PI_FILE="pi.txt",TIME_LIMIT=100,CANDIDATE_FILE = "cand.txt") #
    @printf("%d",v)
    ts[i,:]=t
    results[i]= v
end
println("end")

io = open("route.out","w")
for i in 1:size(ts)[1]
    write(io, string(ts[i,:])*"\n")
end
close(io)

# route是初始路线合集，格式如下图
ts = Array{Int32}(undef, (num,n))
results = Array{Int32}(undef,num)
lines = readlines("route.out")
println("start...")
for i in 1:num
    ts[i,:]=Meta.parse(lines[i]) |> eval
    results[i]= totaldist(ts[i,:])
end
println("end")

在这里插入图片描述
接下来设计一下基本数据结构。上述路径使用Array来存储的，

2. 使用MPS GPU加速的2-opt算法

核函数如下：

using Metal
function twooptkernel(dist, t, x, y, z)
    k = thread_position_in_grid_1d()
    n = Int32(size(z)[1])
    if k <= n
        i = x[k]
        j = y[k]
        z[k] = dist[t[i], t[j]]+ dist[t[i==n ? 1 : i+1], t[j==n ? 1 : j+1]] - dist[t[i], t[i==n ? 1 : i+1]] - dist[t[j], t[j==n ? 1 : j+1]]
    end
    return
end

与cpu版本进行对比：

using Random
m = 4000000
x = rand(1:n, m)
y = rand(1:n, m)
z = similar(x)
mx = MtlArray{Int32}(x);
my = MtlArray{Int32}(y);
mz = similar(mx);
mt = MtlArray{Int32}(t);
@time @metal threads=64 grid=6250  twooptkernel(mdist, mt, mx,my,mz)

function cpu(x,y,z)
    for k in 1:m
        i = x[k]
        j = y[k]
        z[k] = dist[t[i], t[j]]+ dist[t[i==n ? 1 : i+1], t[j==n ? 1 : j+1]] - dist[t[i], t[i==n ? 1 : i+1]] - dist[t[j], t[j==n ? 1 : j+1]]
    end
end
@time cpu(x,y,z)

分别为
0.000656 seconds (254 allocations: 6.031 KiB)
以及
5.414779 seconds (74.26 M allocations: 1.168 GiB, 0.61% compilation time)
m = 40000

完整GPU调用代码如下：

for i in 1:200000
    m = 400000
    mx = MtlArray{Int32}(rand(1:n, m));
    my = MtlArray{Int32}(rand(1:n, m));
    mz = similar(mx);
    mt = MtlArray{Int32}(t);
    @metal threads=64 grid=6250  twooptkernel(mdist, mt, mx,my,mz)
    z = Array(mz)
    v,ind = findmin(z)
    if v < 0
        st,ed = mx[ind], my[ind]
        if st>ed
           st,ed = ed,st 
        end
        reverse!(t, st+1, ed)
    end
    if i % 1000 == 0
        println(i,":",totaldist(t))
    end
end

3. LKH算法包

调用方式极其简单：

using LKH
opt_tour, opt_len = solve_tsp(dist)
opt_tour, opt_len = solve_tsp(x, y; dist="EUC_2D")
opt_tour, opt_len = solve_tsp("gr17.tsp")
opt_tour, opt_len = solve_tsp(dist; INITIAL_TOUR_ALGORITHM="GREEDY", RUNS=5)

可用参数参考这篇：https://blog.csdn.net/kittyzc/article/details/114388836
一些常用参数：

RUNS=1, 重复运行次数
MAX_TRIALS =1000, 每次运行最多尝试交换的次数，默认值等于点数。
MAX_SWAPS=1000, 每次交换最多运行swap的次数。
INITIAL_TOUR_FILE = "temp.out", 初始化路径地址
SEED=inum, 每次运行给一个新的seed
PI_FILE="pi-"*f*".txt", 新的cost，计算一次后重复调用
TIME_LIMIT=100, 时间限制
CANDIDATE_FILE = "cand-"*f*".txt"， 每个点的近邻，计算一次后重复调用即可

4. 简单遗传算法

每次随机选两条路径（t_ids），并随机从第一条路径t1上选两个点（e_ids），将其中的点按照第二条路径t2的顺序重新排列，产生新的路径t3，对此路径使用LKH算法进行优化，如果距离比t1或者t2短，则进行替换。

iter = 500
t_ids = rand(1:size(ts)[1], (iter,2))
e_ids = rand(1:n, (iter,2))
for inum in 1:iter
    t1i,t2i = t_ids[inum,1], t_ids[inum,2] 
    i1, i2 = e_ids[inum,1], e_ids[inum,2] 
    if i1>i2
        i1,i2 = i2,i1
    end 
    t1 = ts[t1i,:]
    t2 = ts[t2i,:]
    t3 = Array{Int32}(undef, n)
    t3[1:i1] = t1[1:i1]
    t3[i2:end] = t1[i2:end]
    t2ind = 1
    for t3ind in i1:i2
        while !(t2[t2ind] in t1[i1:i2])
            t2ind += 1
        end
        t3[t3ind] = t2[t2ind]
        t2ind+=1
    end
    io = open("temp.out","w")
    write(io, "TYPE : TOUR\nDIMENSION : "*string(n)*"\nTOUR_SECTION\n")
    for i in t3
        write(io, string(i)*"\n")
    end
    write(io, "-1\nEOF\n")
    close(io)
    t3, s3 = solve_tsp("/Users/chen/Downloads/vlsi/"*f*".tsp";RUNS=1,MAX_TRIALS =1000, MAX_SWAPS=1000,INITIAL_TOUR_FILE = "temp.out",SEED=inum,PI_FILE="pi.txt",TIME_LIMIT=100,CANDIDATE_FILE = "cand.txt") #
    sk, skind = findmin(results)

    # replace parents
    if s3 < sk
        print("*",inum,": from ",sk," to ",s3,", ")
        io = open("best"*string(s3)*".out","w")
        write(io, "TYPE : TOUR\nDIMENSION : "*string(n)*"\nTOUR_SECTION\n")
        for i in t3
            write(io, string(i)*"\n")
        end
        write(io, "-1\nEOF\n")
        close(io)
    else
        print(" ",inum,":")
    end
    a,b = t1i,t2i
    if s3 < results[a] <= results[b]
        ts[b,:]=t3
        results[b] = s3
        println(" replace ",b)
    elseif s3 < results[b] <= results[a] 
        ts[a,:]=t3
        results[a] = s3
        println(" replace ",a)
    elseif findmin(results)[1]==findmax(results)[1]
        println(" finish")
        break
    else
        println(" pass")
    end
end

5. ALNS算法

核心代码如下，首先生成breaks，打断这些点位连接的边，然后使用solve_tsp_part重新进行组合：

b_try_num = 1000
b_point_num = 400
@assert b_point_num*2 <= tsp.dimension
b_id_lists = sort(rand(0:tsp.dimension - 2*b_point_num, (b_try_num,b_point_num)),dims = 2)

for i in 1:b_try_num
    b = b_id_lists[i,:] + 2*Vector(1:b_point_num)
    t_imp,v_imp = solve_tsp_part(tsp,t,b)
    if v_imp>0
        t = t_imp
    end
    if i % 100 == 0
        println(tour_length(tsp,t_imp))
    end
end

用于重构的solve_tsp_part和恢复路径的restore_tour代码如下：

function restore_tour(t::Vector{Int64},b::Vector{Int64},b_imp::Vector{Int64})
    n = length(t)
    n_b = length(b)
    t_imp = zeros(Int64,n)
    cur_imp = 1
    for i in 1:n_b
        b_i = div(b_imp[2*i]+1,2) # 在b中的下标
        r = false
        if b_imp[2*i]%2==1
            r = true
        end
        st = pos(b,b_i)
        ed = pos(b,b_i+1)-1
        if ed > st
            mov_imp = ed - st
            if r
                t_imp[cur_imp:cur_imp+mov_imp] = reverse(t[st:ed])
                # for j in reverse(st:ed)
                #     print(t[j],",")
                # end
            else   
                t_imp[cur_imp:cur_imp+mov_imp] = t[st:ed]
                # for j in st:ed
                #     print(t[j],",")
                # end
            end
            cur_imp += mov_imp+1
        else
            if r
                mov_imp = ed - 1
                t_imp[cur_imp:cur_imp+mov_imp] = reverse(t[1:ed])
                cur_imp += mov_imp+1
                mov_imp = n - st
                t_imp[cur_imp:cur_imp+mov_imp] = reverse(t[st:n])
                cur_imp += mov_imp+1
                # for j in reverse(1:ed)
                #     print(t[j],",")
                # end
                # for j in reverse(st:n)
                #     print(t[j],",")
                # end
            else
                mov_imp = n - st
                t_imp[cur_imp:cur_imp+mov_imp] = t[st:n]
                cur_imp += mov_imp+1
                mov_imp = ed - 1
                t_imp[cur_imp:cur_imp+mov_imp] = t[1:ed]
                cur_imp += mov_imp+1
                # for j in st:n
                #     print(t[j],",")
                # end
                # for j in 1:ed
                #     print(t[j],",")
                # end
            end
        end
    end
    t_imp
end

function check_bimp(b_imp::Vector{Int64})
    n_b = div(length(b_imp),2)
    for i in 1:n_b
        if abs(b_imp[2*i]-b_imp[2*i-1])!=1 || div(b_imp[2*i]+1,2)!= div(b_imp[2*i-1]+1,2)
            print(i,":",b_imp[2*i],",",b_imp[2*i-1])
            return false
        end
    end
    return true
end

function solve_tsp_part(tsp::TSP,t::Vector{Int64},b::Vector{Int64}) # b for breaks
    # cal distance matrix
    n = length(t)
    n_b = length(b)
    dist_matrix_b = zeros(Int64,2*n_b,2*n_b)
    for i in 1:n_b
        for j in i+1:n_b
            dist_matrix_b[2*i-1,2*j-1] = dist_matrix_b[2*j-1,2*i-1] = Int(round(euclidean(pos(tsp,t,pos(b,i)),pos(tsp,t,pos(b,j)))))
            dist_matrix_b[2*i-1,2*j] = dist_matrix_b[2*j,2*i-1] = Int(round(euclidean(pos(tsp,t,pos(b,i)),pos(tsp,t, pos(b,j+1)-1 ))))
            dist_matrix_b[2*i,2*j-1] = dist_matrix_b[2*j-1,2*i] = Int(round(euclidean(pos(tsp,t,pos(b,i+1)-1  ),pos(tsp,t,b[j]))))
            dist_matrix_b[2*i,2*j] = dist_matrix_b[2*j,2*i] = Int(round(euclidean(pos(tsp,t,pos(b,i+1)-1 ),pos(tsp,t,pos(b,j+1)-1 ))))
        end
    end
    max_element = maximum(dist_matrix_b)
    for i in 1:n_b
        for j in i+1:n_b
            dist_matrix_b[2*i-1,2*j-1] += max_element
            dist_matrix_b[2*j-1,2*i-1] += max_element
            dist_matrix_b[2*i-1,2*j] += max_element
            dist_matrix_b[2*j,2*i-1] += max_element
            dist_matrix_b[2*i,2*j-1] += max_element
            dist_matrix_b[2*j-1,2*i] += max_element
            dist_matrix_b[2*i,2*j] += max_element
            dist_matrix_b[2*j,2*i] += max_element
        end
    end
    # for i in 1:n_b-1
    #     dist_matrix_b[2*i,2*i-1] = dist_matrix_b[2*i-1,2*i] = -max_element
    # end
    b_imp,v_imp = solve_tsp(dist_matrix_b;runs = 1,MAX_TRIALS =50, MAX_SWAPS=400)
    v_imp -= max_element*n_b
    if check_bimp(b_imp)
        v = 0
        for i in 1:n_b
            v +=  Int(round(euclidean(pos(tsp,t,pos(b,i)),pos(tsp,t, pos(b,i)-1 ))))
        end
        imp = v-v_imp
        if imp > 0
            return restore_tour(t,b,b_imp), imp
        else
            return t, 0
        end
    else
        return b_imp, -1
    end
end