r"""
Cells for the L-infinity Delaunay triangulation.
Each cell of the L-infinity Delaunay triangulation can be identified with a
marked triangulation. The marking corresponds to the position of the horizontal
and vertical separatrices. Each triangle hence get one of the following types:
bottom-left, bottom-right, top-left, top-right.
"""
# ****************************************************************************
# This file is part of sage-flatsurf.
#
# Copyright (C) 2016-2019 Vincent Delecroix
# 2016-2019 W. Patrick Hooper
# 2023 Julian Rüth
#
# 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.misc.cachefunc import cached_method
# the types of edges
V_NONE = 0 # start vertex has no horizontal/vertical separatrix
V_LEFT = 1 # horizontal separatrix going right
V_RIGHT = 2 # horizontal separatrix going left
V_BOT = 3 # vertical separatrix going up
V_TOP = 4 # vertical separatrix going down
[docs]
def sign_and_norm_conditions(dim, i, s):
r"""
Inequalities:
if `s=+1`, encode the inequalities `+1 >= x_i >= 0`
if `s=-1`, encode the inequalities `-1 <= x_i <= 0`
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import sign_and_norm_conditions
sage: sorted(Polyhedron(ieqs=sign_and_norm_conditions(1, 0, 1)).vertices_list())
[[0], [1]]
sage: sorted(Polyhedron(ieqs=sign_and_norm_conditions(1, 0, -1)).vertices_list())
[[-1], [0]]
sage: ieqs = []
sage: ieqs.extend(sign_and_norm_conditions(2, 0, 1))
sage: ieqs.extend(sign_and_norm_conditions(2, 1, -1))
sage: sorted(Polyhedron(ieqs=ieqs).vertices_list())
[[0, -1], [0, 0], [1, -1], [1, 0]]
"""
l_sign = [0] * (dim + 1)
l_sign[i + 1] = s
l_norm = [0] * (dim + 1)
l_norm[i + 1] = -s
l_norm[0] = 1
return (l_sign, l_norm)
[docs]
def opposite_condition(dim, i, j):
r"""
encode the equality `x_i = -x_j`
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import \
....: opposite_condition, sign_and_norm_conditions
sage: eq1 = opposite_condition(2, 0, 1)
sage: eq2 = opposite_condition(2, 1, 0)
sage: ieqs1 = sign_and_norm_conditions(2, 0, 1)
sage: ieqs2 = sign_and_norm_conditions(2, 1, -1)
sage: sorted(Polyhedron(eqns=[eq1], ieqs=ieqs1).vertices_list())
[[0, 0], [1, -1]]
sage: sorted(Polyhedron(eqns=[eq1,eq2], ieqs=ieqs1).vertices_list())
[[0, 0], [1, -1]]
sage: sorted(Polyhedron(eqns=[eq1,eq2], ieqs=ieqs1+ieqs2).vertices_list())
[[0, 0], [1, -1]]
"""
linear = [0] * (dim + 1)
linear[i + 1] = 1
linear[j + 1] = 1
return linear
[docs]
def bottom_top_delaunay_condition(dim, p1, e1, p2, e2):
r"""
Delaunay condition for bottom-top pairs of triangles
"""
# re(e2p1) <= im(e1) + im(e2m1)
e2m1 = (e2 + 2) % 3
e2p1 = (e2 + 1) % 3
im_e2m1 = 2 * (3 * p2 + e2m1) + 1
re_e2p1 = 2 * (3 * p2 + e2p1)
im_e1 = 2 * (3 * p1 + e1) + 1
linear = [0] * (dim + 1)
linear[im_e1 + 1] = 1
linear[im_e2m1 + 1] = 1
linear[re_e2p1 + 1] = -1
return linear
[docs]
def right_left_delaunay_condition(dim, p1, e1, p2, e2):
r"""
Delaunay condition for right-left pairs of triangles
"""
# im(e2p1) <= re(e2) + re(e1m1)
e1m1 = (e1 + 2) % 3
e2p1 = (e2 + 1) % 3
im_e2p1 = 3 * (p2 + e2p1) + 1
re_e2 = 3 * (p2 + e2)
re_e1m1 = 3 * (p1 + e1m1)
linear = [0] * (dim + 1)
linear[re_e2 + 1] = 1
linear[re_e1m1 + 1] = 1
linear[im_e2p1 + 1] = -1
return linear
[docs]
class LInfinityMarkedTriangulation:
r"""
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import \
....: V_NONE, V_BOT, V_TOP, V_RIGHT, V_LEFT, LInfinityMarkedTriangulation
sage: gluings = {(0,0):(1,2), (1,2):(0,0), (0,1):(1,0), (1,0):(0,1),
....: (0,2):(1,1), (1,1):(0,2)}
sage: types = [(V_BOT, V_NONE, V_RIGHT), (V_NONE, V_LEFT, V_TOP)]
sage: T = LInfinityMarkedTriangulation(2, gluings, types)
"""
def __init__(self, num_faces, edge_identifications, edge_types, check=True):
from sage.rings.integer_ring import ZZ
self._n = ZZ(num_faces)
self._edge_identifications = edge_identifications
self._edge_types = edge_types
if check:
self._check()
def _check(self):
if self._n % 2:
raise ValueError("the number of faces must be even")
if sorted(self._edge_identifications.keys()) != [
(i, j) for i in range(self._n) for j in range(3)
]:
raise ValueError("should be a triangulation")
if (
not isinstance(self._edge_types, list)
or len(self._edge_types) != self.num_faces()
):
raise ValueError("edge_types invalid")
for i in range(self.num_faces()):
if len(self._edge_types[i]) != 3:
raise ValueError("edge_types invalid")
for j in range(3):
if self._edge_types[i][j] not in [
V_NONE,
V_LEFT,
V_RIGHT,
V_BOT,
V_TOP,
]:
raise ValueError("types[{}] = {} invalid", i, self._edge_types[i])
seen = [False] * self._n
for p in range(self._n):
if seen[p]:
continue
sh = sum(
self._edge_types[p][r] == V_LEFT or self._edge_types[p][r] == V_RIGHT
for r in (0, 1, 2)
)
sv = sum(
self._edge_types[p][r] == V_BOT or self._edge_types[p][r] == V_TOP
for r in (0, 1, 2)
)
if sh != 1 or sv != 1:
raise ValueError(
"triangle {} has invalid types {}".format(p, self._edge_types[p])
)
[docs]
def num_faces(self):
return self._n
[docs]
def num_edges(self):
return 3 * self._n // 2
[docs]
def opposite_edge(self, p, e):
return self._edge_identifications[(p, e)]
def __repr__(self):
return "Marked triangulation made of {} triangles".format(self.num_faces())
[docs]
def bottom_top_pairs(self):
r"""
Return a list ``(p1,e1,p2,e2)``.
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import \
....: V_NONE, V_BOT, V_TOP, V_RIGHT, V_LEFT, LInfinityMarkedTriangulation
sage: gluings = {(0,0):(1,2), (1,2):(0,0), (0,1):(1,0), (1,0):(0,1),
....: (0,2):(1,1), (1,1):(0,2)}
sage: types = [(V_BOT, V_NONE, V_RIGHT), (V_NONE, V_LEFT, V_TOP)]
sage: T = LInfinityMarkedTriangulation(2, gluings, types)
sage: T.bottom_top_pairs()
[(0, 0, 1, 2)]
"""
pairs = []
for p1 in range(self._n):
for e1 in range(3):
if self._edge_types[p1][e1] == V_BOT:
e1p1 = (e1 + 1) % 3
p2, e2 = self.opposite_edge(p1, e1p1)
e2m1 = (e2 - 1) % 3
if self._edge_types[p2][e2m1] == V_TOP:
pairs.append((p1, e1, p2, e2m1))
return pairs
[docs]
def right_left_pairs(self):
r"""
Return a list ``(p1,e1,p2,e2)``
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import \
....: V_NONE, V_BOT, V_TOP, V_RIGHT, V_LEFT, LInfinityMarkedTriangulation
sage: gluings = {(0,0):(1,2), (1,2):(0,0), (0,1):(1,0), (1,0):(0,1),
....: (0,2):(1,1), (1,1):(0,2)}
sage: types = [(V_BOT, V_NONE, V_LEFT), (V_NONE, V_RIGHT, V_TOP)]
sage: T = LInfinityMarkedTriangulation(2, gluings, types)
sage: T.right_left_pairs()
[(1, 1, 0, 2)]
"""
pairs = []
for p1 in range(self._n):
for e1 in range(3):
if self._edge_types[p1][e1] == V_RIGHT:
e1p1 = (e1 + 1) % 3
p2, e2 = self.opposite_edge(p1, e1p1)
e2m1 = (e2 - 1) % 3
if self._edge_types[p2][e2m1] == V_LEFT:
pairs.append((p1, e1, p2, e2m1))
return pairs
[docs]
@cached_method
def polytope(self):
r"""
Each edge correspond to a vector in RR^2 (identified to CC)
We assign the following coordinates
(p,e) -> real part at 2*(3*p + e) and imag part at 2*(3*p + e) + 1
The return polyhedron is compact as we fix each side to be of L-infinity
length less than 1.
"""
dim = 4 * self.num_edges()
eqns = []
ieqs = []
signs = [None] * dim
# edges should sum up to zero
for p in range(self._n):
linear = [0] * (dim + 1)
linear[6 * p + 1] = 1
linear[6 * p + 3] = 1
linear[6 * p + 5] = 1
eqns.append(linear)
linear = [0] * (dim + 1)
linear[6 * p + 2] = 1
linear[6 * p + 4] = 1
linear[6 * p + 6] = 1
eqns.append(linear)
# opposite edges are opposite vectors
for p1 in range(self._n):
for e1 in range(3):
p2, e2 = self.opposite_edge(p1, e1)
re1 = 2 * (3 * p1 + e1)
im1 = 2 * (3 * p1 + e1) + 1
re2 = 2 * (3 * p2 + e2)
im2 = 2 * (3 * p2 + e2) + 1
if re1 < re2:
eqns.append(opposite_condition(dim, re1, re2))
eqns.append(opposite_condition(dim, im1, im2))
# Compute the signs depending on edge types
for p in range(self._n):
for e1, e2 in ((2, 0), (0, 1), (1, 2)):
t = self._edge_types[p][e2]
re1 = 2 * (3 * p + e1)
im1 = 2 * (3 * p + e1) + 1
re2 = 2 * (3 * p + e2)
im2 = 2 * (3 * p + e2) + 1
if t == V_BOT:
signs[re1] = signs[re2] = +1
signs[im1] = -1
signs[im2] = +1
elif t == V_TOP:
signs[re1] = signs[re2] = -1
signs[im1] = +1
signs[im2] = -1
elif t == V_RIGHT:
signs[re1] = +1
signs[re2] = -1
signs[im1] = signs[im2] = +1
elif t == V_LEFT:
signs[re1] = -1
signs[re2] = +1
signs[im1] = signs[im2] = -1
# adding sign conditions
for i in range(dim):
ieqs.extend(sign_and_norm_conditions(dim, i, signs[i]))
# Delaunay conditions
for p1, e1, p2, e2 in self.bottom_top_pairs():
ieqs.append(bottom_top_delaunay_condition(dim, p1, e1, p2, e2))
for p1, e1, p2, e2 in self.right_left_pairs():
ieqs.append(right_left_delaunay_condition(dim, p1, e1, p2, e2))
# return eqns, ieqs
from sage.geometry.polyhedron.constructor import Polyhedron
from sage.rings.rational_field import QQ
return Polyhedron(ieqs=ieqs, eqns=eqns, base_ring=QQ)
[docs]
def barycenter(self):
r"""
Return the translation surface in the barycenter of this polytope.
EXAMPLES::
sage: from flatsurf.geometry.l_infinity_delaunay_cells import \
....: V_NONE, V_BOT, V_TOP, V_RIGHT, V_LEFT, LInfinityMarkedTriangulation
sage: gluings = {(0,0):(1,2), (1,2):(0,0), (0,1):(1,0), (1,0):(0,1),
....: (0,2):(1,1), (1,1):(0,2)}
sage: types = [(V_BOT, V_NONE, V_LEFT), (V_NONE, V_RIGHT, V_TOP)]
sage: T = LInfinityMarkedTriangulation(2, gluings, types)
sage: S = T.barycenter()
sage: S.polygon(0)
Polygon(vertices=[(0, 0), (3/7, 13/21), (-3/7, 11/21)])
sage: S.polygon(1)
Polygon(vertices=[(0, 0), (6/7, 2/21), (3/7, 13/21)])
"""
verts = [v.vector() for v in self.polytope().vertices()]
b = sum(verts) / len(verts)
from flatsurf import Polygon
from sage.rings.rational_field import QQ
from flatsurf import MutableOrientedSimilaritySurface
barycenter = MutableOrientedSimilaritySurface(QQ)
for p in range(self._n):
e1 = (b[6 * p], b[6 * p + 1])
e2 = (b[6 * p + 2], b[6 * p + 3])
e3 = (b[6 * p + 4], b[6 * p + 5])
barycenter.add_polygon(Polygon(edges=[e1, e2, e3], base_ring=QQ))
for gluing in self._edge_identifications.items():
barycenter.glue(*gluing)
return barycenter