| Hosted by CoCalc | Download
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#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()