Ticket #2203: trac_2203-rebase.patch

File trac_2203-rebase.patch, 12.1 KB (added by mvngu, 11 years ago)

rebase of previous patch

  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1274523674 25200
    # Node ID a9fbd55b2bef9829d1205997e5f97de2b5c615d4
    # Parent  35c219593994239c7763d329c2f510487bb0da01
    #2203: a traveling salesman problem solver; is_hamiltonian method
    
    diff --git a/sage/graphs/generic_graph.py b/sage/graphs/generic_graph.py
    a b  
    35903590
    35913591            return val
    35923592
     3593    def traveling_salesman_problem(self, weighted = True):
     3594        r"""
     3595        Solves the traveling salesman problem (TSP)
     3596
     3597        Given a graph (resp. a digraph) `G` with weighted edges,
     3598        the traveling salesman problem consists in finding a
     3599        hamiltonian cycle (resp. circuit) of the graph of
     3600        minimum cost.
     3601
     3602        This TSP is one of the most famous NP-Complete problems,
     3603        this function can thus be expected to take some time
     3604        before returning its result.
     3605
     3606        INPUT:
     3607
     3608        - ``weighted`` (boolean) -- whether to consider the
     3609          weights of the edges.
     3610
     3611              - If set to ``False`` (default), all edges are
     3612                assumed to weight `1`
     3613
     3614              - If set to ``True``, the weights are taken into
     3615                account, and the edges whose weight is ``None``
     3616                are assumed to be set to `1`
     3617
     3618        OUTPUT:
     3619
     3620        A solution to the TSP, as a ``Graph`` object whose vertex
     3621        set is `V(G)`, and whose edges are only those of the
     3622        solution.
     3623
     3624        ALGORITHM:
     3625
     3626        This optimization problem is solved through the use
     3627        of Linear Programming.
     3628       
     3629        NOTE:
     3630
     3631        - This function is correctly defined for both graph and digraphs.
     3632          In the second case, the returned cycle is a circuit of optimal
     3633          cost.
     3634
     3635        EXAMPLES:
     3636
     3637        The Heawood graph is known to be hamiltonian::
     3638
     3639            sage: g = graphs.HeawoodGraph()
     3640            sage: tsp = g.traveling_salesman_problem()   # optional - requires GLPK, CBC, or CPLEX
     3641            sage: tsp                                     # optional - requires GLPK, CBC, or CPLEX
     3642            TSP from Heawood graph: Graph on 14 vertices
     3643
     3644        The solution to the TSP has to be connected ::
     3645
     3646            sage: tsp.is_connected()                      # optional - requires GLPK, CBC, or CPLEX
     3647            True
     3648
     3649        It must also be a `2`-regular graph::
     3650
     3651            sage: tsp.is_regular(k=2)                     # optional - requires GLPK, CBC, or CPLEX
     3652            True
     3653
     3654        And obviously it is a subgraph of the Heawood graph::
     3655
     3656            sage: all([ e in g.edges() for e in tsp.edges()])    # optional - requires GLPK, CBC, or CPLEX
     3657            True
     3658
     3659        On the other hand, the Petersen Graph is known not to
     3660        be hamiltonian::
     3661
     3662            sage: g = graphs.PetersenGraph()                     # optional - requires GLPK, CBC, or CPLEX
     3663            sage: tsp = g.traveling_salesman_problem()          # optional - requires GLPK, CBC, or CPLEX
     3664            Traceback (most recent call last):
     3665            ...
     3666            ValueError: The given graph is not hamiltonian
     3667
     3668        One easy way to change is is obviously to add to this
     3669        graph the edges corresponding to a hamiltonian cycle.
     3670     
     3671        If we do this by setting the cost of these new edges
     3672        to `2`, while the others are set to `1`, we notice
     3673        that not all the edges we added are used in the
     3674        optimal solution ::
     3675
     3676            sage: for u, v in g.edges(labels = None):
     3677            ...      g.set_edge_label(u,v,1)
     3678
     3679            sage: cycle = graphs.CycleGraph(10)
     3680            sage: for u,v in cycle.edges(labels = None):
     3681            ...      if not g.has_edge(u,v):
     3682            ...          g.add_edge(u,v)
     3683            ...      g.set_edge_label(u,v,2)
     3684
     3685            sage: tsp = g.traveling_salesman_problem(weighted = True)   # optional - requires GLPK, CBC, or CPLEX
     3686            sage: sum( tsp.edge_labels() ) < 2*10                       # optional - requires GLPK, CBC, or CPLEX
     3687            True
     3688
     3689        If we pick `1/2` instead of `2` as a cost for these new edges,
     3690        they clearly become the optimal solution
     3691
     3692            sage: for u,v in cycle.edges(labels = None):
     3693            ...      g.set_edge_label(u,v,1/2)
     3694
     3695            sage: tsp = g.traveling_salesman_problem(weighted = True)   # optional - requires GLPK, CBC, or CPLEX
     3696            sage: sum( tsp.edge_labels() ) == (1/2)*10                   # optional - requires GLPK, CBC, or CPLEX
     3697            True
     3698
     3699        """
     3700       
     3701        from sage.numerical.mip import MixedIntegerLinearProgram
     3702
     3703        p = MixedIntegerLinearProgram(maximization = False)
     3704
     3705
     3706        f = p.new_variable()
     3707        r = p.new_variable()
     3708
     3709
     3710
     3711        # If the graph has multiple edges
     3712        if self.has_multiple_edges():
     3713            g = self.copy()
     3714            multi = self.multiple_edges()
     3715            g.delete_edges(multi)
     3716            g.allow_multiple_edges(False)
     3717            if weighted:
     3718                e = {}
     3719
     3720                for u,v,l in multi:
     3721                    u,v = (u,v) if u<v else (v,u)
     3722
     3723                    # The weight of an edge is the minimum over
     3724                    # the weights of the parallel edges
     3725
     3726                    #  new value *if* ( none other        *or*   new==None and last > 1     *else*  change nothing
     3727                    e[(u,v)] = l if (not e.has_key((u,v)) or ( l is None and e[(u,v)] > 1 )) else e[(u,v)]
     3728
     3729                g.add_edges([(u,v) for (u,v),l in e.iteritems()])
     3730
     3731            else:
     3732                from sage.sets.set import Set
     3733                g.add_edges(Set([ (min(u,v),max(u,v)) for u,v,l in multi]))
     3734
     3735        else:
     3736            g = self
     3737
     3738        eps = 1/(2*Integer(g.order()))
     3739        x = g.vertex_iterator().next()
     3740
     3741
     3742        if g.is_directed():
     3743
     3744            # returns the variable corresponding to arc u,v
     3745            E = lambda u,v : f[(u,v)]
     3746
     3747            # All the vertices have in-degree 1
     3748            [p.add_constraint(sum([ f[(u,v)] for u in g.neighbors_in(v)]), min = 1, max = 1) for v in g]
     3749
     3750            # All the vertices have out-degree 1
     3751            [p.add_constraint(sum([ f[(v,u)] for u in g.neighbors_out(v)]), min = 1, max = 1) for v in g]
     3752
     3753
     3754            # r is greater than f
     3755            for u,v in g.edges(labels = None):
     3756                if g.has_edge(v,u):
     3757                    if u < v:
     3758                        p.add_constraint( r[(u,v)] + r[(v,u)]- f[(u,v)] - f[(v,u)], min = 0)
     3759
     3760                        # no 2-cycles
     3761                        p.add_constraint( f[(u,v)] + f[(v,u)], max = 1)
     3762
     3763                else:
     3764                    p.add_constraint( r[(u,v)] + r[(v,u)] - f[(u,v)], min = 0)
     3765                   
     3766
     3767            # defining the answer when g is directed
     3768            from sage.graphs.all import DiGraph
     3769            tsp = DiGraph()
     3770        else:
     3771
     3772            # reorders the edge as they can appear in the two different ways
     3773            R = lambda x,y : (x,y) if x < y else (y,x)
     3774
     3775            # returns the variable corresponding to arc u,v
     3776            E = lambda u,v : f[R(u,v)]
     3777
     3778            # All the vertices have degree 2
     3779            [p.add_constraint(sum([ f[R(u,v)] for u in g.neighbors(v)]), min = 2, max = 2) for v in g]
     3780
     3781            # r is greater than f
     3782            for u,v in g.edges(labels = None):
     3783                p.add_constraint( r[(u,v)] + r[(v,u)] - f[R(u,v)], min = 0)
     3784
     3785            from sage.graphs.all import Graph
     3786
     3787            # defining the answer when g is not directed
     3788            tsp = Graph()
     3789
     3790       
     3791        # no cycle which does not contain x
     3792        for v in g:
     3793            if v != x:
     3794                p.add_constraint( sum([ r[(u,v)] for u in g.neighbors(v)]),max = 1-eps)
     3795
     3796
     3797        weight = lambda u,v : g.edge_label(u,v) if g.edge_label(u,v) is not None else 1
     3798
     3799        if weighted:
     3800            p.set_objective( sum([ weight(u,v)*E(u,v) for u,v in g.edges(labels=None)]) )
     3801        else:
     3802            p.set_objective(None)
     3803
     3804        p.set_binary(f)
     3805
     3806        from sage.numerical.mip import MIPSolverException
     3807
     3808        try:
     3809            obj = p.solve()
     3810            f = p.get_values(f)
     3811            tsp.add_vertices(g.vertices())
     3812            tsp.set_pos(g.get_pos())
     3813            tsp.name("TSP from "+g.name())
     3814            tsp.add_edges([(u,v,l) for u,v,l in g.edges() if E(u,v) == 1])
     3815
     3816            return tsp
     3817       
     3818        except MIPSolverException:
     3819            raise ValueError("The given graph is not hamiltonian")
     3820
     3821
     3822    def hamiltonian_cycle(self):
     3823        r"""
     3824        Returns a hamiltonian cycle/circuit of the current graph/digraph
     3825
     3826        A graph (resp. digraph) is said to be hamiltonian
     3827        if it contains as a subgraph a cycle (resp. a circuit)
     3828        going through all the vertices.
     3829
     3830        Computing a hamiltonian cycle/circuit being NP-Complete,
     3831        this algorithm could run for some time depending on
     3832        the instance.
     3833
     3834        ALGORITHM:
     3835
     3836        See ``Graph.traveling_salesman_problem``.
     3837
     3838        OUTPUT:
     3839
     3840        Returns a hamiltonian cycle/circuit if it exists. Otherwise,
     3841        raises a ``ValueError`` exception.
     3842
     3843        NOTE:
     3844
     3845        This function, as ``is_hamiltonian``, computes a hamiltonian
     3846        cycle if it exists : the user should *NOT* test for
     3847        hamiltonicity using ``is_hamiltonian`` before calling this
     3848        function, as it would result in computing it twice.
     3849
     3850        EXAMPLES:
     3851
     3852        The Heawood Graph is known to be hamiltonian ::
     3853
     3854            sage: g = graphs.HeawoodGraph()
     3855            sage: g.hamiltonian_cycle()         # optional - requires GLPK, CBC, or CPLEX
     3856            TSP from Heawood graph: Graph on 14 vertices
     3857
     3858        The Petergraph, though, is not ::
     3859
     3860            sage: g = graphs.PetersenGraph()
     3861            sage: g.hamiltonian_cycle()         # optional - requires GLPK, CBC, or CPLEX
     3862            Traceback (most recent call last):
     3863            ...
     3864            ValueError: The given graph is not hamiltonian
     3865        """
     3866
     3867        try:
     3868            return self.traveling_salesman_problem(weighted = False)
     3869        except:
     3870            raise ValueError("The given graph is not hamiltonian")
     3871
    35933872    def flow(self, x, y, value_only=True, integer=False, use_edge_labels=True, vertex_bound=False, solver=None, verbose=0):
    35943873        r"""
    35953874        Returns a maximum flow in the graph from ``x`` to ``y``
     
    1140611685                        return False
    1140711686        return True
    1140811687
     11688    def is_hamiltonian(self):
     11689        r"""
     11690        Tests whether the current graph is Hamiltonian
     11691
     11692        A graph (resp. digraph) is said to be hamiltonian
     11693        if it contains as a subgraph a cycle (resp. a circuit)
     11694        going through all the vertices.
     11695
     11696        Testing for hamiltonicity being NP-Complete, this
     11697        algorithm could run for some time depending on
     11698        the instance.
     11699
     11700        ALGORITHM:
     11701
     11702        See ``Graph.traveling_salesman_problem``.
     11703
     11704        OUTPUT:
     11705
     11706        Returns ``True`` if a hamiltonian cycle/circuit exists, and
     11707        ``False`` otherwise.
     11708
     11709        NOTE:
     11710
     11711        This function, as ``hamiltonian_cycle`` and
     11712        ``traveling_salesman_problem``, computes a hamiltonian
     11713        cycle if it exists : the user should *NOT* test for
     11714        hamiltonicity using ``is_hamiltonian`` before calling
     11715        ``hamiltonian_cycle`` or ``traveling_salesman_problem``
     11716        as it would result in computing it twice.
     11717
     11718        EXAMPLES:
     11719
     11720        The Heawood Graph is known to be hamiltonian ::
     11721
     11722            sage: g = graphs.HeawoodGraph()
     11723            sage: g.is_hamiltonian()         # optional - requires GLPK, CBC, or CPLEX
     11724            True
     11725
     11726        The Petergraph, though, is not ::
     11727
     11728            sage: g = graphs.PetersenGraph()
     11729            sage: g.is_hamiltonian()         # optional - requires GLPK, CBC, or CPLEX
     11730            False
     11731
     11732        """
     11733
     11734        try:
     11735            tsp = self.traveling_salesman_problem(weighted = False)
     11736            return True
     11737
     11738        except ValueError:
     11739            return False
     11740
    1140911741    def is_isomorphic(self, other, certify=False, verbosity=0, edge_labels=False):
    1141011742        """
    1141111743        Tests for isomorphism between self and other.