Fix moebius_transform, midpoint and perpendicular_bisector

[kliem]

We fix several related issues in hyperbolic geometry:

• moebius_transform wrong in the orientation-reversing
  (negative determinant) case
• midpoint and perpendicular_bisector
  gave different results for (a, b) and (b, a)

Fuzzing geometry/hyperbolic_space/hyperbolic_geodesic.py
doctests uncovered these issues. See #29935, #29962, #29963.

UHP = HyperbolicPlane().UHP()
sage: g = UHP.random_geodesic()   # Good one
sage: g.start(), g.end()
(Point in UHP -9.26728445125634 + 7.15068223909426*I,
 Point in UHP -8.57515638592863 + 0.982992065975858*I)
sage: g = UHP.random_geodesic()   # Bad one
sage: g.start(), g.end()
(Point in UHP -4.92591958004554 + 7.81075974592370*I,
 Point in UHP -8.44416583993011 + 3.87025183089797*I)

Depends on #29962

CC: @slel @orlitzky @tscrim @greglaun @sagetrac-jhonrubia6

Component: geometry

Keywords: hyperbolic geodesic

Author: Travis Scrimshaw

Branch/Commit: 234604b

Reviewer: Samuel Lelièvre

Issue created by migration from https://trac.sagemath.org/ticket/29936

[kliem]

Branch: public/29936

[kliem]

Author: Jonathan Kliem

[kliem]

comment:1

The first thing is to use the flag # abs tol to obtain more meaningful doctest failures at random states. I get the following (not all during the same run):

File "src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py", line 912, in sage.geometry.hyperbolic_space.hyperbolic_geodesic.HyperbolicGeodesic.perpendicular_bisector
Failed example:
    abs(h.intersection(g)[0].coordinates() - g.midpoint().coordinates())  # abs tol 1e-9
Expected:
    0
Got:
    0.309903204800636
Tolerance exceeded:
    0 vs 0.309903204800636, tolerance 4e-1 > 1e-9
**********************************************************************
File "src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py", line 1376, in sage.geometry.hyperbolic_space.hyperbolic_geodesic.HyperbolicGeodesicUHP.perpendicular_bisector
Failed example:
    abs(c(g.intersection(h)[0]) - c(g.midpoint()))  # abs tol 1e-9
Expected:
    0
Got:
    9.74965482957673
Tolerance exceeded:
    0 vs 9.74965482957673, tolerance 1e1 > 1e-9

**********************************************************************
File "src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py", line 1914, in sage.geometry.hyperbolic_space.hyperbolic_geodesic.HyperbolicGeodesicUHP._crossratio_matrix
Failed example:
    abs(moebius_transform(A, p1))  # abs tol 1e-9
Expected:
    0
Got:
    1.82859072578637
Tolerance exceeded:
    0 vs 1.82859072578637, tolerance 2e0 > 1e-9
**********************************************************************
File "src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py", line 1916, in sage.geometry.hyperbolic_space.hyperbolic_geodesic.HyperbolicGeodesicUHP._crossratio_matrix
Failed example:
    abs(moebius_transform(A, p2) - 1)  # abs tol 1e-9
Expected:
    0
Got:
    0.619923876802936
Tolerance exceeded:
    0 vs 0.619923876802936, tolerance 7e-1 > 1e-9

This tells me that the preciseness that is claimed with the doctests is far from valid.

---

New commits:

23ed583	fix doctest in hyperbolic_space/hyperbolic_point
5283dc4	use abs tol flag

[kliem]

Commit: 5283dc4

[kliem]

comment:2

I know very little about numerics. From what I remember it makes sense to consider the relative error as well:

sage: from sage.geometry.hyperbolic_space.hyperbolic_isometry import moebius_transform
sage: from sage.geometry.hyperbolic_space.hyperbolic_geodesic import HyperbolicGeodesicUHP
sage: UHP = HyperbolicPlane().UHP()
sage: def foo():
....:     (p1, p2, p3) = [UHP.random_point().coordinates() for k in range(3)]
....:     A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3)
....:     a = abs(moebius_transform(A, p1))
....:     b = abs((moebius_transform(A, p2) - 1))
....:     sum_of_norms = A.norm() + p1.norm() + p2.norm() + p3.norm()
....:     return min(a, a/sum_of_norms), min(b, b/sum_of_norms)
....: 
sage: max(foo()[1] for _ in range(1000))
0.157952578757514
sage: max(foo()[0] for _ in range(1000))
0.0615941849084595

The situation is far worse with the other two tests, which are essentially the same, if I understand correctly:

