Refactor for better TSP algorithm.

discrete_optimization.sublime-project

@@ -1,17 +0,0 @@
-{
-	"build_systems":
-	[
-		{
-			"file_regex": "^[ ]*File \"(...*?)\", line ([0-9]*)",
-			"name": "Anaconda Python Builder",
-			"selector": "source.python",
-			"shell_cmd": "\"python\" -u \"$file\""
-		}
-	],
-	"folders":
-	[
-		{
-			"path": "."
-		}
-	]
-}

tsp/geometry.py

@@ -1,81 +0,0 @@
-
-def on_segment(p, q, r):
-    # Given three colinear points p, q, r, the function checks if
-    # point q lies on line segment 'pr'
-    if q.x <= max(p.x, r.x) and q.x >= min(p.x, r.x) and \
-       q.y <= max(p.y, r.y) and q.y >= min(p.y, r.y):
-        return True
-    return False
-
-
-def intersect(a, b, c, d):
-    # Return true if line segments ab and cd intersect
-    def ccw(a, b, c):
-        return (c.y - a.y) * (b.x - a.x) > (b.y - a.y) * (c.x - a.x)
-    return ccw(a, c, d) != ccw(b, c, d) and ccw(a, b, c) != ccw(a, b, d)
-
-
-def orientation(p, q, r):
-    # To find orientation of ordered triplet (p, q, r).
-    # The function returns following values
-    # 0 --> p, q and r are colinear
-    # 1 --> Clockwise
-    # 2 --> Counterclockwise
-
-    # See https://www.geeksforgeeks.org/orientation-3-ordered-points/
-    # for details of below formula.
-    val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y)
-
-    if val == 0:
-        return 0  # colinear
-
-    # clock or counterclock wise
-    if val > 0:
-        return 1
-    return 2
-
-
-def in_square(p1, x_min, x_max, y_min, y_max):
-    if p1.x < x_min or p1.x > x_max:
-        return False
-    if p1.y < y_min or p1.y > y_max:
-        return False
-    return True
-
-
-def in_circle(p1, p_center, radius):
-    return (p1.x - p_center.x)**2 + (p1.y - p_center.y)**2 < radius**2
-
-
-def do_intersect(p1, q1, p2, q2):
-    # The main function that returns true if line segment 'p1q1'
-    # and 'p2q2' intersect.
-    # Find the four orientations needed for general and
-    # special cases
-    o1 = orientation(p1, q1, p2)
-    o2 = orientation(p1, q1, q2)
-    o3 = orientation(p2, q2, p1)
-    o4 = orientation(p2, q2, q1)
-
-    # General case
-    if (o1 != o2 and o3 != o4):
-        return True
-
-    # Special Cases
-    # p1, q1 and p2 are colinear and p2 lies on segment p1q1
-    if (o1 == 0 and on_segment(p1, p2, q1)):
-        return True
-
-    # p1, q1 and q2 are colinear and q2 lies on segment p1q1
-    if (o2 == 0 and on_segment(p1, q2, q1)):
-        return True
-
-    # p2, q2 and p1 are colinear and p1 lies on segment p2q2
-    if (o3 == 0 and on_segment(p2, p1, q2)):
-        return True
-
-    # p2, q2 and q1 are colinear and q1 lies on segment p2q2
-    if (o4 == 0 and on_segment(p2, q1, q2)):
-        return True
-
-    return False

tsp/linkedlist.py

@@ -1,36 +0,0 @@
-
-
-class LinkedList(object):
-    def __init__(self):
-        self.first_node = None
-        self.last_node = None
-
-    def new_node(self, prev_node=None, data=None, next_node=None):
-        return {
-            "prev_node": prev_node,
-            "data": data,
-            "next_node": next_node}
-
-    def append(self, data):
-        new_node = self.new_node(self.last_node, data, None)
-        if not self.first_node:
-            self.first_node = new_node
-            self.last_node = new_node
-        else:
-            self.last_node["next_node"] = new_node
-            self.last_node = new_node
-        return new_node
-
-    def index(self, index):
-        for i, node in enumerate(self.iter()):
-            if i == index:
-                return node["data"]
-
-    def __iter__(self):
-        return self.iter()
-
-    def iter(self):
-        current_node = self.first_node
-        while current_node:
-            yield current_node['data']
-            current_node = current_node['next_node']

tsp/map.py

