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.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


#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)
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))

#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))
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.add_constraint(x_eq[u] + x_eq[v] >= w)
#for v in isolated_vertices:
    #P_eq.add_constraint(x_eq[v] >= 0)
    #P_eq.add_constraint(x_eq[v] == y[v])
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()