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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/diannao/42499.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Groovy中,多种循环方式

1.最常规的for循环 for (def i 0; i < 5; i) { //这个i需要声明println "遍历输出${i}" }2.while循环 def j 0 while (j < 5) {println "遍历输出 ${j}"j }3.for in 循环 def list [0, 1, 2, 3, 4] //这个l无需声明 for (l in list) { printl…

uniapp跨域问题解决

找到menifest文件&#xff0c;在文件的最后添加如下代码&#xff1a; // h5 解决跨域问题"h5":{"devServer": {"proxy": {"/adminapi": {"target": "https://www.demo.com", // 目标访问网址"changeOrigin…

数据库的存储引擎,数据类型,约束条件,严格模式

【一】存储引擎 1.什么是存储引擎存储引擎可以理解为处理数据的不同方式 2.查看存储引擎show engines; 3.须知的引擎MyISAM5.5之前版本MySQL默认的存储引擎特点:存取数据的速度快 但是功能很少 安全性较低速度快因为自带索引InnoDB5.5之后版本MySQL默认的存储引擎特点:有诸多功…

工程施工合同无效但竣工交付,应当参照合同关于工程价款的约定计算折价补偿款

审判实践中&#xff0c;对于建设工程施工合同无效但工程竣工并交付使用的&#xff0c;应以何种标准计算折价补偿款的问题&#xff0c;认识不一致。【法官会议意见】&#xff1a;建设工程施工合同是承包人进行工程建设、交付工作成果即建设工程并由发包人支付价款的合同。建设工…

httpd目录显示乱码问题

vim /etc/httpd/conf/httpd.conf 在<Directory “/var/www/html”>下添加&#xff1a; IndexOptions CharsetUTF-8重启httpd: systemctl restart httpd.service还是不好看&#xff0c;调整下显示宽度&#xff0c;还是这个位置&#xff1a; <Directory “/var/www/ht…

区块链论文速读A会-ISSTA 2023(2/2)如何检测DeFi协议中的价格操纵漏洞

Conference&#xff1a;ACM SIGSOFT International Symposium on Software Testing and Analysis (ISSTA) CCF level&#xff1a;CCF A Categories&#xff1a;Software Engineering/System Software/Programming Languages Year&#xff1a;2023 第1~5篇区块链文章 请点击此…

09视图,触发器,事务,存储过程,函数,流程控制,索引,隔离机制,锁机制,三大范式

【一】视图 (1)视图须知概念 1.什么是视图&#xff1f; 视图就是通过查询得到一张虚拟表&#xff0c;然后保存下来&#xff0c;下次可以直接使用 2.为什么要用视图&#xff1f; 如果要频繁操作一张虚拟表(拼表组成)&#xff0c;就可以制作成视图&#xff0c;后续直接操作 注意…

IDEA 创建springboot项目杂记-更新中

一、工具使用杂记 1、使用maven 创建新springboot项目时&#xff0c;因为https://start.spring.io/ 连接不上项目无法创建。直接把脚手架地址换为国内的 http://start.aliyun.com

田忌赛马 贪心

本题是更难的那道,一场50 最低为o 第一行一个整数 &#x1d45b;n &#xff0c;表示他们各有几匹马&#xff08;两人拥有的马的数目相同&#xff09;。第二行 &#x1d45b;n 个整数&#xff0c;每个整数都代表田忌的某匹马的速度值&#xff08;0≤0≤ 速度值 ≤100≤1…

Python】从文本字符串中提取数字、电话号码、日期、网址的方法

关于从文本字符串中提取数字、电话号码、日期和网址的方法&#xff1a; 提取数字&#xff1a; 在 Python 中&#xff0c;使用正则表达式 \d 来匹配数字。 \d 表示匹配一个数字字符&#xff08;0-9&#xff09;。如果要匹配连续的数字&#xff0c;可以使用 \d 。 import re def …

C++面向对象的常见面试题目(一)

1. 面向对象的三大特征 &#xff08;1&#xff09;封装&#xff1a;隐藏对象的内部状态&#xff0c;只暴露必要的接口。 #include <iostream> #include <string>// 定义一个简单的类 Person class Person { private: // 私有成员&#xff0c;外部不可直接访问std…

Mac OS ssh 连接提示 Permission denied (publickey)

这错误有点奇葩&#xff0c;MacBook的IDE(vscode和pycharm)远程都连不上&#xff0c;terminal能连上&#xff0c;windows的pycharm能连上&#xff0c;见鬼了&#xff0c;所以肯定不是秘钥的问题了&#xff0c;查了好久竟然发现是权限的问题。。 chmod 400 ~/.ssh/id_rsa http…

华为机试HJ37统计每个月兔子的总数

华为机试HJ37统计每个月兔子的总数 题目&#xff1a; 想法&#xff1a; 上述题目实际是一个斐波那契数列&#xff0c;利用斐波那契数列对问题进行求解 input_number int(input())def fib(n):if n < 2:return 1else:n_1 1n_2 1count 2while count < n:n_1, n_2 n_…

【Android】【多屏】多屏异显异触调试技巧总结

这里写目录标题 如何获取多屏IDs获取多屏的size/density如何启动应用到指定DisplayId多屏截屏/录屏screencapscreenrecord发送按键到指定DisplayId 如何获取多屏IDs dumpsys display | grep mDisplayIdtrinket:/ # dumpsys display | grep mDisplayIdmDisplayId0mDisplayId2 t…

【AI资讯】可以媲美GPT-SoVITS的低显存开源文本转语音模型Fish Speech

Fish Speech是一款由fishaudio开发的全新文本转语音工具&#xff0c;支持中英日三种语言&#xff0c;语音处理接近人类水平&#xff0c;使用Flash-Attn算法处理大规模数据&#xff0c;提供高效、准确、稳定的TTS体验。 Fish Audio

区块链技术的应用场景和优势。

区块链技术具有广泛的应用场景和优势。 区块链技术的应用场景&#xff1a; 1. 金融服务&#xff1a;区块链可用于支付、跨境汇款、借贷和结算等金融服务&#xff0c;提高交易效率、降低成本并增强安全性。 2. 物联网&#xff08;IoT&#xff09;&#xff1a;区块链可以用于物…

机器学习Day12:特征选择与稀疏学习

1.子集搜索与评价 相关特征&#xff1a;对当前学习任务有用的特征 无关特征&#xff1a;对当前学习任务没用的特征 特征选择&#xff1a;从给定的特征集合中选择出相关特征子集的过程 为什么要特征选择&#xff1f; 1.任务中经常碰到维数灾难 2.去除不相关的特征能降低学习的…

Git注释规范

主打一个有用 代码的提交规范参考如下&#xff1a; init:初始化项目feat:新功能&#xff08;feature&#xff09;fix:修补bugdocs:文档&#xff08;documentation&#xff09;style:格式&#xff08;不影响代码运行的变动&#xff09;refactor:重构&#xff08;即不是新增功能…

NodeJs获取文件扩展名

path.extname 是 Node.js 路径模块 (path) 中的一个方法&#xff0c;用于获取文件路径的扩展名。扩展名是指文件名中最后一个 .&#xff08;点&#xff09;之后的部分&#xff0c;包括这个 .。 const path require(path);const filename example.txt; const ext path.extna…

计算机网络之令牌环

1.令牌环工作原理 令牌环&#xff08;Token Ring&#xff09;是一种局域网&#xff08;LAN&#xff09;的通信协议&#xff0c;最初由IBM在1984年开发并标准化为IEEE 802.5标准。在令牌环网络中&#xff0c;所有的计算机或工作站被连接成一个逻辑或物理的环形拓扑结构。网络中…