Ticket #9910: trac-9910_reviewer.patch
File trac-9910_reviewer.patch, 21.4 KB (added by , 10 years ago) |
---|
-
sage/graphs/generic_graph.py
# HG changeset patch # User Minh Van Nguyen <nguyenminh2@gmail.com> # Date 1290885462 28800 # Node ID 9a5f569d00d280dfd09ba4f2417ab8745d5c6b3c # Parent 9c9495ebab27b9b708c83c50bf7b3d88c0a903db #9910: reviewer patch: longest path Changes include: * Many typo fixes. * Proper ReST formatting. * Use `range()` instead of `xrange()` to conform to Python 3.x. * Use the function `print()` instead of the statement `print`. * Make the doctests and code conform to Python coding conventions. * When comparing an object with `None`, use the idiom "is not" instead of "!=". * Doctests to take care of the case of the trivial graph. * Remove superfluous trailing white spaces. These are mostly harmless, but in some cases they can result in hard to detect bugs. * Handle the case where the argument `algorithm` is not assigned either "backtrack" or "MILP". * Some general clean-ups. diff --git a/sage/graphs/generic_graph.py b/sage/graphs/generic_graph.py
a b 4374 4374 4375 4375 return val 4376 4376 4377 def longest_path(self, s = None, t = None, weighted = False, algorithm = "MILP", solver = None, verbose =0):4377 def longest_path(self, s=None, t=None, weighted=False, algorithm="MILP", solver=None, verbose=0): 4378 4378 r""" 4379 4379 Returns a longest path of ``self``. 4380 4380 … … 4386 4386 - ``t`` (vertex) -- forces the destination of the path. Set to 4387 4387 ``None`` by default. 4388 4388 4389 - ``weighted`` (boolean) -- whether the label on the edges are4390 to be considered as weights (a label set to ``None`` or4389 - ``weighted`` (boolean) -- whether the labels on the edges are 4390 to be considered as weights (a label set to ``None`` or 4391 4391 ``{}`` being considered as a weight of `1`). Set to 4392 4392 ``False`` by default. 4393 4393 4394 4394 - ``algorithm`` -- one of ``"MILP"`` (default) or 4395 ``"backtrack"``. Two remarks on this respect 4395 ``"backtrack"``. Two remarks on this respect: 4396 4396 4397 4397 * While the MILP formulation returns an exact answer, 4398 4398 the backtrack algorithm is a randomized heuristic. 4399 4399 4400 4400 * As the backtrack algorithm does not support edge 4401 weighting, setting ``weighted `` to ``True`` will force4401 weighting, setting ``weighted=True`` will force 4402 4402 the use of the MILP algorithm. 4403 4403 4404 - ``solver`` -- (default: ``None``) Specify aLinear Program (LP)4404 - ``solver`` -- (default: ``None``) Specify the Linear Program (LP) 4405 4405 solver to be used. If set to ``None``, the default one is used. For 4406 4406 more information on LP solvers and which default solver is used, see 4407 4407 the method … … 4441 4441 The heuristic totally agrees:: 4442 4442 4443 4443 sage: g = graphs.PetersenGraph() 4444 sage: g.longest_path(algorithm ="backtrack").edges()4444 sage: g.longest_path(algorithm="backtrack").edges() 4445 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 4446 4447 Let us compute longest path on random graphs with random4448 weights. We each time ensure the givengraph is indeed a4447 Let us compute longest paths on random graphs with random 4448 weights. Each time, we ensure the resulting graph is indeed a 4449 4449 path:: 4450 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())4451 sage: for i in range(20): 4452 ... g = graphs.RandomGNP(15, 0.3) 4453 ... for u, v in g.edges(labels=False): 4454 ... g.set_edge_label(u, v, random()) 4455 4455 ... lp = g.longest_path() 4456 4456 ... if (not lp.is_forest() or 4457 4457 ... not max(lp.degree()) <= 2 or 4458 4458 ... not lp.is_connected()): 4459 ... print "Error !"4459 ... print("Error!") 4460 4460 ... break 4461 ... 4462 sage: print "Test finished !" 4463 Test finished ! 4464 4465 TESTS:: 4461 4462 TESTS: 4463 4464 The argument ``algorithm`` must be either ``'backtrack'`` or 4465 ``'MILP'``:: 4466 4467 sage: graphs.PetersenGraph().longest_path(algorithm="abc") 4468 Traceback (most recent call last): 4469 ... 4470 ValueError: algorithm must be either 'backtrack' or 'MILP' 4466 4471 4467 4472 Disconnected graphs not weighted:: 4468 4473 … … 4476 4481 Disconnected graphs weighted:: 4477 4482 4478 4483 sage: g1 = graphs.PetersenGraph() 4479 sage: for u,v in g.edges(labels =False):4480 ... g.set_edge_label(u, v,random())4484 sage: for u,v in g.edges(labels=False): 4485 ... g.set_edge_label(u, v, random()) 4481 4486 sage: g2 = 2 * g1 4482 sage: lp1 = g1.longest_path(weighted =True)4483 sage: lp2 = g2.longest_path(weighted =True)4487 sage: lp1 = g1.longest_path(weighted=True) 4488 sage: lp2 = g2.longest_path(weighted=True) 4484 4489 sage: lp1[0] == lp2[0] 4485 4490 True 4486 4491 … … 4488 4493 4489 4494 sage: Graph().longest_path() 4490 4495 Graph on 0 vertices 4491 sage: Graph().longest_path(weighted =True)4496 sage: Graph().longest_path(weighted=True) 4492 4497 [0, Graph on 0 vertices] 4498 sage: graphs.EmptyGraph().longest_path() 4499 Graph on 0 vertices 4500 sage: graphs.EmptyGraph().longest_path(weighted=True) 4501 [0, Graph on 0 vertices] 4502 4503 Trivial graphs:: 4504 4505 sage: G = Graph() 4506 sage: G.add_vertex(0) 4507 sage: G.longest_path() 4508 Graph on 0 vertices 4509 sage: G.longest_path(weighted=True) 4510 [0, Graph on 0 vertices] 4511 sage: graphs.CompleteGraph(1).longest_path() 4512 Graph on 0 vertices 4513 sage: graphs.CompleteGraph(1).longest_path(weighted=True) 4514 [0, Graph on 0 vertices] 4493 4515 4494 4516 Random test for digraphs:: 4495 4517 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())4518 sage: for i in range(20): 4519 ... g = digraphs.RandomDirectedGNP(15, 0.3) 4520 ... for u, v in g.edges(labels=False): 4521 ... g.set_edge_label(u, v, random()) 4500 4522 ... lp = g.longest_path() 4501 4523 ... if (not lp.is_directed_acyclic() or 4502 4524 ... not max(lp.out_degree()) <= 1 or 4503 4525 ... not max(lp.in_degree()) <= 1 or 4504 4526 ... not lp.is_connected()): 4505 ... print "Error !"4527 ... print("Error!") 4506 4528 ... break 4507 ...4508 sage: print "Test finished !"4509 Test finished !4510 4529 """ 4511 4530 if weighted: 4512 4531 algorithm = "MILP" 4513 4532 if algorithm not in ("backtrack", "MILP"): 4533 raise ValueError("algorithm must be either 'backtrack' or 'MILP'") 4514 4534 4515 4535 # Quick improvement 4516 4536 if not self.is_connected(): 4517 4537 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())4538 return max(g.longest_path(s=s, t=t, 4539 weighted=weighted, 4540 algorithm=algorithm) 4541 for g in self.connected_components_subgraphs()) 4542 else: 4543 return max((g.longest_path(s=s, t=t, 4544 weighted=weighted, 4545 algorithm=algorithm) 4546 for g in self.connected_components_subgraphs()), 4547 key=lambda x: x.order()) 4528 4548 4529 4549 # Stupid cases 4530 4531 # - Graph having less than 1 vertex 4550 # - Graph having <= 1 vertex. 4532 4551 # 4533 4552 # - The source has outdegree 0 in a directed graph, or 4534 # degree 0, or is not a vertex of the graph 4553 # degree 0, or is not a vertex of the graph. 4535 4554 # 4536 4555 # - The destination has indegree 0 in a directed graph, or 4537 # degree 0, or is not a vertex of the graph 4556 # degree 0, or is not a vertex of the graph. 4538 4557 # 4539 4558 # - Both s and t are specified, but there is no path between 4540 # the two in a directed graph (the graph is connected) 4541 4559 # the two in a directed graph (the graph is connected). 4542 4560 if (self.order() <= 1 or 4543 (s !=None and (4561 (s is not None and ( 4544 4562 (s not in self) or 4545 4563 (self._directed and self.degree_out(s) == 0) or 4546 4564 (not self._directed and self.degree(s) == 0))) or 4547 4548 (t != None and ( 4565 (t is not None and ( 4549 4566 (t not in self) or 4550 4567 (self._directed and self.degree_in(t) == 0) or 4551 4568 (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 4569 (self._directed and (s is not None) and (t is not None) and 4570 len(self.shortest_path(s, t) == 0))): 4556 4571 if self._directed: 4557 4572 from sage.graphs.all import DiGraph 4558 4573 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 4574 from sage.graphs.all import Graph 4575 return [0, Graph()] if weighted else Graph() 4563 4576 4564 4577 # Calling the backtrack heuristic if asked 4565 4578 if algorithm == "backtrack": 4566 4579 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:])) 4580 x = fh(self, find_path=True)[1] 4581 return self.subgraph(vertices=x, edges=zip(x[:-1], x[1:])) 4570 4582 4571 4583 ################## 4572 4584 # LP Formulation # … … 4578 4590 4579 4591 # Associating a weight to a label 4580 4592 if weighted: 4581 weight = lambda x 4582 else: 4583 weight = lambda x 4593 weight = lambda x: x if (x is not None and x != {}) else 1 4594 else: 4595 weight = lambda x: 1 4584 4596 4585 4597 from sage.numerical.mip import MixedIntegerLinearProgram, Sum 4586 4587 4598 p = MixedIntegerLinearProgram() 4588 4599 4589 4600 # edge_used[(u,v)] == 1 if (u,v) is used 4590 edge_used = p.new_variable(binary =True)4601 edge_used = p.new_variable(binary=True) 4591 4602 4592 4603 # relaxed version of the previous variable, to prevent the 4593 4604 # creation of cycles 4594 4605 r_edge_used = p.new_variable() 4595 4606 4596 4607 # vertex_used[v] == 1 if vertex v is used 4597 vertex_used = p.new_variable(binary = True) 4598 4599 if self._directed: 4600 4608 vertex_used = p.new_variable(binary=True) 4609 4610 if self._directed: 4601 4611 # 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 4612 for u, v in self.edges(labels=False): 4613 if self.has_edge(v, u): 4614 p.add_constraint(edge_used[(u,v)] + edge_used[(v,u)], max=1) 4607 4615 # A vertex is used if one of its incident edges is 4608 4616 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 4617 for e in self.incoming_edges(labels=False): 4618 p.add_constraint(vertex_used[v] - edge_used[e], min=0) 4619 for e in self.outgoing_edges(labels=False): 4620 p.add_constraint(vertex_used[v] - edge_used[e], min=0) 4614 4621 # 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 4622 p.add_constraint( 4623 Sum(vertex_used[v] for v in self) 4624 - Sum(edge_used[e] for e in self.edges(labels=False)), 4625 min=1, max=1) 4619 4626 # A vertex has at most one incoming used edge and at most 4620 4627 # one outgoing used edge 4621 4628 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 4629 p.add_constraint( 4630 Sum(edge_used[(u,v)] for u in self.neighbors_in(v)), 4631 max=1) 4632 p.add_constraint( 4633 Sum(edge_used[(v,u)] for u in self.neighbors_out(v)), 4634 max=1) 4625 4635 # r_edge_used is "more" than edge_used, though it ignores 4626 4636 # 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 4637 for u, v in self.edges(labels=False): 4638 p.add_constraint(r_edge_used[(u,v)] 4639 + r_edge_used[(v,u)] 4640 - edge_used[(u,v)], 4641 min=0) 4633 4642 # No cycles 4634 4643 for v in self: 4635 p.add_constraint( Sum( r_edge_used[(u,v)] for u in self.neighbors(v)), max = 1-epsilon) 4636 4644 p.add_constraint( 4645 Sum(r_edge_used[(u,v)] for u in self.neighbors(v)), 4646 max=1-epsilon) 4637 4647 # Enforcing the source if asked.. If s is set, it has no 4638 4648 # 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 4649 if s is not None: 4650 p.add_constraint( 4651 Sum(edge_used[(u,s)] for u in self.neighbors_in(s)), 4652 max=0, min=0) 4653 p.add_constraint( 4654 Sum(edge_used[(s,u)] for u in self.neighbors_out(s)), 4655 min=1, max=1) 4643 4656 # Enforcing the destination if asked.. If t is set, it has 4644 4657 # 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 4658 if t is not None: 4659 p.add_constraint( 4660 Sum(edge_used[(u,t)] for u in self.neighbors_in(t)), 4661 min=1, max=1) 4662 p.add_constraint( 4663 Sum(edge_used[(t,u)] for u in self.neighbors_out(t)), 4664 max=0, min=0) 4650 4665 # 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 4666 p.set_objective( 4667 Sum(weight(l) * edge_used[(u,v)] for u, v, l in self.edges())) 4668 else: 4655 4669 # f_edge_used calls edge_used through reordering u and v 4656 4670 # 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 4671 f_edge_used = lambda u, v: edge_used[ 4672 (u,v) if hash(u) < hash(v) else (v,u)] 4659 4673 # A vertex is used if one of its incident edges is 4660 4674 for v in self: 4661 4675 for u in self.neighbors(v): 4662 p.add_constraint(vertex_used[v] - f_edge_used(u,v), min = 0) 4663 4676 p.add_constraint(vertex_used[v] - f_edge_used(u,v), min=0) 4664 4677 # 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 4678 p.add_constraint( 4679 Sum(vertex_used[v] for v in self) 4680 - Sum(f_edge_used(u,v) for u, v in self.edges(labels=False)), 4681 min=1, max=1) 4669 4682 # A vertex has at most two incident edges used 4670 4683 for v in self: 4671 p.add_constraint( Sum( f_edge_used(u,v) for u in self.neighbors(v) ), max = 2)4672 4684 p.add_constraint( 4685 Sum(f_edge_used(u,v) for u in self.neighbors(v)), max=2) 4673 4686 # 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 4687 for u, v in self.edges(labels=False): 4688 p.add_constraint(r_edge_used[(u,v)] 4689 + r_edge_used[(v,u)] 4690 - f_edge_used(u,v), 4691 min=0) 4680 4692 # No cycles 4681 4693 for v in self: 4682 p.add_constraint( Sum( r_edge_used[(u,v)] for u in self.neighbors(v)), max = 1-epsilon) 4683 4694 p.add_constraint( 4695 Sum(r_edge_used[(u,v)] for u in self.neighbors(v)), 4696 max=1-epsilon) 4684 4697 # Enforcing the destination if asked.. If s or t are set, 4685 4698 # 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 4699 if s is not None: 4700 p.add_constraint( 4701 Sum(f_edge_used(s,u) for u in self.neighbors(s)), 4702 max=1, min=1) 4703 if t is not None: 4704 p.add_constraint( 4705 Sum(f_edge_used(t,u) for u in self.neighbors(t)), 4706 max=1, min=1) 4691 4707 # Defining the objective 4692 p.set_objective( Sum( weight(l) * f_edge_used(u,v) for u,v,l in self.edges() ))4693 4708 p.set_objective(Sum(weight(l) * f_edge_used(u,v) 4709 for u, v, l in self.edges())) 4694 4710 # Computing the result. No exception has to be raised, as this 4695 4711 # 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 4712 # and a path from s to t if they are specified). 4713 p.solve(solver=solver, log=verbose) 4700 4714 edge_used = p.get_values(edge_used) 4701 4715 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 self._directed: 4717 g = self.subgraph( 4718 vertices=(v for v in self if vertex_used[v] >= 0.5), 4719 edges=((u,v,l) for u, v, l in self.edges() 4720 if edge_used[(u,v)] >= 0.5)) 4721 else: 4722 g = self.subgraph( 4723 vertices=(v for v in self if vertex_used[v] >= 0.5), 4724 edges=((u,v,l) for u, v, l in self.edges() 4725 if f_edge_used(u,v) >= 0.5)) 4716 4726 if weighted: 4717 return sum(map(weight, g.edge_labels())), g4727 return sum(map(weight, g.edge_labels())), g 4718 4728 else: 4719 4729 return g 4720 4730 4721 4722 4731 def traveling_salesman_problem(self, weighted = True, solver = None, verbose = 0): 4723 4732 r""" 4724 4733 Solves the traveling salesman problem (TSP)