Shared2017-03-24-102329.sagewsOpen in CoCalc
from fractions import Fraction
from random import randint
from functools import reduce

def solve(G, P, x, obj):
P.set_objective(obj)
done=False
precision=11
while not done:
s=P.solve()
x_star = P.get_values(x)
H = G.copy()
for (u,v,w) in H.edges():
H.set_edge_label(u,v, w - x_star[u] - x_star[v])
M_H = H.matching(use_edge_labels=True)

nu_H = -1*sum(w for (u,v,w) in M_H)
if round(nu_H-x_star['epsilon'],precision) >= 0:
done=True
return P
else:
P.add_constraint(sum(x[u] + x[v] - G.edge_label(u,v) for (u,v,w) in M_H) - x['epsilon'] >=0)
return P

def maximum_weight_matching(G):
M = G.matching(use_edge_labels=True)
nu = sum(w for (u,v,w) in M)
return (M,nu)

def maximum_weight_fractional_matching(G):
P = MixedIntegerLinearProgram()
x = P.new_variable(real=True, nonnegative=True)
P.set_objective(sum(G.edge_label(e[0],e[1])*x[e] for e in G.edges(labels=False)))
for v in G:
P.add_constraint(sum(x[e] for e in G.edges_incident(v, labels=False)) <= 1)
P.solve()
return P

def print_solution(P,x):
for key, value in P.get_values(x).items():
print(str(key) + ":= " + str(value))
print(" ")

def get_y(G,P,x):
x_bar = [value for key, value in P.get_values(x).items() if key != 'epsilon']

#other_x = [7/5, 22/5, 0, 33/5, 4, 27/5, 22/5, 13/5, 27/5, 27/5, 12/5, 6]
#x_bar = [1/2*a + 1/2*b for (a,b) in zip(x_bar,other_x)]
#x_bar=other_x
#x_bar = [5, 15, 15, 5, 15, 7, 7, 7, 0, 6, 0, 0, 1/10, 0]

excess_edges = [(u,v, w-x_bar[u]-x_bar[v]) for (u,v,w) in G.edges()]
print("Excess Edges")
print(excess_edges)
excess_graph = initialize_graph(excess_edges)
y_bar = optimal_fractional_cover(excess_graph)
print("Excess cover")
print(y_bar)
y = [u+v for (u,v) in zip(x_bar,y_bar)]
return y

def print_all_solutions(G,P,x):
for v in G.vertices():
print("Vertex: " + str(v) + " maximized")
print("Number of constraints " + str(P.number_of_constraints()))
P=solve(G,P,x,x[v])
print_solution(P,x)

def initialize_graph(edge_list):
G = Graph()
G.weighted(true)
return G

def initialize_random_graph(n, p):
G = graphs.RandomGNP(n, p)
G.weighted(true)
for (u,v) in G.edges(labels=False):
w = randint(1,n)
G.set_edge_label(u,v,w)
return G

def initialize_lp(G, delta):
(M,nu) = maximum_weight_matching(G)
print("Max weight of a matching = " + str(nu))
print("Value of delta = " + str(delta))
nu_f=maximum_weight_fractional_matching(G).get_objective_value()
print("Max weight of a fractional matching = " + str(nu_f))
if(nu>delta or nu_f<=delta):
P = MixedIntegerLinearProgram(solver='ppl')
x = P.new_variable(real=True, nonnegative=False)
P.add_constraint(sum( x[v] for v in G) == delta)
for v in G:
P.add_constraint(sum(x[u]+x[v]-w for (u,v,w) in M) - x['epsilon'] >= 0)
return (solve(G,P,x,x['epsilon']),x)

def equality_graph(G,y):
edge_list=[]
for (a,b,w) in G.edges():
if(y[a] + y[b] == w):
edge_list.append((a,b,w))
G_equals = initialize_graph(edge_list)
return G_equals

def optimal_fractional_cover(G):
P = MixedIntegerLinearProgram(solver='ppl')
x = P.new_variable(real=True, nonnegative=True)
for (u,v,w) in G.edges():
P.set_objective(sum(-1*x[v] for v in G.vertices()))
P.solve()
return [value for key,value in P.get_values(x).items()]

def inessential_vertices(G):
nu = G.matching(value_only=True)
inessential_vertices = []
for v in G.vertices():
H = G.copy()
H.delete_vertex(v)
if nu == H.matching(value_only=True):
inessential_vertices.append(v)
return inessential_vertices