@@ -0,0 +1,143 @@
+import math
+
+
+class Map(object):
+    # Create Map. Cluster points into regions.  Calculate distances only to own
+    # and neighbor regions.  We can actually cluster in O(n) when we know how
+    # high and wide the clusters are. Once we have that working we go from
+    # there
+    CLUSTER_SIZE = 3  # How many points we want per cluster.
+
+    def __init__(self):
+        pass
+
+    def calc_corners(self, points):
+        x_min, x_max = float("inf"), float("-inf")
+        y_min, y_max = float("inf"), float("-inf")
+        for p in points:
+            if p.x < x_min:
+                x_min = p.x
+            if p.x > x_max:
+                x_max = p.x
+            if p.y < y_min:
+                y_min = p.y
+            if p.y > y_max:
+                y_max = p.y
+        self.x_min = x_min
+        self.x_max = x_max
+        self.y_min = y_min
+        self.y_max = y_max
+
+    def calc_cluster_dim(self, points):
+        clusters = len(points) // self.CLUSTER_SIZE
+        # Calculate number of clusters to have a square
+        self.clusters_x = math.ceil(math.sqrt(clusters))
+        self.clusters_y = self.clusters_x
+        self.clusters_total = self.clusters_x ** 2
+        self.cluster_x_dim = (self.x_max - self.x_min) / self.clusters_x
+        self.cluster_y_dim = (self.y_max - self.y_min) / self.clusters_y
+
+    def sort_points_into_clusters(self, points):
+        self.clusters = [[[]
+                          for x in range(self.clusters_y)]
+                          for y in range(self.clusters_y)]
+        for p in points:
+            cluster_x = int((p.x - self.x_min) // self.cluster_x_dim)
+            cluster_y = int((p.y - self.y_min) // self.cluster_y_dim)
+
+            # If the point is on the outer edge of the highest cluster
+            # the index will be outside the correct range. We put it
+            # into the closes cluster.
+            if cluster_x == self.clusters_x:
+                cluster_x -= 1
+            if cluster_y == self.clusters_y:
+                cluster_y -= 1
+
+            self.clusters[cluster_x][cluster_y].append(p)
+            p.cluster_x = cluster_x
+            p.cluster_y = cluster_y
+
+    def add_neighbors_to_points(self, points):
+        """ Add all points from the surrounding clusters to each point. """
+        for p in points:
+            clusters_x = [p.cluster_x]
+            clusters_y = [p.cluster_y]
+
+            if p.cluster_x - 1 >= 0:
+                clusters_x.append(p.cluster_x - 1)
+            if p.cluster_x + 1 < self.clusters_x:
+                clusters_x.append(p.cluster_x + 1)
+            if p.cluster_y - 1 >= 0:
+                clusters_y.append(p.cluster_y - 1)
+            if p.cluster_y + 1 < self.clusters_y:
+                clusters_y.append(p.cluster_y + 1)
+
+            clusters = [(x, y)
+                        for x in clusters_x
+                        for y in clusters_y]
+            neighbors = []
+            for x, y in clusters:
+                for p2 in self.clusters[x][y]:
+                    if p is not p2:
+                        neighbors.append(p2)
+            p.add_neighbors(neighbors)
+
+    def cluster(self, points):
+        """ Splits the map into clusters of a size so
+        that each cluster contains CLUSTER_SIZE points on
+        average. Adds all points from the current cluster
+        and the adjacent clusters to each point. """
+        self.calc_corners(points)
+        self.calc_cluster_dim(points)
+        self.sort_points_into_clusters(points)
+        self.add_neighbors_to_points(points)
+        return points
+
+    def plot(self, points):
+        try:
+            import matplotlib.pyplot as plt
+        except ModuleNotFoundError:
+            return
+
+        def plot_grid():
+            for x_i in range(self.clusters_x + 1):
+                x_1 = self.x_min + x_i * self.cluster_x_dim
+                x_2 = x_1
+                y_1 = self.y_min
+                y_2 = self.y_max
+                plt.plot([x_1, x_2], [y_1, y_2], 'b:')
+            for y_i in range(self.clusters_y + 1):
+                x_1 = self.x_min
+                x_2 = self.x_max
+                y_1 = self.y_min + y_i * self.cluster_y_dim
+                y_2 = y_1
+                plt.plot([x_1, x_2], [y_1, y_2], 'b:')
+
+        def plot_arrows():
+            for i in range(len_points):
+                p1 = points[i - 1]
+                p2 = points[i]
+                plot_arrow(p1, p2)
+
+        def plot_arrow(p1, p2):
+            x = p1.x
+            y = p1.y
+            dx = p2.x - x
+            dy = p2.y - y
+            opt = {'head_width': 0.4, 'head_length': 0.4, 'width': 0.05,
+                   'length_includes_head': True}
+            plt.arrow(x, y, dx, dy, **opt)
+
+        def plot_points():
+            for i, p in enumerate(points):
+                plt.plot(p.x, p.y, '')
+                plt.text(p.x, p.y, '  ' + str(p))
+                for nb, _ in p.neighbors:
+                    # plt.plot([p.x, nb.x], [p.y, nb.y], 'r--')
+                    pass
+
+        len_points = len(points)
+        plot_points()
+        plot_grid()
+        plot_arrows()
+        plt.show()

tsp/tsp.py

@@ -2,11 +2,44 @@ import math
 from functools import lru_cache
 from collections import namedtuple
 from random import shuffle
+from map import Map
 import time
 
+@lru_cache(maxsize=1000000)
+def distance(p1, p2):
+    """ Returns the distance between two points. """
+    return math.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)
+
+
+def float_is_equal(a, b):
+    if (a - b) < 0.001:
+        return True
+    return 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()))
+            for i in range(0, node_count)]
+
+
+def prepare_output_data(points):
+    # Basic plausibility checks
+    assert(len(set(points)) == len(points))
+    assert(len(points) > 4)
+    for i, p in enumerate(points):
+        assert(i == p.index)
+    r = Route(points)
+    obj = r.total_distance
+    output_data = '%.2f' % obj + ' ' + str(0) + '\n'
+    output_data += ' '.join(map(lambda p: str(p.id), points))
+    return output_data
+
 
 class Point(object):
     def __init__(self, index, x, y):
