From 4e9a8bd3727e33cf113b62b0ca61628947ee3f99 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Sun, 30 Jul 2023 12:36:05 -0700 Subject: [PATCH 01/33] LatticePolytope: Use experimental method PermutationGroupElement._transpose_left --- src/sage/geometry/lattice_polytope.py | 35 +++++++++---------- .../groups/perm_gps/permgroup_element.pyx | 10 ++++++ 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 26c7d9a3819..a3d17d6cd51 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3215,7 +3215,7 @@ def index_of_max(iterable): for j in range(n_v): m = index_of_max(PM[0, i] for i in range(j, n_v)) if m > 0: - permutations[0][1] = S_v((j + 1, m + j + 1), check=False) * permutations[0][1] + permutations[0][1] = permutations[0][1]._transpose_left(j + 1, m + j + 1) first_row = list(PM[0]) # Arrange other rows one by one and compare with first row @@ -3224,7 +3224,7 @@ def index_of_max(iterable): permutations[n_s] = [S_f.one(), S_v.one()] m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(n_v)) if m > 0: - permutations[n_s][1] = S_v((1, m + 1), check=False) * permutations[n_s][1] + permutations[n_s][1] = permutations[n_s][1]._transpose_left(1, m + 1) d = (PM[k, permutations[n_s][1](1) - 1] - permutations[0][1](first_row)[0]) if d < 0: @@ -3235,8 +3235,7 @@ def index_of_max(iterable): for i in range(1, n_v): m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(i,n_v)) if m > 0: - permutations[n_s][1] = S_v((i + 1, m + i + 1), check=False) \ - * permutations[n_s][1] + permutations[n_s][1] = permutations[n_s][1]._transpose_left(i + 1, m + i + 1) if d == 0: d = (PM[k, permutations[n_s][1](i+1) - 1] -permutations[0][1](first_row)[i]) @@ -3246,7 +3245,7 @@ def index_of_max(iterable): # This row is smaller than 1st row, so nothing to do del permutations[n_s] continue - permutations[n_s][0] = S_f((1, k + 1), check=False) * permutations[n_s][0] + permutations[n_s][0] = permutations[n_s][0]._transpose_left(1, k + 1) if d == 0: # This row is the same, so we have a symmetry! n_s += 1 @@ -3286,7 +3285,7 @@ def index_of_max(iterable): # number of local permutations associated with current global n_p = 0 ccf = cf - permutations_bar = {0: copy(permutations[k])} + permutations_bar = {0: list(permutations[k])} # We look for the line with the maximal entry in the first # subsymmetry block, i.e. we are allowed to swap elements # between 0 and S(0) @@ -3295,14 +3294,14 @@ def index_of_max(iterable): v0 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] vj = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](j+1) - 1] if v0 < vj: - permutations_bar[n_p][1] = S_v((1, j + 1), check=False) * permutations_bar[n_p][1] + permutations_bar[n_p][1] = permutations_bar[n_p][1]._transpose_left(1, j + 1) if ccf == 0: l_r[0] = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] if s != l: - permutations_bar[n_p][0] = S_f((l + 1, s + 1), check=False) * permutations_bar[n_p][0] + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) n_p += 1 ccf = 1 - permutations_bar[n_p] = copy(permutations[k]) + permutations_bar[n_p] = list(permutations[k]) else: d1 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] d = d1 - l_r[0] @@ -3312,20 +3311,20 @@ def index_of_max(iterable): elif d==0: # Maximal values agree, so possible symmetry if s != l: - permutations_bar[n_p][0] = S_f((l + 1, s + 1), check=False) * permutations_bar[n_p][0] + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) n_p += 1 - permutations_bar[n_p] = copy(permutations[k]) + permutations_bar[n_p] = list(permutations[k]) else: # We found a greater maximal value for first entry. # It becomes our new reference: l_r[0] = d1 if s != l: - permutations_bar[n_p][0] = S_f((l + 1, s + 1), check=False) * permutations_bar[n_p][0] + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) # Forget previous work done cf = 0 - permutations_bar = {0:copy(permutations_bar[n_p])} + permutations_bar = {0: list(permutations_bar[n_p])} n_p = 1 - permutations_bar[n_p] = copy(permutations[k]) + permutations_bar[n_p] = list(permutations[k]) n_s = k + 1 # Check if the permutations found just now work # with other elements @@ -3345,7 +3344,7 @@ def index_of_max(iterable): vc = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] vj = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(j+1) - 1] if (vc < vj): - permutations_bar[s][1] = S_v((c + 1, j + 1), check=False) * permutations_bar[s][1] + permutations_bar[s][1] = permutations_bar[s][1]._transpose_left(c + 1, j + 1) if ccf == 0: # Set reference and carry on to next permutation l_r[c] = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] @@ -3356,7 +3355,7 @@ def index_of_max(iterable): if d < 0: n_p -= 1 if s < n_p: - permutations_bar[s] = copy(permutations_bar[n_p]) + permutations_bar[s] = list(permutations_bar[n_p]) elif d > 0: # The current case leads to a smaller matrix, # hence this case becomes our new reference @@ -3366,10 +3365,10 @@ def index_of_max(iterable): n_s = k + 1 # Update permutations if (n_s - 1) > k: - permutations[k] = copy(permutations[n_s - 1]) + permutations[k] = list(permutations[n_s - 1]) n_s -= 1 for s in range(n_p): - permutations[n_s] = copy(permutations_bar[s]) + permutations[n_s] = list(permutations_bar[s]) n_s += 1 cf = n_s permutations = {k: permutations[k] for k in permutations if k < n_s} diff --git a/src/sage/groups/perm_gps/permgroup_element.pyx b/src/sage/groups/perm_gps/permgroup_element.pyx index 2522b5c346d..2b36616b27e 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pyx +++ b/src/sage/groups/perm_gps/permgroup_element.pyx @@ -1298,6 +1298,16 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): return coercion_model.bin_op(left, right, operator.mul) + def _transpose_left(self, j, k): + if j == k: + return self + cdef PermutationGroupElement prod = self._new_c() + cdef int i + for i from 0 <= i < self.n: + prod.perm[i] = self.perm[i] + prod.perm[j - 1], prod.perm[k - 1] = self.perm[k - 1], self.perm[j - 1] + return prod + cpdef _mul_(left, _right): r""" EXAMPLES:: From 684f35968d4c1691c7aaff87e1ac5403ac9cafb2 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 2 Aug 2023 20:24:47 -0700 Subject: [PATCH 02/33] sage.geometry.palp_normal_form: Cythonized _palp_PM_max method --- src/sage/geometry/lattice_polytope.py | 200 +------------------ src/sage/geometry/palp_normal_form.pyx | 262 +++++++++++++++++++++++++ 2 files changed, 264 insertions(+), 198 deletions(-) create mode 100644 src/sage/geometry/palp_normal_form.pyx diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index a3d17d6cd51..1ae8504a306 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3199,204 +3199,8 @@ def _palp_PM_max(self, check=False): sage: all(results) # long time True """ - PM = self.vertex_facet_pairing_matrix() - n_v = PM.ncols() - n_f = PM.nrows() - S_v = SymmetricGroup(n_v) - S_f = SymmetricGroup(n_f) - - # and find all the ways of making the first row of PM_max - def index_of_max(iterable): - # returns the index of max of any iterable - return max(enumerate(iterable), key=lambda x: x[1])[0] - - n_s = 1 - permutations = {0: [S_f.one(), S_v.one()]} - for j in range(n_v): - m = index_of_max(PM[0, i] for i in range(j, n_v)) - if m > 0: - permutations[0][1] = permutations[0][1]._transpose_left(j + 1, m + j + 1) - first_row = list(PM[0]) - - # Arrange other rows one by one and compare with first row - for k in range(1, n_f): - # Error for k == 1 already! - permutations[n_s] = [S_f.one(), S_v.one()] - m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(n_v)) - if m > 0: - permutations[n_s][1] = permutations[n_s][1]._transpose_left(1, m + 1) - d = (PM[k, permutations[n_s][1](1) - 1] - - permutations[0][1](first_row)[0]) - if d < 0: - # The largest elt of this row is smaller than largest elt - # in 1st row, so nothing to do - continue - # otherwise: - for i in range(1, n_v): - m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(i,n_v)) - if m > 0: - permutations[n_s][1] = permutations[n_s][1]._transpose_left(i + 1, m + i + 1) - if d == 0: - d = (PM[k, permutations[n_s][1](i+1) - 1] - -permutations[0][1](first_row)[i]) - if d < 0: - break - if d < 0: - # This row is smaller than 1st row, so nothing to do - del permutations[n_s] - continue - permutations[n_s][0] = permutations[n_s][0]._transpose_left(1, k + 1) - if d == 0: - # This row is the same, so we have a symmetry! - n_s += 1 - else: - # This row is larger, so it becomes the first row and - # the symmetries reset. - first_row = list(PM[k]) - permutations = {0: permutations[n_s]} - n_s = 1 - permutations = {k: permutations[k] for k in permutations if k < n_s} - - b = tuple(PM[permutations[0][0](1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) - # Work out the restrictions the current permutations - # place on other permutations as a automorphisms - # of the first row - # The array is such that: - # S = [i, 1, ..., 1 (ith), j, i+1, ..., i+1 (jth), k ... ] - # describes the "symmetry blocks" - S = list(range(1, n_v + 1)) - for i in range(1, n_v): - if b[i-1] == b[i]: - S[i] = S[i-1] - S[S[i]-1] += 1 - else: - S[i] = i + 1 - - # We determine the other rows of PM_max in turn by use of perms and - # aut on previous rows. - for l in range(1, n_f - 1): - n_s = len(permutations) - n_s_bar = n_s - cf = 0 - l_r = [0]*n_v - # Search for possible local permutations based off previous - # global permutations. - for k in range(n_s_bar - 1, -1, -1): - # number of local permutations associated with current global - n_p = 0 - ccf = cf - permutations_bar = {0: list(permutations[k])} - # We look for the line with the maximal entry in the first - # subsymmetry block, i.e. we are allowed to swap elements - # between 0 and S(0) - for s in range(l, n_f): - for j in range(1, S[0]): - v0 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] - vj = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](j+1) - 1] - if v0 < vj: - permutations_bar[n_p][1] = permutations_bar[n_p][1]._transpose_left(1, j + 1) - if ccf == 0: - l_r[0] = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] - if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) - n_p += 1 - ccf = 1 - permutations_bar[n_p] = list(permutations[k]) - else: - d1 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] - d = d1 - l_r[0] - if d < 0: - # We move to the next line - continue - elif d==0: - # Maximal values agree, so possible symmetry - if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) - n_p += 1 - permutations_bar[n_p] = list(permutations[k]) - else: - # We found a greater maximal value for first entry. - # It becomes our new reference: - l_r[0] = d1 - if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) - # Forget previous work done - cf = 0 - permutations_bar = {0: list(permutations_bar[n_p])} - n_p = 1 - permutations_bar[n_p] = list(permutations[k]) - n_s = k + 1 - # Check if the permutations found just now work - # with other elements - for c in range(1, n_v): - h = S[c] - ccf = cf - # Now let us find out where the end of the - # next symmetry block is: - if h < c + 1: - h = S[h - 1] - s = n_p - # Check through this block for each possible permutation - while s > 0: - s -= 1 - # Find the largest value in this symmetry block - for j in range(c + 1, h): - vc = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] - vj = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(j+1) - 1] - if (vc < vj): - permutations_bar[s][1] = permutations_bar[s][1]._transpose_left(c + 1, j + 1) - if ccf == 0: - # Set reference and carry on to next permutation - l_r[c] = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] - ccf = 1 - else: - d1 = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] - d = d1 - l_r[c] - if d < 0: - n_p -= 1 - if s < n_p: - permutations_bar[s] = list(permutations_bar[n_p]) - elif d > 0: - # The current case leads to a smaller matrix, - # hence this case becomes our new reference - l_r[c] = d1 - cf = 0 - n_p = s + 1 - n_s = k + 1 - # Update permutations - if (n_s - 1) > k: - permutations[k] = list(permutations[n_s - 1]) - n_s -= 1 - for s in range(n_p): - permutations[n_s] = list(permutations_bar[s]) - n_s += 1 - cf = n_s - permutations = {k: permutations[k] for k in permutations if k < n_s} - # If the automorphisms are not already completely restricted, - # update them - if S != list(range(1, n_v + 1)): - # Take the old automorphisms and update by - # the restrictions the last worked out - # row imposes. - c = 0 - M = tuple(PM[permutations[0][0](l+1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) - while c < n_v: - s = S[c] + 1 - S[c] = c + 1 - c += 1 - while c < (s - 1): - if M[c] == M[c - 1]: - S[c] = S[c - 1] - S[S[c] - 1] += 1 - else: - S[c] = c + 1 - c += 1 - # Now we have the perms, we construct PM_max using one of them - PM_max = PM.with_permuted_rows_and_columns(*permutations[0]) - if check: - return (PM_max, permutations) - else: - return PM_max + from .palp_normal_form import _palp_PM_max + return _palp_PM_max(self, check) def npoints(self): r""" diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx new file mode 100644 index 00000000000..77fffa8b761 --- /dev/null +++ b/src/sage/geometry/palp_normal_form.pyx @@ -0,0 +1,262 @@ +from sage.groups.perm_gps.permgroup_named import SymmetricGroup + + +def _palp_PM_max(self, check=False): + r""" + Compute the permutation normal form of the vertex facet pairing + matrix . + + The permutation normal form of a matrix is defined as the lexicographic + maximum under all permutations of its rows and columns. For more + more detail, see also + :meth:`~sage.matrix.matrix2.Matrix.permutation_normal_form`. + + Instead of using the generic method for computing the permutation + normal form, this method uses the PALP algorithm to compute + the permutation normal form and its automorphisms concurrently. + + INPUT: + + - ``check`` -- Boolean (default: ``False``), whether to return + the permutations leaving the maximal vertex-facet pairing + matrix invariant. + + OUTPUT: + + A matrix or a tuple of a matrix and a dict whose values are the + permutation group elements corresponding to the permutations + that permute :meth:`vertices` such that the vertex-facet pairing + matrix is maximal. + + EXAMPLES:: + + sage: o = lattice_polytope.cross_polytope(2) + sage: PM = o.vertex_facet_pairing_matrix() + sage: PM_max = PM.permutation_normal_form() # optional - sage.graphs + sage: PM_max == o._palp_PM_max() # optional - sage.graphs + True + sage: P2 = ReflexivePolytope(2, 0) + sage: PM_max, permutations = P2._palp_PM_max(check=True) + sage: PM_max + [3 0 0] + [0 3 0] + [0 0 3] + sage: list(permutations.values()) + [[(1,2,3), (1,2,3)], + [(1,3,2), (1,3,2)], + [(1,3), (1,3)], + [(1,2), (1,2)], + [(), ()], + [(2,3), (2,3)]] + sage: PM_max.automorphisms_of_rows_and_columns() + [((), ()), + ((1,2,3), (1,2,3)), + ((1,3,2), (1,3,2)), + ((2,3), (2,3)), + ((1,2), (1,2)), + ((1,3), (1,3))] + sage: PMs = ( i._palp_PM_max(check=True) + ....: for i in ReflexivePolytopes(2) ) + sage: results = ( len(i) == len(j.automorphisms_of_rows_and_columns()) + ....: for j, i in PMs ) + sage: all(results) # long time + True + """ + PM = self.vertex_facet_pairing_matrix() + n_v = PM.ncols() + n_f = PM.nrows() + S_v = SymmetricGroup(n_v) + S_f = SymmetricGroup(n_f) + + # and find all the ways of making the first row of PM_max + def index_of_max(iterable): + # returns the index of max of any iterable + return max(enumerate(iterable), key=lambda x: x[1])[0] + + n_s = 1 + permutations = {0: [S_f.one(), S_v.one()]} + for j in range(n_v): + m = index_of_max(PM[0, i] for i in range(j, n_v)) + if m > 0: + permutations[0][1] = permutations[0][1]._transpose_left(j + 1, m + j + 1) + first_row = list(PM[0]) + + # Arrange other rows one by one and compare with first row + for k in range(1, n_f): + # Error for k == 1 already! + permutations[n_s] = [S_f.one(), S_v.one()] + m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(n_v)) + if m > 0: + permutations[n_s][1] = permutations[n_s][1]._transpose_left(1, m + 1) + d = (PM[k, permutations[n_s][1](1) - 1] + - permutations[0][1](first_row)[0]) + if d < 0: + # The largest elt of this row is smaller than largest elt + # in 1st row, so nothing to do + continue + # otherwise: + for i in range(1, n_v): + m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(i,n_v)) + if m > 0: + permutations[n_s][1] = permutations[n_s][1]._transpose_left(i + 1, m + i + 1) + if d == 0: + d = (PM[k, permutations[n_s][1](i+1) - 1] + -permutations[0][1](first_row)[i]) + if d < 0: + break + if d < 0: + # This row is smaller than 1st row, so nothing to do + del permutations[n_s] + continue + permutations[n_s][0] = permutations[n_s][0]._transpose_left(1, k + 1) + if d == 0: + # This row is the same, so we have a symmetry! + n_s += 1 + else: + # This row is larger, so it becomes the first row and + # the symmetries reset. + first_row = list(PM[k]) + permutations = {0: permutations[n_s]} + n_s = 1 + permutations = {k: permutations[k] for k in permutations if k < n_s} + + b = tuple(PM[permutations[0][0](1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) + # Work out the restrictions the current permutations + # place on other permutations as a automorphisms + # of the first row + # The array is such that: + # S = [i, 1, ..., 1 (ith), j, i+1, ..., i+1 (jth), k ... ] + # describes the "symmetry blocks" + S = list(range(1, n_v + 1)) + for i in range(1, n_v): + if b[i-1] == b[i]: + S[i] = S[i-1] + S[S[i]-1] += 1 + else: + S[i] = i + 1 + + # We determine the other rows of PM_max in turn by use of perms and + # aut on previous rows. + for l in range(1, n_f - 1): + n_s = len(permutations) + n_s_bar = n_s + cf = 0 + l_r = [0]*n_v + # Search for possible local permutations based off previous + # global permutations. + for k in range(n_s_bar - 1, -1, -1): + # number of local permutations associated with current global + n_p = 0 + ccf = cf + permutations_bar = {0: list(permutations[k])} + # We look for the line with the maximal entry in the first + # subsymmetry block, i.e. we are allowed to swap elements + # between 0 and S(0) + for s in range(l, n_f): + for j in range(1, S[0]): + v0 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] + vj = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](j+1) - 1] + if v0 < vj: + permutations_bar[n_p][1] = permutations_bar[n_p][1]._transpose_left(1, j + 1) + if ccf == 0: + l_r[0] = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] + if s != l: + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + n_p += 1 + ccf = 1 + permutations_bar[n_p] = list(permutations[k]) + else: + d1 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] + d = d1 - l_r[0] + if d < 0: + # We move to the next line + continue + elif d==0: + # Maximal values agree, so possible symmetry + if s != l: + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + n_p += 1 + permutations_bar[n_p] = list(permutations[k]) + else: + # We found a greater maximal value for first entry. + # It becomes our new reference: + l_r[0] = d1 + if s != l: + permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + # Forget previous work done + cf = 0 + permutations_bar = {0: list(permutations_bar[n_p])} + n_p = 1 + permutations_bar[n_p] = list(permutations[k]) + n_s = k + 1 + # Check if the permutations found just now work + # with other elements + for c in range(1, n_v): + h = S[c] + ccf = cf + # Now let us find out where the end of the + # next symmetry block is: + if h < c + 1: + h = S[h - 1] + s = n_p + # Check through this block for each possible permutation + while s > 0: + s -= 1 + # Find the largest value in this symmetry block + for j in range(c + 1, h): + vc = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] + vj = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(j+1) - 1] + if (vc < vj): + permutations_bar[s][1] = permutations_bar[s][1]._transpose_left(c + 1, j + 1) + if ccf == 0: + # Set reference and carry on to next permutation + l_r[c] = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] + ccf = 1 + else: + d1 = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] + d = d1 - l_r[c] + if d < 0: + n_p -= 1 + if s < n_p: + permutations_bar[s] = list(permutations_bar[n_p]) + elif d > 0: + # The current case leads to a smaller matrix, + # hence this case becomes our new reference + l_r[c] = d1 + cf = 0 + n_p = s + 1 + n_s = k + 1 + # Update permutations + if (n_s - 1) > k: + permutations[k] = list(permutations[n_s - 1]) + n_s -= 1 + for s in range(n_p): + permutations[n_s] = list(permutations_bar[s]) + n_s += 1 + cf = n_s + permutations = {k: permutations[k] for k in permutations if k < n_s} + # If the automorphisms are not already completely restricted, + # update them + if S != list(range(1, n_v + 1)): + # Take the old automorphisms and update by + # the restrictions the last worked out + # row imposes. + c = 0 + M = tuple(PM[permutations[0][0](l+1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) + while c < n_v: + s = S[c] + 1 + S[c] = c + 1 + c += 1 + while c < (s - 1): + if M[c] == M[c - 1]: + S[c] = S[c - 1] + S[S[c] - 1] += 1 + else: + S[c] = c + 1 + c += 1 + # Now we have the perms, we construct PM_max using one of them + PM_max = PM.with_permuted_rows_and_columns(*permutations[0]) + if check: + return (PM_max, permutations) + else: + return PM_max From 26c8607cc4d47bc180537fc45401cd38d6e91e2c Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 2 Aug 2023 21:06:45 -0700 Subject: [PATCH 03/33] src/sage/geometry/palp_normal_form.pyx: More cythonization --- src/sage/geometry/palp_normal_form.pyx | 78 ++++++++++++++++---------- 1 file changed, 47 insertions(+), 31 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 77fffa8b761..ec1a1550e36 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -1,4 +1,6 @@ +from sage.groups.perm_gps.permgroup_element cimport PermutationGroupElement from sage.groups.perm_gps.permgroup_named import SymmetricGroup +from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense def _palp_PM_max(self, check=False): @@ -62,9 +64,9 @@ def _palp_PM_max(self, check=False): sage: all(results) # long time True """ - PM = self.vertex_facet_pairing_matrix() - n_v = PM.ncols() - n_f = PM.nrows() + cdef Matrix_integer_dense PM = self.vertex_facet_pairing_matrix() + cdef int n_v = PM.ncols() + cdef int n_f = PM.nrows() S_v = SymmetricGroup(n_v) S_f = SymmetricGroup(n_f) @@ -73,12 +75,14 @@ def _palp_PM_max(self, check=False): # returns the index of max of any iterable return max(enumerate(iterable), key=lambda x: x[1])[0] - n_s = 1 - permutations = {0: [S_f.one(), S_v.one()]} + cdef int n_s = 1 + cdef dict permutations = {0: [S_f.one(), S_v.one()]} + cdef int j, k, m, d + for j in range(n_v): m = index_of_max(PM[0, i] for i in range(j, n_v)) if m > 0: - permutations[0][1] = permutations[0][1]._transpose_left(j + 1, m + j + 1) + permutations[0][1] = ( permutations[0][1])._transpose_left(j + 1, m + j + 1) first_row = list(PM[0]) # Arrange other rows one by one and compare with first row @@ -87,28 +91,27 @@ def _palp_PM_max(self, check=False): permutations[n_s] = [S_f.one(), S_v.one()] m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(n_v)) if m > 0: - permutations[n_s][1] = permutations[n_s][1]._transpose_left(1, m + 1) - d = (PM[k, permutations[n_s][1](1) - 1] - - permutations[0][1](first_row)[0]) + permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(1, m + 1) + d = PM[k, ( permutations[n_s][1])(1) - 1] - ( permutations[0][1])(first_row)[0] if d < 0: # The largest elt of this row is smaller than largest elt # in 1st row, so nothing to do continue # otherwise: for i in range(1, n_v): - m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(i,n_v)) + m = index_of_max(PM[k, ( permutations[n_s][1](j+1)) - 1] for j in range(i,n_v)) if m > 0: - permutations[n_s][1] = permutations[n_s][1]._transpose_left(i + 1, m + i + 1) + permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m + i + 1) if d == 0: - d = (PM[k, permutations[n_s][1](i+1) - 1] - -permutations[0][1](first_row)[i]) + d = (PM[k, ( permutations[n_s][1])(i+1) - 1] + - ( permutations[0][1])(first_row)[i]) if d < 0: break if d < 0: # This row is smaller than 1st row, so nothing to do del permutations[n_s] continue - permutations[n_s][0] = permutations[n_s][0]._transpose_left(1, k + 1) + permutations[n_s][0] = ( permutations[n_s][0])._transpose_left(1, k + 1) if d == 0: # This row is the same, so we have a symmetry! n_s += 1 @@ -120,14 +123,15 @@ def _palp_PM_max(self, check=False): n_s = 1 permutations = {k: permutations[k] for k in permutations if k < n_s} - b = tuple(PM[permutations[0][0](1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) + cdef tuple b = tuple(PM[( permutations[0][0])(1) - 1, + ( permutations[0][1])(j+1) - 1] for j in range(n_v)) # Work out the restrictions the current permutations # place on other permutations as a automorphisms # of the first row # The array is such that: # S = [i, 1, ..., 1 (ith), j, i+1, ..., i+1 (jth), k ... ] # describes the "symmetry blocks" - S = list(range(1, n_v + 1)) + cdef list S = list(range(1, n_v + 1)) for i in range(1, n_v): if b[i-1] == b[i]: S[i] = S[i-1] @@ -135,6 +139,9 @@ def _palp_PM_max(self, check=False): else: S[i] = i + 1 + cdef int l, np, cf, ccf, n_s_bar, d1, v0, vc, vj + cdef list l_r + # We determine the other rows of PM_max in turn by use of perms and # aut on previous rows. for l in range(1, n_f - 1): @@ -154,19 +161,23 @@ def _palp_PM_max(self, check=False): # between 0 and S(0) for s in range(l, n_f): for j in range(1, S[0]): - v0 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] - vj = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](j+1) - 1] + v0 = PM[( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1] + vj = PM[( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(j+1) - 1] if v0 < vj: - permutations_bar[n_p][1] = permutations_bar[n_p][1]._transpose_left(1, j + 1) + permutations_bar[n_p][1] = ( permutations_bar[n_p][1])._transpose_left(1, j + 1) if ccf == 0: - l_r[0] = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] + l_r[0] = PM[( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1] if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + permutations_bar[n_p][0] = ( permutations_bar[n_p][0])._transpose_left(l + 1, s + 1) n_p += 1 ccf = 1 permutations_bar[n_p] = list(permutations[k]) else: - d1 = PM[permutations_bar[n_p][0](s+1) - 1, permutations_bar[n_p][1](1) - 1] + d1 = PM[( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1] d = d1 - l_r[0] if d < 0: # We move to the next line @@ -174,7 +185,7 @@ def _palp_PM_max(self, check=False): elif d==0: # Maximal values agree, so possible symmetry if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + permutations_bar[n_p][0] = ( permutations_bar[n_p][0])._transpose_left(l + 1, s + 1) n_p += 1 permutations_bar[n_p] = list(permutations[k]) else: @@ -182,7 +193,7 @@ def _palp_PM_max(self, check=False): # It becomes our new reference: l_r[0] = d1 if s != l: - permutations_bar[n_p][0] = permutations_bar[n_p][0]._transpose_left(l + 1, s + 1) + permutations_bar[n_p][0] = ( permutations_bar[n_p][0])._transpose_left(l + 1, s + 1) # Forget previous work done cf = 0 permutations_bar = {0: list(permutations_bar[n_p])} @@ -204,16 +215,20 @@ def _palp_PM_max(self, check=False): s -= 1 # Find the largest value in this symmetry block for j in range(c + 1, h): - vc = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] - vj = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(j+1) - 1] - if (vc < vj): - permutations_bar[s][1] = permutations_bar[s][1]._transpose_left(c + 1, j + 1) + vc = PM[( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1] + vj = PM[( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(j+1) - 1] + if vc < vj: + permutations_bar[s][1] = ( permutations_bar[s][1])._transpose_left(c + 1, j + 1) if ccf == 0: # Set reference and carry on to next permutation - l_r[c] = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] + l_r[c] = PM[( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1] ccf = 1 else: - d1 = PM[(permutations_bar[s][0])(l+1) - 1, (permutations_bar[s][1])(c+1) - 1] + d1 = PM[( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1] d = d1 - l_r[c] if d < 0: n_p -= 1 @@ -242,7 +257,8 @@ def _palp_PM_max(self, check=False): # the restrictions the last worked out # row imposes. c = 0 - M = tuple(PM[permutations[0][0](l+1) - 1, permutations[0][1](j+1) - 1] for j in range(n_v)) + M = tuple(PM[( permutations[0][0])(l+1) - 1, + ( permutations[0][1])(j+1) - 1] for j in range(n_v)) while c < n_v: s = S[c] + 1 S[c] = c + 1 From e5d593549b16152bb6e439f3e71d426ba4f69b1c Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 12:56:17 -0700 Subject: [PATCH 04/33] src/sage/geometry/palp_normal_form.pyx: Cython loop to replace index_of_max --- src/sage/geometry/palp_normal_form.pyx | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index ec1a1550e36..27b847535d5 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -78,6 +78,7 @@ def _palp_PM_max(self, check=False): cdef int n_s = 1 cdef dict permutations = {0: [S_f.one(), S_v.one()]} cdef int j, k, m, d + cdef int element, max_element for j in range(n_v): m = index_of_max(PM[0, i] for i in range(j, n_v)) @@ -89,7 +90,7 @@ def _palp_PM_max(self, check=False): for k in range(1, n_f): # Error for k == 1 already! permutations[n_s] = [S_f.one(), S_v.one()] - m = index_of_max(PM[k, permutations[n_s][1](j+1) - 1] for j in range(n_v)) + m = index_of_max(PM[k, ( permutations[n_s][1])(j+1) - 1] for j in range(n_v)) if m > 0: permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(1, m + 1) d = PM[k, ( permutations[n_s][1])(1) - 1] - ( permutations[0][1])(first_row)[0] @@ -99,9 +100,15 @@ def _palp_PM_max(self, check=False): continue # otherwise: for i in range(1, n_v): - m = index_of_max(PM[k, ( permutations[n_s][1](j+1)) - 1] for j in range(i,n_v)) - if m > 0: - permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m + i + 1) + max_element = PM[k, ( permutations[n_s][1](i + 1)) - 1] + m = i + 1 + for j in range(i + 1, n_v): + element = PM[k, ( permutations[n_s][1](j+1)) - 1] + if element > max_element: + max_element = element + m = j + if m > i + 1: + permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m) if d == 0: d = (PM[k, ( permutations[n_s][1])(i+1) - 1] - ( permutations[0][1])(first_row)[i]) From 1425cb13d28b94dbe5e7c912f76eddf556392bd7 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 14:12:10 -0700 Subject: [PATCH 05/33] src/sage/geometry/palp_normal_form.pyx: Cython loop to replace index_of_max (fixup) --- src/sage/geometry/palp_normal_form.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 27b847535d5..dbb37844f1e 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -100,15 +100,15 @@ def _palp_PM_max(self, check=False): continue # otherwise: for i in range(1, n_v): - max_element = PM[k, ( permutations[n_s][1](i + 1)) - 1] - m = i + 1 + max_element = PM[k, ( permutations[n_s][1])(i + 1) - 1] + m = i for j in range(i + 1, n_v): - element = PM[k, ( permutations[n_s][1](j+1)) - 1] + element = PM[k, ( permutations[n_s][1])(j + 1) - 1] if element > max_element: max_element = element m = j - if m > i + 1: - permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m) + if m > i: + permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m + 1) if d == 0: d = (PM[k, ( permutations[n_s][1])(i+1) - 1] - ( permutations[0][1])(first_row)[i]) From 5fba248fa287813384f4c7482595b43e8d1a41fb Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 16:01:31 -0700 Subject: [PATCH 06/33] src/sage/geometry/palp_normal_form.pyx: Do not keep first_row as a list --- src/sage/geometry/palp_normal_form.pyx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index dbb37844f1e..e97f12c6dc5 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -84,7 +84,8 @@ def _palp_PM_max(self, check=False): m = index_of_max(PM[0, i] for i in range(j, n_v)) if m > 0: permutations[0][1] = ( permutations[0][1])._transpose_left(j + 1, m + j + 1) - first_row = list(PM[0]) + + cdef first_row_index = 0 # Arrange other rows one by one and compare with first row for k in range(1, n_f): @@ -93,7 +94,8 @@ def _palp_PM_max(self, check=False): m = index_of_max(PM[k, ( permutations[n_s][1])(j+1) - 1] for j in range(n_v)) if m > 0: permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(1, m + 1) - d = PM[k, ( permutations[n_s][1])(1) - 1] - ( permutations[0][1])(first_row)[0] + d = (PM[k, ( permutations[n_s][1])(1) - 1] + - PM[first_row_index, ( permutations[0][1])(1) - 1]) if d < 0: # The largest elt of this row is smaller than largest elt # in 1st row, so nothing to do @@ -111,7 +113,7 @@ def _palp_PM_max(self, check=False): permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m + 1) if d == 0: d = (PM[k, ( permutations[n_s][1])(i+1) - 1] - - ( permutations[0][1])(first_row)[i]) + - PM[first_row_index, ( permutations[0][1])(i + 1) - 1]) if d < 0: break if d < 0: @@ -125,7 +127,7 @@ def _palp_PM_max(self, check=False): else: # This row is larger, so it becomes the first row and # the symmetries reset. - first_row = list(PM[k]) + first_row_index = k permutations = {0: permutations[n_s]} n_s = 1 permutations = {k: permutations[k] for k in permutations if k < n_s} From 524c37c767ff501c3bd3485c729be4180489cb02 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 17:38:02 -0700 Subject: [PATCH 07/33] src/sage/geometry/palp_normal_form.pyx: Replace another use of index_of_max --- src/sage/geometry/palp_normal_form.pyx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index e97f12c6dc5..e772cb9455c 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -91,7 +91,13 @@ def _palp_PM_max(self, check=False): for k in range(1, n_f): # Error for k == 1 already! permutations[n_s] = [S_f.one(), S_v.one()] - m = index_of_max(PM[k, ( permutations[n_s][1])(j+1) - 1] for j in range(n_v)) + max_element = PM[k, ( permutations[n_s][1])(1) - 1] + m = 0 + for j in range(1, n_v): + element = PM[k, ( permutations[n_s][1])(j + 1) - 1] + if element > max_element: + max_element = element + m = j if m > 0: permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(1, m + 1) d = (PM[k, ( permutations[n_s][1])(1) - 1] From 011f6989bb74a75968fede4880d46f9f78d25f14 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 17:38:42 -0700 Subject: [PATCH 08/33] src/sage/groups/perm_gps/permgroup_element.pxd: Add cpdef --- src/sage/groups/perm_gps/permgroup_element.pxd | 1 + src/sage/groups/perm_gps/permgroup_element.pyx | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/sage/groups/perm_gps/permgroup_element.pxd b/src/sage/groups/perm_gps/permgroup_element.pxd index 0a584745f96..b013451d5a0 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pxd +++ b/src/sage/groups/perm_gps/permgroup_element.pxd @@ -18,6 +18,7 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): cpdef _set_permutation_group_element(self, PermutationGroupElement p, bint convert) cpdef _mul_(self, other) + cpdef PermutationGroupElement _transpose_left(self, int j, int k) cpdef PermutationGroupElement _generate_new(self, list new_list) cpdef PermutationGroupElement _generate_new_GAP(self, old) cpdef _gap_list(self) diff --git a/src/sage/groups/perm_gps/permgroup_element.pyx b/src/sage/groups/perm_gps/permgroup_element.pyx index 2b36616b27e..d64c01efce8 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pyx +++ b/src/sage/groups/perm_gps/permgroup_element.pyx @@ -1298,7 +1298,7 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): return coercion_model.bin_op(left, right, operator.mul) - def _transpose_left(self, j, k): + cpdef PermutationGroupElement _transpose_left(self, int j, int k): if j == k: return self cdef PermutationGroupElement prod = self._new_c() From b15c25fb8c2fe59435e8c1be277320d5e697f6ad Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 3 Aug 2023 23:03:24 -0700 Subject: [PATCH 09/33] src/sage/geometry/palp_normal_form.pyx: Faster access to elements of PM --- src/sage/geometry/palp_normal_form.pyx | 58 +++++++++++++------------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index e772cb9455c..54dd7239391 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -91,35 +91,35 @@ def _palp_PM_max(self, check=False): for k in range(1, n_f): # Error for k == 1 already! permutations[n_s] = [S_f.one(), S_v.one()] - max_element = PM[k, ( permutations[n_s][1])(1) - 1] + max_element = PM.get_unsafe_int(k, ( permutations[n_s][1])(1) - 1) m = 0 for j in range(1, n_v): - element = PM[k, ( permutations[n_s][1])(j + 1) - 1] + element = PM.get_unsafe_int(k, ( permutations[n_s][1])(j + 1) - 1) if element > max_element: max_element = element m = j if m > 0: permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(1, m + 1) - d = (PM[k, ( permutations[n_s][1])(1) - 1] - - PM[first_row_index, ( permutations[0][1])(1) - 1]) + d = (PM.get_unsafe_int(k, ( permutations[n_s][1])(1) - 1) + - PM.get_unsafe_int(first_row_index, ( permutations[0][1])(1) - 1)) if d < 0: # The largest elt of this row is smaller than largest elt # in 1st row, so nothing to do continue # otherwise: for i in range(1, n_v): - max_element = PM[k, ( permutations[n_s][1])(i + 1) - 1] + max_element = PM.get_unsafe_int(k, ( permutations[n_s][1])(i + 1) - 1) m = i for j in range(i + 1, n_v): - element = PM[k, ( permutations[n_s][1])(j + 1) - 1] + element = PM.get_unsafe_int(k, ( permutations[n_s][1])(j + 1) - 1) if element > max_element: max_element = element m = j if m > i: permutations[n_s][1] = ( permutations[n_s][1])._transpose_left(i + 1, m + 1) if d == 0: - d = (PM[k, ( permutations[n_s][1])(i+1) - 1] - - PM[first_row_index, ( permutations[0][1])(i + 1) - 1]) + d = (PM.get_unsafe_int(k, ( permutations[n_s][1])(i+1) - 1) + - PM.get_unsafe_int(first_row_index, ( permutations[0][1])(i + 1) - 1)) if d < 0: break if d < 0: @@ -138,8 +138,9 @@ def _palp_PM_max(self, check=False): n_s = 1 permutations = {k: permutations[k] for k in permutations if k < n_s} - cdef tuple b = tuple(PM[( permutations[0][0])(1) - 1, - ( permutations[0][1])(j+1) - 1] for j in range(n_v)) + cdef tuple b = tuple(PM.get_unsafe_int(( permutations[0][0])(1) - 1, + ( permutations[0][1])(j+1) - 1) + for j in range(n_v)) # Work out the restrictions the current permutations # place on other permutations as a automorphisms # of the first row @@ -176,23 +177,23 @@ def _palp_PM_max(self, check=False): # between 0 and S(0) for s in range(l, n_f): for j in range(1, S[0]): - v0 = PM[( permutations_bar[n_p][0])(s+1) - 1, - ( permutations_bar[n_p][1])(1) - 1] - vj = PM[( permutations_bar[n_p][0])(s+1) - 1, - ( permutations_bar[n_p][1])(j+1) - 1] + v0 = PM.get_unsafe_int(( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1) + vj = PM.get_unsafe_int(( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(j+1) - 1) if v0 < vj: permutations_bar[n_p][1] = ( permutations_bar[n_p][1])._transpose_left(1, j + 1) if ccf == 0: - l_r[0] = PM[( permutations_bar[n_p][0])(s+1) - 1, - ( permutations_bar[n_p][1])(1) - 1] + l_r[0] = PM.get_unsafe_int(( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1) if s != l: permutations_bar[n_p][0] = ( permutations_bar[n_p][0])._transpose_left(l + 1, s + 1) n_p += 1 ccf = 1 permutations_bar[n_p] = list(permutations[k]) else: - d1 = PM[( permutations_bar[n_p][0])(s+1) - 1, - ( permutations_bar[n_p][1])(1) - 1] + d1 = PM.get_unsafe_int(( permutations_bar[n_p][0])(s+1) - 1, + ( permutations_bar[n_p][1])(1) - 1) d = d1 - l_r[0] if d < 0: # We move to the next line @@ -230,20 +231,20 @@ def _palp_PM_max(self, check=False): s -= 1 # Find the largest value in this symmetry block for j in range(c + 1, h): - vc = PM[( permutations_bar[s][0])(l+1) - 1, - ( permutations_bar[s][1])(c+1) - 1] - vj = PM[( permutations_bar[s][0])(l+1) - 1, - ( permutations_bar[s][1])(j+1) - 1] + vc = PM.get_unsafe_int(( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1) + vj = PM.get_unsafe_int(( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(j+1) - 1) if vc < vj: permutations_bar[s][1] = ( permutations_bar[s][1])._transpose_left(c + 1, j + 1) if ccf == 0: # Set reference and carry on to next permutation - l_r[c] = PM[( permutations_bar[s][0])(l+1) - 1, - ( permutations_bar[s][1])(c+1) - 1] + l_r[c] = PM.get_unsafe_int(( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1) ccf = 1 else: - d1 = PM[( permutations_bar[s][0])(l+1) - 1, - ( permutations_bar[s][1])(c+1) - 1] + d1 = PM.get_unsafe_int(( permutations_bar[s][0])(l+1) - 1, + ( permutations_bar[s][1])(c+1) - 1) d = d1 - l_r[c] if d < 0: n_p -= 1 @@ -272,8 +273,9 @@ def _palp_PM_max(self, check=False): # the restrictions the last worked out # row imposes. c = 0 - M = tuple(PM[( permutations[0][0])(l+1) - 1, - ( permutations[0][1])(j+1) - 1] for j in range(n_v)) + M = tuple(PM.get_unsafe_int(( permutations[0][0])(l+1) - 1, + ( permutations[0][1])(j+1) - 1) + for j in range(n_v)) while c < n_v: s = S[c] + 1 S[c] = c + 1 From 9ba1cd85ab58873dae431aa98f284c1c124f9754 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Fri, 4 Aug 2023 21:22:07 -0700 Subject: [PATCH 10/33] sage.geometry.palp_normal_form: Add _palp_canonical_order, factored out from .lattice_polytope --- src/sage/geometry/lattice_polytope.py | 38 ++++---------- src/sage/geometry/palp_normal_form.pyx | 70 +++++++++++++++++++++++++- 2 files changed, 79 insertions(+), 29 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 1ae8504a306..cd14c9d0ca6 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3155,8 +3155,8 @@ def _palp_PM_max(self, check=False): INPUT: - ``check`` -- Boolean (default: ``False``), whether to return - the permutations leaving the maximal vertex-facet pairing - matrix invariant. + the permutations leaving the maximal vertex-facet pairing + matrix invariant. OUTPUT: @@ -3200,7 +3200,7 @@ def _palp_PM_max(self, check=False): True """ from .palp_normal_form import _palp_PM_max - return _palp_PM_max(self, check) + return _palp_PM_max(self.vertex_facet_pairing_matrix(), check) def npoints(self): r""" @@ -4978,12 +4978,13 @@ def _palp_canonical_order(V, PM_max, permutations): - ``PM_max`` -- the maximal vertex-facet pairing matrix - - ``permutation`` -- the permutations of the vertices yielding + - ``permutations`` -- the permutations of the vertices yielding ``PM_max``. OUTPUT: - The PALP normal form as a :class:`point collection `. + The PALP normal form as a :class:`point collection ` + and a permutation. TESTS:: @@ -4998,32 +4999,15 @@ def _palp_canonical_order(V, PM_max, permutations): M(-1, 0) in 2-d lattice M, (1,3,2,4)) """ - n_v = PM_max.ncols() - S_v = SymmetricGroup(n_v) - p_c = S_v.one() - M_max = [max(row[j] for row in PM_max.rows()) for j in range(n_v)] - S_max = sum(PM_max) - for i in range(n_v): - k = i - for j in range(i + 1, n_v): - if M_max[j] < M_max[k] or \ - (M_max[j] == M_max[k] and S_max[j] < S_max[k]): - k = j - if not k == i: - M_max[i], M_max[k] = M_max[k], M_max[i] - S_max[i], S_max[k] = S_max[k], S_max[i] - p_c = S_v((1 + i, 1 + k), check=False) * p_c - # Create array of possible NFs. - permutations = [p_c * l[1] for l in permutations.values()] + from .palp_normal_form import _palp_canonical_order + Vmatrix = V.column_matrix() - Vs = [(Vmatrix.with_permuted_columns(sig).hermite_form(), sig) - for sig in permutations] - Vmin = min(Vs, key=lambda x:x[0]) Vmodule = V.module() - vertices = [Vmodule(_) for _ in Vmin[0].columns()] + vertices, permutation = _palp_canonical_order(Vmatrix, PM_max, permutations) + vertices = [Vmodule(v) for v in vertices] for v in vertices: v.set_immutable() - return (PointCollection(vertices, Vmodule), Vmin[1]) + return PointCollection(vertices, Vmodule), permutation def _palp_convert_permutation(permutation): diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 54dd7239391..e5efacb6d27 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -1,9 +1,11 @@ from sage.groups.perm_gps.permgroup_element cimport PermutationGroupElement from sage.groups.perm_gps.permgroup_named import SymmetricGroup from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense +from sage.matrix.special import column_matrix +from sage.structure.element import Matrix -def _palp_PM_max(self, check=False): +def _palp_PM_max(Matrix_integer_dense PM, check=False): r""" Compute the permutation normal form of the vertex facet pairing matrix . @@ -64,7 +66,6 @@ def _palp_PM_max(self, check=False): sage: all(results) # long time True """ - cdef Matrix_integer_dense PM = self.vertex_facet_pairing_matrix() cdef int n_v = PM.ncols() cdef int n_f = PM.nrows() S_v = SymmetricGroup(n_v) @@ -293,3 +294,68 @@ def _palp_PM_max(self, check=False): return (PM_max, permutations) else: return PM_max + + +def _palp_canonical_order(vertices, PM_max, permutations): + r""" + Compute the PALP normal form of vertices of a lattice polytope + using auxiliary data computed elsewhere. + + This is a helper function for + :meth:`~sage.geometry.lattice_polytope.LatticePolytopeClass.normal_form` + and should not be called directly. + + Given an iterable of vertices, the maximal vertex-facet pairing matrix + and the permutations realizing this matrix, apply the last part of the + PALP algorithm and return the normal form. + + INPUT: + + - ``vertices`` -- iterable of iterables. The vertices. + + - ``PM_max`` -- the maximal vertex-facet pairing matrix + + - ``permutation`` -- the permutations of the vertices yielding ``PM_max``. + + OUTPUT: + + The PALP normal form as an iterable of integer vectors. + + TESTS:: + + sage: L = lattice_polytope.cross_polytope(2) + sage: V = L.vertices() + sage: PM_max, permutations = L._palp_PM_max(check=True) # optional - sage.groups + sage: from sage.geometry.lattice_polytope import _palp_canonical_order + sage: _palp_canonical_order(V, PM_max, permutations) # optional - sage.groups + (M( 1, 0), + M( 0, 1), + M( 0, -1), + M(-1, 0) + in 2-d lattice M, (1,3,2,4)) + """ + n_v = PM_max.ncols() + S_v = SymmetricGroup(n_v) + p_c = S_v.one() + M_max = [max(row[j] for row in PM_max.rows()) for j in range(n_v)] + S_max = sum(PM_max) + for i in range(n_v): + k = i + for j in range(i + 1, n_v): + if M_max[j] < M_max[k] or \ + (M_max[j] == M_max[k] and S_max[j] < S_max[k]): + k = j + if not k == i: + M_max[i], M_max[k] = M_max[k], M_max[i] + S_max[i], S_max[k] = S_max[k], S_max[i] + p_c = S_v((1 + i, 1 + k), check=False) * p_c + # Create array of possible NFs. + permutations = [p_c * l[1] for l in permutations.values()] + if isinstance(vertices, Matrix): + Vmatrix = vertices + else: + Vmatrix = column_matrix(vertices) + Vs = [(Vmatrix.with_permuted_columns(sig).hermite_form(), sig) + for sig in permutations] + Vmin = min(Vs, key=lambda x: x[0]) + return Vmin[0].columns(), Vmin[1] From a47701b6374e9c894bee01f008267e7af213e3f4 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Fri, 4 Aug 2023 21:43:44 -0700 Subject: [PATCH 11/33] Polyhedron_ZZ.normal_form: New --- src/sage/geometry/lattice_polytope.py | 2 +- src/sage/geometry/polyhedron/base_ZZ.py | 38 +++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index cd14c9d0ca6..fd1bc168bb1 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -5024,7 +5024,7 @@ def _palp_convert_permutation(permutation): OUTPUT: - A :class:`permutation group element `. + A :class:`permutation group element `. EXAMPLES:: diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index 561ed76d70c..f40b88b284b 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -834,3 +834,41 @@ def is_known_summand(poly): decompositions.append((X, Y)) summands += [X, Y] return tuple(decompositions) + + def normal_form(self, algorithm="palp_native", permutation=False): + r""" + EXAMPLES: + + We compute the normal form of the "diamond":: + + sage: d = Polyhedron([(1,0), (0,1), (-1,0), (0,-1)]) + sage: d.vertices() + (A vertex at (-1, 0), + A vertex at (0, -1), + A vertex at (0, 1), + A vertex at (1, 0)) + sage: d.normal_form() + [(1, 0), (0, 1), (0, -1), (-1, 0)] + + It is not possible to compute normal forms for polytopes which do not + span the space:: + + sage: p = Polyhedron([(1,0,0), (0,1,0), (-1,0,0), (0,-1,0)]) + sage: p.normal_form() + Traceback (most recent call last): + ... + ValueError: normal form is not defined for A 2-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices + """ + from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order + + if self.dim() < self.ambient_dim(): + raise ValueError("normal form is not defined for %s" % self) + + PM = self.slack_matrix().transpose() + PM_max, permutations = _palp_PM_max(PM, check=True) + out = _palp_canonical_order(self.vertices(), PM_max, permutations) + + if permutation: + return out + else: + return out[0] From 3411205c29962fac2b1ec7a126157a28bd368cb9 Mon Sep 17 00:00:00 2001 From: Luze Xu Date: Tue, 22 Aug 2023 17:50:32 -0700 Subject: [PATCH 12/33] fix bug --- src/sage/geometry/palp_normal_form.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index e5efacb6d27..1dae6e70e1e 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -82,7 +82,7 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): cdef int element, max_element for j in range(n_v): - m = index_of_max(PM[0, i] for i in range(j, n_v)) + m = index_of_max(PM.get_unsafe_int(0, ( permutations[0][1])(i + 1) - 1) for i in range(j, n_v)) if m > 0: permutations[0][1] = ( permutations[0][1])._transpose_left(j + 1, m + j + 1) From 88ebcf72d3aec71219d55be1687734dd1442066d Mon Sep 17 00:00:00 2001 From: Luze Xu Date: Tue, 22 Aug 2023 18:13:09 -0700 Subject: [PATCH 13/33] add the example --- src/sage/geometry/lattice_polytope.py | 7 +++++++ src/sage/geometry/palp_normal_form.pyx | 12 ++++++++++++ 2 files changed, 19 insertions(+) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index fd1bc168bb1..0138f50a2b3 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3198,6 +3198,13 @@ def _palp_PM_max(self, check=False): ....: for j, i in PMs ) sage: all(results) # long time True + sage: P = Polyhedron([(-4,-6),(-4,-5),(0,0),(1,0),(5,6)]) + sage: P.lattice_polytope()._palp_PM_max() + [9 5 4 0 0] + [6 0 6 5 0] + [1 5 0 0 4] + [0 6 0 1 6] + [0 0 3 5 3] """ from .palp_normal_form import _palp_PM_max return _palp_PM_max(self.vertex_facet_pairing_matrix(), check) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 1dae6e70e1e..d2ffb05c8e1 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -65,6 +65,18 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): ....: for j, i in PMs ) sage: all(results) # long time True + sage: from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order + sage: P = Polyhedron([(-4,-6),(-4,-5),(0,0),(1,0),(5,6)]) + sage: PM = P.slack_matrix().transpose() + sage: _palp_PM_max(PM) + [9 5 4 0 0] + [6 0 6 5 0] + [1 5 0 0 4] + [0 6 0 1 6] + [0 0 3 5 3] + sage: PM_max, permutations = _palp_PM_max(PM, check=True) + sage: _palp_canonical_order(P.vertices(), PM_max, permutations) + ([(1, 0), (0, 0), (2, 4), (1, 5), (-1, 1)], (1,2,3)) """ cdef int n_v = PM.ncols() cdef int n_f = PM.nrows() From 996191cc1c0d3677169cb8bd4d68bd3f7b380767 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 23 Aug 2023 12:33:48 -0700 Subject: [PATCH 14/33] Put new tests in TESTS sections --- src/sage/geometry/lattice_polytope.py | 9 +++++++-- src/sage/geometry/palp_normal_form.pyx | 5 +++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 0138f50a2b3..fc33c78a15f 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3198,8 +3198,13 @@ def _palp_PM_max(self, check=False): ....: for j, i in PMs ) sage: all(results) # long time True - sage: P = Polyhedron([(-4,-6),(-4,-5),(0,0),(1,0),(5,6)]) - sage: P.lattice_polytope()._palp_PM_max() + + TESTS: + + Check that a bug introduced in :issue:`35997` is fixed:: + + sage: P = LatticePolytope([(-4,-6),(-4,-5),(0,0),(1,0),(5,6)]) + sage: P._palp_PM_max() [9 5 4 0 0] [6 0 6 5 0] [1 5 0 0 4] diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index d2ffb05c8e1..a08a9166831 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -65,6 +65,11 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): ....: for j, i in PMs ) sage: all(results) # long time True + + TESTS: + + Check that a bug introduced in :issue:`35997` is fixed:: + sage: from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order sage: P = Polyhedron([(-4,-6),(-4,-5),(0,0),(1,0),(5,6)]) sage: PM = P.slack_matrix().transpose() From d6e433f4eb676110199dd2ffc5639925f35667cc Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 23 Aug 2023 12:58:30 -0700 Subject: [PATCH 15/33] src/sage/geometry/lattice_polytope.py: Update copyright --- src/sage/geometry/lattice_polytope.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index fc33c78a15f..5d4b0a45a41 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -93,8 +93,27 @@ """ # **************************************************************************** -# Copyright (C) 2007-2013 Andrey Novoseltsev -# Copyright (C) 2007-2013 William Stein +# Copyright (C) 2007-2017 Andrey Novoseltsev +# 2007-2013 William Stein +# 2009 Mike Hansen +# 2009-2020 John H. Palmieri +# 2010-2014 Volker Braun +# 2012 Samuel Gonshaw +# 2013 Jan Keitel +# 2014 André Apitzsch +# 2014 Wilfried Luebbe +# 2015-2022 Frédéric Chapoton +# 2015 Ursula Whitcher +# 2016 Jori Mäntysalo +# 2017 Travis Scrimshaw +# 2018 Christian Stump +# 2018 Vincent Klein +# 2019 Vincent Delecroix +# 2019 Jonathan Kliem +# 2020 Samuel Lelièvre +# 2021-2023 Matthias Koeppe +# 2022 David Coudert +# 2023 Luze Xu # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by From e241607216a5a0e4f5f7b2d10995ae3ca31acaed Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 23 Aug 2023 12:58:49 -0700 Subject: [PATCH 16/33] src/sage/geometry/palp_normal_form.pyx: Add copyright --- src/sage/geometry/palp_normal_form.pyx | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index a08a9166831..6eca2a56b9b 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -1,3 +1,24 @@ +r""" +PALP normal form of vertices of a lattice polytope +""" +# **************************************************************************** +# Copyright (C) 2013 Jan Keitel +# 2014 Volker Braun +# 2018 Christian Stump +# 2019 Vincent Delecroix +# 2019 Jonathan Kliem +# 2021 Michael Orlitzky +# 2018-2022 Frédéric Chapoton +# 2023 Luze Xu +# 2023 Matthias Koeppe +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + from sage.groups.perm_gps.permgroup_element cimport PermutationGroupElement from sage.groups.perm_gps.permgroup_named import SymmetricGroup from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense @@ -7,8 +28,7 @@ from sage.structure.element import Matrix def _palp_PM_max(Matrix_integer_dense PM, check=False): r""" - Compute the permutation normal form of the vertex facet pairing - matrix . + Compute the permutation normal form of the vertex facet pairing matrix. The permutation normal form of a matrix is defined as the lexicographic maximum under all permutations of its rows and columns. For more From 6f8769be8e3c3c499d9ef42e20b896c7eaf34435 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Sat, 12 Aug 2023 13:18:33 -0700 Subject: [PATCH 17/33] src/sage/geometry/palp_normal_form.pyx: One more use of _transpose_left --- src/sage/geometry/palp_normal_form.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 6eca2a56b9b..92bbd9a362d 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -385,7 +385,7 @@ def _palp_canonical_order(vertices, PM_max, permutations): if not k == i: M_max[i], M_max[k] = M_max[k], M_max[i] S_max[i], S_max[k] = S_max[k], S_max[i] - p_c = S_v((1 + i, 1 + k), check=False) * p_c + p_c = p_c._transpose_left(1 + i, 1 + k) # Create array of possible NFs. permutations = [p_c * l[1] for l in permutations.values()] if isinstance(vertices, Matrix): From 48348b5a490069825df3fdc9f9c649ab902ff26d Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 23 Aug 2023 13:41:34 -0700 Subject: [PATCH 18/33] src/sage/geometry/palp_normal_form.pyx: Update # needs --- src/sage/geometry/palp_normal_form.pyx | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 92bbd9a362d..212e7736023 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.groups r""" PALP normal form of vertices of a lattice polytope """ @@ -56,8 +57,8 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): sage: o = lattice_polytope.cross_polytope(2) sage: PM = o.vertex_facet_pairing_matrix() - sage: PM_max = PM.permutation_normal_form() # optional - sage.graphs - sage: PM_max == o._palp_PM_max() # optional - sage.graphs + sage: PM_max = PM.permutation_normal_form() + sage: PM_max == o._palp_PM_max() True sage: P2 = ReflexivePolytope(2, 0) sage: PM_max, permutations = P2._palp_PM_max(check=True) @@ -362,9 +363,9 @@ def _palp_canonical_order(vertices, PM_max, permutations): sage: L = lattice_polytope.cross_polytope(2) sage: V = L.vertices() - sage: PM_max, permutations = L._palp_PM_max(check=True) # optional - sage.groups + sage: PM_max, permutations = L._palp_PM_max(check=True) sage: from sage.geometry.lattice_polytope import _palp_canonical_order - sage: _palp_canonical_order(V, PM_max, permutations) # optional - sage.groups + sage: _palp_canonical_order(V, PM_max, permutations) (M( 1, 0), M( 0, 1), M( 0, -1), From 63abe415a738385c65400fa4d9783298a8035ce0 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 23 Aug 2023 15:51:43 -0700 Subject: [PATCH 19/33] PermutationGroupElement._transpose_left: Generalize --- .../groups/perm_gps/permgroup_element.pxd | 2 +- .../groups/perm_gps/permgroup_element.pyx | 25 ++++++++++++++++++- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/sage/groups/perm_gps/permgroup_element.pxd b/src/sage/groups/perm_gps/permgroup_element.pxd index b013451d5a0..21b509f06d1 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pxd +++ b/src/sage/groups/perm_gps/permgroup_element.pxd @@ -18,7 +18,7 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): cpdef _set_permutation_group_element(self, PermutationGroupElement p, bint convert) cpdef _mul_(self, other) - cpdef PermutationGroupElement _transpose_left(self, int j, int k) + cpdef PermutationGroupElement _transpose_left(self, j, k) cpdef PermutationGroupElement _generate_new(self, list new_list) cpdef PermutationGroupElement _generate_new_GAP(self, old) cpdef _gap_list(self) diff --git a/src/sage/groups/perm_gps/permgroup_element.pyx b/src/sage/groups/perm_gps/permgroup_element.pyx index d64c01efce8..f25db93010f 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pyx +++ b/src/sage/groups/perm_gps/permgroup_element.pyx @@ -1298,13 +1298,36 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): return coercion_model.bin_op(left, right, operator.mul) - cpdef PermutationGroupElement _transpose_left(self, int j, int k): + cpdef PermutationGroupElement _transpose_left(self, j, k): + r""" + Return the product of the transposition `(j, k)` with ``self``. + + EXAMPLES:: + + sage: S = SymmetricGroup(5) + sage: s = S([5, 2, 4, 3, 1]) + sage: s._transpose_left(2, 3) + (1,5)(2,4,3) + sage: S((2, 3)) * s + (1,5)(2,4,3) + + sage: S = SymmetricGroup(["a", "b", "c", "d", "e"]) + sage: s = S(["e", "b", "d", "c", "a"]) + sage: s._transpose_left("b", "c") + ('a','e')('b','d','c') + sage: S(("b", "c")) * s + ('a','e')('b','d','c') + """ if j == k: return self cdef PermutationGroupElement prod = self._new_c() cdef int i for i from 0 <= i < self.n: prod.perm[i] = self.perm[i] + if not self._parent._has_natural_domain(): + convert_dict = self._parent._domain_to_gap + j = convert_dict[j] + k = convert_dict[k] prod.perm[j - 1], prod.perm[k - 1] = self.perm[k - 1], self.perm[j - 1] return prod From 38da660c064a64ef1aa09b7c085d1a46ac105f81 Mon Sep 17 00:00:00 2001 From: Luze Xu Date: Wed, 23 Aug 2023 14:08:38 -0700 Subject: [PATCH 20/33] remove index_of_max --- src/sage/geometry/palp_normal_form.pyx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index 212e7736023..fad9f968c06 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -109,18 +109,19 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): S_v = SymmetricGroup(n_v) S_f = SymmetricGroup(n_f) - # and find all the ways of making the first row of PM_max - def index_of_max(iterable): - # returns the index of max of any iterable - return max(enumerate(iterable), key=lambda x: x[1])[0] - cdef int n_s = 1 cdef dict permutations = {0: [S_f.one(), S_v.one()]} cdef int j, k, m, d cdef int element, max_element for j in range(n_v): - m = index_of_max(PM.get_unsafe_int(0, ( permutations[0][1])(i + 1) - 1) for i in range(j, n_v)) + max_element = PM.get_unsafe_int(0, ( permutations[0][1])(j + 1) - 1) + m = 0 + for i in range(j + 1, n_v): + element = PM.get_unsafe_int(0, ( permutations[0][1])(i + 1) - 1) + if element > max_element: + max_element = element + m = i - j if m > 0: permutations[0][1] = ( permutations[0][1])._transpose_left(j + 1, m + j + 1) From e5025aea8fc425abdd2570e3c9195a79a6d79414 Mon Sep 17 00:00:00 2001 From: Luze Xu Date: Tue, 5 Sep 2023 10:38:59 -0700 Subject: [PATCH 21/33] crash example --- src/sage/geometry/lattice_polytope.py | 63 +++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 5d4b0a45a41..4e30dbd242f 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3032,6 +3032,69 @@ def normal_form(self, algorithm="palp", permutation=False): M( 0, -1), M(-1, 0) in 2-d lattice M + + Note that the default algorithm may crash for higher dimensions:: + + sage: P = LatticePolytope([[-3, -3, -6, -6, -1], [3, 3, 6, 6, 1], [-3, -3, -6, -6, 1], + ....: [-3, -3, -3, -6, 0], [-3, -3, -3, 0, 0], [-3, -3, 0, 0, 0], + ....: [-3, 0, -6, -6, 0], [-3, 0, -3, -6, 0], [-3, 0, -3, 0, 0], + ....: [-3, 0, 0, 0, -1], [3, 3, 6, 6, -1], [-3, 0, 0, 0, 1], + ....: [0, -3, -6, -6, 0], [0, -3, -3, -6, 0], [0, -3, -3, 0, 0], + ....: [0, -3, 0, 0, -1], [3, 3, 3, 6, 0], [0, -3, 0, 0, 1], + ....: [0, 0, -6, -6, 0], [0, 0, -3, -6, -1], [3, 3, 3, 0, 0], + ....: [0, 0, -3, -6, 1], [0, 0, -3, 0, -1], [3, 3, 0, 0, 0], + ....: [0, 0, -3, 0, 1], [0, 0, 3, 0, -1], [3, 0, 6, 6, 0], + ....: [0, 0, 3, 0, 1], [0, 0, 3, 6, -1], [3, 0, 3, 6, 0], + ....: [0, 0, 3, 6, 1], [0, 0, 6, 6, 0], [0, 3, 0, 0, -1], + ....: [3, 0, 3, 0, 0], [0, 3, 0, 0, 1], [0, 3, 3, 0, 0], + ....: [0, 3, 3, 6, 0], [0, 3, 6, 6, 0], [3, 0,0, 0, -1], [3, 0, 0, 0, 1]]) + sage: P.normal_form(algorithm="palp") + Traceback (most recent call last): + ... + Output: + b'*** stack smashing detected ***: terminated\nAborted\n' + sage: P.normal_form(algorithm="palp_native") + M( 6, 0, 0, 0, 0), + M( -6, 0, 0, 0, 0), + M( 0, 1, 0, 0, 0), + M( 0, 0, 3, 0, 0), + M( 0, 1, 0, 3, 0), + M( 0, 0, 0, 0, 3), + M( -6, 1, 6, 3, -6), + M( -6, 0, 6, 0, -3), + M(-12, 1, 6, 3, -3), + M( -6, 1, 0, 3, 0), + M( -6, 0, 3, 3, 0), + M( 6, 0, -6, -3, 6), + M(-12, 1, 6, 3, -6), + M(-12, 0, 9, 3, -6), + M( 0, 0, 0, -3, 0), + M(-12, 1, 6, 6, -6), + M(-12, 0, 6, 3, -3), + M( 0, 1, -3, 0, 0), + M( 0, 0, -3, -3, 3), + M( 0, 1, 0, 3, -3), + M( 0, -1, 0, -3, 3), + M( 0, 0, 3, 3, -3), + M( 0, -1, 3, 0, 0), + M( 12, 0, -6, -3, 3), + M( 12, -1, -6, -6, 6), + M( 0, 0, 0, 3, 0), + M( 12, 0, -9, -3, 6), + M( 12, -1, -6, -3, 6), + M( -6, 0, 6, 3, -6), + M( 6, 0, -3, -3, 0), + M( 6, -1, 0, -3, 0), + M(-12, 1, 9, 6, -6), + M( 6, 0, -6, 0, 3), + M( 6, -1, -6, -3, 6), + M( 0, 0, 0, 0, -3), + M( 0, -1, 0, -3, 0), + M( 0, 0, -3, 0, 0), + M( 0, -1, 0, 0, 0), + M( 12, -1, -9, -6, 6), + M( 12, -1, -6, -3, 3) + in 5-d lattice M """ if self.dim() < self.lattice_dim(): raise ValueError("normal form is not defined for %s" % self) From 4e06ea90547c481546092b680ed83d250498df3d Mon Sep 17 00:00:00 2001 From: Luze Xu Date: Tue, 5 Sep 2023 17:18:59 -0700 Subject: [PATCH 22/33] speed tests for palp_native --- src/sage/geometry/lattice_polytope.py | 96 ++++++++++++++++++++++++- src/sage/geometry/polyhedron/base_ZZ.py | 10 +++ 2 files changed, 103 insertions(+), 3 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 4e30dbd242f..4f7e21c4df2 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -2950,10 +2950,11 @@ def normal_form(self, algorithm="palp", permutation=False): - ``algorithm`` -- (default: "palp") The algorithm which is used to compute the normal form. Options are: - * "palp" -- Run external PALP code, usually the fastest option. + * "palp" -- Run external PALP code, usually the fastest option, + but may fail in higher dimensions. * "palp_native" -- The original PALP algorithm implemented - in sage. Currently considerably slower than PALP. + in sage. Currently competitive with PALP in many cases. * "palp_modified" -- A modified version of the PALP algorithm which determines the maximal vertex-facet @@ -3033,7 +3034,49 @@ def normal_form(self, algorithm="palp", permutation=False): M(-1, 0) in 2-d lattice M - Note that the default algorithm may crash for higher dimensions:: + The following examples demonstrate the speed improvement of ``"palp_native"``. + In low dimensions, ``"palp_native"`` is the fastest. + As the dimension increases, ``"palp"`` is relatively faster than ``"palp_native"``. + ``"palp_native"`` is usually much faster than ``"palp_modified"``. + But in some cases when the polytope has high symmetry, however, ``"palp_native"`` is slower:: + + sage: o = lattice_polytope.cross_polytope(2) + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + 625 loops, best of 3: 3.07 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + 625 loops, best of 3: 0.445 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + 625 loops, best of 3: 5.01 ms per loop + + sage: o = lattice_polytope.cross_polytope(3) + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + 625 loops, best of 3: 3.22 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + 625 loops, best of 3: 2.73 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + 625 loops, best of 3: 20.7 ms per loop + + sage: o = lattice_polytope.cross_polytope(4) + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + 625 loops, best of 3: 4.84 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + 625 loops, best of 3: 55.6 ms per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + 625 loops, best of 3: 129 ms per loop + + sage: o = lattice_polytope.cross_polytope(5) + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + 10 loops, best of 3: 0.0364 s per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + 10 loops, best of 3: 1.68 s per loop + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + 10 loops, best of 3: 0.858 s per loop + + Note that the default algorithm ``"palp"`` may crash for higher dimensions because of + the overflow errors as mentioned in :issue:`13525#comment:9`. + Then use ``"palp_native"`` instead, which is usually faster than ``"palp_modified"``. + Below is an example where ``"palp"`` fails, and + ``"palp_native"`` is much faster than ``"palp_modified"``:: sage: P = LatticePolytope([[-3, -3, -6, -6, -1], [3, 3, 6, 6, 1], [-3, -3, -6, -6, 1], ....: [-3, -3, -3, -6, 0], [-3, -3, -3, 0, 0], [-3, -3, 0, 0, 0], @@ -3051,6 +3094,7 @@ def normal_form(self, algorithm="palp", permutation=False): sage: P.normal_form(algorithm="palp") Traceback (most recent call last): ... + RuntimeError: Error executing ... for a polytope sequence! Output: b'*** stack smashing detected ***: terminated\nAborted\n' sage: P.normal_form(algorithm="palp_native") @@ -3095,6 +3139,52 @@ def normal_form(self, algorithm="palp", permutation=False): M( 12, -1, -9, -6, 6), M( 12, -1, -6, -3, 3) in 5-d lattice M + sage: P.normal_form(algorithm="palp_modified") + M( 6, 0, 0, 0, 0), + M( -6, 0, 0, 0, 0), + M( 0, 1, 0, 0, 0), + M( 0, 0, 3, 0, 0), + M( 0, 1, 0, 3, 0), + M( 0, 0, 0, 0, 3), + M( -6, 1, 6, 3, -6), + M( -6, 0, 6, 0, -3), + M(-12, 1, 6, 3, -3), + M( -6, 1, 0, 3, 0), + M( -6, 0, 3, 3, 0), + M( 6, 0, -6, -3, 6), + M(-12, 1, 6, 3, -6), + M(-12, 0, 9, 3, -6), + M( 0, 0, 0, -3, 0), + M(-12, 1, 6, 6, -6), + M(-12, 0, 6, 3, -3), + M( 0, 1, -3, 0, 0), + M( 0, 0, -3, -3, 3), + M( 0, 1, 0, 3, -3), + M( 0, -1, 0, -3, 3), + M( 0, 0, 3, 3, -3), + M( 0, -1, 3, 0, 0), + M( 12, 0, -6, -3, 3), + M( 12, -1, -6, -6, 6), + M( 0, 0, 0, 3, 0), + M( 12, 0, -9, -3, 6), + M( 12, -1, -6, -3, 6), + M( -6, 0, 6, 3, -6), + M( 6, 0, -3, -3, 0), + M( 6, -1, 0, -3, 0), + M(-12, 1, 9, 6, -6), + M( 6, 0, -6, 0, 3), + M( 6, -1, -6, -3, 6), + M( 0, 0, 0, 0, -3), + M( 0, -1, 0, -3, 0), + M( 0, 0, -3, 0, 0), + M( 0, -1, 0, 0, 0), + M( 12, -1, -9, -6, 6), + M( 12, -1, -6, -3, 3) + in 5-d lattice M + sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_native") # not tested + 10 loops, best of 3: 0.137 s per loop + sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_modified") # not tested + 10 loops, best of 3: 22.2 s per loop """ if self.dim() < self.lattice_dim(): raise ValueError("normal form is not defined for %s" % self) diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index f40b88b284b..9ce438105f5 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -837,6 +837,10 @@ def is_known_summand(poly): def normal_form(self, algorithm="palp_native", permutation=False): r""" + Return the normal form of vertices of ``self`` if ``self`` is a lattice polytope, + i.e. all vertices have integer coordinates. For more more detail, + see also :meth:`~sage.geometry.lattice_polytope.LatticePolytopeClass.normal_form`. + EXAMPLES: We compute the normal form of the "diamond":: @@ -849,6 +853,12 @@ def normal_form(self, algorithm="palp_native", permutation=False): A vertex at (1, 0)) sage: d.normal_form() [(1, 0), (0, 1), (0, -1), (-1, 0)] + sage: d.lattice_polytope().normal_form("palp_native") + M( 1, 0), + M( 0, 1), + M( 0, -1), + M(-1, 0) + in 2-d lattice M It is not possible to compute normal forms for polytopes which do not span the space:: From 137c8482020819346985d15a8227d10ae53c27fd Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:06:32 -0700 Subject: [PATCH 23/33] Polyhedron_ZZ.normal_form: Reject unbounded polyhedra --- src/sage/geometry/polyhedron/base_ZZ.py | 26 ++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index 9ce438105f5..0a41dc9371a 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -837,9 +837,10 @@ def is_known_summand(poly): def normal_form(self, algorithm="palp_native", permutation=False): r""" - Return the normal form of vertices of ``self`` if ``self`` is a lattice polytope, - i.e. all vertices have integer coordinates. For more more detail, - see also :meth:`~sage.geometry.lattice_polytope.LatticePolytopeClass.normal_form`. + Return the normal form of vertices of the lattice polytope ``self``. + + For more more detail, + see :meth:`~sage.geometry.lattice_polytope.LatticePolytopeClass.normal_form`. EXAMPLES: @@ -867,12 +868,27 @@ def normal_form(self, algorithm="palp_native", permutation=False): sage: p.normal_form() Traceback (most recent call last): ... - ValueError: normal form is not defined for A 2-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices + ValueError: normal form is not defined for lower-dimensional polyhedra, got + A 2-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices + + The normal form is also not defined for unbounded polyhedra:: + + sage: p = Polyhedron(vertices=[[1, 1]], rays=[[1, 0], [0, 1]], base_ring=ZZ) + sage: p.normal_form() + Traceback (most recent call last): + ... + ValueError: normal form is not defined for unbounded polyhedra, got + A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 2 rays + + See :issue:`15280` for proposed extensions to these cases. """ from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order if self.dim() < self.ambient_dim(): - raise ValueError("normal form is not defined for %s" % self) + raise ValueError("normal form is not defined for lower-dimensional polyhedra, got %s" % self) + + if not self.is_compact(): + raise ValueError("normal form is not defined for unbounded polyhedra, got %s" % self) PM = self.slack_matrix().transpose() PM_max, permutations = _palp_PM_max(PM, check=True) From 26c89cfb17a3d2e1bef63a3c99e8ef74d3f28ee8 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:28:57 -0700 Subject: [PATCH 24/33] Docstring/doctest cosmetics --- src/sage/geometry/lattice_polytope.py | 52 +++++++++++++------------- src/sage/geometry/palp_normal_form.pyx | 12 +++--- 2 files changed, 31 insertions(+), 33 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 4f7e21c4df2..a6c53e03f2d 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -2519,7 +2519,7 @@ def index(self): M(-1, 0) in 2-d lattice M - But they are in the same `GL(Z^n)` orbit and have the same + But they are in the same `GL(\ZZ^n)` orbit and have the same normal form:: sage: d.normal_form() # needs palp @@ -2937,26 +2937,26 @@ def normal_form(self, algorithm="palp", permutation=False): Return the normal form of vertices of ``self``. Two full-dimensional lattice polytopes are in the same - ``GL(\ZZ)``-orbit if and only if their normal forms are the + `GL(\ZZ^n)`-orbit if and only if their normal forms are the same. Normal form is not defined and thus cannot be used for polytopes whose dimension is smaller than the dimension of the ambient space. The original algorithm was presented in [KS1998]_ and implemented in PALP. A modified version of the PALP algorithm is discussed in - [GK2013]_ and available here as "palp_modified". + [GK2013]_ and available here as ``"palp_modified"``. INPUT: - - ``algorithm`` -- (default: "palp") The algorithm which is used + - ``algorithm`` -- (default: ``"palp"``) The algorithm which is used to compute the normal form. Options are: - * "palp" -- Run external PALP code, usually the fastest option, + * ``"palp"`` -- Run external PALP code, usually the fastest option, but may fail in higher dimensions. - * "palp_native" -- The original PALP algorithm implemented + * ``"palp_native"`` -- The original PALP algorithm implemented in sage. Currently competitive with PALP in many cases. - * "palp_modified" -- A modified version of the PALP + * ``"palp_modified"`` -- A modified version of the PALP algorithm which determines the maximal vertex-facet pairing matrix first and then computes its automorphisms, while the PALP algorithm does both things @@ -3040,42 +3040,40 @@ def normal_form(self, algorithm="palp", permutation=False): ``"palp_native"`` is usually much faster than ``"palp_modified"``. But in some cases when the polytope has high symmetry, however, ``"palp_native"`` is slower:: + sage: # not tested sage: o = lattice_polytope.cross_polytope(2) - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") 625 loops, best of 3: 3.07 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") 625 loops, best of 3: 0.445 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") 625 loops, best of 3: 5.01 ms per loop - sage: o = lattice_polytope.cross_polytope(3) - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") 625 loops, best of 3: 3.22 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") 625 loops, best of 3: 2.73 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") 625 loops, best of 3: 20.7 ms per loop - sage: o = lattice_polytope.cross_polytope(4) - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") 625 loops, best of 3: 4.84 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") 625 loops, best of 3: 55.6 ms per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") 625 loops, best of 3: 129 ms per loop - sage: o = lattice_polytope.cross_polytope(5) - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp") 10 loops, best of 3: 0.0364 s per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_native") 10 loops, best of 3: 1.68 s per loop - sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") # not tested + sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") 10 loops, best of 3: 0.858 s per loop Note that the default algorithm ``"palp"`` may crash for higher dimensions because of the overflow errors as mentioned in :issue:`13525#comment:9`. Then use ``"palp_native"`` instead, which is usually faster than ``"palp_modified"``. - Below is an example where ``"palp"`` fails, and + Below is an example where ``"palp"`` fails and ``"palp_native"`` is much faster than ``"palp_modified"``:: sage: P = LatticePolytope([[-3, -3, -6, -6, -1], [3, 3, 6, 6, 1], [-3, -3, -6, -6, 1], @@ -3091,7 +3089,7 @@ def normal_form(self, algorithm="palp", permutation=False): ....: [0, 0, 3, 6, 1], [0, 0, 6, 6, 0], [0, 3, 0, 0, -1], ....: [3, 0, 3, 0, 0], [0, 3, 0, 0, 1], [0, 3, 3, 0, 0], ....: [0, 3, 3, 6, 0], [0, 3, 6, 6, 0], [3, 0,0, 0, -1], [3, 0, 0, 0, 1]]) - sage: P.normal_form(algorithm="palp") + sage: P.normal_form(algorithm="palp") # not tested Traceback (most recent call last): ... RuntimeError: Error executing ... for a polytope sequence! @@ -3139,7 +3137,7 @@ def normal_form(self, algorithm="palp", permutation=False): M( 12, -1, -9, -6, 6), M( 12, -1, -6, -3, 3) in 5-d lattice M - sage: P.normal_form(algorithm="palp_modified") + sage: P.normal_form(algorithm="palp_modified") # long time (22s) M( 6, 0, 0, 0, 0), M( -6, 0, 0, 0, 0), M( 0, 1, 0, 0, 0), @@ -3181,9 +3179,9 @@ def normal_form(self, algorithm="palp", permutation=False): M( 12, -1, -9, -6, 6), M( 12, -1, -6, -3, 3) in 5-d lattice M - sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_native") # not tested + sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_native") # not tested 10 loops, best of 3: 0.137 s per loop - sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_modified") # not tested + sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_modified") # not tested 10 loops, best of 3: 22.2 s per loop """ if self.dim() < self.lattice_dim(): diff --git a/src/sage/geometry/palp_normal_form.pyx b/src/sage/geometry/palp_normal_form.pyx index fad9f968c06..fa756fe65b0 100644 --- a/src/sage/geometry/palp_normal_form.pyx +++ b/src/sage/geometry/palp_normal_form.pyx @@ -43,8 +43,8 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): INPUT: - ``check`` -- Boolean (default: ``False``), whether to return - the permutations leaving the maximal vertex-facet pairing - matrix invariant. + the permutations leaving the maximal vertex-facet pairing + matrix invariant. OUTPUT: @@ -80,10 +80,10 @@ def _palp_PM_max(Matrix_integer_dense PM, check=False): ((2,3), (2,3)), ((1,2), (1,2)), ((1,3), (1,3))] - sage: PMs = ( i._palp_PM_max(check=True) - ....: for i in ReflexivePolytopes(2) ) - sage: results = ( len(i) == len(j.automorphisms_of_rows_and_columns()) - ....: for j, i in PMs ) + sage: PMs = (i._palp_PM_max(check=True) + ....: for i in ReflexivePolytopes(2)) + sage: results = (len(i) == len(j.automorphisms_of_rows_and_columns()) + ....: for j, i in PMs) sage: all(results) # long time True From 7759c528fd989fcd1004ba5f9acf613e8fde77d0 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:40:45 -0700 Subject: [PATCH 25/33] Polyhedron_ZZ.normal_form: Add doc --- src/sage/geometry/polyhedron/base_ZZ.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index 0a41dc9371a..06b93aea513 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -839,6 +839,13 @@ def normal_form(self, algorithm="palp_native", permutation=False): r""" Return the normal form of vertices of the lattice polytope ``self``. + INPUT: + + - ``algorithm`` -- must be ``"palp_native"``, the default. + + - ``permutation`` -- boolean (default: ``False``); if ``True``, the permutation + applied to vertices to obtain the normal form is returned as well. + For more more detail, see :meth:`~sage.geometry.lattice_polytope.LatticePolytopeClass.normal_form`. @@ -884,6 +891,9 @@ def normal_form(self, algorithm="palp_native", permutation=False): """ from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order + if algorithm != "palp_native": + raise ValueError("algorithm must be 'palp_native'") + if self.dim() < self.ambient_dim(): raise ValueError("normal form is not defined for lower-dimensional polyhedra, got %s" % self) From ea4d8571ca257417e612429c61540175b7fe15f3 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:41:37 -0700 Subject: [PATCH 26/33] LatticePolytope.normal_form: Change the default to algorithm='palp_native' --- src/sage/geometry/lattice_polytope.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index a6c53e03f2d..f7306d5b3ff 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -2932,7 +2932,7 @@ def nfacets(self): return len(self.facet_normals()) if self.dim() > 0 else 0 @cached_method - def normal_form(self, algorithm="palp", permutation=False): + def normal_form(self, algorithm="palp_native", permutation=False): r""" Return the normal form of vertices of ``self``. @@ -2947,11 +2947,12 @@ def normal_form(self, algorithm="palp", permutation=False): INPUT: - - ``algorithm`` -- (default: ``"palp"``) The algorithm which is used + - ``algorithm`` -- (default: ``"palp_native"``) The algorithm which is used to compute the normal form. Options are: - * ``"palp"`` -- Run external PALP code, usually the fastest option, - but may fail in higher dimensions. + * ``"palp"`` -- Run external PALP code, usually the fastest option + when it works; but reproducible crashes have been observed in dimension + 5 and higher. * ``"palp_native"`` -- The original PALP algorithm implemented in sage. Currently competitive with PALP in many cases. @@ -2962,7 +2963,7 @@ def normal_form(self, algorithm="palp", permutation=False): automorphisms, while the PALP algorithm does both things concurrently. - - ``permutation`` -- (default: ``False``) If ``True`` the permutation + - ``permutation`` -- boolean (default: ``False``); if ``True``, the permutation applied to vertices to obtain the normal form is returned as well. Note that the different algorithms may return different results that nevertheless lead to the same normal form. @@ -2983,7 +2984,7 @@ def normal_form(self, algorithm="palp", permutation=False): M(-1, 0), M( 0, -1) in 2-d lattice M - sage: d.normal_form() # needs palp + sage: d.normal_form() M( 1, 0), M( 0, 1), M( 0, -1), @@ -3034,11 +3035,11 @@ def normal_form(self, algorithm="palp", permutation=False): M(-1, 0) in 2-d lattice M - The following examples demonstrate the speed improvement of ``"palp_native"``. - In low dimensions, ``"palp_native"`` is the fastest. + The following examples demonstrate the speed of the available algorithms. + In low dimensions, the default algorithm, ``"palp_native"``, is the fastest. As the dimension increases, ``"palp"`` is relatively faster than ``"palp_native"``. ``"palp_native"`` is usually much faster than ``"palp_modified"``. - But in some cases when the polytope has high symmetry, however, ``"palp_native"`` is slower:: + In some cases when the polytope has high symmetry, however, ``"palp_native"`` is slower:: sage: # not tested sage: o = lattice_polytope.cross_polytope(2) @@ -3070,7 +3071,7 @@ def normal_form(self, algorithm="palp", permutation=False): sage: %timeit o.normal_form.clear_cache(); o.normal_form("palp_modified") 10 loops, best of 3: 0.858 s per loop - Note that the default algorithm ``"palp"`` may crash for higher dimensions because of + Note that the algorithm ``"palp"`` may crash for higher dimensions because of the overflow errors as mentioned in :issue:`13525#comment:9`. Then use ``"palp_native"`` instead, which is usually faster than ``"palp_modified"``. Below is an example where ``"palp"`` fails and @@ -3195,8 +3196,7 @@ def normal_form(self, algorithm="palp", permutation=False): elif algorithm == "palp_modified": result = self._palp_modified_normal_form(permutation=permutation) else: - raise ValueError('Algorithm must be palp, ' + - 'palp_native, or palp_modified.') + raise ValueError("algorithm must be 'palp', 'palp_native', or 'palp_modified'") if permutation: vertices, perm = result else: @@ -3219,7 +3219,7 @@ def _palp_modified_normal_form(self, permutation=False): INPUT: - - ``permutation`` -- a Boolean, whether to return the permutation of + - ``permutation`` -- boolean (default: ``False``); whether to return the permutation of the order of the vertices that was applied to obtain this matrix. OUTPUT: From 05f82288e41141de58454f3f68474fa8ab8f99eb Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:45:14 -0700 Subject: [PATCH 27/33] Add # needs sage.groups --- src/sage/geometry/lattice_polytope.py | 6 +++--- src/sage/geometry/polyhedron/base_ZZ.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index f7306d5b3ff..c7dffd27341 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -2984,7 +2984,7 @@ def normal_form(self, algorithm="palp_native", permutation=False): M(-1, 0), M( 0, -1) in 2-d lattice M - sage: d.normal_form() + sage: d.normal_form() # needs sage.groups M( 1, 0), M( 0, 1), M( 0, -1), @@ -3096,7 +3096,7 @@ def normal_form(self, algorithm="palp_native", permutation=False): RuntimeError: Error executing ... for a polytope sequence! Output: b'*** stack smashing detected ***: terminated\nAborted\n' - sage: P.normal_form(algorithm="palp_native") + sage: P.normal_form(algorithm="palp_native") # needs sage.groups M( 6, 0, 0, 0, 0), M( -6, 0, 0, 0, 0), M( 0, 1, 0, 0, 0), @@ -3138,7 +3138,7 @@ def normal_form(self, algorithm="palp_native", permutation=False): M( 12, -1, -9, -6, 6), M( 12, -1, -6, -3, 3) in 5-d lattice M - sage: P.normal_form(algorithm="palp_modified") # long time (22s) + sage: P.normal_form(algorithm="palp_modified") # long time (22s) # needs sage.groups M( 6, 0, 0, 0, 0), M( -6, 0, 0, 0, 0), M( 0, 1, 0, 0, 0), diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index 06b93aea513..cd5dbaa1958 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -859,9 +859,9 @@ def normal_form(self, algorithm="palp_native", permutation=False): A vertex at (0, -1), A vertex at (0, 1), A vertex at (1, 0)) - sage: d.normal_form() + sage: d.normal_form() # needs sage.groups [(1, 0), (0, 1), (0, -1), (-1, 0)] - sage: d.lattice_polytope().normal_form("palp_native") + sage: d.lattice_polytope().normal_form("palp_native") # needs sage.groups M( 1, 0), M( 0, 1), M( 0, -1), From 14a0fe4710f34070bea4b7f6b38babc2c72b21e4 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:54:55 -0700 Subject: [PATCH 28/33] src/sage/geometry/polyhedron/base_ZZ.py: Update copyright --- src/sage/geometry/polyhedron/base_ZZ.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index cd5dbaa1958..2f280a4d74d 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -3,7 +3,14 @@ """ # **************************************************************************** -# Copyright (C) 2011 Volker Braun +# Copyright (C) 2011-2013 Volker Braun +# 2015 Nathann Cohen +# 2015 Vincent Delecroix +# 2017-2018 Frédéric Chapoton +# 2019 Sophia Elia +# 2019-2020 Jonathan Kliem +# 2023 Luze Xu +# 2023 Matthias Koeppe # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by From 13dbc7e56c0f538ed936fd4fe7c83be9413cacb7 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Thu, 7 Sep 2023 11:58:28 -0700 Subject: [PATCH 29/33] src/sage/geometry/lattice_polytope.py: Update # needs --- src/sage/geometry/lattice_polytope.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index c7dffd27341..5874d2d7e9a 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -2522,13 +2522,13 @@ def index(self): But they are in the same `GL(\ZZ^n)` orbit and have the same normal form:: - sage: d.normal_form() # needs palp + sage: d.normal_form() # needs sage.groups M( 1, 0), M( 0, 1), M( 0, -1), M(-1, 0) in 2-d lattice M - sage: lattice_polytope.ReflexivePolytope(2,3).normal_form() + sage: lattice_polytope.ReflexivePolytope(2,3).normal_form() # needs sage.groups M( 1, 0), M( 0, 1), M( 0, -1), From e348fd91a5b042bbdc9e440840beae7cc154f185 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Fri, 3 Nov 2023 11:35:29 -0700 Subject: [PATCH 30/33] src/sage/groups/perm_gps/permgroup_element.pyx: use a standard python-style loop --- src/sage/groups/perm_gps/permgroup_element.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/groups/perm_gps/permgroup_element.pyx b/src/sage/groups/perm_gps/permgroup_element.pyx index cb47bf05033..a459451bd6d 100644 --- a/src/sage/groups/perm_gps/permgroup_element.pyx +++ b/src/sage/groups/perm_gps/permgroup_element.pyx @@ -1322,7 +1322,7 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement): return self cdef PermutationGroupElement prod = self._new_c() cdef int i - for i from 0 <= i < self.n: + for i in range(self.n): prod.perm[i] = self.perm[i] if not self._parent._has_natural_domain(): convert_dict = self._parent._domain_to_gap From 8469a642f1bcc9def53f57874345af014d786abd Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Fri, 3 Nov 2023 12:03:24 -0700 Subject: [PATCH 31/33] src/sage/geometry/polyhedron/base_ZZ.py: Add more tests --- src/sage/geometry/polyhedron/base_ZZ.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/sage/geometry/polyhedron/base_ZZ.py b/src/sage/geometry/polyhedron/base_ZZ.py index 14eadcc7ae9..20b3c95bbcc 100644 --- a/src/sage/geometry/polyhedron/base_ZZ.py +++ b/src/sage/geometry/polyhedron/base_ZZ.py @@ -875,6 +875,11 @@ def normal_form(self, algorithm="palp_native", permutation=False): M(-1, 0) in 2-d lattice M + Using ``permutation=True``:: + + sage: d.normal_form(permutation=True) # needs sage.groups + ([(1, 0), (0, 1), (0, -1), (-1, 0)], ()) + It is not possible to compute normal forms for polytopes which do not span the space:: @@ -895,6 +900,13 @@ def normal_form(self, algorithm="palp_native", permutation=False): A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 2 rays See :issue:`15280` for proposed extensions to these cases. + + TESTS:: + + sage: d.normal_form(algorithm="palp_fiction") + Traceback (most recent call last): + ... + ValueError: algorithm must be 'palp_native' """ from sage.geometry.palp_normal_form import _palp_PM_max, _palp_canonical_order From 0bd2a93f84e43f1190714991767ce8cbdf039443 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Wed, 22 Nov 2023 12:08:31 -0800 Subject: [PATCH 32/33] src/sage/geometry/lattice_polytope.py: Add test --- src/sage/geometry/lattice_polytope.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index 5874d2d7e9a..cdf2c374e9f 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3184,6 +3184,13 @@ def normal_form(self, algorithm="palp_native", permutation=False): 10 loops, best of 3: 0.137 s per loop sage: %timeit P.normal_form.clear_cache(); P.normal_form("palp_modified") # not tested 10 loops, best of 3: 22.2 s per loop + + TESTS:: + + sage: d.normal_form("palp_fiction") + Traceback (most recent call last): + ... + ValueError: algorithm must be 'palp', 'palp_native', or 'palp_modified' """ if self.dim() < self.lattice_dim(): raise ValueError("normal form is not defined for %s" % self) From f9aeca57d76136a88b74a0c3297250070c0227a3 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Mon, 27 Nov 2023 12:59:47 -0800 Subject: [PATCH 33/33] src/sage/geometry/lattice_polytope.py: Mark large example as # not tested --- src/sage/geometry/lattice_polytope.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/geometry/lattice_polytope.py b/src/sage/geometry/lattice_polytope.py index cdf2c374e9f..1cd9f578113 100644 --- a/src/sage/geometry/lattice_polytope.py +++ b/src/sage/geometry/lattice_polytope.py @@ -3138,7 +3138,7 @@ def normal_form(self, algorithm="palp_native", permutation=False): M( 12, -1, -9, -6, 6), M( 12, -1, -6, -3, 3) in 5-d lattice M - sage: P.normal_form(algorithm="palp_modified") # long time (22s) # needs sage.groups + sage: P.normal_form(algorithm="palp_modified") # not tested (22s; MemoryError on 32 bit), needs sage.groups M( 6, 0, 0, 0, 0), M( -6, 0, 0, 0, 0), M( 0, 1, 0, 0, 0),