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
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())