# HG changeset patch
# User Nathann Cohen
# 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/sage/graphs/digraph.py Tue Sep 14 22:59:53 2010 +0200
+++ b/sage/graphs/digraph.py Fri Oct 01 01:01:38 2010 +0200
@@ -1091,7 +1091,7 @@
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,6 +1114,10 @@
- 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,9 +1129,11 @@
- ``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,6 +1157,24 @@
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,7 +1198,25 @@
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
diff -r f30bb70155c9 -r 5c9ed0fbc0ef sage/graphs/generic_graph_pyx.pyx
--- a/sage/graphs/generic_graph_pyx.pyx Tue Sep 14 22:59:53 2010 +0200
+++ b/sage/graphs/generic_graph_pyx.pyx Fri Oct 01 01:01:38 2010 +0200
@@ -949,6 +949,211 @@
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,9 +1262,9 @@
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