sage: def foo():
....:     g = HyperbolicPlane().PD().random_geodesic()
....:     h = g.perpendicular_bisector()
....:     error = abs(h.intersection(g)[0].coordinates() - g.midpoint().coordinates())
....:     sum_of_norms = h.intersection(g)[0].coordinates().norm() + g.midpoint().coordinates().norm()
....:     return min(error, error/sum_of_norms)
....: 
sage: max(foo() for _ in range(10))
0.908342722333419

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

7b244c0

modify doctests to the extend that they hold with fuzz

[sagetrac-git]

Changed commit from 5283dc4 to 7b244c0

[kliem]

comment:4

I tested 10 times and not a single test failed. I hope that is good enough.

[kliem]

Changed keywords from none to hyperbolic geodesic, hyperbolic point

[orlitzky]

comment:6

I've added a few people to CC who may know how these tests are supposed to work. Those errors look pretty large to me.

For something less substantial: I personally always try to limit my use of abs tol to the TESTS block (as opposed to EXAMPLES). The parsing of that comment is an internal detail of our doctest framework, and users who see it in the reference manual generally won't know what it means. Even those users who do understand it have to do a mental calculation when running the example to see if their output matches what is expected. Maybe you agree?

[kliem]

Dependencies: #29962

[kliem]

Changed commit from 7b244c0 to 8292b04

[kliem]

comment:8

Theoretically, there is not need to base this on top of #29962. But basically #29962 discovers those problems. So the doctest works fine until you test it on top of #29962 with a different random seed than 0.

---

Last 10 new commits:

da1c6be	start from a "random" random seed for doctesting
b7b836d	make random seed reproducible
eedbe5e	document random_seed
998b1b9	default random seed 0 for now
1d7b00e	dash instead of underscore for command line options
b31e2d5	Merge branch 'public/29962' of git://trac.sagemath.org/sage into public/29962-reb
2f30dd9	small fixes
b62f781	doctests do not start from a random seed by default yet
1d99129	fix merge conflict
8292b04	Merge branch 'public/29962-reb2' of git://trac.sagemath.org/sage into public/29936-reb

[kliem]

Changed branch from public/29936 to public/29936-reb

[kliem]

comment:9

Replying to @orlitzky:

I've added a few people to CC who may know how these tests are supposed to work. Those errors look pretty large to me.

I agree. That's why I opened #29939 to fix this.

For something less substantial: I personally always try to limit my use of abs tol to the TESTS block (as opposed to EXAMPLES). The parsing of that comment is an internal detail of our doctest framework, and users who see it in the reference manual generally won't know what it means. Even those users who do understand it have to do a mental calculation when running the example to see if their output matches what is expected. Maybe you agree?

