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
    """
    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
    """
    x1, x2 = -inf, -inf  # TODO: something other than -inf (?)
    our_constraints = []
    for (a1, a2), b in constraints:  # TODO: random.shuffle()
        if not a1*x1 + a2*x2 <= b:
            our_constraints.append(((a1, a2), b))

            new_obj = c[0] + (b/a2 - a1/a2)

            constraint_for_1d = []
            for constraint in our_constraints:
                (old_a1, old_a2), old_b = constraint
                new_a1 = old_a1 + ((old_b/old_a2) - (old_a1/old_a2))

                constraint_for_1d.append((new_a1, old_b))

            x1 = solve_1dlp(new_obj, constraint_for_1d)
            x2 = ((b/a2) - (a1/a2))*x1
    return x1, x2


# TODO: 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]

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