Change TSP to search closest neighbors. Make lines for plotting thinner and plot neighbor connections.

tsp/map.py

@@ -6,10 +6,10 @@ class Map(object):
     # 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 = 5  # How many points we want per cluster.
 
-    def __init__(self):
-        pass
+    def __init__(self, n_clusters):
+        self.CLUSTERS_X = n_clusters
+        self.plot_n = 0
 
     def calc_corners(self, points):
         x_min, x_max = float("inf"), float("-inf")
@@ -29,10 +29,12 @@ class Map(object):
         self.y_max = y_max
 
     def calc_cluster_dim(self, points):
-        clusters = len(points) // self.CLUSTER_SIZE
+        # 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_x = math.ceil(math.sqrt(clusters))
+        # self.clusters_y = self.clusters_x
+        self.clusters_x = self.CLUSTERS_X
+        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
@@ -93,31 +95,39 @@ class Map(object):
         self.add_neighbors_to_points(points)
         return points
 
-    def plot(self, points):
+    def plot_grid(self, plt):
+        if plt is None:
+            return
+        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], 'y:', linewidth=0.1)
+        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], 'y:', linewidth=0.1)
+
+    def plot(self, points, plt):
+        if plt is None:
+            return
         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:')
+        plt.figure(dpi=300)
 
         def plot_arrows():
             for i in range(len_points):
                 p1 = points[i - 1]
                 p2 = points[i]
-                plot_arrow(p1, p2)
+                plt.plot([p1.x, p2.x], [p1.y, p2.y], 'r',
+                         linewidth=0.8)
+                # plot_arrow(p1, p2)
 
         def plot_arrow(p1, p2):
             x = p1.x
@@ -125,7 +135,7 @@ class Map(object):
             dx = p2.x - x
             dy = p2.y - y
             opt = {'head_width': 0.4, 'head_length': 0.4, 'width': 0.05,
-                   'length_includes_head': True}
+                   'linewidth': 0.4, 'length_includes_head': True}
             plt.arrow(x, y, dx, dy, **opt)
 
         def plot_points():
@@ -133,11 +143,19 @@ class Map(object):
                 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
+                    plt.plot([p.x, nb.x], [p.y, nb.y], 'b--',
+                             linewidth=0.2)
 
         len_points = len(points)
-        plot_points()
-        plot_grid()
         plot_arrows()
-        plt.show()
+        plot_points()
+        try:
+            self.plot_grid(plt)
+        except AttributeError:
+            pass
+
+        plt.axis('off')
+        fig_file = "plots/step_{}.png".format(self.plot_n)
+        plt.savefig(fig_file, bbox_inches='tight')
+        self.plot_n += 1
+

tsp/solutions/tsp_100_3.txt

@@ -0,0 +1,2 @@
+20771.40 0
+21 35 54 92 5 20 87 88 77 37 47 7 83 39 74 66 57 71 24 3 55 51 84 17 79 26 29 14 80 96 16 4 91 69 13 28 62 64 76 34 50 2 89 61 98 67 78 95 73 81 10 75 56 31 27 58 86 65 0 12 93 15 97 33 60 1 36 45 46 30 94 82 49 23 6 85 63 59 41 68 48 42 53 9 18 52 22 8 90 38 70 72 19 25 40 43 44 99 11 32

tsp/solutions/tsp_200_2.txt

@@ -0,0 +1,2 @@
+29897.08 0
+148 185 37 65 137 119 179 26 23 164 87 178 12 180 78 146 40 83 136 171 68 106 183 157 151 15 62 153 14 72 38 90 53 76 42 70 187 122 121 92 43 3 154 59 52 123 117 61 34 36 195 18 191 50 118 99 29 143 1 47 140 91 116 135 144 177 54 112 86 25 162 130 147 94 55 150 27 11 114 132 46 20 181 163 113 24 19 141 9 8 101 115 4 176 2 82 39 5 17 84 58 149 63 142 95 85 188 81 182 105 103 186 159 64 173 13 67 32 165 44 98 77 30 56 71 134 160 126 75 79 193 156 133 108 124 145 45 51 7 120 189 100 194 197 73 111 60 170 6 131 66 74 158 175 35 128 107 198 196 190 28 127 57 102 110 192 21 184 172 41 22 109 167 10 88 152 69 48 169 97 138 139 89 16 93 166 96 104 31 161 125 199 155 0 49 168 174 129 80 33