I agree. I moved that stuff to TESTS. Especially with the modified doctests, the tests are really hard to read. But I would argue that for some random stuff you cannot expect a good absolute error, but only a good relative error (and the minimum is needed, because you don't want to divide by zero).

[tscrim]

comment:10

Replying to @kliem:

Replying to @orlitzky:

I've added a few people to CC who may know how these tests are supposed to work. Those errors look pretty large to me.

I agree. That's why I opened #29939 to fix this.

It looks like you are comparing using relative tolerance, which if the values are close to 0, could cause large variations. Since you are taking things to be more random, it is quite possible you could divide by something close to 0. So I am not sure this is a bug as instead it could a problem that you have introduced by changing the doctests.

For something less substantial: I personally always try to limit my use of abs tol to the TESTS block (as opposed to EXAMPLES). The parsing of that comment is an internal detail of our doctest framework, and users who see it in the reference manual generally won't know what it means. Even those users who do understand it have to do a mental calculation when running the example to see if their output matches what is expected. Maybe you agree?

I agree. I moved that stuff to TESTS. Especially with the modified doctests, the tests are really hard to read.

I disagree. The problem is numerics and working with inexact fields. You cannot get around this. I don't exactly understand what your modified doctests are doing now other than they are comparing something relatively rather than absolutely.

But I would argue that for some random stuff you cannot expect a good absolute error, but only a good relative error (and the minimum is needed, because you don't want to divide by zero).

You have two competing issues: large values want relative error, small values want absolute error. You have too much randomness now to deal with.

[orlitzky]

comment:11

Replying to @tscrim:

It looks like you are comparing using relative tolerance, which if the values are close to 0, could cause large variations. Since you are taking things to be more random, it is quite possible you could divide by something close to 0. So I am not sure this is a bug as instead it could a problem that you have introduced by changing the doctests.

The original tests also fail if you introduce any randomness at all. I'm mainly concerned about the errors that Jonathan noticed right away, e.g.

File "src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py", line 1914, in sage.geometry.hyperbolic_space.hyperbolic_geodesic.HyperbolicGeodesicUHP._crossratio_matrix
Failed example:
    abs(moebius_transform(A, p1))  # abs tol 1e-9
Expected:
    0
Got:
    1.82859072578637
Tolerance exceeded:
    0 vs 1.82859072578637, tolerance 2e0 > 1e-9

If you have a method that's supposed to return zero but instead returns two, something other than numerical noise might be to blame. Can that test be re-run with a failing random seed, to see what a failing example looks like?

[kliem]

comment:12

Replying to @tscrim:

Replying to @kliem:

Replying to @orlitzky:

I've added a few people to CC who may know how these tests are supposed to work. Those errors look pretty large to me.

I agree. That's why I opened #29939 to fix this.

It looks like you are comparing using relative tolerance, which if the values are close to 0, could cause large variations. Since you are taking things to be more random, it is quite possible you could divide by something close to 0. So I am not sure this is a bug as instead it could a problem that you have introduced by changing the doctests.

For something less substantial: I personally always try to limit my use of abs tol to the TESTS block (as opposed to EXAMPLES). The parsing of that comment is an internal detail of our doctest framework, and users who see it in the reference manual generally won't know what it means. Even those users who do understand it have to do a mental calculation when running the example to see if their output matches what is expected. Maybe you agree?

I agree. I moved that stuff to TESTS. Especially with the modified doctests, the tests are really hard to read.

I disagree. The problem is numerics and working with inexact fields. You cannot get around this. I don't exactly understand what your modified doctests are doing now other than they are comparing something relatively rather than absolutely.

But I would argue that for some random stuff you cannot expect a good absolute error, but only a good relative error (and the minimum is needed, because you don't want to divide by zero).

You have two competing issues: large values want relative error, small values want absolute error. You have too much randomness now to deal with.

I'm taking the minimum of the absolute and the relative error, because I want to be fair. This makes it much better. This is the original test:

sage: from sage.geometry.hyperbolic_space.hyperbolic_isometry import moebius_transform
....: from sage.geometry.hyperbolic_space.hyperbolic_geodesic import HyperbolicGeodesicUHP
....: UHP = HyperbolicPlane().UHP()
....: def foo():
....:     (p1, p2, p3) = [UHP.random_point().coordinates() for k in range(3)]
....:     A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3)
....:     a = abs(moebius_transform(A, p1))
....:     b = abs((moebius_transform(A, p2) - 1))
....:     return (a,b)
....: 
sage: max(foo()[1] for _ in range(1000))
19.6483996640312
sage: max(foo()[0] for _ in range(1000))
76.0110826766808

Here are examples that perform bad. I hope you see, why I chose to change the doctest to minimum of absolute and relative norm:

sage: (p1, p2, p3) = (-2.29920829511711 + 5.52937627459424*I, -2.37141010564321 + 5.57310680132109*I, 8.28221981185745 + 2.99165901130902*I)
sage: A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3)
sage: a = abs(moebius_transform(A, p1)); a
105.706157025922
sage: b = abs((moebius_transform(A, p2) - 1)); b
105.113014884487
sage: sum_of_norms = A.norm() + p1.norm() + p2.norm() + p3.norm(); sum_of_norms
216.645504358830

[kliem]

comment:13

Given those values, I think the original doctests are unacceptable. We know those exactness claims are wrong. Even if it succeeds with random seed 0, we shouldn't leave it as it is.

[tscrim]

comment:14

So I agree with changing the test in hyperbolic_point.py: The equality test is perhaps a bit too narrow and needs to be expanded to allow for some fuzziness (up to some precision set globally).

For the issue in the geodesic, the mathematics is correct:

sage: var('a0,a1,a2,b0,b1,b2')
(a0, a1, a2, b0, b1, b2)
sage: p0 = a0+b0*I
sage: p1 = a1+b1*I
sage: p2 = a2+b2*I
sage: p0 = a0+b0*I
sage: p1 = a1+b1*I
sage: p2 = a2+b2*I
sage: A = HyperbolicGeodesicUHP._crossratio_matrix(p0, p1, p2)
sage: moebius_transform(A, p0)
0
sage: moebius_transform(A, p1).factor()
1
sage: moebius_transform(A, p2)
+Infinity

So I think the issue comes from this in the moebius_transform function:

        if a * d - b * c < 0:
            w = z.conjugate()  # Reverses orientation
        else:
            w = z

At least, I only seem to notice it when the determinant of A has a negative real value. So if I remove the conjugate(), then the issue appears to be fixed.

The problem is not the doctests, which are correct and meaningful (testing the three points are indeed mapped to 0, 1, oo in the UHP model). There is an honest bug and blaming the doctests means you have actually been trying to cover up a bug.

[tscrim]

comment:15

See also: https://en.wikipedia.org/wiki/M%C3%B6bius_transformation#Mapping_first_to_0,_1,_%E2%88%9E

[tscrim]

comment:37

There is an issue with the midpoints:

sage: UHP = HyperbolicPlane().UHP()
sage: g = UHP.get_geodesic(2+I,2+0.5*I)
sage: h = UHP.get_geodesic(2+0.5*I,2+I)
sage: g.midpoint()
Point in UHP 2.00000000000000 + 1.41421356237309*I
sage: h.midpoint()
Point in UHP 2.00000000000000 + 0.707106781186547*I

[sagetrac-git]

Changed commit from 8f00631 to 767a045

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

767a045

Fixing issue with UHP.midpoint() being order dependent.

[tscrim]

comment:39

This now fixes the issue with midpoint(), which was the same issue as in the perp bisector.

[slel]

comment:40

We test a bisector intersects a segment at its midpoint up to rounding error.

Testing that in terms of relative hyperbolic distances would be cleaner.

Proposed replacement for the works_both_ways test:

-            sage: def works_both_ways(a, b):
-            ....:     UHP = HyperbolicPlane().UHP()
-            ....:     g = UHP.get_geodesic(a, b)
-            ....:     h = UHP.get_geodesic(b, a)
-            ....:     p = g.perpendicular_bisector()
-            ....:     q = h.perpendicular_bisector()
-            ....:     c = lambda x: x.coordinates()
-            ....:     error_ab = norm(c(g.intersection(p)[0]) - c(g.midpoint()))
-            ....:     error_ba = norm(c(h.intersection(q)[0]) - c(h.midpoint()))
-            ....:     assert(bool(error_ab < 1e-9) and bool(error_ba < 1e-9)), (error_ab, error_ba)
-            sage: works_both_ways(1 + I, 2 + 0.5*I)
-            sage: works_both_ways(2 + I, 2 + 0.5*I)
+            sage: def bisector_gets_midpoint(a, b):
+            ....:     UHP = HyperbolicPlane().UHP()
+            ....:     g = UHP.get_geodesic(a, b)
+            ....:     p = g.perpendicular_bisector()
+            ....:     x, m = g.intersection(p)[0]), g.midpoint()
+            ....:     return bool(x.dist(m) < 1e-9 * a.dist(b))
+            sage: a, b, c = 1 + I, 2 + I, 2 + 0.5*I
+            sage: pairs = [(a, b), (b, a), (a, c), (c, a), (b, c), (c, b)]
+            sage: all(bisector_gets_midpoint(x, y) for x, y in pairs)

