import math import time from functools import lru_cache 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=100000) 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 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 = "P_{}({}, {})".format(self.index, self.x, self.y) # m = "P_{}({}, {})".format(self.index, self.cluster_x, self.cluster_y) m = "{}".format(self.index) return m def __repr__(self): return self.__str__() def longest_distance(points, ignore_set): """ Returns the point and index of the point with the longest distance to the next point. """ longest_distance = 0 longest_dist_point = None for i in range(len(points)): p1, p2 = points[i - 1], points[i] if p1 in ignore_set: continue current_distance = distance(p1, p2) if current_distance > longest_distance: longest_distance = current_distance longest_dist_point = p1 return longest_dist_point def swap_edges(i, j, points, current_distance=0): """ 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. """ current_distance = total_distance(points) p11, p12 = points[i], points[i + 1] p21, p22 = points[j], points[j + 1] points[i + 1] = p21 points[j] = p12 current_distance -= (distance(p11, p12) + distance(p21, p22)) current_distance += (distance(p11, p21) + distance(p12, p22)) # If we do not correct j = -1 the reverse logic breaks for that case. if j == -1: j = len(points) - 1 # 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 current_distance def k_opt(p1, route): steps = [] ignore_set = set() for _ in range(10): p2 = route.points[(p1.index + 1) % route.len_points] dist_p1p2 = distance(p1, p2) 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 p4 = route.points[(p3.index - 1) % route.len_points] if p4 in ignore_set or p4 is p1: continue if dist_p2p3 < dist_p1p2: break p4 = None if p4 is None: break step = (p1.index, p4.index) new_total = route.swap(p1, p4) steps.append((new_total, step)) return steps 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)) for point in list(route.points): before_k_opt = route.total_distance steps = k_opt(point, route) if not steps: continue # Reverse route to status before k_opt. This is cheaper # than copying the whole route and allows us to search # 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)) new_total = min(steps, key=lambda t: t[0])[0] if new_total + 0.001 < current_total: for total, step in steps: p1, p4 = step current_total = route.swap(p1, p4) if total == new_total: break assert(float_is_equal(route.total_distance, current_total)) no_improvement_iterations = 0 no_improvement_iterations += 1 if no_improvement_iterations > 3: 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 class Route(object): def __init__(self, points): self.points = points self.len_points = len(points) self.total_distance = self.get_total_distance(points) self.point_id_to_point = {p.id: p for p in self.points} def verify_total_distance(self): a = self.total_distance b = self.get_total_distance(self.points) assert(float_is_equal(a, b)) 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 swap(self, p1, p2): """ Swaps two edges. p1 is the first point of the first 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 Afterwards we have to reverse the order of the points [p', p'', p'''] between p2 and p12 while those points themselves are no longer touched. Before swap: [p1, p12, p', p'', p''', p2, p21] After swap: [p1, p2, p', p'', p''', p12, p21] """ if type(p1) is int: p1 = self.points[p1] if type(p2) is int: p2 = self.points[p2] if p1 is p2: return self.total_distance # Handle case when edge goes over the end of the list. p12 = self.points[(p1.index + 1) % self.len_points] p21 = self.points[(p2.index + 1) % self.len_points] # Update self.total_distance. self.total_distance -= (distance(p1, p12) + distance(p2, p21)) self.total_distance += (distance(p1, p2) + distance(p12, p21)) # Swap positions and indices. self.points[p12.index] = p2 self.points[p2.index] = p12 p2.index, p12.index = p12.index, p2.index # Handle case when p2 was before p1 initially. if p12.index > p2.index: len_revers = p12.index - p2.index else: len_revers = (p12.index + self.len_points) - p2.index # Reverse order between p2 and p12. for i in range(1, len_revers // 2 + 1): pa = self.points[(p2.index + i) % self.len_points] pb = self.points[(p12.index - i) % self.len_points] self.points[pa.index] = pb self.points[pb.index] = pa pa.index, pb.index = pb.index, pa.index return self.total_distance def reorder_points_greedy(self): points = list(self.points) current_point, points = points[0], points[1:] self.points = [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() self.points.append(next_point) current_point = next_point self.total_distance = self.get_total_distance(self.points) for i, p in enumerate(self.points): 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)) goal = {51: 429, # 4 100: 20800, # 4 200: 30000, # 8 574: 37600, # 14 # 1889: 323000, # 20 1889: 378069, 33810: 78478868, }[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_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))