Ticket #9923: trac_9923-python.patch

File trac_9923-python.patch, 21.3 KB (added by ncohen, 9 years ago)
  • sage/graphs/digraph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1285887698 -7200
    # Node ID 0476406f40ad3e24b3912e3df228ba0babe8601e
    # Parent  fc908917889694d6310297aa3536c2a0e8ccfa60
    trac 9923 -- Minimum Feedback Arc/Vertex set through constraint generation
    
    diff -r fc9089178896 -r 0476406f40ad sage/graphs/digraph.py
    a b  
    12231223        return sorted(self.out_degree_iterator(), reverse=True)
    12241224
    12251225
    1226     def feedback_edge_set(self, value_only=False, solver=None, verbose=0):
     1226    def feedback_edge_set(self, constraint_generation= True, value_only=False, solver=None, verbose=0):
    12271227        r"""
    1228         Computes the minimum feedback edge set of a digraph
    1229         (also called feedback arc set).
    1230 
    1231         The minimum feedback edge set of a digraph is a set of edges
    1232         that intersect all the circuits of the digraph.
    1233         Equivalently, a minimum feedback arc set of a DiGraph is a set
    1234         `S` of arcs such that the digraph `G-S` is acyclic. For more
    1235         information, see the
    1236         `Wikipedia article on feedback arc sets
    1237         <http://en.wikipedia.org/wiki/Feedback_arc_set>`_.
     1228        Computes the minimum feedback edge set of a digraph (also called
     1229        feedback arc set).
     1230
     1231        The minimum feedback edge set of a digraph is a set of edges that
     1232        intersect all the circuits of the digraph.  Equivalently, a minimum
     1233        feedback arc set of a DiGraph is a set `S` of arcs such that the digraph
     1234        `G-S` is acyclic. For more information, see the `Wikipedia article on
     1235        feedback arc sets <http://en.wikipedia.org/wiki/Feedback_arc_set>`_.
    12381236
    12391237        INPUT:
    12401238
     
    12431241          - When set to ``True``, only the minimum cardinal of a minimum edge
    12441242            set is returned.
    12451243
    1246           - When set to ``False``, the ``Set`` of edges of a minimal edge set
    1247             is returned.
     1244          - When set to ``False``, the ``Set`` of edges of a minimal edge set is
     1245            returned.
     1246
     1247        - ``constraint_generation`` (boolean) -- whether to use constraint
     1248          generation when solving the Mixed Integer Linear Program (default:
     1249          ``True``).
    12481250
    12491251        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
    12501252          solver to be used. If set to ``None``, the default one is used. For
     
    12571259        - ``verbose`` -- integer (default: ``0``). Sets the level of
    12581260          verbosity. Set to 0 by default, which means quiet.
    12591261
    1260         This problem is solved using Linear Programming, which certainly
    1261         is not the best way and will have to be updated. The program to be
    1262         solved is the following:
     1262        ALGORITHM:
     1263
     1264        This problem is solved using Linear Programming, in two different
     1265        ways. The first one is to solve the following formulation:
    12631266
    12641267        .. MATH::
    12651268
     
    12701273                       
    12711274        An explanation:
    12721275
    1273         An acyclic digraph can be seen as a poset, and every poset has
    1274         a linear extension. This means that in any acyclic digraph
    1275         the vertices can be ordered with a total order `<` in such a way
    1276         that if `(u,v)\in G`, then `u<v`.
    1277 
    1278         Thus, this linear program is built in order to assign to each vertex
    1279         `v` a number `d_v\in [0,\dots,n-1]` such that if there exists
    1280         an edge `(u,v)\in G` such that `d_v<d_u`, then the edge `(u,v)` is
    1281         removed.
    1282 
    1283         The number of edges removed is then minimized, which is
    1284         the objective.
     1276        An acyclic digraph can be seen as a poset, and every poset has a linear
     1277        extension. This means that in any acyclic digraph the vertices can be
     1278        ordered with a total order `<` in such a way that if `(u,v)\in G`, then
     1279        `u<v`.
     1280
     1281        Thus, this linear program is built in order to assign to each vertex `v`
     1282        a number `d_v\in [0,\dots,n-1]` such that if there exists an edge
     1283        `(u,v)\in G` such that `d_v<d_u`, then the edge `(u,v)` is removed.
     1284
     1285        The number of edges removed is then minimized, which is the objective.
     1286
     1287        (Constraint Generation)
     1288 
     1289        If the parameter ``constraint_generation`` is enabled, a more efficient
     1290        formulation is used :
     1291 
     1292        .. MATH::
     1293 
     1294            \mbox{Minimize : }&\sum_{(u,v)\in G} b_{(u,v)}\\ 
     1295            \mbox{Such that : }&\\ 
     1296            &\forall C\text{ circuits }\subseteq G, \sum_{uv\in C}b_{(u,v)}\geq 1\\
     1297 
     1298        As the number of circuits contained in a graph is exponential, this LP
     1299        is solved through constraint generation. This means that the solver is
     1300        sequentially asked to solved the problem, knowing only a portion of the
     1301        circuits contained in `G`, each time adding to the list of its
     1302        constraints the circuit which its last answer had left intact.
    12851303
    12861304        EXAMPLES:
    12871305
    1288         If the digraph is created from a graph, and hence is symmetric
    1289         (if `uv` is an edge, then `vu` is an edge too), then
    1290         obviously the cardinality of its feedback arc set is the number
    1291         of edges in the first graph::
     1306        If the digraph is created from a graph, and hence is symmetric (if `uv`
     1307        is an edge, then `vu` is an edge too), then obviously the cardinality of
     1308        its feedback arc set is the number of edges in the first graph::
    12921309
    12931310            sage: cycle=graphs.CycleGraph(5)
    12941311            sage: dcycle=DiGraph(cycle)
     
    12971314            sage: dcycle.feedback_edge_set(value_only=True)
    12981315            5.0
    12991316       
    1300         And in this situation, for any edge `uv` of the first graph, `uv` of `vu`
    1301         is in the returned feedback arc set::
     1317        And in this situation, for any edge `uv` of the first graph, `uv` of
     1318        `vu` is in the returned feedback arc set::
    13021319
    13031320           sage: g = graphs.RandomGNP(5,.3)
    13041321           sage: dg = DiGraph(g)
     
    13061323           sage: (u,v,l) = g.edge_iterator().next()
    13071324           sage: (u,v) in feedback or (v,u) in feedback
    13081325           True
     1326
     1327        TESTS:
     1328 
     1329        Comparing with/without constraint generation::
     1330 
     1331            sage: g = digraphs.RandomDirectedGNP(10,.3)
     1332            sage: x = g.feedback_edge_set(value_only = True)
     1333            sage: y = g.feedback_edge_set(value_only = True,
     1334            ...          constraint_generation = False)
     1335            sage: x == y
     1336            True
    13091337        """
     1338        # It would be a pity to start a LP if the digraph is already acyclic
     1339        if self.is_directed_acyclic():
     1340            return 0 if value_only else []
     1341
     1342        from sage.numerical.mip import MixedIntegerLinearProgram, Sum
     1343
     1344        ########################################
     1345        # Constraint Generation Implementation #
     1346        ########################################
     1347        if constraint_generation:
     1348
     1349            p = MixedIntegerLinearProgram(constraint_generation = True,
     1350                                          maximization = False)
     1351
     1352            # An variable for each edge
     1353            b = p.new_variable(binary = True, dim = 2)
     1354
     1355            # Variables are binary, and their coefficient in the objective is 1
     1356
     1357            p.set_objective( Sum( b[u][v]
     1358                                  for u,v in self.edges(labels = False)))
     1359
     1360            p.solve(log = verbose)
     1361
     1362            # For as long as we do not break because the digraph is
     1363            # acyclic....
     1364            while (1):
     1365
     1366                # Building the graph without the edges removed by the LP
     1367                h = DiGraph()
     1368                for u,v in self.edges(labels = False):
     1369                    if p.get_values(b[u][v]) < .5:
     1370                        h.add_edge(u,v)
     1371
     1372                # Is the digraph acyclic ?
     1373                isok, certificate = h.is_directed_acyclic(certificate = True)
    13101374       
    1311         from sage.numerical.mip import MixedIntegerLinearProgram, Sum
     1375                # If so, we are done !
     1376                if isok:
     1377                    break
     1378
     1379                if verbose:
     1380                    print "Adding a constraint on circuit : ",certificate
     1381
     1382                # There is a circuit left. Let's add the corresponding
     1383                # constraint !
     1384
     1385                p.add_constraint(
     1386                    Sum( b[u][v] for u,v in
     1387                         zip(certificate, certificate[1:] + [certificate[0]])),
     1388                    min = 1)
     1389               
     1390                obj = p.solve(log = verbose)
     1391
     1392            if value_only:
     1393                return obj
     1394           
     1395            else:
     1396           
     1397                # listing the edges contained in the MFAS
     1398                return [(u,v) for u,v in self.edges(labels = False)
     1399                        if p.get_values(b[u][v]) > .5]
     1400
     1401        ######################################
     1402        # Ordering-based MILP Implementation #
     1403        ######################################
     1404        else:
     1405            p=MixedIntegerLinearProgram(maximization=False, solver=solver)
    13121406       
    1313         p=MixedIntegerLinearProgram(maximization=False, solver=solver)
    1314        
    1315         b=p.new_variable(binary = True)
    1316         d=p.new_variable(integer = True)
    1317 
    1318         n=self.order()
    1319 
    1320         for (u,v) in self.edges(labels=None):
    1321             p.add_constraint(d[u]-d[v]+n*(b[(u,v)]),min=1)
    1322 
    1323         for v in self:
    1324             p.add_constraint(d[v],min=n)
    1325 
    1326 
    1327         p.set_objective(Sum([b[(u,v)] for (u,v) in self.edges(labels=None)]))
    1328 
    1329         if value_only:
    1330             return p.solve(objective_only=True, log=verbose)
    1331         else:
    1332             p.solve(log=verbose)
    1333 
    1334             b_sol=p.get_values(b)
    1335 
    1336             from sage.sets.set import Set
    1337             return Set([(u,v) for (u,v) in self.edges(labels=None) if b_sol[(u,v)]==1])
    1338 
    1339     def feedback_vertex_set(self, value_only=False, solver=None, verbose=0):
     1407            b=p.new_variable(binary = True)
     1408            d=p.new_variable(integer = True)
     1409
     1410            n=self.order()
     1411
     1412            for (u,v) in self.edges(labels=None):
     1413                p.add_constraint(d[u]-d[v]+n*(b[(u,v)]),min=1)
     1414
     1415            for v in self:
     1416                p.add_constraint(d[v],min=n)
     1417
     1418
     1419            p.set_objective(Sum([b[(u,v)] for (u,v) in self.edges(labels=None)]))
     1420
     1421            if value_only:
     1422                return p.solve(objective_only=True, log=verbose)
     1423            else:
     1424                p.solve(log=verbose)
     1425               
     1426                b_sol=p.get_values(b)
     1427               
     1428                from sage.sets.set import Set
     1429                return Set([(u,v) for (u,v) in self.edges(labels=None) if b_sol[(u,v)]==1])
     1430
     1431    def feedback_vertex_set(self, value_only=False, solver=None, verbose=0, constraint_generation = True):
    13401432        r"""
    13411433        Computes the minimum feedback vertex set of a digraph.
    13421434
     
    13521444
    13531445        - ``value_only`` -- boolean (default: ``False``)
    13541446
    1355           - When set to ``True``, only the minimum cardinal of a minimum
     1447          - When set to ``True``, only the minimum cardinal of a minimum vertex
     1448            set is returned.
     1449
     1450          - When set to ``False``, the ``Set`` of vertices of a minimal feedback
    13561451            vertex set is returned.
    13571452
    1358           - When set to ``False``, the ``Set`` of vertices of a minimal
    1359             feedback vertex set is returned.
    1360 
    13611453        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
    13621454          solver to be used. If set to ``None``, the default one is used. For
    13631455          more information on LP solvers and which default solver is used,
     
    13691461        - ``verbose`` -- integer (default: ``0``). Sets the level of
    13701462          verbosity. Set to 0 by default, which means quiet.
    13711463
     1464        - ``constraint_generation`` (boolean) -- whether to use constraint
     1465          generation when solving the Mixed Integer Linear Program (default:
     1466          ``True``).
     1467
    13721468        ALGORITHM:
    13731469
    1374         This problem is solved using Linear Programming, which certainly
    1375         is not the best way and will have to be replaced by a better algorithm.
    1376         The program to be solved is the following:
     1470        This problem is solved using Linear Programming, which certainly is not
     1471        the best way and will have to be replaced by a better algorithm.  The
     1472        program to be solved is the following:
    13771473
    13781474        .. MATH::
    13791475
     
    13841480                       
    13851481        A brief explanation:
    13861482
    1387         An acyclic digraph can be seen as a poset, and every poset has
    1388         a linear extension. This means that in any acyclic digraph the
    1389         vertices can be ordered with a total order `<` in such a way
    1390         that if `(u,v)\in G`, then `u<v`.  Thus, this linear program
    1391         is built in order to assign to each vertex `v` a number
    1392         `d_v\in [0,\dots,n-1]` such that if there exists an edge
    1393         `(u,v)\in G` then either `d_v<d_u` or one of `u` or `v` is
    1394         removed.  The number of vertices removed is then minimized,
    1395         which is the objective.
     1483        An acyclic digraph can be seen as a poset, and every poset has a linear
     1484        extension. This means that in any acyclic digraph the vertices can be
     1485        ordered with a total order `<` in such a way that if `(u,v)\in G`, then
     1486        `u<v`.  Thus, this linear program is built in order to assign to each
     1487        vertex `v` a number `d_v\in [0,\dots,n-1]` such that if there exists an
     1488        edge `(u,v)\in G` then either `d_v<d_u` or one of `u` or `v` is removed.
     1489        The number of vertices removed is then minimized, which is the
     1490        objective.
     1491
     1492        (Constraint Generation)
     1493 
     1494        If the parameter ``constraint_generation`` is enabled, a more efficient
     1495        formulation is used :
     1496 
     1497        .. MATH::
     1498 
     1499            \mbox{Minimize : }&\sum_{v\in G} b_{v}\\ 
     1500            \mbox{Such that : }&\\ 
     1501            &\forall C\text{ circuits }\subseteq G, \sum_{v\in C}b_{v}\geq 1\\
     1502 
     1503        As the number of circuits contained in a graph is exponential, this LP
     1504        is solved through constraint generation. This means that the solver is
     1505        sequentially asked to solved the problem, knowing only a portion of the
     1506        circuits contained in `G`, each time adding to the list of its
     1507        constraints the circuit which its last answer had left intact.
    13961508
    13971509        EXAMPLES:
    13981510
    1399         In a digraph built from a graph, any edge is replaced by arcs going
    1400         in the two opposite directions, thus creating a cycle of length two.
    1401         Hence, to remove all the cycles from the graph, each edge must see
    1402         one of its neighbors removed : a feedback vertex set is in this
    1403         situation a vertex cover::
     1511        In a digraph built from a graph, any edge is replaced by arcs going in
     1512        the two opposite directions, thus creating a cycle of length two.
     1513        Hence, to remove all the cycles from the graph, each edge must see one
     1514        of its neighbors removed : a feedback vertex set is in this situation a
     1515        vertex cover::
    14041516
    14051517            sage: cycle=graphs.CycleGraph(5)
    14061518            sage: dcycle=DiGraph(cycle)
     
    14181530            sage: circuit = digraphs.Circuit(5)
    14191531            sage: circuit.feedback_vertex_set(value_only=True) == 1
    14201532            True
    1421         """
    1422        
     1533
     1534        TESTS:
     1535 
     1536        Comparing with/without constraint generation::
     1537 
     1538            sage: g = digraphs.RandomDirectedGNP(10,.3)
     1539            sage: x = g.feedback_vertex_set(value_only = True)
     1540            sage: y = g.feedback_vertex_set(value_only = True,
     1541            ...            constraint_generation = False)
     1542            sage: x == y
     1543            True
     1544         """
     1545
     1546        # It would be a pity to start a LP if the digraph is already acyclic
     1547        if self.is_directed_acyclic():
     1548            if value_only:
     1549                return 0
     1550
     1551            from sage.sets.set import Set
     1552            return Set([0])
     1553
    14231554        from sage.numerical.mip import MixedIntegerLinearProgram, Sum
    1424        
    1425         p=MixedIntegerLinearProgram(maximization=False, solver=solver)
    1426        
    1427         b=p.new_variable(binary = True)
    1428         d=p.new_variable(integer = True)
    1429         n=self.order()
    1430 
    1431         # The removed vertices cover all the back arcs ( third condition )
    1432         for (u, v) in self.edges(labels=False):
    1433             p.add_constraint(d[u] - d[v] + n*(b[u] + b[v]), min=1)
    1434         for u in self:
    1435             p.add_constraint(d[u], max=n)
    1436 
    1437         p.set_objective(Sum([b[v] for v in self]))
    1438 
    1439         if value_only:
    1440             return p.solve(objective_only=True, log=verbose)
     1555
     1556        ########################################
     1557        # Constraint Generation Implementation #
     1558        ########################################
     1559        if constraint_generation:
     1560
     1561            p = MixedIntegerLinearProgram(constraint_generation = True,
     1562                                          maximization = False)
     1563
     1564            # An variable for each vertex
     1565            b = p.new_variable(binary = True)
     1566
     1567            # Variables are binary, and their coefficient in the objective is 1
     1568
     1569            p.set_objective( Sum( b[v] for v in self))
     1570
     1571            p.solve(log = verbose)
     1572
     1573            # For as long as we do not break because the digraph is
     1574            # acyclic....
     1575            while (1):
     1576
     1577                # Building the graph without the edges removed by the LP
     1578                h = self.subgraph(vertices =
     1579                                  [v for v in self if p.get_values(b[v]) < .5])
     1580
     1581                # Is the digraph acyclic ?
     1582                isok, certificate = h.is_directed_acyclic(certificate = True)
     1583
     1584                # If so, we are done !
     1585                if isok:
     1586                    break
     1587
     1588                if verbose:
     1589                    print "Adding a constraint on circuit : ",certificate
     1590
     1591                # There is a circuit left. Let's add the corresponding
     1592                # constraint !
     1593
     1594                p.add_constraint( Sum( b[v] for v in certificate), min = 1)
     1595               
     1596                obj = p.solve(log = verbose)
     1597
     1598            if value_only:
     1599                return obj
     1600           
     1601            else:
     1602           
     1603                # listing the edges contained in the MFAS
     1604                from sage.sets.set import Set
     1605                return Set([v for v in self if p.get_values(b[v]) > .5])
     1606
     1607
    14411608        else:
    1442             p.solve(log=verbose)
    1443             b_sol=p.get_values(b)
    1444 
    1445             from sage.sets.set import Set
    1446             return Set([v for v in self if b_sol[v]==1])
     1609
     1610        ######################################
     1611        # Ordering-based MILP Implementation #
     1612        ######################################
     1613
     1614            p = MixedIntegerLinearProgram(maximization=False, solver=solver)
     1615           
     1616            b = p.new_variable(binary = True)
     1617            d = p.new_variable(integer = True)
     1618            n = self.order()
     1619
     1620            # The removed vertices cover all the back arcs ( third condition )
     1621            for (u,v) in self.edges(labels=None):
     1622                p.add_constraint(d[u]-d[v]+n*(b[u]+b[v]),min=1)
     1623
     1624            for u in self:
     1625                p.add_constraint(d[u],max=n)
     1626
     1627            p.set_objective(Sum([b[v] for v in self]))
     1628
     1629            if value_only:
     1630                return p.solve(objective_only=True, log=verbose)
     1631            else:
     1632                p.solve(log=verbose)
     1633                b_sol=p.get_values(b)
     1634
     1635                from sage.sets.set import Set
     1636                return Set([v for v in self if b_sol[v]==1])
    14471637
    14481638
    14491639    ### Construction
  • sage/graphs/generic_graph_pyx.pyx

    diff -r fc9089178896 -r 0476406f40ad sage/graphs/generic_graph_pyx.pyx
    a b  
    949949        sage_free(u)
    950950        sage_free(v)
    951951
    952 
    953952cpdef tuple find_hamiltonian( G, long max_iter=100000, long reset_bound=30000, long backtrack_bound=1000, find_path=False ):
    954953    r"""
    955954    Randomized backtracking for finding hamiltonian cycles and paths.
  • sage/numerical/mip.pyx

    diff -r fc9089178896 -r 0476406f40ad sage/numerical/mip.pyx
    a b  
    130130      - When set to ``False``, the ``MixedIntegerLinearProgram`` is
    131131        defined as a minimization.
    132132
     133    - ``constraint_generation`` -- whether to require the returned solver to
     134      support constraint generation (excludes Coin). ``False by default``.
     135
    133136    .. SEEALSO::
    134137
    135138     - :func:`default_mip_solver` -- Returns/Sets the default MIP solver.
     
    149152         4.0
    150153    """
    151154   
    152     def __init__(self, solver = None, maximization=True):
     155    def __init__(self, solver = None, maximization=True, constraint_generation = False):
    153156        r"""
    154157        Constructor for the ``MixedIntegerLinearProgram`` class.
    155158
     
    178181          - When set to ``False``, the ``MixedIntegerLinearProgram`` is
    179182            defined as a minimization.
    180183
     184        - ``constraint_generation`` -- whether to require the returned solver to
     185          support constraint generation (excludes Coin). ``False by default``.
     186
    181187        .. SEEALSO::
    182188
    183189        - :meth:`default_mip_solver` -- Returns/Sets the default MIP solver.
     
    189195        """
    190196
    191197        from sage.numerical.backends.generic_backend import get_solver
    192         self._backend = get_solver(solver = solver)
     198
     199        self._backend = get_solver(solver=solver,
     200                                   constraint_generation=constraint_generation)
    193201
    194202        if not maximization:
    195203            self._backend.set_sense(-1)