This also seems more readable to me. Sorry for the slow iterations.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

bd2d756

Doctest to check the distance in the hyperbolic metric.

[sagetrac-git]

Changed commit from 767a045 to bd2d756

[tscrim]

comment:42

I agree that that version is more readable. I tweaked the input points to have floating point numbers as otherwise the code wants to do a really horrible (but exact) symbolic computation. I also made the distance check absolute instead of relative since they should be the same point.

[tscrim]

comment:43

The gap between iterations is no problem. Thank you for taking the time to review this.

[slel]

comment:44

Replying to @tscrim:

made the distance check absolute instead of relative
since they should be the same point.

For short segments (say of length 1e-12) we might
prefer relative, but the current doctest uses segments
of length on the order of one anyway, so I don't mind.

I agree with you about doctesting with floating-point values.
How about this then:

-            sage: a, b, c = 1.0 + I, 2.0 + I, 2.0 + 0.5*I
-            sage: pairs = [(a, b), (b, a), (a, c), (c, a), (b, c), (c, b)]
-            sage: all(bisector_gets_midpoint(x, y) for x, y in pairs)
+            sage: c, d, e = CC(1, 1), CC(2, 1), CC(2, 0.5)
+            sage: pairs = [(c, d), (d, c), (c, e), (e, c), (d, e), (e, d)]
+            sage: all(bisector_gets_midpoint(a, b) for a, b in pairs)

In hyperbolic_isometry.py, pyflakes complains:
'sage.rings.all.RR' imported but unused.

[tscrim]

comment:46

I will make that change tomorrow.

Ah. right, I forgot to take care of that pyflakes issue...

[sagetrac-git]

Changed commit from bd2d756 to 234604b

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

234604b

Last details of #29936.

[tscrim]

comment:48

Fixed.

[tscrim]

comment:49

Green patchbot.

[slel]

comment:50

Good to see these issues solved.

[tscrim]

comment:51

Thank you for the review.

At some point, when I have more free time, I want to implement higher dimensional hyperbolic space. Unfortunately other directly-applicable-to-my-research-code (and paper writing) takes priority...

[vbraun]

Changed branch from public/geometry/fix_hyperbolic_plane-29936 to 234604b