def excess_weights_graph(G, x):
H = G.copy()
for (u,v,w) in H.edges():
H.set_edge_label(u,v, w-x[u]-x[v])
return H
def positive_excess_weights_graph(G,x):
H = excess_weights_graph(G,x)
for (u,v,w) in H.edges():
if (w<=0):
H.delete_edge(u,v)
return H
def optimal_face(G, opt_epsilon):
F = MixedIntegerLinearProgram(solver="ppl")
b = F.new_variable(integer=True, nonnegative=True)

F.set_objective(sum(b[e] for e in G.edges(labels=False)))
for v in G:
F.add_constraint(sum(b[e] for e in G.edges_incident(v, labels=False)) <= 1)
F.add_constraint(sum(G.edge_label(e[0],e[1])*b[e] for e in G.edges(labels=False)) == opt_epsilon)
return (F,b)

def matched(M):
return [u for (u,v) in M] + [v for (u,v) in M]

def normal_objective(G,K):
L = MixedIntegerLinearProgram(solver="ppl")
l = L.new_variable(real=True, nonnegative=False)
L.set_objective(sum(l[e] for e in G.edges(labels=False)))
for e in G.edges(labels=False):
for M in K:
L.add_constraint(sum(l[e] for e in M) == 0)
L.solve()
return L.get_values(l)

def universal_allocation(G, P, x):
x_star = [value for key, value in P.get_values(x).items() if key != 'epsilon']
opt_epsilon = -1*P.get_values(x)['epsilon']
(F,b) = optimal_face(excess_weights_graph(G,x_star), opt_epsilon)
F.solve()
M = [key  for (key,value) in F.get_values(b).items() if(value == 1)]
K = []
k = normal_objective(G,K)
flag = True
while flag:
P.set_objective(sum(x[u] + x[v] for (u,v) in M))
P.solve()
y = [value for key, value in P.get_values(x).items() if key != 'epsilon']
yM = sum(y[u] + y[v] for (u,v) in M)
xM = sum(x_star[u] + x_star[v] for (u,v) in M)
if yM > xM:
for v in G.vertices():
x_star[v] = 1/2*(x_star[v] + y[v])
(F,b) = optimal_face(excess_weights_graph(G,x_star), opt_epsilon)
F.set_objective(sum(k[e]*b[e] for e in G.edges(labels=False)))
k1 = F.solve()
M1 = [key  for (key,value) in F.get_values(b).items() if(value == 1)]
if k1 > 0:
M = M1
else:
F.set_objective(-1*sum(k[e]*b[e] for e in G.edges(labels=False)))
k2=F.solve()
M2 = [key  for (key,value) in F.get_values(b).items() if(value == 1)]
if k2 > 0:
M = M2
else:
flag = False
else:
K = K + [M]
k = normal_objective(G,K)
F.set_objective(sum(k[e]*b[e] for e in G.edges(labels=False)))
k1 = F.solve()
M1 = [key  for (key,value) in F.get_values(b).items() if(value == 1)]
if k1 > 0:
M = M1
else:
F.set_objective(-1*sum(k[e]*b[e] for e in G.edges(labels=False)))
k2=F.solve()
M2 = [key  for (key,value) in F.get_values(b).items() if(value == 1)]
if k2 > 0:
M = M2
else:
flag = False
return x_star

