# Ticket #11429: trac_11429_cythonize_lattice_points.patch

File trac_11429_cythonize_lattice_points.patch, 52.0 KB (added by vbraun, 11 years ago)

Updated patch

• ## module_list.py

# HG changeset patch
# User Volker Braun <vbraun@stp.dias.ie>
# Date 1323988670 0
# Node ID 34cc23b585e0eb2958a554f928cd32aed04df6be
Trac #11429: Count integral points without PALP

This patch implements the basic algorithms in Cython
and optimizes the heck out of them. Its fast but sometimes
not easy to read. No functionality visible to the
end-user is added, the purpose of this patch is purely
speed.

diff --git a/module_list.py b/module_list.py
 a sources = ['sage/geometry/toric_lattice_element.pyx'], libraries=['gmp']), Extension('sage.geometry.integral_points', sources = ['sage/geometry/integral_points.pyx']), Extension('sage.geometry.triangulation.base', sources = ['sage/geometry/triangulation/functions.cc', 'sage/geometry/triangulation/data.cc',
• ## sage/geometry/cone.py

diff --git a/sage/geometry/cone.py b/sage/geometry/cone.py
 a from sage.libs.ppl import C_Polyhedron, Generator_System, Constraint_System, \ Linear_Expression, ray as PPL_ray, point as PPL_point, \ Poly_Con_Relation from sage.combinat.cartesian_product import CartesianProduct from sage.geometry.integral_points import parallelotope_points def is_Cone(x): gens.update(cone.semigroup_generators()) return tuple(gens) if self.is_trivial(): return tuple() gens = list(parallelotope_points(self.rays())) + list(self.rays()) gens = list(parallelotope_points(self.rays(), N)) + list(self.rays()) gens = filter(lambda v:gcd(v)==1, gens) return tuple(gens) return vector(ZZ, p.get_values(x)) ################################################################# def parallelotope_points(spanning_points): r""" Return integral points in the parallelotope spanned by a the rays. See :meth:~ConvexRationalPolyhedralCone.semigroup_generators for a description of the algorithm. INPUT: - spanning_points -- a non-empty list of linearly independent rays (\ZZ-vectors or :class:toric lattice  elements), not necessarily primitive lattice points. OUTPUT: The tuple of all lattice points in the half-open parallelotope spanned by the rays r_i, .. math:: \mathop{par}(\{r_i\}) = \sum_{0\leq a_i < 1} a_i r_i By half-open parallelotope, we mean that the points in the facets not meeting the origin are omitted. EXAMPLES: Note how the points on the outward-facing factes are omitted:: sage: from sage.geometry.cone import parallelotope_points sage: rays = map(vector, [(2,0), (0,2)]) sage: parallelotope_points(rays) ((0, 0), (1, 0), (0, 1), (1, 1)) The rays can also be toric lattice points:: sage: rays = map(ToricLattice(2), [(2,0), (0,2)]) sage: parallelotope_points(rays) (N(0, 0), N(1, 0), N(0, 1), N(1, 1)) A non-smooth cone:: sage: c = Cone([ (1,0), (1,2) ]) sage: parallelotope_points(c.rays()) (N(0, 0), N(1, 1)) A ValueError is raised if the spanning_points are not linearly independent:: sage: rays = map(ToricLattice(2), [(1,1)]*2) sage: parallelotope_points(rays) Traceback (most recent call last): ... ValueError: The spanning points are not linearly independent! """ N = spanning_points[0].parent()  # this is why there needs to be at least one ray # compute points in the open parallelotope, see [BrunsKoch] R = matrix(spanning_points).transpose() D,U,V = R.smith_form() e = D.diagonal()          # the elementary divisors d = prod(e)               # the determinant if d==0: raise ValueError('The spanning points are not linearly independent!') u = U.inverse().columns() # generators for gp(semigroup) gens = [] # "inverse" of the ray matrix as far as possible over ZZ # R*Rinv == diagonal_matrix([d]*D.ncols() + [0]*(D.nrows()-D.ncols())) # If R is full rank, this is Rinv = matrix(ZZ, R.inverse() * d) Dinv = D.transpose() for i in range(0,D.ncols()): Dinv[i,i] = d/D[i,i] Rinv = V * Dinv * U for b in CartesianProduct(*[ range(0,i) for i in e ]): # this is our generator modulo the lattice spanned by the rays gen_mod_rays = sum( b_i*u_i for b_i, u_i in zip(b,u) ) q_times_d = Rinv * gen_mod_rays q_times_d = vector(ZZ,[ q_i % d  for q_i in q_times_d ]) gen = N(R*q_times_d / d) gen.set_immutable() gens.append(gen) assert(len(gens) == d) return tuple(gens)
• ## new file sage/geometry/integral_points.pyx

diff --git a/sage/geometry/integral_points.pyx b/sage/geometry/integral_points.pyx
new file mode 100644
 - # cython: profile=True r""" Cython helper methods to compute integral points in polyhedra. """ #***************************************************************************** #       Copyright (C) 2010 Volker Braun # #  Distributed under the terms of the GNU General Public License (GPL) # #                  http://www.gnu.org/licenses/ #***************************************************************************** from sage.matrix.constructor import matrix, column_matrix, vector, diagonal_matrix from sage.rings.all import QQ, RR, ZZ, gcd, lcm from sage.combinat.permutation import Permutation from sage.combinat.cartesian_product import CartesianProduct from sage.misc.all import prod import copy ############################################################################## # The basic idea to enumerate the lattice points in the parallelotope # is to use the Smith normal form of the ray coordinate matrix to # bring the lattice into a "nice" basis. Here is the straightforward # implementation. Note that you do not need to reduce to the # full-dimensional case, the Smith normal form takes care of that for # you. # ## def parallelotope_points(spanning_points, lattice): ##     # compute points in the open parallelotope, see [BrunsKoch] ##     R = matrix(spanning_points).transpose() ##     D,U,V = R.smith_form() ##     e = D.diagonal()          # the elementary divisors ##     d = prod(e)               # the determinant ##     u = U.inverse().columns() # generators for gp(semigroup) ## ##     # "inverse" of the ray matrix as far as possible over ZZ ##     # R*Rinv == diagonal_matrix([d]*D.ncols() + [0]*(D.nrows()-D.ncols())) ##     # If R is full rank, this is Rinv = matrix(ZZ, R.inverse() * d) ##     Dinv = D.transpose() ##     for i in range(0,D.ncols()): ##         Dinv[i,i] = d/D[i,i] ##     Rinv = V * Dinv * U ## ##     gens = [] ##     for b in CartesianProduct(*[ range(0,i) for i in e ]): ##         # this is our generator modulo the lattice spanned by the rays ##         gen_mod_rays = sum( b_i*u_i for b_i, u_i in zip(b,u) ) ##         q_times_d = Rinv * gen_mod_rays ##         q_times_d = vector(ZZ,[ q_i % d  for q_i in q_times_d ]) ##         gen = lattice(R*q_times_d / d) ##         gen.set_immutable() ##         gens.append(gen) ##     assert(len(gens) == d) ##     return tuple(gens) # # The problem with the naive implementation is that it is slow: # #   1. You can simplify some of the matrix multiplications # #   2. The inner loop keeps creating new matrices and vectors, which #      is slow. It needs to recycle objects. Instead of creating a new #      lattice point again and again, change the entries of an #      existing lattice point and then copy it! def parallelotope_points(spanning_points, lattice): r""" Return integral points in the parallelotope starting at the origin and spanned by the spanning_points. See :meth:~ConvexRationalPolyhedralCone.semigroup_generators for a description of the algorithm. INPUT: - spanning_points -- a non-empty list of linearly independent rays (\ZZ-vectors or :class:toric lattice  elements), not necessarily primitive lattice points. OUTPUT: The tuple of all lattice points in the half-open parallelotope spanned by the rays r_i, .. math:: \mathop{par}(\{r_i\}) = \sum_{0\leq a_i < 1} a_i r_i By half-open parallelotope, we mean that the points in the facets not meeting the origin are omitted. EXAMPLES: Note how the points on the outward-facing factes are omitted:: sage: from sage.geometry.integral_points import parallelotope_points sage: rays = map(vector, [(2,0), (0,2)]) sage: parallelotope_points(rays, ZZ^2) ((0, 0), (1, 0), (0, 1), (1, 1)) The rays can also be toric lattice points:: sage: rays = map(ToricLattice(2), [(2,0), (0,2)]) sage: parallelotope_points(rays, ToricLattice(2)) (N(0, 0), N(1, 0), N(0, 1), N(1, 1)) A non-smooth cone:: sage: c = Cone([ (1,0), (1,2) ]) sage: parallelotope_points(c.rays(), c.lattice()) (N(0, 0), N(1, 1)) A ValueError is raised if the spanning_points are not linearly independent:: sage: rays = map(ToricLattice(2), [(1,1)]*2) sage: parallelotope_points(rays, ToricLattice(2)) Traceback (most recent call last): ... ValueError: The spanning points are not linearly independent! TESTS:: sage: rays = map(vector,[(-3, -2, -3, -2), (-2, -1, -8, 5), (1, 9, -7, -4), (-3, -1, -2, 2)]) sage: len(parallelotope_points(rays, ZZ^4)) 967 """ R = matrix(spanning_points).transpose() e, d, VDinv = ray_matrix_normal_form(R) points = loop_over_parallelotope_points(e, d, VDinv, R, lattice) for p in points: p.set_immutable() return points def ray_matrix_normal_form(R): r""" Compute the Smith normal form of the ray matrix for :func:parallelotope_points. INPUT: - R -- \ZZ-matrix whose columns are the rays spanning the parallelotope. OUTPUT: A tuple containing e, d, and VDinv. EXAMPLES:: sage: from sage.geometry.integral_points import ray_matrix_normal_form sage: R = column_matrix(ZZ,[3,3,3]) sage: ray_matrix_normal_form(R) ([3], 3, [1]) """ D,U,V = R.smith_form() e = D.diagonal()            # the elementary divisors d = prod(e)                 # the determinant if d==0: raise ValueError('The spanning points are not linearly independent!') Dinv = diagonal_matrix(ZZ,[ d/D[i,i] for i in range(0,D.ncols()) ]) VDinv = V * Dinv return (e,d,VDinv) # The optimized version avoids constructing new matrices, vectors, and lattice points cpdef loop_over_parallelotope_points(e, d, VDinv, R, lattice, A=None, b=None): r""" The inner loop of :func:parallelotope_points. INPUT: See :meth:parallelotope_points for e, d, VDinv, R, lattice. - A, b: Either both None or a vector and number. If present, only the parallelotope points satisfying A x \leq b are returned. OUTPUT: The points of the half-open parallelotope as a tuple of lattice points. EXAMPLES: sage: e = [3] sage: d = prod(e) sage: VDinv = matrix(ZZ, [[1]]) sage: R = column_matrix(ZZ,[3,3,3]) sage: lattice = ZZ^3 sage: from sage.geometry.integral_points import loop_over_parallelotope_points sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice) ((0, 0, 0), (1, 1, 1), (2, 2, 2)) sage: A = vector(ZZ, [1,0,0]) sage: b = 1 sage: loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b) ((0, 0, 0), (1, 1, 1)) """ cdef int i, j cdef int dim = VDinv.nrows() cdef int ambient_dim = R.nrows() gens = [] s = ZZ.zero() # summation variable gen = lattice(0) q_times_d = vector(ZZ, dim) for base in CartesianProduct(*[ range(0,i) for i in e ]): for i in range(0, dim): s = ZZ.zero() for j in range(0, dim): s += VDinv[i,j] * base[j] q_times_d[i] = s % d for i in range(0, ambient_dim): s = ZZ.zero() for j in range(0, dim): s += R[i,j] * q_times_d[j] gen[i] = s / d if A is not None: s = ZZ.zero() for i in range(0, ambient_dim): s += A[i] * gen[i] if s>b: continue gens.append(copy.copy(gen)) if A is None: assert(len(gens) == d) return tuple(gens) ############################################################################## def simplex_points(vertices): r""" Return the integral points in a lattice simplex. INPUT: - vertices -- an iterable of integer coordinate vectors. The indices of vertices that span the simplex under consideration. OUTPUT: A tuple containing the integral point coordinates as \ZZ-vectors. EXAMPLES:: sage: from sage.geometry.integral_points import simplex_points sage: simplex_points([(1,2,3), (2,3,7), (-2,-3,-11)]) ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7)) The simplex need not be full-dimensional:: sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)]) sage: simplex_points(simplex.Vrepresentation()) ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5)) sage: simplex_points([(2,3,7)]) ((2, 3, 7),) TESTS:: sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)] sage: simplex = Polyhedron(v); simplex A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices. sage: pts = simplex_points(simplex.Vrepresentation()) sage: len(pts) 49 sage: for p in pts: p.set_immutable() sage: len(set(pts)) 49 sage: all(simplex.contains(p) for p in pts) True sage: v = [(4,-1,-1,-1), (-1,4,-1,-1), (-1,-1,4,-1), (-1,-1,-1,4), (-1,-1,-1,-1)] sage: P4mirror = Polyhedron(v); P4mirror A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices. sage: len( simplex_points(P4mirror.Vrepresentation()) ) 126 """ rays = [ vector(ZZ, v) for v in vertices ] if len(rays)==0: return tuple() origin = rays.pop() origin.set_immutable() if len(rays)==0: return tuple([origin]) translate_points(rays, origin) # Find equation Ax<=b that cuts out simplex from parallelotope Rt = matrix(rays) R = Rt.transpose() if R.is_square(): b = abs(R.det()) A = R.solve_left(vector([b]*len(rays))) else: RtR = Rt * R b = abs(RtR.det()) A = RtR.solve_left(vector([b]*len(rays))) * Rt # e, d, VDinv = ray_matrix_normal_form(R) #    print origin #    print rays #    print parallelotope_points(rays, origin.parent()) #    print 'A = ', A #    print 'b = ', b e, d, VDinv = ray_matrix_normal_form(R) lattice = origin.parent() points = loop_over_parallelotope_points(e, d, VDinv, R, lattice, A, b) + tuple(rays) translate_points(points, -origin) return points cdef translate_points(v_list, delta): r""" Add delta to each vector in v_list. """ cdef int dim = delta.degree() for v in v_list: for i in range(0,dim): v[i] -= delta[i] ############################################################################## # For points with "small" coordinates (that is, fitting into a small # rectangular bounding box) it is faster to naively enumerate the # points. This saves the overhead of triangulating the polytope etc. def rectangular_box_points(box_min, box_max, polyhedron=None): r""" Return the integral points in the lattice bounding box that are also contained in the given polyhedron. INPUT: - box_min -- A list of integers. The minimal value for each coordinate of the rectangular bounding box. - box_max -- A list of integers. The maximal value for each coordinate of the rectangular bounding box. - polyhedron -- A :class:~sage.geometry.polyhedra.Polyhedron or None (default). OUTPUT: A tuple containing the integral points of the rectangular box spanned by box_min and box_max and that lie inside the polyhedron. For sufficiently large bounding boxes, this are all integral points of the polyhedron. If no polyhedron is specified, all integral points of the rectangular box are returned. ALGORITHM: This function implements the naive algorithm towards counting integral points. Given min and max of vertex coordinates, it iterates over all points in the bounding box and checks whether they lie in the polyhedron. The following optimizations are implemented: * Cython: Use machine integers and optimizing C/C++ compiler where possible, arbitrary precision integers where necessary. Bounds checking, no compile time limits. * Unwind inner loop (and next-to-inner loop): .. math:: Ax\leq b \quad \Leftrightarrow \quad a_1 x_1 ~\leq~ b - \sum_{i=2}^d a_i x_i so we only have to evaluate a_1 * x_1 in the inner loop. * Coordinates are permuted to make the longest box edge the inner loop. The inner loop is optimized to run very fast, so its best to do as much work as possible there. * Continuously reorder inequalities and test the most restrictive inequalities first. * Use convexity and only find first and last allowed point in the inner loop. The points in-between must be points of the polyhedron, too. EXAMPLES:: sage: from sage.geometry.integral_points import rectangular_box_points sage: rectangular_box_points([0,0,0],[1,2,3]) ((0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3), (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 1, 3), (0, 2, 0), (0, 2, 1), (0, 2, 2), (0, 2, 3), (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 0, 3), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 1, 3), (1, 2, 0), (1, 2, 1), (1, 2, 2), (1, 2, 3)) sage: cell24 = polytopes.twenty_four_cell() sage: rectangular_box_points([-1]*4, [1]*4, cell24) ((-1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0), (0, 0, 0, -1), (0, 0, 0, 0), (0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0)) sage: d = 3 sage: dilated_cell24 = d*cell24 sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) ) 305 sage: d = 6 sage: dilated_cell24 = d*cell24 sage: len( rectangular_box_points([-d]*4, [d]*4, dilated_cell24) ) 3625 sage: polytope = Polyhedron([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4)]) sage: pts = rectangular_box_points([-4]*4, [4]*4, polytope); pts ((-4, -3, -2, -1), (-1, 0, 0, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 1, 3, 0), (1, 2, 1, 1), (1, 2, 2, 2), (1, 3, 2, 4), (2, 1, 1, 1), (3, 1, 1, 1)) sage: all(polytope.contains(p) for p in pts) True sage: set(map(tuple,pts)) == \ ...   set([(-4,-3,-2,-1),(3,1,1,1),(1,2,1,1),(1,1,3,0),(1,3,2,4), ...        (0,1,1,1),(1,2,2,2),(-1,0,0,1),(1,1,1,1),(2,1,1,1)])   # computed with PALP True Long ints and non-integral polyhedra are explictly allowed:: sage: polytope = Polyhedron([[1], [10*pi.n()]], field=RDF) sage: len( rectangular_box_points([-100], [100], polytope) ) 31 sage: halfplane = Polyhedron(ieqs=[(-1,1,0)]) sage: rectangular_box_points([0,-1+10^50], [0,1+10^50]) ((0, 99999999999999999999999999999999999999999999999999), (0, 100000000000000000000000000000000000000000000000000), (0, 100000000000000000000000000000000000000000000000001)) sage: len( rectangular_box_points([0,-100+10^50], [1,100+10^50], halfplane) ) 201 """ assert len(box_min)==len(box_max) cdef int d = len(box_min) diameter = sorted([ (box_max[i]-box_min[i], i) for i in range(0,d) ], reverse=True) diameter_value = [ x[0] for x in diameter ] diameter_index = [ x[1] for x in diameter ] sort_perm = Permutation([ i+1 for i in diameter_index]) orig_perm = sort_perm.inverse() orig_perm_list = [ i-1 for i in orig_perm ] box_min = sort_perm.action(box_min) box_max = sort_perm.action(box_max) inequalities = InequalityCollection(polyhedron, sort_perm, box_min, box_max) v = vector(ZZ, d) points = [] for p in loop_over_rectangular_box_points(box_min, box_max, inequalities, d): #  v = vector(ZZ, orig_perm.action(p))   # too slow for i in range(0,d): v[i] = p[orig_perm_list[i]] points.append(copy.copy(v)) return tuple(points) cdef loop_over_rectangular_box_points(box_min, box_max, inequalities, int d): """ The inner loop of :func:rectangular_box_points. INPUT: - box_min, box_max -- the bounding box. - inequalities -- a :class:InequalityCollection containing the inequalities defining the polyhedron. - d -- the ambient space dimension. OUTPUT: The integral points in the bounding box satisfying all inequalities. """ cdef int inc points = [] p = copy.copy(box_min) inequalities.prepare_next_to_inner_loop(p) while True: inequalities.prepare_inner_loop(p) i_min = box_min[0] i_max = box_max[0] # Find the lower bound for the allowed region while i_min <= i_max: if inequalities.are_satisfied(i_min): break i_min += 1 # Find the upper bound for the allowed region while i_min <= i_max: if inequalities.are_satisfied(i_max): break i_max -= 1 # The points i_min .. i_max are contained in the polyhedron i = i_min while i <= i_max: p[0] = i points.append(tuple(p)) i += 1 # finally increment the other entries in p to move on to next inner loop inc = 1 if d==1: return points while True: if p[inc]==box_max[inc]: p[inc] = box_min[inc] inc += 1 if inc==d: return points else: p[inc] += 1 break if inc>1: inequalities.prepare_next_to_inner_loop(p) cdef class Inequality_generic: """ An inequality whose coefficients are arbitrary Python/Sage objects INPUT: - A -- list of integers. - b -- integer OUTPUT: Inequality A x + b \geq 0. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_generic sage: Inequality_generic([2*pi,sqrt(3),7/2], -5.5) generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0 """ cdef object A cdef object b cdef object coeff cdef object cache def __cinit__(self, A, b): """ The Cython constructor INPUT: See :class:Inequality_generic. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_generic sage: Inequality_generic([2*pi,sqrt(3),7/2], -5.5) generic: (2*pi, sqrt(3), 7/2) x + -5.50000000000000 >= 0 """ self.A = A self.b = b self.coeff = 0 self.cache = 0 def __repr__(self): """ Return a string representation. OUTPUT: String. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_generic sage: Inequality_generic([2,3,7], -5).__repr__() 'generic: (2, 3, 7) x + -5 >= 0' """ s = 'generic: (' s += ', '.join([str(self.A[i]) for i in range(0,len(self.A))]) s += ') x + ' + str(self.b) + ' >= 0' return s cdef prepare_next_to_inner_loop(self, p): """ In :class:Inequality_int this method is used to peel of the next-to-inner loop. See :meth:InequalityCollection.prepare_inner_loop for more details. """ pass cdef prepare_inner_loop(self, p): """ Peel off the inner loop. See :meth:InequalityCollection.prepare_inner_loop for more details. """ cdef int j self.coeff = self.A[0] self.cache = self.b for j in range(1,len(self.A)): self.cache += self.A[j] * p[j] cdef bint is_not_satisfied(self, inner_loop_variable): r""" Test the inequality, using the cached value from :meth:prepare_inner_loop """ return inner_loop_variable*self.coeff + self.cache < 0 DEF INEQ_INT_MAX_DIM = 20 cdef class Inequality_int: """ Fast version of inequality in the case that all coefficient fit into machine ints. INPUT: - A -- list of integers. - b -- integer - max_abs_coordinates -- the maximum of the coordinates that one wants to evalate the coordinates on. Used for overflow checking. OUTPUT: Inequality A x + b \geq 0. A OverflowError is raised if a machine integer is not long enough to hold the results. A ValueError is raised if some of the input is not integral. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_int sage: Inequality_int([2,3,7], -5, [10]*3) integer: (2, 3, 7) x + -5 >= 0 sage: Inequality_int([1]*21, -5, [10]*21) Traceback (most recent call last): ... OverflowError: Dimension limit exceeded. sage: Inequality_int([2,3/2,7], -5, [10]*3) Traceback (most recent call last): ... ValueError: Not integral. sage: Inequality_int([2,3,7], -5.2, [10]*3) Traceback (most recent call last): ... ValueError: Not integral. sage: Inequality_int([2,3,7], -5*10^50, [10]*3)  # actual error message can differ between 32 and 64 bit Traceback (most recent call last): ... OverflowError: ... """ cdef int A[INEQ_INT_MAX_DIM] cdef int b cdef int dim # the innermost coefficient cdef int coeff cdef int cache # the next-to-innermost coefficient cdef int coeff_next cdef int cache_next cdef int _to_int(self, x) except? -999: if not x in ZZ: raise ValueError('Not integral.') return int(x)  # raises OverflowError in Cython if necessary def __cinit__(self, A, b, max_abs_coordinates): """ The Cython constructor See :class:Inequality_int for input. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_int sage: Inequality_int([2,3,7], -5, [10]*3) integer: (2, 3, 7) x + -5 >= 0 """ cdef int i self.dim = self._to_int(len(A)) if self.dim<1 or self.dim>INEQ_INT_MAX_DIM: raise OverflowError('Dimension limit exceeded.') for i in range(0,self.dim): self.A[i] = self._to_int(A[i]) self.b = self._to_int(b) self.coeff = self.A[0] if self.dim>0: self.coeff_next = self.A[1] # finally, make sure that there cannot be any overflow during the enumeration self._to_int(sum( [ZZ(b)]+[ZZ(A[i])*ZZ(max_abs_coordinates[i]) for i in range(0,self.dim)] )) def __repr__(self): """ Return a string representation. OUTPUT: String. EXAMPLES:: sage: from sage.geometry.integral_points import Inequality_int sage: Inequality_int([2,3,7], -5, [10]*3).__repr__() 'integer: (2, 3, 7) x + -5 >= 0' """ s = 'integer: (' s += ', '.join([str(self.A[i]) for i in range(0,self.dim)]) s += ') x + ' + str(self.b) + ' >= 0' return s cdef prepare_next_to_inner_loop(Inequality_int self, p): """ Peel off the next-to-inner loop. See :meth:InequalityCollection.prepare_inner_loop for more details. """ cdef int j self.cache_next = self.b for j in range(2,self.dim): self.cache_next += self.A[j] * p[j] cdef prepare_inner_loop(Inequality_int self, p): """ Peel off the inner loop. See :meth:InequalityCollection.prepare_inner_loop for more details. """ cdef int j if self.dim>1: self.cache = self.cache_next + self.coeff_next*p[1] else: self.cache = self.cache_next cdef bint is_not_satisfied(Inequality_int self, int inner_loop_variable): return inner_loop_variable*self.coeff + self.cache < 0 cdef class InequalityCollection: """ A collection of inequalities. INPUT: - polyhedron -- a polyhedron defining the inequalities. - permutation -- a permution of the coordinates. Will be used to permute the coordinates of the inequality. - box_min, box_max -- the (not permuted) minimal and maximal coordinates of the bounding box. Used for bounds checking. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], field=QQ) sage: ieq = InequalityCollection(P_QQ, Permutation([1,2,3]), [0]*3,[1]*3); ieq The collection of inequalities integer: (-1, -1, -1) x + 1 >= 0 integer: (-1, -1, 4) x + 1 >= 0 integer: (-1, 4, -1) x + 1 >= 0 integer: (3, -2, -2) x + 2 >= 0 sage: P_RR = Polyhedron(identity_matrix(2).columns() + [(-2.7, -1)], field=RDF) sage: InequalityCollection(P_RR, Permutation([1,2]), [0]*2, [1]*2) The collection of inequalities integer: (-1, -1) x + 1 >= 0 generic: (-1.0, 3.7) x + 1.0 >= 0 generic: (1.0, -1.35) x + 1.35 >= 0 sage: line = Polyhedron(eqns=[(2,3,7)]) sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ) The collection of inequalities integer: (3, 7) x + 2 >= 0 integer: (-3, -7) x + -2 >= 0 TESTS:: sage: TestSuite(ieq).run(skip='_test_pickling') """ cdef object ineqs_int cdef object ineqs_generic def __repr__(self): r""" Return a string representation. OUTPUT: String. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: line = Polyhedron(eqns=[(2,3,7)]) sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ).__repr__() 'The collection of inequalities\ninteger: (3, 7) x + 2 >= 0\ninteger: (-3, -7) x + -2 >= 0' """ s = 'The collection of inequalities\n' for ineq in self.ineqs_int: s += str(ineq) + '\n' for ineq in self.ineqs_generic: s += str(ineq) + '\n' return s.strip() def _make_A_b(self, Hrep_obj, permutation): r""" Return the coefficients and constant of the H-representation object. INPUT: - Hrep_obj -- a H-representation object of the polyhedron. - permutation -- the permutation of the coordinates to apply to A. OUTPUT: A pair (A,b). EXAXMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: line = Polyhedron(eqns=[(2,3,7)]) sage: ieq = InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ) sage: ieq._make_A_b(line.Hrepresentation(0), Permutation([1,2])) ([3, 7], 2) sage: ieq._make_A_b(line.Hrepresentation(0), Permutation([2,1])) ([7, 3], 2) """ v = list(Hrep_obj) A = permutation.action(v[1:]) b = v[0] try: x = lcm([a.denominator() for a in A] + [b.denominator()]) A = [ a*x for a in A ] b = b*x except AttributeError: pass return (A,b) def __cinit__(self, polyhedron, permutation, box_min, box_max): """ The Cython constructor See the class documentation for the desrciption of the arguments. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: line = Polyhedron(eqns=[(2,3,7)]) sage: InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ) The collection of inequalities integer: (3, 7) x + 2 >= 0 integer: (-3, -7) x + -2 >= 0 """ max_abs_coordinates = [ max(abs(c_min), abs(c_max)) for c_min, c_max in zip(box_min, box_max) ] max_abs_coordinates = permutation.action(max_abs_coordinates) self.ineqs_int = [] self.ineqs_generic = [] if not polyhedron: return for Hrep_obj in polyhedron.inequality_generator(): A, b = self._make_A_b(Hrep_obj, permutation) try: H = Inequality_int(A, b, max_abs_coordinates) self.ineqs_int.append(H) except (OverflowError, ValueError): H = Inequality_generic(A, b) self.ineqs_generic.append(H) for Hrep_obj in polyhedron.equation_generator(): A, b = self._make_A_b(Hrep_obj, permutation) # add inequality try: H = Inequality_int(A, b, max_abs_coordinates) self.ineqs_int.append(H) except (OverflowError, ValueError): H = Inequality_generic(A, b) self.ineqs_generic.append(H) # add sign-reversed inequality A = [ -a for a in A ] b = -b try: H = Inequality_int(A, b, max_abs_coordinates) self.ineqs_int.append(H) except (OverflowError, ValueError): H = Inequality_generic(A, b) self.ineqs_generic.append(H) def prepare_next_to_inner_loop(self, p): r""" Peel off the next-to-inner loop. In the next-to-inner loop of :func:rectangular_box_points, we have to repeatedly evaluate A x-A_0 x_0+b. To speed up computation, we pre-evaluate .. math:: c = b + sum_{i=2} A_i x_i and only compute A x-A_0 x_0+b = A_1 x_1 +c \geq 0 in the next-to-inner loop. INPUT: - p -- the point coordinates. Only p[2:] coordinates are potentially used by this method. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection, print_cache sage: P = Polyhedron(ieqs=[(2,3,7,11)]) sage: ieq = InequalityCollection(P, Permutation([1,2,3]), [0]*3,[1]*3); ieq The collection of inequalities integer: (3, 7, 11) x + 2 >= 0 sage: ieq.prepare_next_to_inner_loop([2,1,3]) sage: ieq.prepare_inner_loop([2,1,3]) sage: print_cache(ieq) Cached inner loop: 3 * x_0 + 42 >= 0 Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0 """ for ineq in self.ineqs_int: (ineq).prepare_next_to_inner_loop(p) for ineq in self.ineqs_generic: (ineq).prepare_next_to_inner_loop(p) def prepare_inner_loop(self, p): r""" Peel off the inner loop. In the inner loop of :func:rectangular_box_points, we have to repeatedly evaluate A x+b\geq 0. To speed up computation, we pre-evaluate .. math:: c = A x - A_0 x_0 +b = b + sum_{i=1} A_i x_i and only test A_0 x_0 +c \geq 0 in the inner loop. You must call :meth:prepare_next_to_inner_loop before calling this method. INPUT: - p -- the coordinates of the point to loop over. Only the p[1:] entries are used. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection, print_cache sage: P = Polyhedron(ieqs=[(2,3,7,11)]) sage: ieq = InequalityCollection(P, Permutation([1,2,3]), [0]*3,[1]*3); ieq The collection of inequalities integer: (3, 7, 11) x + 2 >= 0 sage: ieq.prepare_next_to_inner_loop([2,1,3]) sage: ieq.prepare_inner_loop([2,1,3]) sage: print_cache(ieq) Cached inner loop: 3 * x_0 + 42 >= 0 Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 35 >= 0 """ for ineq in self.ineqs_int: (ineq).prepare_inner_loop(p) for ineq in self.ineqs_generic: (ineq).prepare_inner_loop(p) def swap_ineq_to_front(self, i): r""" Swap the i-th entry of the list to the front of the list of inequalities. INPUT: - i -- Integer. The :class:Inequality_int to swap to the beginnig of the list of integral inequalities. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: P_QQ = Polyhedron(identity_matrix(3).columns() + [(-2, -1,-1)], field=QQ) sage: iec = InequalityCollection(P_QQ, Permutation([1,2,3]), [0]*3,[1]*3) sage: iec The collection of inequalities integer: (-1, -1, -1) x + 1 >= 0 integer: (-1, -1, 4) x + 1 >= 0 integer: (-1, 4, -1) x + 1 >= 0 integer: (3, -2, -2) x + 2 >= 0 sage: iec.swap_ineq_to_front(3) sage: iec The collection of inequalities integer: (3, -2, -2) x + 2 >= 0 integer: (-1, -1, -1) x + 1 >= 0 integer: (-1, -1, 4) x + 1 >= 0 integer: (-1, 4, -1) x + 1 >= 0 """ i_th_entry = self.ineqs_int.pop(i) self.ineqs_int.insert(0, i_th_entry) def are_satisfied(self, inner_loop_variable): r""" Return whether all inequalities are satisfied. You must call :meth:prepare_inner_loop before calling this method. INPUT: - inner_loop_variable -- Integer. the 0-th coordinate of the lattice point. OUTPUT: Boolean. Whether the lattice point is in the polyhedron. EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection sage: line = Polyhedron(eqns=[(2,3,7)]) sage: ieq = InequalityCollection(line, Permutation([1,2]), [0]*2, [1]*2 ) sage: ieq.prepare_next_to_inner_loop([3,4]) sage: ieq.prepare_inner_loop([3,4]) sage: ieq.are_satisfied(3) False """ cdef int i for i in range(0,len(self.ineqs_int)): ineq = self.ineqs_int[i] if (ineq).is_not_satisfied(inner_loop_variable): if i>0: self.swap_ineq_to_front(i) return False for i in range(0,len(self.ineqs_generic)): ineq = self.ineqs_generic[i] if (ineq).is_not_satisfied(inner_loop_variable): return False return True cpdef print_cache(InequalityCollection inequality_collection): r""" Print the cached values in :class:Inequality_int (for debugging/doctesting only). EXAMPLES:: sage: from sage.geometry.integral_points import InequalityCollection, print_cache sage: P = Polyhedron(ieqs=[(2,3,7)]) sage: ieq = InequalityCollection(P, Permutation([1,2]), [0]*2,[1]*2); ieq The collection of inequalities integer: (3, 7) x + 2 >= 0 sage: ieq.prepare_next_to_inner_loop([3,5]) sage: ieq.prepare_inner_loop([3,5]) sage: print_cache(ieq) Cached inner loop: 3 * x_0 + 37 >= 0 Cached next-to-inner loop: 3 * x_0 + 7 * x_1 + 2 >= 0 """ cdef Inequality_int ieq = (inequality_collection.ineqs_int[0]) print 'Cached inner loop: ' + \ str(ieq.coeff) + ' * x_0 + ' + str(ieq.cache) + ' >= 0' print 'Cached next-to-inner loop: ' + \ str(ieq.coeff) + ' * x_0 + ' + \ str(ieq.coeff_next) + ' * x_1 + ' + str(ieq.cache_next) + ' >= 0'
• ## sage/geometry/polyhedra.py

diff --git a/sage/geometry/polyhedra.py b/sage/geometry/polyhedra.py
 a from sage.categories.objects import Objects from subprocess import Popen, PIPE from sage.misc.all import tmp_filename from sage.misc.all import tmp_filename, cached_method, prod from sage.misc.functional import norm from sage.misc.package import is_package_installed """ Returns a string representation of the vertex. OUTPUT: String. TESTS:: sage: p = Polyhedron(ieqs = [[0,0,1],[0,1,0],[1,-1,0]]) """ if not self.is_compact(): raise ValueError, 'Can only enumerate points in a compact polyhedron.' lp = self.lattice_polytope(True) # remove cached values to get accurate timings try: del lp._points del lp._npoints except AttributeError: pass if self.is_lattice_polytope(): return lp.points().columns() points = filter(lambda p: self.contains(p), lp.points().columns()) return points def _integral_points_lattice_simplex(self, vertices=None): @cached_method def bounding_box(self, integral=False): r""" Return the integral points in a lattice simplex. Return the coordinates of a rectangular box containing the non-empty polytope. INPUT: - vertices -- None (default) or a list of integers. The indices of vertices that span the simplex under consideration. By default, it is assumed that the polytope itself is a simplex and all vertices are used. OUTPUT: A tuple containing the integral point coordinates as \ZZ-vectors. EXAMPLES:: - integral -- Boolean (default: False). Whether to only allow integral coordinates in the bounding box. OUTPUT: A pair of tuples (box_min, box_max) where box_min are the coordinates of a point bounding the coordinates of the polytope from below and box_max bounds the coordinates from above. EXAMPLES:: sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box() ((1/3, 1/3), (2/3, 2/3)) sage: Polyhedron([ (1/3,2/3), (2/3, 1/3) ]).bounding_box(integral=True) ((0, 0), (1, 1)) sage: polytopes.buckyball().bounding_box() ((-1059/1309, -1059/1309, -1059/1309), (1059/1309, 1059/1309, 1059/1309)) """ box_min = [] box_max = [] if self.n_vertices==0: raise ValueError('Empty polytope is not allowed') for i in range(0,self.ambient_dim()): coords = [ v[i] for v in self.Vrep_generator() ] max_coord = max(coords) min_coord = min(coords) if integral: box_max.append(ceil(max_coord)) box_min.append(floor(min_coord)) else: box_max.append(max_coord) box_min.append(min_coord) return (tuple(box_min), tuple(box_max)) def integral_points(self, threshold=100000): r""" Return the integral points in the polyhedron. Uses either the naive algorithm (iterate over a rectangular bounding box) or triangulation + Smith form. INPUT: - threshold -- integer (default: 100000). Use the naive algorith as long as the bounding box is smaller than this. OUTPUT: The list of integral points in the polyhedron. If the polyhedron is not compact, a ValueError is raised. EXAMPLES:: sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]).integral_points() ((-1, -1), (0, 0), (0, 1), (1, 0), (1, 1)) sage: simplex = Polyhedron([(1,2,3), (2,3,7), (-2,-3,-11)]) sage: simplex._integral_points_lattice_simplex() sage: simplex.integral_points() ((-2, -3, -11), (0, 0, -2), (1, 2, 3), (2, 3, 7)) The simplex need not be full-dimensional:: The polyhedron need not be full-dimensional:: sage: simplex = Polyhedron([(1,2,3,5), (2,3,7,5), (-2,-3,-11,5)]) sage: simplex._integral_points_lattice_simplex() sage: simplex.integral_points() ((-2, -3, -11, 5), (0, 0, -2, 5), (1, 2, 3, 5), (2, 3, 7, 5)) sage: point = Polyhedron([(2,3,7)]) sage: point._integral_points_lattice_simplex() sage: point.integral_points() ((2, 3, 7),) TESTS:: sage: empty = Polyhedron() sage: empty.integral_points() () Here is a simplex where the naive algorithm of running over all points in a rectangular bounding box no longer works fast enough:: sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)] sage: simplex = Polyhedron(v); simplex A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 5 vertices. sage: len(simplex._integral_points_lattice_simplex()) sage: len(simplex.integral_points()) 49 sage: v = [(4,-1,-1,-1), (-1,4,-1,-1), (-1,-1,4,-1), (-1,-1,-1,4), (-1,-1,-1,-1)] sage: P4mirror = Polyhedron(v) sage: len( P4mirror._integral_points_lattice_simplex() ) 126 """ if not vertices: assert self.is_simplex() and self.is_lattice_polytope() vertices = [ vector(ZZ, v) for v in self.Vrepresentation() ] else: vertices = [ vector(ZZ, self.Vrepresentation(i)) for i in vertices ] origin = vertices.pop() origin.set_immutable() if len(vertices)==0: return tuple([origin]) rays = [ v-origin for v in vertices ] # Find equation Ax<=b that cuts out simplex from parallelotope R = matrix(rays) if R.is_square(): b = 1 A = R.solve_right(vector([1]*len(rays))) else: b = 1 RRt = R * R.transpose() A = RRt.solve_right(vector([1]*len(rays))) * R from sage.geometry.cone import parallelotope_points gens = list(parallelotope_points(rays)) + list(rays) gens = [ origin+x for x in list(parallelotope_points(rays)) + list(rays) if A*x<=b ] for x in gens: x.set_immutable() return tuple(gens) def integral_points(self): r""" Return the integral points in the polyhedron. OUTPUT: The list of integral points in the polyhedron. If the polyhedron is not compact, a ValueError is raised. EXAMPLES:: sage: Polyhedron(vertices=[(-1,-1),(1,0),(1,1),(0,1)]).integral_points() ((0, 1), (1, 0), (0, 0), (1, 1), (-1, -1)) Finally, 3-d reflexive polytope number 4078:: Finally, the 3-d reflexive polytope number 4078:: sage: v = [(1,0,0), (0,1,0), (0,0,1), (0,0,-1), (0,-2,1), ...        (-1,2,-1), (-1,2,-2), (-1,1,-2), (-1,-1,2), (-1,-3,2)] sage: pts1 = set(Polyhedron(v).integral_points())    # Sage's own code sage: P = Polyhedron(v) sage: pts1 = P.integral_points()                     # Sage's own code sage: all(P.contains(p) for p in pts1) True sage: pts2 = LatticePolytope(v).points().columns()   # PALP sage: for p in pts2: ...       p.set_immutable() sage: pts2 = set(pts2) sage: pts1 == pts2 sage: for p in pts1: p.set_immutable() sage: for p in pts2: p.set_immutable() sage: set(pts1) == set(pts2) True TODO: profile and make faster:: sage: timeit('Polyhedron(v).integral_points()')   # random output sage: timeit('LatticePolytope(v).points()')       # random output """ if not self.is_compact(): raise ValueError('Can only enumerate points in a compact polyhedron.') if not self.is_lattice_polytope(): raise NotImplementedError('Only implemented for polyhedra with integral vertices') if self.n_vertices() == 0: return tuple() # for small bounding boxes, it is faster to naively iterate over the points of the box box_min, box_max = self.bounding_box(integral=True) box_points = prod(max_coord-min_coord+1 for min_coord, max_coord in zip(box_min, box_max)) if  not self.is_lattice_polytope() or \ (self.is_simplex() and box_points<1000) or \ box_points