r"""
The category of Euclidean polygons defined in the real plane.
See :mod:`flatsurf.geometry.categories` for a general description of the
category framework in sage-flatsurf.
Normally, you won't create this (or any other) category directly. The correct
category of a polygon is automatically determined.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: C = EuclideanPolygons(QQ)
sage: from flatsurf import polygons
sage: polygons.square() in C
True
.. jupyter-execute::
:hide-code:
# Allow jupyter-execute blocks in this module to contain doctests
import jupyter_doctest_tweaks
"""
# ****************************************************************************
# This file is part of sage-flatsurf.
#
# Copyright (C) 2016-2020 Vincent Delecroix
# 2020-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.categories.category_types import Category_over_base_ring
from sage.categories.category_with_axiom import CategoryWithAxiom_over_base_ring
from sage.misc.cachefunc import cached_method
from sage.all import FreeModule
from sage.misc.abstract_method import abstract_method
from sage.structure.element import get_coercion_model
from flatsurf.geometry.categories.polygons import Polygons
from flatsurf.geometry.euclidean import ccw
cm = get_coercion_model()
[docs]class EuclideanPolygons(Category_over_base_ring):
r"""
The category of Euclidean polygons defined in the real plane
over a fixed base ring.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons(QQ)
Category of euclidean polygons over Rational Field
"""
[docs] def super_categories(self):
r"""
Return the categories Euclidean polygons are also contained in, namely
the polygons.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons(QQ).super_categories()
[Category of polygons over Rational Field]
"""
return [Polygons(self.base_ring())]
[docs] class ParentMethods:
r"""
Provides methods available to all Euclidean polygons in the real plane.
If you want to add functionality to such polygons, you probably want to
put it here.
"""
[docs] def vector_space(self):
r"""
Return the vector space of dimension 2 in which this polygon embeds.
EXAMPLES::
sage: from flatsurf import Polygons
sage: C = Polygons(QQ)
sage: C.vector_space()
doctest:warning
...
UserWarning: vector_space() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring().fraction_field()**2 instead
Vector space of dimension 2 over Rational Field
"""
import warnings
warnings.warn(
"vector_space() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring().fraction_field()**2 instead"
)
return self.base_ring().fraction_field() ** 2
[docs] def module(self):
r"""
Return the free module of rank 2 in which this polygon embeds.
EXAMPLES::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S.module()
doctest:warning
...
UserWarning: module() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring()**2 instead
Vector space of dimension 2 over Rational Field
"""
import warnings
warnings.warn(
"module() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring()**2 instead"
)
return self.base_ring() ** 2
[docs] def field(self):
r"""
EXAMPLES::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S.field()
doctest:warning
...
UserWarning: field() has been deprecated and will be removed from a future version of sage-flatsurf; use base_ring() or base_ring().fraction_field() instead
Rational Field
"""
import warnings
warnings.warn(
"field() has been deprecated and will be removed from a future version of sage-flatsurf; use base_ring() or base_ring().fraction_field() instead"
)
return self.base_ring().fraction_field()
def _mul_(self, g, switch_sides=None):
r"""
Apply the 2x2 matrix `g` to this polygon.
The matrix must have non-zero determinant. If the determinant is
negative, then the vertices and edges are relabeled according to the
involutions `v \mapsto (n-v)%n` and `e \mapsto n-1-e` respectively.
EXAMPLES::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices = [(1,0),(0,1),(-1,-1)])
sage: p
Polygon(vertices=[(1, 0), (0, 1), (-1, -1)])
sage: matrix(ZZ,[[0, 1], [1, 0]]) * p
Polygon(vertices=[(0, 1), (-1, -1), (1, 0)])
sage: matrix(ZZ,[[2, 0], [0, 1]]) * p
Polygon(vertices=[(2, 0), (0, 1), (-2, -1)])
"""
from flatsurf import Polygon
if g in self.base_ring():
from sage.all import MatrixSpace
g = MatrixSpace(self.base_ring(), 2)(g)
det = g.det()
if det == 0:
raise ValueError(
"Can not act on a polygon with matrix with zero determinant"
)
if det < 0:
# Note that in this case we reverse the order
vertices = [g * self.vertex(0)]
for i in range(len(self.vertices()) - 1, 0, -1):
vertices.append(g * self.vertex(i))
return Polygon(vertices=vertices, check=False)
return Polygon(
vertices=[g * v for v in self.vertices()],
check=False,
category=self.category(),
)
[docs] @cached_method
def is_rational(self):
r"""
Return whether this is a rational polygon, i.e., all its
:meth:`angles` are rational multiples of π.
EXAMPLES::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices = [(1, 0), (0, 1), (-1, -1)])
sage: p.is_rational()
False
Note that determining rationality is somewhat costly. Once
established, this refines the category of the triangle::
sage: p = Polygon(vertices = [(0, 0), (1, 0), (0, 1)])
sage: p.category()
Category of convex simple euclidean polygons over Rational Field
sage: p.is_rational()
True
sage: p.category()
Category of rational convex simple euclidean polygons over Rational Field
"""
for e in range(len(self.vertices())):
u = self.edge(e)
v = -self.edge((e - 1) % len(self.vertices()))
cos = u.dot_product(v)
sin = u[0] * v[1] - u[1] * v[0]
from flatsurf.geometry.euclidean import is_cosine_sine_of_rational
if not is_cosine_sine_of_rational(cos, sin, scaled=True):
return False
self._refine_category_(self.category().Rational())
return True
[docs] def is_simple(self):
r"""
Return whether this is a simple polygon, i.e., without
self-intersection.
EXAMPLES:
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.is_simple()
True
"""
n = len(self.vertices())
for i in range(n):
ei = (self.vertex(i), self.vertex(i + 1))
for j in range(i + 2, n + 1):
if (i - j) % n in [-1, 0, 1]:
continue
ej = (self.vertex(j), self.vertex(j + 1))
from flatsurf.geometry.euclidean import is_segment_intersecting
if is_segment_intersecting(ei, ej):
return False
return True
[docs] @abstract_method
def vertices(self, marked_vertices=True):
r"""
Return the vertices of this polygon in counterclockwise order as
vectors in the real plane.
INPUT:
- ``marked_vertices`` -- a boolean (default: ``True``); whether to
include marked vertices that are not actually corners of the
polygon.
EXAMPLES::
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.vertices()
((0, 0), (1, 0), (1, 1), (0, 1))
"""
[docs] def vertex(self, i):
r"""
Return coordinates for the ``i``-th vertex of this polygon.
EXAMPLES:
The ``i`` wraps around if it is negative or exceeds the number of
vertices in this polygon::
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.vertex(-1)
(0, 1)
sage: s.vertex(0)
(0, 0)
sage: s.vertex(1)
(1, 0)
sage: s.vertex(2)
(1, 1)
sage: s.vertex(3)
(0, 1)
sage: s.vertex(4)
(0, 0)
"""
vertices = self.vertices()
return vertices[i % len(vertices)]
[docs] def edges(self):
r"""
Return the edges of this polygon as vectors in the plane going from
one vertex to the next one.
EXAMPLES::
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.edges()
[(1, 0), (0, 1), (-1, 0), (0, -1)]
"""
return [self.edge(i) for i in range(len(self.vertices()))]
[docs] def edge(self, i):
r"""
Return the vector going from vertex ``i`` to the following vertex
in counter-clockwise order.
EXAMPLES:
Note that this wraps around if ``i`` is negative or exceeds the
number of vertices::
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.edge(-1)
(0, -1)
sage: s.edge(0)
(1, 0)
sage: s.edge(1)
(0, 1)
sage: s.edge(2)
(-1, 0)
sage: s.edge(3)
(0, -1)
sage: s.edge(4)
(1, 0)
"""
return self.vertex(i + 1) - self.vertex(i)
[docs] def is_convex(self, strict=False):
r"""
Return whether this is a convex polygon.
INPUT:
- ``strict`` -- whether to check for strict convexity, i.e., a
polygon with a π angle is not considered convex.
EXAMPLES::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S.is_convex()
True
sage: S.is_convex(strict=True)
True
"""
from flatsurf.geometry.euclidean import ccw
for i in range(len(self.vertices())):
consecutive_ccw = ccw(self.edge(i), self.edge(i + 1))
if strict:
if consecutive_ccw <= 0:
return False
else:
if consecutive_ccw < 0:
return False
return True
def _test_marked_vertices(self, **options):
r"""
Verify that :meth:`vertices` and :meth:`is_convex` are compatible.
EXAMPLES::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S._test_marked_vertices()
"""
tester = self._tester(**options)
if self.is_convex():
tester.assertEqual(
self.is_convex(strict=True),
self.vertices() == self.vertices(marked_vertices=False),
)
[docs] def is_degenerate(self):
r"""
Return whether this polygon is considered degenerate.
This implements
:meth:`flatsurf.geometry.categories.polygons.Polygons.ParentMethods.is_degenerate`.
EXAMPLES:
Polygons with zero area are considered degenerate::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (1, 0)], check=False)
sage: p.is_degenerate()
True
Polygons with marked vertices are considered degenerate::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (4, 0), (2, 2)])
sage: p.is_degenerate()
True
"""
if self.area() == 0:
return True
if self.vertices() != self.vertices(marked_vertices=False):
return True
return False
[docs] def slopes(self, relative=False):
r"""
Return the slopes of this polygon as vectors in the plane.
INPUT:
- ``relative`` -- a boolean (default: ``False``); whether to return the
slopes not as absolute vectors parallel to the edges but relative to
the previous edge, i.e., after turning the previous edge to be
parallel to the x axis.
EXAMPLES::
sage: from flatsurf import polygons
sage: s = polygons.square()
sage: s.slopes()
[(1, 0), (0, 1), (-1, 0), (0, -1)]
sage: s.slopes(relative=True)
[(0, 1), (0, 1), (0, 1), (0, 1)]
A polygon with a marked point::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (4, 0), (2, 2)])
sage: p.slopes()
[(2, 0), (2, 0), (-2, 2), (-2, -2)]
sage: p.slopes(relative=True)
[(-4, 4), (4, 0), (-4, 4), (0, 8)]
"""
if not relative:
return self.edges()
edges = [
(self.edge((e - 1) % len(self.vertices())), self.edge(e))
for e in range(len(self.vertices()))
]
cos = [u.dot_product(v) for (u, v) in edges]
sin = [u[0] * v[1] - u[1] * v[0] for (u, v) in edges]
from sage.all import vector
return [vector((c, s)) for (c, s) in zip(cos, sin)]
[docs] def erase_marked_vertices(self):
r"""
Return a copy of this polygon without marked vertices.
EXAMPLES::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (4, 0), (2, 2)])
sage: p.erase_marked_vertices()
Polygon(vertices=[(0, 0), (4, 0), (2, 2)])
"""
from flatsurf import Polygon
return Polygon(vertices=self.vertices(marked_vertices=False))
[docs] def is_equilateral(self):
r"""
Return whether all sides of this polygon have the same length.
EXAMPLES::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (2, 2), (0, 2)])
sage: p.is_equilateral()
True
"""
return len({edge[0] ** 2 + edge[1] ** 2 for edge in self.edges()}) == 1
[docs] def is_equiangular(self):
r"""
Return whether all sides of this polygon meet at the same angle.
EXAMPLES::
sage: from flatsurf import Polygon
sage: p = Polygon(vertices=[(0, 0), (2, 0), (2, 2), (0, 2)])
sage: p.is_equiangular()
True
"""
slopes = self.slopes(relative=True)
from flatsurf.geometry.euclidean import is_parallel
return all(
is_parallel(slopes[i - 1], slopes[i]) for i in range(len(slopes))
)
[docs] def plot(
self,
translation=None,
polygon_options={},
edge_options={},
vertex_options={},
):
r"""
Return a plot of this polygon with the origin at ``translation``.
EXAMPLES:
.. jupyter-execute::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S.plot()
...Graphics object consisting of 3 graphics primitives
We can specify an explicit ``zorder`` to render edges and vertices on
top of the axes which are rendered at z-order 3:
.. jupyter-execute::
sage: S.plot(edge_options={'zorder': 3}, vertex_options={'zorder': 3})
...Graphics object consisting of 3 graphics primitives
We can control the colors, e.g., we can render transparent polygons,
with red edges and blue vertices:
.. jupyter-execute::
sage: S.plot(polygon_options={'fill': None}, edge_options={'color': 'red'}, vertex_options={'color': 'blue'})
...Graphics object consisting of 3 graphics primitives
"""
from sage.plot.point import point2d
from sage.plot.line import line2d
from sage.plot.polygon import polygon2d
P = self.vertices(translation)
polygon_options = {"alpha": 0.3, "zorder": 1, **polygon_options}
edge_options = {"color": "orange", "zorder": 2, **edge_options}
vertex_options = {"color": "red", "zorder": 2, **vertex_options}
return (
polygon2d(P, **polygon_options)
+ line2d(P + (P[0],), **edge_options)
+ point2d(P, **vertex_options)
)
[docs] def angles(self, numerical=None, assume_rational=None):
r"""
Return the list of angles of this polygon (divided by `2 \pi`).
EXAMPLES::
sage: from flatsurf import Polygon
sage: T = Polygon(angles=[1, 2, 3])
sage: [T.angle(i) for i in range(3)]
[1/12, 1/6, 1/4]
sage: T.angles()
(1/12, 1/6, 1/4)
sage: sum(T.angle(i) for i in range(3))
1/2
"""
if assume_rational is not None:
import warnings
warnings.warn(
"assume_rational has been deprecated as a keyword to angles() and will be removed from a future version of sage-flatsurf"
)
angles = tuple(
self.angle(i, numerical=numerical) for i in range(len(self.vertices()))
)
if not numerical:
self._refine_category_(self.category().WithAngles(angles))
return angles
[docs] def angle(self, e, numerical=None, assume_rational=None):
r"""
Return the angle at the beginning of the start point of the edge ``e``.
EXAMPLES::
sage: from flatsurf.geometry.polygon import polygons
sage: polygons.square().angle(0)
1/4
sage: polygons.regular_ngon(8).angle(0)
3/8
sage: from flatsurf import Polygon
sage: T = Polygon(vertices=[(0,0), (3,1), (1,5)])
sage: [T.angle(i, numerical=True) for i in range(3)] # abs tol 1e-13
[0.16737532973071603, 0.22741638234956674, 0.10520828791971722]
sage: sum(T.angle(i, numerical=True) for i in range(3)) # abs tol 1e-13
0.5
"""
if assume_rational is not None:
import warnings
warnings.warn(
"assume_rational has been deprecated as a keyword to angle() and will be removed from a future version of sage-flatsurf"
)
if numerical is None:
numerical = not self.is_rational()
if numerical:
import warnings
warnings.warn(
"the behavior of angle() has been changed in recent versions of sage-flatsurf; for non-rational polygons, numerical=True must be set explicitly to get a numerical approximation of the angle"
)
from flatsurf.geometry.euclidean import angle
return angle(
self.edge(e),
-self.edge((e - 1) % len(self.vertices())),
numerical=numerical,
)
[docs] def area(self):
r"""
Return the area of this polygon.
EXAMPLES::
sage: from flatsurf.geometry.polygon import polygons
sage: polygons.regular_ngon(8).area()
2*a + 2
sage: _ == 2*AA(2).sqrt() + 2
True
sage: AA(polygons.regular_ngon(11).area())
9.36563990694544?
sage: polygons.square().area()
1
sage: (2*polygons.square()).area()
4
"""
# Will use an area formula obtainable from Green's theorem. See for instance:
# http://math.blogoverflow.com/2014/06/04/greens-theorem-and-area-of-polygons/
total = self.base_ring().zero()
for i in range(len(self.vertices())):
total += (self.vertex(i)[0] + self.vertex(i + 1)[0]) * self.edge(i)[1]
from sage.all import ZZ
return total / ZZ(2)
[docs] def centroid(self):
r"""
Return the coordinates of the centroid of this polygon.
ALGORITHM:
We use the customary formula of the centroid of polygons, see
https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
EXAMPLES::
sage: from flatsurf.geometry.polygon import polygons
sage: P = polygons.regular_ngon(4)
sage: P
Polygon(vertices=[(0, 0), (1, 0), (1, 1), (0, 1)])
sage: P.centroid()
(1/2, 1/2)
sage: P = polygons.regular_ngon(8); P
Polygon(vertices=[(0, 0), (1, 0), (1/2*a + 1, 1/2*a), (1/2*a + 1, 1/2*a + 1), (1, a + 1), (0, a + 1), (-1/2*a, 1/2*a + 1), (-1/2*a, 1/2*a)])
sage: P.centroid()
(1/2, 1/2*a + 1/2)
sage: P = polygons.regular_ngon(11)
sage: C = P.centroid()
sage: P = P.translate(-C)
sage: P.centroid()
(0, 0)
"""
x, y = list(zip(*self.vertices()))
nvertices = len(x)
A = self.area()
from sage.all import vector
return vector(
(
~(6 * A)
* sum(
[
(x[i - 1] + x[i]) * (x[i - 1] * y[i] - x[i] * y[i - 1])
for i in range(nvertices)
]
),
~(6 * A)
* sum(
[
(y[i - 1] + y[i]) * (x[i - 1] * y[i] - x[i] * y[i - 1])
for i in range(nvertices)
]
),
)
)
[docs] def get_point_position(self, point, translation=None):
r"""
Return the combinatorial classification of a point and a polygon.
INPUT:
- ``point`` -- a point in the plane as a SageMath vector or pair of
numbers
OUTPUT:
a :class:`.geometry.polygon.PolygonPosition` object
ALGORITHM:
We use a winding number algorithm see
https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm
EXAMPLES::
sage: from flatsurf import polygons, Polygon
sage: S = polygons.square()
sage: S.get_point_position((1/2, 1/2))
point positioned in interior of polygon
sage: S.get_point_position((1, 0))
point positioned on vertex 1 of polygon
sage: S.get_point_position((1, 1/2))
point positioned on interior of edge 1 of polygon
sage: S.get_point_position((1, 3/2))
point positioned outside polygon
sage: p = Polygon(edges=[(1, 0), (1, 0), (1, 0), (0, 1), (-3, 0), (0, -1)])
sage: p.get_point_position([10, 0])
point positioned outside polygon
sage: p.get_point_position([1/2, 0])
point positioned on interior of edge 0 of polygon
sage: p.get_point_position([3/2, 0])
point positioned on interior of edge 1 of polygon
sage: p.get_point_position([2,0])
point positioned on vertex 2 of polygon
sage: p.get_point_position([5/2, 0])
point positioned on interior of edge 2 of polygon
sage: p.get_point_position([5/2, 1/4])
point positioned in interior of polygon
"""
from sage.all import vector
point = vector(point)
if translation is not None:
import warnings
warnings.warn(
"the translation keyword argument to get_point_position() has been deprecated and will be removed in a future version of sage-flatsurf; shift the point instead with the - operator"
)
from sage.all import vector
return self.get_point_position(point - vector(translation))
from flatsurf.geometry.euclidean import ccw
from flatsurf.geometry.polygon import PolygonPosition
# Determine whether the point is a vertex of the polygon.
for i, v in enumerate(self.vertices()):
if point == v:
return PolygonPosition(PolygonPosition.VERTEX, vertex=i)
# Determine whether the point is on an edge of the polygon.
for i, (v, e) in enumerate(zip(self.vertices(), self.edges())):
if ccw(e, point - v) == 0:
# The point lies on the line through this edge.
if 0 < e.dot_product(point - v) < e.dot_product(e):
return PolygonPosition(PolygonPosition.EDGE_INTERIOR, edge=i)
# Determine whether the point is inside or outside by computing the
# winding number of the polygon.
winding_number = 0
for v, w in zip(self.vertices(), self.vertices()[1:] + self.vertices()[:1]):
if v[1] < point[1] and w[1] >= point[1] and ccw(w - v, point - v) > 0:
winding_number += 1
if v[1] >= point[1] and w[1] < point[1] and ccw(w - v, point - v) < 0:
winding_number -= 1
if winding_number % 2:
return PolygonPosition(PolygonPosition.INTERIOR)
return PolygonPosition(PolygonPosition.OUTSIDE)
[docs] class Rational(CategoryWithAxiom_over_base_ring):
r"""
The category of rational Euclidean polygons.
.. NOTE::
This category must be defined here to make SageMath's test suite
pass. Otherwise we get "The super categories of a category over
base should be a category over base (or the related Bimodules) or a
singleton category"; we did not investigate what exactly is going on
here.
"""
[docs] class SubcategoryMethods:
[docs] @cached_method
def module(self):
r"""
Return the free module of rank 2 in which these polygons embed.
EXAMPLES::
sage: from flatsurf import Polygons
sage: C = Polygons(QQ)
sage: C.module()
doctest:warning
...
UserWarning: module() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring()**2 instead
Vector space of dimension 2 over Rational Field
::
sage: from flatsurf import EuclideanPolygonsWithAngles
sage: C = EuclideanPolygonsWithAngles(1, 2, 3)
sage: C.module()
Vector space of dimension 2 over Number Field in c with defining polynomial x^2 - 3 with c = 1.732050807568878?
"""
import warnings
warnings.warn(
"module() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring()**2 instead"
)
return FreeModule(self.base_ring(), 2)
[docs] @cached_method
def vector_space(self):
r"""
Return the vector space of dimension 2 in which these polygons embed.
EXAMPLES::
sage: from flatsurf import Polygons
sage: C = Polygons(QQ)
sage: C.vector_space()
Vector space of dimension 2 over Rational Field
::
sage: from flatsurf import EuclideanPolygonsWithAngles
sage: C = EuclideanPolygonsWithAngles(1, 2, 3)
sage: C.vector_space()
doctest:warning
...
UserWarning: vector_space() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring().fraction_field()**2 instead
Vector space of dimension 2 over Number Field in c with defining polynomial x^2 - 3 with c = 1.732050807568878?
"""
import warnings
warnings.warn(
"vector_space() has been deprecated and will be removed in a future version of sage-flatsurf; use base_ring().fraction_field()**2 instead"
)
from sage.all import VectorSpace
return VectorSpace(self.base_ring().fraction_field(), 2)
[docs] def WithAngles(self, angles):
r"""
Return the subcategory of polygons with fixed ``angles``.
INPUT:
- ``angles`` -- a finite sequence of numbers, the inner angles of
the polygon; the angles are automatically normalized to sum to
`(n-2)π`.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons(AA).WithAngles([1, 2, 3])
Category of euclidean triangles with angles (1/12, 1/6, 1/4) over Algebraic Real Field
"""
from flatsurf.geometry.categories.euclidean_polygons_with_angles import (
EuclideanPolygonsWithAngles,
)
angles = EuclideanPolygonsWithAngles._normalize_angles(angles)
return EuclideanPolygonsWithAngles(self.base_ring(), angles) & self
def __call__(self, *args, **kwds):
r"""
TESTS::
sage: from flatsurf import Polygons, ConvexPolygons
sage: C = Polygons(QQ)
sage: p = C(vertices=[(0,0),(1,0),(2,0),(1,1)])
doctest:warning
...
UserWarning: Polygons(…)(…) has been deprecated and will be removed in a future version of sage-flatsurf; use Polygon() instead
sage: p
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: C(p) is p
False
sage: C(p) == p
True
sage: C((1,0), (0,1), (-1, 1))
Traceback (most recent call last):
...
ValueError: the polygon does not close up
sage: D = ConvexPolygons(QQbar)
doctest:warning
...
UserWarning: ConvexPolygons() has been deprecated and will be removed from a future version of sage-flatsurf; use Polygon() to create polygons.
If you really need the category of convex polygons over a ring use EuclideanPolygons(ring).Simple().Convex() instead.
sage: D(p)
doctest:warning
...
UserWarning: ConvexPolygons(…)(…) has been deprecated and will be removed in a future version of sage-flatsurf; use Polygon() instead
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: D(vertices=p.vertices())
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: D(edges=p.edges())
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
"""
# We cannot have a __call__() in SubcategoryMethods so there is no good
# way to support this in the category framework. Also, this code is
# duplicated in several places and the Polygon() helper seems to be
# much more versatile.
import warnings
warnings.warn(
"Polygons(…)(…) has been deprecated and will be removed in a future version of sage-flatsurf; use Polygon() instead"
)
check = kwds.pop("check", True)
from flatsurf.geometry.polygon import EuclideanPolygon
if len(args) == 1 and isinstance(args[0], EuclideanPolygon):
if args[0].category() is self:
return args[0]
vertices = [self.vector_space()(v) for v in args[0].vertices()]
args = ()
else:
vertices = kwds.pop("vertices", None)
edges = kwds.pop("edges", None)
base_point = kwds.pop("base_point", (0, 0))
if (vertices is None) and (edges is None):
if len(args) == 1:
edges = args[0]
elif args:
edges = args
else:
raise ValueError(
"exactly one of 'vertices' or 'edges' must be provided"
)
if kwds:
raise ValueError("invalid keyword {!r}".format(next(iter(kwds))))
if edges is not None:
v = self.vector_space()(base_point)
vertices = []
for e in map(self.vector_space(), edges):
vertices.append(v)
v += e
if v != vertices[0]:
raise ValueError("the polygon does not close up")
from flatsurf.geometry.polygon import Polygon
return Polygon(
base_ring=self.base(), vertices=vertices, category=self, check=check
)
[docs] class Convex(CategoryWithAxiom_over_base_ring):
r"""
The subcategory of convex Euclidean polygons in the real plane.
EXAMPLES:
For historic reasons, there is the shortcut ``ConvexPolygons`` to get
the Euclidean convex polygons::
sage: from flatsurf import ConvexPolygons
sage: C = ConvexPolygons(QQ)
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: C is EuclideanPolygons(QQ).Convex().Simple()
True
sage: C(vertices=[(0,0), (2,0), (1,1)])
Polygon(vertices=[(0, 0), (2, 0), (1, 1)])
sage: C(edges=[(1,0), (0,1), (-1,0), (0,-1)])
Polygon(vertices=[(0, 0), (1, 0), (1, 1), (0, 1)])
This axiom can also be created over non-fields::
sage: ConvexPolygons(ZZ)
Category of convex simple euclidean polygons over Integer Ring
TESTS::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: TestSuite(EuclideanPolygons(QQ).Convex()).run()
sage: TestSuite(EuclideanPolygons(QQbar).Convex()).run()
sage: TestSuite(EuclideanPolygons(ZZ).Convex()).run()
"""
[docs] class ParentMethods:
r"""
Provides methods available to all convex Euclidean polygons in the
real plane.
If you want to add functionality to such polygons, you probably
want to put it here.
"""
[docs] def is_convex(self, strict=False):
r"""
Return whether this is a convex polygon.
INPUT:
- ``strict`` -- whether to check for strict convexity, i.e., a
polygon with a π angle is not considered convex.
EXAMPLES::
sage: from flatsurf import polygons
sage: S = polygons.square()
sage: S.is_convex()
True
sage: S.is_convex(strict=True)
True
"""
if not strict:
return True
return EuclideanPolygons.ParentMethods.is_convex(self, strict=strict)
[docs] class Simple(CategoryWithAxiom_over_base_ring):
r"""
The subcategory of Euclidean polygons without self-intersection in the
real plane.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons(QQ).Simple()
Category of simple euclidean polygons over Rational Field
"""
[docs] class ParentMethods:
r"""
Provides methods available to all simple Euclidean polygons.
If you want to add functionality to all polygons, independent of
implementation, you probably want to put it here.
"""
[docs] def triangulation(self):
r"""
Return a list of pairs of indices of vertices that together with the boundary
form a triangulation.
EXAMPLES:
We triangulate a non-convex polygon::
sage: from flatsurf import Polygon
sage: P = Polygon(vertices=[(0,0), (1,0), (1,1), (0,1), (0,2), (-1,2), (-1,1), (-2,1),
....: (-2,0), (-1,0), (-1,-1), (0,-1)])
sage: P.triangulation()
[(0, 2), (2, 8), (3, 5), (6, 8), (8, 3), (3, 6), (9, 11), (0, 9), (2, 9)]
TESTS::
sage: Polygon(vertices=[(0,0), (1,0), (1,1), (0,1)]).triangulation()
[(0, 2)]
sage: quad = [(0,0), (1,-1), (0,1), (-1,-1)]
sage: for i in range(4):
....: Polygon(vertices=quad[i:] + quad[:i]).triangulation()
[(0, 2)]
[(1, 3)]
[(0, 2)]
[(1, 3)]
sage: poly = [(0,0),(1,1),(2,0),(3,1),(4,0),(4,2),
....: (-4,2),(-4,0),(-3,1),(-2,0),(-1,1)]
sage: Polygon(vertices=poly).triangulation()
[(1, 3), (3, 5), (5, 8), (6, 8), (8, 10), (10, 1), (1, 5), (5, 10)]
sage: for i in range(len(poly)):
....: Polygon(vertices=poly[i:] + poly[:i]).triangulation()
[(1, 3), (3, 5), (5, 8), (6, 8), (8, 10), (10, 1), (1, 5), (5, 10)]
[(0, 2), (2, 4), (4, 7), (5, 7), (7, 9), (9, 0), (0, 4), (4, 9)]
[(1, 3), (3, 6), (4, 6), (6, 8), (8, 10), (10, 1), (3, 8), (10, 3)]
[(0, 2), (2, 5), (3, 5), (5, 7), (7, 9), (9, 0), (2, 7), (9, 2)]
[(1, 4), (2, 4), (4, 6), (6, 8), (8, 10), (10, 1), (1, 6), (8, 1)]
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 9), (9, 0), (0, 5), (7, 0)]
[(0, 2), (2, 4), (4, 6), (6, 8), (8, 10), (10, 2), (4, 10), (6, 10)]
[(1, 3), (3, 5), (5, 7), (7, 9), (9, 1), (10, 1), (3, 9), (5, 9)]
[(0, 2), (2, 4), (4, 6), (6, 8), (8, 0), (9, 0), (2, 8), (4, 8)]
[(1, 3), (3, 5), (5, 7), (7, 10), (8, 10), (10, 1), (1, 7), (3, 7)]
[(0, 2), (2, 4), (4, 6), (6, 9), (7, 9), (9, 0), (0, 6), (2, 6)]
sage: poly = [(0,0), (1,0), (2,0), (2,1), (2,2), (1,2), (0,2), (0,1)]
sage: Polygon(vertices=poly).triangulation()
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
sage: for i in range(len(poly)):
....: Polygon(vertices=poly[i:] + poly[:i]).triangulation()
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
sage: poly = [(0,0), (1,2), (3,3), (1,4), (0,6), (-1,4), (-3,-3), (-1,2)]
sage: Polygon(vertices=poly).triangulation()
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
sage: for i in range(len(poly)):
....: Polygon(vertices=poly[i:] + poly[:i]).triangulation()
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 2), (3, 5), (5, 7), (7, 3), (0, 3)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
[(0, 6), (1, 3), (3, 5), (5, 1), (6, 1)]
[(0, 2), (2, 4), (4, 6), (6, 0), (0, 4)]
sage: x = polygen(QQ)
sage: p = x^4 - 5*x^2 + 5
sage: r = AA.polynomial_root(p, RIF(1.17,1.18))
sage: K.<a> = NumberField(p, embedding=r)
sage: poly = [(1/2*a^2 - 3/2, 1/2*a),
....: (-a^3 + 2*a^2 + 2*a - 4, 0),
....: (1/2*a^2 - 3/2, -1/2*a),
....: (1/2*a^3 - a^2 - 1/2*a + 1, 1/2*a^2 - a),
....: (-1/2*a^2 + 1, 1/2*a^3 - 3/2*a),
....: (-1/2*a + 1, a^3 - 3/2*a^2 - 2*a + 5/2),
....: (1, 0),
....: (-1/2*a + 1, -a^3 + 3/2*a^2 + 2*a - 5/2),
....: (-1/2*a^2 + 1, -1/2*a^3 + 3/2*a),
....: (1/2*a^3 - a^2 - 1/2*a + 1, -1/2*a^2 + a)]
sage: Polygon(vertices=poly).triangulation()
[(0, 3), (1, 3), (3, 5), (5, 7), (7, 9), (9, 3), (3, 7)]
sage: z = QQbar.zeta(24)
sage: pts = [(1+i%2) * z**i for i in range(24)]
sage: pts = [vector(AA, (x.real(), x.imag())) for x in pts]
sage: Polygon(vertices=pts).triangulation()
[(0, 2), ..., (16, 0)]
This is https://github.com/flatsurf/sage-flatsurf/issues/87 ::
sage: x = polygen(QQ)
sage: K.<c> = NumberField(x^2 - 3, embedding=AA(3).sqrt())
sage: Polygon(vertices=[(0, 0), (1, 0), (1/2*c + 1, -1/2), (c + 1, 0), (-3/2*c + 1, 5/2), (0, c - 2)]).triangulation()
[(0, 4), (1, 3), (4, 1)]
"""
vertices = self.vertices()
n = len(vertices)
if n < 3:
raise ValueError
if n == 3:
return []
# NOTE: The algorithm is naive. We look at all possible chords between
# the i-th and j-th vertices. If the chord does not intersect any edge
# then we cut the polygon along this edge and call recursively
# triangulate on the two pieces.
for i in range(n - 1):
eiright = vertices[(i + 1) % n] - vertices[i]
eileft = vertices[(i - 1) % n] - vertices[i]
for j in range(i + 2, (n if i else n - 1)):
ejright = vertices[(j + 1) % n] - vertices[j]
ejleft = vertices[(j - 1) % n] - vertices[j]
chord = vertices[j] - vertices[i]
from flatsurf.geometry.euclidean import is_between
# check angles with neighbouring edges
if not (
is_between(eiright, eileft, chord)
and is_between(ejright, ejleft, -chord)
):
continue
# check intersection with other edges
e = (vertices[i], vertices[j])
good = True
for k in range(n):
if k == i or k == j or k == (i - 1) % n or k == (j - 1) % n:
continue
f = (vertices[k], vertices[(k + 1) % n])
from flatsurf.geometry.euclidean import (
is_segment_intersecting,
)
res = is_segment_intersecting(e, f)
assert res != 1
if res == 2:
good = False
break
if good:
from flatsurf import Polygon
part0 = [
(s + i, t + i)
for s, t in Polygon(
vertices=vertices[i : j + 1], check=False
).triangulation()
]
part1 = []
for s, t in Polygon(
vertices=vertices[j:] + vertices[: i + 1], check=False
).triangulation():
if s < n - j:
s += j
else:
s -= n - j
if t < n - j:
t += j
else:
t -= n - j
part1.append((s, t))
return [(i, j)] + part0 + part1
assert False
[docs] class Convex(CategoryWithAxiom_over_base_ring):
r"""
The subcategory of the simple convex Euclidean polygons.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons(QQ).Simple().Convex()
Category of convex simple euclidean polygons over Rational Field
"""
def __call__(self, *args, **kwds):
r"""
TESTS::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: C = EuclideanPolygons(QQ).Convex().Simple()
sage: p = C(vertices=[(0,0),(1,0),(2,0),(1,1)])
doctest:warning
...
UserWarning: ConvexPolygons(…)(…) has been deprecated and will be removed in a future version of sage-flatsurf; use Polygon() instead
sage: p
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: C(p) is p
True
sage: C((1,0), (0,1), (-1, 1))
Traceback (most recent call last):
...
ValueError: the polygon does not close up
sage: D = EuclideanPolygons(QQbar).Convex().Simple()
sage: D(p)
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: D(vertices=p.vertices())
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
sage: D(edges=p.edges())
Polygon(vertices=[(0, 0), (1, 0), (2, 0), (1, 1)])
"""
# We cannot have a __call__() in SubcategoryMethods so there is no good
# way to support this in the category framework. Also, this code is
# duplicated in several places and the Polygon() helper seems to be
# much more versatile.
import warnings
warnings.warn(
"ConvexPolygons(…)(…) has been deprecated and will be removed in a future version of sage-flatsurf; use Polygon() instead"
)
check = kwds.pop("check", True)
from flatsurf.geometry.polygon import EuclideanPolygon
if len(args) == 1 and isinstance(args[0], EuclideanPolygon):
if args[0].category() is self:
return args[0]
vertices = [self.vector_space()(v) for v in args[0].vertices()]
args = ()
else:
vertices = kwds.pop("vertices", None)
edges = kwds.pop("edges", None)
base_point = kwds.pop("base_point", (0, 0))
if (vertices is None) and (edges is None):
if len(args) == 1:
edges = args[0]
elif args:
edges = args
else:
raise ValueError(
"exactly one of 'vertices' or 'edges' must be provided"
)
if kwds:
raise ValueError(
"invalid keyword {!r}".format(next(iter(kwds)))
)
if edges is not None:
v = (self.base_ring() ** 2)(base_point)
vertices = []
for e in map(self.base_ring() ** 2, edges):
vertices.append(v)
v += e
if v != vertices[0]:
raise ValueError("the polygon does not close up")
from flatsurf.geometry.polygon import Polygon
return Polygon(
base_ring=self.base(), vertices=vertices, category=self, check=check
)
[docs] class ParentMethods:
r"""
Provides methods available to all simple convex Euclidean polygons.
If you want to add functionality to all polygons, independent of
implementation, you probably want to put it here.
"""
[docs] def find_separatrix(self, direction=None, start_vertex=0):
r"""
Returns a pair (v,same) where v is a vertex and same is a boolean.
The provided parameter "direction" should be a non-zero vector with
two entries, or by default direction=(0,1).
A separatrix is a ray leaving a vertex and entering the polygon.
The vertex v will have a separatrix leaving it which is parallel to
direction. The returned value "same" answers the question if this separatrix
points in the same direction as "direction". There is a boundary case:
we allow the separatrix to be an edge if and only if traveling along
the sepatrix from the vertex would travel in a counter-clockwise
direction about the polygon.
The vertex returned is uniquely defined from the above if the polygon
is a triangle. Otherwise, we return the first vertex with this property
obtained by inspecting starting at start_vertex (defaults to 0) and
then moving in the counter-clockwise direction.
EXAMPLES::
sage: from flatsurf import polygons
sage: p=polygons.square()
sage: print(p.find_separatrix())
(1, True)
sage: print(p.find_separatrix(start_vertex=2))
(3, False)
"""
if direction is None:
direction = (self.base_ring() ** 2)(
(self.base_ring().zero(), self.base_ring().one())
)
else:
assert not direction.is_zero()
v = start_vertex
n = len(self.vertices())
for i in range(len(self.vertices())):
if (
ccw(self.edge(v), direction) >= 0
and ccw(self.edge(v + n - 1), direction) > 0
):
return v, True
if (
ccw(self.edge(v), direction) <= 0
and ccw(self.edge(v + n - 1), direction) < 0
):
return v, False
v = v + 1 % n
raise RuntimeError("Failed to find a separatrix")
[docs] def contains_point(self, point, translation=None):
r"""
Return whether the point is within the polygon (after the polygon is possibly translated)
"""
return self.get_point_position(
point, translation=translation
).is_inside()
[docs] def flow_to_exit(self, point, direction):
r"""
Flow a point in the direction of holonomy until the point leaves the
polygon. Note that ValueErrors may be thrown if the point is not in the
polygon, or if it is on the boundary and the holonomy does not point
into the polygon.
INPUT:
- ``point`` -- a point in the closure of the polygon (as a vector)
- ``holonomy`` -- direction of motion (a vector of non-zero length)
OUTPUT:
- The point in the boundary of the polygon where the trajectory exits
- a PolygonPosition object representing the combinatorial position of the stopping point
"""
from flatsurf.geometry.polygon import PolygonPosition
V = self.base_ring().fraction_field() ** 2
if direction == V.zero():
raise ValueError("Zero vector provided as direction.")
v0 = self.vertex(0)
for i in range(len(self.vertices())):
e = self.edge(i)
from sage.all import matrix
m = matrix([[e[0], -direction[0]], [e[1], -direction[1]]])
try:
ret = m.inverse() * (point - v0)
s = ret[0]
t = ret[1]
# What if the matrix is non-invertible?
# Answer: You'll get a ZeroDivisionError which means that the edge is parallel
# to the direction.
# s is location it intersects on edge, t is the portion of the direction to reach this intersection
if t > 0 and 0 <= s and s <= 1:
# The ray passes through edge i.
if s == 1:
# exits through vertex i+1
v0 = v0 + e
return v0, PolygonPosition(
PolygonPosition.VERTEX,
vertex=(i + 1) % len(self.vertices()),
)
if s == 0:
# exits through vertex i
return v0, PolygonPosition(
PolygonPosition.VERTEX, vertex=i
)
# exits through vertex i
# exits through interior of edge i
prod = t * direction
return point + prod, PolygonPosition(
PolygonPosition.EDGE_INTERIOR, edge=i
)
except ZeroDivisionError:
# Here we know the edge and the direction are parallel
if ccw(e, point - v0) == 0:
# In this case point lies on the edge.
# We need to work out which direction to move in.
from flatsurf.geometry.euclidean import is_parallel
if (point - v0).is_zero() or is_parallel(e, point - v0):
# exits through vertex i+1
return self.vertex(i + 1), PolygonPosition(
PolygonPosition.VERTEX,
vertex=(i + 1) % len(self.vertices()),
)
else:
# exits through vertex i
return v0, PolygonPosition(
PolygonPosition.VERTEX, vertex=i
)
pass
v0 = v0 + e
# Our loop has terminated. This can mean one of several errors...
pos = self.get_point_position(point)
if pos.is_outside():
raise ValueError("Started with point outside polygon")
raise ValueError(
"Point on boundary of polygon and direction not pointed into the polygon."
)
[docs] def flow_map(self, direction):
r"""
Return a polygonal map associated to the flow in ``direction`` in this
polygon.
EXAMPLES::
sage: from flatsurf.geometry.polygon import Polygon
sage: S = Polygon(vertices=[(0,0),(2,0),(2,2),(1,2),(0,2),(0,1)])
sage: S.flow_map((0,1))
Flow polygon map:
3 2
0
top lengths: [1, 1]
bot lengths: [2]
sage: S.flow_map((1,1))
Flow polygon map:
3 2 1
4 5 0
top lengths: [1, 1, 2]
bot lengths: [1, 1, 2]
sage: S.flow_map((-1,-1))
Flow polygon map:
0 5 4
1 2 3
top lengths: [2, 1, 1]
bot lengths: [2, 1, 1]
sage: K.<sqrt2> = NumberField(x^2 - 2, embedding=AA(2).sqrt())
sage: S.flow_map((sqrt2,1))
Flow polygon map:
3 2 1
4 5 0
top lengths: [1, 1, 2*sqrt2]
bot lengths: [sqrt2, sqrt2, 2]
"""
from sage.all import vector
direction = vector(direction)
DP = direction.parent()
P = self.base_ring().fraction_field() ** 2
if DP != P:
P = cm.common_parent(DP, P)
ring = P.base_ring()
direction = direction.change_ring(ring)
else:
ring = P.base_ring()
# first compute the transversal length of each edge
t = P([direction[1], -direction[0]])
lengths = [t.dot_product(e) for e in self.edges()]
n = len(lengths)
for i in range(n):
j = (i + 1) % len(lengths)
l0 = lengths[i]
l1 = lengths[j]
if l0 >= 0 and l1 < 0:
rt = j
if l0 > 0 and l1 <= 0:
rb = j
if l0 <= 0 and l1 > 0:
lb = j
if l0 < 0 and l1 >= 0:
lt = j
if rt < lt:
top_lengths = lengths[rt:lt]
top_labels = list(range(rt, lt))
else:
top_lengths = lengths[rt:] + lengths[:lt]
top_labels = list(range(rt, n)) + list(range(lt))
top_lengths = [-x for x in reversed(top_lengths)]
top_labels.reverse()
if lb < rb:
bot_lengths = lengths[lb:rb]
bot_labels = list(range(lb, rb))
else:
bot_lengths = lengths[lb:] + lengths[:rb]
bot_labels = list(range(lb, n)) + list(range(rb))
from flatsurf.geometry.interval_exchange_transformation import (
FlowPolygonMap,
)
return FlowPolygonMap(
ring, bot_labels, bot_lengths, top_labels, top_lengths
)
[docs] def flow(self, point, holonomy, translation=None):
r"""
Flow a point in the direction of holonomy for the length of the
holonomy, or until the point leaves the polygon. Note that ValueErrors
may be thrown if the point is not in the polygon, or if it is on the
boundary and the holonomy does not point into the polygon.
INPUT:
- ``point`` -- a point in the closure of the polygon (vector over the underlying base_ring())
- ``holonomy`` -- direction and magnitude of motion (vector over the underlying base_ring())
- ``translation`` -- optional translation to applied to the polygon (vector over the underlying base_ring())
OUTPUT:
- The point within the polygon where the motion stops (or leaves the polygon)
- The amount of holonomy left to flow
- a PolygonPosition object representing the combinatorial position of the stopping point
EXAMPLES::
sage: from flatsurf.geometry.polygon import polygons
sage: s = polygons.square()
sage: V = QQ**2
sage: p = V((1/2, 1/2))
sage: w = V((2, 0))
sage: s.flow(p, w)
((1, 1/2), (3/2, 0), point positioned on interior of edge 1 of polygon)
"""
from flatsurf.geometry.polygon import PolygonPosition
V = self.base_ring().fraction_field() ** 2
if holonomy == V.zero():
# not flowing at all!
return (
point,
V.zero(),
self.get_point_position(point, translation=translation),
)
if translation is None:
v0 = self.vertex(0)
else:
v0 = self.vertex(0) + translation
for i in range(len(self.vertices())):
e = self.edge(i)
from sage.all import matrix
m = matrix([[e[0], -holonomy[0]], [e[1], -holonomy[1]]])
try:
ret = m.inverse() * (point - v0)
s = ret[0]
t = ret[1]
# What if the matrix is non-invertible?
# s is location it intersects on edge, t is the portion of the holonomy to reach this intersection
if t > 0 and 0 <= s and s <= 1:
# The ray passes through edge i.
if t > 1:
# the segment from point with the given holonomy stays within the polygon
return (
point + holonomy,
V.zero(),
PolygonPosition(PolygonPosition.INTERIOR),
)
if s == 1:
# exits through vertex i+1
v0 = v0 + e
return (
v0,
point + holonomy - v0,
PolygonPosition(
PolygonPosition.VERTEX,
vertex=(i + 1) % len(self.vertices()),
),
)
if s == 0:
# exits through vertex i
return (
v0,
point + holonomy - v0,
PolygonPosition(
PolygonPosition.VERTEX, vertex=i
),
)
# exits through vertex i
# exits through interior of edge i
prod = t * holonomy
return (
point + prod,
holonomy - prod,
PolygonPosition(
PolygonPosition.EDGE_INTERIOR, edge=i
),
)
except ZeroDivisionError:
# can safely ignore this error. It means that the edge and the holonomy are parallel.
pass
v0 = v0 + e
# Our loop has terminated. This can mean one of several errors...
pos = self.get_point_position(point, translation=translation)
if pos.is_outside():
raise ValueError("Started with point outside polygon")
raise ValueError(
"Point on boundary of polygon and holonomy not pointed into the polygon."
)
[docs] def circumscribing_circle(self):
r"""
Returns the circle which circumscribes this polygon.
Raises a ValueError if the polygon is not circumscribed by a circle.
EXAMPLES::
sage: from flatsurf import Polygon
sage: P = Polygon(vertices=[(0,0),(1,0),(2,1),(-1,1)])
sage: P.circumscribing_circle()
Circle((1/2, 3/2), 5/2)
"""
from flatsurf.geometry.circle import circle_from_three_points
circle = circle_from_three_points(
self.vertex(0), self.vertex(1), self.vertex(2), self.base_ring()
)
for i in range(3, len(self.vertices())):
if not circle.point_position(self.vertex(i)) == 0:
raise ValueError(
"Vertex " + str(i) + " is not on the circle."
)
return circle
[docs] def subdivide(self):
r"""
Return a list of triangles that partition this polygon.
For each edge of the polygon one triangle is created that joins this
edge to the
:meth:`~.EuclideanPolygons.ParentMethods.centroid` of
this polygon.
EXAMPLES::
sage: from flatsurf import polygons
sage: P = polygons.regular_ngon(3); P
Polygon(vertices=[(0, 0), (1, 0), (1/2, 1/2*a)])
sage: print(P.subdivide())
[Polygon(vertices=[(0, 0), (1, 0), (1/2, 1/6*a)]),
Polygon(vertices=[(1, 0), (1/2, 1/2*a), (1/2, 1/6*a)]),
Polygon(vertices=[(1/2, 1/2*a), (0, 0), (1/2, 1/6*a)])]
::
sage: P = polygons.regular_ngon(4)
sage: print(P.subdivide())
[Polygon(vertices=[(0, 0), (1, 0), (1/2, 1/2)]),
Polygon(vertices=[(1, 0), (1, 1), (1/2, 1/2)]),
Polygon(vertices=[(1, 1), (0, 1), (1/2, 1/2)]),
Polygon(vertices=[(0, 1), (0, 0), (1/2, 1/2)])]
Sometimes alternating with :meth:`subdivide_edges` can produce a more
uniform subdivision::
sage: P = polygons.regular_ngon(4)
sage: print(P.subdivide_edges(2).subdivide())
[Polygon(vertices=[(0, 0), (1/2, 0), (1/2, 1/2)]),
Polygon(vertices=[(1/2, 0), (1, 0), (1/2, 1/2)]),
Polygon(vertices=[(1, 0), (1, 1/2), (1/2, 1/2)]),
Polygon(vertices=[(1, 1/2), (1, 1), (1/2, 1/2)]),
Polygon(vertices=[(1, 1), (1/2, 1), (1/2, 1/2)]),
Polygon(vertices=[(1/2, 1), (0, 1), (1/2, 1/2)]),
Polygon(vertices=[(0, 1), (0, 1/2), (1/2, 1/2)]),
Polygon(vertices=[(0, 1/2), (0, 0), (1/2, 1/2)])]
"""
vertices = self.vertices()
center = self.centroid()
from flatsurf import Polygon
return [
Polygon(
vertices=(
vertices[i],
vertices[(i + 1) % len(vertices)],
center,
),
)
for i in range(len(vertices))
]
[docs] def subdivide_edges(self, parts=2):
r"""
Return a copy of this polygon whose edges have been split into
``parts`` equal parts each.
INPUT:
- ``parts`` -- a positive integer (default: 2)
EXAMPLES::
sage: from flatsurf import polygons
sage: P = polygons.regular_ngon(3); P
Polygon(vertices=[(0, 0), (1, 0), (1/2, 1/2*a)])
sage: P.subdivide_edges(1) == P
True
sage: P.subdivide_edges(2)
Polygon(vertices=[(0, 0), (1/2, 0), (1, 0), (3/4, 1/4*a), (1/2, 1/2*a), (1/4, 1/4*a)])
sage: P.subdivide_edges(3)
Polygon(vertices=[(0, 0), (1/3, 0), (2/3, 0), (1, 0), (5/6, 1/6*a), (2/3, 1/3*a), (1/2, 1/2*a), (1/3, 1/3*a), (1/6, 1/6*a)])
"""
if parts < 1:
raise ValueError("parts must be a positive integer")
steps = [e / parts for e in self.edges()]
from flatsurf import Polygon
return Polygon(edges=[e for e in steps for p in range(parts)])
[docs] def j_invariant(self):
r"""
Return the Kenyon-Smille J-invariant of this polygon.
The base ring of the polygon must be a number field.
The output is a triple ``(Jxx, Jyy, Jxy)`` that corresponds
respectively to the Sah-Arnoux-Fathi invariant of the vertical flow,
the Sah-Arnoux-Fathi invariant of the horizontal flow and the `xy`-matrix.
EXAMPLES::
sage: from flatsurf import polygons
sage: polygons.right_triangle(1/3,1).j_invariant()
(
[0 0]
(0), (0), [1 0]
)
The regular 8-gon::
sage: polygons.regular_ngon(8).j_invariant()
(
[2 2]
(0), (0), [2 1]
)
(
[ 0 3/2]
(1/2), (-1/2), [3/2 0]
)
Some extra debugging::
sage: K.<a> = NumberField(x^3 - 2, embedding=AA(2)**(1/3))
sage: ux = 1 + a + a**2
sage: uy = -2/3 + a
sage: vx = 1/5 - a**2
sage: vy = a + 7/13*a**2
sage: from flatsurf import Polygon
sage: p = Polygon(edges=[(ux, uy), (vx,vy), (-ux-vx,-uy-vy)], base_ring=K)
sage: Jxx, Jyy, Jxy = p.j_invariant()
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons.Simple.Convex.ParentMethods._wedge_product(ux.vector(), vx.vector()) == Jxx
True
sage: EuclideanPolygons.Simple.Convex.ParentMethods._wedge_product(uy.vector(), vy.vector()) == Jyy
True
"""
from sage.all import QQ, matrix
if self.base_ring() is QQ:
raise NotImplementedError
K = self.base_ring()
try:
V, from_V, to_V = K.vector_space()
except (AttributeError, ValueError):
raise ValueError(
"the surface needs to be define over a number field"
)
dim = K.degree()
M = K ** (dim * (dim - 1) // 2)
Jxx = Jyy = M.zero()
Jxy = matrix(K, dim)
vertices = list(self.vertices())
vertices.append(vertices[0])
for i in range(len(vertices) - 1):
a = to_V(vertices[i][0])
b = to_V(vertices[i][1])
c = to_V(vertices[i + 1][0])
d = to_V(vertices[i + 1][1])
Jxx += self._wedge_product(a, c)
Jyy += self._wedge_product(b, d)
Jxy += self._tensor_product(a, d)
Jxy -= self._tensor_product(c, b)
return (Jxx, Jyy, Jxy)
@staticmethod
def _wedge_product(v, w):
r"""
Return the wedge product of ``v`` and ``w``.
This is a helper method for :meth:`j_invariant`.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons.Simple.Convex.ParentMethods._wedge_product(vector((1, 2)), vector((1, 2)))
(0)
sage: EuclideanPolygons.Simple.Convex.ParentMethods._wedge_product(vector((1, 2)), vector((2, 1)))
(-3)
sage: EuclideanPolygons.Simple.Convex.ParentMethods._wedge_product(vector((1, 2, 3)), vector((2, 3, 4)))
(-1, -2, -1)
"""
d = len(v)
assert len(w) == d
R = v.base_ring()
from sage.all import free_module_element
return free_module_element(
R,
d * (d - 1) // 2,
[
(v[i] * w[j] - v[j] * w[i])
for i in range(d - 1)
for j in range(i + 1, d)
],
)
@staticmethod
def _tensor_product(u, v):
r"""
Return the tensor product of ``u`` and ``v``.
This is a helper method for :meth:`j_invariant`.
EXAMPLES::
sage: from flatsurf.geometry.categories import EuclideanPolygons
sage: EuclideanPolygons.Simple.Convex.ParentMethods._tensor_product(vector((2, 3, 5)), vector((7, 11, 13)))
[14 21 35]
[22 33 55]
[26 39 65]
"""
from sage.all import vector
u = vector(u)
v = vector(v)
d = len(u)
R = u.base_ring()
assert len(u) == len(v) and v.base_ring() == R
from sage.all import matrix
return matrix(
R, d, [u[i] * v[j] for j in range(d) for i in range(d)]
)
[docs] def is_isometric(self, other, certificate=False):
r"""
Return whether ``self`` and ``other`` are isometric convex polygons via an orientation
preserving isometry.
If ``certificate`` is set to ``True`` return also a pair ``(index, rotation)``
of an integer ``index`` and a matrix ``rotation`` such that the given rotation
matrix identifies this polygon with the other and the edges 0 in this polygon
is mapped to the edge ``index`` in the other.
.. TODO::
Implement ``is_linearly_equivalent`` and ``is_similar``.
EXAMPLES::
sage: from flatsurf import Polygon, polygons
sage: S = polygons.square()
sage: S.is_isometric(S)
True
sage: U = matrix(2,[0,-1,1,0]) * S
sage: U.is_isometric(S)
True
sage: x = polygen(QQ)
sage: K.<sqrt2> = NumberField(x^2 - 2, embedding=AA(2)**(1/2))
sage: S = S.change_ring(K)
sage: U = matrix(2, [sqrt2/2, -sqrt2/2, sqrt2/2, sqrt2/2]) * S
sage: U.is_isometric(S)
True
sage: U2 = Polygon(edges=[(1,0), (sqrt2/2, sqrt2/2), (-1,0), (-sqrt2/2, -sqrt2/2)])
sage: U2.is_isometric(U)
False
sage: U2.is_isometric(U, certificate=True)
(False, None)
sage: S = Polygon(edges=[(1,0), (sqrt2/2, 3), (-2,3), (-sqrt2/2+1, -6)])
sage: T = Polygon(edges=[(sqrt2/2,3), (-2,3), (-sqrt2/2+1, -6), (1,0)])
sage: isometric, cert = S.is_isometric(T, certificate=True)
sage: assert isometric
sage: shift, rot = cert
sage: Polygon(edges=[rot * S.edge((k + shift) % 4) for k in range(4)]).translate(T.vertex(0)) == T
True
sage: T = (matrix(2, [sqrt2/2, -sqrt2/2, sqrt2/2, sqrt2/2]) * S).translate((3,2))
sage: isometric, cert = S.is_isometric(T, certificate=True)
sage: assert isometric
sage: shift, rot = cert
sage: Polygon(edges=[rot * S.edge(k + shift) for k in range(4)]).translate(T.vertex(0)) == T
True
"""
from flatsurf.geometry.polygon import EuclideanPolygon
if not isinstance(other, EuclideanPolygon):
raise TypeError("other must be a polygon")
if not other.is_convex():
raise TypeError("other must be convex")
n = len(self.vertices())
if len(other.vertices()) != n:
return False
sedges = self.edges()
oedges = other.edges()
slengths = [x**2 + y**2 for x, y in sedges]
olengths = [x**2 + y**2 for x, y in oedges]
for i in range(n):
if slengths == olengths:
# we have a match of lengths after a shift by i
xs, ys = sedges[0]
xo, yo = oedges[0]
from sage.all import matrix
ms = matrix(2, [xs, -ys, ys, xs])
mo = matrix(2, [xo, -yo, yo, xo])
rot = mo * ~ms
assert rot.det() == 1 and (rot * rot.transpose()).is_one()
assert oedges[0] == rot * sedges[0]
if all(oedges[i] == rot * sedges[i] for i in range(1, n)):
return (
(True, (0 if i == 0 else n - i, rot))
if certificate
else True
)
olengths.append(olengths.pop(0))
oedges.append(oedges.pop(0))
return (False, None) if certificate else False
[docs] def is_translate(self, other, certificate=False):
r"""
Return whether ``other`` is a translate of ``self``.
EXAMPLES::
sage: from flatsurf import Polygon
sage: S = Polygon(vertices=[(0,0), (3,0), (1,1)])
sage: T1 = S.translate((2,3))
sage: S.is_translate(T1)
True
sage: T2 = Polygon(vertices=[(-1,1), (1,0), (2,1)])
sage: S.is_translate(T2)
False
sage: T3 = Polygon(vertices=[(0,0), (3,0), (2,1)])
sage: S.is_translate(T3)
False
sage: S.is_translate(T1, certificate=True)
(True, (0, 1))
sage: S.is_translate(T2, certificate=True)
(False, None)
sage: S.is_translate(T3, certificate=True)
(False, None)
"""
if type(self) is not type(other):
raise TypeError
n = len(self.vertices())
if len(other.vertices()) != n:
return False
sedges = self.edges()
oedges = other.edges()
for i in range(n):
if sedges == oedges:
return (True, (i, 1)) if certificate else True
oedges.append(oedges.pop(0))
return (False, None) if certificate else False
[docs] def is_half_translate(self, other, certificate=False):
r"""
Return whether ``other`` is a translate or half-translate of ``self``.
If ``certificate`` is set to ``True`` then return also a pair ``(orientation, index)``.
EXAMPLES::
sage: from flatsurf import Polygon
sage: S = Polygon(vertices=[(0,0), (3,0), (1,1)])
sage: T1 = S.translate((2,3))
sage: S.is_half_translate(T1)
True
sage: T2 = Polygon(vertices=[(-1,1), (1,0), (2,1)])
sage: S.is_half_translate(T2)
True
sage: T3 = Polygon(vertices=[(0,0), (3,0), (2,1)])
sage: S.is_half_translate(T3)
False
sage: S.is_half_translate(T1, certificate=True)
(True, (0, 1))
sage: half_translate, cert = S.is_half_translate(T2, certificate=True)
sage: assert half_translate
sage: shift, rot = cert
sage: Polygon(edges=[rot * S.edge(k + shift) for k in range(3)]).translate(T2.vertex(0)) == T2
True
sage: S.is_half_translate(T3, certificate=True)
(False, None)
"""
if type(self) is not type(other):
raise TypeError
n = len(self.vertices())
if len(other.vertices()) != n:
return False
sedges = self.edges()
oedges = other.edges()
for i in range(n):
if sedges == oedges:
return (True, (i, 1)) if certificate else True
oedges.append(oedges.pop(0))
assert oedges == other.edges()
oedges = [-e for e in oedges]
for i in range(n):
if sedges == oedges:
return (
(True, (0 if i == 0 else n - i, -1))
if certificate
else True
)
oedges.append(oedges.pop(0))
return (False, None) if certificate else False