Opened 2 years ago

Closed 19 months ago

# Fix moebius_transform, midpoint and perpendicular_bisector

Reported by: Owned by: gh-kliem major sage-9.4 geometry hyperbolic geodesic Samuel Lelièvre, Michael Orlitzky, Travis Scrimshaw, Greg Laun, Javier Honrubia González Travis Scrimshaw Samuel Lelièvre N/A 234604b 234604b9faddb714d10b588325f258a0241304e2 #29962

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)


### comment:1 Changed 2 years ago by gh-kliem

Authors: → Jonathan Kliem → public/29936 → 5283dc49042f5bfd45e2cfa92c567e6843b22f6d

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

### comment:2 Changed 2 years ago by gh-kliem

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


### comment:3 Changed 2 years ago by git

Commit: 5283dc49042f5bfd45e2cfa92c567e6843b22f6d → 7b244c03bb763bde8a3b95fa2ce5c129562c3260

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

 ​7b244c0 modify doctests to the extend that they hold with fuzz

### comment:4 Changed 2 years ago by gh-kliem

Description: modified (diff) hyperbolic geodesic hyperbolic point added new → needs_review Fix precicion errors in geometry → Fix wrong precicion claims in geometry doctests

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

### comment:5 Changed 2 years ago by Samuel Lelièvre

Cc: Samuel Lelièvre added modified (diff) Fix wrong precicion claims in geometry doctests → Fix wrong precision claims in geometry doctests

### comment:6 follow-up:  9 Changed 2 years ago by Michael Orlitzky

Cc: Michael Orlitzky Travis Scrimshaw Greg Laun Javier Honrubia González added

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?

### comment:7 Changed 2 years ago by gh-kliem

Dependencies: → #29962

### comment:8 Changed 2 years ago by gh-kliem

Branch: public/29936 → public/29936-reb 7b244c03bb763bde8a3b95fa2ce5c129562c3260 → 8292b04158fcc2f3e8fed87670026cbdfb418744

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

### comment:9 in reply to:  6 ; follow-up:  10 Changed 2 years ago by gh-kliem

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

### comment:10 in reply to:  9 ; follow-ups:  11  12 Changed 2 years ago by Travis Scrimshaw

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.

### comment:11 in reply to:  10 Changed 2 years ago by Michael Orlitzky

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?

### comment:12 in reply to:  10 Changed 2 years ago by gh-kliem

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


### comment:13 Changed 2 years ago by gh-kliem

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.

### comment:14 follow-up:  20 Changed 2 years ago by Travis Scrimshaw

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.

Last edited 2 years ago by Travis Scrimshaw (previous) (diff)

### comment:16 follow-up:  18 Changed 2 years ago by Travis Scrimshaw

Note that removing the conjugate() makes the test RP*gP == gP in HyperbolicGeodesic.reflection_involution() fail.

I still don't know what is going on with the perpendicular bisector.

### comment:17 Changed 2 years ago by Travis Scrimshaw

With the bisector, sometimes it is not intersecting the geodesic (but the completed versions do intersect perpendicularly). So that is why the test is failing, but I don't know why it is doing this.

### comment:18 in reply to:  16 Changed 2 years ago by gh-kliem

Note that removing the conjugate() makes the test RP*gP == gP in HyperbolicGeodesic.reflection_involution() fail.

I still don't know what is going on with the perpendicular bisector.

You're right. It is a sign issue and not a precision issue:

....: 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)
....:     a1 = abs(moebius_transform(A, p1))
....:     a2 = abs(moebius_transform(A, p1.conjugate()))
....:     a = min(a1, a2)
....:     b1 = abs((moebius_transform(A, p2) - 1))
....:     b2 = abs((moebius_transform(A, p2.conjugate()) - 1))
....:     b = min(b1, b2)
....:     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))
3.32270874593069e-17
sage: max(foo()[0] for _ in range(1000))
0.000000000000000


