No description has been provided for this image

Programming large-scale parallel systems¶

Traveling sales person¶

Contents¶

In this notebook, we will learn

  • How to parallelize the solution of the traveling sales person problem
  • The concept of search overhead

The traveling sales person (TSP) problem¶

Let us start by presenting the traveling sales person (ASP) problem and its solution using a branch and bound algorithm.

Problem statement¶

  • Given a graph $G$ with a distance table $C$ and an initial node (i.e. a city) in the graph
  • Compute the shortest route that visits all cities exactly once, without returning to the initial city.

Note that there we consider a version of the problem in which do not return to the initial city. However, the classic version of the problem includes returning to the initial city.

As for the ASP problem we represent the distance table as a matrix, where $C_{ij}$ is the distance from node $i$ to node $j$. Next figure shows the input and solution (output) of the TSP problem for a simple 4-node graph.

No description has been provided for this image

Branch and bound algorithm¶

We consider a branch and bound strategy to solve the ASP problem. The first ingredient of this strategy is to enumerate all potential solutions of the problem using a search tree. The tree is defined as follows:

  • The root is the initial city
  • The children of a node are the neighbor cities not visited so-far
  • We sort the children using the nearest city first heuristic
  • A node will be a leaf of the tree when all neighbor cities are already visited

The nearest city first heuristic is a way to speed up the search as we increase the possibility of finding the solution in the first paths enumerated in the tree.

No description has been provided for this image

Pruning the search tree¶

We look for a solution by going over the branches of the search tree via first-depth search. Each time we hit a lead in the tree we have a route between all cities and this a possible solution of the TSP problem (called a partial solution). In this process, we keep track of which is the shortest partial solution found so far. Note that we don't need to completely traverse all branches. In some cases, the distance at an intermediate point in a route will be already equal or greater than the best distance found so-far. In this case, we don't need to continue moving within this branch and we can jump to the next one. This is the "bound" criterion in our branch and bound strategy, which is also referred to as "pruning" the search tree.

In the example below, the first route has distance 6. At an intermediate point of the 2nd route we have already covered a distance equal to 6 which guarantees that the solution will not be in this part of the tree.

No description has been provided for this image

Thus, there are only 2 more complete routes that we need to visit. Using this strategy we have significantly reduced the routes visited: from 6 possible routes to 3 visited routes.

No description has been provided for this image

Serial implementation¶

Nearest-city first heuristic¶

In [ ]:
function sort_neighbors(C)
    n = size(C,1)
    map(1:n) do i
        Ci = C[i,:]
        cities = sortperm(Ci)
        distances = Ci[cities]
        collect(zip(cities,distances))
    end
end
In [ ]:
C = [
    0 2 3 2
    2 0 4 1
    3 4 0 3
    2 1 3 0    
]
In [ ]:
C_sorted = sort_neighbors(C)
In [ ]:
city = 3
C_sorted[city]
No description has been provided for this image

Loop over all paths¶

In [ ]:
function visital_all_paths(C_sorted,city)
    num_cities = length(C_sorted)
    path=zeros(Int,num_cities)
    hops = 1
    path[hops] = city
    visital_all_paths_recursive!(C_sorted,hops,path)
end
function visital_all_paths_recursive!(C_sorted,hops,path)
    num_cities = length(C_sorted)
    if hops != num_cities
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            if ! already_visited
                path[next_hops] = next_city
                visital_all_paths_recursive!(C_sorted,next_hops,path)
            end
        end
        return nothing
    else
        @show path
        return nothing
    end
end
In [ ]:
city = 1
visital_all_paths(C_sorted,city)
No description has been provided for this image

Serial implementation without pruning¶

In [ ]:
function tsp_serial_no_prune(C_sorted,city)
    num_cities = length(C_sorted)
    path=zeros(Int,num_cities)
    hops = 1
    path[hops] = city
    distance=0
    min_distance = typemax(Int)
    tsp_serial_no_prune_recursive!(C_sorted,hops,path,distance,min_distance)
end
function tsp_serial_no_prune_recursive!(C_sorted,hops,path,distance,min_distance)
    num_cities = length(C_sorted)
    if hops != num_cities
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            if ! already_visited
                next_distance = distance + distance_increment
                path[next_hops] = next_city
                min_distance = tsp_serial_no_prune_recursive!(
                    C_sorted,next_hops,path,next_distance,min_distance)
            end
        end
        return min_distance
    else
        min_distance = min(distance,min_distance)
        #@show path, distance, min_distance
        return min_distance
    end
end
In [ ]:
city = 1
min_distance = tsp_serial_no_prune(C_sorted,city)

Final serial implementation¶

In [ ]:
function tsp_serial(C_sorted,city)
    num_cities = length(C_sorted)
    path=zeros(Int,num_cities)
    hops = 1
    path[hops] = city
    distance=0
    min_distance = typemax(Int)
    tsp_serial_recursive!(C_sorted,hops,path,distance,min_distance)
end
function tsp_serial_recursive!(C_sorted,hops,path,distance,min_distance)
    if distance >= min_distance
        return min_distance
    end
    num_cities = length(C_sorted)
    if hops != num_cities
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            if ! already_visited
                next_distance = distance + distance_increment
                path[next_hops] = next_city
                min_distance = tsp_serial_recursive!(
                    C_sorted,next_hops,path,next_distance,min_distance)
            end
        end
        return min_distance
    else
        min_distance = min(distance,min_distance)
        #@show path, distance, min_distance
        return min_distance
    end
end
In [ ]:
city = 1
min_distance = tsp_serial(C_sorted,city)

Performance¶

