From 26dfb6d305c32f141682c3320314deeadced3965 Mon Sep 17 00:00:00 2001 From: "Casper V. Kristensen" Date: Tue, 9 Oct 2018 19:26:55 +0200 Subject: [PATCH] Cleanup. --- h2/gift_wrapper.py | 31 +----- h2/{sidedness.py => graham.py} | 160 ++++++++++++------------------ h2/mbc.py | 172 +++++---------------------------- h2/quick_hull.py | 36 +------ h2/util.py | 58 +++++++++++ 5 files changed, 146 insertions(+), 311 deletions(-) rename h2/{sidedness.py => graham.py} (60%) create mode 100644 h2/util.py diff --git a/h2/gift_wrapper.py b/h2/gift_wrapper.py index edc28dd..e07a2c0 100644 --- a/h2/gift_wrapper.py +++ b/h2/gift_wrapper.py @@ -1,21 +1,7 @@ # Use atan2 instead of acos to calc angle; atan2(x,y) of the point we potentially want to add -import random -from collections import namedtuple -import matplotlib.pyplot as plt from math import acos, sqrt -Point = namedtuple('Point', 'x y') -Vector = namedtuple('Vector', 'x y') - - -def gen_point(): - a = random.uniform(1, 5) - b = random.uniform(1, 5) - - x_i = random.uniform(1, 5) - p_i = Point(x_i, a * x_i + b) - - return p_i +from util import Vector, Point, gen_point, display def calc_angle(v1: Vector, v2: Vector) -> float: @@ -31,25 +17,10 @@ def calc_angle(v1: Vector, v2: Vector) -> float: return acos(hmmmmmmm) - - def calc_vector(p1: Point, p2: Point) -> Vector: return Vector((p2.x - p1.x), (p2.y - p1.y)) -def display(points, hull): - x = [point.x for point in points] - y = [point.y for point in points] - - h_x = [point.x for point in hull] - h_y = [point.y for point in hull] - - plt.plot(h_x, h_y, 'ro-') - - plt.scatter(x, y) - plt.show() - - def rapper(points: set): min_pt = min(points) hull = [min_pt] diff --git a/h2/sidedness.py b/h2/graham.py similarity index 60% rename from h2/sidedness.py rename to h2/graham.py index 40ddc25..41f1f0a 100644 --- a/h2/sidedness.py +++ b/h2/graham.py @@ -1,97 +1,63 @@ -import random -from collections import namedtuple -from enum import Enum, auto -import matplotlib.pyplot as plt - -Point = namedtuple('Point', ['x', 'y']) - - -class Side(Enum): - ON = auto() - ABOVE = auto() - BELOW = auto() - - -def sidedness(p1, p2, p3, eps=0.0000001): - - # Find line from p1 to p2, ask where p3 is in relation to this - - y = p3.y * (p2.x - p1.x) - x = p3.x - - a = (p2.y - p1.y) - b = p2.y * (p2.x - p1.x) - a * p2.x - - if y - eps < a * x + b < y + eps: - return Side.ON - elif y > a * x + b: - return Side.ABOVE - return Side.BELOW - - -# test - - -p1 = Point(4, 4) -p2 = Point(0, 0) -p3 = Point(5, 2) - - -# print(sidedness(p1, p2, p3)) - - -def genPoint(): - a = random.uniform(1, 5) - b = random.uniform(1, 5) - - x_i = random.uniform(1, 5) - p_i = Point(x_i, a * x_i + b) - - return p_i - - -def graham_scan(points): - - # A funky issue where both a and b become negative in the sidedness test causes us to have to use - # Side.ABOVE for both tests, regardless of UH or LH. - - sorted_points = sorted(points) - - UH = sorted_points[:2] - #del sorted_points[0] - - for s in sorted_points[2:]: - while len(UH) > 1 and (sidedness(UH[-2], UH[-1], s) != Side.ABOVE): - del UH[-1] - UH.append(s) - - reversed_list = list(reversed(sorted_points)) - reversed_list.append(UH[0]) - LH = reversed_list[:2] - #del reversed_list[0] - - for s in reversed_list[2:]: - while len(LH) > 1 and (sidedness(LH[-2], LH[-1], s) != Side.ABOVE): - del LH[-1] - LH.append(s) - - return UH, LH - - -p = [genPoint() for _ in range(30)] -UH, LH = graham_scan(p) - -x = [point.x for point in p] -y = [point.y for point in p] - -UH_x = [point.x for point in UH] -UH_y = [point.y for point in UH] - -LH_x = [point.x for point in LH] -LH_y = [point.y for point in LH] - -plt.plot(UH_x, UH_y, 'ro-') -plt.plot(LH_x, LH_y, 'ro-') - -plt.scatter(x,y) -plt.show() +from util import gen_point, Side, Point, display + + +def sidedness(p1, p2, p3, eps=0.0000001): + + # Find line from p1 to p2, ask where p3 is in relation to this + + y = p3.y * (p2.x - p1.x) + x = p3.x + + a = (p2.y - p1.y) + b = p2.y * (p2.x - p1.x) - a * p2.x + + if y - eps < a * x + b < y + eps: + return Side.ON + elif y > a * x + b: + return Side.ABOVE + return Side.BELOW + + +# test + + +p1 = Point(4, 4) +p2 = Point(0, 0) +p3 = Point(5, 2) + + +# print(sidedness(p1, p2, p3)) + + +def graham_scan(points): + + # A funky issue where both a and b become negative in the sidedness test causes us to have to use + # Side.ABOVE for both tests, regardless of UH or LH. + + sorted_points = sorted(points) + + UH = sorted_points[:2] + #del sorted_points[0] + + for s in sorted_points[2:]: + while len(UH) > 1 and (sidedness(UH[-2], UH[-1], s) != Side.ABOVE): + del UH[-1] + UH.append(s) + + reversed_list = list(reversed(sorted_points)) + reversed_list.append(UH[0]) + LH = reversed_list[:2] + #del reversed_list[0] + + for s in reversed_list[2:]: + while len(LH) > 1 and (sidedness(LH[-2], LH[-1], s) != Side.ABOVE): + del LH[-1] + LH.append(s) + + return UH, LH + + +p = [gen_point() for _ in range(30)] +UH, LH = graham_scan(p) + +display(p, {*UH, *LH}) diff --git a/h2/mbc.py b/h2/mbc.py index 71bc11e..16f2335 100644 --- a/h2/mbc.py +++ b/h2/mbc.py @@ -1,67 +1,12 @@ -import random import statistics -from collections import namedtuple -from enum import Enum, auto from math import inf from typing import Set, List, Tuple -import matplotlib.pyplot as plt -import numpy as np -from scipy.optimize import linprog - -Point = namedtuple('Point', 'x y') +from util import Side, Point, gen_point, display -def gen_point(lower: int, upper: int) -> Point: - a = random.uniform(lower, upper) - b = random.uniform(lower, upper) - - x_i = random.uniform(lower, upper) - p_i = Point(x_i, a * x_i + b) - p_i = Point(a, b) - return p_i - - -def display(points: Set[Point], hull: Set[Point]): - x = [point.x for point in points] - y = [point.y for point in points] - - h_x = [point.x for point in hull] - h_y = [point.y for point in hull] - - plt.plot(h_x, h_y, 'ro') - - plt.scatter(x, y) - plt.show() - - -def display_line_only(points: Set[Point], slope: int, intercept: int, line_points: Set[Point]): - x = [point.x for point in points] - y = [point.y for point in points] - - plt.scatter(x, y) - - # Plot a line from slope and intercept - axes = plt.gca() - x_vals = np.array(axes.get_xlim()) - y_vals = intercept + slope * x_vals - - for point in line_points: - plt.plot(point.x, point.y, 'go') - - plt.plot(x_vals, y_vals, '--') - - plt.show() - - -class Side(Enum): - ON = auto() - ABOVE = auto() - BELOW = auto() - - -def sidedness(slope: int, intersection: int, p3: Point, linprog_flipper: callable, eps=0.0000001) -> Side: - # finds where a point is in regards to a line +def sidedness(slope: float, intersection: float, p3: Point, linprog_flipper: callable, eps=0.0000001) -> Side: + # finds where a point is in regards to a line if linprog_flipper(p3.y) - eps <= linprog_flipper(slope * p3.x + intersection) <= linprog_flipper(p3.y) + eps: return Side.ON elif p3.y > slope * p3.x + intersection: @@ -75,38 +20,19 @@ def solve_1dlp(c: float, constraints: List[Tuple[float, float]]): :param constraints: [(ai, bi), ...] :return: x1 """ - - min_ = -10000 - max_ = 10000 - - for constraint in constraints: - (a, b) = constraint - - if a == 0: - assert(b >= 0) - - if a > 0: - max_ = min(b/a, max_) - else: - min_ = max(b/a, min_) + try: + if c > 0: + return max(b/a for a, b in constraints if a < 0) + return min(b/a for a, b in constraints if a > 0) + except ValueError: # unbounded + return -inf if c > 0 else inf - if c > 0: - return min_ - else: - return max_ - - - # if c > 0: - # return max(b/a for a, b in constraints if a < 0) - # return min(b/a for a, b in constraints if a > 0) - - -#assert solve_1dlp(1, [(-1, -2)]) == 2 -#assert solve_1dlp(1, [(-1, -2), (-1, -3)]) == 3 -#assert solve_1dlp(1, [(-1, -3), (-1, -2)]) == 3 -#assert solve_1dlp(-1, [(1, 3), (1, 2)]) == 2 -#assert solve_1dlp(1, [(-1, 3), (-1, 2)]) == -2 +assert solve_1dlp(1, [(-1, -2)]) == 2 +assert solve_1dlp(1, [(-1, -2), (-1, -3)]) == 3 +assert solve_1dlp(1, [(-1, -3), (-1, -2)]) == 3 +assert solve_1dlp(-1, [(1, 3), (1, 2)]) == 2 +assert solve_1dlp(1, [(-1, 3), (-1, 2)]) == -2 def solve_2dlp(c: Tuple[float, float], constraints: List[Tuple[Tuple[float, float], float]]): @@ -115,69 +41,20 @@ def solve_2dlp(c: Tuple[float, float], constraints: List[Tuple[Tuple[float, floa :param constraints: [(ai1, ai2, bi), ...] :return: x1, x2 """ + c1, c2 = c + x1, x2 = (-inf, -inf) if c1 > 0 else (inf, inf) + #random.shuffle(constraints) - if c[0] > 0: - x1 = -10000 - else: - x1 = 10000 - - if c[1] > 0: - x2 = -10000 - else: - x2 = 10000 - - our_constraints = [] - for (a1, a2), b in constraints: # TODO: random.shuffle() - - print("x1 and x2", x1, x2) - - if len(our_constraints) == 0: - our_constraints.append(((a1, a2), b)) + for i, ((a1, a2), b) in enumerate(constraints[1:], start=1): + if a1*x1 + a2*x2 <= b: continue - print("{} + {} <= {}".format(a1*x1, a2*x2, b)) - if not (a1*x1 + a2*x2 <= b): - - constraint_for_1d = [] - - new_obj = c[0] - ((c[1]*a1)/a2) - - for constraint in our_constraints: - (a_i1, a_i2), b_i = constraint - - a_prime = a_i1 - ((a_i2*a1)/a2) - b_prime = b_i - ((a_i2*b)/a2) - constraint_for_1d.append((a_prime, b_prime)) - - print("obj", new_obj) - print("const", constraint_for_1d) - - print("lol",[cons[0] for cons in constraint_for_1d]) - # res = linprog([new_obj], [[cons[0]] for cons in constraint_for_1d], [[cons[1]] for cons in constraint_for_1d], bounds=[(None, None)]) - x1 = solve_1dlp(new_obj, constraint_for_1d) - # x1 = res.x[0] - x2 = ((b/a2) - (a1/a2)*x1) - - our_constraints.append(((a1, a2), b)) + x1 = solve_1dlp(c1 - c2*a1/a2, [(ai1 - ai2*a1 / a2, bi - ai2*b / a2) for (ai1, ai2), bi in constraints[:i]]) + x2 = (b - a1*x1) / a2 return x1, x2 -# assert solve_2dlp((-3, -2), [((-1, 3), 12), ((1, 1), 8), ((2, -1), 10)]) == (6, 2) - -c = (-3, -2) -constraints = [ - ((-1, 3), 12), - ((1, 1), 8), - ((2, -1), 10), -] -# = (6.0, 2.0) - -result = solve_2dlp(c, constraints) -print(result) -#exit() - - def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]: if len(points) <= 2: return points @@ -190,8 +67,8 @@ def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]: pr = {p for p in points if p.x >= med_x} # Find the bridge over the vertical line in pm - c = [linprog_flipper(med_x), linprog_flipper(1)] - A = [[linprog_flipper(-p.x), linprog_flipper(-1)] for p in points] + c = (linprog_flipper(med_x), linprog_flipper(1)) + A = [(linprog_flipper(-p.x), linprog_flipper(-1)) for p in points] b = [linprog_flipper(-p.y) for p in points] #result = linprog(c, A, b, bounds=[(None, None), (None, None)], options={"tol": 0.00001}) @@ -205,9 +82,6 @@ def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]: left_point = next(p for p in pl if sidedness(slope, intercept, p, linprog_flipper) == Side.ON) right_point = next(p for p in pr if sidedness(slope, intercept, p, linprog_flipper) == Side.ON) - - - # Prune the points between the two line points pl = {p for p in pl if p.x <= left_point.x} pr = {p for p in pr if p.x >= right_point.x} diff --git a/h2/quick_hull.py b/h2/quick_hull.py index bd5f346..7b1c105 100644 --- a/h2/quick_hull.py +++ b/h2/quick_hull.py @@ -1,41 +1,7 @@ -import random -from collections import namedtuple -from enum import Enum, auto from math import sqrt from typing import Set -import matplotlib.pyplot as plt - -Point = namedtuple('Point', 'x y') - - -class Side(Enum): - ON = auto() - ABOVE = auto() - BELOW = auto() - - -def display(points, hull): - x = [point.x for point in points] - y = [point.y for point in points] - - h_x = [point.x for point in hull] - h_y = [point.y for point in hull] - - plt.plot(h_x, h_y, 'ro') - - plt.scatter(x, y) - plt.show() - - -def gen_point(): - a = random.uniform(1, 5) - b = random.uniform(1, 5) - - x_i = random.uniform(1, 5) - p_i = Point(x_i, a * x_i + b) - - return p_i +from util import Point, gen_point, display def distance(a, b, c): diff --git a/h2/util.py b/h2/util.py new file mode 100644 index 0000000..837d826 --- /dev/null +++ b/h2/util.py @@ -0,0 +1,58 @@ +import random +from collections import namedtuple +from enum import Enum, auto +from typing import Set + +import matplotlib.pyplot as plt +import numpy as np + +Point = namedtuple('Point', 'x y') +Vector = namedtuple('Vector', 'x y') + + +def gen_point(lower: int = 0, upper: int = 10) -> Point: + a = random.uniform(lower, upper) + b = random.uniform(lower, upper) + + x_i = random.uniform(lower, upper) + p_i = Point(x_i, a * x_i + b) + p_i = Point(a, b) + return p_i + + +def display(points: Set[Point], hull: Set[Point]): + x = [point.x for point in points] + y = [point.y for point in points] + + h_x = [point.x for point in hull] + h_y = [point.y for point in hull] + + plt.plot(h_x, h_y, 'ro') + + plt.scatter(x, y) + plt.show() + + +def display_line_only(points: Set[Point], slope: int, intercept: int, line_points: Set[Point]): + x = [point.x for point in points] + y = [point.y for point in points] + + plt.scatter(x, y) + + # Plot a line from slope and intercept + axes = plt.gca() + x_vals = np.array(axes.get_xlim()) + y_vals = intercept + slope * x_vals + + for point in line_points: + plt.plot(point.x, point.y, 'go') + + plt.plot(x_vals, y_vals, '--') + + plt.show() + + +class Side(Enum): + ON = auto() + ABOVE = auto() + BELOW = auto()