Simple implementation (I don't remember where I found it, but it was in C/C++) of Lin-Kernighan heuristic for TSP.

51_1.tsp

51
27 68
30 48
43 67
58 48
58 27
37 69
38 46
46 10
61 33
62 63
63 69
32 22
45 35
59 15
5 6
10 17
21 10
5 64
30 15
39 10
32 39
25 32
25 55
48 28
56 37
30 40
37 52
49 49
52 64
20 26
40 30
21 47
17 63
31 62
52 33
51 21
42 41
31 32
5 25
12 42
36 16
52 41
27 23
17 33
13 13
57 58
62 42
42 57
16 57
8 52
7 38

lk.py

import os
import signal

import numpy as np
from scipy.spatial.distance import cdist


class ObjectiveFunction:

    def __init__(self, points, use_precomputed=True):
        self.points = points
        precomputed_file = '{}.npy'.format(points.shape[0])
        if use_precomputed:
            if os.path.isfile(precomputed_file):
                print('Loading distance matrix from pre-computed file.')
                self._distances = np.load(precomputed_file)
            else:
                print('Pre-computing distance matrix from scratch.')
                self._distances = cdist(points, points, metric='euclidean')
                print('Saving distance matrix for later use.')
                np.save(precomputed_file, self._distances)
        else:
            self._distances = cdist(points, points, metric='euclidean')

    @property
    def node_count(self):
        return self._distances.shape[0]

    def __getitem__(self, from_to):
        (x, y) = from_to
        return self._distances[x, y]

    def __call__(self, s):
        start, end = s[0], s[-1]
        is_feasible = len(set(s)) == len(s) and len(s) == self.node_count
        distance = self._distances[end, start] + sum(self._distances[x, y] for x, y in zip(s, s[1:]))
        return distance, is_feasible

class GracefulKiller:
    kill_now = False

    def __init__(self):
        signal.signal(signal.SIGINT, self.exit_gracefully)
        signal.signal(signal.SIGTERM, self.exit_gracefully)

    def exit_gracefully(self, signum, frame):
        self.kill_now = True

def random_solution(f):
    points = range(f.node_count)
    return list(np.random.permutation(points))

def lin_kernighan(f, **kwargs):

    sorted_pair = lambda x, y: (x, y) if x < y else (y, x)

    def reverse(tour, start, end):
        current = start
        next = tour[start]
        while True:
            nextNext = tour[next]
            tour[next] = current
            current = next
            next = nextNext
            if current == end:
                break

    def is_tour(t):
        count = 1
        start = t[0]
        while start != 0:
            start = t[start]
            count += 1
        return count == len(t)

    def tour_ids(t):
        current = 0
        tour = []
        while True:
            tour.append(ids[current])
            current = t[current]
            if current == 0:
                break
        return tour

    def lkmove(tour, tour_start):
        broken_set = set()
        joined_set = set()
        tour_opt = tour.copy()
        g_opt = 0
        g = 0
        last_next_v = tour_start
        last_possible_next_v = 0
        from_v = tour[last_next_v]
        while True:
            next_v = -1
            broken_edge = sorted_pair(last_next_v, from_v)
            broken_edge_length = f[broken_edge]
            print('Breaking {} {}'.format(last_next_v, from_v))
            if broken_edge in joined_set:
                break

            # emulate for-cycle
            possible_next_v = tour[from_v]
            while next_v == - 1 and possible_next_v != tour_start:
                g_local = broken_edge_length - f[from_v, possible_next_v]
                print('Testing {}'.format(possible_next_v))
                if not (sorted_pair(from_v, possible_next_v) not in broken_set and g + g_local > 0 and sorted_pair(last_possible_next_v, possible_next_v) not in joined_set and tour[possible_next_v] != 0 and possible_next_v != tour[from_v]):
                    last_possible_next_v = possible_next_v
                    # iteration step
                    possible_next_v = tour[possible_next_v]
                    continue
                next_v = possible_next_v
                print('Moving to {}'.format(next_v))
                # iteration step
                possible_next_v = tour[possible_next_v]
            if next_v != -1:
                broken_set.add(broken_edge)
                joined_set.add(sorted_pair(from_v, next_v))
                y_opt_length = f[from_v, tour_start]
                g_opt_local = g + (broken_edge_length - y_opt_length)
                if g_opt_local > g_opt:
                    g_opt = g_opt_local
                    tour_opt = tour.copy()
                    tour_opt[tour_start] = from_v
                print('Joining of distance: {}'.format(f[from_v, next_v]))
                g += broken_edge_length - f[from_v, next_v]
                reverse(tour, from_v, last_possible_next_v)
                next_from = last_possible_next_v
                tour[from_v] = next_v
                last_next_v = next_v
                from_v = next_from

            # end of do-while
            if next_v == -1:
                break
        tour = tour_opt
        assert is_tour(tour)
        return tour

    t = random_solution(f)
    kill_gracefully = kwargs.get('kill_gracefully', True)
    killer = GracefulKiller() if kill_gracefully else None

    # don't know what going on here - but for some reason I wrote it like this because
    # they were using linked list represented as array.
    t = [(i + 1) % f.node_count for i in range(f.node_count)]
    assert is_tour(t)
    ids = list(range(f.node_count))
    old_distance = 0
    for j in range(100):
        for i in range(f.node_count):
            t = lkmove(t, i)
            s = tour_ids(t)
            distance, is_feasible = f(s)
            print('Iteration: {}, Solution: {:.2f}, {}'.format(j, distance, is_feasible))
            if kill_gracefully and killer.kill_now:
                break
        diff = old_distance - distance
        if j != 0 and diff == 0:
            print('Converged after {} iterations'.format(j))
            break
        old_distance = distance
        if kill_gracefully and killer.kill_now:
            break

    # Wrap into Solution class
    return tour_ids(t)

if __name__ == '__main__':
    with open('tsp_51_1', 'r') as input_data_file:
        input_data = input_data_file.read()
    lines = input_data.split('\n')
    node_count = int(lines[0])
    X = np.array([[float(p) for p in l.split(' ')] for l in lines[1:node_count + 1]])
    f = ObjectiveFunction(X, use_precomputed=False)
    lin_kernighan(f)