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 vbraun, 12 years ago)

Updated patch

  • 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  
    12341234        """
    12351235        return Hobj.A() * self._representation_vector + Hobj.b()
    12361236
     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)
    12371252
    12381253
    12391254class Ray(Vrepresentation):
     
    34023417
    34033418        return True
    34043419
    3405 
    34063420    def gale_transform(self):
    34073421        """
    34083422        Return the Gale transform of a polytope as described in the
     
    34273441
    34283442            Lectures in Geometric Combinatorics, R.R.Thomas, 2006, AMS Press.
    34293443        """
    3430         if not self.is_compact(): raise ValueError, "Not a polytope"
     3444        if not self.is_compact(): raise ValueError('Not a polytope.')
    34313445
    34323446        A = matrix(self.n_vertices(),
    34333447                   [ [1]+list(x) for x in self.vertex_generator()])
     
    34353449        A_ker = A.right_kernel()
    34363450        return A_ker.basis_matrix().transpose().rows()
    34373451
     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()
    34383525
    34393526    def triangulated_facial_incidences(self):
    34403527        """
     
    34543541        If the figure is already composed of triangles, then all is well::
    34553542
    34563543            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.
    34573546            [[0, [0, 2, 3]], [1, [0, 1, 2]], [2, [1, 2, 3]], [3, [0, 1, 3]]]
    34583547                             
    34593548        Otherwise some faces get split up to triangles::
     
    34623551            sage: Polyhedron(vertices = [[2,0,0],[4,1,0],[0,5,0],[5,5,0],[1,1,0],[0,0,1]]).triangulated_facial_incidences()
    34633552            [[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]]]
    34643553        """
     3554        from sage.misc.misc import deprecation
     3555        deprecation('This method is deprecated. Use triangulate() instead.', 'Sage Version 4.7.1')
    34653556        try:
    34663557            return self._triangulated_facial_incidences
    34673558        except AttributeError:
     
    35283619
    35293620            sage: p = polytopes.cuboctahedron()
    35303621            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.
    35313626            sage: sc   
    35323627            Simplicial complex with 13 vertices and 20 facets
    35333628        """
     3629        from sage.misc.misc import deprecation
     3630        deprecation('This method is deprecated. Use triangulate().simplicial_complex() instead.',
     3631                    'Sage Version 4.7.1')
    35343632        from sage.homology.simplicial_complex import SimplicialComplex
    35353633        return SimplicialComplex(vertex_set = self.n_vertices(),
    35363634                                 maximal_faces = [x[1] for x in self.triangulated_facial_incidences()])
     
    43874485
    43884486        return True
    43894487
     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
    43904503    def is_lattice_polytope(self):
    43914504        r"""
    43924505        Return whether the polyhedron is a lattice polytope.
     
    44074520            return self._is_lattice_polytope
    44084521        except AttributeError:
    44094522            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())
    44194525        return self._is_lattice_polytope
    44204526       
    44214527    def lattice_polytope(self, envelope=False):
     
    45184624        self._lattice_polytope = LatticePolytope(vertices)
    45194625        return self._lattice_polytope
    45204626
    4521     def integral_points(self):
     4627    def _integral_points_PALP(self):
    45224628        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.
    45244632
    45254633        OUTPUT:
    45264634       
     
    45294637
    45304638        EXAMPLES::
    45314639       
    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()
    45334641            [(-1, -1), (1, 0), (1, 1), (0, 1), (0, 0)]
    45344642            sage: Polyhedron(vertices=[(-1/2,-1/2),(1,0),(1,1),(0,1)]).lattice_polytope(True).points()
    45354643            [ 0 -1 -1  1  1  0  0]
    45364644            [-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()
    45384646            [(1, 0), (1, 1), (0, 1), (0, 0)]
    45394647        """
    45404648        if not self.is_compact():
    45414649            raise ValueError, 'Can only enumerate points in a compact polyhedron.'
    45424650
    45434651        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
    45444659
    45454660        if self.is_lattice_polytope():
    45464661            return lp.points().columns()
     
    45494664                        lp.points().columns())
    45504665        return points
    45514666
     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
    45524789   
    45534790    def combinatorial_automorphism_group(self):
    45544791        """
  • sage/geometry/triangulation/point_configuration.py

    diff --git a/sage/geometry/triangulation/point_configuration.py b/sage/geometry/triangulation/point_configuration.py
    a b  
    690690        return Fan(self, (vector(R, p) - origin for p in points))
    691691
    692692
     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
    693715
    694716########################################################################
    695717class PointConfiguration(UniqueRepresentation, PointConfiguration_base):