+        self.id = index
         self.index = index
         self.x = x
         self.y = y
@@ -23,48 +56,93 @@ class Point(object):
     def __str__(self):
         # m = "P_{}({}, {})".format(self.index, self.x, self.y)
         # m = "P_{}({}, {})".format(self.index, self.cluster_x, self.cluster_y)
-        m = "P({})".format(self.index)
+        id = self.id
+        index = self.index
+        m = f"P({id=}, {index=})"
         return m
 
     def __repr__(self):
         return self.__str__()
 
 
-def parse_input_data(input_data):
+class Route(object):
-    lines = input_data.split('\n')
+    def __init__(self, points):
-    node_count = int(lines[0])
+        self.points = points
-    return [Point(str(i), *map(float, lines[i + 1].split()))
+        self.len_points = len(points)
-            for i in range(0, node_count)]
+        self.total_distance = self.get_total_distance(self.points)
 
+    def get_total_distance(self, points):
+        """ Calculate the total distance of the point sequence. """
+        # Use negative indexing to get the distance from last to first point
+        return sum([distance(points[i - 1], points[i])
+                   for i in range(self.len_points)])
 
-def float_is_equal(a, b):
+    def swap(self, p1, p2):
-    if (a - b) < 0.001:
+        """
-        return True
+        Swaps two edges. p1 is the first point of the first
-    return False
+        edge and p2 is the first point of the second edge.
+        The first point of edge 1 (p1) points to the first point
+        of edge two (p2) after the swap, while the second point
+        of edge 1 (p12) points to the second point of edge two (p22).
+        This means we swap p12 and p2 and update their indices.
 
+            Before: p1 -> p12 and p2 -> p22
+            After: p1 -> p2 and p12 -> p22
 
-def prepare_output_data(points):
+        Afterwards we have to reverse the order of the points between
-    # Basic plausibility checks
+        p2 and p12 while those points themselves are no longer touched.
-    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: str(p.index), points))
-    return output_data
 
+        """
+        p12 = self.points[(p1.index + 1) % self.len_points]
+        p22 = self.points[(p2.index + 1) % self.len_points]
 
-@lru_cache(maxsize=1000000)
+        # Swap positions in route.
-def distance(point_1, point_2):
+        self.points[p12.index] = p2
-    """ Calculate the distance between two points. """
+        self.points[p2.index] = p12
-    p1, p2 = point_1, point_2
+        # Swap indices.
-    return math.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)
+        p2.index, p12.index = p12.index, p2.index
 
+        # TODO(felixm): Update self.total_distance.
+        # TODO(felixm): Reverse order between p2 and p12.
 
