Ticket #12553: trac_12553_ppl_count_points.patch

File trac_12553_ppl_count_points.patch, 11.3 KB (added by vbraun, 10 years ago)

Rediffed patch

  • sage/geometry/integral_points.pyx

    # HG changeset patch
    # User Volker Braun <vbraun@stp.dias.ie>
    # Date 1341159365 -3600
    # Node ID 5f313ba8af428f45db2b59cbb9ef04caf40fa1f6
    # Parent  0f6f5a5dc03bb8e98bd3581806fefe0d9c9a24cd
    Trac #12553: Add interface for PALP polytope databases
    
    This patch implement counting of integral points (without enumerating them)
    
    diff --git a/sage/geometry/integral_points.pyx b/sage/geometry/integral_points.pyx
    a b  
    1414from sage.rings.all import QQ, RR, ZZ, gcd, lcm
    1515from sage.combinat.permutation import Permutation
    1616from sage.combinat.cartesian_product import CartesianProduct
    17 from sage.misc.all import prod
     17from sage.misc.all import prod, uniq
    1818import copy
    1919
    2020##############################################################################
     
    346346    - ``box_max`` -- A list of integers. The maximal value for each
    347347      coordinate of the rectangular bounding box.
    348348
    349     - ``polyhedron`` -- A :class:`~sage.geometry.polyhedra.Polyhedron`
    350       or ``None`` (default).
     349    - ``polyhedron`` -- A
     350      :class:`~sage.geometry.polyhedron.base.Polyhedron_base`, a PPL
     351      :class:`~sage.libs.ppl.C_Polyhedron`, or ``None`` (default).
    351352
    352353    OUTPUT:
    353354
     
    443444         (0, 100000000000000000000000000000000000000000000000001))
    444445        sage: len( rectangular_box_points([0,-100+10^50], [1,100+10^50], halfplane) )
    445446        201
     447
     448    Using a PPL polyhedron::
     449   
     450        sage: from sage.libs.ppl import Variable, Generator_System, C_Polyhedron, point
     451        sage: gs = Generator_System()
     452        sage: x = Variable(0); y = Variable(1); z = Variable(2)
     453        sage: gs.insert(point(0*x + 1*y + 0*z))
     454        sage: gs.insert(point(0*x + 1*y + 3*z))
     455        sage: gs.insert(point(3*x + 1*y + 0*z))
     456        sage: gs.insert(point(3*x + 1*y + 3*z))
     457        sage: poly = C_Polyhedron(gs)
     458        sage: rectangular_box_points([0]*3, [3]*3, poly)
     459        ((0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3),
     460         (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 1, 3), (3, 1, 0), (3, 1, 1), (3, 1, 2), (3, 1, 3))
    446461    """
    447462    assert len(box_min)==len(box_max)
    448463    cdef int d = len(box_min)
     
    551566    cdef object b
    552567    cdef object coeff
    553568    cdef object cache
     569    # The index of the inequality in the polyhedron H-representation
     570    cdef int index
    554571
    555     def __cinit__(self, A, b):
     572    def __cinit__(self, A, b, index=-1):
    556573        """
    557574        The Cython constructor
    558575
     
    570587        self.b = b
    571588        self.coeff = 0
    572589        self.cache = 0
     590        self.index = int(index)
    573591
    574592    def __repr__(self):
    575593        """
     
    614632    cdef bint is_not_satisfied(self, inner_loop_variable):
    615633        r"""
    616634        Test the inequality, using the cached value from :meth:`prepare_inner_loop`
     635
     636        OUTPUT:
     637
     638        Boolean. Whether the inequality is not satisfied.
    617639        """
    618640        return inner_loop_variable*self.coeff + self.cache < 0       
    619641
     642    cdef bint is_equality(Inequality_generic self, int inner_loop_variable):
     643        r"""
     644        Test the inequality, using the cached value from :meth:`prepare_inner_loop`
     645
     646        OUTPUT:
     647
     648        Boolean. Given the inequality `Ax+b\geq 0`, this method
     649        returns whether the equality `Ax+b=0` is satisfied.
     650        """
     651        return inner_loop_variable*self.coeff + self.cache == 0       
     652
    620653
    621654
    622655DEF INEQ_INT_MAX_DIM = 20
     
    677710    # the next-to-innermost coefficient
    678711    cdef int coeff_next
    679712    cdef int cache_next
     713    # The index of the inequality in the polyhedron H-representation
     714    cdef int index
    680715
    681716    cdef int _to_int(self, x) except? -999:
    682717        if not x in ZZ: raise ValueError('Not integral.')
    683718        return int(x)  # raises OverflowError in Cython if necessary
    684719
    685     def __cinit__(self, A, b, max_abs_coordinates):
     720    def __cinit__(self, A, b, max_abs_coordinates, index=-1):
    686721        """
    687722        The Cython constructor
    688723
     
    696731        """
    697732        cdef int i
    698733        self.dim = self._to_int(len(A))
     734        self.index = int(index)
    699735        if self.dim<1 or self.dim>INEQ_INT_MAX_DIM:
    700736            raise OverflowError('Dimension limit exceeded.')
    701737        for i in range(0,self.dim):
     
    752788    cdef bint is_not_satisfied(Inequality_int self, int inner_loop_variable):
    753789        return inner_loop_variable*self.coeff + self.cache < 0       
    754790
     791    cdef bint is_equality(Inequality_int self, int inner_loop_variable):
     792        return inner_loop_variable*self.coeff + self.cache == 0       
     793
    755794
    756795
    757796cdef class InequalityCollection:
     
    798837    """
    799838    cdef object ineqs_int
    800839    cdef object ineqs_generic
    801 
     840   
    802841    def __repr__(self):
    803842        r"""
    804843        Return a string representation.
     
    880919        self.ineqs_generic = []
    881920        if polyhedron is None:
    882921            return
     922
     923        try:
     924            # polyhedron is a PPL C_Polyhedron class?
     925            self._cinit_from_PPL(max_abs_coordinates, permutation, polyhedron)
     926        except AttributeError:
     927            try:
     928                # polyhedron is a Polyhedron class?
     929                self._cinit_from_Polyhedron(max_abs_coordinates, permutation, polyhedron)
     930            except AttributeError:
     931                raise TypeError('Cannot extract Hrepresentation data from polyhedron.')
     932
     933    cdef _cinit_from_PPL(self, max_abs_coordinates, permutation, polyhedron):
     934        """
     935        Initialize the inequalities from a PPL C_Polyhedron
     936
     937        See __cinit__ for a description of the arguments.
     938
     939        EXAMPLES::
     940       
     941            sage: from sage.libs.ppl import Variable, Generator_System, C_Polyhedron, point
     942            sage: gs = Generator_System()
     943            sage: x = Variable(0); y = Variable(1); z = Variable(2)
     944            sage: gs.insert(point(0*x + 0*y + 1*z))
     945            sage: gs.insert(point(0*x + 3*y + 1*z))
     946            sage: gs.insert(point(3*x + 0*y + 1*z))
     947            sage: gs.insert(point(3*x + 3*y + 1*z))
     948            sage: poly = C_Polyhedron(gs)
     949            sage: from sage.geometry.integral_points import InequalityCollection
     950            sage: InequalityCollection(poly, Permutation([1,3,2]), [0]*3, [3]*3 )
     951            The collection of inequalities
     952            integer: (0, 1, 0) x + -1 >= 0
     953            integer: (0, -1, 0) x + 1 >= 0
     954            integer: (1, 0, 0) x + 0 >= 0
     955            integer: (0, 0, 1) x + 0 >= 0
     956            integer: (-1, 0, 0) x + 3 >= 0
     957            integer: (0, 0, -1) x + 3 >= 0
     958        """
     959        for index,c in enumerate(polyhedron.minimized_constraints()):
     960            A = permutation.action(c.coefficients())
     961            b = c.inhomogeneous_term()
     962            try:
     963                H = Inequality_int(A, b, max_abs_coordinates, index)
     964                self.ineqs_int.append(H)
     965            except (OverflowError, ValueError):
     966                H = Inequality_generic(A, b, index)
     967                self.ineqs_generic.append(H)               
     968            if c.is_equality():
     969                A = [ -a for a in A ]
     970                b = -b
     971                try:
     972                    H = Inequality_int(A, b, max_abs_coordinates, index)
     973                    self.ineqs_int.append(H)
     974                except (OverflowError, ValueError):
     975                    H = Inequality_generic(A, b, index)
     976                    self.ineqs_generic.append(H)
     977
     978    cdef _cinit_from_Polyhedron(self, max_abs_coordinates, permutation, polyhedron):
     979        """
     980        Initialize the inequalities from a Sage Polyhedron
     981
     982        See __cinit__ for a description of the arguments.
     983
     984        EXAMPLES::
     985       
     986            sage: from sage.geometry.integral_points import InequalityCollection
     987            sage: line = Polyhedron(eqns=[(2,3,7)])
     988            sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 )
     989            The collection of inequalities
     990            integer: (3, 7) x + 2 >= 0
     991            integer: (-3, -7) x + -2 >= 0
     992        """
    883993        for Hrep_obj in polyhedron.inequality_generator():
    884994            A, b = self._make_A_b(Hrep_obj, permutation)
    885995            try:
    886                 H = Inequality_int(A, b, max_abs_coordinates)
     996                H = Inequality_int(A, b, max_abs_coordinates, Hrep_obj.index())
    887997                self.ineqs_int.append(H)
    888998            except (OverflowError, ValueError):
    889                 H = Inequality_generic(A, b)
     999                H = Inequality_generic(A, b, Hrep_obj.index())
    8901000                self.ineqs_generic.append(H)
    8911001        for Hrep_obj in polyhedron.equation_generator():
    8921002            A, b = self._make_A_b(Hrep_obj, permutation)
    8931003            # add inequality
    8941004            try:
    895                 H = Inequality_int(A, b, max_abs_coordinates)
     1005                H = Inequality_int(A, b, max_abs_coordinates, Hrep_obj.index())
    8961006                self.ineqs_int.append(H)
    8971007            except (OverflowError, ValueError):
    898                 H = Inequality_generic(A, b)
     1008                H = Inequality_generic(A, b, Hrep_obj.index())
    8991009                self.ineqs_generic.append(H)
    9001010            # add sign-reversed inequality
    9011011            A = [ -a for a in A ]
    9021012            b = -b
    9031013            try:
    904                 H = Inequality_int(A, b, max_abs_coordinates)
     1014                H = Inequality_int(A, b, max_abs_coordinates, Hrep_obj.index())
    9051015                self.ineqs_int.append(H)
    9061016            except (OverflowError, ValueError):
    907                 H = Inequality_generic(A, b)
     1017                H = Inequality_generic(A, b, Hrep_obj.index())
    9081018                self.ineqs_generic.append(H)
    9091019       
    9101020    def prepare_next_to_inner_loop(self, p):
     
    10541164                return False
    10551165        return True
    10561166   
     1167    def satisfied_as_equalities(self, inner_loop_variable):
     1168        """
     1169        Return the inequalities (by their index) that are satisfied as
     1170        equalities.
    10571171
     1172        INPUT:
     1173
     1174        - ``inner_loop_variable`` -- Integer. the 0-th coordinate of
     1175          the lattice point.
     1176
     1177        OUTPUT:
     1178       
     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. 
     1183       
     1184        EXAMPLES::
     1185       
     1186            sage: from sage.geometry.integral_points import InequalityCollection
     1187            sage: quadrant = Polyhedron(rays=[(1,0), (0,1)])
     1188            sage: ieqs = InequalityCollection(quadrant, Permutation([1,2]), [-1]*2, [1]*2 )
     1189            sage: ieqs.prepare_next_to_inner_loop([-1,0])
     1190            sage: ieqs.prepare_inner_loop([-1,0])
     1191            sage: ieqs.satisfied_as_equalities(-1)
     1192            (1,)
     1193            sage: ieqs.satisfied_as_equalities(0)
     1194            (0, 1)
     1195            sage: ieqs.satisfied_as_equalities(1)
     1196            (1,)
     1197        """
     1198        cdef int i
     1199        result = []
     1200        for i in range(0,len(self.ineqs_int)):
     1201            ineq = self.ineqs_int[i]
     1202            if (<Inequality_int>ineq).is_equality(inner_loop_variable):
     1203                result.append( (<Inequality_int>ineq).index )
     1204        for i in range(0,len(self.ineqs_generic)):
     1205            ineq = self.ineqs_generic[i]
     1206            if (<Inequality_generic>ineq).is_equality(inner_loop_variable):
     1207                result.append( (<Inequality_generic>ineq).index )
     1208        return tuple(uniq(result))
     1209   
    10581210
    10591211
    10601212cpdef print_cache(InequalityCollection inequality_collection):