felixm/discrete_optimization
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Eod commit.

Browse Source
This commit is contained in:
Felix Martin 2019-12-23 17:29:44 -05:00
parent 7481ee3f0a
commit 19c1515c0e
3 changed files with 169 additions and 123 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
tsp/geometry.py
View File

@ -8,6 +8,13 @@ def on_segment(p, q, r):
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

36
tsp/linkedlist.py Normal file
View File

@ -0,0 +1,36 @@
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']

249
tsp/tsp.py
View File

@ -1,7 +1,7 @@
import math
from functools import lru_cache
from collections import namedtuple
from geometry import do_intersect
from geometry import intersect
Point = namedtuple("Point", ['index', 'x', 'y'])
@ -9,33 +9,32 @@ Point = namedtuple("Point", ['index', 'x', 'y'])
def parse_input_data(input_data):
lines = input_data.split('\n')
node_count = int(lines[0])
return [Point(i, *map(float, lines[i].split()))
for i in range(1, node_count + 1)]
return [Point(i, *map(float, lines[i + 1].split()))
for i in range(0, node_count)]
def plot_solution(solution, points):
def plot_graph(points):
import matplotlib.pyplot as plt
def plot_arrows():
for i in range(len(solution)):
p_1_index = solution[i - 1]
p_2_index = solution[i]
plot_arrow(p_1_index, p_2_index)
for i in range(len(points)):
p1 = points[i - 1]
p2 = points[i]
plot_arrow(p1, p2)
def plot_arrow(p_1_index, p_2_index):
p_1 = points[p_1_index]
p_2 = points[p_2_index]
x = p_1.x
y = p_1.y
dx = p_2.x - x
dy = p_2.y - y
plt.arrow(x, y, dx, dy,
head_width=0.5, head_length=0.5)
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):
for p in points:
plt.plot(p.x, p.y, '')
plt.text(p.x, p.y, ' ' + str(i))
plt.text(p.x, p.y, ' ' + str(p.index))
plot_points()
plot_arrows()
@ -45,78 +44,31 @@ def plot_solution(solution, points):
def solve_it(input_data):
@lru_cache(maxsize=1000000)
def length(point_1_index, point_2_index):
p_1 = points[point_1_index]
p_2 = points[point_2_index]
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(solution):
obj = calculate_length(solution)
def prepare_output_data(points):
obj = total_length(points)
output_data = '%.2f' % obj + ' ' + str(0) + '\n'
output_data += ' '.join(map(str, solution))
output_data += ' '.join(map(lambda p: str(p.index), points))
return output_data
def is_valid(solution):
points = set(range(len(solution)))
assert(set(solution) == points)
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 calculate_length(solution):
obj = 0
for i in range(0, len(solution)):
point_1_index = solution[i - 1]
point_2_index = solution[i]
obj += length(point_1_index, point_2_index)
return obj
def total_length(points):
return sum([length(points[i - 1], points[i])
for i in range(len(points))])
def initial_solution_naiv():
return list(range(0, node_count))
return []
def does_edge_cause_intersection(edge, edges):
p1 = points[edge[0]]
p2 = points[edge[1]]
for existing_edge in edges:
q1 = points[existing_edge[0]]
q2 = points[existing_edge[1]]
print(p1, p2, q1, q2)
if do_intersect(p1, q1, p2, q2):
return True
return False
def get_dimensions(point_indices):
p = points[0]
x_min_p = p
x_max_p = p
y_min_p = p
y_max_p = p
x_min = p.x
x_max = p.x
y_min = p.y
y_max = p.y
for p in points:
if p.x < x_min:
x_min = p.x
x_min_p = p
if p.y < y_min:
y_min = p.y
y_min_p = p
if p.x > x_max:
x_max = p.x
x_max_p = p
if p.y > y_max:
y_max = p.y
y_max_p = p
return (x_min_p, x_max_p, y_min_p, y_max_p)
def initial_solution_greedy(point_indices):
current_point = get_dimensions(point_indices)[0].index
xs = [current_point]
points = set(point_indices)
points.remove(current_point)
def reorder_points_greedy(points):
current_point = points[0]
solution = [current_point]
points = set(points[1:])
while points:
min_length = 999999
@ -126,50 +78,101 @@ def solve_it(input_data):
if new_length < min_length:
min_length = new_length
min_point = next_point
xs.append(min_point)
points.remove(min_point)
current_point = min_point
return xs
def local_search(solution):
# Find longest edges to swap
max_len_1 = 0
max_len_2 = 0
edge_1 = None
edge_2 = None
for i in range(node_count):
new_len = length(solution[i - 1], solution[i])
if new_len > max_len_1:
max_len_1 = new_len
edge_1 = (i - 1, i)
for i in range(node_count):
new_len = length(solution[i - 1], solution[i])
if new_len > max_len_2 and new_len != max_len_1:
max_len_2 = new_len
edge_2 = (i - 1, i)
a_1, a_2 = edge_1
print(edge_1, edge_2)
n_a_1, n_a_2 = [solution[a_1], solution[a_2]]
print(n_a_1, n_a_2)
b_1, b_2 = edge_2
n_b_1, n_b_2 = [solution[b_1], solution[b_2]]
print(n_b_1, n_b_2)
solution[a_2] = n_b_2
solution[b_2] = n_a_2
solution.append(current_point)
points.remove(current_point)
return solution
points = parse_input_data(input_data)
node_count = len(points)
solution = initial_solution_greedy(list(range(node_count)))
local_search(solution)
# solution = initial_solution_naiv()
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)
ignore_list = []
while True:
points_before_change = list(points)
try:
points = local_search(points, ignore_list)
except UnboundLocalError:
break
is_valid(points, num_nodes)
value = total_length(points)
# plot_graph(points_before_change)
is_valid(points, num_nodes)
return prepare_output_data(points)
if __name__ == "__main__":
file_location = "data/tsp_200_2"
with open(file_location, 'r') as input_data_file:
input_data = input_data_file.read()
print(solve_it(input_data))
plot_solution(solution, points)
is_valid(solution)
return prepare_output_data(solution)