Ticket #11429: trac_11429_native_enumeration_of_lattice_polytope_points.patch
| File trac_11429_native_enumeration_of_lattice_polytope_points.patch, 16.0 KB (added by , 12 years ago) |
|---|
-
sage/geometry/polyhedra.py
# HG changeset patch # User Volker Braun <vbraun@stp.dias.ie> # Date 1307937431 25200 # Node ID ab54908876df0724035fdd72e0a31ac29de52331 # Parent 5dbd209f6246214ffb6e1603b90df7d87be2c5bc Trac #11429: Count integral points without PALP We want our own code to enumerate lattice points in polyhedra because: * Going through the PALP pexpect interface is annoyingly slow. * no more compile-time bounds * It seems like PALP uses a very unsophisticated algorithm diff --git a/sage/geometry/polyhedra.py b/sage/geometry/polyhedra.py
a b 1234 1234 """ 1235 1235 return Hobj.A() * self._representation_vector + Hobj.b() 1236 1236 1237 def is_integral(self): 1238 r""" 1239 Return whether the coordinates of the vertex are all integral. 1240 1241 OUTPUT: 1242 1243 Boolean. 1244 1245 EXAMPLES:: 1246 1247 sage: p = Polyhedron([(1/2,3,5), (0,0,0), (2,3,7)]) 1248 sage: [ v.is_integral() for v in p.vertex_generator() ] 1249 [False, True, True] 1250 """ 1251 return all(x in ZZ for x in self._representation_vector) 1237 1252 1238 1253 1239 1254 class Ray(Vrepresentation): … … 3402 3417 3403 3418 return True 3404 3419 3405 3406 3420 def gale_transform(self): 3407 3421 """ 3408 3422 Return the Gale transform of a polytope as described in the … … 3427 3441 3428 3442 Lectures in Geometric Combinatorics, R.R.Thomas, 2006, AMS Press. 3429 3443 """ 3430 if not self.is_compact(): raise ValueError , "Not a polytope"3444 if not self.is_compact(): raise ValueError('Not a polytope.') 3431 3445 3432 3446 A = matrix(self.n_vertices(), 3433 3447 [ [1]+list(x) for x in self.vertex_generator()]) … … 3435 3449 A_ker = A.right_kernel() 3436 3450 return A_ker.basis_matrix().transpose().rows() 3437 3451 3452 def triangulate(self, engine='auto', connected=True, fine=False, regular=None, star=None): 3453 r""" 3454 Returns a triangulation of the polytope. 3455 3456 INPUT: 3457 3458 - ``engine`` -- either 'auto' (default), 'internal', or 3459 'TOPCOM'. The latter two instruct this package to always 3460 use its own triangulation algorithms or TOPCOM's algorithms, 3461 respectively. By default ('auto'), TOPCOM is used if it is 3462 available and internal routines otherwise. 3463 3464 The remaining keyword parameters are passed through to the 3465 :class:`~sage.geometry.triangulation.point_configuration.PointConfiguration` 3466 constructor: 3467 3468 - ``connected`` -- boolean (default: ``True``). Whether the 3469 triangulations should be connected to the regular 3470 triangulations via bistellar flips. These are much easier to 3471 compute than all triangulations. 3472 3473 - ``fine`` -- boolean (default: ``False``). Whether the 3474 triangulations must be fine, that is, make use of all points 3475 of the configuration. 3476 3477 - ``regular`` -- boolean or ``None`` (default: 3478 ``None``). Whether the triangulations must be regular. A 3479 regular triangulation is one that is induced by a 3480 piecewise-linear convex support function. In other words, 3481 the shadows of the faces of a polyhedron in one higher 3482 dimension. 3483 3484 * ``True``: Only regular triangulations. 3485 3486 * ``False``: Only non-regular triangulations. 3487 3488 * ``None`` (default): Both kinds of triangulation. 3489 3490 - ``star`` -- either ``None`` (default) or a point. Whether 3491 the triangulations must be star. A triangulation is star if 3492 all maximal simplices contain a common point. The central 3493 point can be specified by its index (an integer) in the 3494 given points or by its coordinates (anything iterable.) 3495 3496 OUTPUT: 3497 3498 A triangulation of the convex hull of the vertices as a 3499 :class:`~sage.geometry.triangulation.point_configuration.Triangulation`. The 3500 indices in the triangulation correspond to the 3501 :meth:`Vrepresentation` objects. 3502 3503 EXAMPLES:: 3504 3505 sage: cube = polytopes.n_cube(3) 3506 sage: triangulation = cube.triangulate(engine='internal') # to make doctest independent of TOPCOM 3507 sage: triangulation 3508 (<0,1,2,4>, <1,2,3,4>, <1,3,4,5>, <2,3,4,6>, <3,4,5,6>, <3,5,6,7>) 3509 sage: simplex_indices = triangulation[0]; simplex_indices 3510 (0, 1, 2, 4) 3511 sage: simplex_vertices = [ cube.Vrepresentation(i) for i in simplex_indices ] 3512 sage: simplex_vertices 3513 [A vertex at (1, 1, 1), A vertex at (-1, 1, 1), 3514 A vertex at (1, -1, 1), A vertex at (1, 1, -1)] 3515 sage: Polyhedron(simplex_vertices) 3516 A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 4 vertices. 3517 """ 3518 if not self.is_compact(): 3519 raise NotImplementedError('I can only triangulate compact polytopes.') 3520 from sage.geometry.triangulation.point_configuration import PointConfiguration 3521 pc = PointConfiguration((v.vector() for v in self.vertex_generator()), 3522 connected=connected, fine=fine, regular=regular, star=star) 3523 pc.set_engine(engine) 3524 return pc.triangulate() 3438 3525 3439 3526 def triangulated_facial_incidences(self): 3440 3527 """ … … 3454 3541 If the figure is already composed of triangles, then all is well:: 3455 3542 3456 3543 sage: Polyhedron(vertices = [[5,0,0],[0,5,0],[5,5,0],[2,2,5]]).triangulated_facial_incidences() 3544 doctest:...: DeprecationWarning: (Since Sage Version 4.7.1) 3545 This method is deprecated. Use triangulate() instead. 3457 3546 [[0, [0, 2, 3]], [1, [0, 1, 2]], [2, [1, 2, 3]], [3, [0, 1, 3]]] 3458 3547 3459 3548 Otherwise some faces get split up to triangles:: … … 3462 3551 sage: Polyhedron(vertices = [[2,0,0],[4,1,0],[0,5,0],[5,5,0],[1,1,0],[0,0,1]]).triangulated_facial_incidences() 3463 3552 [[0, [0, 1, 5]], [1, [0, 4, 5]], [2, [2, 4, 5]], [3, [2, 3, 5]], [4, [1, 3, 5]], [5, [0, 1, 4]], [5, [1, 4, 3]], [5, [4, 3, 2]]] 3464 3553 """ 3554 from sage.misc.misc import deprecation 3555 deprecation('This method is deprecated. Use triangulate() instead.', 'Sage Version 4.7.1') 3465 3556 try: 3466 3557 return self._triangulated_facial_incidences 3467 3558 except AttributeError: … … 3528 3619 3529 3620 sage: p = polytopes.cuboctahedron() 3530 3621 sage: sc = p.simplicial_complex() 3622 doctest:...: DeprecationWarning: (Since Sage Version 4.7.1) 3623 This method is deprecated. Use triangulate().simplicial_complex() instead. 3624 doctest:...: DeprecationWarning: (Since Sage Version 4.7.1) 3625 This method is deprecated. Use triangulate() instead. 3531 3626 sage: sc 3532 3627 Simplicial complex with 13 vertices and 20 facets 3533 3628 """ 3629 from sage.misc.misc import deprecation 3630 deprecation('This method is deprecated. Use triangulate().simplicial_complex() instead.', 3631 'Sage Version 4.7.1') 3534 3632 from sage.homology.simplicial_complex import SimplicialComplex 3535 3633 return SimplicialComplex(vertex_set = self.n_vertices(), 3536 3634 maximal_faces = [x[1] for x in self.triangulated_facial_incidences()]) … … 4387 4485 4388 4486 return True 4389 4487 4488 def is_simplex(self): 4489 r""" 4490 Return whether the polyhedron is a simplex. 4491 4492 EXAMPLES:: 4493 4494 sage: Polyhedron([(0,0,0), (1,0,0), (0,1,0)]).is_simplex() 4495 True 4496 sage: polytopes.n_simplex(3).is_simplex() 4497 True 4498 sage: polytopes.n_cube(3).is_simplex() 4499 False 4500 """ 4501 return self.is_compact() and (self.dim()+1==self.n_vertices()) 4502 4390 4503 def is_lattice_polytope(self): 4391 4504 r""" 4392 4505 Return whether the polyhedron is a lattice polytope. … … 4407 4520 return self._is_lattice_polytope 4408 4521 except AttributeError: 4409 4522 pass 4410 4411 self._is_lattice_polytope = False 4412 if self.is_compact(): 4413 try: 4414 matrix(ZZ, self.vertices()) 4415 self._is_lattice_polytope = True 4416 except TypeError: 4417 pass 4418 4523 self._is_lattice_polytope = self.is_compact() and \ 4524 all(v.is_integral() for v in self.vertex_generator()) 4419 4525 return self._is_lattice_polytope 4420 4526 4421 4527 def lattice_polytope(self, envelope=False): … … 4518 4624 self._lattice_polytope = LatticePolytope(vertices) 4519 4625 return self._lattice_polytope 4520 4626 4521 def integral_points(self):4627 def _integral_points_PALP(self): 4522 4628 r""" 4523 Return the integral points in the polyhedron. 4629 Return the integral points in the polyhedron using PALP. 4630 4631 This method is for testing purposes and will eventually be removed. 4524 4632 4525 4633 OUTPUT: 4526 4634 … … 4529 4637 4530 4638 EXAMPLES:: 4531 4639 4532 sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]). integral_points()4640 sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)])._integral_points_PALP() 4533 4641 [(-1, -1), (1, 0), (1, 1), (0, 1), (0, 0)] 4534 4642 sage: Polyhedron(vertices=[(-1/2,-1/2),(1,0),(1,1),(0,1)]).lattice_polytope(True).points() 4535 4643 [ 0 -1 -1 1 1 0 0] 4536 4644 [-1 0 -1 0 1 1 0] 4537 sage: Polyhedron(vertices=[(-1/2,-1/2),(1,0),(1,1),(0,1)]). integral_points()4645 sage: Polyhedron(vertices=[(-1/2,-1/2),(1,0),(1,1),(0,1)])._integral_points_PALP() 4538 4646 [(1, 0), (1, 1), (0, 1), (0, 0)] 4539 4647 """ 4540 4648 if not self.is_compact(): 4541 4649 raise ValueError, 'Can only enumerate points in a compact polyhedron.' 4542 4650 4543 4651 lp = self.lattice_polytope(True) 4652 4653 # remove cached values to get accurate timings 4654 try: 4655 del lp._points 4656 del lp._npoints 4657 except AttributeError: 4658 pass 4544 4659 4545 4660 if self.is_lattice_polytope(): 4546 4661 return lp.points().columns() … … 4549 4664 lp.points().columns()) 4550 4665 return points 4551 4666 4667 def _integral_points_lattice_simplex(self, vertices=None): 4668 r""" 4669 Return the integral points in a lattice simplex. 4670 4671 INPUT: 4672 4673 - ``vertices`` -- ``None`` (default) or a list of 4674 integers. The indices of vertices that span the simplex 4675 under consideration. By default, it is assumed that the 4676 polytope itself is a simplex and all vertices are used. 4677 4678 OUTPUT: 4679 4680 A tuple containing the integral point coordinates as `\ZZ`-vectors. 4681 4682 EXAMPLES:: 4683 4684 sage: simplex = Polyhedron([(1,2,3), (2,3,7), (-2,-3,-11)]) 4685 sage: simplex._integral_points_lattice_simplex() 4686 ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7)) 4687 4688 The simplex need not be full-dimensional:: 4689 4690 sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)]) 4691 sage: simplex._integral_points_lattice_simplex() 4692 ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5)) 4693 4694 sage: point = Polyhedron([(2,3,7)]) 4695 sage: point._integral_points_lattice_simplex() 4696 ((2, 3, 7),) 4697 4698 TESTS:: 4699 4700 sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)] 4701 sage: simplex = Polyhedron(v); simplex 4702 A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices. 4703 sage: len(simplex._integral_points_lattice_simplex()) 4704 49 4705 4706 sage: v = [(4,-1,-1,-1), (-1,4,-1,-1), (-1,-1,4,-1), (-1,-1,-1,4), (-1,-1,-1,-1)] 4707 sage: P4mirror = Polyhedron(v) 4708 sage: len( P4mirror._integral_points_lattice_simplex() ) 4709 126 4710 """ 4711 if not vertices: 4712 assert self.is_simplex() and self.is_lattice_polytope() 4713 vertices = [ vector(ZZ, v) for v in self.Vrepresentation() ] 4714 else: 4715 vertices = [ vector(ZZ, self.Vrepresentation(i)) for i in vertices ] 4716 4717 origin = vertices.pop() 4718 origin.set_immutable() 4719 if len(vertices)==0: 4720 return tuple([origin]) 4721 rays = [ v-origin for v in vertices ] 4722 4723 # Find equation Ax<=b that cuts out simplex from parallelotope 4724 R = matrix(rays) 4725 if R.is_square(): 4726 b = 1 4727 A = R.solve_right(vector([1]*len(rays))) 4728 else: 4729 b = 1 4730 RRt = R * R.transpose() 4731 A = RRt.solve_right(vector([1]*len(rays))) * R 4732 4733 from sage.geometry.cone import parallelotope_points 4734 gens = list(parallelotope_points(rays)) + list(rays) 4735 gens = [ origin+x for x in 4736 list(parallelotope_points(rays)) + list(rays) 4737 if A*x<=b ] 4738 for x in gens: 4739 x.set_immutable() 4740 return tuple(gens) 4741 4742 def integral_points(self): 4743 r""" 4744 Return the integral points in the polyhedron. 4745 4746 OUTPUT: 4747 4748 The list of integral points in the polyhedron. If the 4749 polyhedron is not compact, a ``ValueError`` is raised. 4750 4751 EXAMPLES:: 4752 4753 sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]).integral_points() 4754 ((0, 1), (1, 0), (0, 0), (1, 1), (-1, -1)) 4755 4756 Finally, 3-d reflexive polytope number 4078:: 4757 4758 sage: v = [(1,0,0), (0,1,0), (0,0,1), (0,0,-1), (0,-2,1), 4759 ... (-1,2,-1), (-1,2,-2), (-1,1,-2), (-1,-1,2), (-1,-3,2)] 4760 sage: pts1 = set(Polyhedron(v).integral_points()) # Sage's own code 4761 sage: pts2 = LatticePolytope(v).points().columns() # PALP 4762 sage: for p in pts2: 4763 ... p.set_immutable() 4764 sage: pts2 = set(pts2) 4765 sage: pts1 == pts2 4766 True 4767 4768 TODO: profile and make faster:: 4769 4770 sage: timeit('Polyhedron(v).integral_points()') # random output 4771 sage: timeit('LatticePolytope(v).points()') # random output 4772 """ 4773 if not self.is_compact(): 4774 raise ValueError('Can only enumerate points in a compact polyhedron.') 4775 if not self.is_lattice_polytope(): 4776 raise NotImplementedError('Only implemented for polyhedra with integral vertices') 4777 4778 if self.is_simplex(): 4779 points = self._integral_points_lattice_simplex() 4780 else: 4781 triangulation = self.triangulate() 4782 points = set() 4783 for simplex in triangulation: 4784 points.update(self._integral_points_lattice_simplex(simplex)) 4785 points = tuple(points) 4786 4787 # assert all(self.contains(p) for p in points) # slow 4788 return points 4552 4789 4553 4790 def combinatorial_automorphism_group(self): 4554 4791 """ -
sage/geometry/triangulation/point_configuration.py
diff --git a/sage/geometry/triangulation/point_configuration.py b/sage/geometry/triangulation/point_configuration.py
a b 690 690 return Fan(self, (vector(R, p) - origin for p in points)) 691 691 692 692 693 def simplicial_complex(self): 694 r""" 695 Return a simplicial complex from a triangulation of the point 696 configuration. 697 698 OUTPUT: 699 700 A :class:`~sage.homology.simplicial_complex.SimplicialComplex`. 701 702 EXAMPLES:: 703 704 sage: p = polytopes.cuboctahedron() 705 sage: sc = p.triangulate().simplicial_complex() 706 sage: sc 707 Simplicial complex with 12 vertices and 16 facets 708 """ 709 from sage.homology.simplicial_complex import SimplicialComplex 710 from sage.misc.all import flatten 711 vertex_set = set(flatten(self)) 712 return SimplicialComplex(vertex_set = vertex_set, 713 maximal_faces = self) 714 693 715 694 716 ######################################################################## 695 717 class PointConfiguration(UniqueRepresentation, PointConfiguration_base):
