diff --git a/veerer/layout.py b/veerer/layout.py index 67373f8..b2c5cc9 100644 --- a/veerer/layout.py +++ b/veerer/layout.py @@ -785,10 +785,15 @@ def _plot_edge_label(self, a, tilde=None, **opts): if tilde: if a > ep[a]: - lab = "%s=~%s" % (a, ep[a]) + #lab = "%s=~%s" % (a, ep[a]) + # nicer shorter version + lab = "~" + str(ep[a]) else: lab = str(a) else: + if a > ep[a]: + from sage.plot.graphics import Graphics + return Graphics() lab = str(a) x, y = self._triangulation._holonomies[a] diff --git a/veerer/linear_family.py b/veerer/linear_family.py index 90c393d..dd3f893 100644 --- a/veerer/linear_family.py +++ b/veerer/linear_family.py @@ -310,6 +310,17 @@ def _horizontal_subspace(self): mat[i, j] *= -1 return mat + def switch_subspace_generators_matrix(self, slope=VERTICAL): + if slope == VERTICAL: + if not self._mutable: + return self._subspace + else: + return self._subspace.__copy__() + elif slope == HORIZONTAL: + return self._horizontal_subspace + else: + raise ValueError + def as_linear_family(self): return self diff --git a/veerer/permutation.pyx b/veerer/permutation.pyx index f24a85d..83336ac 100644 --- a/veerer/permutation.pyx +++ b/veerer/permutation.pyx @@ -374,6 +374,40 @@ def str_to_cycles(s): return r +def str_to_cycles_and_data(s): + r""" + Return a list of cycles and data from a string. + + EXAMPLES:: + + sage: from veerer.permutation import str_to_cycles_and_data + sage: str_to_cycles_and_data('(0:1,1:2)') + ([[0, 1]], {0: '1', 1: '2'}) + sage: str_to_cycles_and_data('(0:1,1:2)(3:0)') + ([[0, 1], [3]], {0: 1, 1: 2, 3: 0}) + """ + r = [] + data = {} + for c_str in s[1:-1].split(')('): + if not c_str: + continue + cycle = [] + for c in c_str.replace(' ', '').split(','): + i = c.find(':') + if i == -1 or c.find(':', i + 1) != -1: + raise ValueError('invalid input string') + j = c[i+1:] + c = c[:i] + if c[0] == '~': + i = ~int(c[1:]) + else: + i = int(c) + cycle.append(i) + data[i] = int(j) + r.append(cycle) + return r, data + + def perm_random(int n): r""" Return a random permutation. @@ -788,6 +822,31 @@ def perm_cycle_type(array.array p, int n=-1): return c +def perm_cycles_to_string(list cycles, involution=None, data=None): + r""" + Return a string representing a list of cycles. + + INPUT: + + - ``cycles`` -- list of cycles + + - ``involution`` -- optional involution (possibly with fixed points) + + - ``data`` -- optional data + """ + if involution: + if data: + elt = lambda e: ('%d' % e if e <= involution[e] else '~%d' % involution[e]) + (':%d' % data[e]) + else: + elt = lambda e: ('%d' % e if e <= involution[e] else '~%d' % involution[e]) + elif data: + elt = lambda e: ('%d' % e) + (':%d' % data[e]) + else: + elt = str + + return ''.join(map(lambda x: '(' + ','.join(map(elt, x)) + ')', cycles)) + + def perm_cycle_string(array.array p, singletons=True, n=-1, involution=None): r""" Return a string representing the cycle decomposition of `p` @@ -802,13 +861,7 @@ def perm_cycle_string(array.array p, singletons=True, n=-1, involution=None): sage: perm_cycle_string(array('i', [0,2,1]), False) '(1,2)' """ - if involution: - elt = lambda e: '%d'%e if e <= involution[e] else '~%d'%involution[e] - else: - elt = str - - return ''.join(map(lambda x: '('+','.join(map(elt, x))+')', - perm_cycles(p, singletons, n))) + return perm_cycles_to_string(perm_cycles(p, singletons, n), involution) def perm_orbit(array.array p, int i): diff --git a/veerer/polyhedron/linear_algebra.py b/veerer/polyhedron/linear_algebra.py index a7b0577..c408b4e 100644 --- a/veerer/polyhedron/linear_algebra.py +++ b/veerer/polyhedron/linear_algebra.py @@ -99,3 +99,6 @@ def linear_form_normalize(base_ring, linear_form): linear_form[k] /= c return linear_form + + +vector_normalize = linear_form_normalize diff --git a/veerer/triangulation.py b/veerer/triangulation.py index a278dab..b86877f 100644 --- a/veerer/triangulation.py +++ b/veerer/triangulation.py @@ -30,13 +30,16 @@ from .permutation import (perm_init, perm_check, perm_cycles, perm_dense_cycles, perm_invert, perm_conjugate, perm_cycle_string, perm_cycles_lengths, - perm_num_cycles, str_to_cycles, perm_compose, perm_from_base64_str, + perm_cycles_to_string, perm_on_list, + perm_num_cycles, str_to_cycles, str_to_cycles_and_data, perm_compose, perm_from_base64_str, uint_base64_str, uint_from_base64_str, perm_base64_str, perms_are_transitive, triangulation_relabelling_from) -def face_edge_perms_init(data): +def face_edge_perms_init(faces): r""" + INPUT:: ``faces`` - a list or a string encoding a permutation + EXAMPLES:: sage: from veerer.triangulation import face_edge_perms_init # random output due to deprecation warnings from realalg @@ -47,6 +50,9 @@ def face_edge_perms_init(data): sage: face_edge_perms_init('(0,1,2)') (array('i', [1, 2, 0]), array('i', [0, 1, 2])) + sage: face_edge_perms_init('(0,1,2)(~0)(~1)(~2)') + (array('i', [1, 2, 0, 3, 4, 5]), array('i', [5, 4, 3, 2, 1, 0])) + TESTS: Check that edge permutation do not depend on the details of faces:: @@ -56,10 +62,10 @@ def face_edge_perms_init(data): sage: f3 = "(6,4,3)(~6,~5,1)(5,0,2)" sage: assert face_edge_perms_init(f1)[1] == face_edge_perms_init(f2)[1] == face_edge_perms_init(f3)[1] """ - if isinstance(data, str): - l = str_to_cycles(data) + if isinstance(faces, str): + l = str_to_cycles(faces) else: - l = [[int(i) for i in c] for c in data] + l = [[int(i) for i in c] for c in faces] pos = [] neg = [] @@ -117,6 +123,37 @@ def face_edge_perms_init(data): return perm_init(fp), perm_init(ep) + +def boundary_init(fp, ep, boundary): + n = len(ep) + + if boundary is None: + return array('i', [0] * n) + elif isinstance(boundary, (array, tuple, list)): + if len(boundary) != n: + raise ValueError('invalid input argument') + return array('i', boundary) + elif isinstance(boundary, dict): + output = array('i', [0] * n) + for e, v in boundary.items(): + if isinstance(e, str): + if not e: + raise ValueError('keys must be valid edges, got {!r}'.format(e)) + elif e[0] == '~': + e = ~int(e[1:]) + else: + e = int(e) + elif isinstance(e, numbers.Integral): + e = int(e) + if e > n: + raise ValueError('keys must be valid edges, got {!r}'.format(e)) + output[e] = v + return output + else: + raise TypeError('invalid boundary data') + + + # NOTE: we don't really care that we have a triangulation here. When # there is no restriction on the cycle decomposition of _fp we got # a coding for any cell decomposition (with possibly folded edges). @@ -142,6 +179,17 @@ def face_edge_perms_init(data): # https://github.com/flatsurf/sage-flatsurf +# TODO: surface_dynamics compatible datastructure +# { +# _n: number of standard edges +# _m: number of self-glued edges +# _vp, _fp: array of size 2 _n + _m +# _boundary: bitset +# } +# boundary requirement: boundary is necessarily glued to a non-boundary +# edges are 0, 1, ..., n+m-1 and half-edges are +# 0, ~0, 1, ~1, ..., (n-1), ~(n-1), n, n+1, ..., n+m-1 + class Triangulation(object): r""" A triangulation of an oriented surface @@ -215,15 +263,40 @@ class Triangulation(object): Traceback (most recent call last): ... ValueError: (fp, ep, vp) do not generate a transitive group + + Examples with boundaries:: + + sage: Triangulation("(0,1,2)(~0)(~1)(~2)", boundary={"~0": 1, "~1": 1, "~2": 1}) + Triangulation("(0,1,2)", boundary="(~2:1)(~1:1)(~0:1)") + sage: Triangulation("(0,1,2)", boundary="(~0:1)(~1:1,~2:1)") + Triangulation("(0,1,2)", boundary="(~2:1,~1:1)(~0:1)") + sage: Triangulation("(0,1,2)(~0,~1,~2)", boundary={"~0": 0}) + Triangulation("(0,1,2)(~2,~0,~1)") + + Examples with invalid boundaries:: + + sage: Triangulation("(0,1,2)(~0,~1,~2)", boundary={"~0": 1, "~1": 1, "~2": 0}) + Traceback (most recent call last): + ... + ValueError: invalid boundary data: 0 at half-edge i=~2 and 1 at half-edge fp[i]=~0 """ - __slots__ = ['_mutable', '_n', '_fp', '_ep', '_vp'] + __slots__ = ['_mutable', '_n', '_fp', '_ep', '_vp', '_bdry'] - def __init__(self, triangles, mutable=False, check=True): + def __init__(self, triangles, boundary=None, mutable=False, check=True): if isinstance(triangles, Triangulation): self._fp = triangles.face_permutation(copy=True) self._ep = triangles.edge_permutation(copy=True) + self._bdry = triangles.boundary_vector(copy=True) else: + if boundary is not None and isinstance(boundary, str): + boundary_cycles, boundary = str_to_cycles_and_data(boundary) + if isinstance(triangles, str): + triangles = str_to_cycles(triangles) + else: + triangles = list(triangles) + triangles.extend(boundary_cycles) self._fp, self._ep = face_edge_perms_init(triangles) + self._bdry = boundary_init(self._fp, self._ep, boundary) self._mutable = mutable @@ -231,24 +304,10 @@ def __init__(self, triangles, mutable=False, check=True): ep = self._ep n = self._n = len(fp) - # TODO: face labels are disabled for now. We should actually - # make it a *partition* of the faces (two faces are allowed to - # have the same label) The trivial partition {0,1,2,...,F-1} - # would correspond to no labels and the atomic - # {{0},{1},...,{F-1}} would correspond to everybody being - # labelled - # fl = self._fl = [None] * F # face labels - - # TODO: edge labels are disabled for now. - # el = self._fl = [None] * E # edge labels - vp = self._vp = array('i', [-1] * n) for i in range(n): vp[fp[ep[i]]] = i - # TODO: vertex labels are disabled for now - # vl = self._vl = perm_cycles(vp)[1].... - if check: self._check(ValueError) @@ -259,10 +318,11 @@ def __getstate__(self): sage: from veerer import Triangulation sage: t = Triangulation("(0,1,2)") sage: dumps(t) # indirect doctest - b'x\x9ck`J.KM-J-\xd2+)\xcaL\xccK/\xcdI,\xc9\xcc\xcf\xe3\nA\xe1\x152h6\x162\xc6\x162ix3z3y3\x00!\x90\xeeLM\xd2\x03\x00\xe7\x1e\x14^' + b'x\x9ck`J.KM-J-\xd2+)\xcaL\xccK/\xcdI,\xc9\xcc\xcf\xe3\nA\xe1\x152h6\x162\xc6\x162ix3z3y3\x00!\x8cf\xe8LM\xd2\x03\x00_u\x15?' """ a = list(self._fp) a.extend(self._ep) + a.extend(self._bdry) a.append(self._mutable) return a @@ -283,10 +343,11 @@ def __setstate__(self, arg): True sage: s1._check() """ - n = (len(arg) - 1) // 2 + n = (len(arg) - 1) // 3 self._n = n self._fp = array('i', arg[:n]) self._ep = array('i', arg[n:2*n]) + self._bdry = array('i', arg[2*n:3*n]) self._mutable = arg[-1] self._vp = array('i', [-1] * n) for i in range(n): @@ -356,7 +417,7 @@ def __hash__(self): return x @staticmethod - def from_face_edge_perms(fp, ep, vp=None, mutable=False, check=True): + def from_face_edge_perms(fp, ep, vp=None, boundary=None, mutable=False, check=True): r""" INPUT: @@ -387,6 +448,10 @@ def from_face_edge_perms(fp, ep, vp=None, mutable=False, check=True): for i in range(n): vp[fp[ep[i]]] = i T._vp = vp + if boundary is None: + T._bdry = array('i', [0] * n) + else: + T._bdry = array('i', boundary) T._mutable = mutable @@ -431,7 +496,7 @@ def _check(self, error=RuntimeError): sage: Triangulation("(0)") Traceback (most recent call last): ... - ValueError: broken face permutation + ValueError: non-triangle internal face starting at half-edge i=0 sage: from array import array sage: fp = array('i', [1,2,0]) @@ -440,7 +505,7 @@ def _check(self, error=RuntimeError): sage: Triangulation.from_face_edge_perms(fp, ep, vp) Traceback (most recent call last): ... - ValueError: fev relation not satisfied + ValueError: fev relation not satisfied at half-edge i=0 sage: fp = array('i', [1,2,0]) sage: ep = array('i', [1,2,0]) @@ -448,10 +513,12 @@ def _check(self, error=RuntimeError): sage: Triangulation.from_face_edge_perms(fp, ep, vp) Traceback (most recent call last): ... - ValueError: broken edge permutation + ValueError: invalid edge permutation at half-edge i=0 """ n = self._n + if not (hasattr(self, '_vp') and hasattr(self, '_ep') and hasattr(self, '_fp') and hasattr(self, '_bdry')): + raise error('missing attributes: these must be _vp, _ep, _fp, _bdry') if not perm_check(self._fp, n): raise error('fp is not a permutation: {}'.format(self._fp)) if not perm_check(self._ep, n): @@ -460,6 +527,8 @@ def _check(self, error=RuntimeError): raise error('vp is not a permutation: {}'.format(self._vp)) if not perms_are_transitive([self._fp, self._ep, self._vp]): raise error('(fp, ep, vp) do not generate a transitive group') + if not isinstance(self._bdry, array) or self._bdry.typecode != 'i' or len(self._bdry) != n: + raise error('bdry must be an integer array of same length as the underlying permutations') # The face, edge, and vertex permutations fp, ep, and vp must # satisfy the relations fp^3 = ep^2 = fp.ep.vp = Id. Note @@ -467,14 +536,19 @@ def _check(self, error=RuntimeError): # \infty) into the symmetric group on the half edges. for i in range(n): - # NOTE: this first test is relevant only if we enforce triangulations - # when we generalize to CW complex we can simply remove it - if self._fp[i] == i or self._fp[self._fp[self._fp[i]]] != i: - raise error('broken face permutation') + j = self._fp[i] + if (self._bdry[i] == 0) != (self._bdry[j] == 0): + raise error('invalid boundary data: {} at half-edge i={} and {} at half-edge fp[i]={}'.format( + self._bdry[i], self._half_edge_string(i), + self._bdry[j], self._half_edge_string(j))) + + for i in range(n): + if not self._bdry[i] and (self._fp[i] == i or self._fp[self._fp[self._fp[i]]] != i): + raise error('non-triangle internal face starting at half-edge i={}'.format(self._half_edge_string(i))) if self._ep[self._ep[i]] != i: - raise error('broken edge permutation') + raise error('invalid edge permutation at half-edge i={}'.format(self._half_edge_string(i))) if self._fp[self._ep[self._vp[i]]] != i: - raise error('fev relation not satisfied') + raise error('fev relation not satisfied at half-edge i={}'.format(self._half_edge_string(i))) def _check_half_edge(self, e): if not isinstance(e, numbers.Integral): @@ -561,6 +635,7 @@ def copy(self, mutable=None): T._fp = self._fp T._ep = self._ep T._vp = self._vp + T._bdry = self._bdry T._mutable = mutable return T @@ -570,6 +645,7 @@ def copy(self, mutable=None): T._fp = self._fp[:] T._ep = self._ep[:] T._vp = self._vp[:] + T._bdry = self._bdry[:] T._mutable = mutable return T @@ -666,6 +742,12 @@ def vertex_permutation(self, copy=True): else: return self._vp + def boundary_vector(self, copy=True): + if copy: + return self._bdry[:] + else: + return self._bdry + def next_at_vertex(self, i, check=True): r""" EXAMPLES:: @@ -720,7 +802,7 @@ def previous_at_vertex(self, i, check=True): i = self._check_half_edge(i) return self._fp[self._ep[i]] - def previous_in_face(self, i): + def previous_in_face(self, i, check=True): r""" EXAMPLES:: @@ -734,6 +816,8 @@ def previous_in_face(self, i): sage: T.previous_in_face(3) 5 """ + if check: + i = self._check_half_edge(i) return self._ep[self._vp[i]] def num_half_edges(self): @@ -763,22 +847,52 @@ def _norm(self, e): f = self._ep[e] return f if f < e else e + def _half_edge_string(self, e): + f = self._ep[e] + return '~%d' % f if f < e else '%d' % e + def num_faces(self): - return perm_num_cycles(self._fp) + r""" + Return the number of triangles. + """ + return sum(self._bdry[c[0]] == 0 for c in perm_cycles(self._fp)) + + def num_boundaries(self): + r""" + Return the number of boundaries. + """ + return sum(self._bdry[c[0]] == 1 for c in perm_cycles(self._fp)) def faces(self): r""" - Return the list of faces as triples of half-edges + Return the list of internal faces as triples of half-edges EXAMPLES:: - sage: from veerer import * + sage: from veerer import Triangulation sage: T = Triangulation("(0,1,2)(3,4,5)(~0,~3,6)") sage: T.faces() [[0, 1, 2], [3, 4, 5], [6, 8, 7]] """ - return perm_cycles(self._fp, True, self._n) + return [c for c in perm_cycles(self._fp, True, self._n) if self._bdry[c[0]] == 0] + + def boundaries(self): + r""" + Return the list of boundaries as lists of half-edges. + + EXAMPLES:: + + sage: from veerer import Triangulation + + sage: T = Triangulation("(0,1,2)(3,4,5)(~0,~3,6)") + sage: T.boundaries() + [] + sage: T = Triangulation("(0,1,2)(~0)(~1,~2)", {"~0": 1, "~1": 1, "~2": 1}) + sage: T.boundaries() + [[3, 4], [5]] + """ + return [c for c in perm_cycles(self._fp, True, self._n) if self._bdry[c[0]]] def edges(self): r""" @@ -809,6 +923,9 @@ def vertices(self): return perm_cycles(self._vp, True, self._n) def num_vertices(self): + r""" + Return the number of vertices. + """ return perm_num_cycles(self._vp) def euler_characteristic(self): @@ -836,6 +953,18 @@ def euler_characteristic(self): sage: T = Triangulation("(0,1,2)(~2,3,4)(~4,5,6)(~6,~0,7)(~7,~1,8)(~8,~3,~5)") sage: T.euler_characteristic() -2 + + A cylinder:: + + sage: T = Triangulation("(0,1,2)(~0,3,4)(~1,~2)(~3,~4)", {"~1": 1, "~2": 1, "~3": 1, "~4": 1}) + sage: T.euler_characteristic() + 0 + + A pair of pants:: + + sage: T = Triangulation("(0,1,2)(~0)(~1)(~2)", {"~0": 1, "~1": 1, "~2": 1}) + sage: T.euler_characteristic() + -1 """ return self.num_faces() - self.num_edges() + (self.num_vertices() + self.num_folded_edges()) @@ -854,9 +983,21 @@ def genus(self): sage: T = Triangulation([(-3,1,-1), (-2,0,2)]) sage: T.genus() 1 + + A cylinder:: + + sage: T = Triangulation("(0,1,2)(~0,3,4)(~1,~2)(~3,~4)", {"~1": 1, "~2": 1, "~3": 1, "~4": 1}) + sage: T.genus() + 0 + + A pair of pants:: + + sage: T = Triangulation("(0,1,2)(~0)(~1)(~2)", {"~0": 1, "~1": 1, "~2": 1}) + sage: T.genus() + 0 """ - # chi = 2 - 2g - return (2 - self.euler_characteristic()) // 2 + # chi = 2 - 2g - n + return (2 - self.euler_characteristic() - self.num_boundaries()) // 2 def __str__(self): r""" @@ -867,7 +1008,13 @@ def __str__(self): sage: str(T) 'Triangulation("(0,1,2)(~2,~0,~1)")' """ - return 'Triangulation("%s")' % perm_cycle_string(self._fp, n=self._n, involution=self._ep) + cycles = perm_cycles(self._fp, n=self._n) + face_cycles = perm_cycles_to_string([c for c in cycles if not self._bdry[c[0]]], involution=self._ep) + bdry_cycles = perm_cycles_to_string([c for c in cycles if self._bdry[c[0]]], involution=self._ep, data=self._bdry) + if bdry_cycles: + return 'Triangulation("%s", boundary="%s")' % (face_cycles, bdry_cycles) + else: + return 'Triangulation("%s")' % face_cycles def __repr__(self): return str(self) @@ -1176,7 +1323,7 @@ def is_flippable(self, e): E = self._ep[e] a = self._fp[e] b = self._fp[a] - return a != E and b != E + return a != E and b != E and self._bdry[e] == 0 and self._bdry[E] == 0 def flippable_edges(self): r""" @@ -1190,12 +1337,16 @@ def flippable_edges(self): sage: V = VeeringTriangulation(T, [RED, RED, BLUE]) sage: V.flippable_edges() [0, 1] + + sage: T = Triangulation("(0,1,2)(~0,3,4)(~1,~2)(~3,~4)", {"~1": 1, "~2": 1, "~3": 1, "~4": 1}) + sage: T.flippable_edges() + [0] """ n = self._n ep = self._ep return [e for e in range(n) if e <= ep[e] and self.is_flippable(e)] - def square_about_edge(self, e): + def square_about_edge(self, e, check=True): r""" Return the four edges that makes ``e`` the diagonal of a quadrilateral. @@ -1225,8 +1376,12 @@ def square_about_edge(self, e): # v/ c | # x---------->x - e = int(e) + if check: + e = self._check_half_edge(e) + E = self._ep[e] + if self._bdry[e] or self._bdry[E]: + raise ValueError('non internal edge') a = self._fp[e] b = self._fp[a] @@ -1235,7 +1390,7 @@ def square_about_edge(self, e): return a, b, c, d - def swap(self, e): + def swap(self, e, check=True): r""" Change the orientation of the edge ``e``. @@ -1288,6 +1443,9 @@ def swap(self, e): if not self._mutable: raise ValueError('immutable triangulation; use a mutable copy instead') + if check: + e = self._check_half_edge(e) + vp = self._vp ep = self._ep fp = self._fp @@ -1321,6 +1479,8 @@ def swap(self, e): vp[e] = E_vp vp[E] = e_vp + self._bdry[e], self._bdry[E] = self._bdry[E], self._bdry[e] + def relabel(self, p, check=True): r""" Relabel this triangulation inplace according to the permutation ``p``. @@ -1354,6 +1514,13 @@ def relabel(self, p, check=True): sage: T.relabel("(0,2)(1,3)") sage: T == T0 True + + An example with boundary:: + + sage: t = Triangulation("(0,1,2)(~0,3,4)(~4,~3,~2,~1)", {"~1": 1, "~2": 1, "~3": 1, "~4": 1}, mutable=True) + sage: t.relabel("(0,3)(1,~2)") + sage: t + Triangulation("(0,4,~3)(3,~2,~1)", boundary="(1:1,2:1,~4:1,~0:1)") """ if not self._mutable: raise ValueError('immutable triangulation; use a mutable copy instead') @@ -1369,6 +1536,7 @@ def relabel(self, p, check=True): self._vp = perm_conjugate(self._vp, p) self._ep = perm_conjugate(self._ep, p) self._fp = perm_conjugate(self._fp, p) + perm_on_list(p, self._bdry, n) def flip(self, e, check=True): r""" @@ -1392,6 +1560,17 @@ def flip(self, e, check=True): sage: T == Triangulation("(0,1,2)") True + + An example with boundaries:: + + sage: t = Triangulation("(0,1,2)(~0,3,4)", boundary="(~4:1,~3:1,~2:1,~1:1)", mutable=True) + sage: t.flip(0) + sage: t + Triangulation("(0,2,3)(1,~0,4)", boundary="(~4:1,~3:1,~2:1,~1:1)") + sage: t.flip(2) + Traceback (most recent call last): + ... + ValueError: can not flip non internal edge 2 """ # v<----------u v<----------u # | a ^^ |^ a ^ @@ -1414,6 +1593,8 @@ def flip(self, e, check=True): e = self._check_half_edge(e) E = self._ep[e] + if self._bdry[e] or self._bdry[E]: + raise ValueError('can not flip non internal edge %s' % self._norm(e)) a = self._fp[e] b = self._fp[a] @@ -1508,6 +1689,9 @@ def flip_back(self, e, check=True): E = self._ep[e] + if self._bdry[e] or self._bdry[E]: + raise ValueError('can not flip non internal edge %s' % self._norm(e)) + a = self._fp[e] b = self._fp[a] if a == E or b == E: @@ -1677,6 +1861,7 @@ def _relabelling_from(self, start_edge): return triangulation_relabelling_from(self._vp, self._ep, start_edge) def _automorphism_good_starts(self): + # TODO: we should also discriminate based on boundaries! # we discriminate based on the lengths of cycles in # vp, ep, fp and vp[ep[fp]] n = self._n @@ -1689,7 +1874,7 @@ def _automorphism_good_starts(self): d = {} for i in range(n): - s = ((lv[i] * n + le[i]) * n + lf[i]) * n + lvef[i] + s = (((self._bdry[i] * n + lv[i]) * n + le[i]) * n + lf[i]) * n + lvef[i] if s in d: d[s].append(i) else: @@ -1716,6 +1901,9 @@ def automorphisms(self): The output is a list of arrays that are permutations acting on the set of half edges. + For triangulations with boundaries, we allow automorphism to permute + boundaries. Though, boundary edge have to be mapped on boundary edge. + EXAMPLES:: sage: from veerer import * @@ -1745,9 +1933,22 @@ def automorphisms(self): sage: for a in A: ....: S.relabel(a) ....: assert S == V + + Examples with boundaries:: + + sage: t = Triangulation("(0,1,2)", boundary="(~0:1)(~1:1)(~2:1)") + sage: len(t.automorphisms()) + 3 + sage: t = Triangulation("(0,1,2)", boundary="(~0:1,~1:1,~2:1)") + sage: len(t.automorphisms()) + 3 + sage: t = Triangulation("(0,1,2)", boundary="(~0:1,~1:1,~2:2)") + sage: len(t.automorphisms()) + 1 """ fp = self._fp ep = self._ep + bdry = self._bdry best = None best_relabellings = [] @@ -1757,8 +1958,9 @@ def automorphisms(self): fp_new = perm_conjugate(fp, relabelling) ep_new = perm_conjugate(ep, relabelling) + bdry_new = perm_on_list(relabelling, self._bdry, self._n) - T = (fp_new, ep_new) + T = (fp_new, ep_new, bdry_new) if best is None or T == best: best_relabellings.append(relabelling) best = T @@ -1785,8 +1987,10 @@ def best_relabelling(self, all=False): fp_new = perm_conjugate(fp, relabelling) ep_new = perm_conjugate(ep, relabelling) + bdry_new = self._bdry[:] + perm_on_list(relabelling, bdry_new, self._n) - T = (fp_new, ep_new) + T = (fp_new, ep_new, bdry_new) if best is None or T < best: best_relabelling = relabelling best = T @@ -1841,6 +2045,20 @@ def set_canonical_labels(self): ....: L.relabel(p) ....: L.set_canonical_labels() ....: assert L == L0, (L, L0) + + Examples with boundaries:: + + sage: t = Triangulation("(0,1,2)(~0,3,4)", boundary="(~1:1)(~2:2)(~3:3)(~4:4)", mutable=True) + sage: t.set_canonical_labels() + sage: t + Triangulation("(1,~3,~2)(~4,~1,~0)", boundary="(0:1)(2:4)(3:3)(4:2)") + + sage: t0 = t.copy(mutable=False) + sage: for _ in range(10): + ....: p = perm_random_centralizer(t.edge_permutation(copy=False)) + ....: t.relabel(p) + ....: t.set_canonical_labels() + ....: assert t == t0, (t, t0) """ if not self._mutable: raise ValueError('immutable triangulation; use a mutable copy instead') diff --git a/veerer/veering_triangulation.py b/veerer/veering_triangulation.py index 694267b..6dd6474 100644 --- a/veerer/veering_triangulation.py +++ b/veerer/veering_triangulation.py @@ -35,13 +35,15 @@ import ppl from sage.structure.richcmp import op_LT, op_LE, op_EQ, op_NE, op_GT, op_GE, rich_to_bool +from sage.modules.free_module import FreeModule +from sage.matrix.constructor import matrix from .constants import * from .permutation import * from .misc import det2 from .triangulation import Triangulation from .polyhedron import LinearExpressions, ConstraintSystem -from .polyhedron.linear_algebra import linear_form_project, linear_form_normalize +from .polyhedron.linear_algebra import linear_form_project, vector_normalize class VeeringTriangulation(Triangulation): @@ -88,11 +90,34 @@ class VeeringTriangulation(Triangulation): sage: VeeringTriangulation("(0,~6,~3)(1,7,~2)(2,~1,~0)(3,5,~4)(4,8,~5)(6,~8,~7)", "RBPBRBPRB") VeeringTriangulation("(0,~6,~3)(1,7,~2)(2,~1,~0)(3,5,~4)(4,8,~5)(6,~8,~7)", "RBPBRBPRB") + + Triangulations with boundary:: + + sage: VeeringTriangulation("(0,1,2)(~0,3,4)", "(~1:1)(~2:1)(~3:1)(~4:1)", "RBRBR") + VeeringTriangulation("(0,1,2)(3,4,~0)", boundary="(~4)(~3)(~2)(~1)", colouring="RBRBR") + sage: VeeringTriangulation("(0,1,2)(~0,3,4)", boundary="(~1:1)(~2:1)(~3:1)(~4:1)", colouring="RBRBR") + VeeringTriangulation("(0,1,2)(3,4,~0)", boundary="(~4)(~3)(~2)(~1)", colouring="RBRBR") """ __slots__ = ['_colouring'] - def __init__(self, triangulation, colouring, mutable=False, check=True): - Triangulation.__init__(self, triangulation, mutable=mutable, check=False) + def __init__(self, *args, triangulation=None, boundary=None, colouring=None, mutable=False, check=True): + if len(args) == 3: + if triangulation is not None or boundary is not None or colouring is not None: + raise ValueError('invalid data in constructor') + # triangulation_data boundary colouring + triangulation, boundary, colouring = args + elif len(args) == 2: + # either (triangulation, colouring) + # or (triangulation_data, colouring) + if triangulation is not None or colouring is not None: + raise ValueError('invalid data in constructor') + triangulation, colouring = args + elif len(args) == 1: + if triangulation is not None: + raise ValueError('invalid data in constructor') + triangulation, = args + + Triangulation.__init__(self, triangulation, boundary=boundary, mutable=mutable, check=False) n = self._n # number of half edges (get initialized by the triangulation) if isinstance(colouring, str): @@ -137,6 +162,8 @@ def _check(self, error=RuntimeError): # 1-degenerate: PBR (PURPLE), GRB (GREEN) # 2-degenerate: BPG (BLUE|PURPLE|GREEN), RGP (RED|GREEN|PURPLE) for a in range(n): + if self._bdry[a]: + continue col, a, b, c = self.triangle(a) good = False if col == BLUE: @@ -157,9 +184,11 @@ def _check(self, error=RuntimeError): # no monochromatic vertex for v in self.vertices(): + if self._bdry[v[0]]: + continue col = cols[v[0]] i = 1 - while i < len(v) and cols[v[i]] == col: + while i < len(v) and cols[v[i]] == col and not self._bdry[v[i]]: i += 1 if i == len(v): raise error('monochromatic vertex {} of colour {}'.format(v, colour_to_string(cols[v[0]]))) @@ -171,10 +200,15 @@ def __getstate__(self): sage: from veerer import VeeringTriangulation sage: t = VeeringTriangulation("(0,1,2)", "BBR") sage: dumps(t) # indirect doctest - b'x\x9ck`J.KM-J-\xd2\x03Q\x99y\xe9\xf1%E\x99\x89y\xe9\xa59\x89%\x99\xf9y\\a\x10\xd1\x10\x14\xc1B\x06\xcd\xc6B\xc6\xd8B&\roFo&o\x06 \x04\xd1 \xc8\xd8\x99\x9a\xa4\x07\x005\xe2\x1bc' + b'x\x9ck`J.KM-J-\xd2\x03Q\x99y\xe9\xf1%E\x99\x89y\xe9\xa59\x89%\x99\xf9y\\a\x10\xd1\x10\x14\xc1B\x06\xcd\xc6B\xc6\xd8B&\roFo&o\x06 \x84\xd1\x0c@\x9a\xc9\x9b\xb135I\x0f\x00\xd8*\x1cD' + sage: t = VeeringTriangulation("(0,1,2)(~0,3,4)", "(~1:1)(~2:1)(~3:1)(~4:1)", "RBRBR") + sage: dumps(t) # indirect doctest + b'x\x9ck`J.KM-J-\xd2\x03Q\x99y\xe9\xf1%E\x99\x89y\xe9\xa59\x89%\x99\xf9y\\a\x10\xd1\x10\x14\xc1B\x06\xcd\xc6B\xc6\xd8B&\roFo&o\x06o\x16oNoVo6ovo\x0eof \x9b\x03\xc8b\x03\x8a\xb0\x00yL@5\x0cH\x90\x11\n\x19\xc0z!\x18\xce\xeaLM\xd2\x03\x00\xd7\x88$\xd9' """ + # TODO: add boundary if needed a = list(self._fp) a.extend(self._ep) + a.extend(self._bdry) a.extend(self._colouring) a.append(self._mutable) return a @@ -186,6 +220,8 @@ def __setstate__(self, arg): sage: from veerer import VeeringTriangulation sage: t0 = VeeringTriangulation("(0,1,2)", "BBR", mutable=False) sage: t1 = VeeringTriangulation("(0,1,2)", "BBR", mutable=True) + sage: t2 = VeeringTriangulation("(0,1,2)(~0,3,4)", "(~1:1)(~2:1)(~3:1)(~4:1)", "RBRBR") + sage: s0 = loads(dumps(t0)) # indirect doctest sage: assert s0 == t0 sage: s0._mutable @@ -197,12 +233,16 @@ def __setstate__(self, arg): sage: s1._mutable True sage: s1._check() + + sage: s2 = loads(dumps(t2)) + sage: assert s2 == t2 """ - n = (len(arg) - 1) // 3 + n = (len(arg) - 1) // 4 self._n = n self._fp = array('i', arg[:n]) self._ep = array('i', arg[n:2*n]) - self._colouring = array('i', arg[2*n:3*n]) + self._bdry = array('i', arg[2*n:3*n]) + self._colouring = array('i', arg[3*n:4*n]) self._mutable = arg[-1] self._vp = array('i', [-1] * n) for i in range(n): @@ -311,7 +351,7 @@ def as_linear_family(self, mutable=False): P = self.train_track_linear_space() return VeeringTriangulationLinearFamily(self, self.train_track_linear_space().lines(), mutable=mutable) - def triangle(self, a): + def triangle(self, a, check=True): r""" Return a quadruple ``(colour, e0, e1, e2)`` in canonical form for the triangle with half edge ``a``. @@ -345,6 +385,10 @@ def triangle(self, a): # non-dgenerate: BBR (BLUE), RRB (RED) # 1-degenerate: PBR (PURPLE), GRB (GREEN) # 2-degenerate: RGP (RED|GREEN|PURPLE), BPG (BLUE|GREEN|PURPLE) + if check: + a = self._check_half_edge(a) + if self._bdry[a]: + raise ValueError('boundary edge') b = self._fp[a] c = self._fp[b] cols = self._colouring @@ -604,6 +648,8 @@ def from_face_edge_perms(self, colouring, fp, ep, vp=None, mutable=False, check= vp = array('i', [-1] * n) for i in range(n): vp[fp[ep[i]]] = i + + T._bdry = array('i', [0] * n) T._vp = vp T._colouring = colouring T._mutable = mutable @@ -777,6 +823,7 @@ def copy(self, mutable=None): T._vp = self._vp T._ep = self._ep T._fp = self._fp + T._bdry = self._bdry T._colouring = self._colouring T._mutable = mutable @@ -787,6 +834,7 @@ def copy(self, mutable=None): T._vp = self._vp[:] T._ep = self._ep[:] T._fp = self._fp[:] + T._bdry = self._bdry[:] T._colouring = self._colouring[:] T._mutable = mutable @@ -821,10 +869,22 @@ def __str__(self): sage: VeeringTriangulation("(0,1,2)", [RED, RED, BLUE]) VeeringTriangulation("(0,1,2)", "RRB") + sage: VeeringTriangulation("(0,1,2)(3,4,~0)", "(~4:1,~3:1,~2:1,~1:1)", "RRBRB") + VeeringTriangulation("(0,1,2)(3,4,~0)", boundary="(~4,~3,~2,~1)", colouring="RRBRB") """ - return "VeeringTriangulation(\"{}\", \"{}\")".format( - perm_cycle_string(self._fp, False, self._n, self._ep), - self._colouring_string(short=True)) + cycles = perm_cycles(self._fp, n=self._n) + face_cycles = perm_cycles_to_string([c for c in cycles if self._bdry[c[0]] == 0], involution=self._ep) + bdry_cycles = perm_cycles_to_string([c for c in cycles if self._bdry[c[0]] == 1], involution=self._ep) + colouring = self._colouring_string(short=True) + if bdry_cycles: + return "VeeringTriangulation(\"{}\", boundary=\"{}\", colouring=\"{}\")".format( + face_cycles, + bdry_cycles, + colouring) + else: + return "VeeringTriangulation(\"{}\", \"{}\")".format( + face_cycles, + colouring) def __repr__(self): return str(self) @@ -834,10 +894,84 @@ def colour(self, e): e = int(e) return self._colouring[e] - def angles(self): + def vertex_angle(self, e): + r""" + Return the angle at the vertex the half-edge ``e`` is adjacent to. + + EXAMPLES:: + + sage: from veerer import VeeringTriangulation + sage: T = VeeringTriangulation([(-12, 4, -4), (-11, -1, 11), (-10, 0, 10), + ....: (-9, 9, 1), (-8, 8, -2), (-7, 7, 2), + ....: (-6, 6, -3), (-5, 5, 3)], + ....: [RED, RED, RED, RED, BLUE, BLUE, BLUE, BLUE, BLUE, BLUE, BLUE, BLUE]) + sage: T.vertex_angle(0) + 1 + sage: T.vertex_angle(1) + 3 + """ + e = self._check_half_edge(e) + vp = self.vertex_permutation(copy=False) + + a = 0 + e0 = e + col = self._colouring[e] + # NOTE: we count by multiples of pi/2 and then divide by 2 + ee = vp[e] + while True: + ee = vp[e] + ccol = self._colouring[ee] + switch = bool(col != ccol and (((col & (BLUE|RED)) and (ccol & (BLUE|RED))) or (col & (GREEN|PURPLE)))) + a += switch + 2 * self._bdry[e] + e = ee + col = ccol + if e == e0: + break + assert a % 2 == 0, a + return a // 2 + + def face_angle(self, e): + r""" + Return the angle associated to the pole in the middle of the face ``e`` is adjacent to. + + EXAMPLES:: + + sage: from veerer import Triangulation, VeeringTriangulation + sage: t = Triangulation("(0,1,2)", boundary="(~2:2,~1:1,~0:1)") + sage: vt = VeeringTriangulation(t, "BRB") + sage: vt.face_angle(3) + -2 + """ + e = self._check_half_edge(e) + fp = self.face_permutation(copy=False) + + if not self._bdry[e]: + raise ValueError('e={} not on a boundary face'.format(self._half_edge_string(e))) + + cum_angle = 0 + alternations = 0 + num_sides = 0 + col = self._colouring[e] + e0 = e + while True: + ee = fp[e] + ccol = self._colouring[ee] + cum_angle += self._bdry[e] + alternations += (col != ccol) + num_sides += 1 + col = ccol + e = ee + if e == e0: + break + return (num_sides - cum_angle - alternations // 2) + + def angles(self, half_edge_representatives=False): r""" Return the list of angles (divided by \pi). + There is a positive angle for each vertex and a non-positive angle for + each boundary face (if any). + EXAMPLES:: sage: from veerer import * @@ -871,30 +1005,88 @@ def angles(self): sage: t.forgot_forward_flippable_colour() sage: t.angles() [3, 3, 3, 3] + + The plane with 4 marked points:: + + sage: t = Triangulation("(0,1,2)(~0,3,4)", boundary="(~4:1,~3:1,~2:1,~1:1)") + sage: VeeringTriangulation(t, "RRBRB").angles() + [2, 2, 2, 2, -2] + + A bi-infinite cylinder with 2 marked points:: + + sage: t = Triangulation("(0,1,2)(~1,~2,3)", boundary="(~0:1)(~3:1)") + sage: VeeringTriangulation(t, "RBRR").angles() + [2, 2, 0, 0] + + A bi-infinite cylinder with a single marked points:: + + sage: t = Triangulation("", boundary="(0:1)(~0:1)") + sage: VeeringTriangulation(t, "R").angles() + [2, 0, 0] + + Complement of a segment in the plane:: + + sage: t = Triangulation("", boundary="(0:2,~0:2)") + sage: VeeringTriangulation(t, "R").angles() + [2, 2, -2] + + Complement of a triangle in the plane:: + + sage: t = Triangulation("(0,1,2)", boundary="(~2:2,~1:1,~0:1)") + sage: VeeringTriangulation(t, "BRB").angles() + [2, 2, 2, -2] """ + # TODO: handle non-connected stratum correctly! n = self.num_half_edges() angles = [] + half_edges = [] seen = [False] * n vp = self.vertex_permutation(copy=False) - for e in range(n): - if seen[e]: continue + fp = self.face_permutation(copy=False) + # zeros (and simple poles) from vertices + for e in range(n): + if seen[e]: + continue a = 0 col = self._colouring[e] # NOTE: we count by multiples of pi/2 and then divide by 2 + half_edges.append(e) while not seen[e]: seen[e] = True ee = vp[e] ccol = self._colouring[ee] - a += bool(col != ccol and (((col & (BLUE|RED)) and (ccol & (BLUE|RED))) or - (col & (GREEN|PURPLE)))) + switch = bool(col != ccol and (((col & (BLUE|RED)) and (ccol & (BLUE|RED))) or (col & (GREEN|PURPLE)))) + a += switch + 2 * self._bdry[e] e = ee col = ccol - a += col != ccol - + assert a % 2 == 0, a angles.append(a // 2) - return angles + [1] * self.num_folded_edges() + # angles at infinity from boundary faces (meromorphic differentials) + seen = [False] * n + for e in range(n): + if seen[e] or not self._bdry[e]: + continue + + cum_angle = 0 + alternations = 0 + num_sides = 0 + col = self._colouring[e] + half_edges.append(e) + while not seen[e]: + seen[e] = True + ee = fp[e] + ccol = self._colouring[ee] + cum_angle += self._bdry[e] + alternations += (col != ccol) + num_sides += 1 + col = ccol + e = ee + angles.append(num_sides - cum_angle - alternations // 2) + angles.extend([1] * self.num_folded_edges()) + + return (angles, half_edges) if half_edge_representatives else angles def is_abelian(self, certificate=False): r""" @@ -949,14 +1141,26 @@ def is_abelian(self, certificate=False): sage: t.forgot_backward_flippable_colour() sage: t.is_abelian() False + + Examples with boundaries:: + + sage: t = Triangulation("", "(0:1)(~0:1)") + sage: VeeringTriangulation(t, "R").is_abelian() + True + + sage: t = Triangulation("", "(0:1)(1:1)(~0:1,2:1,~1:1,~2:1)") + sage: VeeringTriangulation(t, "RRR").is_abelian() + False """ ep = self._ep vp = self._vp cols = self._colouring - # Perform a BFS. We store an orientation for each edge that specifies - # its orientation in R^2: True if x > 0 or (x == 0 and y > 0) - # False if x < 0 or (x == 0 and y < 0) + # Try to choose holonomy with signs for each vector. It is enough + # to store a boolean for each edge: + # True: (+, +) or (+,-) + # False: (-,-) or (-,+) + # There is a coherence around vertices and faces oris = [None] * self._n oris[0] = True q = [0] @@ -965,7 +1169,10 @@ def is_abelian(self, certificate=False): o = oris[e0] f = vp[e] while f != e0: - if (cols[e] == RED and cols[f] == BLUE) or (cols[e] == GREEN): + if self._bdry[f]: + if not (cols[e] == RED and cols[f] == BLUE): + o = not o + elif (cols[e] == RED and cols[f] == BLUE) or (cols[e] == GREEN): o = not o oris[f] = o @@ -976,8 +1183,12 @@ def is_abelian(self, certificate=False): return (False, None) if certificate else False e, f = f, vp[f] - if (cols[e] == RED and cols[f] == BLUE) or cols[e] == GREEN: + if self._bdry[f]: + if not (cols[e] == RED and cols[f] == BLUE): + o = not o + elif (cols[e] == RED and cols[f] == BLUE) or (cols[e] == GREEN): o = not o + if oris[e0] != o: return (False, None) if certificate else False @@ -1092,17 +1303,23 @@ def stratum(self): sage: t.forgot_forward_flippable_colour() # optional - surface_dynamics sage: t.stratum() # optional - surface_dynamics Q_2(1^4) + + Some examples with boundaries:: + + sage: t = Triangulation("", "(0:1)(~0:1)") + sage: VeeringTriangulation(t, "R").stratum() # optional - surface_dynamics + H_0(0, -1^2) """ from .features import surface_dynamics_feature surface_dynamics_feature.require() + from surface_dynamics.flat_surfaces.strata import Stratum + A = self.angles() - if any(a%2 for a in A) or not self.is_abelian(): - from surface_dynamics import QuadraticStratum - return QuadraticStratum([(a-2) for a in A]) + if any(a % 2 for a in A) or not self.is_abelian(): + return Stratum([(a - 2) for a in A], 2) else: - from surface_dynamics import AbelianStratum - return AbelianStratum([(a-2)/2 for a in A]) + return Stratum([(a - 2) // 2 for a in A], 1) def dimension(self): r""" @@ -1565,6 +1782,7 @@ def _automorphism_good_starts(self, has_purple=None, has_green=None): return starts + # TODO: check the argument!!! def edge_colour(self, e): return self._colouring[e] @@ -2879,6 +3097,13 @@ def train_track_polytope(self, slope=VERTICAL, low_bound=0, backend=None): [[1, 0, 1], [1, 1, 0]] sage: sorted(T.train_track_polytope(HORIZONTAL, backend='sage').rays()) [[0, 1, 1]] + + Examples with boundaries (meromorphic differentials):: + + sage: VeeringTriangulation("", "(0:1)(~0:1)", "R").train_track_polytope() + Cone of dimension 1 in ambient dimension 1 made of 1 facets (backend=ppl) + sage: VeeringTriangulation("(0,1,2)(~1,~2,3)", "(~0:1)(~3:1)", "RBRR").train_track_polytope() + Cone of dimension 2 in ambient dimension 4 made of 2 facets (backend=ppl) """ from sage.rings.integer_ring import ZZ L = LinearExpressions(ZZ) @@ -3750,7 +3975,7 @@ def geometric_flips(self, backend=None): constraint = x[self._norm(e)] == y[self._norm(a)] + y[self._norm(d)] constraint = constraint.coefficients(dim=2*ne, homogeneous=True) linear_form_project(eqns, constraint) - linear_form_normalize(base_ring, constraint) + vector_normalize(base_ring, constraint) constraint = tuple(constraint) if constraint in delaunay_facets: delaunay_facets[constraint].append(e) @@ -3831,6 +4056,124 @@ def random_forward_flip(self, repeat=1): else: self.flip_back(e, old_col) + def switch_subspace_generators_matrix(self, slope=VERTICAL): + r""" + EXAMPLES:: + + sage: from veerer import VeeringTriangulation + sage: vt = VeeringTriangulation("(0,1,2)(3,4,5)(6,7,8)(~0,~7,~5)(~3,~4,~2)(~6,~1,~8)", "RRBRRBRRB") + sage: vt.switch_subspace_generators_matrix() + [ 1 1 0 0 0 0 1 1 0] + [ 0 0 0 1 1 0 0 0 0] + [ 1 0 -1 1 0 -1 0 0 0] + [ 0 0 0 0 0 0 1 0 -1] + """ + return matrix(self.base_ring(), self.train_track_switch_constraints(slope).cone().lines()) + + def parallel_cylinders(self, col=RED): + r""" + Return the L-parallel cylinders. + + EXAMPLES:: + + sage: from veerer.linear_family import VeeringTriangulationLinearFamilies + sage: X9 = VeeringTriangulationLinearFamilies.prototype_H1_1(0, 2, 1, -1) + sage: X9.parallel_cylinders() + [([(0, 0, 0, 1, 1, 0, 0, 1, 1), (0, 0, 0, 0, 0, 1, 1, 0, 0)], [1, 1])] + """ + cylinders = list(self.cylinders(col)) + if not cylinders: + [] + + B = [] # boundary edges + C = [] # middle edges + for cyl in cylinders: + b = [0] * self.num_edges() # indicatrix of bottom edges + t = [0] * self.num_edges() # indicatrix of top edges + c = [0] * self.num_edges() # indicatrix of the middle edges + for e in cyl[0]: + c[self._norm(e)] = 1 + for e in cyl[1]: + b[self._norm(e)] = 1 + for e in cyl[2]: + t[self._norm(e)] = 1 + C.append(c) + B.append((b, t)) + + # take intersection of the cylinder twists in the tangent space + F = FreeModule(self.base_ring(), self.num_edges()) + U = F.submodule(self.switch_subspace_generators_matrix()) + T = F.submodule(C) + + # what would be even nicer are the equations on the coefficients of C + basis = U.intersection(T).basis_matrix() + twist_coeffs = [] + C = matrix(C) + seen = set() + parallel_families = [] + for b in basis.rows(): + # conjecture: the reduced echelon form is a basis with the following property + # * the nonzero positions for elements of the basis are disjoint (correspond to a weighted union of cylinders) + # * it is >= 0 + # In other words, the equivalence classes of Benirschke are tight by a single equation. + # Moreover, the coefficients give the circumference ratios + u = C.solve_left(b) + u = vector_normalize(self.base_ring(), u) + pos = u.nonzero_positions() + assert not any(i in seen for i in pos) + seen.update(pos) + parallel_cylinders = [C[i] for i in pos] + parallel_families.append((parallel_cylinders, [u[i] for i in pos])) + + return parallel_families + + def contraction(self, edges): + r""" + Return the degeneration obtained by contracting the edges in ``contractions`` and + blowing-up the ones in blowups. + """ + + def blowup(self, edges): + n = self._n + edges = set(edges) + edges.update([self._ep[e] for e in edges]) + + kept_half_edges = [e for e in range(self._n) if e not in edges] + m = len(kept_half_edges) + relabelling = {e: j for j, e in enumerate(kept_half_edges)} + + fp = array('i', [0] * m) + ep = array('i', [0] * m) + boundary = array('i', [0] * m) + colouring = array('i', [0] * m) + for e in half_edges: + ep[relabelling[e]] = relabelling[self._ep[e]] + f = self._fp[e] + while f not in relabelling: + f = self._fp[f] + fp[relabelling[e]] = relabelling[f] + # TODO: should be updated + boundary[relabelling[e]] = self._bdry[e] + colouring[relabelling[e]] = self._colouring[e] + + t = Triangulation.from_face_edge_perms(fp, ep, boundary=boundary) + return VeeringTriangulation(t, colouring), edge_map + + def codimension_one_horizontal_degenerations(self, i=None): + r""" + Return codimension one horizontal degenerations. + """ + if i is not None and i != 0: + raise ValueError('invalid level specification') + + n = self._n + for cyl_family in self.parallel_cylinders(): + edges = [] + for cyl in cyl_family[0]: + for i, j in enumerate(cyl): + if j: + edges.append(i) + yield self.blowup(edges) class VeeringTriangulations: @staticmethod