-def total_distance(points):
+        return points
-    """ Calculate the total distance of the point sequence. """
+
-    # Use negative indexing to get the distance from last to first point
+    def reorder_points_greedy(self):
-    return sum([distance(points[i - 1], points[i])
+        best_distance = float("inf")
-                for i in range(len(points))])
+        best_solution = None
+        points = self.points
+
+        for i in range(1000):
+            shuffle(points)
+            current_point, points = points[0], points[1:]
+            solution = [current_point]
+
+            while points:
+                next_point = None
+                # Select the closest point as the following one.
+                for neighbor, _ in current_point.neighbors:
+                    if neighbor in points:
+                        next_point = neighbor
+                        points.remove(next_point)
+                        break
+
+                # If none of the neighbors could be selected use any point.
+                if next_point is None:
+                    next_point = points.pop()
+
+                solution.append(next_point)
+                current_point = next_point
+
+            total_distance = self.get_total_distance(solution)
+            points = solution
+            if total_distance < best_distance:
+                best_distance = total_distance
+                best_solution = solution.copy()
+
+        self.points = best_solution
+        for i, p in enumerate(self.points):
+            p.index = i
+        return self.points
 
 
 def longest_distance(points, ignore_set):
@@ -117,40 +195,6 @@ def swap_edges(i, j, points, current_distance=0):
     return current_distance
 
 
-def reorder_points_greedy(points):
-    best_length = float("inf")
-    best_solution = None
-
-    for i in range(1000):
-        shuffle(points)
-        current_point, points = points[0], points[1:]
-        solution = [current_point]
-
-        while points:
-            next_point = None
-            # Select the closest point as the following one.
-            for neighbor, _ in current_point.neighbors:
-                if neighbor in points:
-                    next_point = neighbor
-                    points.remove(next_point)
-                    break
-
-            # If none of the neighbors could be selected use any point.
-            if next_point is None:
-                next_point = points.pop()
-
-            solution.append(next_point)
-            current_point = next_point
-
-        total_length = total_distance(solution)
-        points = solution
-        if total_length < best_length:
-            best_length = total_length
-            best_solution = solution.copy()
-
-    return best_solution
-
-
 def k_opt(p1_index, points, steps):
     ignore_set = set()
 
@@ -213,154 +257,17 @@ def local_search_k_opt(points):
     return points
 
 
