Ticket #9923: trac_9923.patch

sage/graphs/digraph.py

@@ -1091 +1091 @@
         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).
@@ -1114 +1114 @@
           - 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
@@ -1125 +1129 @@
         - ``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::
 
@@ -1151 +1157 @@
         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 
@@ -1174 +1198 @@
            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

@@ -949 +949 @@
         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"""
@@ -1057 +1262 @@
     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