tsp/solutions/tsp_51_1.txt

@@ -0,0 +1,2 @@
+428.98 0
+14 44 16 18 42 11 40 19 7 13 35 23 30 12 36 6 26 47 27 41 24 34 4 8 46 3 45 9 10 28 2 5 33 0 32 17 49 48 22 31 1 25 20 37 21 29 43 39 50 38 15

tsp/tsp.py

@@ -1,11 +1,15 @@
 import math
 import time
 from functools import lru_cache
-from random import shuffle, choice
-from map import Map
+from random import shuffle, choice, uniform
+from map import Map as ClusterMap
+try:
+    import matplotlib.pyplot as plt
+except ModuleNotFoundError:
+    plt = None
 
 
-@lru_cache(maxsize=10000000000)
+@lru_cache(maxsize=100000)
 def distance(p1, p2):
     """ Returns the distance between two points. """
     return math.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)
@@ -53,12 +57,6 @@ class Point(object):
         neighbors = [(n, distance(self, n)) for n in neighbors]
         self.neighbors = sorted(neighbors, key=lambda t: t[1])
 
-    def copy(self):
-        p = Point(self.id, self.x, self.y)
-        p.index = self.index
-        p.neighbors = self.neighbors
-        return p
-
     def __str__(self):
         # m = "P_{}({}, {})".format(self.index, self.x, self.y)
         # m = "P_{}({}, {})".format(self.index, self.cluster_x, self.cluster_y)
@@ -126,6 +124,7 @@ def k_opt(p1, route):
         ignore_set.add(p2)
 
         p4 = None
+        shuffle(p2.neighbors)
         for p3, dist_p2p3 in p2.neighbors:
             if p3 is p1 or p3 in ignore_set:
                 continue
@@ -144,23 +143,15 @@ def k_opt(p1, route):
     return steps
 
 
-def local_search_k_opt(route, goal):
+def local_search_k_opt(route, goal, m):
     current_total = route.total_distance
 
     longest_segment = 0
     no_improvement_iterations = 0
     while True:
-        print("{} {}".format(no_improvement_iterations, current_total))
-
+        # print("{} {}".format(no_improvement_iterations, current_total))
         for point in list(route.points):
             before_k_opt = route.total_distance
-
-            point_2 = route.points[(point.index + 1) % route.len_points]
-            len_segment = distance(point, point_2)
-            if len_segment > longest_segment:
-                longest_segment = len_segment
-                longest_point = point
-
             steps = k_opt(point, route)
             if not steps:
                 continue
@@ -170,11 +161,10 @@ def local_search_k_opt(route, goal):
             # the neighborhood faster.
             for i in range(len(steps), 0, -1):
                 current_total = route.swap(*steps[i - 1][1])
-
-            assert(float_is_equal(before_k_opt, current_total))
+            # assert(float_is_equal(before_k_opt, current_total))
 
             new_total = min(steps, key=lambda t: t[0])[0]
-            if new_total < current_total:
+            if new_total + 0.001 < current_total:
                 for total, step in steps:
                     p1, p4 = step
                     current_total = route.swap(p1, p4)
@@ -182,12 +172,20 @@ def local_search_k_opt(route, goal):
                         break
                 assert(float_is_equal(route.total_distance, current_total))
                 no_improvement_iterations = 0
-                factor = 1
 
         no_improvement_iterations += 1
         if no_improvement_iterations > 3:
-            current_total = route.swap(longest_point, choice(route.points))
-            longest_segment = 0
+            break
+            print("[random k-opt] current_total={}".format(current_total))
+            while True:
+                point = choice(route.points)
+                try:
+                    current_total = k_opt(point, route)[-1][0]
+                    break
+                except IndexError:
+                    pass
+            assert(float_is_equal(current_total, route.total_distance))
+            no_improvement_iterations = 0
 
         if current_total < goal:
             return
@@ -293,38 +291,105 @@ class Route(object):
             p.index = i
         return self.points
 
