Ticket #9910: trac_9910.2.patch

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