Implement k-opt.

tsp/data/tsp_lec

@@ -0,0 +1,13 @@
+12
+0.0 0.0
+0.6 4.2
+2.7 3.7
+2.1 4.3
+7.2 4.6
+8.3 3.2
+7.0 0.5
+4.7 2.8
+5.2 3.6
+3.7 3.7
+3.6 2.8
+2.1 0.0

tsp/tsp.py

@@ -3,18 +3,21 @@ from functools import lru_cache
 from collections import namedtuple
 from geometry import intersect
 
-Point = namedtuple("Point", ['index', 'x', 'y'])
+Point = namedtuple("P", ['name', 'x', 'y'])
+DEBUG = False
 
 
 def parse_input_data(input_data):
     lines = input_data.split('\n')
     node_count = int(lines[0])
-    return [Point(i, *map(float, lines[i + 1].split()))
+    return [Point(str(i), *map(float, lines[i + 1].split()))
             for i in range(0, node_count)]
 
 
 def plot_graph(points):
     import matplotlib.pyplot as plt
+    if not DEBUG:
+        return
 
     def plot_arrows():
         for i in range(len(points)):
@@ -34,144 +37,228 @@ def plot_graph(points):
     def plot_points():
         for p in points:
             plt.plot(p.x, p.y, '')
-            plt.text(p.x, p.y, '  ' + str(p.index))
+            plt.text(p.x, p.y, '  ' + p.name)
 
     plot_points()
     plot_arrows()
     plt.show()
 
 
-def solve_it(input_data):
+def prepare_output_data(points):
+    # Basic plausibility checks
+    assert(len(set(points)) == len(points))
+    assert(len(points) > 4)
+    obj = total_distance(points)
+    output_data = '%.2f' % obj + ' ' + str(0) + '\n'
+    output_data += ' '.join(map(lambda p: p.name, points))
+    return output_data
 
-    @lru_cache(maxsize=1000000)
-    def length(p_1, p_2):
-        return math.sqrt((p_1.x - p_2.x)**2 + (p_1.y - p_2.y)**2)
 
-    def prepare_output_data(points):
+@lru_cache(maxsize=1000000)
-        obj = total_length(points)
+def distance(point_1, point_2):
-        output_data = '%.2f' % obj + ' ' + str(0) + '\n'
+    """ Calculate the distance between two points. """
-        output_data += ' '.join(map(lambda p: str(p.index), points))
+    p1, p2 = point_1, point_2
-        return output_data
+    return math.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)
 
-    def is_valid(points, num_nodes):
-        expected_points = set(range(num_nodes))
-        assert(set([p.index for p in points]) == expected_points)
-        return True
 
-    def total_length(points):
+def total_distance(points):
-        return sum([length(points[i - 1], points[i])
+    """ Calculate the total distance of the point sequence. """
-                    for i in range(len(points))])
+    # Use negative indexing to get the distance from last to first point
+    return sum([distance(points[i - 1], points[i])
+                for i in range(len(points))])
 