In [ ]:
n = 11
using Random
Random.seed!(1)
C = rand(1:10,n,n)
C_sorted = sort_neighbors(C)
city = 1
@time tsp_serial_no_prune(C_sorted,city)
@time tsp_serial(C_sorted,city)

Parallel implementation¶

In [ ]:
using Distributed
In [ ]:
if workers() == procs()
    addprocs(3)
end
In [ ]:
function enumerate_paths_dist(C_sorted,city,max_hops)
    T = typeof((0,Int[]))
    jobs_chnl = RemoteChannel(()->Channel{T}(1))
    @sync begin
        for w in workers()
            @spawnat w consume_jobs(C_sorted,jobs_chnl)
        end
        generate_jobs(C_sorted,city,max_hops,jobs_chnl)
    end
end
function generate_jobs(C_sorted,city,max_hops,jobs_chnl)
    num_cities = length(C_sorted)
    path=zeros(Int,num_cities)
    hops = 1
    path[hops] = city
    generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl)
    for w in workers()
        put!(jobs_chnl,(0,Int[]))
    end
    close(jobs_chnl)
end
function generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl)
    num_cities = length(C_sorted)
    if hops == max_hops
        @show path
        put!(jobs_chnl,(hops,copy(path)))
        return nothing
    else
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            if ! already_visited
                path[next_hops] = next_city
                generate_jobs_recursive(C_sorted,next_hops,path,max_hops,jobs_chnl)
            end
        end
        return nothing
    end
end
@everywhere function consume_jobs(C_sorted,jobs_chnl)
    while true
        hops,path = take!(jobs_chnl)
        if hops == 0
            println("Done!")
            break
        end
        consume_jobs_recursive(C_sorted,hops,path)
    end
end
@everywhere function consume_jobs_recursive(C_sorted,hops,path)
    num_cities = length(C_sorted)
    if hops != num_cities
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            if ! already_visited
                path[next_hops] = next_city
                consume_jobs_recursive(C_sorted,next_hops,path)
            end
        end
        return nothing
    else
        @show path
        return nothing
    end
end
In [ ]:
C = [
    0 2 3 2
    2 0 4 1
    3 4 0 3
    2 1 3 0    
]
C_sorted = sort_neighbors(C)
city = 1
max_hops = 2
min_distance = enumerate_paths_dist(C_sorted,city,max_hops)

How to track the global minimum distance?¶

In [ ]:
buffer = 1 # very important
min_distance_chnl = RemoteChannel(()->Channel{Int}(buffer))
put!(min_distance_chnl,typemax(Int))
@sync for w in workers()
    @spawnat w begin
         sleep(rand(1:3))
         min_distance = take!(min_distance_chnl)
         @show min_distance
         distance = rand(5:10)
         min_distance = min(distance,min_distance)
         @show distance
         put!(min_distance_chnl,min_distance)
    end
end
min_distance = take!(min_distance_chnl)
@show min_distance
close(min_distance_chnl)

Final parallel implementation¶

In [ ]:
function tsp_dist(C_sorted,city,max_hops)
    T = typeof((0,Int[],0))
    jobs_chnl = RemoteChannel(()->Channel{T}(1))
    min_distance_chnl = RemoteChannel(()->Channel{Int}(1))
    put!(min_distance_chnl,typemax(Int))
    @sync begin
        for w in workers()
            @spawnat w consume_jobs(C_sorted,jobs_chnl,min_distance_chnl)
        end
        generate_jobs(C_sorted,city,max_hops,jobs_chnl)
    end
    min_distance = take!(min_distance_chnl)
    close(min_distance_chnl)
    return min_distance
end
function generate_jobs(C_sorted,city,max_hops,jobs_chnl)
    num_cities = length(C_sorted)
    path=zeros(Int,num_cities)
    hops = 1
    path[hops] = city
    distance = 0
    generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl,distance)
    for w in workers()
        put!(jobs_chnl,(0,Int[],0))
    end
    close(jobs_chnl)
end
function generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl,distance)
    num_cities = length(C_sorted)
    if hops == max_hops
        @show path, distance
        put!(jobs_chnl,(hops,copy(path),distance))
        return nothing
    else
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            next_distance = distance + distance_increment
            if ! already_visited
                path[next_hops] = next_city
                generate_jobs_recursive(
                        C_sorted,next_hops,path,max_hops,jobs_chnl,next_distance)
            end
        end
        return nothing
    end
end
@everywhere function consume_jobs(C_sorted,jobs_chnl,min_distance_chnl)
    while true
        hops,path,distance = take!(jobs_chnl)
        if hops == 0
            println("Done!")
            break
        end
        min_distance = take!(min_distance_chnl)
        put!(min_distance_chnl,min_distance)
        consume_jobs_recursive(C_sorted,hops,path,min_distance_chnl,distance,min_distance)
    end
end
@everywhere function consume_jobs_recursive(C_sorted,hops,path,min_distance_chnl,distance,min_distance)
    if distance >= min_distance
        return min_distance
    end
    num_cities = length(C_sorted)
    if hops != num_cities
        city = path[hops]
        connections = C_sorted[city]
        next_hops = hops + 1
        for (next_city,distance_increment) in connections
            already_visited = (next_city in view(path,1:hops))
            next_distance = distance + distance_increment
            if ! already_visited
                path[next_hops] = next_city
                min_distance = consume_jobs_recursive(
                        C_sorted,next_hops,path,min_distance_chnl,next_distance,min_distance)
            end
        end
        return min_distance
    else
        min_distance = take!(min_distance_chnl)
        min_distance = min(min_distance,distance)
        put!(min_distance_chnl,min_distance)
        @show path, distance, min_distance
        return min_distance
    end
end
In [ ]:
city = 1
max_hops = 2
min_distance = tsp_dist(C_sorted,city,max_hops)