2-opt algorithm to solve the Travelling Salesman Problem in Python

Question

I couldn't find any complete implementation of the 2-opt algorithm in Python so I am trying to add the missing parts to the code found here, which I present below.

def two_opt(route):
     best = route
     improved = True
     while improved:
          improved = False
          for i in range(1, len(route)-2):
               for j in range(i+1, len(route)):
                    if j-i == 1: continue # changes nothing, skip then
                    new_route = route[:]
                    new_route[i:j] = route[j-1:i-1:-1] # this is the 2woptSwap
                    if cost(new_route) < cost(best):  # what should cost be?
                         best = new_route
                         improved = True
          route = best
     return best

In order to complete this code, I made a small program to extract long/lat co-ords from a text file and fill in an adjacency matrix with the cost for each point. Full code, including samples of input co-ordinates and adjacency matrix may be found on Code Review.

Since I do not know what the cost function is from the code above, my idea was to work out all the costs from one point to another and placed in an adjacency matrix: adj_matrix. This represents how far each point is from the others.

I tried passing my cost/adjacency matrix to the function to use that, however I am unable to calculate the cost given my adjacency matrix.

def main():
    # code to read from file
    # code to append co-ordinates to points and calculate the haversine distance between each point
    route = random.sample(range(10), 10)
    best = two_opt(route, adj_matrix)  # passing my adjacency matrix
    print(best)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Another Python 2-opt question: Generate all neighbors for 2OPT in python

Any suggestions on how I can find the correct cost from the adjacency matrix would be appreciated.

Answer 1

First of all, an adjacency matrix is typically a (0, 1)-matrix. What you have here is variously referred to as the cost, weight or distance matrix.

Now to your question.

The cost function can be as simple as:

def cost(cost_mat, route):
   return cost_mat[np.roll(route, 1), route].sum()

Here, np.roll() "rotates" the route by one position to make it easy to use it with route to index into the cost matrix. The sum() simply adds up the costs of the individual segment to compute the total cost of the route.

(If at some point you decide to look at the Asymmetric TSP, you'll need to make sure the row/column order matches how cost_mat is constructed; for the Euclidean TSP this doesn't matter as the cost matrix is symmetric.)

Example of use:

cost_mat = np.array([
   [0, 1, 2, 3],
   [1, 0, 4, 5],
   [2, 4, 0, 7],
   [3, 5, 7, 0],
])

route = np.array([2, 1, 3, 0])

print(cost(cost_mat, route))

Answer 2

2-opt deletes two edges and creates two new ones (assuming the cost matrix is symmetric), so the cost function can be simplified to consider only the edges that change. For large arrays this is much faster than enumerating over the whole route.

import numpy as np

def cost_change(cost_mat, n1, n2, n3, n4):
    return cost_mat[n1][n3] + cost_mat[n2][n4] - cost_mat[n1][n2] - cost_mat[n3][n4]


def two_opt(route, cost_mat):
    best = route
    improved = True
    while improved:
        improved = False
        for i in range(1, len(route) - 2):
            for j in range(i + 1, len(route)):
                if j - i == 1: continue
                if cost_change(cost_mat, best[i - 1], best[i], best[j - 1], best[j]) < 0:
                    best[i:j] = best[j - 1:i - 1:-1]
                    improved = True
        route = best
    return best


if __name__ == '__main__':
    nodes = 1000
    init_route = list(range(nodes))
    print(init_route)
    cost_mat = np.random.randint(100, size=(nodes, nodes))
    cost_mat += cost_mat.T
    np.fill_diagonal(cost_mat, 0)
    cost_mat = list(cost_mat)
    best_route = two_opt(init_route, cost_mat)
    print(best_route)