Naive local search

tsp_ls_naive.jl

mutable struct TSP
    name::String
    num::Int
    coord::Array{Float64, 2}
end

function read(fn::String = "./data/rat783.tsp")
    tspname = ""
    tspn = 0
    coords = nothing
    flag = false
    for line in readlines(fn)
        line = strip(line)
        if contains(line, ":")
            data = strip.(split(line, ":"))
            if data[1] == "NAME"
                tspname = data[2]
            elseif data[1] == "TYPE"
                if data[2] != "TSP"
                    exit(0)
                end
            elseif data[1] == "DIMENSION"
                tspn = parse(Int, data[2])
                coords = zeros(Int, tspn, 2)
            elseif data[1] == "EDGE_WEIGHT_TYPE"
                if data[2] != "EUC_2D"
                    exit(0)
                end
            end
        end

        if line == "NODE_COORD_SECTION"
            flag = true
            continue
        end

        if flag && line != "EOF"
            n, x, y = parse.(Int, strip.(split(line)))
            coords[n, 1] = x
            coords[n, 2] = y
        end
    end

    TSP(tspname, tspn, coords)
end


function distance(tsp::TSP, i::Int, j::Int)
    @assert 1 <= i <= tsp.num
    @assert 1 <= j <= tsp.num
    xd = tsp.coord[i, 1] - tsp.coord[j, 1]
    yd = tsp.coord[i, 2] - tsp.coord[j, 2]
    return floor(Int, sqrt(xd * xd + yd * yd) + 0.5)
end


mutable struct Work
    tour::Vector{Int}
    obj::Int
end


function calculate_objective(tsp::TSP, tour)
    length = 0.0
    for i in 1:tsp.num
        if i < tsp.num
            length += distance(tsp, tour[i], tour[i + 1])
        else
            length += distance(tsp, tour[tsp.num], tour[1])
        end
    end
    return length
end


function initialize(tsp::TSP)
    tour = collect(1:tsp.num)
    obj = calculate_objective(tsp, tour)
    Work(tour, obj)
end


function nearest_neighbor!(tsp::TSP, work)
    for i in 2:tsp.num
        mindist = Inf
        arg_mindist = nothing
        for j in i:tsp.num
            dist = distance(tsp, work.tour[i - 1], work.tour[j])
            if dist < mindist
                mindist = dist
                arg_mindist = j
            end
        end
        work.tour[i], work.tour[arg_mindist] = work.tour[arg_mindist], work.tour[i]
    end
    work.obj = calculate_objective(tsp, work.tour)
end


function local_search!(tsp::TSP, work; size=3)
    println("\n[local search algorithm]")

    count = -1
    while true
        count += 1
        # 2-opt
        two_opt_search!(tsp, work)
        println("$(count) 2-opt  $(work.obj)")
        # or-opt
        if or_opt_search!(tsp, work; size=size)
            println("$(count) Or-opt $(work.obj)")
            continue
        end
        
        # 3-opt
        if three_opt_search!(tsp, work)
            println("$(count) 3-opt  $(work.obj)")
            continue
        end
        break
    end

    println("length=$(work.obj)")
end


function two_opt_search!(tsp::TSP, work)
    function eval_diff(w, i, j)
        ni = i + 1 > tsp.num ? 1 : i + 1
        nj = j + 1 > tsp.num ? 1 : j + 1
        u, next_u = w.tour[i], w.tour[ni]
        v, next_v = w.tour[j], w.tour[nj]
        cur = distance(tsp, u, next_u) + distance(tsp, v, next_v)
        new = distance(tsp, u, v) + distance(tsp, next_u, next_v)
        return new - cur
    end

    function change_tour!(w, i, j)
        w.obj += eval_diff(w, i, j)
        w.tour[i+1:j] = reverse(w.tour[i+1:j])
    end

    improved = false
    restart = true
    while restart
        restart = false
        nbhd = ((i, j) for i in 1:length(work.tour) for j in (i+2):length(work.tour))
        for (i, j) in nbhd
            delta = eval_diff(work, i, j)
            if delta < 0
                change_tour!(work, i, j)
                improved = true
                restart = true
                break
            end
        end
    end
    return improved
end