If you return just a, b (the absolute errors), you get:

sage: max(foo()[1] for _ in range(1000))
2.79547014692442e-15
sage: max(foo()[0] for _ in range(1000))
0.000000000000000


I guess the absolute error is well again, because The points are uniformly distributed over the rectangle [-10, 10] \times [0, 10i].

Version 0, edited 2 years ago by gh-kliem (next)

### comment:19 Changed 2 years ago by gh-kliem

So the original doctests are fine up to a sign issue.

### comment:20 in reply to:  14 Changed 2 years ago by gh-kliem

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.

I opened #29939 to fix this problem. But if we can fix this here that is fine with me to. If it would have been a precision problem, I don't think I would have been able to fix this. So this is why I opened a separate ticket.

Anyway, I'm glad you figured out the actual problem.

### comment:21 Changed 2 years ago by Travis Scrimshaw

Branch: public/29936-reb → public/geometry/fix_hyperbolic_plane-29936 8292b04158fcc2f3e8fed87670026cbdfb418744 → 2028478e2e1ffd9294b7a5b0e96626479ce7f546 needs_review → needs_work

So the issue with the perp bisector comes from this:

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)


In particular, the completed geodesics always have the endpoints in order. I am not sure if this is the best fix, but it does fix that issue.

However, I don't yet understand enough of the mathematics to fix the issue with the crossratio matrix.

New commits:

 ​2028478 Fixing perp bisector when coordinates out of order.

### comment:22 Changed 2 years ago by gh-kliem

Authors: Jonathan Kliem modified (diff) point removed Fix wrong precision claims in geometry doctests → Sign error in hyperbolic geodesic

### comment:23 Changed 2 years ago by Samuel Lelièvre

I would rewrite the test along these lines:

+        Check the result is independent of the order::
+
+            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()
+            ....:     assert bool(c(g.intersection(p)[0]) - c(g.midpoint()) < 1e-9)
+            ....:     assert bool(c(h.intersection(q)[0]) - c(g.midpoint()) < 1e-9)
+            sage: works_both_ways(1 + I, 2 + 0.5*I)
+            sage: works_both_ways(2 + I, 2 + 0.5*I)


### comment:24 Changed 2 years ago by Travis Scrimshaw

That is a better way to do the test. I will change that next chance I get (although it might not be until the end of the week).

### comment:25 Changed 2 years ago by Samuel Lelièvre

Authors: → Travis Scrimshaw → Samuel Lelièvre

Some more suggestions.

Add trac ticket number to test:

+        Check the result is independent of the order (:trac:29936)::


Rephrase comment?

-        # Since the complete geodesic p1 -> p2 always returns p1 < p2,
-        #   we need to swap start and end when that happens
+        # Completed geodesics always have endpoints in order,
+        # so we swap start and end if needed.


Define S immediately after the swap, since that's where .complete() happens.

+        S = self.complete()._to_std_geod(start)
d = self._model._dist_points(start, self._end.coordinates()) / 2
-        S = self.complete()._to_std_geod(start)


Avoid detour to the symbolic ring.

-        s2 = sqrt(2) * 0.5
+        s2 = sqrt(0.5)


### comment:26 Changed 2 years ago by git

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

 ​ad8102e Addressing reviewer comments.

### comment:27 Changed 2 years ago by Travis Scrimshaw

I managed to find a bit of time today. Here are the changes you requested although I used a different version of the comment.

Now to figure out how to fix the issue in comment:14 as some parts of the code that use moebius_transformation seem to use that conjugation...

### comment:28 Changed 2 years ago by Matthias Köppe

Milestone: sage-9.2 → sage-9.3

### comment:29 Changed 22 months ago by Matthias Köppe

Milestone: sage-9.3 → sage-9.4

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

### comment:30 Changed 22 months ago by Samuel Lelièvre

