Source code for flatsurf.geometry.gl2r_orbit_closure

r"""
GL(2,R)-orbit closure of translation surfaces

.. TODO::

    - Theorem 1.9 of Alex Wright: the field of definition is contained in the field generated by
      the ratio of circumferences. We should provide a method, .reset_field_of_definition or
      something similar

    - flow decompositions should be accessible on surfaces rather than on ``GL2ROrbitClosure``

EXAMPLES:

Let us first construct a Veech surface in the stratum H(2)::

    sage: from flatsurf import translation_surfaces
    sage: from flatsurf import GL2ROrbitClosure

    sage: x = polygen(QQ)
    sage: K.<a> = NumberField(x^3 - 2, embedding=AA(2)**(1/3))
    sage: S = translation_surfaces.mcmullen_L(1,1,1,a)
    sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf  # random output due to matplotlib warnings with some combinations of setuptools and matplotlib
    sage: O.decomposition((1,2)).cylinders() # optional: pyflatsurf
    [Cylinder with perimeter [...]]

The following is also a Veech surface. However the flow decomposition
in directions with long cylinders might not discover them if a limit
is set::

    sage: S = translation_surfaces.mcmullen_genus2_prototype(4,2,1,1,1/4)
    sage: l = S.base_ring().gen()
    sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf
    sage: dec = O.decomposition((8*l - 25, 16), 9) # optional: pyflatsurf
    sage: dec.undeterminedComponents() # optional: pyflatsurf
    [Component with perimeter [...]]

Further refinement might change the status of undetermined components::

    sage: import pyflatsurf # optional: pyflatsurf
    sage: dec.decompose(10r) # optional: pyflatsurf
    True
    sage: dec.undeterminedComponents() # optional: pyflatsurf
    []

Veech surfaces have the property that any saddle connection direction is
parabolic::

    sage: S = translation_surfaces.veech_double_n_gon(5)
    sage: O = GL2ROrbitClosure(S)  # optional: pyflatsurf
    sage: all(d.parabolic() for d in O.decompositions_depth_first(3))  # optional: pyflatsurf
    True

For surfaces in rank one loci, even though they are completely periodic,
they are generally not parabolic::

    sage: S = translation_surfaces.mcmullen_genus2_prototype(4,2,1,1,1/4)
    sage: O = GL2ROrbitClosure(S)  # optional: pyflatsurf
    sage: all((d.hasCylinder() == False) or d.parabolic() for d in O.decompositions(6))  # optional: pyflatsurf
    False
    sage: all((d.completelyPeriodic() == True) or (d.hasCylinder() == False) for d in O.decompositions(6))  # optional: pyflatsurf
    True
"""
# ****************************************************************************
#  This file is part of sage-flatsurf.
#
#        Copyright (C) 2019-2022 Julian Rüth
#                      2020      Vincent Delecroix
#
#  sage-flatsurf 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.
#
#  sage-flatsurf is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with sage-flatsurf. If not, see <https://www.gnu.org/licenses/>.
# ****************************************************************************

from sage.all import FreeModule, matrix, identity_matrix, ZZ, QQ, Unknown, vector, prod


