2748-Mine Optimisation & Steiner Trees.sagews

%html
<h1>Mine Optimisation & Steiner Trees</h1>
<p>The following is program to compute approximate steiner trees for a gradient restricted network, such as a mine</p>

Mine Optimisation & Steiner Trees

The following is program to compute approximate steiner trees for a gradient restricted network, such as a mine

## The terminals of the steiner tree (e.g. ore bodies, surface outlets)
nodes=[(0,0,0),(0,1,-0.25),(0,-1,-2),(-1,0,-0.75)]#,(1,-1,-1),(0,0,-0.3)]

## The restriction on the gradient of edges of the steiner tree
gradient_restriction=1/7.0

## 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,args={},args2={'color':'red'}): 
    e = g.edges() 
    #print e 
    g = Graphics() 
    points = [] 
    for i,j,w in e: 
        g += line3d([pos[i],pos[j]],**args) 
        if i not in points: 
            points.append(i) 
            if i in special: 
                g += point3d(pos[i],**args2) 
            else: 
                g += point3d(pos[i]) 
        if j not in points: 
            points.append(j) 
            if j in special: 
                g += point3d(pos[j],**args2) 
            else: g += point3d(pos[j]) 
    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 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 3 vertices 
# returning True, p iff p is the steiner point and is 
# different to p1 p2 p3
def find_steiner_point_3(m,p1,p2,p3):
    if not ( len(p1) == len(p2) and len(p1) == len(p3) ):
        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/3,add_points(p1,p2,p3))
    centroid_dist = total_distance(m,centroid,p1,p2,p3)
    
    
    # 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,p1,p2,p3)
            dist2 = total_distance(m,guess2,p1,p2,p3)
            
            # 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 
    if metric(m,guess,p1) < threshold:
        return False,p1
    elif metric(m,guess,p2) < threshold:
        return False,p2
    elif metric(m,guess,p3) < threshold:
        return False,p3
    else:
        return True,guess



def get_all_3_steiner_points(ps):
    return [find_steiner_point_3(m, *p) for p in k_subsets(ps,3)]

## In this cell we compute the candidates for the steiner points
## and then only keep the useful/valid ones

m = gradient_restriction

# make sure all the nodes are floating point (not symbolic) for efficiency
fnodes = [ [ float(x) for x in p ] for p in nodes ] 

# Get a few candidates for steiner points
new_nodes = [p for is_new_point, p in get_all_3_steiner_points(fnodes) if is_new_point ]

# compute the (normal) convex hull of the points (the field=RDF
#   is to stop the nodes being rounded to integral points)
convex_hull = Polyhedron(fnodes,field=RDF)

# filter out the steiner points outside this hull
valid_nodes = []
invalid_nodes = []
for n in new_nodes:
    if convex_hull.contains(n):
        if n not in valid_nodes:
            valid_nodes.append(n)
    else:
        invalid_nodes.append(n)
        
print len(valid_nodes),"of the",len(new_nodes),"new nodes are valid"

# Make a plot of the interesting bits
G = Graphics()
if len(invalid_nodes) > 0:
    G += points(invalid_nodes,color='orange',size=3)
if len(valid_nodes) > 0:
    G += points(valid_nodes)
if len(fnodes) > 0:
    G += points(fnodes,color='red',size=10)
    G += convex_hull.show(opacity=0.1)
    
print 
print "====================================="
print "  The candidates for steiner points"
G.show(aspect_ratio=1)

2 of the 4 new nodes are valid ===================================== The candidates for steiner points

## Now we construct the complete graph on these vertex, with the appropriate weights

all_nodes = fnodes+valid_nodes
n_o_nodes = len(all_nodes)

# a dictionary of the node numbers => euclidean position
node_dict = {}
for i,node in enumerate(all_nodes):
    node_dict[i] = node

# the dictionary required to make the graph, i.e:
#    each node number => (dictionary of neighbour number => weight of edge)
graph_dict = {}
for start in range(n_o_nodes):
    start_point = node_dict[start]
    ends_dict = {}
    for end in range(start+1,n_o_nodes):
        ends_dict[end] = metric(m,start_point,node_dict[end])
    graph_dict[start] = ends_dict

complete_graph=Graph(graph_dict,weighted=True)
graph_plot(complete_graph,node_dict,range(len(nodes)))

## compute the steiner tree

smt_of_complete = complete_graph.steiner_tree(range(len(nodes)),weighted=True)

## print it out
graph_plot(smt_of_complete,node_dict,range(len(nodes)))