| 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 | |