# Ticket #9923: trac_9923.patch

File trac_9923.patch, 10.8 KB (added by ncohen, 10 years ago)
• ## sage/graphs/digraph.py

# HG changeset patch
# User Nathann Cohen <nathann.cohen@gmail.com>
# Date 1285887698 -7200
# Node ID 5c9ed0fbc0efa3956f509f6cfccc6039eeb8ac1a
# Parent  f30bb70155c938e77cdea1da6c71e688fc3502b5
trac 9923 -- Minimum Feedback Arc/Edge set through constraint generation

diff -r f30bb70155c9 -r 5c9ed0fbc0ef sage/graphs/digraph.py
 a return sorted(self.out_degree_iterator(), reverse=True) def feedback_edge_set(self, value_only=False, solver=None, verbose=0): def feedback_edge_set(self, constraint_generation= True, value_only=False, solver=None, verbose=0): r""" Computes the minimum feedback edge set of a digraph (also called feedback arc set). - When set to False, the Set of edges of a minimal edge set is returned. - constraint_generation (boolean) -- whether to use constraint generation when solving the Mixed Integer Linear Program (default: True). - solver -- (default: None) Specify a Linear Program (LP) solver to be used. If set to None, the default one is used. For more information on LP solvers and which default solver is used, see - verbose -- integer (default: 0). Sets the level of verbosity. Set to 0 by default, which means quiet. This problem is solved using Linear Programming, which certainly is not the best way and will have to be updated. The program to be solved is the following: ALGORITHM: This problem is solved using Linear Programming, in two different ways. The first one is to solve the following formulation: .. MATH:: The number of edges removed is then minimized, which is the objective. (Constraint Generation) If the parameter constraint_generation is enabled, a more efficient formulation is used : .. MATH:: \mbox{Minimize : }&\sum_{(u,v)\in G} b_{(u,v)}\\ \mbox{Such that : }&\\ &\forall C\text{ circuits }\subseteq\in G, \sum_{uv\in C}b_{(u,v)}\geq 1\\ As the number of circuits contained in a graph is exponential, this LP is solved through constraint generation. This means that the solver is sequentially asked to solved the problem, knowing only a portion of the circuits contained in G, each time adding to the list of its constraints the circuit which its last answer had left intact. EXAMPLES: If the digraph is created from a graph, and hence is symmetric sage: (u,v,l) = g.edge_iterator().next() sage: (u,v) in feedback or (v,u) in feedback True TESTS: Comparing with/without constraint generation:: sage: g = digraphs.RandomDirectedGNP(10,.3) sage: x = g.feedback_edge_set(value_only = True) sage: y = g.feedback_edge_set(value_only = True, constraint_generation = False) sage: x == y True """ if constraint_generation: from sage.graphs.generic_graph_pyx import minimum_feedback_arc_set obj,mfas = minimum_feedback_arc_set(self, constraint_log = verbose, solver_log = verbose) if value_only: return obj else: return mfas from sage.numerical.mip import MixedIntegerLinearProgram, Sum
• ## sage/graphs/generic_graph_pyx.pyx