The patchbot's pyflakes plugin complains that moebius_transform is undefined.

### comment:31 Changed 20 months ago by git

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

 ​3e6410b Fixing perp bisector when coordinates out of order. ​d0ee202 Addressing reviewer comments. ​3ead571 Moebius transformations should not do conjugation.

### comment:32 Changed 20 months ago by Travis Scrimshaw

Status: needs_work → needs_review

Here is a rebased branch that fixes the other problem. It is indeed related with the moebius_tranformation. That is really meant for isomorphisms of CP1 with maps in PGL(2, C), so we should not do complex conjugation. Moreover, we cannot expect the determinant to be a real number unless we choose a good representative, but it is more numerically stable to not do more divisions than necessary. The orientation preserving maps of H2 in the UHP model are for PGL(2, R), so there is a realness assumption of the determinant that we can apply, which means we can compare with 0. Because of this, I moved the orientation reversing code directly into the model computations. This should now fix all of the issues noted above.

### comment:33 Changed 20 months ago by Travis Scrimshaw

We also need the detour to the symbolic ring because sometimes a computation could be wanted in the symbolic ring and it makes checking for valid isometries much more stable.

### comment:34 Changed 20 months ago by Samuel Lelièvre

I missed one thing: to check a complex number is small, use its norm:

             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()
-            ....:     assert bool(c(g.intersection(p)[0]) - c(g.midpoint()) < 1e-9)
-            ....:     assert bool(c(h.intersection(q)[0]) - c(g.midpoint()) < 1e-9)
+            ....:     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))


Indeed, c(g.intersection(p)[0]) - c(g.midpoint()) is a complex number, and:

sage: I < 1e-9
True

Last edited 20 months ago by Samuel Lelièvre (previous) (diff)

### comment:35 Changed 20 months ago by git

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

 ​8f00631 Better check for trac #29936.

### comment:36 Changed 20 months ago by Travis Scrimshaw

Status: needs_review → needs_work

Good point, which uncovered that there is still an error in the second test...

### comment:37 Changed 20 months ago by Travis Scrimshaw

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


### comment:38 Changed 20 months ago by git

Commit: 8f00631626fa636d1979842a1e0b6761f1e2d399 → 767a045b913ac1858d8785f2f78a3aedcc0718fe

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

 ​767a045 Fixing issue with UHP.midpoint() being order dependent.

### comment:39 Changed 20 months ago by Travis Scrimshaw

Status: needs_work → needs_review

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

### comment:40 Changed 20 months ago by Samuel Lelièvre

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.

### comment:41 Changed 20 months ago by git

Commit: 767a045b913ac1858d8785f2f78a3aedcc0718fe → bd2d756c73c0f6d2f9677dfec9338752fcaeb6c8

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

 ​bd2d756 Doctest to check the distance in the hyperbolic metric.

### comment:42 follow-up:  44 Changed 20 months ago by Travis Scrimshaw

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.

### comment:43 Changed 20 months ago by Travis Scrimshaw

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

### comment:44 in reply to:  42 Changed 20 months ago by Samuel Lelièvre

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.

-            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.

### comment:45 Changed 20 months ago by Samuel Lelièvre

Description: modified (diff) Sign error in hyperbolic geodesic → Fix moebius_transform, midpoint and perpendicular_bisector

### comment:46 Changed 20 months ago by Travis Scrimshaw

I will make that change tomorrow.

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

### comment:47 Changed 20 months ago by git

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

 ​234604b Last details of #29936.

Fixed.

Green patchbot.

### comment:50 Changed 19 months ago by Samuel Lelièvre

Status: needs_review → positive_review

Good to see these issues solved.

### comment:51 Changed 19 months ago by Travis Scrimshaw

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...

### comment:52 Changed 19 months ago by Volker Braun

Branch: public/geometry/fix_hyperbolic_plane-29936 → 234604b9faddb714d10b588325f258a0241304e2 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.