2746.sagews

nodes=[(0.0,0.0,0.0),(0,0,-1.0),(0,0.83,-0.5)]
m=1

# helpers
def get_sets_of_n(list,n):
    if len(list) < n or n <= 0:
        return []
    elif n == 1:
        return map(lambda x:[x],list)
    newlist = []
    for index,i in enumerate(list[:-n+1]):
        for j in get_sets_of_n(list[index+1:],n-1):
            newlist.append([i]+j)
    return newlist
    
def metric(m,a,b):
    horiz = (a[0]-b[0])^2+(a[1]-b[1])^2
    vert = abs(a[2]-b[2])
    #return '',sqrt(horiz+vert^2)
    if vert^2 > m*horiz:
        return 'b',float(sqrt(1+m^-2)*vert)
    elif vert^2/horiz == m:
        return 'm',float(sqrt(horiz + vert^2))
    else:
         return 'f',float(sqrt(horiz + vert^2))
def nub(list):
    last = None
    newlist = []
    for i in list:
        if i != last:
            newlist.append(i)
            last = i
    return newlist
def average(list):
    return sum(list)/len(list)
def trace(a):
    print a
    return a

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 scale_point(s,p):
    return map(lambda x:x*s,p)
def add_points(p,q):
    return map(lambda x,y:x+y,p,q)    
    
def compute_metric_steiner2(l):
    if len(l) != 2:
        return None
    first_guess=map(average,(zip(*l)))
    return first_guess
def compute_metric_steiner3(l,p=False):
    if len(l) != 3:
        return None
    guess=map(average,(zip(*l)))
    def compute_dist(xx):
        return metric(m,xx,l[0])[1]+metric(m,xx,l[1])[1]+metric(m,xx,l[2])[1]
    dist = compute_dist(guess)
    step = dist/2
    while step > 1e-3:
        if p:
            print "** ", guess, dist
        best = guess
        bestd = dist
        for perm in [(step,0,0),(0,step,0),(0,0,step)]:
            nperm = scale_point(-1,perm)
            guess1 = add_points(guess,perm)
            guess2 = add_points(guess,nperm)
            d1 = compute_dist(guess1)
            d2 = compute_dist(guess2)
            if p:
                print ">>> ",guess1,d1,"\n    ", guess2,d2
            if d1 < bestd:
                best = guess1
                bestd = d2
            if d2 < bestd:
                best = guess2
                bestd = d2
        guess = best
        dist = bestd
        step /= 2 
    return guess
    
def compute_metric_steiner4(l):
    if len(l) != 4:
        return None
    first_guess=map(average,(zip(*l)))
    return first_guess
def compute_metric_steiner5(l):
    if len(l) != 5:
        return None
    first_guess=map(average,(zip(*l)))
    return first_guess

# we compute the average of all the subsets of vertices of a given size
# we first get the sets, i.e. a list of lists of points, and then transpose 
# each list in the big list, to a list of lists of length 3 of values, each of these sublists
# corresponds to a coordinate (i.e. it contains list of all the x values, then y values and z values)
# and so just average these
len_nodes = len(nodes)
node_numbers = range(len_nodes)
newnodes=map(compute_metric_steiner3, get_sets_of_n(nodes,3))#+map(compute_metric_steiner4, get_sets_of_n(nodes,4))#+map(compute_metric_steiner5,get_sets_of_n(nodes,5))
allnodes = nodes+nub(sorted(newnodes))

pos = {}
for index,i in enumerate(allnodes):
    pos[index] = i

c_hull = Polyhedron(vertices=nodes,field=RDF)
useful_v=node_numbers[:]
for i,p in pos.items():
    if c_hull.contains(p) and i not in useful_v: #and list(l).count(0) < len(l) - 2: # or list(l)[:len_nodes].count(0) < len_nodes:
        useful_v.append(i)    
print "Useful: ",useful_v, "\t#:",len(useful_v)
diction = {}
for i,j in get_sets_of_n(list(enumerate(allnodes)),2):
    if i[0] not in diction:
        diction[i[0]] = {j[0]:metric(m,i[1],j[1])[1]}
    else:
        diction[i[0]][j[0]] = metric(m,i[1],j[1])[1]
g=Graph(diction)
print len(allnodes)

# compute the spanning tree using those ^ vertices
mst=g.min_spanning_tree(lambda a:a[2])
mst_dict = {}
for i,j,d in mst:
    if i not in mst_dict:
        mst_dict[i] = {j:d}
    else:
        mst_dict[i][j] = d
mstg = Graph(mst_dict,weighted=True)
mstgwam = mstg.weighted_adjacency_matrix()

#print allnodes

gwam = g.weighted_adjacency_matrix()
new_wam = []
for i in range(len(allnodes)):
    row = gwam[i]
    newrow = [0]*len(row)
    if i in useful_v:
        for j in useful_v:
            newrow[j] = row[j]
    new_wam.append(newrow)

trimmed_mstg = Graph(data=Matrix(new_wam),format='weighted_adjacency_matrix')
G=Graphics()
G+=c_hull.show(opacity=0.1)
G+=graph_plot(mstg,pos,useful_v,{'opacity':0.5},{'size':7,'color':'red'})
G+=graph_plot(trimmed_mstg,pos,node_numbers,{'thickness':1},{'size':10,'color':'red'})
print "The big vertices are the given points"
print "The shaded solid is their convex hull"
G.show(viewer='tachyon',aspect_ratio=1)

Useful: [0, 1, 2, 3] #: 4 4 The big vertices are the given points The shaded solid is their convex hull

## actually compute steiner tree

#smt=g.steiner_tree(node_numbers,weighted=True)
trimmed_smt=trimmed_mstg.steiner_tree(node_numbers,weighted=True)
G=Graphics()
#G+=graph_plot(mstg,pos,node_numbers,{'color':'black'})
G+=graph_plot(trimmed_smt,pos,node_numbers)#,{'thickness':'2'},{'size':8,'color':'red'})
#G+=graph_plot(smt,pos,node_numbers,{'thickness':'2','color':'black'})
print "Steiner Tree of",nodes, "with gradient constraint",m
G.show(viewer='tachyon')

Steiner Tree of [(0.000000000000000, 0.000000000000000, 0.000000000000000), (0, 0, -1.00000000000000), (0, 0.830000000000000, -0.500000000000000)] with gradient constraint 1