Ticket #14434: trac_14434-ALL.patch

File trac_14434-ALL.patch, 22.8 KB (added by ncohen, 6 years ago)
  • sage/graphs/digraph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1365581838 -7200
    # Node ID 0dce8ee365602d8f7af7653d95b6d765ee7117c3
    # Parent  d0b9d29d21159047139b49edc8beaf6641361e3d
    Move digraph.feedback_vertex_set to generic_digraph.py
    * * *
    implements GenericGraph.feedback_vertex_set for undirected graphs
    * * *
    Implement feedback_vertex_set for graphs -- broken doctest
    * * *
    trac 14334: clean documentation for feedback_vertex_set
    * * *
    Implement feedback_vertex_set for graphs -- a bug in is_tree
    
    diff --git a/sage/graphs/digraph.py b/sage/graphs/digraph.py
    a b  
    8585    :delim: |
    8686
    8787    :meth:`~DiGraph.feedback_edge_set` | Computes the minimum feedback edge (arc) set of a digraph
    88     :meth:`~DiGraph.feedback_vertex_set` | Computes the minimum feedback vertex set of a digraph.
    89 
    9088
    9189Methods
    9290-------
     
    16041602
    16051603                return [(u,v) for (u,v) in self.edges(labels=None) if b_sol[(u,v)]==1]
    16061604
    1607     def feedback_vertex_set(self, value_only=False, solver=None, verbose=0, constraint_generation = True):
    1608         r"""
    1609         Computes the minimum feedback vertex set of a digraph.
    1610 
    1611         The minimum feedback vertex set of a digraph is a set of vertices
    1612         that intersect all the circuits of the digraph.
    1613         Equivalently, a minimum feedback vertex set of a DiGraph is a set
    1614         `S` of vertices such that the digraph `G-S` is acyclic. For more
    1615         information, see the
    1616         `Wikipedia article on feedback vertex sets
    1617         <http://en.wikipedia.org/wiki/Feedback_vertex_set>`_.
    1618 
    1619         INPUT:
    1620 
    1621         - ``value_only`` -- boolean (default: ``False``)
    1622 
    1623           - When set to ``True``, only the minimum cardinal of a minimum vertex
    1624             set is returned.
    1625 
    1626           - When set to ``False``, the ``Set`` of vertices of a minimal feedback
    1627             vertex set is returned.
    1628 
    1629         - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
    1630           solver to be used. If set to ``None``, the default one is used. For
    1631           more information on LP solvers and which default solver is used,
    1632           see the method
    1633           :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
    1634           of the class
    1635           :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.
    1636 
    1637         - ``verbose`` -- integer (default: ``0``). Sets the level of
    1638           verbosity. Set to 0 by default, which means quiet.
    1639 
    1640         - ``constraint_generation`` (boolean) -- whether to use constraint
    1641           generation when solving the Mixed Integer Linear Program (default:
    1642           ``True``).
    1643 
    1644         ALGORITHM:
    1645 
    1646         This problem is solved using Linear Programming, which certainly is not
    1647         the best way and will have to be replaced by a better algorithm.  The
    1648         program to be solved is the following:
    1649 
    1650         .. MATH::
    1651 
    1652             \mbox{Minimize : }&\sum_{v\in G} b_v\\
    1653             \mbox{Such that : }&\\
    1654             &\forall (u,v)\in G, d_u-d_v+nb_u+nb_v\geq 0\\
    1655             &\forall u\in G, 0\leq d_u\leq |G|\\
    1656                        
    1657         A brief explanation:
    1658 
    1659         An acyclic digraph can be seen as a poset, and every poset has a linear
    1660         extension. This means that in any acyclic digraph the vertices can be
    1661         ordered with a total order `<` in such a way that if `(u,v)\in G`, then
    1662         `u<v`.  Thus, this linear program is built in order to assign to each
    1663         vertex `v` a number `d_v\in [0,\dots,n-1]` such that if there exists an
    1664         edge `(u,v)\in G` then either `d_v<d_u` or one of `u` or `v` is removed.
    1665         The number of vertices removed is then minimized, which is the
    1666         objective.
    1667 
    1668         (Constraint Generation)
    1669  
    1670         If the parameter ``constraint_generation`` is enabled, a more efficient
    1671         formulation is used :
    1672  
    1673         .. MATH::
    1674  
    1675             \mbox{Minimize : }&\sum_{v\in G} b_{v}\\ 
    1676             \mbox{Such that : }&\\ 
    1677             &\forall C\text{ circuits }\subseteq G, \sum_{v\in C}b_{v}\geq 1\\
    1678  
    1679         As the number of circuits contained in a graph is exponential, this LP
    1680         is solved through constraint generation. This means that the solver is
    1681         sequentially asked to solved the problem, knowing only a portion of the
    1682         circuits contained in `G`, each time adding to the list of its
    1683         constraints the circuit which its last answer had left intact.
    1684 
    1685         EXAMPLES:
    1686 
    1687         In a digraph built from a graph, any edge is replaced by arcs going in
    1688         the two opposite directions, thus creating a cycle of length two.
    1689         Hence, to remove all the cycles from the graph, each edge must see one
    1690         of its neighbors removed : a feedback vertex set is in this situation a
    1691         vertex cover::
    1692 
    1693             sage: cycle=graphs.CycleGraph(5)
    1694             sage: dcycle=DiGraph(cycle)
    1695             sage: cycle.vertex_cover(value_only=True)
    1696             3
    1697             sage: feedback = dcycle.feedback_vertex_set()
    1698             sage: len(feedback)
    1699             3
    1700             sage: (u,v,l) = cycle.edge_iterator().next()
    1701             sage: u in feedback or v in feedback
    1702             True
    1703            
    1704         For a circuit, the minimum feedback arc set is clearly `1`::
    1705 
    1706             sage: circuit = digraphs.Circuit(5)
    1707             sage: circuit.feedback_vertex_set(value_only=True) == 1
    1708             True
    1709 
    1710         TESTS:
    1711  
    1712         Comparing with/without constraint generation::
    1713  
    1714             sage: g = digraphs.RandomDirectedGNP(10,.3)
    1715             sage: x = g.feedback_vertex_set(value_only = True)
    1716             sage: y = g.feedback_vertex_set(value_only = True,
    1717             ...            constraint_generation = False)
    1718             sage: x == y
    1719             True
    1720          """
    1721 
    1722         # It would be a pity to start a LP if the digraph is already acyclic
    1723         if self.is_directed_acyclic():
    1724             if value_only:
    1725                 return 0
    1726             return []
    1727 
    1728         from sage.numerical.mip import MixedIntegerLinearProgram
    1729 
    1730         ########################################
    1731         # Constraint Generation Implementation #
    1732         ########################################
    1733         if constraint_generation:
    1734 
    1735             p = MixedIntegerLinearProgram(constraint_generation = True,
    1736                                           maximization = False)
    1737 
    1738             # An variable for each vertex
    1739             b = p.new_variable(binary = True)
    1740 
    1741             # Variables are binary, and their coefficient in the objective is 1
    1742 
    1743             p.set_objective( p.sum( b[v] for v in self))
    1744 
    1745             p.solve(log = verbose)
    1746 
    1747             # For as long as we do not break because the digraph is
    1748             # acyclic....
    1749             while (1):
    1750 
    1751                 # Building the graph without the edges removed by the LP
    1752                 h = self.subgraph(vertices =
    1753                                   [v for v in self if p.get_values(b[v]) < .5])
    1754 
    1755                 # Is the digraph acyclic ?
    1756                 isok, certificate = h.is_directed_acyclic(certificate = True)
    1757 
    1758                 # If so, we are done !
    1759                 if isok:
    1760                     break
    1761 
    1762                 if verbose:
    1763                     print "Adding a constraint on circuit : ",certificate
    1764 
    1765                 # There is a circuit left. Let's add the corresponding
    1766                 # constraint !
    1767 
    1768                 p.add_constraint( p.sum( b[v] for v in certificate), min = 1)
    1769                
    1770                 obj = p.solve(log = verbose)
    1771 
    1772             if value_only:
    1773                 return obj
    1774            
    1775             else:
    1776            
    1777                 # listing the edges contained in the MFAS
    1778                 return [v for v in self if p.get_values(b[v]) > .5]
    1779 
    1780 
    1781         else:
    1782 
    1783         ######################################
    1784         # Ordering-based MILP Implementation #
    1785         ######################################
    1786 
    1787             p = MixedIntegerLinearProgram(maximization=False, solver=solver)
    1788            
    1789             b = p.new_variable(binary = True)
    1790             d = p.new_variable(integer = True)
    1791             n = self.order()
    1792 
    1793             # The removed vertices cover all the back arcs ( third condition )
    1794             for (u,v) in self.edges(labels=None):
    1795                 p.add_constraint(d[u]-d[v]+n*(b[u]+b[v]),min=1)
    1796 
    1797             for u in self:
    1798                 p.add_constraint(d[u],max=n)
    1799 
    1800             p.set_objective(p.sum([b[v] for v in self]))
    1801 
    1802             if value_only:
    1803                 return Integer(round(p.solve(objective_only=True, log=verbose)))
    1804             else:
    1805                 p.solve(log=verbose)
    1806                 b_sol=p.get_values(b)
    1807 
    1808                 return [v for v in self if b_sol[v]==1]
    1809 
    1810 
    18111605    ### Construction
    18121606
    18131607    def reverse(self):
  • sage/graphs/generic_graph.py

    diff --git a/sage/graphs/generic_graph.py b/sage/graphs/generic_graph.py
    a b  
    281281
    282282    :meth:`~GenericGraph.steiner_tree` | Returns a tree of minimum weight connecting the given set of vertices.
    283283    :meth:`~GenericGraph.edge_disjoint_spanning_trees` | Returns the desired number of edge-disjoint spanning trees/arborescences.
     284    :meth:`~GenericGraph.feedback_vertex_set` | Computes the minimum feedback vertex set of a (di)graph.
    284285    :meth:`~GenericGraph.multiway_cut` | Returns a minimum edge multiway cut
    285286    :meth:`~GenericGraph.max_cut` | Returns a maximum edge cut of the graph.
    286287    :meth:`~GenericGraph.longest_path` | Returns a longest path of ``self``.
     
    49984999        height = p.new_variable(dim = 2)
    49995000
    50005001        # cut[e] represents whether e is in the cut
    5001         cut = p.new_variable()
     5002        cut = p.new_variable(binary = True)
    50025003
    50035004        # Reorder
    50045005        R = lambda x,y : (x,y) if x<y else (y,x)
     
    50405041                    p.add_constraint( height[(s,t)][u] - height[(s,t)][v] - cut[R(u,v)], max = 0)
    50415042                    p.add_constraint( height[(s,t)][v] - height[(s,t)][u] - cut[R(u,v)], max = 0)
    50425043
    5043 
    5044         p.set_binary(cut)       
    50455044        if value_only:
    50465045            if use_edge_labels:
    50475046                return p.solve(objective_only = True, log = verbose)
     
    50535052        cut = p.get_values(cut)
    50545053
    50555054        if self.is_directed():
    5056             return filter(lambda (u,v,l) : cut[u,v] > .5, self.edge_iterator())
    5057 
    5058         else:
    5059             return filter(lambda (u,v,l) : cut[R(u,v)] > .5, self.edge_iterator())
     5055            return filter(lambda (u,v,l) : cut[u,v] == 1, self.edge_iterator())
     5056
     5057        else:
     5058            return filter(lambda (u,v,l) : cut[R(u,v)] ==1, self.edge_iterator())
    50605059
    50615060
    50625061    def max_cut(self, value_only=True, use_edge_labels=False, vertices=False, solver=None, verbose=0):
     
    55705569        vertex_used = p.get_values(vertex_used)
    55715570        if self._directed:
    55725571            g = self.subgraph(
    5573                 vertices=(v for v in self if vertex_used[v] >= 0.5),
     5572                vertices=(v for v in self if vertex_used[v] == 1),
    55745573                edges=((u,v,l) for u, v, l in self.edges()
    5575                        if edge_used[(u,v)] >= 0.5))
     5574                       if edge_used[(u,v)] == 1))
    55765575        else:
    55775576            g = self.subgraph(
    5578                 vertices=(v for v in self if vertex_used[v] >= 0.5),
     5577                vertices=(v for v in self if vertex_used[v] == 1),
    55795578                edges=((u,v,l) for u, v, l in self.edges()
    5580                        if f_edge_used(u,v) >= 0.5))
     5579                       if f_edge_used(u,v) == 1))
    55815580        if use_edge_labels:
    55825581            return sum(map(weight, g.edge_labels())), g
    55835582        else:
     
    58565855                    # We build the DiGraph representing the current solution
    58575856                    h = DiGraph()
    58585857                    for u,v,l in g.edges():
    5859                         if p.get_values(b[u][v]) > .5:
     5858                        if p.get_values(b[u][v]) == 1:
    58605859                            h.add_edge(u,v,l)
    58615860
    58625861                    # If there is only one circuit, we are done !
     
    59075906                    # We build the DiGraph representing the current solution
    59085907                    h = Graph()
    59095908                    for u,v,l in g.edges():
    5910                         if p.get_values(B(u,v)) > .5:
     5909                        if p.get_values(B(u,v)) == 1:
    59115910                            h.add_edge(u,v,l)
    59125911
    59135912                    # If there is only one circuit, we are done !
     
    61316130        else:
    61326131            raise ValueError("``algorithm`` (%s) should be 'tsp' or 'backtrack'."%(algorithm))
    61336132
     6133    def feedback_vertex_set(self, value_only=False, solver=None, verbose=0, constraint_generation = True):
     6134        r"""
     6135        Computes the minimum feedback vertex set of a (di)graph.
     6136
     6137        The minimum feedback vertex set of a (di)graph is a set of vertices that
     6138        intersect all of its cycles.  Equivalently, a minimum feedback vertex
     6139        set of a (di)graph is a set `S` of vertices such that the digraph `G-S`
     6140        is acyclic. For more information, see :wikipedia:`Feedback_vertex_set`.
     6141
     6142        INPUT:
     6143
     6144        - ``value_only`` -- boolean (default: ``False``)
     6145
     6146          - When set to ``True``, only the minimum cardinal of a minimum vertex
     6147            set is returned.
     6148
     6149          - When set to ``False``, the ``Set`` of vertices of a minimal feedback
     6150            vertex set is returned.
     6151
     6152        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
     6153          solver to be used. If set to ``None``, the default one is used. For
     6154          more information on LP solvers and which default solver is used,
     6155          see the method
     6156          :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
     6157          of the class
     6158          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.
     6159
     6160        - ``verbose`` -- integer (default: ``0``). Sets the level of
     6161          verbosity. Set to 0 by default, which means quiet.
     6162
     6163        - ``constraint_generation`` (boolean) -- whether to use constraint
     6164          generation when solving the Mixed Integer Linear Program (default:
     6165          ``True``).
     6166
     6167        ALGORITHMS:
     6168
     6169        (Constraints generation)
     6170
     6171        When the parameter ``constraint_generation`` is enabled (default) the
     6172        following MILP formulation is used to solve the problem:
     6173
     6174        .. MATH::
     6175
     6176            \mbox{Minimize : }&\sum_{v\in G} b_{v}\\
     6177            \mbox{Such that : }&\\
     6178            &\forall C\text{ circuits }\subseteq G, \sum_{v\in C}b_{v}\geq 1\\
     6179
     6180        As the number of circuits contained in a graph is exponential, this LP
     6181        is solved through constraint generation. This means that the solver is
     6182        sequentially asked to solve the problem, knowing only a portion of the
     6183        circuits contained in `G`, each time adding to the list of its
     6184        constraints the circuit which its last answer had left intact.
     6185
     6186        (Another formulation based on an ordering of the vertices)
     6187
     6188        When the graph is directed, a second (and very slow) formulation is
     6189        available, which should only be used to check the result of the first
     6190        implementation in case of doubt.
     6191
     6192        .. MATH::
     6193
     6194            \mbox{Minimize : }&\sum_{v\in G} b_v\\
     6195            \mbox{Such that : }&\\
     6196            &\forall (u,v)\in G, d_u-d_v+nb_u+nb_v\geq 0\\
     6197            &\forall u\in G, 0\leq d_u\leq |G|\\
     6198
     6199        A brief explanation:
     6200
     6201        An acyclic digraph can be seen as a poset, and every poset has a linear
     6202        extension. This means that in any acyclic digraph the vertices can be
     6203        ordered with a total order `<` in such a way that if `(u,v)\in G`, then
     6204        `u<v`.  Thus, this linear program is built in order to assign to each
     6205        vertex `v` a number `d_v\in [0,\dots,n-1]` such that if there exists an
     6206        edge `(u,v)\in G` then either `d_v<d_u` or one of `u` or `v` is removed.
     6207        The number of vertices removed is then minimized, which is the
     6208        objective.
     6209
     6210        EXAMPLES:
     6211
     6212        The necessary example::
     6213
     6214            sage: g = graphs.PetersenGraph()
     6215            sage: fvs = g.feedback_vertex_set()
     6216            sage: len(fvs)
     6217            3
     6218            sage: g.delete_vertices(fvs)
     6219            sage: g.is_forest()
     6220            True
     6221
     6222        In a digraph built from a graph, any edge is replaced by arcs going in
     6223        the two opposite directions, thus creating a cycle of length two.
     6224        Hence, to remove all the cycles from the graph, each edge must see one
     6225        of its neighbors removed: a feedback vertex set is in this situation a
     6226        vertex cover::
     6227
     6228            sage: cycle = graphs.CycleGraph(5)
     6229            sage: dcycle = DiGraph(cycle)
     6230            sage: cycle.vertex_cover(value_only=True)
     6231            3
     6232            sage: feedback = dcycle.feedback_vertex_set()
     6233            sage: len(feedback)
     6234            3
     6235            sage: (u,v,l) = cycle.edge_iterator().next()
     6236            sage: u in feedback or v in feedback
     6237            True
     6238
     6239        For a circuit, the minimum feedback arc set is clearly `1`::
     6240
     6241            sage: circuit = digraphs.Circuit(5)
     6242            sage: circuit.feedback_vertex_set(value_only=True) == 1
     6243            True
     6244
     6245        TESTS:
     6246
     6247        Comparing with/without constraint generation::
     6248
     6249            sage: g = digraphs.RandomDirectedGNP(10,.3)
     6250            sage: x = g.feedback_vertex_set(value_only = True)
     6251            sage: y = g.feedback_vertex_set(value_only = True,
     6252            ....:          constraint_generation = False)
     6253            sage: x == y
     6254            True
     6255
     6256        Bad algorithm::
     6257
     6258            sage: g = graphs.PetersenGraph()
     6259            sage: g.feedback_vertex_set(constraint_generation = False)
     6260            Traceback (most recent call last):
     6261            ...
     6262            ValueError: The only implementation available for undirected graphs is with constraint_generation set to True.
     6263        """
     6264        if not constraint_generation and not self.is_directed():
     6265            raise ValueError("The only implementation available for "
     6266                             "undirected graphs is with constraint_generation "
     6267                             "set to True.")
     6268
     6269        # It would be a pity to start a LP if the graph is already acyclic
     6270        if ((not self.is_directed() and self.is_forest()) or
     6271            (    self.is_directed() and self.is_directed_acyclic())):
     6272            if value_only:
     6273                return 0
     6274            return []
     6275
     6276        from sage.numerical.mip import MixedIntegerLinearProgram
     6277
     6278        ########################################
     6279        # Constraint Generation Implementation #
     6280        ########################################
     6281        if constraint_generation:
     6282
     6283            p = MixedIntegerLinearProgram(constraint_generation = True,
     6284                                          maximization = False)
     6285
     6286            # A variable for each vertex
     6287            b = p.new_variable(binary = True)
     6288
     6289            # Variables are binary, and their coefficient in the objective is 1
     6290
     6291            p.set_objective(p.sum( b[v] for v in self))
     6292
     6293            p.solve(log = verbose)
     6294
     6295            # For as long as we do not break because the digraph is
     6296            # acyclic....
     6297            while True:
     6298
     6299                # Building the graph without the edges removed by the LP
     6300                h = self.subgraph(vertices =
     6301                                  [v for v in self if p.get_values(b[v]) == 0])
     6302
     6303                # Is the graph acyclic ?
     6304                if self.is_directed():
     6305                    isok, certificate = h.is_directed_acyclic(certificate = True)
     6306                else:
     6307                    isok, certificate = h.is_forest(certificate = True)
     6308
     6309                # If so, we are done !
     6310                if isok:
     6311                    break
     6312
     6313                if verbose:
     6314                    print "Adding a constraint on circuit : ",certificate
     6315
     6316                # There is a circuit left. Let's add the corresponding
     6317                # constraint !
     6318
     6319                p.add_constraint(p.sum(b[v] for v in certificate), min = 1)
     6320
     6321                obj = p.solve(log = verbose)
     6322
     6323            if value_only:
     6324                return obj
     6325
     6326            else:
     6327
     6328                # listing the edges contained in the MFVS
     6329                return [v for v in self if p.get_values(b[v]) == 1]
     6330
     6331        else:
     6332
     6333        ######################################
     6334        # Ordering-based MILP Implementation #
     6335        ######################################
     6336
     6337            p = MixedIntegerLinearProgram(maximization = False, solver = solver)
     6338
     6339            b = p.new_variable(binary = True)
     6340            d = p.new_variable(integer = True)
     6341            n = self.order()
     6342
     6343            # The removed vertices cover all the back arcs ( third condition )
     6344            for (u,v) in self.edges(labels = None):
     6345                p.add_constraint(d[u]-d[v]+n*(b[u]+b[v]), min = 1)
     6346
     6347            for u in self:
     6348                p.add_constraint(d[u], max = n)
     6349
     6350            p.set_objective(p.sum([b[v] for v in self]))
     6351
     6352            if value_only:
     6353                return Integer(round(p.solve(objective_only = True, log = verbose)))
     6354            else:
     6355                p.solve(log=verbose)
     6356                b_sol = p.get_values(b)
     6357
     6358                return [v for v in self if b_sol[v] == 1]
     6359
    61346360    def flow(self, x, y, value_only=True, integer=False, use_edge_labels=True, vertex_bound=False, method = None, solver=None, verbose=0):
    61356361        r"""
    61366362        Returns a maximum flow in the graph from ``x`` to ``y``
     
    61486374            &\forall x\in G, b_x\mbox{ is a binary variable}
    61496375
    61506376        INPUT:
    6151    
     6377
    61526378        - ``x`` -- Source vertex
    61536379
    61546380        - ``y`` -- Sink vertex
  • sage/graphs/graph.py

    diff --git a/sage/graphs/graph.py b/sage/graphs/graph.py
    a b  
    16901690            sage: -1 in cycle
    16911691            True
    16921692
     1693        TESTS:
     1694
     1695        :trac:`14434` is fixed::
     1696
     1697            sage: g = Graph({0:[1,4,5],3:[4,8,9],4:[9],5:[7,8],7:[9]})
     1698            sage: _,cycle = g.is_tree(certificate=True)
     1699            sage: g.size()
     1700            10
     1701            sage: g.add_cycle(cycle)
     1702            sage: g.size()
     1703            10
    16931704        """
    16941705
    16951706        if self.order() == 0:
     
    17071718            n = self.order()
    17081719            seen = {}
    17091720            u = self.vertex_iterator().next()
    1710             v = self.neighbor_iterator(u).next()
    17111721            seen[u] = u
    1712             stack = [(u,v)]
     1722            stack = [(u,v) for v in self.neighbor_iterator(u)]
    17131723            while stack:
    17141724                u,v = stack.pop(-1)
    17151725                if v in seen:
     
    17541764            sage: g.is_forest(certificate = True)
    17551765            (True, None)
    17561766            sage: (2*g + graphs.PetersenGraph() + g).is_forest(certificate = True)
    1757             (False, [60, 61, 62, 63, 64])
     1767            (False, [63, 62, 61, 60, 64])
    17581768        """
    17591769        number_of_connected_components = len(self.connected_components())
    17601770        isit = (self.num_verts() ==