Add geometry in the hyperbolic plane

saraedum commented 2 years ago

TODO

[x] Added news/
[x] Everything is doctested.
[x] ~~The intersection algorithm is unit tested with (more complicated) random inputs.~~ (let's do this when there is a problem.)
[x] Generalize, see https://github.com/flatsurf/sage-flatsurf/pull/158#issuecomment-1103917076 (let's put this into a followup issue; should probably all go into a Geometry class.)
[x] Fix https://github.com/conda-forge/sagelib-feedstock/issues/136 with a patch to https://github.com/conda-forge/conda-forge-repodata-patches-feedstock
[x] perpendicular to a given geodesic passing through a given point (wikipedia)
[x] perpendicular to two geodesics (wikipedia)
[x] midpoint between two points (wikipedia)
[x] perpendicular bisector of a line segment
[x] ~~bisect an angle (wikipedia)~~ (will do later.)
[x] ~~centroid of a polygon~~ (will do later.)
[x] Make I in HyperbolicPlane() work.
[x] ~~RealSet.unbounded_below_closed cannot be constructed from number field elements even if embedded.~~ (Left a note in the implementation. Let's wait for it to become an actual issue.)
[x] Check that unoriented objects have no leading signs in their _repr_.
[x] Add section to the documentation.

Random plots

from flatsurf.geometry.hyperbolic import HyperbolicPlane
H = HyperbolicPlane(QQ)

a = H(-4)
b = H(5/4*I)
c = H.half_circle(0, 1).intersection(H.geodesic(1/3, 3))
d = H(1)
e = H(1/2 * I)
f = H(0)
g = H(-1/3)
h = H(-3/2)

H.polygon([H.geodesic(a, b).left_half_space(), H.geodesic(b, c).left_half_space(), H.geodesic(c, d).left_half_space()]).plot(model="klein") +\
H.polygon([H.geodesic(d, c).left_half_space(), H.geodesic(c, e).left_half_space(), H.geodesic(e, f).left_half_space()]).plot(model="klein", color='#ffefff') +\
H.polygon([H.geodesic(g, e).left_half_space(), H.geodesic(e, c).left_half_space(), H.geodesic(c, b).left_half_space(), H.geodesic(b, h).left_half_space()]).plot(model="klein", color="#ffffef") +\
H.polygon([H.geodesic(h, b).left_half_space(), H.geodesic(b, a).left_half_space()]).plot(model="klein", color="#ffefff") +\
H.polygon([H.geodesic(f, e).left_half_space(), H.geodesic(e, g).left_half_space()]).plot(model="klein", color="#efefff")

And the same in the upper half plane

References

Cannon, Floyd, Kenyon, Parry "Hyperbolic geometry" (1997)

videlec commented 2 years ago

Just for reference, the interface for convex polygons in R^d is as follows

Polyhedron(vertices=XXX)    or    Polyhedron(eqns=XXX, ieqs=YYY)

Where eqns is a list of triples (a, b, c) meaning a + bx + cy = 0 and ieqs is a list of triples (a, b, c) meaning a + bx + cy >= 0.

videlec commented 2 years ago

Isometries in the Klein model are particularly nice : they are induced from the isometries of the ambient R^{1,2} (that is the group SO(1,2)) in particular projective (x, y) -> ((ax + by + c) / (gx + hy + i) , (dx + ey + f) / (gx + hy + i)). There is an actual isomorphism of Lie groupes between PSL(2,R) and SO(1, 2). The direction PSL(2,R) -> SO(1,2) is nice and simple.

Do you intend to have a dedicated Parent/Element for isometries ?

saraedum commented 2 years ago

Isometries in the Klein model are particularly nice : they are induced from the isometries of the ambient R^{1,2} (that is the group SO(1,2)) in particular projective (x, y) -> ((ax + by + c) / (gx + hy + i) , (dx + ey + f) / (gx + hy + i)). There is an actual isomorphism of Lie groupes between PSL(2,R) and SO(1, 2). The direction PSL(2,R) -> SO(1,2) is nice and simple.

@slel I think you were wondering about this part.

Do you intend to have a dedicated Parent/Element for isometries ?

No, I think for the time being I just want to be able to apply a matrix to a convex set.

videlec commented 2 years ago

Do you intend to have a dedicated Parent/Element for isometries ?

No, I think for the time being I just want to be able to apply a matrix to a convex set.

A 2x2 matrix in SL(2,R) or a 3x3 matrix in SO(1,2)?

videlec commented 2 years ago

Generators of SO(1, 2) are

[ cos(theta)  -sin(theta)  0 ]
[ sin(theta)   cos(theta)  0 ]
[ 0            0           1 ]

and

[ cosh(t) 0 sinh(t) ]
[ 0       1 0       ]
[ sinh(t) 0 cosh(t) ]

More precisely any matrix in SO(1, 2) can be written as a composition k1 a k2 with k1, k2 of the first type and a of the second type. This is the so-called KAK decomposition.

videlec commented 2 years ago

I wrote down the lift SL(2,R) -> SO(1, 2). The normalization of the Klein model we use is not very nice as this lift does not preserve integrality (there is a division by 2 in the formula).

videlec commented 2 years ago

@saraedum This comment is not directly related to flatsurf (which is very much SL(2, R)-focused). It turns out that working with a given quadratic form Q of signature (1, 2) in R^3 is very convenient (eg when dealing with quaternion algebras or when defining Coxeter groups such as triangle groups). In this branch we fix the canonical Q(x,y,z) = -x^2 - y^2 + z^2 but there is not much reason to do so.

videlec commented 2 years ago

@saraedum hmmm

sage: H.geodesic(p0, p1).apply_isometry(m)
{2*(x^2 + y^2) - 5*x + 3 = 0}
sage: H.geodesic(q0, q1)
{2*(x^2 + y^2) - 5*x + 3 = 0}
sage:  H.geodesic(p0, p1).apply_isometry(m) == H.geodesic(q0, q1)
False

videlec commented 2 years ago

@saraedum The underlying reason is

sage: H.geodesic(5, -5, -1, model="half_plane") == H.geodesic(5/13, -5/13, -1/13, model="half_plane")
False

videlec commented 2 years ago

Indeed, the following is not a projective comparison

return rich_to_bool(op, (self._a - other._a).sign())

saraedum commented 2 years ago

Indeed, the following is not a projective comparison
return rich_to_bool(op, (self._a - other._a).sign())

Yes, I noticed that as well now. I think richcmp should just be equality and not conflate this with anything else.

I now introduced a separate sorting key for the half plane intersection algorithm and a separate is_subset (which then also ignores the orientation of a geodesic.)

saraedum commented 2 years ago

@saraedum This comment is not directly related to flatsurf (which is very much SL(2, R)-focused). It turns out that working with a given quadratic form Q of signature (1, 2) in R^3 is very convenient (eg when dealing with quaternion algebras or when defining Coxeter groups such as triangle groups). In this branch we fix the canonical Q(x,y,z) = -x^2 - y^2 + z^2 but there is not much reason to do so.

Sorry, but I don't understand what you are trying to say here.

videlec commented 2 years ago

@saraedum This comment is not directly related to flatsurf (which is very much SL(2, R)-focused). It turns out that working with a given quadratic form Q of signature (1, 2) in R^3 is very convenient (eg when dealing with quaternion algebras or when defining Coxeter groups such as triangle groups). In this branch we fix the canonical Q(x,y,z) = -x^2 - y^2 + z^2 but there is not much reason to do so.

Sorry, but I don't understand what you are trying to say here.

The Klein model is the geometry of lines intersected with x^2 + y^2 < 1 (which is secretely the projectivization of -x^2 - y^2 + z^2 > 0). But it is sometimes more natural to consider an other quadratic form in R^3 with signature (1, 2) and the associated Klein model P ( {(x, y, z): Q(x, y, z) > 0} ). Especially when dealing with subgroups of isometries coming from Coxeter groups. I was wondering whether one could make the code flexible enough so that no specific feature of -x^2 - y^2 + z^2 is used.

saraedum commented 2 years ago

@saraedum This comment is not directly related to flatsurf (which is very much SL(2, R)-focused). It turns out that working with a given quadratic form Q of signature (1, 2) in R^3 is very convenient (eg when dealing with quaternion algebras or when defining Coxeter groups such as triangle groups). In this branch we fix the canonical Q(x,y,z) = -x^2 - y^2 + z^2 but there is not much reason to do so.

Sorry, but I don't understand what you are trying to say here.

The Klein model is the geometry of lines intersected with x^2 + y^2 < 1 (which is secretely the projectivization of -x^2 - y^2 + z^2 > 0). But it is sometimes more natural to consider an other quadratic form in R^3 with signature (1, 2) and the associated Klein model P ( {(x, y, z): Q(x, y, z) > 0} ). Especially when dealing with subgroups of isometries coming from Coxeter groups. I was wondering whether one could make the code flexible enough so that no specific feature of -x^2 - y^2 + z^2 is used.

We should discuss this in our weekly call. I understand what you're saying but I would like to understand the use case a bit better.

videlec commented 2 years ago

@saraedum This comment is not directly related to flatsurf (which is very much SL(2, R)-focused). It turns out that working with a given quadratic form Q of signature (1, 2) in R^3 is very convenient (eg when dealing with quaternion algebras or when defining Coxeter groups such as triangle groups). In this branch we fix the canonical Q(x,y,z) = -x^2 - y^2 + z^2 but there is not much reason to do so.

Sorry, but I don't understand what you are trying to say here.

The Klein model is the geometry of lines intersected with x^2 + y^2 < 1 (which is secretely the projectivization of -x^2 - y^2 + z^2 > 0). But it is sometimes more natural to consider an other quadratic form in R^3 with signature (1, 2) and the associated Klein model P ( {(x, y, z): Q(x, y, z) > 0} ). Especially when dealing with subgroups of isometries coming from Coxeter groups. I was wondering whether one could make the code flexible enough so that no specific feature of -x^2 - y^2 + z^2 is used.

We should discuss this in our weekly call. I understand what you're saying but I would like to understand the use case a bit better.

One example comes from Coxeter groups that are generated by reflections. The canonical quadratic form depends on the choice of reflections. For Triangle(2,3,7) I got the quadratic form

[                      1                       0 1/2*E(7)^3 + 1/2*E(7)^4]
[                      0                       1                    -1/2]
[1/2*E(7)^3 + 1/2*E(7)^4                    -1/2                       1]

where E(7) is the seventh root of unity.

videlec commented 2 years ago

Some additional useful methods (maybe for a later PR)

perpendicular to a given geodesic passing through a given point
perpendicular to two geodesics
midpoint between two points
mediatrix of two points

The constructions are on wikiedia.

videlec commented 2 years ago

Two remarks about segments.

I would like to build hyperbolic segments between pairs of points. It turns out that I must do

from flatsurf.geometry.hyperbolic import HyperbolicPlane
H = HyperbolicPlane(QQ)
A = H.point(1/2, 1/2, model='klein')
B = H.point(-1/3, 3/5, model='klein')
s = H.segment(H.geodesic(A, B), A, B)  # !?

I find the first argument a bit weird.

Also, plotting via s.plot() does plot the region between the oriented segment and the boundary. Not just the segment.

saraedum commented 2 years ago

Also, plotting via s.plot() does plot the region between the oriented segment and the boundary. Not just the segment.

Should be fixed now.

saraedum commented 2 years ago

@videlec thanks for the usability feedback. I am a bit unsure how to address it.

I could do HyperbolicPlane.segment() accept three arguments, at most one can be a geodesic the other two can be points, if there is only one point, a geodesic must be present, if that geodesic is oriented, then …

It gets a bit complicated (to document at least) so I now added HyperbolicPoint.segment() so you can say A.segment(B) to get the segment from A to B in the above example. Do you think that's good enough?

videlec commented 2 years ago

It gets a bit complicated (to document at least) so I now added HyperbolicPoint.segment() so you can say A.segment(B) to get the segment from A to B in the above example. Do you think that's good enough?

A.segment(B) looks like a perfect alternative interface.

videlec commented 2 years ago

Still about geodesic segments, when playing with inexact ring the constructor is very unhappy

sage: from flatsurf.geometry.hyperbolic import HyperbolicPlane
sage: H = HyperbolicPlane(RDF)
sage: p0 = H.point(0, 1, 'half_plane')
sage: p1 = H.point(1, 1, 'half_plane')
sage: p0.segment(p1)
Traceback (most recent call last)
...
NotImplementedError: cannot decide containment for null sets

Is there a way to plot geodesic segments between inexact points?

saraedum commented 2 years ago

p0.segment(p1, check=False, assume_normalized=True)

should work.

saraedum commented 2 years ago

@videlec could you clarify what's the mediatrix of two points? Wikipedia says

In Catholic Mariology, the title Mediatrix refers to the intercessory role of the Blessed Virgin Mary as a mediator in the salvific redemption by her son Jesus Christ and that he bestows graces through her.

Do you mean bisecting an angle? That seems like the more common term to me.

slel commented 2 years ago

French: médiatrice English: perpendicular bisector

videlec commented 2 years ago

@videlec could you clarify what's the mediatrix of two points? Wikipedia says

In Catholic Mariology, the title Mediatrix refers to the intercessory role of the Blessed Virgin Mary as a mediator in the salvific redemption by her son Jesus Christ and that he bestows graces through her.

Do you mean bisecting an angle? That seems like the more common term to me.

It is the french "médiatrice" which seems to be translated as "perpendicular bisector". Given two points, the perpendicular bisector is the orthogonal to the line between these two points and passing by their middle.

saraedum commented 1 year ago

Btw. flake8 --select RST --rst-roles meth,class,func --rst-directives SEEALSO,TODO flatsurf/**/*.py is awesome to find problems in documentation (since Sphinx error messages are often useless, see https://github.com/sphinx-doc/sphinx/issues/6036) You need to mamba install flake8-rst-docstrings for this.

github-actions[bot] commented 1 year ago

Documentation preview for this PR is ready! :tada: Built with commit: 52a1773e61cd1c45fb0f2c0325772d1ca77e9830

flatsurf / sage-flatsurf