Ticket #9910: trac_9910.patch

File trac_9910.patch, 14.4 KB (added by ncohen, 10 years ago)
  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1284491936 -7200
    # Node ID 9c9495ebab27b9b708c83c50bf7b3d88c0a903db
    # Parent  e369cc8bfa30fbab80056966a4d484c1eb2766f6
    trac 9910 -- Longest Path method for graphs/digraphs
    
    diff -r e369cc8bfa30 -r 9c9495ebab27 sage/graphs/generic_graph.py
    a b  
    43744374
    43754375            return val
    43764376
     4377    def longest_path(self, s = None, t = None, weighted = False, algorithm = "MILP", solver = None, verbose = 0):
     4378        r"""
     4379        Returns a longest path of ``self``.
     4380
     4381        INPUT:
     4382
     4383        - ``s`` (vertex) -- forces the source of the path. Set to
     4384          ``None`` by default.
     4385
     4386        - ``t`` (vertex) -- forces the destination of the path. Set to
     4387          ``None`` by default.
     4388
     4389        - ``weighted`` (boolean) -- whether the label on the edges are
     4390          tobe considered as weights (a label set to ``None`` or
     4391          ``{}`` being considered as a weight of `1`). Set to
     4392          ``False`` by default.
     4393
     4394        - ``algorithm`` -- one of ``"MILP"`` (default) or
     4395          ``"backtrack"``. Two remarks on this respect :
     4396
     4397              * While the MILP formulation returns an exact answer,
     4398                the backtrack algorithm is a randomized heuristic.
     4399
     4400              * As the backtrack algorithm does not support edge
     4401                weighting, setting ``weighted`` to ``True`` will force
     4402                the use of the MILP algorithm.
     4403
     4404        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
     4405          solver to be used. If set to ``None``, the default one is used. For
     4406          more information on LP solvers and which default solver is used, see
     4407          the method
     4408          :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
     4409          of the class
     4410          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.
     4411
     4412        - ``verbose`` -- integer (default: ``0``). Sets the level of
     4413          verbosity. Set to 0 by default, which means quiet.
     4414
     4415        .. NOTE::
     4416
     4417            The length of a path is assumed to be the number of its
     4418            edges, or the sum of their labels.
     4419
     4420        OUTPUT:
     4421
     4422        A subgraph of ``self`` corresponding to a (directed if
     4423        ``self`` is directed) longest path. If ``weighted == True``, a
     4424        pair ``weight, path`` is returned.
     4425
     4426        ALGORITHM:
     4427
     4428        Mixed Integer Linear Programming. (This problem is known to be
     4429        NP-Hard).
     4430
     4431        EXAMPLES:
     4432
     4433        Petersen's graph being hypohamiltonian, it has a longest path
     4434        of length `n-2`::
     4435
     4436            sage: g = graphs.PetersenGraph()
     4437            sage: lp = g.longest_path()
     4438            sage: lp.order() >= g.order() - 2
     4439            True
     4440
     4441        The heuristic totally agrees::
     4442
     4443            sage: g = graphs.PetersenGraph()
     4444            sage: g.longest_path(algorithm = "backtrack").edges()
     4445            [(0, 1, None), (1, 2, None), (2, 3, None), (3, 4, None), (4, 9, None), (5, 7, None), (5, 8, None), (6, 8, None), (6, 9, None)]
     4446
     4447        Let us compute longest path on random graphs with random
     4448        weights. We each time ensure the given graph is indeed a
     4449        path::
     4450
     4451            sage: for i in xrange(20):
     4452            ...       g = graphs.RandomGNP(15,.3)
     4453            ...       for u,v in g.edges(labels = False):
     4454            ...          g.set_edge_label(u,v,random())
     4455            ...       lp = g.longest_path()
     4456            ...       if (not lp.is_forest() or
     4457            ...           not max(lp.degree()) <= 2 or
     4458            ...           not lp.is_connected()):
     4459            ...           print "Error !"
     4460            ...           break
     4461            ...
     4462            sage: print "Test finished !"
     4463            Test finished !
     4464
     4465        TESTS::
     4466
     4467        Disconnected graphs not weighted::
     4468       
     4469            sage: g1 = graphs.PetersenGraph()
     4470            sage: g2 = 2 * g1
     4471            sage: lp1 = g1.longest_path()
     4472            sage: lp2 = g2.longest_path()
     4473            sage: len(lp1) == len(lp2)
     4474            True
     4475
     4476        Disconnected graphs weighted::
     4477
     4478            sage: g1 = graphs.PetersenGraph()
     4479            sage: for u,v in g.edges(labels = False):
     4480            ...       g.set_edge_label(u,v,random())
     4481            sage: g2 = 2 * g1
     4482            sage: lp1 = g1.longest_path(weighted = True)
     4483            sage: lp2 = g2.longest_path(weighted = True)
     4484            sage: lp1[0] == lp2[0]
     4485            True
     4486
     4487        Empty graphs::
     4488
     4489            sage: Graph().longest_path()
     4490            Graph on 0 vertices
     4491            sage: Graph().longest_path(weighted = True)
     4492            [0, Graph on 0 vertices]
     4493
     4494        Random test for digraphs::
     4495
     4496            sage: for i in xrange(20):
     4497            ...       g = digraphs.RandomDirectedGNP(15,.3)
     4498            ...       for u,v in g.edges(labels = False):
     4499            ...          g.set_edge_label(u,v,random())
     4500            ...       lp = g.longest_path()
     4501            ...       if (not lp.is_directed_acyclic() or
     4502            ...           not max(lp.out_degree()) <= 1 or
     4503            ...           not max(lp.in_degree()) <= 1 or
     4504            ...           not lp.is_connected()):
     4505            ...           print "Error !"
     4506            ...           break
     4507            ...
     4508            sage: print "Test finished !"
     4509            Test finished !
     4510        """
     4511        if weighted:
     4512            algorithm = "MILP"
     4513
     4514
     4515        # Quick improvement
     4516        if not self.is_connected():
     4517            if weighted:
     4518                return max( g.longest_path(s = s, t = t,
     4519                                           weighted = weighted,
     4520                                           algorithm = algorithm)
     4521                            for g in self.connected_components_subgraphs() )
     4522            else:
     4523                return max( (g.longest_path(s = s, t = t,
     4524                                            weighted = weighted,
     4525                                            algorithm = algorithm)
     4526                            for g in self.connected_components_subgraphs() ),
     4527                            key = lambda x : x.order() )
     4528
     4529        # Stupid cases
     4530
     4531        # - Graph having less than 1 vertex
     4532        #
     4533        # - The source has outdegree 0 in a directed graph, or
     4534        #   degree 0, or is not a vertex of the graph
     4535        #
     4536        # - The destination has indegree 0 in a directed graph, or
     4537        #   degree 0, or is not a vertex of the graph
     4538        #
     4539        # - Both s and t are specified, but there is no path between
     4540        #   the two in a directed graph (the graph is connected)
     4541
     4542        if (self.order() <= 1 or
     4543            (s != None and (
     4544                    (s not in self) or
     4545                    (self._directed and self.degree_out(s) == 0) or
     4546                    (not self._directed and self.degree(s) == 0))) or
     4547
     4548            (t != None and (
     4549                    (t not in self) or
     4550                    (self._directed and self.degree_in(t) == 0) or
     4551                    (not self._directed and self.degree(t) == 0))) or
     4552
     4553            (self._directed and s != None and t != None and
     4554             len(self.shortest_path(s,t) == 0))):
     4555
     4556            if self._directed:
     4557                from sage.graphs.all import DiGraph
     4558                return [0, DiGraph()] if weighted else DiGraph()
     4559            else:
     4560                from sage.graphs.all import Graph
     4561                return [0, Graph()] if weighted else Graph()
     4562
     4563
     4564        # Calling the backtrack heuristic if asked
     4565        if algorithm == "backtrack":
     4566            from sage.graphs.generic_graph_pyx import find_hamiltonian as fh
     4567            x = fh( self, find_path = True )[1]
     4568            return self.subgraph(vertices = x,
     4569                                 edges = zip(x[:-1], x[1:]))
     4570
     4571        ##################
     4572        # LP Formulation #
     4573        ##################
     4574
     4575        # Epsilon... Must be less than 1/(n+1), but we want to avoid
     4576        # numerical problems...
     4577        epsilon = 1/(6*float(self.order()))
     4578
     4579        # Associating a weight to a label
     4580        if weighted:
     4581            weight = lambda x : x if (x is not None and x != {}) else 1
     4582        else:
     4583            weight = lambda x : 1
     4584
     4585        from sage.numerical.mip import MixedIntegerLinearProgram, Sum
     4586
     4587        p = MixedIntegerLinearProgram()
     4588       
     4589        # edge_used[(u,v)] == 1 if (u,v) is used
     4590        edge_used = p.new_variable(binary = True)
     4591
     4592        # relaxed version of the previous variable, to prevent the
     4593        # creation of cycles
     4594        r_edge_used = p.new_variable()
     4595
     4596        # vertex_used[v] == 1 if vertex v is used
     4597        vertex_used = p.new_variable(binary = True)
     4598
     4599        if self._directed:
     4600
     4601            # if edge uv is used, vu can not be
     4602            for u,v in self.edges(labels = False):
     4603                if self.has_edge(v,u):
     4604                    p.add_constraint(edge_used[(u,v)] + edge_used[(v,u)], max = 1)
     4605           
     4606           
     4607            # A vertex is used if one of its incident edges is
     4608            for v in self:
     4609                for e in self.incoming_edges(labels = False):
     4610                    p.add_constraint(vertex_used[v] - edge_used[e], min = 0)
     4611                for e in self.outgoing_edges(labels = False):
     4612                    p.add_constraint(vertex_used[v] - edge_used[e], min = 0)
     4613
     4614            # A path is a tree. If n vertices are used, at most n-1 edges are
     4615            p.add_constraint(Sum( vertex_used[v] for v in self)
     4616                             - Sum( edge_used[e] for e in self.edges(labels = False)),
     4617                             min = 1, max = 1)
     4618
     4619            # A vertex has at most one incoming used edge and at most
     4620            # one outgoing used edge
     4621            for v in self:
     4622                p.add_constraint( Sum( edge_used[(u,v)] for u in self.neighbors_in(v) ), max = 1)
     4623                p.add_constraint( Sum( edge_used[(v,u)] for u in self.neighbors_out(v) ), max = 1)
     4624
     4625            # r_edge_used is "more" than edge_used, though it ignores
     4626            # the direction
     4627            for u,v in self.edges(labels = False):
     4628                p.add_constraint(  r_edge_used[(u,v)]
     4629                                 + r_edge_used[(v,u)]
     4630                                 - edge_used[(u,v)],
     4631                                   min = 0)
     4632
     4633            # No cycles
     4634            for v in self:
     4635                p.add_constraint( Sum( r_edge_used[(u,v)] for u in self.neighbors(v)), max = 1-epsilon)
     4636
     4637            # Enforcing the source if asked.. If s is set, it has no
     4638            # incoming edge and exactly one son
     4639            if s!=None:
     4640                p.add_constraint( Sum( edge_used[(u,s)] for u in self.neighbors_in(s)), max = 0, min = 0)
     4641                p.add_constraint( Sum( edge_used[(s,u)] for u in self.neighbors_out(s)), min = 1, max = 1)
     4642
     4643            # Enforcing the destination if asked.. If t is set, it has
     4644            # no outgoing edge and exactly one parent
     4645            if t!=None:
     4646                p.add_constraint( Sum( edge_used[(u,t)] for u in self.neighbors_in(t)), min = 1, max = 1)
     4647                p.add_constraint( Sum( edge_used[(t,u)] for u in self.neighbors_out(t)), max = 0, min = 0)
     4648
     4649               
     4650            # Defining the objective
     4651            p.set_objective( Sum( weight(l) * edge_used[(u,v)] for u,v,l in self.edges() ) )
     4652           
     4653        else:
     4654
     4655            # f_edge_used calls edge_used through reordering u and v
     4656            # to avoid having two different variables
     4657            f_edge_used = lambda u,v : edge_used[ (u,v) if hash(u) < hash(v) else (v,u) ]
     4658
     4659            # A vertex is used if one of its incident edges is
     4660            for v in self:
     4661                for u in self.neighbors(v):
     4662                    p.add_constraint(vertex_used[v] - f_edge_used(u,v), min = 0)
     4663
     4664            # A path is a tree. If n vertices are used, at most n-1 edges are
     4665            p.add_constraint( Sum( vertex_used[v] for v in self)
     4666                              - Sum( f_edge_used(u,v) for u,v in self.edges(labels = False)),
     4667                              min = 1, max = 1)
     4668
     4669            # A vertex has at most two incident edges used
     4670            for v in self:
     4671                p.add_constraint( Sum( f_edge_used(u,v) for u in self.neighbors(v) ), max = 2)
     4672
     4673            # r_edge_used is "more" than edge_used
     4674            for u,v in self.edges(labels = False):
     4675                p.add_constraint(  r_edge_used[(u,v)]
     4676                                 + r_edge_used[(v,u)]
     4677                                 - f_edge_used(u,v),
     4678                                   min = 0)
     4679
     4680            # No cycles
     4681            for v in self:
     4682                p.add_constraint( Sum( r_edge_used[(u,v)] for u in self.neighbors(v)), max = 1-epsilon)
     4683
     4684            # Enforcing the destination if asked.. If s or t are set,
     4685            # they have exactly one incident edge
     4686            if s != None:
     4687                p.add_constraint( Sum( f_edge_used(s,u) for u in self.neighbors(s) ), max = 1, min = 1)
     4688            if t != None:
     4689                p.add_constraint( Sum( f_edge_used(t,u) for u in self.neighbors(t) ), max = 1, min = 1)
     4690
     4691            # Defining the objective
     4692            p.set_objective( Sum( weight(l) * f_edge_used(u,v) for u,v,l in self.edges() ) )
     4693
     4694        # Computing the result. No exception has to be raised, as this
     4695        # problem always has a solution (there is at least one edge,
     4696        # and a path from s to t if they are specified
     4697
     4698        p.solve(solver = solver, log = verbose)
     4699
     4700        edge_used = p.get_values(edge_used)
     4701        vertex_used = p.get_values(vertex_used)
     4702
     4703        if self._directed:
     4704            g = self.subgraph( vertices =
     4705                                  (v for v in self if vertex_used[v] >= .5),
     4706                               edges =
     4707                                  ( (u,v,l) for u,v,l in self.edges()
     4708                                    if edge_used[(u,v)] >= .5 ))
     4709        else:
     4710            g = self.subgraph( vertices =
     4711                                  (v for v in self if vertex_used[v] >= .5),
     4712                               edges =
     4713                                  ( (u,v,l) for u,v,l in self.edges()
     4714                                    if f_edge_used(u,v) >= .5 ))
     4715
     4716        if weighted:
     4717            return sum(map(weight,g.edge_labels())), g
     4718        else:
     4719            return g
     4720
     4721
    43774722    def traveling_salesman_problem(self, weighted = True, solver = None, verbose = 0):
    43784723        r"""
    43794724        Solves the traveling salesman problem (TSP)