...
This commit is contained in:
parent
d580ea4c89
commit
56c0948ba7
|
|
@ -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()
|
|
||||||
174
h2/mbc.py
174
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
|
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:
|
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
|
return Side.BELOW
|
||||||
|
|
||||||
|
|
||||||
def solve_1dlp(c: float, constraints: List[Tuple[float, float]]):
|
def solve_1dlp(c, constraints):
|
||||||
"""
|
|
||||||
: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):
|
|
||||||
c1, c2 = c
|
c1, c2 = c
|
||||||
((a1, a2), b) = constraints[-1]
|
((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)
|
bj, aj = (lel_b - lel_a2 * q), (lel_a1 - lel_a2 * p)
|
||||||
if aj < 0 and bj / aj > interval[0]:
|
if aj < 0 and bj / aj > interval[0]:
|
||||||
interval[0] = bj / aj
|
interval[0] = bj / aj
|
||||||
elif aj > 0 and bj / aj < interval[1]:
|
elif aj > 0 and bj / aj < interval[1]:
|
||||||
interval[1] = bj/aj
|
interval[1] = bj / aj
|
||||||
|
|
||||||
c = -(c1 - c2 * p)
|
c = -(c1 - c2 * p)
|
||||||
if c < 0:
|
if c < 0:
|
||||||
|
|
@ -74,65 +36,14 @@ def new_solve_1dlp(c, constraints, idx):
|
||||||
return interval[1], q - (p * interval[1])
|
return interval[1], q - (p * interval[1])
|
||||||
|
|
||||||
|
|
||||||
|
def solve_2dlp(c, constraints):
|
||||||
def new_solve_2dlp(c, constraints):
|
|
||||||
c1, c2 = c
|
c1, c2 = c
|
||||||
x1 = -10000 if c1 > 0 else 10000
|
x1 = -10_000 if c1 > 0 else 10_000
|
||||||
x2 = -10000 if c2 > 0 else 10000
|
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):
|
if not (a1*x1 + a2*x2 <= b):
|
||||||
x1,x2 = new_solve_1dlp(c, constraints[:idx+1], idx)
|
x1, x2 = solve_1dlp(c, constraints[:i])
|
||||||
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)
|
|
||||||
|
|
||||||
|
|
||||||
return x1, x2
|
return x1, x2
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -147,72 +58,27 @@ def find_median(points):
|
||||||
|
|
||||||
return median[0]
|
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:
|
if len(points) < 2:
|
||||||
return points
|
return points
|
||||||
|
|
||||||
# Find the point with median x-coordinate, and partition the points on this point
|
# Find the point with median x-coordinate, and partition the points on this point
|
||||||
med_x = find_median(points)
|
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
|
# Find left and right points in regards to median
|
||||||
pl = {p for p in points if p.x < med_x}
|
pl = {p for p in points if p.x < med_x}
|
||||||
pr = {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
|
# Find the bridge over the vertical line in pm
|
||||||
#slope, intercept = solve_2dlp((flipper(med_x), flipper(1)),
|
slope, intercept = solve_2dlp((flipper(med_x), flipper(1)),
|
||||||
# [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points]) # confirmed correct
|
[((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points])
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
|
# 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)
|
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)
|
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
|
# Prune the points between the two line points
|
||||||
pl = {p for p in pl if p.x <= left_point.x}
|
pl = {p for p in pl if p.x <= left_point.x}
|
||||||
pr = {p for p in pr if p.x >= right_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__':
|
if __name__ == '__main__':
|
||||||
random.seed(1337_420)
|
random.seed(1337_420)
|
||||||
points = {gen_point(0, 20) for _ in range(20)}
|
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 = {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=-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)}
|
#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)
|
upper_hull_points = mbc_ch(points, lambda x: x)
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue
Block a user