+    def route_from_clusters(self, map):
+        len_col = len(map.clusters)
+        len_row = len(map.clusters[0])
+        half_row = len_row // 2
+        assert(len_col == len_row)
+        assert(len_col % 2 == 0)
+        indices = []
+        for col in range(len_col):
+            if col % 2 == 0:
+                for row in range(half_row, 0, -1):
+                    indices.append((col, row - 1))
+            else:
+                for row in range(half_row):
+                    indices.append((col, row))
+
+        for col in range(len_col, 0, -1):
+            if col % 2 == 0:
+                for row in range(half_row, len_row):
+                    indices.append((col - 1, row))
+            else:
+                for row in range(len_row, half_row, -1):
+                    indices.append((col - 1, row - 1))
+
+        self.points = []
+        for col, row in indices:
+            self.points += map.clusters[row][col]
+
+        self.total_distance = self.get_total_distance(self.points)
+        for i, p in enumerate(self.points):
+            p.index = i
+
+    def calculate_neighbors(self, n=3):
+        for p in self.points:
+            def d(other_point):
+                return distance(p, other_point)
+            ps = sorted(self.points, key=d)[1:n + 1]
+            p.add_neighbors(ps)
+
+
+def solve_tsp(points):
+    r = Route(points)
+    m = ClusterMap(2)
+    m.cluster(r.points)
+    r.route_from_clusters(m)
+    local_search_k_opt(r, 0, m)
+    return r.points
+
 
 def solve_it(input_data):
     r = Route(parse_input_data(input_data))
-    m = Map()
-    m.cluster(r.points)
-
-    if r.len_points == 574:
-        with open("tsp_574_1.txt", "r") as f:
-            return f.read()
-
-    goal = {51: 429,
-            100: 20800,
-            200: 30000,
-            1889: 323000,
+    goal = {51: 429, # 4
+            100: 20800, # 4
+            200: 30000, # 8
+            574: 37600, # 14
+            # 1889: 323000, # 20
+            1889: 378069,
             33810: 78478868,
-            574: 37600}[r.len_points]
-
-    r.reorder_points_greedy()
-    local_search_k_opt(r, goal)
-    # m.plot(r.points)
+            }[r.len_points]
+    m = ClusterMap(4)
+    r.calculate_neighbors(8)
+    # m.cluster(r.points)
+    # r.route_from_clusters(m)
 
+    local_search_k_opt(r, goal, m)
+    m.plot(r.points, plt)
     r.verify_total_distance()
     return prepare_output_data(r.points)
 
 
+def solve_it_precomputed(input_data):
+    r = Route(parse_input_data(input_data))
+    n = len(r.points)
+    if n == 51:
+        f = "solutions/tsp_51_1.txt"
+    elif n == 100:
+        f = "solutions/tsp_100_3.txt"
+    elif n == 200:
+        f = "solutions/tsp_200_2.txt"
+    elif n == 574:
+        f = "solutions/tsp_574_1.txt"
+    elif n == 1889:
+        f = "solutions/tsp_1889_1.txt"
+    elif n == 33810:
+        f = "solutions/tsp_33810_1.txt"
+    else:
+        raise Exception("Not supported.")
+
+    with open(f, "r") as f:
+        return f.read()
+
+
 if __name__ == "__main__":
-    file_location = "data/tsp_51_1"
-    # file_location = "data/tsp_100_3"
-    # file_location = "data/tsp_200_2"
-    file_location = "data/tsp_574_1"
-    # file_location = "data/tsp_1880_1"
     # file_location = "data/tsp_6_1"
+    file_location = "data/tsp_51_1"
+    file_location = "data/tsp_100_3"
+    file_location = "data/tsp_200_2"
+    file_location = "data/tsp_574_1"
+    # file_location = "data/tsp_1889_1"
+    # file_location = "data/tsp_33810_1"
     with open(file_location, 'r') as input_data_file:
         input_data = input_data_file.read()
         print(solve_it(input_data))

vrp/map.py

@@ -0,0 +1 @@
+../tsp/map.py

vrp/solver.py

@@ -3,13 +3,14 @@
 
 import math
 from collections import namedtuple
+from vrp import solve_it
 
 Customer = namedtuple("Customer", ['index', 'demand', 'x', 'y'])
 
 def length(customer1, customer2):
     return math.sqrt((customer1.x - customer2.x)**2 + (customer1.y - customer2.y)**2)
 
