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

Fix branches.

Browse Source
This commit is contained in:
Casper 2018-10-08 15:17:02 +02:00
parent 0553d80379
commit 08c4b1f790
No known key found for this signature in database
GPG Key ID: B1156723DB3BDDA8
2 changed files with 293 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

387
h2/MBC_carlsen.py
View File

@ -1,193 +1,194 @@
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import collections import collections
from enum import Enum, auto from enum import Enum, auto
from random import randint from random import randint
Point = collections.namedtuple("Point", "x y") Point = collections.namedtuple("Point", "x y")
class Side(Enum): class Side(Enum):
ON = auto() ON = auto()
ABOVE = auto() ABOVE = auto()
BELOW = auto() BELOW = auto()
def sidedness(slope, intersection, p3, eps=0.0000001): def sidedness(slope, intersection, p3, eps=0.0000001):
print(slope * p3[0] + intersection ) print(slope * p3[0] + intersection )
""" finds where a point is in regards to a line """ """ 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: 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 return Side.ON
elif p3[1] > slope * p3[0] + intersection: elif p3[1] > slope * p3[0] + intersection:
return Side.ABOVE return Side.ABOVE
return Side.BELOW return Side.BELOW
def diplay_prune_points(points, p1, p2): def diplay_prune_points(points, p1, p2):
xs = [p[0] for p in points] xs = [p[0] for p in points]
ys = [p[1] for p in points] ys = [p[1] for p in points]
plt.plot(xs, ys, 'ro') plt.plot(xs, ys, 'ro')
plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) plt.plot([p1[0], p2[0]], [p1[1], p2[1]])
plt.show() plt.show()
def solve1D(points, xm, iteration_num): def solve1D(points, xm, iteration_num):
#print("iter:", iteration_num) #print("iter:", iteration_num)
point = points[iteration_num] point = points[iteration_num]
if iteration_num == 0: if iteration_num == 0:
return -float('Inf'), -float('Inf') return -float('Inf'), -float('Inf')
# lad point[1] = point[0] * a + b <=> y = x * a + b # lad point[1] = point[0] * a + b <=> y = x * a + b
# isolere b og sæt ind i constraints # isolere b og sæt ind i constraints
a = None a = None
b = point[1] - point[0] # * a b = point[1] - point[0] # * a
# minimere xm*a + b, hvor b har ny værdi # minimere xm*a + b, hvor b har ny værdi
# Vi regner kun med koefficienterne # Vi regner kun med koefficienterne
# obj_fun = (xm - point[0]) + points[1] = (xm*a - xi*a) + y # obj_fun = (xm - point[0]) + points[1] = (xm*a - xi*a) + y
# max eller min # max eller min
#print("XM og point[0]:", xm, point[0]) #print("XM og point[0]:", xm, point[0])
a = xm - point[0] a = xm - point[0]
#print("a lige her:", a) #print("a lige her:", a)
a_constraint_list = [] a_constraint_list = []
# looping over the i constraints # looping over the i constraints
for p in points[:iteration_num]: for p in points[:iteration_num]:
# we can't make a straigt vertical line # we can't make a straigt vertical line
if p[0] == point[0]: if p[0] == point[0]:
if p[1] > point[1]: if p[1] > point[1]:
p[0] = p[0] + 0.0001 p[0] = p[0] + 0.0001
else: else:
point[0] = point[0] + 0.0001 point[0] = point[0] + 0.0001
print(p, point) print(p, point)
# Spring over den i'te constraint # Spring over den i'te constraint
if p != point: if p != point:
# y_j - yi # y_j - yi
c = p[1] - point[1] c = p[1] - point[1]
# x_j * a - xi * a # x_j * a - xi * a
a_diff = p[0] - point[0] a_diff = p[0] - point[0]
print(c, a_diff, c/a_diff) print(c, a_diff, c/a_diff)
# det her er forkert og skal fikses # det her er forkert og skal fikses
if a >= 0 and a_diff < 0: if a >= 0 and a_diff < 0:
a_constraint_list.append(-float('Inf')) a_constraint_list.append(-float('Inf'))
else: else:
a_constraint_list.append(c/a_diff) a_constraint_list.append(c/a_diff)
# hvis a > 0 så min # hvis a > 0 så min
# hvis a < 0 så max. # hvis a < 0 så max.
if a >= 0: if a >= 0:
v1 = max(a_constraint_list) v1 = max(a_constraint_list)
v2 = point[1] - point[0] * v1 v2 = point[1] - point[0] * v1
v = (v1, v2) v = (v1, v2)
elif a < 0: elif a < 0:
v1 = min(a_constraint_list) v1 = min(a_constraint_list)
v2 = point[1] - point[0] * v1 v2 = point[1] - point[0] * v1
v = (v1, v2) v = (v1, v2)
return v return v
def findBridge(points, xm): def findBridge(points, xm):
if xm > 0: if xm > 0:
v = (-float('Inf'), -float('Inf')) v = (-float('Inf'), -float('Inf'))
elif xm < 0: elif xm < 0:
4 # noget 4 # noget
else: else:
# TODO: xm == 0 # TODO: xm == 0
4 4
# looping over constraints # looping over constraints
for point in points: for point in points:
# checking for violation # checking for violation
if not point[1] <= point[0] * v[0] + v[1]: if not point[1] <= point[0] * v[0] + v[1]:
v = solve1D(points, xm, points.index(point)) v = solve1D(points, xm, points.index(point))
#print("HER OVER: ", v) #print("HER OVER: ", v)
slope = v[0] slope = v[0]
intercept = v[1] intercept = v[1]
line_points = [p for p in points if sidedness(slope, intercept, p) == Side.ON] line_points = [p for p in points if sidedness(slope, intercept, p) == Side.ON]
if len(line_points) > 3: if len(line_points) > 3:
print("Halli HAllo") print("Halli HAllo")
print(line_points) print(line_points)
return line_points[0], line_points[1] return line_points[0], line_points[1]
def find_median(points): def find_median(points):
if len(points) % 2 == 0: if len(points) % 2 == 0:
first_med_idx = int(len(points) / 2 - 1) first_med_idx = int(len(points) / 2 - 1)
second_med_idx = int(len(points) / 2) second_med_idx = int(len(points) / 2)
return (points[first_med_idx][0] + points[second_med_idx][0]) / 2 return (points[first_med_idx][0] + points[second_med_idx][0]) / 2
else: else:
idx = int((len(points)-1) / 2) idx = int((len(points)-1) / 2)
return points[idx][0] return points[idx][0]
def upperHull(points, all_points): def upperHull(points, all_points):
print("Punkter:", points) print("Punkter:", points)
if len(points) < 2: if len(points) < 2:
return [] return []
xm = find_median(points) xm = find_median(points)
# end-points of bridge # end-points of bridge
(xi, yi), (xj, yj) = findBridge(points, xm) (xi, yi), (xj, yj) = findBridge(points, xm)
#print(xm) #print(xm)
print((xi, yi), (xj, yj)) print((xi, yi), (xj, yj))
prune_points = [p for p in points if p[0] < xi or xj < p[0]] + [(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 # 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)] 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)) diplay_prune_points(prune_all_points, (xi, yi), (xj, yj))
Pl = [p for p in prune_points if p[0] < xm] Pl = [p for p in prune_points if p[0] < xm]
Pr = [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("Pl:", Pl, "Pr", Pr)
print("\n") print("\n")
# recurse results and return # recurse results and return
ret = [(xi, yi), (xj, yj)] 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(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] ret = ret + [p for p in upperHull(Pr, prune_all_points) if p not in ret]
return ret return ret
p1 = (2, 1) p1 = (2, 1)
p2 = (3, 4) p2 = (3, 4)
p3 = (5, 2) p3 = (5, 2)
p4 = (6, 5) p4 = (6, 5)
list_of_points = [] list_of_points = []
list_of_points.append(p1) list_of_points.append(p1)
list_of_points.append(p2) list_of_points.append(p2)
list_of_points.append(p3) list_of_points.append(p3)
list_of_points.append(p4) list_of_points.append(p4)
xs = [p[0] for p in list_of_points] xs = [p[0] for p in list_of_points]
ys = [p[1] for p in list_of_points] ys = [p[1] for p in list_of_points]
plt.plot(xs, ys, 'ro') plt.plot(xs, ys, 'ro')
plt.show() plt.show()
result = list(sorted(upperHull(list_of_points, list_of_points))) result = list(sorted(upperHull(list_of_points, list_of_points)))
result result
res_xs = [p[0] for p in result] res_xs = [p[0] for p in result]
res_ys = [p[1] for p in result] res_ys = [p[1] for p in result]
#print("result", result) #print("result", result)
plt.plot(res_xs, res_ys) plt.plot(res_xs, res_ys)
plt.plot(xs, ys, 'ro') plt.plot(xs, ys, 'ro')
plt.show() plt.show()

110
h2/mbc.py
View File

@ -1,21 +1,109 @@
import random
import statistics import statistics
from collections import namedtuple from collections import namedtuple
from typing import List, Set from enum import Enum, auto
from typing import Set
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import linprog
Point = namedtuple('Point', 'x y') Point = namedtuple('Point', 'x y')
def mbc_ch(points: Set[Point]): 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 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 # Find the point with median x-coordinate, and partition the points on this point
pm = statistics.median_high(points) med_x = statistics.median(p.x for p in points)
pl = {point # Find left and right points in regards to median
for point in points pl = {p for p in points if p.x < med_x}
if point.x < pm.x} pr = {p for p in points if p.x >= med_x}
pr = {point
for point in points
if point.x >= pm.x}
# Find the bridge over the vertical line in pm # 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))