function or_opt_search!(tsp::TSP, work; size=OR_OPT_SIZE)
    N = tsp.num

    function eval_diff(w, s, i, j)
        @assert s >= 1
        head_p = w.tour[mod1(i, N)]
        tail_p = w.tour[mod1(i + s - 1, N)]
        prev_p = w.tour[mod1(i - 1, N)]
        next_p = w.tour[mod1(i + s, N)]
        v = w.tour[mod1(j, N)]
        next_v = w.tour[mod1(j + 1, N)]
        cur = distance(tsp, prev_p, head_p) + distance(tsp, tail_p, next_p) + distance(tsp, v, next_v)
        new_f = distance(tsp, prev_p, next_p) + distance(tsp, v, head_p) + distance(tsp, tail_p, next_v)
        new_b = distance(tsp, prev_p, next_p) + distance(tsp, v, tail_p) + distance(tsp, head_p, next_v)
        if new_f <= new_b
            return new_f - cur, "fwd"
        else
            return new_b - cur, "bak"
        end
    end

    function change_tour!(w, s, i, j, oper)
        # subpath [i, ..., i + s - 1]
        subpath = []
        for h in 0:s-1
            push!(subpath, w.tour[mod1(i + h, N)])
        end
        if oper == "bak"
            reverse!(subpath)
        end

        # move subpath [i, ..., i + 1 - 1] to j + 1
        for h in i+s:j
            w.tour[mod1(N + h - s, N)] = w.tour[mod1(h, N)]
        end
        for h in 0:s-1
            w.tour[mod1(N + j + 1 - s + h, N)] = subpath[h + 1]
        end
        work.obj = calculate_objective(tsp, work.tour)
    end

    improved = false
    restart = true
    while restart
        restart = false
        nbhd = ((s, i, j) for s in 1:(size + 1) for i in 1:length(work.tour) for j in (i+s):(i+length(work.tour)-2))
        for (s, i, j) in nbhd
            delta, oper = eval_diff(work, s, i, j)
            if delta < 0
                change_tour!(work, s, i, j, oper)
                improved = true
                restart = true
                break
            end
        end
    end
    return improved
end

function three_opt_search!(tsp::TSP, work)
    N = tsp.num

    function eval_diff(work, i, j, k)
        best, arg_best = Inf, nothing
        u, next_u = work.tour[i], work.tour[mod1(i + 1, N)]
        v, next_v = work.tour[j], work.tour[mod1(j + 1, N)]
        w, next_w = work.tour[k], work.tour[mod1(k + 1, N)]
        cur = distance(tsp, u, next_u) + distance(tsp, v, next_v) + distance(tsp, w, next_w)
        new = distance(tsp, u, next_v) + distance(tsp, v, next_w) + distance(tsp, w, next_u)
        if new - cur < best
            best, arg_best = new - cur, "type1"
        end
        new = distance(tsp, u, w) + distance(tsp, next_v, next_u) + distance(tsp, v, next_w)
        if new - cur < best
            best, arg_best = new - cur, "type2"
        end
        new = distance(tsp, u, next_v) + distance(tsp, w, v) + distance(tsp, next_u, next_w)
        if new - cur < best
            best, arg_best = new - cur, "type3"
        end
        new = distance(tsp, v, u) + distance(tsp, next_w, next_v) + distance(tsp, w, next_u)
        if new - cur < best
            best, arg_best = new - cur, "type4"
        end
        return best, arg_best
    end

    function change_tour!(w, i, j, k, oper)
        @assert i < j
        @assert j < k
        subpath = []
        if oper == "type1"
            for z in (j + 1):k
                push!(subpath, w.tour[z])
            end
            for z in (i + 1):j
                push!(subpath, w.tour[z])
            end
        elseif oper == "type2"
            for z in range(k, j + 1; step=-1)
                push!(subpath, w.tour[z])
            end
            for z in (i + 1):j
                push!(subpath, w.tour[z])
            end
        elseif oper == "type3"
            for z in (j + 1):k
                push!(subpath, w.tour[z])
            end
            for z in range(j, i + 1; step=-1)
                push!(subpath, w.tour[z])
            end
        elseif oper == "type4"
            for z in range(j, i + 1; step=-1)
                push!(subpath, w.tour[z])
            end
            for z in range(k, j + 1; step=-1)
                push!(subpath, w.tour[z])
            end
        end
        w.tour[i + 1:k] .= subpath
        work.obj = calculate_objective(tsp, work.tour)
    end

    improved = false
    restart = true
    
    while restart
        restart = false
        nbhd = ((i, j, k) 
                    for i in 1:N
                    for j in (i + 2):N
                    for k in (j + 2):N)
        for (idx, (i, j, k)) in enumerate(nbhd)
            delta, oper = eval_diff(work, i, j, k)
            if delta < 0
                change_tour!(work, i, j, k, oper)
                improved = true
                restart = true
                break
            end
        end
    end
    return improved
end



function main()
    OR_OPT_SIZE = 3
    start_time = time()
    # tsp = read("./data/rat783.tsp")
    tsp = read("./data/a280.tsp")
    work = initialize(tsp)
    println("seq length=$(work.obj)")
    nearest_neighbor!(tsp, work)
    local_search!(tsp, work; size=OR_OPT_SIZE)
    end_time = time()

    println("\nTotal time:\t$(end_time - start_time) sec")
end