BerGeo/h2/mbc.py

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


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 
    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:
        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
    """

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


    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


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

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

    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

    # 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
    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, [])

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




    # 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, linprog_flipper), {left_point, right_point}, mbc_ch(pr, linprog_flipper))


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