From c6bfb7d30ecea137655a24569ad54349736f3a24 Mon Sep 17 00:00:00 2001 From: "Casper V. Kristensen" Date: Thu, 11 Oct 2018 15:54:52 +0200 Subject: [PATCH] Lol --- h2/gift_wrapper.py | 5 ++--- h2/mbc.py | 44 ++++++++++++++++++++------------------------ h2/tmptest.py | 26 ++++++++++++++++++-------- h2/util.py | 2 +- 4 files changed, 41 insertions(+), 36 deletions(-) diff --git a/h2/gift_wrapper.py b/h2/gift_wrapper.py index 541fd3c..8670a98 100644 --- a/h2/gift_wrapper.py +++ b/h2/gift_wrapper.py @@ -1,7 +1,7 @@ # Use atan2 instead of acos to calc angle; atan2(x,y) of the point we potentially want to add from math import acos, sqrt -from util import Vector, Point, gen_point, display +from util import Vector, Point, gen_point def calc_angle(v1: Vector, v2: Vector) -> float: @@ -10,8 +10,7 @@ def calc_angle(v1: Vector, v2: Vector) -> float: len_2 = sqrt(v2.x**2 + v2.y**2) tmp = dot / (len_1 * len_2) - - return acos(round(tmp, 6)) + return acos(max(min(tmp, 1), -1)) # acos is only defined in [-1,1] def calc_vector(p1: Point, p2: Point) -> Vector: diff --git a/h2/mbc.py b/h2/mbc.py index 4a50ace..b0758a6 100644 --- a/h2/mbc.py +++ b/h2/mbc.py @@ -1,13 +1,17 @@ +import random import statistics from math import inf from typing import Set, List, Tuple +from scipy.optimize import linprog + +import util from util import Side, Point, gen_point, display -def sidedness(slope: float, intersection: float, p3: Point, linprog_flipper: callable, eps=0.0000001) -> Side: +def sidedness(slope: float, intersection: float, p3: Point, 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: + if flipper(p3.y) - eps <= flipper(slope * p3.x + intersection) <= flipper(p3.y) + eps: return Side.ON elif p3.y > slope * p3.x + intersection: return Side.ABOVE @@ -45,17 +49,16 @@ def solve_2dlp(c: Tuple[float, float], constraints: List[Tuple[Tuple[float, floa x1, x2 = (-inf, -inf) if c1 > 0 else (inf, inf) #random.shuffle(constraints) - for i, ((a1, a2), b) in enumerate(constraints[1:], start=1): - if a1*x1 + a2*x2 <= b: - continue - - 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 + for i, ((a1, a2), b) in enumerate(constraints[1:], start=0): + if not a1*x1 + a2*x2 <= 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 -def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]: +def mbc_ch(points: Set[Point], flipper: callable) -> Set[Point]: if len(points) <= 2: return points @@ -67,26 +70,19 @@ 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] - b = [linprog_flipper(-p.y) for p in points] - - #result = linprog(c, A, b, bounds=[(None, None), (None, None)], options={"tol": 0.00001}) - #slope, intercept = result.x[0], result.x[1] - - slope, intercept = solve_2dlp(c, list(zip(A, b))) - - #display_line_only(points, slope, intercept, []) + slope, intercept = solve_2dlp((flipper(med_x), flipper(1)), + [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) # Find the two points which are on the line, should work - 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) + on = {p for p in points if sidedness(slope, intercept, p, flipper) == Side.ON} + left_point = min(on) + right_point = max(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} + pl = {p for p in points if p.x <= left_point.x} + pr = {p for p in points if p.x >= right_point.x} - return set.union(mbc_ch(pl, linprog_flipper), {left_point, right_point}, mbc_ch(pr, linprog_flipper)) + return set.union(mbc_ch(pl, flipper), {left_point, right_point}, mbc_ch(pr, flipper)) def mbc(points: Set[Point]) -> Set[Point]: diff --git a/h2/tmptest.py b/h2/tmptest.py index fecd4ee..e75124c 100644 --- a/h2/tmptest.py +++ b/h2/tmptest.py @@ -7,16 +7,26 @@ from graham import graham_scan from mbc import mbc from quick_hull import quick_hull -random.seed(1337_420) +#random.seed(1337_420) -points = {util.gen_point(-100, 100) for i in range(5_000)} +while True: + print("="*100) + points = {util.gen_point(-100, 100) for i in range(100)} + print("Points:", sorted(points)) -# Sanity check -graham = set(graham_scan(points)) -gift = set(rapper(points)) -quick = quick_hull(points) -mbch = set.union(mbc(points)) -assert gift == graham == quick == mbch + # Sanity check + graham = set(graham_scan(points)) + gift = set(rapper(points)) + quick = quick_hull(points) + mbch = set.union(mbc(points)) + + print("graham:", sorted(graham)) + print("gift :", sorted(gift)) + print("quick :", sorted(quick)) + print("mbch :", sorted(mbch)) + + if not gift == graham == quick == mbch: + util.display(points=points, hull=mbch) def time_it(f: callable, args: tuple = (), iterations=20): diff --git a/h2/util.py b/h2/util.py index ed72b21..3fff2e8 100644 --- a/h2/util.py +++ b/h2/util.py @@ -32,7 +32,7 @@ def display(points: Set[Point], hull: Set[Point]): plt.show() -def display_line_only(points: Set[Point], slope: int, intercept: int, line_points: Set[Point]): +def display_line_only(points: Set[Point], slope: float, intercept: float, line_points: Set[Point]): x = [point.x for point in points] y = [point.y for point in points]