## helper functions
# return a list of all the k subsets of a given list
def k_subsets(list,k):
if len(list) < k or k <= 0:
return []
elif k == 1: # treating this separately makes things simpler
return [[x] for x in list]
else:
subsets = []
# select a first element, and then find all the k-1 subsets
# from the later elements, and prepend your chosen first element
for i,elem in enumerate(list):
for smaller_subset in k_subsets(list[i+1:], k - 1):
subsets.append([elem] + smaller_subset)
return subsets
# compute the given metric for two points (in any dimension > 1)
def metric(m,p1,p2,label_too=False):
if len(p1) != len(p2):
return 0
# square of the horizontal distance
horiz = sum( [ (a-b)^2 for a,b in zip(p1,p2)[:-1] ] )
# the vertical distance
vert = abs(p1[-1] - p2[-1])
# rephrasing of the condition for distinguishing the types of edges
test = vert^2 - m^2 * horiz
if test > 0:
if label_too:
return 'b', sqrt(1+1/m^2)*vert
else:
return sqrt(1+1/m^2)*vert
elif label_too:
if test == 0:
return 'm', sqrt(horiz + vert^2)
else:
return 'f', sqrt(horiz + vert^2)
else:
return sqrt(horiz + vert^2)
### plots the graph nicely
## this needs to be rewritten
def graph_plot(g,pos,special,m,edges={'b':{'color':'red'},'m':{'color':'orange'},'f':{'color':'green'}},special_vertex={'color':'red','size':8},normal_vertex={'size':6}):
e = g.edges()
default = {}
if 'all' in edges:
default = edges['all']
if 'b' not in edges:
edges['b'] = default
if 'm' not in edges:
edges['m'] = default
if 'f' not in edges:
edges['f'] = default
g = Graphics()
points = []
for i,j,w in e:
label,d = metric(m,pos[i],pos[j],label_too=True)
g += line3d([pos[i],pos[j]],**edges[label])
if i not in points:
points.append(i)
if i in special:
g += point3d(pos[i],**special_vertex)
else:
g += point3d(pos[i],**normal_vertex)
if j not in points:
points.append(j)
if j in special:
g += point3d(pos[j],**special_vertex)
else:
g += point3d(pos[j],**normal_vertex)
return g
def average(list):
return sum(list)/len(list)
# can add an arbitrary number of points
# does it by grouping the coefficients by index
# and then summing them.
def add_points(*points):
return [sum(row) for row in zip(*points)]
def scale_point(s,p1):
return [s*x for x in p1]
# computes the total distance of one point to any number of others
def total_distance(m,p,*other_points):
return sum([ metric(m,p,other_point,label_too=False) for other_point in other_points ])
# find/approximate a steiner point of some vertices returning (True, p)
# iff p is the steiner point and is different to all of them
def find_steiner_point(m,*ps):
if len(ps) == 0:
return False, None
l = len(ps[0])
if not all([len(x) == l for x in ps]):
return False,None
# Finds an alright guess, i.e. the geometric centriod of the points,
# and then optimises that guess by adding a small amount here and there to try to improve it.
centroid = scale_point(1/len(ps),add_points(*ps))
centroid_dist = total_distance(m,centroid,*ps)
# initial conditions
guess = centroid
dist = centroid_dist
step = dist / 10 # 10 is just for fun
while step > 5e-4 * centroid_dist:
best = guess
best_dist = dist
# this is our permuting stage, so it adds a small amount in each direction (6 of them)
# and checks to see if these new points are better
for twiddler in [(step,0,0),(0,step,0),(0,0,step)]:
neg_twiddler = scale_point(-1,twiddler)
guess1 = add_points(guess, twiddler)
guess2 = add_points(guess, neg_twiddler)
dist1 = total_distance(m,guess1,*ps)
dist2 = total_distance(m,guess2,*ps)
# do the checks
if dist1 < best_dist:
best = guess1
best_dist = dist1
if dist2 < best_dist:
best = guess2
best_dist = dist2
# we've found our new winners, so remember them for next time
guess = best
dist = best_dist
step /= 2
# now check that the steiner point isn't too close to one of the terminals
# if it is, we will assume that the terminal is the steiner point
threshold = 1e-2 * centroid_dist
for p in ps:
if metric(m,p,guess) < threshold:
return False,p
# if i'm here, then guess passed all the tests, so return it
return True,guess
def get_all_3_steiner_points(ps):
v = []
# do the 3 points
for set_of_points in k_subsets(ps,3):
is_point,real_point = find_steiner_point(m, *set_of_points)
if is_point:
v.append(real_point)
return v
def get_all_4_steiner_points(ps):
v = []
for set_of_points in k_subsets(ps,4):
is_point,real_point = find_steiner_point(m, *set_of_points)
if is_point:
v.append(real_point)
return v