Ticket #12553: trac_12553_ppl_lattice_polytope.patch

File trac_12553_ppl_lattice_polytope.patch, 52.8 KB (added by vbraun, 9 years ago)

Updated patch

  • sage/geometry/integral_points.pyx

    # HG changeset patch
    # User Volker Braun <vbraun.name@gmail.com>
    # Date 1367768319 -3600
    # Node ID b92dcc9715e6f6c0d65b6d7987309468dd3d6620
    # Parent  ea8f149ec0e88183d31b0f76b95758fda2df3ad6
    Trac #12553: Add interface for PALP polytope databases
    
    A PPL-based lattice polytope class
    
    diff --git a/sage/geometry/integral_points.pyx b/sage/geometry/integral_points.pyx
    a b  
    333333# rectangular bounding box) it is faster to naively enumerate the
    334334# points. This saves the overhead of triangulating the polytope etc.
    335335
    336 def rectangular_box_points(box_min, box_max, polyhedron=None):
     336def rectangular_box_points(box_min, box_max, polyhedron=None,
     337                           count_only=False, return_saturated=False):
    337338    r"""
    338339    Return the integral points in the lattice bounding box that are
    339340    also contained in the given polyhedron.
     
    350351      :class:`~sage.geometry.polyhedron.base.Polyhedron_base`, a PPL
    351352      :class:`~sage.libs.ppl.C_Polyhedron`, or ``None`` (default).
    352353
     354    - ``count_only`` -- Boolean (default: ``False``). Whether to
     355      return only the total number of vertices, and not their
     356      coordinates. Enabling this option speeds up the
     357      enumeration. Cannot be combined with the ``return_saturated``
     358      option.
     359
     360    - ``return_saturated`` -- Boolean (default: ``False``. Whether to
     361      also return which inequalities are saturated for each point of
     362      the polyhedron. Enabling this slows down the enumeration. Cannot
     363      be combined with the ``count_only`` option.
     364
    353365    OUTPUT:
    354366
    355     A tuple containing the integral points of the rectangular box
    356     spanned by `box_min` and `box_max` and that lie inside the
    357     ``polyhedron``. For sufficiently large bounding boxes, this are
    358     all integral points of the polyhedron.
     367    By default, this function returns a tuple containing the integral
     368    points of the rectangular box spanned by `box_min` and `box_max`
     369    and that lie inside the ``polyhedron``. For sufficiently large
     370    bounding boxes, this are all integral points of the polyhedron.
    359371
    360372    If no polyhedron is specified, all integral points of the
    361373    rectangular box are returned.
    362374
     375    If ``count_only`` is specified, only the total number (an integer)
     376    of found lattice points is returned.
     377
     378    If ``return_saturated`` is enabled, then for each integral point a
     379    pair ``(point, Hrep)`` is returned where ``point`` is the point
     380    and ``Hrep`` is the set of indices of the H-representation objects
     381    that are saturated at the point.
     382
    363383    ALGORITHM:
    364384
    365385    This function implements the naive algorithm towards counting
     
    404424         (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
    405425         (1, 2, 0), (1, 2, 1), (1, 2, 2), (1, 2, 3))
    406426
     427        sage: from sage.geometry.integral_points import rectangular_box_points
     428        sage: rectangular_box_points([0,0,0],[1,2,3], count_only=True)
     429        24
     430
    407431        sage: cell24 = polytopes.twenty_four_cell()
    408432        sage: rectangular_box_points([-1]*4, [1]*4, cell24)
    409433        ((-1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0), (0, 0, 0, -1),
     
    419443        sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) )
    420444        3625
    421445
     446        sage: rectangular_box_points([-d]*4, [d]*4, dilated_cell24, count_only=True)
     447        3625
     448
    422449        sage: polytope = Polyhedron([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4)])
    423450        sage: pts = rectangular_box_points([-4]*4, [4]*4, polytope); pts
    424451        ((-4, -3, -2, -1), (-1, 0, 0, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 1, 3, 0),
     
    446473        201
    447474
    448475    Using a PPL polyhedron::
    449    
     476
    450477        sage: from sage.libs.ppl import Variable, Generator_System, C_Polyhedron, point
    451478        sage: gs = Generator_System()
    452479        sage: x = Variable(0); y = Variable(1); z = Variable(2)
     
    458485        sage: rectangular_box_points([0]*3, [3]*3, poly)
    459486        ((0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
    460487         (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 1, 3), (3, 1, 0), (3, 1, 1), (3, 1, 2), (3, 1, 3))
     488
     489    Optionally, return the information about the saturated inequalities as well::
     490
     491        sage: cube = polytopes.n_cube(3)
     492        sage: cube.Hrepresentation(0)             
     493        An inequality (0, 0, -1) x + 1 >= 0
     494        sage: cube.Hrepresentation(1)
     495        An inequality (0, -1, 0) x + 1 >= 0
     496        sage: cube.Hrepresentation(2)
     497        An inequality (-1, 0, 0) x + 1 >= 0
     498        sage: rectangular_box_points([0]*3, [1]*3, cube, return_saturated=True)
     499        (((0, 0, 0), frozenset([])),
     500         ((0, 0, 1), frozenset([0])),
     501         ((0, 1, 0), frozenset([1])),
     502         ((0, 1, 1), frozenset([0, 1])),
     503         ((1, 0, 0), frozenset([2])),
     504         ((1, 0, 1), frozenset([0, 2])),
     505         ((1, 1, 0), frozenset([1, 2])),
     506         ((1, 1, 1), frozenset([0, 1, 2])))
    461507    """
    462508    assert len(box_min)==len(box_max)
     509    assert not (count_only and return_saturated)
    463510    cdef int d = len(box_min)
    464511    diameter = sorted([ (box_max[i]-box_min[i], i) for i in range(0,d) ], reverse=True)
    465512    diameter_value = [ x[0] for x in diameter ]
     
    472519    box_max = sort_perm.action(box_max)
    473520    inequalities = InequalityCollection(polyhedron, sort_perm, box_min, box_max)
    474521
     522    if count_only:
     523        return loop_over_rectangular_box_points(box_min, box_max, inequalities, d, count_only)
     524
     525    points = []
    475526    v = vector(ZZ, d)
    476     points = []
    477     for p in loop_over_rectangular_box_points(box_min, box_max, inequalities, d):
    478         #  v = vector(ZZ, orig_perm.action(p))   # too slow
    479         for i in range(0,d):
    480             v[i] = p[orig_perm_list[i]]
    481         points.append(copy.copy(v))
     527    if not return_saturated:
     528        for p in loop_over_rectangular_box_points(box_min, box_max, inequalities, d, count_only):
     529            #  v = vector(ZZ, orig_perm.action(p))   # too slow
     530            for i in range(0,d):
     531                v.set(i, p[orig_perm_list[i]])
     532            v_copy = copy.copy(v)
     533            v_copy.set_immutable()
     534            points.append(v_copy)
     535    else:
     536        for p, saturated in \
     537                loop_over_rectangular_box_points_saturated(box_min, box_max, inequalities, d):
     538            for i in range(0,d):
     539                v.set(i, p[orig_perm_list[i]])
     540            v_copy = copy.copy(v)
     541            v_copy.set_immutable()
     542            points.append( (v_copy, saturated) )
     543
    482544    return tuple(points)
    483545
    484546
    485 cdef loop_over_rectangular_box_points(box_min, box_max, inequalities, int d):
     547cdef loop_over_rectangular_box_points(box_min, box_max, inequalities, int d, bint count_only):
    486548    """
    487549    The inner loop of :func:`rectangular_box_points`.
    488550   
     
    495557
    496558    - ``d`` -- the ambient space dimension.
    497559
     560    - ``count_only`` -- whether to only return the total number of
     561      lattice points.
     562
    498563    OUTPUT:
    499564
    500565    The integral points in the bounding box satisfying all
    501566    inequalities.
    502567    """
    503568    cdef int inc
     569    if count_only:
     570        points = 0
     571    else:
     572        points = []
     573    p = copy.copy(box_min)
     574    inequalities.prepare_next_to_inner_loop(p)
     575    while True:
     576        inequalities.prepare_inner_loop(p)
     577        i_min = box_min[0]
     578        i_max = box_max[0]
     579        # Find the lower bound for the allowed region
     580        while i_min <= i_max:
     581            if inequalities.are_satisfied(i_min):
     582                break
     583            i_min += 1
     584        # Find the upper bound for the allowed region
     585        while i_min <= i_max:
     586            if inequalities.are_satisfied(i_max):
     587                break
     588            i_max -= 1
     589        # The points i_min .. i_max are contained in the polyhedron
     590        if count_only:
     591            if i_max>=i_min:
     592                points += i_max-i_min+1
     593        else:
     594            i = i_min
     595            while i <= i_max:
     596                p[0] = i
     597                points.append(tuple(p))
     598                i += 1
     599        # finally increment the other entries in p to move on to next inner loop
     600        inc = 1
     601        if d==1: return points
     602        while True:
     603            if p[inc]==box_max[inc]:
     604                p[inc] = box_min[inc]
     605                inc += 1
     606                if inc==d:
     607                    return points
     608            else:
     609                p[inc] += 1
     610                break
     611        if inc>1:
     612            inequalities.prepare_next_to_inner_loop(p)
     613
     614
     615
     616cdef loop_over_rectangular_box_points_saturated(box_min, box_max, inequalities, int d):
     617    """
     618    The analog of :func:`rectangular_box_points` except that it keeps
     619    track of which inequalities are saturated.
     620   
     621    INPUT:
     622
     623    See :func:`rectangular_box_points`.
     624
     625    OUTPUT:
     626
     627    The integral points in the bounding box satisfying all
     628    inequalities, each point being returned by a coordinate vector and
     629    a set of H-representation object indices.
     630    """
     631    cdef int inc
    504632    points = []
    505633    p = copy.copy(box_min)
    506634    inequalities.prepare_next_to_inner_loop(p)
     
    522650        i = i_min
    523651        while i <= i_max:
    524652            p[0] = i
    525             points.append(tuple(p))
     653            saturated = inequalities.satisfied_as_equalities(i)
     654            points.append( (tuple(p), saturated) )
    526655            i += 1
    527656        # finally increment the other entries in p to move on to next inner loop
    528657        inc = 1
     
    651780        return inner_loop_variable*self.coeff + self.cache == 0       
    652781
    653782
    654 
     783# if dim>20 then we always use the generic inequalities (Inequality_generic)
    655784DEF INEQ_INT_MAX_DIM = 20
    656785
    657786cdef class Inequality_int:
     
    11761305
    11771306        OUTPUT:
    11781307       
    1179         A tuple of integers in ascending order. Each integer is the
    1180         index of a H-representation object of the
    1181         polyhedron. Equalities are treated as a pair of opposite
    1182         inequalities. 
     1308        A set of integers in ascending order. Each integer is the
     1309        index of a H-representation object of the polyhedron (either a
     1310        inequality or an equation).
    11831311       
    11841312        EXAMPLES::
    11851313       
     
    11891317            sage: ieqs.prepare_next_to_inner_loop([-1,0])
    11901318            sage: ieqs.prepare_inner_loop([-1,0])
    11911319            sage: ieqs.satisfied_as_equalities(-1)
    1192             (1,)
     1320            frozenset([1])
    11931321            sage: ieqs.satisfied_as_equalities(0)
    1194             (0, 1)
     1322            frozenset([0, 1])
    11951323            sage: ieqs.satisfied_as_equalities(1)
    1196             (1,)
     1324            frozenset([1])
    11971325        """
    11981326        cdef int i
    11991327        result = []
     
    12051333            ineq = self.ineqs_generic[i]
    12061334            if (<Inequality_generic>ineq).is_equality(inner_loop_variable):
    12071335                result.append( (<Inequality_generic>ineq).index )
    1208         return tuple(uniq(result))
     1336        return frozenset(result)
    12091337   
    12101338
    12111339
  • new file sage/geometry/polyhedron/ppl_lattice_polytope.py

    diff --git a/sage/geometry/polyhedron/ppl_lattice_polytope.py b/sage/geometry/polyhedron/ppl_lattice_polytope.py
    new file mode 100644
    - +  
     1"""
     2Fast Lattice Polytopes using PPL.
     3
     4The :func:`LatticePolytope_PPL` class is a thin wrapper around PPL
     5polyhedra. Its main purpose is to be fast to construct, at the cost of
     6being much less full-featured than the usual polyhedra. This makes it
     7possible to iterate with it over the list of all 473800776 reflexive
     8polytopes in 4 dimensions.
     9
     10.. NOTE::
     11
     12    For general lattice polyhedra you should use
     13    :func:`~sage.geometry.polyhedon.constructor.Polyhedron` with
     14    `base_ring=ZZ`.
     15
     16EXAMPLES::
     17
     18    sage: vertices = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1), (-9, -6, -1, -1)]
     19    sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     20    sage: P = LatticePolytope_PPL(vertices);  P
     21    A 4-dimensional lattice polytope in ZZ^4 with 5 vertices
     22    sage: P.integral_points()
     23    ((-9, -6, -1, -1), (-3, -2, 0, 0), (-2, -1, 0, 0), (-1, -1, 0, 0),
     24     (-1, 0, 0, 0), (0, 0, 0, 0), (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 0, 1), (0, 0, 1, 0))
     25    sage: P.integral_points_not_interior_to_facets()
     26    ((-9, -6, -1, -1), (-3, -2, 0, 0), (0, 0, 0, 0), (1, 0, 0, 0),
     27     (0, 1, 0, 0), (0, 0, 0, 1), (0, 0, 1, 0))
     28
     29Fibrations of the lattice polytopes are defined as lattice
     30sub-polytopes and give rise to fibrations of toric varieties for
     31suitable fan refinements. We can compute them using
     32:meth:`~LatticePolytope_PPL.fibration_generator` ::
     33
     34    sage: F = P.fibration_generator(2).next()
     35    sage: F.vertices()
     36    ((1, 0, 0, 0), (0, 1, 0, 0), (-3, -2, 0, 0))
     37
     38Finally, we can compute automorphisms and identify fibrations that
     39only differ by a lattice automorphism::
     40
     41    sage: square = LatticePolytope_PPL((-1,-1),(-1,1),(1,-1),(1,1))
     42    sage: fibers = [ f.vertices() for f in square.fibration_generator(1) ];  fibers
     43    [((1, 0), (-1, 0)), ((0, 1), (0, -1)), ((-1, -1), (1, 1)), ((-1, 1), (1, -1))]
     44    sage: square.pointsets_mod_automorphism(fibers)
     45    (frozenset([(0, 1), (0, -1)]), frozenset([(1, 1), (-1, -1)]))
     46
     47AUTHORS:
     48
     49    - Volker Braun: initial version, 2012
     50"""
     51
     52########################################################################
     53#       Copyright (C) 2012 Volker Braun <vbraun.name@gmail.com>
     54#
     55#  Distributed under the terms of the GNU General Public License (GPL)
     56#
     57#                  http://www.gnu.org/licenses/
     58########################################################################
     59
     60
     61
     62import copy
     63from sage.rings.integer import GCD_list
     64from sage.rings.integer_ring import ZZ
     65from sage.misc.all import union, cached_method, prod, uniq
     66from sage.matrix.constructor import matrix, column_matrix, vector, diagonal_matrix
     67from sage.libs.ppl import (
     68    C_Polyhedron, Linear_Expression, Variable,
     69    point, ray, line, Generator, Generator_System,
     70    Constraint_System,
     71    Poly_Con_Relation)
     72
     73
     74
     75
     76########################################################################
     77def LatticePolytope_PPL(*args):
     78    """
     79    Construct a new instance of the PPL-based lattice polytope class.
     80
     81    EXAMPLES::
     82
     83        sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     84        sage: LatticePolytope_PPL((0,0),(1,0),(0,1))
     85        A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     86
     87        sage: from sage.libs.ppl import point, Generator_System, C_Polyhedron, Linear_Expression, Variable
     88        sage: p = point(Linear_Expression([2,3],0));  p
     89        point(2/1, 3/1)
     90        sage: LatticePolytope_PPL(p)
     91        A 0-dimensional lattice polytope in ZZ^2 with 1 vertex
     92
     93        sage: P = C_Polyhedron(Generator_System(p));  P
     94        A 0-dimensional polyhedron in QQ^2 defined as the convex hull of 1 point
     95        sage: LatticePolytope_PPL(P)
     96        A 0-dimensional lattice polytope in ZZ^2 with 1 vertex
     97
     98    A ``TypeError`` is raised if the arguments do not specify a lattice polytope::
     99
     100        sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     101        sage: LatticePolytope_PPL((0,0),(1/2,1))
     102        Traceback (most recent call last):
     103        ...
     104        TypeError: no conversion of this rational to integer
     105
     106        sage: from sage.libs.ppl import point, Generator_System, C_Polyhedron, Linear_Expression, Variable
     107        sage: p = point(Linear_Expression([2,3],0), 5);  p
     108        point(2/5, 3/5)
     109        sage: LatticePolytope_PPL(p)
     110        Traceback (most recent call last):
     111        ...
     112        TypeError: The generator is not a lattice polytope generator.
     113
     114        sage: P = C_Polyhedron(Generator_System(p));  P
     115        A 0-dimensional polyhedron in QQ^2 defined as the convex hull of 1 point
     116        sage: LatticePolytope_PPL(P)
     117        Traceback (most recent call last):
     118        ...
     119        TypeError: The polyhedron has non-integral generators.
     120    """
     121    if len(args)==1 and isinstance(args[0], C_Polyhedron):
     122        polyhedron = args[0]
     123        #if not polyhedron.is_bounded():
     124        #    raise TypeError('The polyhedron is unbounded.')
     125        if not all(p.is_point() and p.divisor().is_one() for p in polyhedron.generators()):
     126            raise TypeError('The polyhedron has non-integral generators.')
     127        return LatticePolytope_PPL_class(polyhedron)
     128    if len(args)==1:
     129        vertices = args[0]
     130        try:
     131            return LatticePolytope_PPL(*vertices)
     132        except TypeError:
     133            pass
     134    vertices = args
     135    gs = Generator_System()
     136    for v in vertices:
     137        if isinstance(v, Generator):
     138            if (not v.is_point()) or (not v.divisor().is_one()):
     139                raise TypeError('The generator is not a lattice polytope generator.')
     140            gs.insert(v)
     141        else:
     142            gs.insert(point(Linear_Expression(v, 0)))
     143    return LatticePolytope_PPL_class(gs)
     144
     145
     146########################################################################
     147class LatticePolytope_PPL_class(C_Polyhedron):
     148    """
     149    The lattice polytope class.
     150
     151    You should use :func:LatticePolytope_PPL` to construct instances.
     152
     153    EXAMPLES::
     154
     155        sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     156        sage: LatticePolytope_PPL((0,0),(1,0),(0,1))
     157        A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     158    """
     159
     160    def _repr_(self):
     161        """
     162        Return the string representation
     163
     164        OUTPUT:
     165
     166        A string.
     167
     168        EXAMPLES::
     169
     170            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     171            sage: P = LatticePolytope_PPL((0,0),(1,0),(0,1))
     172            sage: P
     173            A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     174            sage: P._repr_()
     175            'A 2-dimensional lattice polytope in ZZ^2 with 3 vertices'
     176
     177            sage: LatticePolytope_PPL()
     178            The empty lattice polytope in ZZ^0
     179        """
     180        desc = ''
     181        if self.n_vertices()==0:
     182            desc += 'The empty lattice polytope'
     183        else:
     184            desc += 'A ' + repr(self.affine_dimension()) + '-dimensional lattice polytope'
     185        desc += ' in ZZ^' + repr(self.space_dimension())
     186
     187        if self.n_vertices()>0:
     188            desc += ' with '
     189            desc += repr(self.n_vertices())
     190            if self.n_vertices()==1: desc += ' vertex'
     191            else:                    desc += ' vertices'
     192        return desc
     193
     194
     195
     196    def is_bounded(self):
     197        """
     198        Return whether the lattice polytope is compact.
     199
     200        OUTPUT:
     201
     202        Always ``True``, since polytopes are by definition compact.
     203
     204        EXAMPLES::
     205
     206            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     207            sage: LatticePolytope_PPL((0,0),(1,0),(0,1)).is_bounded()
     208            True
     209        """
     210        return True
     211
     212    @cached_method
     213    def n_vertices(self):
     214        """
     215        Return the number of vertices.
     216
     217        OUTPUT:
     218
     219        An integer, the number of vertices.
     220
     221        EXAMPLES::
     222
     223            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     224            sage: LatticePolytope_PPL((0,0,0), (1,0,0), (0,1,0)).n_vertices()
     225            3
     226        """
     227        return len(self.minimized_generators())
     228
     229    @cached_method
     230    def is_simplex(self):
     231        r"""
     232        Return whether the polyhedron is a simplex.
     233
     234        OUTPUT:
     235
     236        Boolean, whether the polyhedron is a simplex (possibly of
     237        strictly smaller dimension than the ambient space).
     238
     239        EXAMPLES::
     240
     241            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     242            sage: LatticePolytope_PPL((0,0,0), (1,0,0), (0,1,0)).is_simplex()
     243            True
     244        """
     245        return self.affine_dimension()+1==self.n_vertices()
     246
     247    @cached_method
     248    def bounding_box(self):
     249        r"""
     250        Return the coordinates of a rectangular box containing the non-empty polytope.
     251
     252        OUTPUT:
     253
     254        A pair of tuples ``(box_min, box_max)`` where ``box_min`` are
     255        the coordinates of a point bounding the coordinates of the
     256        polytope from below and ``box_max`` bounds the coordinates
     257        from above.
     258
     259        EXAMPLES::
     260
     261            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     262            sage: LatticePolytope_PPL((0,0),(1,0),(0,1)).bounding_box()
     263            ((0, 0), (1, 1))
     264        """
     265        box_min = []
     266        box_max = []
     267        if self.is_empty():
     268            raise ValueError('Empty polytope is not allowed')
     269        for i in range(0, self.space_dimension()):
     270            x = Variable(i)
     271            coords = [ v.coefficient(x) for v in self.generators() ]
     272            max_coord = max(coords)
     273            min_coord = min(coords)
     274            box_max.append(max_coord)
     275            box_min.append(min_coord)
     276        return (tuple(box_min), tuple(box_max))
     277
     278    @cached_method
     279    def n_integral_points(self):
     280        """
     281        Return the number of integral points.
     282
     283        OUTPUT:
     284
     285        Integer. The number of integral points contained in the
     286        lattice polytope.
     287
     288        EXAMPLES::
     289
     290            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     291            sage: LatticePolytope_PPL((0,0),(1,0),(0,1)).n_integral_points()
     292            3
     293        """
     294        if self.is_empty():
     295            return tuple()
     296        box_min, box_max = self.bounding_box()
     297        from sage.geometry.integral_points import rectangular_box_points
     298        return rectangular_box_points(box_min, box_max, self, count_only=True)
     299
     300    @cached_method
     301    def integral_points(self):
     302        r"""
     303        Return the integral points in the polyhedron.
     304
     305        Uses the naive algorithm (iterate over a rectangular bounding
     306        box).
     307
     308        OUTPUT:
     309
     310        The list of integral points in the polyhedron. If the
     311        polyhedron is not compact, a ``ValueError`` is raised.
     312
     313        EXAMPLES::
     314
     315            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     316            sage: LatticePolytope_PPL((-1,-1),(1,0),(1,1),(0,1)).integral_points()
     317            ((-1, -1), (0, 0), (0, 1), (1, 0), (1, 1))
     318
     319            sage: simplex = LatticePolytope_PPL((1,2,3), (2,3,7), (-2,-3,-11))
     320            sage: simplex.integral_points()
     321            ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7))
     322
     323        The polyhedron need not be full-dimensional::
     324
     325            sage: simplex = LatticePolytope_PPL((1,2,3,5), (2,3,7,5), (-2,-3,-11,5))
     326            sage: simplex.integral_points()
     327            ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5))
     328
     329            sage: point = LatticePolytope_PPL((2,3,7))
     330            sage: point.integral_points()
     331            ((2, 3, 7),)
     332
     333            sage: empty = LatticePolytope_PPL()
     334            sage: empty.integral_points()
     335            ()
     336
     337        Here is a simplex where the naive algorithm of running over
     338        all points in a rectangular bounding box no longer works fast
     339        enough::
     340
     341            sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)]
     342            sage: simplex = LatticePolytope_PPL(v); simplex
     343            A 4-dimensional lattice polytope in ZZ^4 with 5 vertices
     344            sage: len(simplex.integral_points())
     345            49
     346
     347        Finally, the 3-d reflexive polytope number 4078::
     348
     349            sage: v = [(1,0,0), (0,1,0), (0,0,1), (0,0,-1), (0,-2,1),
     350            ...        (-1,2,-1), (-1,2,-2), (-1,1,-2), (-1,-1,2), (-1,-3,2)]
     351            sage: P = LatticePolytope_PPL(*v)
     352            sage: pts1 = P.integral_points()                     # Sage's own code
     353            sage: pts2 = LatticePolytope(v).points().columns()   # PALP
     354            sage: for p in pts1: p.set_immutable()
     355            sage: for p in pts2: p.set_immutable()
     356            sage: set(pts1) == set(pts2)
     357            True
     358
     359            sage: timeit('Polyhedron(v).integral_points()')   # random output
     360            sage: timeit('LatticePolytope(v).points()')       # random output
     361            sage: timeit('LatticePolytope_PPL(*v).integral_points()')       # random output
     362        """
     363        if self.is_empty():
     364            return tuple()
     365        box_min, box_max = self.bounding_box()
     366        from sage.geometry.integral_points import rectangular_box_points
     367        points = rectangular_box_points(box_min, box_max, self)
     368        if not self.n_integral_points.is_in_cache():
     369            self.n_integral_points.set_cache(len(points))
     370        return points
     371
     372    @cached_method
     373    def _integral_points_saturating(self):
     374        """
     375        Return the integral points together with information about
     376        which inequalities are saturated.
     377
     378        See :func:`~sage.geometry.integral_points.rectangular_box_points`.
     379
     380        OUTPUT:
     381
     382        A tuple of pairs (one for each integral point) consisting of a
     383        pair ``(point, Hrep)``, where ``point`` is the coordinate
     384        vector of the intgeral point and ``Hrep`` is the set of
     385        indices of the :meth:`minimized_constraints` that are
     386        saturated at the point.
     387
     388        EXAMPLES::
     389
     390            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     391            sage: quad = LatticePolytope_PPL((-1,-1),(0,1),(1,0),(1,1))
     392            sage: quad._integral_points_saturating()
     393            (((-1, -1), frozenset([0, 1])),
     394             ((0, 0), frozenset([])),
     395             ((0, 1), frozenset([0, 3])),
     396             ((1, 0), frozenset([1, 2])),
     397             ((1, 1), frozenset([2, 3])))
     398        """
     399        if self.is_empty():
     400            return tuple()
     401        box_min, box_max = self.bounding_box()
     402        from sage.geometry.integral_points import rectangular_box_points
     403        points= rectangular_box_points(box_min, box_max, self, return_saturated=True)
     404        if not self.n_integral_points.is_in_cache():
     405            self.n_integral_points.set_cache(len(points))
     406        if not self.integral_points.is_in_cache():
     407            self.integral_points.set_cache(tuple(p[0] for p in points))
     408        return points
     409
     410    @cached_method
     411    def integral_points_not_interior_to_facets(self):
     412        """
     413        Return the integral points not interior to facets
     414
     415        OUTPUT:
     416
     417        A tuple whose entries are the coordinate vectors of integral
     418        points not interior to facets (codimension one faces) of the
     419        lattice polytope.
     420
     421        EXAMPLES::
     422
     423            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     424            sage: square = LatticePolytope_PPL((-1,-1),(-1,1),(1,-1),(1,1))
     425            sage: square.n_integral_points()
     426            9
     427            sage: square.integral_points_not_interior_to_facets()
     428            ((-1, -1), (-1, 1), (0, 0), (1, -1), (1, 1))
     429        """
     430        n = 1 + self.space_dimension() - self.affine_dimension()
     431        return tuple(p[0] for p in self._integral_points_saturating() if len(p[1])!=n)
     432
     433    @cached_method
     434    def vertices(self):
     435        r"""
     436        Return the vertices as a tuple of `\ZZ`-vectors.
     437
     438        OUTPUT:
     439
     440        A tuple of `\ZZ`-vectors. Each entry is the coordinate vector
     441        of an integral points of the lattice polytope.
     442
     443        EXAMPLES::
     444
     445            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     446            sage: p = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
     447            sage: p.vertices()
     448            ((-9, -6, -1, -1), (0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0))
     449            sage: p.minimized_generators()
     450            Generator_System {point(-9/1, -6/1, -1/1, -1/1), point(0/1, 0/1, 0/1, 1/1),
     451            point(0/1, 0/1, 1/1, 0/1), point(0/1, 1/1, 0/1, 0/1), point(1/1, 0/1, 0/1, 0/1)}
     452        """
     453        d = self.space_dimension()
     454        v = vector(ZZ, d)
     455        points = []
     456        for g in self.minimized_generators():
     457            for i in range(0,d):
     458                v[i] = g.coefficient(Variable(i))
     459            v_copy = copy.copy(v)
     460            v_copy.set_immutable()
     461            points.append(v_copy)
     462        return tuple(points)
     463
     464    @cached_method
     465    def is_full_dimensional(self):
     466        """
     467        Return whether the lattice polytope is full dimensional.
     468
     469        OUTPUT:
     470
     471        Boolean. Whether the :meth:`affine_dimension` equals the
     472        ambient space dimension.
     473
     474        EXAMPLES::
     475
     476            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     477            sage: p = LatticePolytope_PPL((0,0),(0,1))
     478            sage: p.is_full_dimensional()
     479            False
     480            sage: q = LatticePolytope_PPL((0,0),(0,1),(1,0))
     481            sage: q.is_full_dimensional()
     482            True
     483        """
     484
     485        return self.affine_dimension() == self.space_dimension()
     486
     487    def fibration_generator(self, dim):
     488        """
     489        Generate the lattice polytope fibrations.
     490
     491        For the purposes of this function, a lattice polytope fiber is
     492        a sub-lattice polytope. Projecting the plane spanned by the
     493        subpolytope to a point yields another lattice polytope, the
     494        base of the fibration.
     495
     496        INPUT:
     497
     498        - ``dim`` -- integer. The dimension of the lattice polytope
     499          fiber.
     500
     501        OUTPUT:
     502
     503        A generator yielding the distinct lattice polytope fibers of
     504        given dimension.
     505
     506        EXAMPLES::
     507
     508            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     509            sage: p = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
     510            sage: list( p.fibration_generator(2) )
     511            [A 2-dimensional lattice polytope in ZZ^4 with 3 vertices]
     512        """
     513        assert self.is_full_dimensional()
     514        codim = self.space_dimension() - dim
     515        # "points" are the potential vertices of the fiber. They are
     516        # in the $codim$-skeleton of the polytope, which is contained
     517        # in the points that saturate at least $dim$ equations.
     518        points = [ p for p in self._integral_points_saturating() if len(p[1])>=dim ]
     519        points = sorted(points, key=lambda x:len(x[1]))
     520
     521        # iterate over point combinations subject to all points being on one facet.
     522        def point_combinations_iterator(n, i0=0, saturated=None):
     523            for i in range(i0, len(points)):
     524                p, ieqs = points[i]
     525                if saturated is None:
     526                    saturated_ieqs = ieqs
     527                else:
     528                    saturated_ieqs = saturated.intersection(ieqs)
     529                if len(saturated_ieqs)==0:
     530                    continue
     531                if n == 1:
     532                    yield [i]
     533                else:
     534                    for c in point_combinations_iterator(n-1, i+1, saturated_ieqs):
     535                        yield [i] + c
     536
     537        point_lines = [ line(Linear_Expression(p[0].list(),0)) for p in points ]
     538        origin = point()
     539        fibers = set()
     540        gs = Generator_System()
     541        for indices in point_combinations_iterator(dim):
     542            gs.clear()
     543            gs.insert(origin)
     544            for i in indices:
     545                gs.insert(point_lines[i])
     546            plane = C_Polyhedron(gs)
     547            if plane.affine_dimension() != dim:
     548                continue
     549            plane.intersection_assign(self)
     550            if (not self.is_full_dimensional()) and (plane.affine_dimension() != dim):
     551                continue
     552            try:
     553                fiber = LatticePolytope_PPL(plane)
     554            except TypeError:   # not a lattice polytope
     555                continue
     556            fiber_vertices = tuple(sorted(fiber.vertices()))
     557            if fiber_vertices not in fibers:
     558                yield fiber
     559                fibers.update([fiber_vertices])
     560
     561    def pointsets_mod_automorphism(self, pointsets):
     562        """
     563        Return ``pointsets`` modulo the automorphisms of ``self``.
     564
     565        INPUT:
     566
     567        - ``polytopes`` a tuple/list/iterable of subsets of the
     568          integral points of ``self``.
     569
     570        OUTPUT:
     571
     572        Representatives of the point sets modulo the
     573        :meth:`lattice_automorphism_group` of ``self``.
     574
     575        EXAMPLES::
     576
     577            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     578            sage: square = LatticePolytope_PPL((-1,-1),(-1,1),(1,-1),(1,1))
     579            sage: fibers = [ f.vertices() for f in square.fibration_generator(1) ]
     580            sage: square.pointsets_mod_automorphism(fibers)
     581            (frozenset([(0, 1), (0, -1)]), frozenset([(1, 1), (-1, -1)]))
     582
     583            sage: cell24 = LatticePolytope_PPL(
     584            ...   (1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,-1,-1,1),(0,0,-1,1),
     585            ...   (0,-1,0,1),(-1,0,0,1),(1,0,0,-1),(0,1,0,-1),(0,0,1,-1),(-1,1,1,-1),
     586            ...   (1,-1,-1,0),(0,0,-1,0),(0,-1,0,0),(-1,0,0,0),(1,-1,0,0),(1,0,-1,0),
     587            ...   (0,1,1,-1),(-1,1,1,0),(-1,1,0,0),(-1,0,1,0),(0,-1,-1,1),(0,0,0,-1))
     588            sage: fibers = [ f.vertices() for f in cell24.fibration_generator(2) ]
     589            sage: cell24.pointsets_mod_automorphism(fibers)   # long time
     590            (frozenset([(1, 0, -1, 0), (-1, 0, 1, 0), (0, -1, -1, 1), (0, 1, 1, -1)]),
     591             frozenset([(-1, 0, 0, 0), (1, 0, 0, 0), (0, 0, 0, 1),
     592                        (1, 0, 0, -1), (0, 0, 0, -1), (-1, 0, 0, 1)]))
     593        """
     594        points = set()
     595        for ps in pointsets:
     596            points.update(ps)
     597        points = tuple(points)
     598        Aut = self.lattice_automorphism_group(points,
     599                                              point_labels=tuple(range(len(points))))
     600        indexsets = set([ frozenset([points.index(p) for p in ps]) for ps in pointsets ])
     601        orbits = []
     602        while len(indexsets)>0:
     603            idx = indexsets.pop()
     604            orbits.append(frozenset([points[i] for i in idx]))
     605            for g in Aut:
     606                g_idx = frozenset([g(i) for i in idx])
     607                indexsets.difference_update([g_idx])
     608        return tuple(orbits)
     609
     610    @cached_method
     611    def ambient_space(self):
     612        r"""
     613        Return the ambient space.
     614
     615        OUTPUT:
     616
     617        The free module `\ZZ^d`, where `d` is the ambient space
     618        dimension.
     619
     620        EXAMPLES::
     621
     622            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     623            sage: point = LatticePolytope_PPL((1,2,3))
     624            sage: point.ambient_space()
     625            Ambient free module of rank 3 over the principal ideal domain Integer Ring
     626        """
     627        from sage.modules.free_module import FreeModule
     628        return FreeModule(ZZ, self.space_dimension())
     629
     630    def contains(self, point_coordinates):
     631        r"""
     632        Test whether point is contained in the polytope.
     633       
     634        INPUT:
     635 
     636        - ``point_coordinates`` -- a list/tuple/iterable of rational
     637          numbers. The coordinates of the point.
     638
     639        OUTPUT:
     640
     641        Boolean.
     642       
     643        EXAMPLES::
     644
     645            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     646            sage: line = LatticePolytope_PPL((1,2,3), (-1,-2,-3))
     647            sage: line.contains([0,0,0])
     648            True
     649            sage: line.contains([1,0,0])
     650            False
     651        """
     652        p = C_Polyhedron(point(Linear_Expression(list(point_coordinates), 1)))
     653        is_included = Poly_Con_Relation.is_included()
     654        for c in self.constraints():
     655            if not p.relation_with(c).implies(is_included):
     656                return False
     657        return True
     658       
     659    @cached_method
     660    def contains_origin(self):
     661        """
     662        Test whether the polytope contains the origin
     663
     664        OUTPUT:
     665       
     666        Boolean.
     667
     668        EXAMPLES::
     669
     670            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     671            sage: LatticePolytope_PPL((1,2,3), (-1,-2,-3)).contains_origin()
     672            True
     673            sage: LatticePolytope_PPL((1,2,5), (-1,-2,-3)).contains_origin()
     674            False
     675        """
     676        return self.contains(self.ambient_space().zero())
     677
     678    @cached_method
     679    def affine_space(self):
     680        r"""
     681        Return the affine space spanned by the polytope.
     682
     683        OUTPUT:
     684
     685        The free module `\ZZ^n`, where `n` is the dimension of the
     686        affine space spanned by the points of the polytope.
     687
     688        EXAMPLES::
     689
     690            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     691            sage: point = LatticePolytope_PPL((1,2,3))
     692            sage: point.affine_space()
     693            Free module of degree 3 and rank 0 over Integer Ring
     694            Echelon basis matrix:
     695            []
     696            sage: line = LatticePolytope_PPL((1,1,1), (1,2,3))
     697            sage: line.affine_space()
     698            Free module of degree 3 and rank 1 over Integer Ring
     699            Echelon basis matrix:
     700            [0 1 2]
     701        """
     702        vertices = self.vertices()
     703        if not self.contains_origin():
     704            v0 = vertices[0]
     705            vertices = [v-v0 for v in vertices]
     706        return self.ambient_space().span(vertices).saturation()
     707
     708    def affine_lattice_polytope(self):
     709        """
     710        Return the lattice polytope restricted to
     711        :meth:`affine_space`.
     712
     713        OUTPUT:
     714
     715        A new, full-dimensional lattice polytope.
     716
     717        EXAMPLES::
     718
     719            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     720            sage: poly_4d = LatticePolytope_PPL((-9,-6,0,0),(0,1,0,0),(1,0,0,0));  poly_4d
     721            A 2-dimensional lattice polytope in ZZ^4 with 3 vertices
     722            sage: poly_4d.space_dimension()
     723            4
     724            sage: poly_2d = poly_4d.affine_lattice_polytope();  poly_2d
     725            A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     726            sage: poly_2d.space_dimension()
     727            2
     728        """
     729        V = self.affine_space()
     730        if self.contains_origin():
     731            vertices = [ V.coordinates(v) for v in self.vertices() ]
     732        else:
     733            v0 = vertices[0]
     734            vertices = [ V.coordinates(v-v0) for v in self.vertices() ]
     735        return LatticePolytope_PPL(*vertices)
     736
     737    @cached_method
     738    def base_projection(self, fiber):
     739        """
     740        The projection that maps the sub-polytope ``fiber`` to a
     741        single point.
     742
     743        OUTPUT:
     744
     745        The quotient module of the ambient space modulo the
     746        :meth:`affine_space` spanned by the fiber.
     747
     748        EXAMPLES::
     749
     750            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     751            sage: poly = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
     752            sage: fiber = poly.fibration_generator(2).next()
     753            sage: poly.base_projection(fiber)
     754            Finitely generated module V/W over Integer Ring with invariants (0, 0)
     755        """
     756        return self.ambient_space().quotient(fiber.affine_space())
     757
     758    @cached_method
     759    def base_projection_matrix(self, fiber):
     760        """
     761        The projection that maps the sub-polytope ``fiber`` to a
     762        single point.
     763
     764        OUTPUT:
     765
     766        An integer matrix that represents the projection to the
     767        base.
     768
     769        .. SEEALSO::
     770       
     771            The :meth:`base_projection` yields equivalent information,
     772            and is easier to use. However, just returning the matrix
     773            has lower overhead.
     774
     775        EXAMPLES::
     776
     777            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     778            sage: poly = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
     779            sage: fiber = poly.fibration_generator(2).next()
     780            sage: poly.base_projection_matrix(fiber)
     781            [0 0 1 0]
     782            [0 0 0 1]
     783
     784        Note that the basis choice in :meth:`base_projection` for the
     785        quotient is usually different::
     786
     787            sage: proj = poly.base_projection(fiber)
     788            sage: proj_matrix = poly.base_projection_matrix(fiber)
     789            sage: [ proj(p) for p in poly.integral_points() ]
     790            [(-1, -1), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (1, 0), (0, 1)]
     791            sage: [ proj_matrix*p for p in poly.integral_points() ]
     792            [(-1, -1), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (1, 0)]
     793        """
     794        return matrix(ZZ, fiber.vertices()).right_kernel_matrix()
     795
     796    def base_rays(self, fiber, points):
     797        """
     798        Return the primitive lattice vectors that generate the
     799        direction given by the base projection of points.
     800
     801        INPUT:
     802
     803        - ``fiber`` -- a sub-lattice polytope defining the
     804          :meth:`base_projection`.
     805
     806        - ``points`` -- the points to project to the base.
     807
     808        OUTPUT:
     809
     810        A tuple of primitive `\ZZ`-vectors.
     811
     812        EXAMPLES::
     813
     814            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     815            sage: poly = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
     816            sage: fiber = poly.fibration_generator(2).next()
     817            sage: poly.base_rays(fiber, poly.integral_points_not_interior_to_facets())
     818            ((-1, -1), (0, 1), (1, 0))
     819
     820            sage: p = LatticePolytope_PPL((1,0),(1,2),(-1,0))
     821            sage: f = LatticePolytope_PPL((1,0),(-1,0))
     822            sage: p.base_rays(f, p.integral_points())
     823            ((1),)
     824        """
     825        quo = self.base_projection(fiber)
     826        vertices = []
     827        for p in points:
     828            v = vector(ZZ,quo(p))
     829            if v.is_zero():
     830                continue
     831            d =  GCD_list(v.list())
     832            if d>1:
     833                for i in range(0,v.degree()):
     834                    v[i] /= d
     835            v.set_immutable()
     836            vertices.append(v)
     837        return tuple(uniq(vertices))
     838
     839    @cached_method
     840    def has_IP_property(self):
     841        """
     842        Whether the lattice polytope has the IP property.
     843
     844        That is, the polytope is full-dimensional and the origin is a
     845        interior point not on the boundary.
     846
     847        OUTPUT:
     848
     849        Boolean.
     850
     851        EXAMPLES::
     852
     853            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     854            sage: LatticePolytope_PPL((-1,-1),(0,1),(1,0)).has_IP_property()
     855            True
     856            sage: LatticePolytope_PPL((-1,-1),(1,1)).has_IP_property()
     857            False
     858        """
     859        origin = C_Polyhedron(point(0*Variable(self.space_dimension())))
     860        is_included = Poly_Con_Relation.is_included()
     861        saturates = Poly_Con_Relation.saturates()
     862        for c in self.constraints():
     863            rel = origin.relation_with(c)
     864            if (not rel.implies(is_included)) or rel.implies(saturates):
     865                return False
     866        return True
     867
     868    @cached_method
     869    def restricted_automorphism_group(self, vertex_labels=None):
     870        r"""
     871        Return the restricted automorphism group.
     872
     873        First, let the linear automorphism group be the subgroup of
     874        the Euclidean group `E(d) = GL(d,\RR) \ltimes \RR^d`
     875        preserving the `d`-dimensional polyhedron. The Euclidean group
     876        acts in the usual way `\vec{x}\mapsto A\vec{x}+b` on the
     877        ambient space. The restricted automorphism group is the
     878        subgroup of the linear automorphism group generated by
     879        permutations of vertices. If the polytope is full-dimensional,
     880        it is equal to the full (unrestricted) automorphism group.
     881
     882        INPUT:
     883
     884        - ``vertex_labels`` -- a tuple or ``None`` (default). The
     885          labels of the vertices that will be used in the output
     886          permutation group. By default, the vertices are used
     887          themselves.
     888
     889        OUTPUT:
     890
     891        A
     892        :class:`PermutationGroup<sage.groups.perm_gps.permgroup.PermutationGroup_generic>`
     893        acting on the vertices (or the ``vertex_labels``, if
     894        specified).
     895
     896        REFERENCES:
     897
     898        ..  [BSS]
     899            David Bremner, Mathieu Dutour Sikiric, Achill Schuermann:
     900            Polyhedral representation conversion up to symmetries.
     901            http://arxiv.org/abs/math/0702239
     902
     903        EXAMPLES::
     904
     905            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     906            sage: Z3square = LatticePolytope_PPL((0,0), (1,2), (2,1), (3,3))
     907            sage: Z3square.restricted_automorphism_group(vertex_labels=(1,2,3,4))
     908            Permutation Group with generators [(2,3), (1,2)(3,4), (1,4)]
     909            sage: G = Z3square.restricted_automorphism_group(); G
     910            Permutation Group with generators [((1,2),(2,1)),
     911            ((0,0),(1,2))((2,1),(3,3)), ((0,0),(3,3))]
     912            sage: tuple(G.domain()) == Z3square.vertices()
     913            True
     914            sage: G.orbit(Z3square.vertices()[0])
     915            ((0, 0), (1, 2), (3, 3), (2, 1))
     916
     917            sage: cell24 = LatticePolytope_PPL(
     918            ...   (1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,-1,-1,1),(0,0,-1,1),
     919            ...   (0,-1,0,1),(-1,0,0,1),(1,0,0,-1),(0,1,0,-1),(0,0,1,-1),(-1,1,1,-1),
     920            ...   (1,-1,-1,0),(0,0,-1,0),(0,-1,0,0),(-1,0,0,0),(1,-1,0,0),(1,0,-1,0),
     921            ...   (0,1,1,-1),(-1,1,1,0),(-1,1,0,0),(-1,0,1,0),(0,-1,-1,1),(0,0,0,-1))
     922            sage: cell24.restricted_automorphism_group().cardinality()
     923            1152
     924        """
     925        if not self.is_full_dimensional():
     926            return self.affine_lattice_polytope().\
     927                restricted_automorphism_group(vertex_labels=vertex_labels)
     928        if vertex_labels is None:
     929            vertex_labels = self.vertices()
     930        from sage.groups.perm_gps.permgroup import PermutationGroup
     931        from sage.graphs.graph import Graph
     932        # good coordinates for the vertices
     933        v_list = []
     934        for v in self.minimized_generators():
     935            assert v.divisor().is_one()
     936            v_coords = (1,) + v.coefficients()
     937            v_list.append(vector(v_coords))
     938
     939        # Finally, construct the graph
     940        Qinv = sum( v.column() * v.row() for v in v_list ).inverse()
     941        G = Graph()
     942        for i in range(0,len(v_list)):
     943            for j in range(i+1,len(v_list)):
     944                v_i = v_list[i]
     945                v_j = v_list[j]
     946                G.add_edge(vertex_labels[i], vertex_labels[j], v_i * Qinv * v_j)
     947        return G.automorphism_group(edge_labels=True)
     948
     949    @cached_method
     950    def lattice_automorphism_group(self, points=None, point_labels=None):
     951        """
     952        The integral subgroup of the restricted automorphism group.
     953
     954        INPUT:
     955
     956        - ``points`` -- A tuple of coordinate vectors or ``None``
     957          (default). If specified, the points must form complete
     958          orbits under the lattice automorphism group. If ``None`` all
     959          vertices are used.
     960
     961        - ``point_labels`` -- A tuple of labels for the ``points`` or
     962          ``None`` (default). These will be used as labels for the do
     963          permutation group. If ``None`` the ``points`` will be used
     964          themselves.
     965
     966        OUTPUT:
     967
     968        The integral subgroup of the restricted automorphism group
     969        acting on the given ``points``, or all vertices if not
     970        specified.
     971
     972        EXAMPLES::
     973
     974            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     975            sage: Z3square = LatticePolytope_PPL((0,0), (1,2), (2,1), (3,3))
     976            sage: Z3square.lattice_automorphism_group()
     977            Permutation Group with generators [(), ((1,2),(2,1)),
     978            ((0,0),(3,3)), ((0,0),(3,3))((1,2),(2,1))]
     979
     980            sage: G1 = Z3square.lattice_automorphism_group(point_labels=(1,2,3,4));  G1
     981            Permutation Group with generators [(), (2,3), (1,4), (1,4)(2,3)]
     982            sage: G1.cardinality()
     983            4
     984
     985            sage: G2 = Z3square.restricted_automorphism_group(vertex_labels=(1,2,3,4)); G2
     986            Permutation Group with generators [(2,3), (1,2)(3,4), (1,4)]
     987            sage: G2.cardinality()
     988            8
     989
     990            sage: points = Z3square.integral_points();  points
     991            ((0, 0), (1, 1), (1, 2), (2, 1), (2, 2), (3, 3))
     992            sage: Z3square.lattice_automorphism_group(points, point_labels=(1,2,3,4,5,6))
     993            Permutation Group with generators [(), (3,4), (1,6)(2,5), (1,6)(2,5)(3,4)]
     994        """
     995        if not self.is_full_dimensional():
     996            return self.affine_lattice_polytope().lattice_automorphism_group()
     997
     998        if points is None:
     999            points = self.vertices()
     1000        if point_labels is None:
     1001            point_labels = tuple(points)
     1002        points = [ vector(ZZ, [1]+v.list()) for v in points ]
     1003        map(lambda x:x.set_immutable(), points)
     1004       
     1005        vertices = [ vector(ZZ, [1]+v.list()) for v in self.vertices() ]
     1006        pivots = matrix(ZZ, vertices).pivot_rows()
     1007        basis = matrix(ZZ, [ vertices[i] for i in pivots ])
     1008        Mat_ZZ = basis.parent()
     1009        basis_inverse = basis.inverse()
     1010
     1011        from sage.groups.perm_gps.permgroup import PermutationGroup
     1012        from sage.groups.perm_gps.permgroup_element import PermutationGroupElement
     1013        lattice_gens = []
     1014        G = self.restricted_automorphism_group(
     1015            vertex_labels=tuple(range(len(vertices))))
     1016        for g in G:
     1017            image = matrix(ZZ, [ vertices[g(i)] for i in pivots ])
     1018            m = basis_inverse*image
     1019            if m not in Mat_ZZ:
     1020                continue
     1021            perm_list = [ point_labels[points.index(p*m)]
     1022                          for p in points ]
     1023            lattice_gens.append(perm_list)
     1024        return PermutationGroup(lattice_gens, domain=point_labels)
     1025
     1026
     1027
     1028
  • sage/libs/ppl.pyx

    diff --git a/sage/libs/ppl.pyx b/sage/libs/ppl.pyx
    a b  
    148148#                  http://www.gnu.org/licenses/
    149149#*****************************************************************************
    150150
     151from sage.structure.sage_object cimport SageObject
    151152from sage.libs.gmp.mpz cimport mpz_t, mpz_set
    152153from sage.rings.integer cimport Integer
    153154from sage.rings.rational import Rational
     
    164165#  - solve linear program
    165166# These can only be triggered by methods in the Polyhedron class
    166167# they need to be wrapped in sig_on() / sig_off()
    167 
    168168####################################################
    169169cdef extern from "gmpxx.h":
    170170    cdef cppclass mpz_class:
     
    466466
    467467
    468468### Forward declarations ###########################
    469 cdef class _mutable_or_immutable(object)
     469cdef class _mutable_or_immutable(SageObject)
    470470cdef class Variable(object)
    471471cdef class Linear_Expression(object)
    472472cdef class Generator(object)
     
    486486####################################################
    487487### _mutable_or_immutable ##########################
    488488####################################################
    489 cdef class _mutable_or_immutable(object):
     489cdef class _mutable_or_immutable(SageObject):
    490490    r"""
    491491    A base class for mutable or immutable objects.
    492492
     
    11901190        raise NotImplementedError, 'The Polyhedron class is abstract, you must not instantiate it.'
    11911191
    11921192
    1193     def __repr__(self):
     1193    def _repr_(self):
    11941194        """
    11951195        Return a string representation.
    11961196
     
    12031203            sage: from sage.libs.ppl import Variable, C_Polyhedron
    12041204            sage: x = Variable(0)
    12051205            sage: y = Variable(1)
    1206             sage: C_Polyhedron( 5*x-2*y >=  x+y-1 ).__repr__()
     1206            sage: C_Polyhedron( 5*x-2*y >=  x+y-1 )._repr_()
    12071207            'A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 1 point, 1 ray, 1 line'
    12081208
    12091209        Special cases::
    12101210
    1211             sage: C_Polyhedron(3, 'empty').__repr__()
     1211            sage: C_Polyhedron(3, 'empty')._repr_()
    12121212            'The empty polyhedron in QQ^3'
    1213             sage: C_Polyhedron(3, 'universe').__repr__()
     1213            sage: C_Polyhedron(3, 'universe')._repr_()
    12141214            'The space-filling polyhedron in QQ^3'
    12151215        """
    12161216        dim = self.affine_dimension()
     
    51505150        """
    51515151        return self.thisptr.OK()
    51525152
     5153   
     5154    def __len__(self):
     5155        """
     5156        Return the number of generators in the system.
     5157
     5158            sage: from sage.libs.ppl import Variable, Generator_System, point
     5159            sage: x = Variable(0)
     5160            sage: y = Variable(1)
     5161            sage: gs = Generator_System()
     5162            sage: gs.insert(point(3*x+2*y))
     5163            sage: gs.insert(point(x))
     5164            sage: gs.insert(point(y))
     5165            sage: len(gs)
     5166            3
     5167        """
     5168        return sum([1 for g in self])
     5169
    51535170
    51545171    def __iter__(self):
    51555172        """
     
    61266143        """
    61276144        return self.thisptr.OK()
    61286145
     6146   
     6147    def __len__(self):
     6148        """
     6149        Return the number of constraints in the system.
     6150
     6151        EXAMPLES::
     6152
     6153            sage: from sage.libs.ppl import Variable, Constraint_System
     6154            sage: x = Variable(0)
     6155            sage: cs = Constraint_System( x>0 )
     6156            sage: cs.insert( x<1 )
     6157            sage: len(cs)
     6158            2
     6159        """
     6160        return sum([1 for c in self])
     6161
    61296162
    61306163    def __iter__(self):
    61316164        """