-    def initial_solution_naiv():
+
+def longest_distance(points, ignore_list):
+    """ Returns the point and index of the
+    point with the longest distance to the next point. """
+    longest_distance = 0
+    longest_dist_point = None
+    longest_dist_index = None
+    for i in range(len(points)):
+        p1, p2 = points[i - 1], points[i]
+        if p1 in ignore_list:
+            continue
+        current_distance = distance(p1, p2)
+        if current_distance > longest_distance:
+            longest_distance = current_distance
+            longest_dist_point = p1
+            longest_dist_index = i - 1
+    return longest_dist_point, longest_dist_index
+
+
+def swap_edges(i, j, points):
+    """
+    Swaps edges in-place. Also returns result.
+
+    :param i: Index of first point of first edge.
+    :param j: Index if first point of second edge.
+    """
+    assert(i != j)
+    _, p12 = points[i], points[i + 1]
+    p21, _ = points[j], points[j + 1]
+    points[i + 1] = p21
+    points[j] = p12
+
+    # Reverse order of points between swapped lines.
+    if i < j:
+        points[i + 2:j] = points[i + 2:j][::-1]
+    else:
+        # List goes over boundaries
+        len_points = len(points)
+        segment = points[i + 2:] + points[:j]
+        segment.reverse()
+        points[i + 2:] = segment[:len_points - i - 2]
+        points[:j] = segment[len_points - i - 2:]
+    return points
+
+
+def local_search(points, ignore_list):
+    #print("-" * 80)
+    #print("Local search")
+    #print("ignore_list", ignore_list)
+
+    max_len = 0
+    max_index = None
+    for i in range(len(points)):
+        if points[i - 1] in ignore_list:
+            continue
+        new_len = length(points[i - 1], points[i])
+        if new_len > max_len:
+            max_len = new_len
+            p_i = i - 1
+    p1 = points[p_i]
+    p2 = points[p_i + 1]
+    #print("Found max_len for ", edge(p1, p2))
+
+    current_length = total_distance(points)
+    for p_j in range(len(points)):
+        if p_j in [p_i, p_i + 1, p_i + 2]:
+            continue
+        q1 = points[p_j - 1]
+        q2 = points[p_j]
+
+        new_points = list(points)
+        swap_edges(p_i, p_j - 1, new_points)
+        new_length = total_distance(new_points)
+        if new_length < current_length:
+            #print("Swaping", edge(points[p_i], points[p_i + 1]), "and", edge(points[p_j - 1], points[p_j]))
+            #print("Better new_points", new_length, "smaller", current_length)
+            ignore_list.clear()
+            return new_points
+
+    #print("Did not find an intersection that provides better results.")
+    ignore_list.append(p1)
+    return points
+
+
+def reorder_points_greedy(points):
+    current_point = points[0]
+    solution = [current_point]
+    points = points[1:]
+
+    while points:
+        min_length = 999999
+        min_point = None
+        for next_point in points:
+            new_length = distance(current_point, next_point)
+            if new_length < min_length:
+                min_length = new_length
+                min_point = next_point
+        current_point = min_point
+        solution.append(current_point)
+        points.remove(current_point)
+
+    return solution
+
+
+def print_swap(i, j, points):
+    if not DEBUG:
+        return
+    print("Swap:", points[i].name, " <-> ", points[j].name)
+
+
+def k_opt(p1_index, points, ignore_list, swaps):
+    print("k_opt ignore_list len", len(ignore_list))
+    i = p1_index
+    p1, p2 = points[i], points[i + 1]
+    dist_p1p2 = distance(p1, p2)
+    ignore_list.append(p2)
+
+    p4_index = None
+    for p3_index in range(len(points)):
+        p3 = points[p3_index]
+        p4 = points[p3_index - 1]
+        if p4 in ignore_list or p4 is p1:
+            continue
+        dist_p2p3 = distance(p2, p3)
+        if dist_p2p3 < dist_p1p2:
+            p4_index = p3_index - 1
+            dist_p1p2 = dist_p2p3
+
+    if not p4_index:
         return []
 
