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

news/2024/7/16 6:15:22 标签: julia, 算法, 开发语言

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

http://www.niftyadmin.cn/n/5542847.html

相关文章

并发编程工具集——读写锁-ReadWriteLock（上篇）（十六）

并发编程工具集——读写锁-ReadWriteLock（上篇）（十六）

什么是读写锁 基本原则&#xff1a;&#xff08;读读不互斥、读写互斥、写写互斥&#xff09; 允许多个线程同时读共享变量&#xff1b;只允许一个线程写共享变量&#xff1b;如果一个写线程正在执行写操作&#xff0c;此时禁止读线程读共享变量。读写锁与互斥锁的一个重要区别…
阅读更多...
唤醒知识循环，共筑绿色阅读梦——探索旧书回收小程序的无限可能

唤醒知识循环，共筑绿色阅读梦——探索旧书回收小程序的无限可能

在这个信息爆炸的时代&#xff0c;书籍作为知识与智慧的载体&#xff0c;其重要性不言而喻。然而&#xff0c;随着电子阅读的兴起和书籍更新换代的加速&#xff0c;大量旧书被束之高阁&#xff0c;甚至面临被遗弃的命运。这不仅是对宝贵文化资源的浪费&#xff0c;也是对环境保…
阅读更多...
Webpack: Loader开发 (1)

Webpack: Loader开发 (1)

概述 如何扩展 Webpack&#xff1f;有两种主流方式&#xff0c;一是 Loader —— 主要负责将资源内容翻译成 Webpack 能够理解、处理的 JavaScript 代码&#xff1b;二是 Plugin —— 深度介入 Webpack 构建过程&#xff0c;重塑 构建逻辑。 相对而言&#xff0c;Loader 的职责…
阅读更多...
最新整理的机器人相关数据合集（1993-2022年不等 具体看数据类型）

最新整理的机器人相关数据合集（1993-2022年不等 具体看数据类型）

机器人安装数据是指记录全球或特定区域内工业机器人新安装数量的信息&#xff0c;这一数据由国际机器人联合会(IFR)等权威机构定期发布。这些数据不仅揭示了机器人技术的市场需求趋势&#xff0c;还反映了各国和地区自动化水平及产业升级的步伐。例如&#xff0c;数据显示中国在…
阅读更多...
【MySQL基础篇】函数及约束

【MySQL基础篇】函数及约束

1、函数 函数是指一段可以直接被另一段程序程序调用的程序或代码。 函数 - 字符串函数 MySQL中内置了很多字符串函数&#xff0c;常用的几个如下&#xff1a; 函数功能CONCAT(S1,S2,...,Sn)字符串拼接&#xff0c;将S1,S2,...,Sn拼接成一个字符串LOWER(str)将字符串str全部…
阅读更多...
人工智能--循环神经网络

人工智能--循环神经网络

个人主页&#xff1a;欢迎来到 Papicatch的博客 课设专栏 &#xff1a;学生成绩管理系统 专业知识专栏&#xff1a; 专业知识 文章目录 &#x1f349;引言 &#x1f349;概述 &#x1f348;基本概念 &#x1f34d;定义 &#x1f34d;结构 &#x1f34c;输入层 &#…
阅读更多...
解释如何在使用Bitmap时进行优化，以减少内存占用和提高性能。

解释如何在使用Bitmap时进行优化，以减少内存占用和提高性能。

在使用Android开发中的Bitmap时&#xff0c;优化其使用以减少内存占用和提高性能是一个重要且复杂的任务。Bitmap作为图像处理的核心&#xff0c;其处理不当往往会导致内存溢出&#xff08;OutOfMemoryError&#xff09;或应用性能下降。下面从技术难点、面试官关注点、回答吸引…
阅读更多...
【Unity小技巧】Unity字典序列化

【Unity小技巧】Unity字典序列化

字典序列化 在 Unity 中&#xff0c;标准的 C# 字典&#xff08;Dictionary<TKey, TValue>&#xff09;是不能直接序列化的&#xff0c;因为 Unity 的序列化系统不支持非 Unity 序列化的集合类型。可以通过手写字典实现 效果&#xff1a; 实现步骤&#xff1a; 继承ISe…
阅读更多...
最新文章

Copyright @ 2022~2023