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: print("rrgeerg"*30) constraint_for_1d = [] # a_prime = a1 - (a2*a1)/a2 # b_prime = b - (a2*b)/a2 # print("a,b prime", a_prime, b_prime) # constraint_for_1d.append((a_prime, b_prime)) # Fix this new_obj = c[0] - (c[1]*a1)/a2 for constraint in our_constraints: (a_i1, a_i2), b_i = constraint print("lel", a_i1, a_i2, b_i) 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("lal", a1, a2, b) print(new_obj) print(constraint_for_1d) x1 = solve_1dlp(new_obj, constraint_for_1d) 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] # 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))