diff -r f30bb70155c9 -r 5c9ed0fbc0ef sage/graphs/generic_graph_pyx.pyx
 a sage_free(u) sage_free(v) ############################## # Minimum feedback arc set   # ############################## from sage.numerical.backends.generic_backend cimport GenericBackend cpdef minimum_feedback_arc_set(g, constraint_log = 0,solver_log = 0): r""" Returns a minimum feedback arc set of the given digraph. INPUT: - g -- a digraph - constraint_log (boolean) -- whether to prin informations about the constraints generation. 0 by default. - solver_log (integer) -- log lever for GLPK. 0 by default. ALGORITHM: Mixed Integer Linear Program with constraint generation. Each edge has a corresponding binary variable, set to 1 if the edge is included in the MFAS, and set to 0 otherwise. We want to minimize their number. We first solve the problem without any constraint, which means all the variables have a value of 0, and that no edges are to be removed. We are likely to find a circuit in this graph, so we look for it using the method find_cycle, then add to our LP that the sum of the variables must be at least 1 over this cycle (it has to be broken !) EXAMPLES: Disjoint union of cycles (this method is tested through DiGraph.feedback_edge_set too):: sage: from sage.graphs.generic_graph_pyx import minimum_feedback_arc_set sage: g = digraphs.Circuit(5)*10 sage: minimum_feedback_arc_set(g)[0] 10.0 """ from sage.numerical.backends.generic_backend import get_solver cdef int n = g.order() cdef int m = g.size() cdef list indices cdef list values # Associating its variable number to each edge of the graph cdef dict edge_id = dict( [(x,y) for (y,x) in enumerate(g.edges(labels = False))] ) cdef GenericBackend p p = get_solver(constraint_generation = True) p.add_variables(m) p.set_sense(-1) # Variables are binary, and their coefficient in the objective is # 1 for 0 <= i < m: p.set_variable_type(i,0) p.set_objective([1]*m) p.set_verbosity(solver_log) p.solve() # For as long as we do not break because the digraph is # acyclic.... while (1): # Find a cycle c = find_cycle(p, g, edge_id) # If there is none, we are done ! if not c: break if constraint_log: print "Adding a constraint on cycle : ",c indices = [] for 0<= i < len(c)-1: indices.append(edge_id[(c[i],c[i+1])]) indices.append(edge_id[(c[len(c)-1],c[0])]) values = [1] * len(c) p.add_linear_constraint(indices, values, -1, 1) p.solve() # Getting the objective back obj = p.get_objective_value() # listing the edges contained in the MFAS edges = [(u,v) for u,v in g.edges(labels = False) if p.get_variable_value(edge_id[(u,v)]) > .5] # Accomplished ! return obj,edges # Find a cycle in the digraph given by the previous solution returned cdef list find_cycle(GenericBackend lp, g, dict edges_id): r""" Returns a cycle if there exists one in the digraph where the edges removed by the current LP solution have been removed. If there is none, an empty list is returned. """ # Vertices that have already been seen cdef dict seen = dict( [(v,False) for v in g] ) # Activated vertices cdef dict activated = dict( [(v,True) for v in g] ) # Vertices whose neighbors have already been added to the stack cdef dict tried = dict( [(v,False) for v in g] ) # Parent of a vertex in the discovery tree cdef dict parent = {} # The vertices left to be visited cdef list stack = [] # We try any vertex as the source of the exploration tree for v in g: # We are not interested in trying de-activated vertices if not activated[v]: continue stack = [v] seen[v] = True # For as long as some vertices are to be visited while stack: # We take the last one (depth-first search) u = stack.pop(-1) # If we never tried it, we put it at the end of the stack # again, but remember we tried it once. if not tried[u]: tried[u] = True seen[u] = True stack.append(u) # If we tried all the path starting from this vertex # already, it means it leads to no circuit. We can remove # it, and go to the next one else: seen[u] = False activated[u] = False continue # For each out-neighbor of the current vertex for uu in g.neighbors_out(u): # We are not interested in using the edges removed by # the previous solution if lp.get_variable_value(edges_id[(u,uu)]) > .5: continue # If we have found a new vertex, we put it at the end # of the stack. We ignored de-activated vertices. if not seen[uu]: if activated[uu]: parent[uu] = u stack.append(uu) # If we have already met this vertex, it means we have # found a circuit ! else: # We build it, then return it answer = [u] tmp = u while u != uu: u = parent.get(u,uu) answer.append(u) answer.reverse() return answer # Each time we have explored all the descendants of a vertex v # and met no circuit, we de-activate all these vertices, as we # know we do not need them anyway. for v in g: if seen[v]: seen[v] = False activated[v] = False tried[v] = True # No Cycle... Good news ! Let's return it. return [] cpdef tuple find_hamiltonian( G, long max_iter=100000, long reset_bound=30000, long backtrack_bound=1000, find_path=False ): r""" answers are still satisfiable path (the algorithm would raise an exception otherwise):: sage: for i in range(200): ...      g = graphs.RandomGNP(20,.1) ...      _ = fh(G,find_path=True) sage: for i in range(50):              # long ...      g = graphs.RandomGNP(20,.4)   # long ...      _ = fh(G,find_path=True)      # long """ from sage.misc.prandom import randint