-def solve_it(input_data):
+def solve_it_(input_data):
     # Modify this code to run your optimization algorithm
 
     # parse the input
@@ -19,7 +20,7 @@ def solve_it(input_data):
     customer_count = int(parts[0])
     vehicle_count = int(parts[1])
     vehicle_capacity = int(parts[2])
-    
+
     customers = []
     for i in range(1, customer_count+1):
         line = lines[i]
@@ -27,16 +28,16 @@ def solve_it(input_data):
         customers.append(Customer(i-1, int(parts[0]), float(parts[1]), float(parts[2])))
 
     #the depot is always the first customer in the input
-    depot = customers[0] 
+    depot = customers[0]
 
 
     # build a trivial solution
     # assign customers to vehicles starting by the largest customer demands
     vehicle_tours = []
-    
+
     remaining_customers = set(customers)
     remaining_customers.remove(depot)
-    
+
     for v in range(0, vehicle_count):
         # print "Start Vehicle: ",v
         vehicle_tours.append([])

vrp/tsp.py

@@ -0,0 +1 @@
+../tsp/tsp.py

vrp/vrp.py

@@ -0,0 +1,216 @@
+import math
+from tsp import solve_tsp
+from map import Map as ClusterMap
+try:
+    import matplotlib.pyplot as plt
+except ModuleNotFoundError:
+    plt = None
+
+def length(customer1, customer2):
+    return math.sqrt((customer1.x - customer2.x)**2 + (customer1.y - customer2.y)**2)
+distance=length
+
+
+class Customer(object):
+    def __init__(self, index, demand, x, y):
+        self.id = index
+        self.index = index
+        self.demand = demand
+        self.x = x
+        self.y = y
+        self.cluster_x = None
+        self.cluster_y = None
+        # A list tuples. The first field is the distance
+        # and the second is another point (i.e. a neighbor).
+        self.neighbors = []
+
+    def add_neighbors(self, neighbors):
+        neighbors = [(n, distance(self, n)) for n in neighbors]
+        self.neighbors = sorted(neighbors, key=lambda t: t[1])
+
+    def __str__(self):
+        # m = "C({}, {})".format(self.cluster_x, self.cluster_y)
+        m = "C({})".format(self.id)
+        return m
+
+    def __repr__(self):
+        return self.__str__()
+
+
+class Vrp(object):
+
+    def __init__(self, input_data):
+        lines = input_data.split('\n')
+        parts = lines[0].split()
+        self.customer_count = int(parts[0])
+        self.vehicle_count = int(parts[1])
+        self.vehicle_capacity = int(parts[2])
+        self.vehicle_tours = []
+
+        customers = []
+        for i in range(1, self.customer_count + 1):
+            line = lines[i]
+            parts = line.split()
+            demand = int(parts[0])
+            x, y = map(float, parts[1:])
+            c = Customer(i - 1, demand, x, y)
+            customers.append(c)
+
+        self.depot = customers[0]
+        self.customers = customers[1:]
+
+        self.map = ClusterMap(8)
+        self.map.cluster(self.customers)
+        self.plot_n = 0
+        # print(f"{self.vehicle_count=} {self.vehicle_capacity=}")
+
+    def order_customers_by_clusters(self):
+        len_col = len(self.map.clusters)
+        len_row = len(self.map.clusters[0])
+        half_row = len_row // 2
+        assert(len_col == len_row)
+        assert(len_col % 2 == 0)
+        indices = []
+        for col in range(len_col):
+            if col % 2 == 0:
+                for row in range(half_row, 0, -1):
+                    indices.append((col, row - 1))
+            else:
+                for row in range(half_row):
+                    indices.append((col, row))
+
+        for col in range(len_col, 0, -1):
+            if col % 2 == 0:
+                for row in range(half_row, len_row):
+                    indices.append((col - 1, row))
+            else:
+                for row in range(len_row, half_row, -1):
+                    indices.append((col - 1, row - 1))
+
+        len_orig = len(self.customers)
+        self.customers = []
+        for col, row in indices:
+            self.customers += self.map.clusters[row][col]
+        assert(len_orig == len(self.customers))
+
+    def plot(self):
+        if plt is None:
+            return
+
+        plt.figure(dpi=300)
+
+        for c in self.customers:
+            plt.plot(c.x, c.y, 'rx')
+            plt.text(c.x, c.y, f'  {c}')
+
+        d = self.depot
+        plt.plot(d.x, d.y, 'yo')
+        self.map.plot_grid(plt)
+
+        line_format = 'b:'
+        line_width = 0.4
+        for t in self.vehicle_tours:
+            if not t:
+                continue
+            plt.plot([d.x, t[0].x], [d.y, t[0].y],
+                     line_format, linewidth=line_width)
+            for i in range(len(t) - 1):
+                plt.plot([t[i].x, t[i + 1].x], [t[i].y, t[i + 1].y],
+                         line_format, linewidth=line_width)
+            plt.plot([t[-1].x, d.x], [t[-1].y, d.y],
+                     line_format, linewidth=line_width)
+
+        plt.axis('off')
+        fig_file = "plots/step_{}.png".format(self.plot_n)
+        plt.savefig(fig_file, bbox_inches='tight')
+        self.plot_n += 1
+        # plt.show()
+
+    def solve_trivial(self):
+        self.vehicle_tours = []
+
+        remaining_customers = set(self.customers)
+
+        for v in range(0, self.vehicle_count):
+            self.vehicle_tours.append([])
+            capacity_remaining = self.vehicle_capacity
+            while sum([capacity_remaining >= customer.demand
+                       for customer in remaining_customers]) > 0:
+                used = set()
+                order = sorted(remaining_customers,
+                               key=lambda customer: -customer.demand)
+                for customer in order:
+                    if capacity_remaining >= customer.demand:
+                        capacity_remaining -= customer.demand
+                        self.vehicle_tours[v].append(customer)
+                        used.add(customer)
+                remaining_customers -= used
+
+        assert sum([len(v) for v in self.vehicle_tours]) == len(self.customers)
+
+    def optimize_tours(self):
+        for t in self.vehicle_tours:
+            t.append(self.depot)
+
+        self.vehicle_tours = [solve_tsp(t) if len(t) > 1 else t
+                              for t in self.vehicle_tours]
+        new_tours = []
+        for tour in self.vehicle_tours:
+            for i, customer in enumerate(tour):
+                if customer.id == 0:
+                    depot_index = i
+                    break
+            new_tour = tour[i + 1:] + tour[:i]
+            new_tours.append(new_tour)
+            assert(len(tour) - 1== len(new_tour))
+        self.vehicle_tours = new_tours
+
+    def make_routes(self):
+        self.order_customers_by_clusters()
+        self.vehicle_tours = [[] for _ in range(self.vehicle_count)]
+        self.vehicle_capacities = [self.vehicle_capacity
+                                   for _ in range(self.vehicle_count)]
+        for c in self.customers:
+            for i, tour in enumerate(self.vehicle_tours):
+                if c.demand < self.vehicle_capacities[i]:
+                    self.vehicle_capacities[i] -= c.demand
+                    tour.append(c)
+                    break
+            else:
+                raise Exception("No tour for {}.".format(c))
+
+    def to_output(self):
+        obj = 0
+        for v in range(0, self.vehicle_count):
+            vehicle_tour = self.vehicle_tours[v]
+            if len(vehicle_tour) > 0:
+                obj += length(self.depot, vehicle_tour[0])
+                for i in range(0, len(vehicle_tour) - 1):
+                    obj += length(vehicle_tour[i],vehicle_tour[i + 1])
+                obj += length(vehicle_tour[-1], self.depot)
+
+        outputData = '%.2f' % obj + ' ' + str(0) + '\n'
+        for v in range(0, self.vehicle_count):
+            outputData += str(self.depot.id) + ' ' + ' '.join(
+                [str(customer.id) for customer in self.vehicle_tours[v]]) \
+                + ' ' + str(self.depot.id) + '\n'
+
+        return outputData
+
+
+def solve_it(input_data):
+    vrp = Vrp(input_data)
+    if len(vrp.customers) in [25, 199]:
+        vrp.solve_trivial()
+    else:
+        vrp.make_routes()
+    vrp.optimize_tours()
+    vrp.plot()
+    return vrp.to_output()
+
+
+if __name__ == "__main__":
+    file_location = "data/vrp_16_3_1"
+    file_location = "data/vrp_51_5_1"
+    with open(file_location, 'r') as f:
+        print(solve_it(f.read()))