-    def reorder_points_greedy(points):
+    print_swap(p1_index, p4_index, points)
-        current_point = points[0]
+    plot_graph(points)
-        solution = [current_point]
+    swap_edges(p1_index, p4_index, points)
-        points = set(points[1:])
+    swaps.append([p1_index, p4_index])
-
+    new_total = total_distance(points)
-        while points:
+    print("Current distance", new_total)
-            min_length = 999999
+    r = k_opt(p1_index, points, ignore_list, list(swaps))
-            min_point = None
+    r.append((new_total, swaps))
-            for next_point in points:
+    return r
-                new_length = length(current_point, next_point)
-                if new_length < min_length:
-                    min_length = new_length
-                    min_point = next_point
-            current_point = min_point
-            solution.append(current_point)
-            points.remove(current_point)
-
-        return solution
-
-    def swap_edges(i, j, points):
-        _, p12 = points[i], points[i + 1]
-        p21, _ = points[j], points[j + 1]
-
-        points[i + 1] = p21
-        points[j] = p12
-        points[i + 2:j] = points[i + 2:j][::-1]
-        return points
-
-    def edge(p1, p2):
-        return "{} -> {}".format(p1.index, p2.index)
-
-    def local_search(points, ignore_list):
-        # Find longest edges to swap
-        #print("-" * 80)
-        #print("Local search")
-        #print("ignore_list", ignore_list)
-
-        max_len = 0
-        max_index = None
-        for i in range(len(points)):
-            if points[i - 1] in ignore_list:
-                continue
-            new_len = length(points[i - 1], points[i])
-            if new_len > max_len:
-                max_len = new_len
-                p_i = i - 1
-        p1 = points[p_i]
-        p2 = points[p_i + 1]
-        #print("Found max_len for ", edge(p1, p2))
-
-        current_length = total_length(points)
-        for p_j in range(len(points)):
-            if p_j in [p_i, p_i + 1, p_i + 2]:
-                continue
-            q1 = points[p_j - 1]
-            q2 = points[p_j]
-
-            new_points = list(points)
-            swap_edges(p_i, p_j - 1, new_points)
-            new_length = total_length(new_points)
-            if new_length < current_length:
-                #print("Swaping", edge(points[p_i], points[p_i + 1]),
-                #      "and", edge(points[p_j - 1], points[p_j]))
-                #print("Test new_points", new_length, "smaller", current_length)
-                #print("Better, return new_points.")
-                ignore_list.clear()
-                return new_points
-
-            #if intersect(p1, p2, q1, q2):
-            #    print(edge(p1, p2), "intersects", edge(q1, q2))
-            #    new_points = list(points)
-            #    swap_edges(p_i, p_j - 1, new_points)
-            #    new_length = total_length(new_points)
-            #    print("Test new_points", new_length, "smaller", current_length)
-            #    if new_length < current_length:
-            #        print("Better, return new_points.")
-            #        ignore_list.append(p1)
-            #        return new_points
-            #    else:
-            #        print("Worse, find better intersection.")
-
-        #print("Did not find an intersection that provides better results.")
-        ignore_list.append(p1)
-        return points
-
-    points = parse_input_data(input_data)
-    points = reorder_points_greedy(points)
-    num_nodes = len(points)
 
 
+def local_search_k_opt(points):
+    current_total = total_distance(points)
     ignore_list = []
     while True:
-        points_before_change = list(points)
+        print()
-        try:
+        print("--- new iteration ---")
-            points = local_search(points, ignore_list)
+        print("Ignored points", [p.name for p in ignore_list])
-        except UnboundLocalError:
+        point, index = longest_distance(points, ignore_list)
+        if not point:
+            print("No more points")
             break
-        is_valid(points, num_nodes)
-        value = total_length(points)
 
-    # plot_graph(points_before_change)
+        ignore_list.append(point)
-    is_valid(points, num_nodes)
+        print("Next point (longest_distance)", point)
+
+        r = k_opt(index, list(points), [], [])
+        print("k-opt", len(r))
+        if not r:
+            print("Found no better solution.")
+            continue
+
+        new_total, steps = min(r)
+
+        print("new_total", new_total, "current_total", current_total)
+        if new_total < current_total:
+            print("Improvment. Apply steps.")
+            for step in steps:
+                swap_edges(*step, points)
+            assert(total_distance(points) == new_total)
+            current_total = new_total
+            ignore_list = []
+        else:
+            print("No changes.")
+        plot_graph(points)
+    return points
+
+
+def solve_it(input_data):
+    points = parse_input_data(input_data)
+    num_points = len(points)
+    #points = reorder_points_greedy(points)
+    local_search_k_opt(points)
+
+    # plot_graph(points)
     return prepare_output_data(points)
 
 
 if __name__ == "__main__":
-    file_location = "data/tsp_200_2"
+    file_location = "data/tsp_51_1"
+    # DEBUG = True
     with open(file_location, 'r') as input_data_file:
         input_data = input_data_file.read()
         print(solve_it(input_data))