discrete_optimization/tsp/tsp.py

391 lines
12 KiB
Python

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 > 10:
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 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 = ClusterMap(4)
m.cluster(r.points)
goal = {51: 429, # 4
100: 20800, # 4
200: 30000, # 8
574: 37600, # 14
# 1889: 323000, # 20
1889: 378069,
33810: 78478868,
}[r.len_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_(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))