#Specify your graph here as a list of triples (u,v,weight) specifying the edge endpoints and the edge weight
#If you are working on my sheet (instead of cloning) I just ask that you comment out my previous graph so it's preserved.
#When I have more time I'll add more features for usability, and file io so we can save and load graphs as desired.
size = 12
probability = 0.35
G = initialize_random_graph(size, probability)
while(maximum_weight_matching(G)[1] > maximum_weight_fractional_matching(G).get_objective_value()-1):
G = initialize_random_graph(size,probability)
#G = initialize_graph([(0, 3, 5), (0, 5, 8), (0, 8, 8), (0, 9, 6), (1, 3, 3), (1, 5, 8), (1, 8, 10), (1, 9, 9), (1, 10, 7), (2, 5, 3), (2, 11, 6), (3, 6, 11), (3, 8, 12), (3, 9, 8), (3, 11, 7), (4, 5, 7), (4, 9, 1), (4, 11, 10), (5, 8, 12), (5, 10, 8), (5, 11, 10), (7, 9, 8), (7, 10, 5), (7, 11, 8), (8, 10, 7)])
#G = initialize_graph([(0, 2, 8), (0, 3, 8), (1, 3, 10), (1, 4, 7), (2, 3, 12), (2, 4, 8)])
#G = initialize_graph([(0, 4, 8), (0, 1, 8), (1, 2, 9), (2, 3, 7), (4, 1, 11), (3, 4, 8)])
#G = initialize_graph([(0, 4, 8), (0, 1, 8), (1, 2, 9), (2, 3, 7), (4, 1, 11), (3, 4, 8), (5,6,7), (6,7,7), (7,8,7), (8,9,7), (9,5,7)])
#G = initialize_graph([(u,v,70*w) for (u,v,w) in [(0,1,1),(0,2,1), (1,2,1), (2,3,6/7), (3,4,1), (4,5,1), (5,6,1), (6,4,1), (7,8,1), (8,9,1), (9,10,1), (10,11,1), (11,7,1), (3,8,6/7)]])
#G = initialize_graph([(u,v,70*w) for (u,v,w) in [(0,1,1),(0,2,1), (1,2,1), (2,3,6/7), (3,4,1), (4,5,1), (5,6,1), (6,7,1), (7,8,1), (8,4,1)]])
#G = initialize_graph([(u,v,70*w) for (u,v,w) in [(0,1,1),(0,2,1), (1,2,1), (2,3,6/7), (3,4,1), (4,5,1), (5,6,1), (6,4,1), (7,8,1), (8,9,1), (9,7,1), (3,8,6/7)]])
#G = initialize_graph([(u,v,7*w) for (u,v,w) in [(0,1,1),(0,2,1), (1,2,1), (4,5,1), (5,6,1), (6,4,1), (7,8,1), (8,9,1), (9,10,1), (10,11,1), (11,7,1)]])
#G = initialize_graph([(u,v,7*w) for (u,v,w) in [(0,1,19/21),(0,2,19/21), (1,2,19/21), (2,3,6/7), (3,4,1), (4,5,1), (5,6,1), (6,4,1), (7,8,1), (8,9,1), (9,10,1), (10,11,1), (11,7,1)]])
#G = initialize_graph([(u,v,7*w) for (u,v,w) in [(0,1,1),(0,2,10), (1,2,1), (2,3,6/7), (3,8,1), (8,4,1), (4,5,1), (5,6,1), (6,4,1)]])
G = initialize_graph([(u,v,70*w) for (u,v,w) in [(0,1,1),(0,2,1), (1,2,1), (2,3,6/7), (3,4,1),  (4,5,1), (5,6,1), (6,4,1)]])
#G = initialize_graph([(0, 7, 2), (0, 8, 8), (0, 11, 12), (1, 3, 10), (1, 4, 10), (1, 5, 11), (1, 11, 9), (2, 5, 7), (2, 10, 1), (3, 4, 8), (3, 5, 10), (3, 11, 9), (4, 9, 11), (5, 6, 5), (5, 7, 6), (5, 8, 9), (5, 10, 11), (6, 9, 10), (6, 11, 5), (7, 8, 9)])
#G = initialize_graph([(0,1,20), (1,2,20), (2,3,20), (3,4,20), (4,0,20), (5,6,20), (6,7,20), (7,5,20), (8,9,1), (9,10,1), (11,12, 1/10), (12,13,1/10), (0,9,11), (1,9,11), (2,9,11), (3,9,11), (4,9,11), (5,12, 11/10), (6,12,11/10), (7,12,11/10)])
#G = initialize_graph([(0, 1, 18/5), (0, 2, 18/5), (1, 2, 18/5), (2, 3, 0), (3, 4, 0), (4, 5, 8/5), (4, 6, 8/5), (5, 6, 8/5), (7, 8, 0), (7, 11, 0), (8, 9, 0), (9, 10, 0), (10, 11, 0)])
G.plot(edge_labels='true')
#G.add_edges([ (1,2,1), (2,3,1), (3,4,1), (4,5,1), (5,1,1), (2,4,2)])
#Run this block of code to compile and plot your graph

()
#graph = initialize_graph([(0,1,1), (1,2,1), (2,3,1), (3,4,1), (4,0,1), (5,6,1), (6,7,1), (0,6,1), (1,6,1), (2,6,1), (3,6,1), (4,6,1), (8,9,1), (9,10, 1), (10,8,1), (11, 12,1), (12,13,1), (8,12,1), (9,12,1), (10,12,1)])
#view(graph)
#opts = graph.latex_options()
#opts.set_option('edge_labels', True)
#opts.set_option('graphic_size', (15,7))
#print(opts.latex())

print(G.edges())
(P,x) = initialize_lp(G, maximum_weight_matching(G)[1])
print("Initial  Solution")
print_solution(P,x)
print_all_solutions(G,P,x)