[docs] class GL2ROrbitClosure: r""" Lower bound approximation to the tangent space of a GL(2,R)-orbit closure of a linear family of translation surfaces. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(3, 3, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: GL2ROrbitClosure(S) # optional: pyflatsurf GL(2,R)-orbit closure of dimension at least 2 in H_5(4, 2^2) (ambient dimension 12) Computing an orbit closure over an exact real ring with transcendental elements:: sage: from flatsurf import Polygon, EuclideanPolygonsWithAngles sage: from pyexactreal import ExactReals # optional: pyexactreal # random output due to matplotlib warnings with some combinations of setuptools and matplotlib sage: E = EuclideanPolygonsWithAngles((1, 5, 5, 5)) sage: R = ExactReals(E.base_ring()) # optional: pyexactreal sage: slopes = E.slopes() sage: T = Polygon(angles=(1, 5, 5, 5), edges=[slopes[0], R.random_element(1/4) * slopes[1]]) # optional: pyexactreal sage: S = similarity_surfaces.billiard(T) # optional: pyexactreal sage: S = S.minimal_cover(cover_type="translation") # optional: pyexactreal sage: O = GL2ROrbitClosure(S); O # optional: pyflatsurf, optional: pyexactreal GL(2,R)-orbit closure of dimension at least 4 in H_7(4^3, 0) (ambient dimension 17) sage: bound = E.billiard_unfolding_stratum('half-translation', marked_points=True).dimension() sage: for decomposition in O.decompositions(1): # long time, optional: pyflatsurf, optional: pyexactreal ....: O.update_tangent_space_from_flow_decomposition(decomposition) ....: if O.dimension() == bound: break sage: O # long time, optional: pyflatsurf, optional: pyexactreal GL(2,R)-orbit closure of dimension at least 8 in H_7(4^3, 0) (ambient dimension 17) TESTS:: sage: from flatsurf import translation_surfaces sage: x = polygen(QQ) sage: K = NumberField(x**3 - 2, 'a', embedding=AA(2)**QQ((1,3))) sage: a = K.gen() sage: S = translation_surfaces.mcmullen_genus2_prototype(2, 1, 0, -1, a/4) sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: O._U.base_ring() # optional: pyflatsurf Number Field in a with defining polynomial x^3 - 2 with a = 1.259921049894873? We now illustrate the projection:: sage: try: # optional: pyflatsurf ....: V = O.proj.row_ambient_module() ....: except AttributeError: ....: V = O.proj._row_ambient_module() sage: try: # optional: pyflatsurf ....: H = O.proj.column_ambient_module() ....: except AttributeError: ....: H = O.proj._column_ambient_module() sage: assert O.proj.rank() == H.dimension() # optional: pyflatsurf sage: for b in O.boundaries(): # optional: pyflatsurf ....: assert (O.proj * b).is_zero() sage: for i, e in enumerate(O.spanning_set): # optional: pyflatsurf ....: assert (O.proj * V.gen(e.index())) == H.gen(i) """ def __init__(self, surface): from flatsurf.geometry.categories import TranslationSurfaces from flatsurf.geometry.surface import Surface_base if isinstance(surface, Surface_base): if surface not in TranslationSurfaces(): raise NotImplementedError( "cannot compute orbit closure of a non-translation surface" ) base_ring = surface.base_ring() self._surface = surface.pyflatsurf().codomain().flat_triangulation() else: from flatsurf.geometry.pyflatsurf.conversion import sage_ring base_ring = sage_ring(surface) self._surface = surface # A model of the vector space R² in libflatsurf, e.g., to represent the # vector associated to a saddle connection. from flatsurf.features import pyflatsurf_feature pyflatsurf_feature.require() import pyflatsurf.vector self.V2 = pyflatsurf.vector.Vectors(base_ring) # We construct a spanning set of edges, that is a subset of the # edges that form a basis of H_1(S, Sigma; Z) # It comes together with a projection matrix t, m = self._spanning_tree() assert set(t.keys()) == {f[2] for f in self._surface.faces()} self.spanning_set = [] v = set(t.values()) for e in self._surface.edges(): if e.positive() not in v and e.negative() not in v: self.spanning_set.append(e) self.d = len(self.spanning_set) assert 3 * self.d - 3 == self._surface.size() assert m.rank() == self.d m = m.transpose() # projection matrix from Z^E to H_1(S, Sigma; Z) in the basis # of spanning edges self.proj = matrix(ZZ, [r for r in m.rows() if not r.is_zero()]) self.Omega = self._intersection_matrix(t, self.spanning_set) self.V = FreeModule(self.V2.base_ring(), self.d) self.H = matrix(self.V2.base_ring(), self.d, 2) for i in range(self.d): s = self._surface.fromHalfEdge(self.spanning_set[i].positive()) self.H[i] = self.V2._isomorphic_vector_space(self.V2(s)) self.Hdual = self.Omega * self.H # Note that we don't use Sage vector spaces because they are usually # way too slow (in particular we avoid calling .echelonize()) self._U = matrix(self.V2._algebraic_ring(), self.d) # Computing the rank of _U can be very expensive. Therefore, we use # that rank _Ubar ≤ rank _U where _Ubar = _U mod _p for some prime # ideal's valuation _p, see _rank(). # Since _p is a bit expensive to compute, it's initialized on demand. self._p = None self._Ubar = None self._U_rank = 0 self.update_tangent_space_from_vector(self.H.transpose()[0]) self.update_tangent_space_from_vector(self.H.transpose()[1])
[docs] def dimension(self): r""" Return the current complex dimension of the GL(2,R)-orbit closure. Note that this is not the dimension of the orbit closure but only a lower bound. It is always at least 2 (coming from a GL(2,R)-orbit). The current tangent space could be refined via :meth:`update_tangent_space_from_flow_decomposition`. EXAMPLES:: sage: from flatsurf import Polygon, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = Polygon(angles=(1, 3, 5)) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: O.dimension() # optional: pyflatsurf 2 sage: bound = T.category().billiard_unfolding_stratum('half-translation', marked_points=True).dimension() sage: for decomposition in O.decompositions(1): # long time, optional: pyflatsurf ....: if O.dimension() == bound: break ....: O.update_tangent_space_from_flow_decomposition(decomposition) ....: print(O.dimension()) 3 5 5 7 9 10 """ return self._U_rank
[docs] def ambient_stratum(self): r""" Return the stratum of Abelian differentials this surface belongs to. EXAMPLES:: sage: from flatsurf import Polygon, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = Polygon(angles=(1, 3, 5)) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: O.ambient_stratum() # optional: pyflatsurf H_3(4, 0^4) """ from surface_dynamics import Stratum surface = self._surface angles = [surface.angle(v) for v in surface.vertices()] return Stratum([a - 1 for a in angles], 1)
[docs] def base_ring(self): r""" Return the underlying base ring """ return self._U.base_ring()
[docs] def field_of_definition(self): r""" Return the field of definition of the current subspace. .. WARNING:: This involves the computation of the echelon form of the matrix. It might be rather expensive if the computation of the tangent space is not terminated. EXAMPLES:: sage: from flatsurf import Polygon, similarity_surfaces, EuclideanPolygonsWithAngles sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: from pyexactreal import ExactReals # optional: pyexactreal sage: E = EuclideanPolygonsWithAngles((1, 5, 5, 5)) sage: R = ExactReals(E.base_ring()) # optional: pyexactreal sage: slopes = E.slopes() sage: T = Polygon(angles=(1, 5, 5, 5), edges=[slopes[0], R.random_element(1/4) * slopes[1]]) # optional: pyexactreal sage: S = similarity_surfaces.billiard(T) # optional: pyexactreal sage: S = S.minimal_cover(cover_type="translation") # optional: pyexactreal sage: O = GL2ROrbitClosure(S); O # optional: pyflatsurf, optional: pyexactreal GL(2,R)-orbit closure of dimension at least 4 in H_7(4^3, 0) (ambient dimension 17) sage: O.field_of_definition() # optional: pyflatsurf, optional: pyexactreal Number Field in c0 with defining polynomial x^2 - 2 with c0 = 1.414213562373095? sage: bound = E.billiard_unfolding_stratum('half-translation', marked_points=True).dimension() sage: for decomposition in O.decompositions(1): # long time, optional: pyflatsurf, optional: pyexactreal ....: if O.dimension() == bound: break ....: O.update_tangent_space_from_flow_decomposition(decomposition) sage: O.field_of_definition() # long time, optional: pyflatsurf, optional: pyexactreal Rational Field sage: T = Polygon(angles=(1, 3, 5)) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: O.field_of_definition() # optional: pyflatsurf Number Field in c0 with defining polynomial x^3 - 3*x - 1 with c0 = 1.879385241571817? sage: bound = T.category().billiard_unfolding_stratum('half-translation', marked_points=True).dimension() sage: for decomposition in O.decompositions(1): # long time, optional: pyflatsurf ....: if O.dimension() == bound: break ....: O.update_tangent_space_from_flow_decomposition(decomposition) sage: O.field_of_definition() # long time, optional: pyflatsurf Rational Field """ M = self._U.echelon_form() from flatsurf.geometry.subfield import subfield_from_elements L, elts, phi = subfield_from_elements(M.base_ring(), M[: self._U_rank].list()) return L
def _half_edge_to_face(self, h): r""" Return a canonical half-edge encoding the face bounded by ``h``. """ surface = self._surface h1 = h h2 = surface.nextInFace(h1) h3 = surface.nextInFace(h2) return min([h1, h2, h3], key=lambda x: x.index()) def __repr__(self): return ( "GL(2,R)-orbit closure of dimension at least %d in %s (ambient dimension %d)" % (self._U_rank, self.ambient_stratum(), self.d) )
[docs] def holonomy(self, v): r""" Return the holonomy of ``v`` (with respect to the chosen homology basis) EXAMPLES:: sage: from flatsurf import translation_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: x = polygen(QQ) sage: K.<a> = NumberField(x^3 - 2, embedding=AA(2)**(1/3)) sage: S = translation_surfaces.mcmullen_L(1,1,1,a) sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: edges = O._surface.edges() # optional: pyflatsurf sage: F = FreeModule(ZZ, len(edges)) # optional: pyflatsurf sage: all(O.V2(O.holonomy(O.proj * F.gen(i))).vector == O.V2(O._surface.fromHalfEdge(e.positive())).vector for i, e in enumerate(edges)) # optional: pyflatsurf True """ return self.V(v) * self.H
[docs] def holonomy_dual(self, v): r""" Return the holonomy of the dual of ``v`` """ return self.V(v) * self.Hdual
[docs] def tangent_space_basis(self): return self._U[: self._U_rank].rows()
[docs] def lift(self, v): r""" Given a vector in the "spanning set basis" return a vector on the full basis of edges. EXAMPLES:: sage: from flatsurf import polygons, translation_surfaces, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: S = translation_surfaces.mcmullen_genus2_prototype(4,2,1,1,0) sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: u0,u1 = O.tangent_space_basis() # optional: pyflatsurf sage: v0 = O.lift(u0) # optional: pyflatsurf sage: v1 = O.lift(u1) # optional: pyflatsurf sage: span([v0, v1]) # optional: pyflatsurf Vector space of degree 9 and dimension 2 over Real Embedded Number Field in l with defining polynomial x^2 - x - 8 with l = 3.372281323269015? Basis matrix: [ 1 0 -1 (1/8*l+7/8 ~ 1.2965352) (-1/8*l+1/8 ~ -0.29653517) -1 (5/8*l-5/8 ~ 1.4826758) (-1/2*l+3/2 ~ -0.18614066) (-5/8*l+13/8 ~ -0.48267583)] [ 0 1 -1 (1/4*l-1/4 ~ 0.59307033) (-1/4*l+1/4 ~ -0.59307033) 0 (1/4*l-1/4 ~ 0.59307033) 0 (-1/4*l+1/4 ~ -0.59307033)] This can be used to deform the surface:: sage: T = polygons.triangle(3,4,13) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover("translation").erase_marked_points() # long time (3s, #122), optional: pyflatsurf sage: O = GL2ROrbitClosure(S) # long time (above), optional: pyflatsurf sage: for d in O.decompositions(4, 20): # long time (2s, #124), optional: pyflatsurf ....: O.update_tangent_space_from_flow_decomposition(d) ....: if O.dimension() == 4: ....: break sage: d1,d2,d3,d4 = [O.lift(b) for b in O.tangent_space_basis()] # long time (above), optional: pyflatsurf sage: dreal = d1/132 + d2/227 + d3/1280 - d4/13201 # long time (above), optional: pyflatsurf sage: dimag = d1/141 - d2/233 + d4/1230 + d4/14250 # long time (above), optional: pyflatsurf sage: d = [O.V2((x,y)).vector for x,y in zip(dreal,dimag)] # long time (above), optional: pyflatsurf sage: S2 = O._surface + d # long time (6s), optional: pyflatsurf sage: O2 = GL2ROrbitClosure(S2.surface()) # long time (above), optional: pyflatsurf sage: for d in O2.decompositions(1, 20): # long time (25s, #124), optional: pyflatsurf ....: O2.update_tangent_space_from_flow_decomposition(d) """ # given the values on the spanning edges we reconstruct the unique vector that # vanishes on the boundary bdry = self.boundaries() n = self._surface.edges().size() k = len(self.spanning_set) assert k + len(bdry) == n + 1 A = matrix(QQ, n + 1, n) for i, e in enumerate(self.spanning_set): A[i, e.index()] = 1 for i, b in enumerate(bdry): A[k + i, :] = b u = vector(self.V2.base_ring(), n + 1) u[:k] = v from sage.all import Fields if not self.V2.base_ring() in Fields(): assert all(uu._backend.coefficients().size() == 1 for uu in u) u = u.parent().change_ring(self.V2.base_ring().base_ring())( [uu._backend.coefficients()[0] for uu in u] ) return A.solve_right(u)
[docs] def absolute_homology(self): vert_index = {v: i for i, v in enumerate(self._surface.vertices())} m = len(vert_index) if m == 1: return self.V rows = [] from flatsurf.features import pyflatsurf_feature pyflatsurf_feature.require() import pyflatsurf for e in self.spanning_set: r = [0] * m i = vert_index[ pyflatsurf.flatsurf.Vertex.target( e.positive(), self._surface.combinatorial() ) ] j = vert_index[ pyflatsurf.flatsurf.Vertex.source( e.positive(), self._surface.combinatorial() ) ] if i != j: r[i] = 1 r[j] = -1 rows.append(r) return matrix(rows).left_kernel()
[docs] def absolute_dimension(self): r""" EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1,3,4) # Veech octagon sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover("translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: O.absolute_dimension() # optional: pyflatsurf 2 The triangular billiard (5,6,7) belongs to the canonical double cover of the stratum Q(5,3,0^3) in genus 3. The orbit is dense and we can check that the absolute dimension is indeed `6 = 2 rank`:: sage: T = polygons.triangle(5,6,7) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover("translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: for d in O.decompositions(5, 100): # long time (3s) # optional: pyflatsurf ....: O.update_tangent_space_from_flow_decomposition(d) ....: if O.dimension() == 9: ....: break sage: O.absolute_dimension() # long time (above) # optional: pyflatsurf 6 """ return ( self.absolute_homology().matrix() * self._U[: self._U_rank].transpose() ).rank()
def _spanning_tree(self, root=None): r""" Return a pair ``(tree, proj)`` where - ``tree`` is a tree encoded in a dictionary. Its keys are the faces (coded by their minimal adjacent half-edge) and the corresponding value is the half-edge to cross to go toward the root face. - ``proj`` a projection matrix : for a vector ``v``, the vector ``v * proj`` is cohomologous to ``v`` and only takes values on the spanning set. EXAMPLES: The following example illustrate the projection matrix:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1,3,4) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover("translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: num_edges = O._surface.edges().size() # optional: pyflatsurf sage: V = VectorSpace(QQ, num_edges) # optional: pyflatsurf sage: tree, proj = O._spanning_tree() # optional: pyflatsurf The matrix ``proj`` is indeed a projection:: sage: proj * proj == proj # optional: pyflatsurf True We now check that the matrix maps vectors to cohomologous vectors and takes values only on the chosen spanning set of edges:: sage: values = tree.values() # optional: pyflatsurf sage: indices = set(e.index() for e in O._surface.edges() if e.positive() not in values and e.negative() not in values) # optional: pyflatsurf sage: B = V.subspace(O.boundaries()) # optional: pyflatsurf sage: for e in range(num_edges): # optional: pyflatsurf ....: v = V.gen(e) ....: assert (v * proj - v) in B ....: if e in indices: ....: assert v * proj == v ....: else: ....: assert (proj * v).is_zero() """ if root is None: root = next(iter(self._surface.edges())).positive() root = self._half_edge_to_face(root) t = {root: None} # face -> half edge to take to go to the root todo = [root] edges = [] # store edges in topological order to perform Gauss reduction while todo: f = todo.pop() for _ in range(3): f1 = -f g = self._half_edge_to_face(f1) if g not in t: t[g] = f1 todo.append(g) edges.append(f1) f = self._surface.nextInFace(f) # gauss reduction n = self._surface.size() proj = identity_matrix(ZZ, n) edges.reverse() for f1 in edges: f2 = self._surface.nextInFace(f1) f3 = self._surface.nextInFace(f2) assert self._surface.nextInFace(f3) == f1 i1 = f1.index() s1 = -1 if i1 % 2 else 1 i2 = f2.index() s2 = -1 if i2 % 2 else 1 i3 = f3.index() s3 = -1 if i3 % 2 else 1 i1 = f1.edge().index() i2 = f2.edge().index() i3 = f3.edge().index() proj[i1] = -s1 * (s2 * proj[i2] + s3 * proj[i3]) for j in range(n): assert proj[j, i1] == 0 return (t, proj) def _intersection_matrix(self, t, spanning_set): r""" Given a spanning tree, compute the associated intersection matrix. It can be used to compute holonomies. (we can be off by a - sign) """ d = len(spanning_set) h = spanning_set[0].positive() all_edges = {e.positive() for e in spanning_set} all_edges.update([e.negative() for e in spanning_set]) contour = [] contour_inv = {} # half edge -> position in contour while h not in contour_inv: contour_inv[h] = len(contour) contour.append(h) h = self._surface.nextAtVertex(-h) while h not in all_edges: h = self._surface.nextAtVertex(h) assert len(contour) == len(all_edges) # two curves intersect when their relative position in the contour # are x y x y or y x y x Omega = matrix(ZZ, d) for i in range(len(spanning_set)): ei = spanning_set[i] pi1 = contour_inv[ei.positive()] pi2 = contour_inv[ei.negative()] if pi1 > pi2: si = -1 pi1, pi2 = pi2, pi1 else: si = 1 for j in range(i + 1, len(spanning_set)): ej = spanning_set[j] pj1 = contour_inv[ej.positive()] pj2 = contour_inv[ej.negative()] if pj1 > pj2: sj = -1 pj1, pj2 = pj2, pj1 else: sj = 1 # pj1 pj2 pi1 pi2: pj2 < pi1 # pi1 pi2 pj1 pj2: pi2 < pj1 # pi1 pj1 pj2 pi2: pi1 < pj1 and pj2 < pi2 # pj1 pi1 pi2 pj2: pj1 < pi1 and pi2 < pj2 if ( (pj2 < pi1) or (pi2 < pj1) or (pj1 > pi1 and pj2 < pi2) or (pj1 < pi1 and pj2 > pi2) ): # no intersection continue if pi1 < pj1 < pi2: # one sign Omega[i, j] = si * sj else: # other sign assert pi1 < pj2 < pi2, (pi1, pi2, pj1, pj2) Omega[i, j] = -si * sj Omega[j, i] = -Omega[i, j] return Omega
[docs] def boundaries(self): r""" Return the list of boundaries (ie sum of edges around a triangular face). These are elements of H_1(S, Sigma; Z). TESTS:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: from itertools import product sage: for a in range(1,5): # long time (1s), optional: pyflatsurf ....: for b in range(a, 5): ....: for c in range(b, 5): ....: if gcd([a, b, c]) != 1 or (a,b,c) == (1,1,2): ....: continue ....: T = polygons.triangle(a, b, c) ....: S = similarity_surfaces.billiard(T) ....: S = S.minimal_cover(cover_type="translation") ....: O = GL2ROrbitClosure(S) ....: for b in O.boundaries(): ....: assert (O.proj * b).is_zero() """ n = self._surface.size() V = FreeModule(ZZ, n) B = [] for f1, f2, f3 in self._surface.faces(): i1 = f1.index() s1 = -1 if i1 % 2 else 1 i2 = f2.index() s2 = -1 if i2 % 2 else 1 i3 = f3.index() s3 = -1 if i3 % 2 else 1 i1 = f1.edge().index() i2 = f2.edge().index() i3 = f3.edge().index() v = [0] * n v[i1] = 1 v[i2] = s1 * s2 v[i3] = s1 * s3 B.append(V(v)) B[-1].set_immutable() return B
[docs] def decomposition(self, v, limit=-1): v = self.V2(v) from flatsurf.features import pyflatsurf_feature pyflatsurf_feature.require() import pyflatsurf decomposition = pyflatsurf.flatsurf.makeFlowDecomposition( self._surface, v.vector ) if limit != 0: decomposition.decompose(int(limit)) return decomposition
[docs] def decompositions(self, bound, limit=-1, bfs=False): limit = int(limit) connections = self._surface.connections().bound(int(bound)) if bfs: connections = connections.byLength() slopes = None from flatsurf.features import cppyy_feature cppyy_feature.require() import cppyy for connection in connections: direction = connection.vector() if slopes is None: slopes = cppyy.gbl.std.set[ type(direction), type(direction).CompareSlope ]() if slopes.find(direction) != slopes.end(): continue slopes.insert(direction) yield self.decomposition(direction, limit)
[docs] def decompositions_depth_first(self, bound, limit=-1): return self.decompositions(bound, bfs=False, limit=limit)
[docs] def decompositions_breadth_first(self, bound, limit=-1): return self.decompositions(bound, bfs=True, limit=limit)
[docs] def is_teichmueller_curve(self, bound, limit=-1): r""" Return ``False`` when the program can find a direction which is either completely periodic with incomensurable moduli or a direction with at least one cylinder and at least one minimal component. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: for a in range(1,6): # long time (8s), optional: pyflatsurf ....: for b in range(a,6): ....: for c in range(b,6): ....: if a + b + c > 7 or gcd([a,b,c]) != 1: ....: continue ....: T = polygons.triangle(a, b, c) ....: S = similarity_surfaces.billiard(T) ....: S = S.minimal_cover(cover_type="translation") ....: O = GL2ROrbitClosure(S) ....: if O.is_teichmueller_curve(3, 50) != False: ....: print(a,b,c) 1 1 1 1 1 2 1 1 4 1 2 2 1 2 3 1 3 3 """ if self.V2.base_ring in [ZZ, QQ]: # square tiled surface return True if self.dimension() > 2: # too large orbit closure return False k = self.field_of_definition() if k in [ZZ, QQ]: # square tiled return True nv = len(self._surface.vertices()) ne = len(self._surface.edges()) nf = len(self._surface.faces()) genus = (ne - nv - nf) // 2 + 1 if k.degree() > genus or not k.is_totally_real(): return False for decomposition in self.decompositions_depth_first(bound, limit): if ( decomposition.parabolic() == False ): # noqa, we are comparing to a boost tribool so this cannot be replaced by "is False" return False return Unknown
[docs] def cylinder_circumference(self, component, A, sc_index, proj): r""" Return the circumference of the cylinder ``component`` in the homology of the underlying surface. INPUT: - ``component`` -- a cylinder - ``A``, ``sc_index``, ``proj`` -- the output of ``flow_decomposition_kontsevich_zorich_cocycle`` EXAMPLES:: sage: from flatsurf import translation_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: S = translation_surfaces.veech_double_n_gon(5) sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: for dec in O.decompositions_depth_first(1): # optional: pyflatsurf ....: break sage: c0, c1 = dec.components() # optional: pyflatsurf sage: kz = O.flow_decomposition_kontsevich_zorich_cocycle(dec) # optional: pyflatsurf sage: O.cylinder_circumference(c0, *kz) # optional: pyflatsurf (1, 1, 1, -1) sage: O.cylinder_circumference(c1, *kz) # optional: pyflatsurf (0, 0, 1, 0) """ if ( component.cylinder() != True ): # noqa, we are comparing to a boost tribool so this cannot be replaced by "is not True" raise ValueError perimeters = list(component.perimeter()) per = perimeters[0] assert not per.vertical() sc = per.saddleConnection() i = sc_index[sc] if i < 0: s = -1 i = -i - 1 else: s = 1 v = s * proj.column(i) circumference = -A.solve_right(v) # check hol = self.holonomy_dual(circumference) holbis = component.circumferenceHolonomy() holbis = self.V2._isomorphic_vector_space(self.V2(holbis)) assert hol == holbis, (hol, holbis) return circumference
[docs] def cylinder_deformation_subspace(self, decomposition): r""" Return a subspace included in the tangent space to the GL(2,R)-orbit closure computed from the flow decomposition ``decomposition``. From A. Wright cylinder deformation Theorem. """ def eliminate_denominators(fractions): r""" Given a list of ``fractions``, pairs of numerators `n_i` and denominators `d_i`, return a list of fractions `c n_i/d_i` scaled uniformly such that the value can be represented in the underlying ring. """ fractions = list(fractions) try: return [x.parent()(x / y) for x, y in fractions] except (ValueError, ArithmeticError, NotImplementedError): denominators = {denominator for numerator, denominator in fractions} return [ numerator * prod([d for d in denominators if denominator != d]) for (numerator, denominator) in fractions ] module_fractions = [] vcyls = [] kz = self.flow_decomposition_kontsevich_zorich_cocycle(decomposition) for component in decomposition.components(): if ( component.cylinder() == False ): # noqa, we are comparing to a boost tribool so this cannot be replaced by "is False" continue elif ( component.cylinder() == True ): # noqa, we are comparing to a boost tribool so this cannot be replaced with "is True" vcyls.append(self.cylinder_circumference(component, *kz)) width = self.V2._isomorphic_vector_space.base_ring()( self.V2.base_ring()(component.width()) ) height = self.V2._isomorphic_vector_space.base_ring()( self.V2.base_ring()( component.vertical().project(component.circumferenceHolonomy()) ) ) module_fractions.append((width, height)) else: return [] if not module_fractions: return [] modules = eliminate_denominators(module_fractions) if hasattr(modules[0], "_backend"): # Make sure all modules live in the same K-Module so that .coefficients() below produces coefficient lists of the same length. from functools import reduce parent = reduce( lambda m, n: m.span(m, n), [module._backend.module() for module in modules], modules[0]._backend.module(), ) modules = [ module.parent()(module._backend.promote(parent)) for module in modules ] def to_rational_vector(x): r""" Return the rational coefficients of `x` over its implicit basis, e.g., if `x` is in a number field K, return the coefficients of x in `K` as a vector space over the rationals. """ if x.parent() in [ZZ, QQ]: ret = [QQ(x)] elif hasattr(x, "vector"): ret = x.vector() elif hasattr(x, "renf_elem"): ret = x.parent().number_field(x).vector() elif hasattr(x, "_backend"): from itertools import chain ret = list( chain( *[ to_rational_vector( self.V2.base_ring().base_ring()( self.V2.base_ring().base_ring()(c) ) ) for c in x._backend.coefficients() ] ) ) else: raise NotImplementedError( "cannot turn {}, i.e., a {}, into a rational vector yet".format( x, type(x) ) ) assert all(y in QQ for y in ret) return ret assert all(module.parent() is modules[0].parent() for module in modules) M = matrix([to_rational_vector(module) for module in modules]) assert M.base_ring() is QQ relations = self._left_kernel_matrix(M) assert len(vcyls) == len(module_fractions) == relations.ncols() vectors = [ sum( t * vcyl * module[1] for (t, vcyl, module) in zip(relation, vcyls, module_fractions) ) for relation in self._right_kernel_matrix(relations).rows() ] assert all( v.base_ring() is self.V2._isomorphic_vector_space.base_ring() for v in vectors ) return vectors
def _flow_decomposition_spanning_tree(self, decomposition, sc_index, sc_comp): r""" Helper for :meth:`flow_decomposition_kontsevich_zorich_cocycle`. The method is similar to `_spanning_tree` but adapted to the basis used in the flow decomposition ``decomposition``. It returns a pair ``(tree, proj)`` where ``tree`` is a tree encoded in a dictionary and ``proj`` is a projection matrix. INPUT: - ``decomposition`` -- a flow decomposition - ``sc_index`` -- a dictionary whose keys are the saddle connections on the boundary of ``decomposition`` and the corresponding values are the indices used in the matrix - ``sc_comp`` -- a dictionary whose keys are the saddle connections on the boundary of ``decomposition`` and the corresponding values are the indices of the components of ``decomposition`` """ components = list(decomposition.components()) n = len(sc_index) assert n % 2 == 0 n //= 2 for p in components[0].perimeter(): break t = {0: None} # face -> half edge to take to go to the root todo = [0] edges = [] # store edges in topological order to perform Gauss reduction while todo: i = todo.pop() c = components[i] for sc in c.perimeter(): sc1 = -sc.saddleConnection() j = sc_comp[sc1] if j not in t: t[j] = sc1 todo.append(j) edges.append(sc1) # gauss reduction spanning_set = set(range(n)) proj = identity_matrix(ZZ, n) edges.reverse() for sc1 in edges: i1 = sc_index[sc1] if i1 < 0: s1 = -1 i1 = -i1 - 1 else: s1 = 1 comp = components[sc_comp[sc1]] proj[i1] = 0 for p in comp.perimeter(): sc = p.saddleConnection() if sc == sc1: continue j = sc_index[sc] if j < 0: s = -1 j = -j - 1 else: s = 1 proj[i1] = -s1 * s * proj[j] spanning_set.remove(i1) for j in range(n): assert proj[j, i1] == 0 return (t, sorted(spanning_set), proj)
[docs] def flow_decomposition_kontsevich_zorich_cocycle(self, decomposition): r""" Base change from the homology of ``decomposition`` to the underlying surface. Return a triple ``A``, ``sc_index``, ``proj`` where - ``A`` is a `d \times d` matrix where `d` is the dimension of the ambient stratum (ie the dimension of the homology relative to the singularities) - ``sc_index`` is a dictionary whose keys are saddle connections appearing in the various components of this decomposition and the corresponding value is an integral index - ``proj`` is a projection matrix from vectors expressed as a sum of saddle connections appearing in this component to the chosen representative of cohomology with respect to this component EXAMPLES:: sage: from flatsurf import translation_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: S = translation_surfaces.square_torus() sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: for dec in O.decompositions_depth_first(3): # optional: pyflatsurf ....: kz = O.flow_decomposition_kontsevich_zorich_cocycle(dec) # optional: pyflatsurf ....: print(kz[0]) [ 0 1] [-1 -1] [ 0 1] [-1 -2] [-1 -1] [ 1 0] [-1 -1] [ 2 1] [-2 -1] [ 3 1] [-1 -1] [ 3 2] [1 0] [0 1] [ 1 0] [-1 1] """ sc_pos = ( [] ) # list of positive boundary saddle connections (store only one orientation for each) sc_index = ( {} ) # inverse of sc_pos (and also reverse orientation with negative index) sc_comp = {} n_saddles = 0 components = list(decomposition.components()) n_components = len(components) for i, comp in enumerate(components): for p in comp.perimeter(): sc = p.saddleConnection() sc_comp[sc] = i if sc not in sc_index: sc_index[sc] = n_saddles sc_pos.append(sc) sc_index[-sc] = -n_saddles - 1 n_saddles += 1 t, spanning_set, proj = self._flow_decomposition_spanning_tree( decomposition, sc_index, sc_comp ) assert proj.rank() == len(spanning_set) == n_saddles - n_components + 1 proj = proj.transpose() proj = matrix(ZZ, [r for r in proj.rows() if not r.is_zero()]) assert proj.nrows() == self.proj.nrows() # Write the base change V^*(T') -> V^*(T) relative to our bases A = matrix(ZZ, self.d) for i, sc in enumerate(spanning_set): sc = sc_pos[sc] c = sc.chain() for edge in self._surface.edges(): A[i] += ZZ(str(c[edge])) * self.proj.column(edge.index()) assert A.det().is_unit() return A, sc_index, proj
[docs] def update_tangent_space_from_flow_decomposition(self, decomposition): r""" Update the current tangent space by using the cylinder deformation vectors from ``decomposition``. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1, 2, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: for d in O.decompositions(1): # long time (1s), optional: pyflatsurf ....: O.update_tangent_space_from_flow_decomposition(d) sage: assert O.dimension() == 2 # long time (above), optional: pyflatsurf TESTS: A regression test for https://github.com/flatsurf/sage-flatsurf/pull/69:: sage: from itertools import islice sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: for (a,b,c,dim) in [(3,2,2,7), (4,2,1,8), (4,4,3,11), (5,3,3,11), (5,4,4,13), (5,5,3,13)]: # long time, optional: pyflatsurf ....: T = polygons.triangle(a, b, c) ....: S = similarity_surfaces.billiard(T) ....: S = S.minimal_cover(cover_type="translation") ....: O = GL2ROrbitClosure(S) # optional: pyflatsurf ....: nsteps = 0 ....: for d in islice(O.decompositions(3,100), 10): # optional: pyflatsurf ....: O.update_tangent_space_from_flow_decomposition(d) ....: nsteps += 1 ....: if O.dimension() == dim: ....: break ....: print("(%d, %d, %d): %d" % (a, b, c, nsteps)) (3, 2, 2): 5 (4, 2, 1): 4 (4, 4, 3): 4 (5, 3, 3): 4 (5, 4, 4): 7 (5, 5, 3): 4 """ if self._U_rank == self._U.nrows(): return for v in self.cylinder_deformation_subspace(decomposition): self.update_tangent_space_from_vector(v) if self._U_rank == self._U.nrows(): return
def _rank(self): r""" Return a lower bound for the rank of the matrix _U. """ while True: if self._p is not None: try: self._Ubar = self._U.apply_map( phi=self._p.reduce, R=self._p.residue_field() ) break except ValueError: # The matrix _U cannot be reduced mod p because some entries have negative valuation. pass p = 2**30 if self._p is None else self._p.p() from sage.all import next_prime self._p = ZZ.valuation(next_prime(p)).extensions(self._U.base_ring())[0] return self._Ubar.rank()
[docs] def update_tangent_space_from_vector(self, v): if self._U_rank == self._U.nrows(): return v = vector(v) if v.base_ring() is not self.V2._algebraic_ring(): for gen, p in self.V2.decomposition(v): self.update_tangent_space_from_vector(p) return self._U[self._U_rank] = v rank = self._rank() if rank > self._U_rank: assert rank == self._U_rank + 1 self._U_rank += 1
def __eq__(self, other): r""" Return whether ``other`` was built starting from the same surface than this orbit closure. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1, 2, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: GL2ROrbitClosure(S) == GL2ROrbitClosure(S) # optional: pyflatsurf True """ return self._surface == other._surface def __ne__(self, other): r""" Return whether ``other`` was not built starting from the same surface than this orbit closure. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1, 2, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: GL2ROrbitClosure(S) != GL2ROrbitClosure(S) # optional: pyflatsurf False """ return not (self == other) def __hash__(self): r""" Return a hash value of this object compatible with equality operators. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1, 2, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: hash(O) == hash(O) # optional: pyflatsurf True """ return hash(self._surface) def __reduce__(self): r""" Return a serializable representation of this Orbit Closure. EXAMPLES:: sage: from flatsurf import polygons, similarity_surfaces sage: from flatsurf import GL2ROrbitClosure # optional: pyflatsurf sage: T = polygons.triangle(1, 2, 5) sage: S = similarity_surfaces.billiard(T) sage: S = S.minimal_cover(cover_type="translation") sage: O = GL2ROrbitClosure(S) # optional: pyflatsurf sage: loads(dumps(O)) == O # optional: pyflatsurf # long time (6s, #123) True """ return ( GL2ROrbitClosure, (self._surface,), {"_U": self._U, "_U_rank": self._U_rank}, ) @classmethod def _right_kernel_matrix(cls, M): r""" Compute the right kernel of the rational matrix `M`. See https://github.com/flatsurf/sage-flatsurf/issues/100. """ M = M._clear_denom()[0] rows = M.nrows() columns = M.ncols() # The backends cannot handle these trivial cases if M.ncols() == 0: return M.new_matrix(nrows=0, ncols=M.ncols()) if M.nrows() == 0: return M.matrix_space(M.ncols(), M.ncols()).identity_matrix() height = M.height() if height >= 2**256: algorithm = "flint" elif columns >= 64 and rows >= 64: algorithm = "padic" elif rows * columns <= 32: algorithm = "flint" else: algorithm = "pari" if algorithm == "pari": ker = M.__pari__().matker().mattranspose().sage() elif algorithm == "flint": ker = M._rational_kernel_flint().transpose() else: ker = M._rational_kernel_iml().transpose() if ker.nrows() == 0 and ker.ncols() != M.ncols(): ker = ker.new_matrix(nrows=0, ncols=M.ncols()) return ker @classmethod def _left_kernel_matrix(cls, M): r""" Compute the left kernel of the rational matrix `M`. See https://github.com/flatsurf/sage-flatsurf/issues/100. """ return cls._right_kernel_matrix(M.transpose())