BerGeo/h2/mbc.py

import random
from typing import Set

from profile import Profiler
from util import Side, Point, gen_point, display, gen_circular_point


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, constraints):
    c1, c2 = c
    ((a1, a2), b) = constraints[-1]
    q, p = b / a2, a1 / a2

    interval = [-10_000, 10_000]

    for (lel_a1, lel_a2), lel_b in constraints:
        bj, aj = (lel_b - lel_a2 * q), (lel_a1 - lel_a2 * p)
        if aj < 0 and bj / aj > interval[0]:
            interval[0] = bj / aj
        elif aj > 0 and bj / aj < interval[1]:
            interval[1] = bj / aj

    c = -(c1 - c2 * p)
    if c < 0:
        return interval[0], q - (p * interval[0])
    elif c >= 0:
        return interval[1], q - (p * interval[1])


@Profiler("solving LP")
def solve_2dlp(c, constraints):
    c1, c2 = c
    x1 = -10_000 if c1 > 0 else 10_000
    x2 = -10_000 if c2 > 0 else 10_000

    for i, ((a1, a2), b) in enumerate(constraints, start=1):
        if not (a1*x1 + a2*x2 <= b):
            x1, x2 = solve_1dlp(c, constraints[:i])
    return x1, x2


@Profiler("finding median")
def find_median(points):
    num_candidates = min(5, len(points))
    candidates = random.sample(points, num_candidates)

    candidates.sort(key=lambda p: p.x)
    median_i = num_candidates // 2
    median = ((candidates[median_i - 1].x + candidates[median_i].x)/2,
              (candidates[median_i - 1].y + candidates[median_i].y)/2)

    return median[0]


def is_left(a: Point, b: Point, c: Point):
    return ((b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)) > 0


def mbc_ch(points: Set[Point], flipper: callable, extra_prune=False, shuffle=True) -> Set[Point]:
    if len(points) < 2:
        return points

    # Extra pruning step
    if extra_prune:
        with Profiler("extra pruning step"):
            left_point = min(points, key=lambda p: (p.x, -p.y))
            right_point = max(points, key=lambda p: (p.x, -p.y))
            if flipper(1) == 1:
                points = {p for p in points if is_left(left_point, right_point, p)}.union({left_point, right_point})
            else:
                points = {p for p in points if not is_left(left_point, right_point, p)}.union({left_point, right_point})

    # Find the point with median x-coordinate, and partition the points on this point
    med_x = find_median(points)

    # Find left and right points in regards to median
    with Profiler("partitioning set"):
        pl = {p for p in points if p.x < med_x}
        pr = {p for p in points if p.x >= med_x}

    # Shuffle
    with Profiler("flipping constraints"):
        constraints = [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]
    if shuffle:
        with Profiler("shuffling constraints"):
            random.shuffle(constraints)

    # Find the bridge over the vertical line in pm
    slope, intercept = solve_2dlp((flipper(med_x), flipper(1)), constraints)

    with Profiler("finding bridge points"):
        left_point = next(p for p in pl if sidedness(slope, intercept, p, flipper) == Side.ON)
        right_point = next(p for p in pr if sidedness(slope, intercept, p, flipper) == Side.ON)

    # Prune the points between the two line points
    with Profiler("pruning between 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, extra_prune=extra_prune, shuffle=shuffle),
                     {left_point, right_point},
                     mbc_ch(pr, flipper, extra_prune=extra_prune, shuffle=shuffle))


@Profiler("mbc")
def mbc(points: Set[Point]) -> Set[Point]:
    return set.union(mbc_ch(points, lambda x: x, extra_prune=False, shuffle=True),
                     mbc_ch(points, lambda x: -x, extra_prune=False, shuffle=True))


@Profiler("mbc2")
def mbc2(points: Set[Point]) -> Set[Point]:
    return set.union(mbc_ch(points, lambda x: x, extra_prune=True, shuffle=True),
                     mbc_ch(points, lambda x: -x, extra_prune=True, shuffle=True))


@Profiler("mbc_no_shuffle")
def mbc_no_shuffle(points: Set[Point]) -> Set[Point]:
    return set.union(mbc_ch(points, lambda x: x, extra_prune=False, shuffle=False),
                     mbc_ch(points, lambda x: -x, extra_prune=False, shuffle=False))


@Profiler("mbc2_no_shuffle")
def mbc2_no_shuffle(points: Set[Point]) -> Set[Point]:
    return set.union(mbc_ch(points, lambda x: x, extra_prune=True, shuffle=False),
                     mbc_ch(points, lambda x: -x, extra_prune=True, shuffle=False))


if __name__ == '__main__':
    random.seed(1337_420)
    points = {gen_point(0, 20) for _ in range(20)}
    points = {gen_circular_point(1, 100, 50) for _ in range(200)}
    #points = {gen_triangular_point(Point(1,1), Point(51,1), Point(26, 30)) for _ in range(200)}
    #points = {Point(x=-33.11091053638924, y=38.88967778961347), Point(x=61.20269947424177, y=-78.96305419217254), Point(x=99.44053842147957, y=-89.11579172297581), Point(x=-92.40380889537532, y=84.33904351991652), Point(x=-90.63139185545595, y=-91.13793846505985)}
    #points = {Point(x=5.2, y=9.7), Point(x=5.3, y=4.9), Point(x=3.3, y=3.6), Point(x=9.2, y=4.8), Point(x=9.7, y=5.7), Point(x=5.6, y=8.7)}

    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))