diff --git a/h2/MBC_carlsen.py b/h2/MBC_carlsen.py deleted file mode 100644 index beb7e55..0000000 --- a/h2/MBC_carlsen.py +++ /dev/null @@ -1,192 +0,0 @@ -import matplotlib.pyplot as plt -import collections -from enum import Enum, auto -from random import randint - -Point = collections.namedtuple("Point", ("x", "y")) - - -class Side(Enum): - ON = auto() - ABOVE = auto() - BELOW = auto() - - -def sidedness(slope, intersection, p3, eps=0.0000001): - print(slope * p3[0] + intersection ) - # finds where a point is in regards to a line - if p3[1] - eps <= slope * p3[0] + intersection <= p3[1] + eps or p3[0] - eps <= (p3[1] - intersection)/slope <= p3[0] + eps: - return Side.ON - elif p3[1] > slope * p3[0] + intersection: - return Side.ABOVE - return Side.BELOW - - -def diplay_prune_points(points, p1, p2): - xs = [p[0] for p in points] - ys = [p[1] for p in points] - - plt.plot(xs, ys, 'ro') - plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) - plt.show() - - -def solve1D(points, xm, iteration_num): - #print("iter:", iteration_num) - point = points[iteration_num] - - if iteration_num == 0: - return -float('Inf'), -float('Inf') - - # lad point[1] = point[0] * a + b <=> y = x * a + b - # isolere b og sæt ind i constraints - a = None - b = point[1] - point[0] # * a - - # minimere xm*a + b, hvor b har ny værdi - # Vi regner kun med koefficienterne - # obj_fun = (xm - point[0]) + points[1] = (xm*a - xi*a) + y - - # max eller min - #print("XM og point[0]:", xm, point[0]) - a = xm - point[0] - #print("a lige her:", a) - - a_constraint_list = [] - # looping over the i constraints - for p in points[:iteration_num]: - # we can't make a straigt vertical line - if p[0] == point[0]: - if p[1] > point[1]: - p[0] = p[0] + 0.0001 - else: - point[0] = point[0] + 0.0001 - - print(p, point) - # Spring over den i'te constraint - if p != point: - # y_j - yi - c = p[1] - point[1] - # x_j * a - xi * a - a_diff = p[0] - point[0] - - print(c, a_diff, c/a_diff) - # det her er forkert og skal fikses - if a >= 0 and a_diff < 0: - a_constraint_list.append(-float('Inf')) - else: - a_constraint_list.append(c/a_diff) - - # hvis a > 0 så min - # hvis a < 0 så max. - if a >= 0: - v1 = max(a_constraint_list) - v2 = point[1] - point[0] * v1 - v = (v1, v2) - elif a < 0: - v1 = min(a_constraint_list) - v2 = point[1] - point[0] * v1 - v = (v1, v2) - - return v - - -def findBridge(points, xm): - if xm > 0: - v = (-float('Inf'), -float('Inf')) - elif xm < 0: - 4 # noget - else: - # TODO: xm == 0 - 4 - # looping over constraints - for point in points: - # checking for violation - if not point[1] <= point[0] * v[0] + v[1]: - v = solve1D(points, xm, points.index(point)) - #print("HER OVER: ", v) - - slope = v[0] - intercept = v[1] - - line_points = [p for p in points if sidedness(slope, intercept, p) == Side.ON] - - if len(line_points) > 3: - print("Halli HAllo") - print(line_points) - - return line_points[0], line_points[1] - - -def find_median(points): - if len(points) % 2 == 0: - first_med_idx = int(len(points) / 2 - 1) - second_med_idx = int(len(points) / 2) - return (points[first_med_idx][0] + points[second_med_idx][0]) / 2 - else: - idx = int((len(points)-1) / 2) - return points[idx][0] - - -def upperHull(points, all_points): - print("Punkter:", points) - if len(points) < 2: - return [] - xm = find_median(points) - - # end-points of bridge - (xi, yi), (xj, yj) = findBridge(points, xm) - - #print(xm) - print((xi, yi), (xj, yj)) - prune_points = [p for p in points if p[0] < xi or xj < p[0]] + [(xi, yi), (xj, yj)] - - # Neden for er mest for visualisering - prune_all_points = [p for p in all_points if p[0] < xi or xj < p[0]] + [(xi, yi), (xj, yj)] - diplay_prune_points(prune_all_points, (xi, yi), (xj, yj)) - - Pl = [p for p in prune_points if p[0] < xm] - Pr = [p for p in prune_points if p[0] >= xm] - #print("Pl:", Pl, "Pr", Pr) - - print("\n") - # recurse results and return - ret = [(xi, yi), (xj, yj)] - ret = ret + [p for p in upperHull(Pl, prune_all_points) if p not in ret] - ret = ret + [p for p in upperHull(Pr, prune_all_points) if p not in ret] - return ret - - - -p1 = (2, 1) -p2 = (3, 4) -p3 = (5, 2) -p4 = (6, 5) - -list_of_points = [] - -list_of_points.append(p1) -list_of_points.append(p2) -list_of_points.append(p3) -list_of_points.append(p4) - - - - -xs = [p[0] for p in list_of_points] -ys = [p[1] for p in list_of_points] - -plt.plot(xs, ys, 'ro') -plt.show() - -result = list(sorted(upperHull(list_of_points, list_of_points))) - -result -res_xs = [p[0] for p in result] -res_ys = [p[1] for p in result] - - -#print("result", result) -plt.plot(res_xs, res_ys) -plt.plot(xs, ys, 'ro') -plt.show() diff --git a/h2/mbc.py b/h2/mbc.py index 51234d4..5526a6b 100644 --- a/h2/mbc.py +++ b/h2/mbc.py @@ -1,12 +1,8 @@ -import statistics -from math import inf, isnan, sqrt -from typing import Set, List, Tuple -import util -from scipy.optimize import linprog -import scipy import random +from math import sqrt +from typing import Set -from util import Side, Point, gen_point, display, display_line_only, gen_circular_point, gen_weird_point, gen_triangular_point +from util import Side, Point, gen_point, display, gen_circular_point, gen_triangular_point def sidedness(slope: float, intersection: float, p3: Point, flipper: callable, eps=0.0000001) -> Side: @@ -18,54 +14,20 @@ def sidedness(slope: float, intersection: float, p3: Point, flipper: callable, e 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_) - elif a < 0: - min_ = max(b/a, min_) - - - if c > 0: - return min_ - else: - return max_ - - - -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 new_solve_1dlp(c, constraints, idx): +def solve_1dlp(c, constraints): c1, c2 = c ((a1, a2), b) = constraints[-1] - q, p = b / a2, a1 /a2 + q, p = b / a2, a1 / a2 - interval = [-9999999, 9999999] + interval = [-10_000, 10_000] - for new_idx, ((lel_a1, lel_a2), lel_b) in enumerate(constraints): + 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 + interval[1] = bj / aj c = -(c1 - c2 * p) if c < 0: @@ -74,65 +36,14 @@ def new_solve_1dlp(c, constraints, idx): return interval[1], q - (p * interval[1]) - -def new_solve_2dlp(c, constraints): +def solve_2dlp(c, constraints): c1, c2 = c - x1 = -10000 if c1 > 0 else 10000 - x2 = -10000 if c2 > 0 else 10000 + x1 = -10_000 if c1 > 0 else 10_000 + x2 = -10_000 if c2 > 0 else 10_000 - for idx, ((a1, a2), b) in enumerate(constraints): + for i, ((a1, a2), b) in enumerate(constraints, start=1): if not (a1*x1 + a2*x2 <= b): - x1,x2 = new_solve_1dlp(c, constraints[:idx+1], idx) - return x1,x2 - - - -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 = -10000 if c1 > 0 else 10000 - x2 = -10000 if c2 > 0 else 10000 - - #random.shuffle(constraints) - - - for idx, ((a1, a2), b) in enumerate(constraints[1:]): - - #print("x1 and x2", x1, x2) - - - #print("{} + {} <= {}".format(a1*x1, a2*x2, b)) - #print("pls",a1*x1 + a2*x2) - - #print("yes"*10) if isnan(a1*x1+a2*x2) else print("no"*10) - - if not (a1*x1 + a2*x2 <= b): - constraint_for_1d = [] - - new_obj = c[0] - ((c[1]*a1)/a2) - - for constraint in constraints[:idx]: - (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 - #print(res) - x2 = ((b/a2) - (a1/a2)*x1) - - + x1, x2 = solve_1dlp(c, constraints[:i]) return x1, x2 @@ -147,72 +58,27 @@ def find_median(points): return median[0] -def mbc_ch(points: Set[Point], flipper: callable) -> Set[Point]: +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 = find_median(points) - #med_x = statistics.median(p.x for p in points) - #print(med_x) # 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} - #print("pl\n",pl) - #print("pr\n",pr) - - c = [flipper(med_x), flipper(1)] - A = [[flipper(-p.x), flipper(-1)] for p in points] - b = [flipper(-p.y) for p in points] # Find the bridge over the vertical line in pm - #slope, intercept = solve_2dlp((flipper(med_x), flipper(1)), - # [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) # confirmed correct - - slope, intercept = new_solve_2dlp((flipper(med_x), flipper(1)), - [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) - - - #print("slope, intercept:",slope, intercept) - - res = linprog(c, A, b, bounds=[[None, None], [None, None]], options={"tol": 0.01}) - #print("res0, res1:",res.x[0], res.x[1]) - - #slope, intercept = res.x[0], res.x[1] - - #display_line_only(points, slope, intercept, []) - - # 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} - #print("On Points:",on) - #left_point = min(on) - #right_point = max(on) + slope, intercept = solve_2dlp((flipper(med_x), flipper(1)), + [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) + # Find the two points which are on the line dist_to_line = lambda p: abs(intercept + slope * p.x - p.y)/sqrt(1 + slope**2) - left_point = min(pl, key = dist_to_line) + left_point = min(pl, key=dist_to_line) right_point = min(pr, key=dist_to_line) - #display_line_only(points, slope, intercept, [left_point, right_point]) - - - #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) - - - - - def find_med_point(points, med_x): - for p in points: - if med_x+0.001 >= p.x >= med_x-0.001: - return {p} - return {} - - #print("med point:",find_med_point(points, med_x)) - - #display_line_only(points, slope, intercept, {left_point, right_point}) - # 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} @@ -227,11 +93,9 @@ def mbc(points: Set[Point]) -> Set[Point]: 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_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)