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']
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.add_edges(edge_list)
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):
print("Bad 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(x[v] >= 0)
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.add_constraint(x[u] + x[v] >= w)
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):
L.add_constraint(l[e] <= 1)
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
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([(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.plot(edge_labels='true')
()
print(G.edges())
(P,x) = initialize_lp(G, maximum_weight_matching(G)[1])
print("Initial Solution")
print_solution(P,x)
P.add_constraint(x['epsilon'] == P.get_values(x)['epsilon'])
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))
y = get_y(G,P,x)
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')
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))
P_eq.add_constraint(x_eq['epsilon'] == P_eq.get_values(x_eq)['epsilon'])
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)))
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.add_constraint(x_Q[v] >= 0)
Q.set_objective(x_Q['epsilon'])
for v in G_eq_D.vertices():
Q.add_constraint(x_Q[v] <= y[v])
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())):
Q.add_constraint(x_Q[u] +x_Q[v] >= w)
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)
Q.add_constraint(x_Q['epsilon'] == epsilon_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()