Pownie/BerGeo
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merged

Browse Source
This commit is contained in:
Alexander Munch-Hansen 2018-10-18 11:56:16 +02:00
commit 86e2b4c8f2
2 changed files with 21 additions and 344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

192
h2/MBC_carlsen.py
View File

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

173
h2/mbc.py
View File

@ -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,71 +58,31 @@ 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)
#dist_to_line = lambda p: abs(intercept + slope * p.x - p.y)/sqrt(1 + slope**2)
#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])
slope, intercept = solve_2dlp((flipper(med_x), flipper(1)),
[((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points])
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})
# 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)
#right_point = min(pr, key=dist_to_line)
# Prune the points between the two line points
pl = {p for p in pl if p.x <= left_point.x}
@ -227,11 +98,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)