97 lines
3.1 KiB
Python
97 lines
3.1 KiB
Python
import statistics
|
|
from math import inf
|
|
from typing import Set, List, Tuple
|
|
|
|
from util import Side, Point, gen_point, display
|
|
|
|
|
|
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 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
|
|
return Side.BELOW
|
|
|
|
|
|
def solve_1dlp(c: float, constraints: List[Tuple[float, float]]):
|
|
"""
|
|
:param c: c1
|
|
:param constraints: [(ai, bi), ...]
|
|
:return: x1
|
|
"""
|
|
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
|
|
|
|
|
|
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]]):
|
|
"""
|
|
:param c: (c1, c2)
|
|
:param constraints: [(ai1, ai2, bi), ...]
|
|
:return: x1, x2
|
|
"""
|
|
c1, c2 = c
|
|
x1 = -inf if c1 > 0 else inf
|
|
x2 = -inf if c2 > 0 else inf
|
|
|
|
#random.shuffle(constraints)
|
|
|
|
for i, ((a1, a2), b) in enumerate(constraints[1:]):
|
|
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], flipper: callable) -> Set[Point]:
|
|
if len(points) <= 2:
|
|
return points
|
|
|
|
# Find the point with median x-coordinate, and partition the points on this point
|
|
med_x = statistics.median(p.x for p in points)
|
|
|
|
# Find left and right points in regards to median
|
|
pl = {p for p in points if p.x < med_x}
|
|
pr = {p for p in points if p.x >= med_x}
|
|
|
|
# Find the bridge over the vertical line in pm
|
|
slope, intercept = solve_2dlp((flipper(med_x), flipper(1)), # confirmed correct
|
|
[((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) # confirmed correct
|
|
|
|
# Find the two points which are on the line, should work
|
|
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}
|
|
|
|
return set.union(mbc_ch(pl, flipper), {left_point, right_point}, mbc_ch(pr, flipper))
|
|
|
|
|
|
def mbc(points: Set[Point]) -> Set[Point]:
|
|
return set.union(mbc_ch(points, lambda x: x), mbc_ch(points, lambda x: -x))
|
|
|
|
|
|
if __name__ == '__main__':
|
|
points = {gen_point(1, 10) for _ in range(20)}
|
|
|
|
upper_hull_points = mbc_ch(points, lambda x: x)
|
|
lower_hull_points = mbc_ch(points, lambda x: -x)
|
|
|
|
display(points, upper_hull_points.union(lower_hull_points))
|