x_star = universal_allocation(G,P,x)
print(x_star)
H = positive_excess_weights_graph(G,x_star)
for v in H.vertices():
if H.neighbors(v) == []:
H.delete_vertex(v)
H.plot(edge_labels=True)
print(H.edges(labels=True))
inessential_H = inessential_vertices(H)
factor_critical = True
if sorted(inessential_H) != sorted(H.vertices()):
factor_critical = False
excess_balanced = True
for cc in H.connected_components():
C = H.subgraph(cc)
(u,v,w0) = C.edges(labels=True)[0]
for (u,v,w) in C.edges(labels=True):
if w != w0:
excess_balanced = False
break
print("Factor Critical: " + str(factor_critical))

print("Excess Balanced: " + str(excess_balanced))

#Run this block of code for leastcore calculations

[(0, 1, 70), (0, 2, 70), (1, 2, 70), (2, 3, 60), (3, 4, 70), (4, 5, 70), (4, 6, 70), (5, 6, 70)] Max weight of a matching = 210 Value of delta = 210 Max weight of a fractional matching = 245.0 Initial Solution 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 0 maximized Number of constraints 16 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 1 maximized Number of constraints 16 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 2 maximized Number of constraints 16 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 3 maximized Number of constraints 16 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 4 maximized Number of constraints 17 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 5 maximized Number of constraints 17 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 Vertex: 6 maximized Number of constraints 17 0:= 24 1:= 24 2:= 24 3:= 36 4:= 34 5:= 34 6:= 34 epsilon:= -24 [24, 24, 24, 36, 34, 34, 34]
[(0, 1, 22), (0, 2, 22), (1, 2, 22), (4, 5, 2), (4, 6, 2), (5, 6, 2)] Factor Critical: True Excess Balanced: True
 #y = optimal_fractional_cover(G)
y = get_y(G,P,x)
#y=[3.5,3.5,66.5, 3.5, 3.5, 3.5, 3.5, 3.5,3.5]
#y=[3.5,3.5,3.5, 43/10, 3.5, 3.5, 3.5, 3.5,3.5,3.5,3.5,3.5]
print(y)
G_eq = equality_graph(G,y)
print(G_eq.matching())
G_eq.plot(edge_labels='true')
D = inessential_vertices(G_eq)
print(D)
G_eq_D = G_eq.subgraph(D)
G_eq_D.plot(edge_labels='true')
#TODO Figure out how to obtain non-chords:
chords = []
for (u,v,w) in G.edges():
for cc in G_eq_D.connected_components():
if u in cc and v in cc:
chords.append((u,v,w))
non_chords = list(set(G.edges()) - set(chords) - set(G_eq.edges()))
print(non_chords)
isolated_vertices = list(set(G.vertices()) - set(G_eq.vertices()))
print(isolated_vertices)
(P_eq,x_eq) = initialize_lp(G_eq, maximum_weight_matching(G)[1]-1/2*(7/2 + 43/10))#sum(y[v] for v in isolated_vertices))
slack_edges = list(set(G.edges())-set(G_eq.edges()))
print(slack_edges)
for (u,v,w) in non_chords:
print("Not chord " + str((u,v,w)))
#for v in isolated_vertices:
P_eq=solve(G_eq,P_eq,x_eq, x_eq['epsilon'])
print_all_solutions(G_eq,P_eq,x_eq)

Q = MixedIntegerLinearProgram(solver='ppl')
x_Q = Q.new_variable(real=True, nonnegative=False)
delta=sum(w for (u,v,w) in G.matching(use_edge_labels=True))
Q.add_constraint(sum(x_Q[v] for v in G) == delta)
for v in G:
Q.set_objective(x_Q['epsilon'])
for v in G_eq_D.vertices():
for cc in G_eq_D.connected_components():
for u in cc:
for v in cc:
Q.add_constraint(x_Q[v] -y[v] == x_Q[u] - y[u])
for (u,v,w) in list(set(G.edges()) - set(G_eq_D.edges())):
M_Q = G_eq_D.matching(use_edge_labels='True')
Q.add_constraint(sum(x_Q[u] + x_Q[v] -w for (u,v,w) in M_Q) >= x_Q['epsilon'])
epsilon_Q = Q.solve()
print_solution(Q,x_Q)
for u in G:
print(str(u) + " maximized")
Q.set_objective(x_Q[u])
Q.solve()
print_solution(Q,x_Q)

G = initialize_graph([(0,1,11/5), (1,2,11/5), (2,0, 11/5), (2,3,0), (3,4,0), (4,5,1/5), (5,6,1/5), (6,4,1/5)])
G.plot(edge_labels=True)
y = optimal_fractional_cover(G)
print(y)

P.constraints()