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;endreturn 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);endreturn 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 viewsmst_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] = MaxNumfor _ in 1:(n-1) # need to add n - 1 other verts to treecost, newvert = findmin(tree_dists)treevert = closest_tree_verts[newvert];mst_cost += costpush!(mst_edges, tuple(sort([treevert, newvert])...))tree_dists[newvert] = MaxNumfor i in 1:nif tree_dists[i] >= MaxNum;continue;endif tree_dists[i] > dm[i, newvert];tree_dists[i] = dm[i, newvert];closest_tree_verts[i] = newvert;endendendreturn mst_edges, mst_cost
end
這里是老版的函數
using TSPLIB
using LKH
using Distances
using Randomfunction pos(tsp::TSP,t::Vector{Int},i::Int)n = length(t)i = mod(i, n)if i == 0return tsp.nodes[t[n],1:2]elsereturn tsp.nodes[t[i],1:2]end
endfunction pos(b::Vector{Int},i::Int)n_b = length(b)i = mod(i, n_b)if i == 0return b[n_b]elsereturn b[i]end
endfunction tour_length(tsp::TSP,t::Vector{Int64})l = 0for i in 1:length(t)l += Int.(round.(euclidean(pos(tsp,t,i),pos(tsp,t,i+1))))endl
endfunction greedyConstruct(tsp)n = tsp.dimensionuv = Set{Int}(2:n)t = Array{Int64}(undef, n)t[1] = 1for i in 2:ncur = t[i-1]mindist = Infind = 1for j in collect(uv)dis = round.(Int, euclidean(tsp.nodes[cur,1:2],tsp.nodes[j,1:2]))if dis < mindistmindist = disind = jendendt[i] = indpop!(uv,ind)endreturn t
endnum = 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,:]=tresults[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:numts[i,:]=Meta.parse(lines[i]) |> evalresults[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 <= ni = 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]]endreturn
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:mi = 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:200000m = 400000mx = 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 < 0st,ed = mx[ind], my[ind]if st>edst,ed = ed,st endreverse!(t, st+1, ed)endif i % 1000 == 0println(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:itert1i,t2i = t_ids[inum,1], t_ids[inum,2] i1, i2 = e_ids[inum,1], e_ids[inum,2] if i1>i2i1,i2 = i2,i1end t1 = ts[t1i,:]t2 = ts[t2i,:]t3 = Array{Int32}(undef, n)t3[1:i1] = t1[1:i1]t3[i2:end] = t1[i2:end]t2ind = 1for t3ind in i1:i2while !(t2[t2ind] in t1[i1:i2])t2ind += 1endt3[t3ind] = t2[t2ind]t2ind+=1endio = open("temp.out","w")write(io, "TYPE : TOUR\nDIMENSION : "*string(n)*"\nTOUR_SECTION\n")for i in t3write(io, string(i)*"\n")endwrite(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 parentsif s3 < skprint("*",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 t3write(io, string(i)*"\n")endwrite(io, "-1\nEOF\n")close(io)elseprint(" ",inum,":")enda,b = t1i,t2iif s3 < results[a] <= results[b]ts[b,:]=t3results[b] = s3println(" replace ",b)elseif s3 < results[b] <= results[a] ts[a,:]=t3results[a] = s3println(" replace ",a)elseif findmin(results)[1]==findmax(results)[1]println(" finish")breakelseprintln(" 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_numb = b_id_lists[i,:] + 2*Vector(1:b_point_num)t_imp,v_imp = solve_tsp_part(tsp,t,b)if v_imp>0t = t_impendif i % 100 == 0println(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 = 1for i in 1:n_bb_i = div(b_imp[2*i]+1,2) # 在b中的下標r = falseif b_imp[2*i]%2==1r = trueendst = pos(b,b_i)ed = pos(b,b_i+1)-1if ed > stmov_imp = ed - stif rt_imp[cur_imp:cur_imp+mov_imp] = reverse(t[st:ed])# for j in reverse(st:ed)# print(t[j],",")# endelse t_imp[cur_imp:cur_imp+mov_imp] = t[st:ed]# for j in st:ed# print(t[j],",")# endendcur_imp += mov_imp+1elseif rmov_imp = ed - 1t_imp[cur_imp:cur_imp+mov_imp] = reverse(t[1:ed])cur_imp += mov_imp+1mov_imp = n - stt_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],",")# endelsemov_imp = n - stt_imp[cur_imp:cur_imp+mov_imp] = t[st:n]cur_imp += mov_imp+1mov_imp = ed - 1t_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],",")# endendendendt_imp
endfunction check_bimp(b_imp::Vector{Int64})n_b = div(length(b_imp),2)for i in 1:n_bif 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 falseendendreturn true
endfunction solve_tsp_part(tsp::TSP,t::Vector{Int64},b::Vector{Int64}) # b for breaks# cal distance matrixn = length(t)n_b = length(b)dist_matrix_b = zeros(Int64,2*n_b,2*n_b)for i in 1:n_bfor j in i+1:n_bdist_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 ))))endendmax_element = maximum(dist_matrix_b)for i in 1:n_bfor j in i+1:n_bdist_matrix_b[2*i-1,2*j-1] += max_elementdist_matrix_b[2*j-1,2*i-1] += max_elementdist_matrix_b[2*i-1,2*j] += max_elementdist_matrix_b[2*j,2*i-1] += max_elementdist_matrix_b[2*i,2*j-1] += max_elementdist_matrix_b[2*j-1,2*i] += max_elementdist_matrix_b[2*i,2*j] += max_elementdist_matrix_b[2*j,2*i] += max_elementendend# 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# endb_imp,v_imp = solve_tsp(dist_matrix_b;runs = 1,MAX_TRIALS =50, MAX_SWAPS=400)v_imp -= max_element*n_bif check_bimp(b_imp)v = 0for i in 1:n_bv += Int(round(euclidean(pos(tsp,t,pos(b,i)),pos(tsp,t, pos(b,i)-1 ))))endimp = v-v_impif imp > 0return restore_tour(t,b,b_imp), impelsereturn t, 0endelsereturn b_imp, -1end
end
)
)
