diff --git a/veerer/automaton.py b/veerer/automaton.py index eb964a5..c1732d6 100644 --- a/veerer/automaton.py +++ b/veerer/automaton.py @@ -603,7 +603,7 @@ def run(self, max_size=None): if is_target_valid: # = T is a valid vertex r, _ = T.best_relabelling() rinv = perm_invert(r) - T.relabel(r) + T.relabel(r, check=False) new_state = T.copy(mutable=False) new_state.set_immutable() graph[state_branch[-1]].append((new_state, new_flip)) @@ -635,7 +635,7 @@ def run(self, max_size=None): if self._verbosity >= 2: print('[automaton] known vertex') sys.stdout.flush() - T.relabel(rinv) + T.relabel(rinv, check=False) else: # = T is not a valid vertex if self._verbosity >= 2: @@ -652,7 +652,7 @@ def run(self, max_size=None): while flips and not flip_branch[-1]: rinv = relabellings.pop() - T.relabel(rinv) + T.relabel(rinv, check=False) flip_back_data = flips.pop() self._flip_back(flip_back_data) state_branch.pop() @@ -845,6 +845,10 @@ class GeometricAutomaton(Automaton): sage: GeometricAutomaton(vt) Geometric veering automaton with 54 vertices + sage: from surface_dynamics import AbelianStratum + sage: for comp in AbelianStratum(4).components(): + ....: print(comp, GeometricAutomaton(VeeringTriangulation.from_stratum(comp))) + One can check that the cardinality is indeed correct:: sage: sum(x.is_geometric() for x in CoreAutomaton(vt)) @@ -869,6 +873,8 @@ def _flip(self, flip_data): flip_back_data = tuple((e, self._state.colour(e)) for e, _ in flip_data) for e, col in flip_data: self._state.flip(e, col) + if env.CHECK and not self._state.is_geometric(): + raise RuntimeError('that was indeed possible!') return True, flip_back_data def _flip_back(self, flip_back_data): @@ -896,6 +902,18 @@ class GeometricAutomatonSubspace(Automaton): sage: sum(x.is_geometric() for x in CoreAutomaton(vt)) 54 + + Less trivial examples:: + + sage: vt, s, t = VeeringTriangulations.L_shaped_surface(1, 1, 1, 1) + sage: f = VeeringTriangulationLinearFamily(vt, [s, t]) + sage: GeometricAutomatonSubspace(f) + Geometric veering linear constraint automaton with 6 vertices + + sage: vt, s, t = VeeringTriangulations.L_shaped_surface(2, 3, 4, 5, 1, 2) + sage: f = VeeringTriangulationLinearFamily(vt, [s, t]) + sage: GeometricAutomatonSubspace(f) + """ _name = 'geometric veering linear constraint' @@ -916,6 +934,10 @@ def _flip(self, flip_data): flip_back_data = tuple((e, self._state.colour(e)) for e, _ in flip_data) for e, col in flip_data: self._state.flip(e, col) + if env.CHECK: + self._state._check(RuntimeError) + if not self._state.is_geometric(): + raise RuntimeError return True, flip_back_data def _flip_back(self, flip_back_data): diff --git a/veerer/env.py b/veerer/env.py index 29f2bb9..bc7b11b 100644 --- a/veerer/env.py +++ b/veerer/env.py @@ -6,7 +6,7 @@ sage: from veerer.env import sage, flipper, surface_dynamics, ppl """ -CHECK = True +CHECK = False try: import sage.all diff --git a/veerer/linear_family.py b/veerer/linear_family.py index 3e772ba..872dc58 100644 --- a/veerer/linear_family.py +++ b/veerer/linear_family.py @@ -23,9 +23,11 @@ from array import array from copy import copy +import collections +import itertools from .env import sage, ppl -from .constants import VERTICAL, HORIZONTAL +from .constants import VERTICAL, HORIZONTAL, BLUE, RED from .permutation import perm_cycle_string, perm_cycles, perm_check, perm_conjugate, perm_on_list @@ -274,7 +276,6 @@ def _check(self, error=ValueError): self._set_switch_conditions(self._tt_check, v, VERTICAL) if subspace != subspace.echelon_form(): raise error('subspace not in echelon form') - if self._mutable != self._subspace.is_mutable(): raise error('incoherent mutability states') @@ -503,19 +504,140 @@ def flip_back(self, e, col, check=True): VeeringTriangulation.flip_back(self, e, col, Gx=self._subspace, check=check) self._subspace.echelonize() + def _conjugate_subspace(self): + mat = copy(self._subspace) + ne = self.num_edges() + ep = self._ep + for j in range(ne): + if ep[j] < j: + raise ValueError('not in standard form') + if self._colouring[j] == BLUE: + for i in range(mat.nrows()): + mat[i, j] *= -1 + return mat + def geometric_polytope(self, x_low_bound=0, y_low_bound=0, hw_bound=0): r""" - Return the geometric polyopte + Return the geometric polytope. + + EXAMPLES:: + + sage: from veerer import * + + sage: T = VeeringTriangulation("(0,1,2)(~0,~1,~2)", "RRB") + sage: T.geometric_polytope() + sage: T.as_linear_family().geometric_polytope() """ - return VeeringTriangulation.geometric_polytope(self, Gx=self._subspace, x_low_bound=x_low_bound, y_low_bound=y_low_bound, hw_bound=hw_bound) + dim = self._subspace.ncols() + gens = [tuple(r) + (0,) * dim for r in self._subspace] + gens.extend((0,) * dim + tuple(r) for r in self._conjugate_subspace()) + subspace = Polyhedron(lines=gens) + + ieqs = [] + eqns = [] + for g in VeeringTriangulation.geometric_polytope(self).constraints(): + l = (g.inhomogeneous_term(),) + g.coefficients() + if g.is_equality(): + eqns.append(l) + elif g.is_inequality(): + ieqs.append(l) + return subspace.intersection(Polyhedron(ieqs=ieqs, eqns=eqns)) def geometric_flips(self): r""" - Return the list of geometric neighbours. + Return the list of geometric flips. + + A flip is geometric if it arises generically as a flip of L^oo-Delaunay + triangulations along Teichmueller geodesics. These correspond to a subset + of facets of the geometric polytope. OUTPUT: a list of pairs (edge number, new colour) + + EXAMPLES: + + L-shaped square tiled surface with 3 squares (given as a sphere with + 3 triangles). It has two geometric neighbors corresponding to simultaneous + flipping of the diagonals 3, 4 and 5:: + + sage: from veerer import * + sage: T, s, t = VeeringTriangulations.L_shaped_surface(1, 1, 1, 1) + sage: f = VeeringTriangulationLinearFamily(T, [s, t]) + sage: sorted(T.geometric_flips()) + [[(3, 1), (4, 1), (5, 1)], [(3, 2), (4, 2), (5, 2)]] + + To be compared with the geometric flips in the ambient stratum:: + + sage: sorted(T.as_linear_family().geometric_flips()) + [[(3, 1)], [(3, 2)], [(4, 1)], [(4, 2)], [(5, 1)], [(5, 2)]] + + A more complicated example:: + + sage: T, s, t = VeeringTriangulations.L_shaped_surface(2, 3, 5, 2, 1, 1) + sage: f = VeeringTriangulationLinearFamily(T, [s, t]) + sage: sorted(T.geometric_flips()) + [[(4, 2)], [(5, 1)], [(5, 2)]] + + TESTS:: + + sage: from veerer import VeeringTriangulation + sage: fp = "(0,~8,~7)(1,3,~2)(2,7,~3)(4,6,~5)(5,8,~6)(~4,~1,~0)" + sage: cols = "RBRRRRBBR" + sage: vt = VeeringTriangulation(fp, cols) + sage: sorted(vt.geometric_flips()) + [[(2, 1)], [(2, 2)], [(4, 1), (8, 1)], [(4, 2), (8, 2)]] + sage: sorted(vt.as_linear_family().geometric_flips()) + [[(2, 1)], [(2, 2)], [(4, 1), (8, 1)], [(4, 2), (8, 2)]] """ - return VeeringTriangulation.geometric_flips(self, Gx=self._subspace) + ambient_dim = self._subspace.ncols() + dim = self._subspace.nrows() + P = self.geometric_polytope() + if P.dimension() != 2 * dim: + raise ValueError('not geometric P.dimension() = {} while 2 * dim = {}'.format(P.dimension(), 2 * dim)) + + eqns = P.equations_list() + ieqs = P.inequalities_list() + + delaunay_facets = collections.defaultdict(list) + for e in self.forward_flippable_edges(): + a, b, c, d = self.square_about_edge(e) + + hyperplane = [0] * (1 + 2 * ambient_dim) + hyperplane[1 + self._norm(e)] = 1 + hyperplane[1 + ambient_dim + self._norm(a)] = -1 + hyperplane[1 + ambient_dim + self._norm(d)] = -1 + Q = Polyhedron(base_ring=P.base_ring(), + eqns=eqns + [hyperplane], + ieqs=ieqs) + if Q.dimension() == 2 * dim - 1: + delaunay_facets[Q].append(e) + + neighbours = [] + for Q, edges in delaunay_facets.items(): + eqns = Q.equations_list() + ieqs = Q.inequalities_list() + for cols in itertools.product([BLUE, RED], repeat=len(edges)): + Z = list(zip(edges, cols)) + half_spaces = [] + for e, col in Z: + a, b, c, d = self.square_about_edge(e) + half_space = [0] * (1 + 2 * ambient_dim) + if col == RED: + # x[a] <= x[d] + half_space[1 + self._norm(d)] = +1 + half_space[1 + self._norm(a)] = -1 + else: + # x[a] >= x[d] + half_space[1 + self._norm(d)] = -1 + half_space[1 + self._norm(a)] = +1 + half_spaces.append(half_space) + S = Polyhedron(base_ring=Q.base_ring(), + eqns=eqns, + ieqs=ieqs + half_spaces) + if S.dimension() == 2 * dim - 1: + neighbours.append(Z) + + return neighbours + class VeeringTriangulationLinearFamilies: diff --git a/veerer/permutation.py b/veerer/permutation.py index 27d1015..59165bc 100644 --- a/veerer/permutation.py +++ b/veerer/permutation.py @@ -626,6 +626,48 @@ def perm_cycles(p, singletons=True, n=None): return res +def perm_cycles_lengths(p, n=None): + r""" + Return the array of orbit sizes. + + INPUT: + + - ``p`` -- the permutation + + - ``n`` -- (optional) only use the first ``n`` elements of the permutation ``p`` + + EXAMPLES:: + + sage: from veerer.permutation import perm_cycles_lengths + + sage: perm_cycles_lengths([0,2,1]) + [1, 2, 2] + sage: perm_cycles_lengths([2,-1,0]) + [2, -1, 2] + sage: perm_cycles_lengths([2,0,1,None,None], n=3) + [3, 3, 3] + """ + if n is None: + n = len(p) + elif n < 0 or n > len(p): + raise ValueError + + lengths = [-1] * n + for i in range(n): + if lengths[i] != -1 or p[i] == -1: + continue + j = i + m = 0 + while lengths[j] == -1: + lengths[j] = 0 + j = p[j] + m += 1 + while lengths[j] == 0: + lengths[j] = m + j = p[j] + + return lengths + def perm_num_cycles(p, n=None): r""" Return the number of cycles of the permutation ``p``. diff --git a/veerer/triangulation.py b/veerer/triangulation.py index 56221eb..c14c888 100644 --- a/veerer/triangulation.py +++ b/veerer/triangulation.py @@ -908,7 +908,7 @@ def flip_homological_action(self, e, m, twist=False): else: m[e] = (m[a] if a < A else -m[A]) + (m[d] if d < D else -m[D]) - def relabel_homological_action(self, p, m, twist=False): + def relabel_homological_action(self, p, m, twist=False, check=True): r""" Acts on the homological vectors ``m`` by the relabeling permutation ``p``. @@ -946,7 +946,7 @@ def relabel_homological_action(self, p, m, twist=False): """ n = self._n ne = self.num_edges() - if not perm_check(p, n): + if check and not perm_check(p, n): p = perm_init(p, self._n, self._ep) if not perm_check(p, n): raise ValueError('invalid relabeling permutation') @@ -1159,10 +1159,7 @@ def swap(self, e): vp[e] = E_vp vp[E] = e_vp - # TODO: remove check - self._check() - - def relabel(self, p): + def relabel(self, p, check=True): r""" Relabel this triangulation inplace according to the permutation ``p``. @@ -1200,7 +1197,7 @@ def relabel(self, p): raise ValueError('immutable triangulation; use a mutable copy instead') n = self._n - if not perm_check(p, n): + if check and not perm_check(p, n): # if the input is not a valid permutation, we assume that half-edges # are not separated p = perm_init(p, n, self._ep) @@ -1211,7 +1208,6 @@ def relabel(self, p): self._ep = perm_conjugate(self._ep, p) self._fp = perm_conjugate(self._fp, p) - def flip(self, e): r""" Flip the edge ``e``. @@ -1551,30 +1547,11 @@ def _automorphism_good_starts(self): # vp, ep, fp and vp[ep[fp]] n = self._n - lv = [-1] * n - for c in perm_cycles(self._vp): - m = len(c) - for i in c: - lv[i] = m - - le = [-1] * n - for c in perm_cycles(self._ep): - m = len(c) - for i in c: - le[i] = m - - lf = [-1] * n - for c in perm_cycles(self._fp): - m = len(c) - for i in c: - lf[i] = m - - lvef = [-1] * n + lv = perm_cycles_lengths(self._vp, n) + le = perm_cycles_lengths(self._ep, n) + lf = perm_cycles_lengths(self._fp, n) vef = perm_compose(perm_compose(self._fp, self._ep), self._vp) - for c in perm_cycles(vef): - m = len(c) - for i in c: - lvef[i] = m + lvef = perm_cycles_lengths(vef, n) d = {} for i in range(n): @@ -1584,7 +1561,7 @@ def _automorphism_good_starts(self): else: d[s] = [i] m = min(len(x) for x in d.values()) - candidates = [s for s,v in d.items() if len(v) == m] + candidates = [s for s, v in d.items() if len(v) == m] winner = min(candidates) return d[winner] @@ -1726,7 +1703,7 @@ def set_canonical_labels(self): raise ValueError('immutable triangulation; use a mutable copy instead') r, _ = self.best_relabelling() - self.relabel(r) + self.relabel(r, check=False) def iso_sig(self): r""" diff --git a/veerer/veering_triangulation.py b/veerer/veering_triangulation.py index 47ec2c9..dc3080e 100644 --- a/veerer/veering_triangulation.py +++ b/veerer/veering_triangulation.py @@ -1027,7 +1027,7 @@ def abelian_cover(self, mutable=None): colouring_cov = self._colouring * 2 # TODO: remove the check argument vt = self.from_face_edge_perms(colouring_cov, array('l', fp_cov), array('l', ep_cov), mutable=True, check=True) - vt.relabel(vt._relabelling_from(0)) + vt.relabel(vt._relabelling_from(0), check=False) if not mutable: vt.set_immutable() return vt @@ -1396,7 +1396,7 @@ def relabel(self, p, check=True): if not perm_check(p, n): raise ValueError('invalid relabelling permutation') - Triangulation.relabel(self, p) + Triangulation.relabel(self, p, check=False) perm_on_list(p, self._colouring) def _automorphism_good_starts(self, has_purple=None, has_green=None): @@ -1612,6 +1612,7 @@ def best_relabelling(self, all=False): best = None if all: relabellings = [] + for start_edge in self._automorphism_good_starts(): relabelling = self._relabelling_from(start_edge) @@ -2518,7 +2519,6 @@ def switch_matrix(self, slope=VERTICAL): EXAMPLES:: sage: from veerer import VeeringTriangulation, HORIZONTAL - sage: import ppl sage: vt = VeeringTriangulation("(0,1,2)(3,4,5)(6,7,8)(~8,~0,~7)(~6,~1,~5)(~4,~2,~3)", "RRBRRBRRB") sage: vt.switch_matrix() [ 1 -1 1 0 0 0 0 0 0] @@ -2535,7 +2535,6 @@ def switch_matrix(self, slope=VERTICAL): [ 0 1 0 0 0 1 -1 0 0] [ 0 0 1 -1 1 0 0 0 0] """ - require_package('ppl', 'switch_matrix') require_package('sage', 'switch_matrix') from sage.matrix.constructor import matrix @@ -2928,7 +2927,6 @@ def _linear_subspace_generator_system(self, x, y, Gx): gs.insert(ppl.line(g)) return gs - # TODO: ppl only allows *rational* polytope which is not enough for our purpose def geometric_polytope(self, Lx=None, Gx=None, x_low_bound=0, y_low_bound=0, hw_bound=0): r""" Return the geometric polytope of this veering triangulation. @@ -2975,9 +2973,10 @@ def geometric_polytope(self, Lx=None, Gx=None, x_low_bound=0, y_low_bound=0, hw_ x = [ppl.Variable(e) for e in range(ne)] y = [ppl.Variable(ne+e) for e in range(ne)] + if Lx is None and Gx is None: - P = self.train_track_polytope(VERTICAL, low_bound=x_low_bound) - P.concatenate_assign(self.train_track_polytope(HORIZONTAL, low_bound=y_low_bound)) + P = VeeringTriangulation.train_track_polytope(self, VERTICAL, low_bound=x_low_bound) + P.concatenate_assign(VeeringTriangulation.train_track_polytope(self, HORIZONTAL, low_bound=y_low_bound)) self._set_geometric_constraints(P.add_constraint, x, y, hw_bound=hw_bound) return P @@ -3642,7 +3641,8 @@ def geometric_flips(self, Lx=None, Gx=None): else: dim = self.stratum_dimension() - P = self.geometric_polytope(Lx=Lx, Gx=Gx) + # NOTE: geometric_polytope might be redefined in subclasses + P = VeeringTriangulation.geometric_polytope(self, Lx=Lx, Gx=Gx) if P.affine_dimension() != 2 * dim: raise ValueError('not geometric or invalid constraints') @@ -3654,19 +3654,23 @@ def geometric_flips(self, Lx=None, Gx=None): Q = ppl.C_Polyhedron(P) Q.add_constraint(x[self._norm(e)] == y[self._norm(a)] + y[self._norm(d)]) - if Q.affine_dimension() == 2*dim - 1: + facet_dim = Q.affine_dimension() + assert facet_dim < 2 * dim + if facet_dim == 2*dim - 1: hQ = ppl_cone_to_hashable(Q) if hQ not in delaunay_facets: delaunay_facets[hQ] = [Q, []] delaunay_facets[hQ][1].append(e) # computing colouring of the new triangulations + # NOTE: is it really enough? We need the neighbour to also be geometric + # (and not 2*dim - 1 dimensional) neighbours = [] for Q, edges in delaunay_facets.values(): for cols in product([BLUE, RED], repeat=len(edges)): S = ppl.C_Polyhedron(Q) Z = list(zip(edges, cols)) - for e,col in Z: + for e, col in Z: a,b,c,d = self.square_about_edge(e) if col == RED: S.add_constraint(x[self._norm(a)] <= x[self._norm(d)])