Ticket #2203: trac_2203.patch

File trac_2203.patch, 12.1 KB (added by ncohen, 11 years ago)
  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1265713488 -3600
    # Node ID 59390123f8f939884d08c7565c32812a31099a3b
    # Parent  eb27a39a6df43f557bdc8bc838c30f3b5c41e925
    trac 2203 : Graph.travelling_salesman_problem + is_hamiltonian
    
    diff -r eb27a39a6df4 -r 59390123f8f9 sage/graphs/generic_graph.py
    a b  
    36383638
    36393639            return val
    36403640
     3641    def traveling_salesman_problem(self, weighted = True):
     3642        r"""
     3643        Solves the traveling salesman problem (TSP)
     3644
     3645        Given a graph (resp. a digraph) `G` with weighted edges,
     3646        the traveling salesman problem consists in finding a
     3647        hamiltonian cycle (resp. circuit) of the graph of
     3648        minimum cost.
     3649
     3650        This TSP is one of the most famous NP-Complete problems,
     3651        this function can thus be expected to take some time
     3652        before returning its result.
     3653
     3654        INPUT:
     3655
     3656        - ``weighted`` (boolean) -- whether to consider the
     3657          weights of the edges.
     3658
     3659              - If set to ``False`` (default), all edges are
     3660                assumed to weight `1`
     3661
     3662              - If set to ``True``, the weights are taken into
     3663                account, and the edges whose weight is ``None``
     3664                are assumed to be set to `1`
     3665
     3666        OUTPUT:
     3667
     3668        A solution to the TSP, as a ``Graph`` object whose vertex
     3669        set is `V(G)`, and whose edges are only those of the
     3670        solution.
     3671
     3672        ALGORITHM:
     3673
     3674        This optimization problem is solved through the use
     3675        of Linear Programming.
     3676       
     3677        NOTE:
     3678
     3679        - This function is correctly defined for both graph and digraphs.
     3680          In the second case, the returned cycle is a circuit of optimal
     3681          cost.
     3682
     3683        EXAMPLES:
     3684
     3685        The Heawood graph is known to be hamiltonian::
     3686
     3687            sage: g = graphs.HeawoodGraph()
     3688            sage: tsp = g.traveling_salesman_problem()   # optional - requires GLPK, CBC, or CPLEX
     3689            sage: tsp                                     # optional - requires GLPK, CBC, or CPLEX
     3690            TSP from Heawood graph: Graph on 14 vertices
     3691
     3692        The solution to the TSP has to be connected ::
     3693
     3694            sage: tsp.is_connected()                      # optional - requires GLPK, CBC, or CPLEX
     3695            True
     3696
     3697        It must also be a `2`-regular graph::
     3698
     3699            sage: tsp.is_regular(k=2)                     # optional - requires GLPK, CBC, or CPLEX
     3700            True
     3701
     3702        And obviously it is a subgraph of the Heawood graph::
     3703
     3704            sage: all([ e in g.edges() for e in tsp.edges()])    # optional - requires GLPK, CBC, or CPLEX
     3705            True
     3706
     3707        On the other hand, the Petersen Graph is known not to
     3708        be hamiltonian::
     3709
     3710            sage: g = graphs.PetersenGraph()                     # optional - requires GLPK, CBC, or CPLEX
     3711            sage: tsp = g.traveling_salesman_problem()          # optional - requires GLPK, CBC, or CPLEX
     3712            Traceback (most recent call last):
     3713            ...
     3714            ValueError: The given graph is not hamiltonian
     3715
     3716        One easy way to change is is obviously to add to this
     3717        graph the edges corresponding to a hamiltonian cycle.
     3718     
     3719        If we do this by setting the cost of these new edges
     3720        to `2`, while the others are set to `1`, we notice
     3721        that not all the edges we added are used in the
     3722        optimal solution ::
     3723
     3724            sage: for u, v in g.edges(labels = None):
     3725            ...      g.set_edge_label(u,v,1)
     3726
     3727            sage: cycle = graphs.CycleGraph(10)
     3728            sage: for u,v in cycle.edges(labels = None):
     3729            ...      if not g.has_edge(u,v):
     3730            ...          g.add_edge(u,v)
     3731            ...      g.set_edge_label(u,v,2)
     3732
     3733            sage: tsp = g.traveling_salesman_problem(weighted = True)   # optional - requires GLPK, CBC, or CPLEX
     3734            sage: sum( tsp.edge_labels() ) < 2*10                       # optional - requires GLPK, CBC, or CPLEX
     3735            True
     3736
     3737        If we pick `1/2` instead of `2` as a cost for these new edges,
     3738        they clearly become the optimal solution
     3739
     3740            sage: for u,v in cycle.edges(labels = None):
     3741            ...      g.set_edge_label(u,v,1/2)
     3742
     3743            sage: tsp = g.traveling_salesman_problem(weighted = True)   # optional - requires GLPK, CBC, or CPLEX
     3744            sage: sum( tsp.edge_labels() ) == (1/2)*10                   # optional - requires GLPK, CBC, or CPLEX
     3745            True
     3746
     3747        """
     3748       
     3749        from sage.numerical.mip import MixedIntegerLinearProgram
     3750
     3751        p = MixedIntegerLinearProgram(maximization = False)
     3752
     3753
     3754        f = p.new_variable()
     3755        r = p.new_variable()
     3756
     3757
     3758
     3759        # If the graph has multiple edges
     3760        if self.has_multiple_edges():
     3761            g = self.copy()
     3762            multi = self.multiple_edges()
     3763            g.delete_edges(multi)
     3764            g.allow_multiple_edges(False)
     3765            if weighted:
     3766                e = {}
     3767
     3768                for u,v,l in multi:
     3769                    u,v = (u,v) if u<v else (v,u)
     3770
     3771                    # The weight of an edge is the minimum over
     3772                    # the weights of the parallel edges
     3773
     3774                    #  new value *if* ( none other        *or*   new==None and last > 1     *else*  change nothing
     3775                    e[(u,v)] = l if (not e.has_key((u,v)) or ( l is None and e[(u,v)] > 1 )) else e[(u,v)]
     3776
     3777                g.add_edges([(u,v) for (u,v),l in e.iteritems()])
     3778
     3779            else:
     3780                from sage.sets.set import Set
     3781                g.add_edges(Set([ (min(u,v),max(u,v)) for u,v,l in multi]))
     3782
     3783        else:
     3784            g = self
     3785
     3786        eps = 1/(2*Integer(g.order()))
     3787        x = g.vertex_iterator().next()
     3788
     3789
     3790        if g.is_directed():
     3791
     3792            # returns the variable corresponding to arc u,v
     3793            E = lambda u,v : f[(u,v)]
     3794
     3795            # All the vertices have in-degree 1
     3796            [p.add_constraint(sum([ f[(u,v)] for u in g.neighbors_in(v)]), min = 1, max = 1) for v in g]
     3797
     3798            # All the vertices have out-degree 1
     3799            [p.add_constraint(sum([ f[(v,u)] for u in g.neighbors_out(v)]), min = 1, max = 1) for v in g]
     3800
     3801
     3802            # r is greater than f
     3803            for u,v in g.edges(labels = None):
     3804                if g.has_edge(v,u):
     3805                    if u < v:
     3806                        p.add_constraint( r[(u,v)] + r[(v,u)]- f[(u,v)] - f[(v,u)], min = 0)
     3807
     3808                        # no 2-cycles
     3809                        p.add_constraint( f[(u,v)] + f[(v,u)], max = 1)
     3810
     3811                else:
     3812                    p.add_constraint( r[(u,v)] + r[(v,u)] - f[(u,v)], min = 0)
     3813                   
     3814
     3815            # defining the answer when g is directed
     3816            from sage.graphs.all import DiGraph
     3817            tsp = DiGraph()
     3818        else:
     3819
     3820            # reorders the edge as they can appear in the two different ways
     3821            R = lambda x,y : (x,y) if x < y else (y,x)
     3822
     3823            # returns the variable corresponding to arc u,v
     3824            E = lambda u,v : f[R(u,v)]
     3825
     3826            # All the vertices have degree 2
     3827            [p.add_constraint(sum([ f[R(u,v)] for u in g.neighbors(v)]), min = 2, max = 2) for v in g]
     3828
     3829            # r is greater than f
     3830            for u,v in g.edges(labels = None):
     3831                p.add_constraint( r[(u,v)] + r[(v,u)] - f[R(u,v)], min = 0)
     3832
     3833            from sage.graphs.all import Graph
     3834
     3835            # defining the answer when g is not directed
     3836            tsp = Graph()
     3837
     3838               
     3839        # no cycle which does not contain x
     3840        for v in g:
     3841            if v != x:
     3842                p.add_constraint( sum([ r[(u,v)] for u in g.neighbors(v)]),max = 1-eps)
     3843
     3844
     3845        weight = lambda u,v : g.edge_label(u,v) if g.edge_label(u,v) is not None else 1
     3846
     3847        if weighted:
     3848            p.set_objective( sum([ weight(u,v)*E(u,v) for u,v in g.edges(labels=None)]) )
     3849        else:
     3850            p.set_objective(None)
     3851
     3852        p.set_binary(f)
     3853
     3854        from sage.numerical.mip import MIPSolverException
     3855
     3856        try:
     3857            obj = p.solve()
     3858            f = p.get_values(f)
     3859            tsp.add_vertices(g.vertices())
     3860            tsp.set_pos(g.get_pos())
     3861            tsp.name("TSP from "+g.name())
     3862            tsp.add_edges([(u,v,l) for u,v,l in g.edges() if E(u,v) == 1])
     3863
     3864            return tsp
     3865       
     3866        except MIPSolverException:
     3867            raise ValueError("The given graph is not hamiltonian")
     3868
     3869
     3870    def hamiltonian_cycle(self):
     3871        r"""
     3872        Returns a hamiltonian cycle/circuit of the current graph/digraph
     3873
     3874        A graph (resp. digraph) is said to be hamiltonian
     3875        if it contains as a subgraph a cycle (resp. a circuit)
     3876        going through all the vertices.
     3877
     3878        Computing a hamiltonian cycle/circuit being NP-Complete,
     3879        this algorithm could run for some time depending on
     3880        the instance.
     3881
     3882        ALGORITHM:
     3883
     3884        See ``Graph.traveling_salesman_problem``.
     3885
     3886        OUTPUT:
     3887
     3888        Returns a hamiltonian cycle/circuit if it exists. Otherwise,
     3889        raises a ``ValueError`` exception.
     3890
     3891        NOTE:
     3892
     3893        This function, as ``is_hamiltonian``, computes a hamiltonian
     3894        cycle if it exists : the user should *NOT* test for
     3895        hamiltonicity using ``is_hamiltonian`` before calling this
     3896        function, as it would result in computing it twice.
     3897
     3898        EXAMPLES:
     3899
     3900        The Heawood Graph is known to be hamiltonian ::
     3901
     3902            sage: g = graphs.HeawoodGraph()
     3903            sage: g.hamiltonian_cycle()         # optional - requires GLPK, CBC, or CPLEX
     3904            TSP from Heawood graph: Graph on 14 vertices
     3905
     3906        The Petergraph, though, is not ::
     3907
     3908            sage: g = graphs.PetersenGraph()
     3909            sage: g.hamiltonian_cycle()         # optional - requires GLPK, CBC, or CPLEX
     3910            Traceback (most recent call last):
     3911            ...
     3912            ValueError: The given graph is not hamiltonian
     3913        """
     3914
     3915        try:
     3916            return self.traveling_salesman_problem(weighted = False)
     3917        except:
     3918            raise ValueError("The given graph is not hamiltonian")
     3919
    36413920    def flow(self,x,y,value_only=True,integer=False, use_edge_labels=True,vertex_bound=False):
    36423921        r"""
    36433922        Returns a maximum flow in the graph from ``x`` to ``y``
     
    1048510764                        return False
    1048610765        return True
    1048710766
     10767    def is_hamiltonian(self):
     10768        r"""
     10769        Tests whether the current graph is Hamiltonian
     10770
     10771        A graph (resp. digraph) is said to be hamiltonian
     10772        if it contains as a subgraph a cycle (resp. a circuit)
     10773        going through all the vertices.
     10774
     10775        Testing for hamiltonicity being NP-Complete, this
     10776        algorithm could run for some time depending on
     10777        the instance.
     10778
     10779        ALGORITHM:
     10780
     10781        See ``Graph.traveling_salesman_problem``.
     10782
     10783        OUTPUT:
     10784
     10785        Returns ``True`` if a hamiltonian cycle/circuit exists, and
     10786        ``False`` otherwise.
     10787
     10788        NOTE:
     10789
     10790        This function, as ``hamiltonian_cycle`` and
     10791        ``traveling_salesman_problem``, computes a hamiltonian
     10792        cycle if it exists : the user should *NOT* test for
     10793        hamiltonicity using ``is_hamiltonian`` before calling
     10794        ``hamiltonian_cycle`` or ``traveling_salesman_problem``
     10795        as it would result in computing it twice.
     10796
     10797        EXAMPLES:
     10798
     10799        The Heawood Graph is known to be hamiltonian ::
     10800
     10801            sage: g = graphs.HeawoodGraph()
     10802            sage: g.is_hamiltonian()         # optional - requires GLPK, CBC, or CPLEX
     10803            True
     10804
     10805        The Petergraph, though, is not ::
     10806
     10807            sage: g = graphs.PetersenGraph()
     10808            sage: g.is_hamiltonian()         # optional - requires GLPK, CBC, or CPLEX
     10809            False
     10810
     10811        """
     10812
     10813        try:
     10814            tsp = self.traveling_salesman_problem(weighted = False)
     10815            return True
     10816
     10817        except ValueError:
     10818            return False
     10819
    1048810820    def is_isomorphic(self, other, certify=False, verbosity=0, edge_labels=False):
    1048910821        """
    1049010822        Tests for isomorphism between self and other.