-
-class Map(object):
-    # Create Map. Cluster points into regions.  Calculate distances only to own
-    # and neighbor regions.  We can actually cluster in O(n) when we know how
-    # high and wide the clusters are. Once we have that working we go from
-    # there
-    CLUSTER_SIZE = 3  # How many points we want per cluster.
-
-    def __init__(self, points):
-        self.points = points
-        self.num_points = len(points)
-
-        self.calc_corners()
-        self.calc_cluster_dim()
-        self.sort_points_into_clusters()
-        self.add_neighbors_to_points()
-
-    def calc_cluster_dim(self):
-        clusters = self.num_points // self.CLUSTER_SIZE
-        # Calculate number of clusters to have a square
-        self.clusters_x = math.ceil(math.sqrt(clusters))
-        self.clusters_y = self.clusters_x
-        self.clusters_total = self.clusters_x ** 2
-        self.cluster_x_dim = (self.x_max - self.x_min) / self.clusters_x
-        self.cluster_y_dim = (self.y_max - self.y_min) / self.clusters_y
-
-    def add_neighbors_to_points(self):
-        """ Add all points from the surrounding clusters to each point. """
-        for p in self.points:
-            clusters_x = [p.cluster_x]
-            clusters_y = [p.cluster_y]
-
-            if p.cluster_x - 1 >= 0:
-                clusters_x.append(p.cluster_x - 1)
-            if p.cluster_x + 1 < self.clusters_x:
-                clusters_x.append(p.cluster_x + 1)
-            if p.cluster_y - 1 >= 0:
-                clusters_y.append(p.cluster_y - 1)
-            if p.cluster_y + 1 < self.clusters_y:
-                clusters_y.append(p.cluster_y + 1)
-
-            clusters = [(x, y)
-                        for x in clusters_x
-                        for y in clusters_y]
-            neighbors = []
-            for x, y in clusters:
-                for p2 in self.clusters[x][y]:
-                    if p is not p2:
-                        neighbors.append(p2)
-            p.add_neighbors(neighbors)
-
-    def sort_points_into_clusters(self):
-        self.clusters = [[[]
-                          for x in range(self.clusters_y)]
-                          for y in range(self.clusters_y)]
-        for p in self.points:
-            cluster_x = int((p.x - self.x_min) // self.cluster_x_dim)
-            cluster_y = int((p.y - self.y_min) // self.cluster_y_dim)
-
-            # If the point is on the outer edge of the highest cluster
-            # the index will be outside the correct range. We put it
-            # into the closes cluster.
-            if cluster_x == self.clusters_x:
-                cluster_x -= 1
-            if cluster_y == self.clusters_y:
-                cluster_y -= 1
-
-            self.clusters[cluster_x][cluster_y].append(p)
-            p.cluster_x = cluster_x
-            p.cluster_y = cluster_y
-
-    def calc_corners(self):
-        x_min, x_max = float("inf"), float("-inf")
-        y_min, y_max = float("inf"), float("-inf")
-        for p in self.points:
-            if p.x < x_min:
-                x_min = p.x
-            if p.x > x_max:
-                x_max = p.x
-            if p.y < y_min:
-                y_min = p.y
-            if p.y > y_max:
-                y_max = p.y
-        self.x_min = x_min
-        self.x_max = x_max
-        self.y_min = y_min
-        self.y_max = y_max
-
-    def plot(self):
-        try:
-            import matplotlib.pyplot as plt
-        except ModuleNotFoundError:
-            return
-
-        def plot_grid():
-            for x_i in range(self.clusters_x + 1):
-                x_1 = self.x_min + x_i * self.cluster_x_dim
-                x_2 = x_1
-                y_1 = self.y_min
-                y_2 = self.y_max
-                plt.plot([x_1, x_2], [y_1, y_2], 'b:')
-            for y_i in range(self.clusters_y + 1):
-                x_1 = self.x_min
-                x_2 = self.x_max
-                y_1 = self.y_min + y_i * self.cluster_y_dim
-                y_2 = y_1
-                plt.plot([x_1, x_2], [y_1, y_2], 'b:')
-
-        def plot_arrows():
-            for i in range(self.num_points):
-                p1 = self.points[i - 1]
-                p2 = self.points[i]
-                plot_arrow(p1, p2)
-
-        def plot_arrow(p1, p2):
-            x = p1.x
-            y = p1.y
-            dx = p2.x - x
-            dy = p2.y - y
-            opt = {'head_width': 0.4, 'head_length': 0.4, 'width': 0.05,
-                   'length_includes_head': True}
-            plt.arrow(x, y, dx, dy, **opt)
-
-        def plot_points():
-            for i, p in enumerate(self.points):
-                plt.plot(p.x, p.y, '')
-                plt.text(p.x, p.y, '  ' + str(p))
-                for nb, _ in p.neighbors:
-                    # plt.plot([p.x, nb.x], [p.y, nb.y], 'r--')
-                    pass
-
-        plot_points()
-        plot_grid()
-        plot_arrows()
-        plt.show()
-
-
 def solve_it(input_data):
-    points = parse_input_data(input_data)
+    r = Route(parse_input_data(input_data))
-    # Initialiaze map before algorithm because it clusters the points
+    m = Map()
-    # and adds the neighbors to each point.
+
-    m = Map(points)
+    m.cluster(r.points)
-    m.points = reorder_points_greedy(points)
+    r.reorder_points_greedy()
-    # FIXME(felixm): Don't do this here.
+    # m.plot(r.points)
-    m.points = local_search_k_opt(m.points)
+
-    m.plot()
+    return prepare_output_data(r.points)
+
 
-    return prepare_output_data(m.points)
 
 
 if __name__ == "__main__":
@@ -369,32 +276,3 @@ if __name__ == "__main__":
         input_data = input_data_file.read()
         print(solve_it(input_data))
 
-
-def local_search_2_opt(points):
-    current_total = total_distance(points)
-    ignore_set = set()
-    while True:
-        pi, i = longest_distance(points, ignore_set)
-        ignore_set.add(pi)
-        if not pi:
-            break
-
-        best_new_total = current_total
-        best_points = None
-        swap = None
-        for j in range(len(points)):
-            if j in [i, i + 1, i + 2]:
-                continue
-            new_points = list(points)
-            swap_edges(i, j - 1, new_points)
-            new_total = total_distance(new_points)
-            if new_total < best_new_total:
-                swap = (points[i], points[j - 1])
-                best_new_total = new_total
-                best_points = new_points
-
-        if best_new_total < current_total:
-            current_total = best_new_total
-            points = best_points
-            ignore_set = set()
-    return points