Opened 2 years ago

Closed 19 months ago

#29936 closed defect (fixed)

Fix moebius_transform, midpoint and perpendicular_bisector

Reported by: gh-kliem Owned by:
Priority: major Milestone: sage-9.4
Component: geometry Keywords: hyperbolic geodesic
Cc: Samuel Lelièvre, Michael Orlitzky, Travis Scrimshaw, Greg Laun, Javier Honrubia González Merged in:
Authors: Travis Scrimshaw Reviewers: Samuel Lelièvre
Report Upstream: N/A Work issues:
Branch: 234604b (Commits, GitHub, GitLab) Commit: 234604b9faddb714d10b588325f258a0241304e2
Dependencies: #29962 Stopgaps:

Status badges

Description (last modified by Samuel Lelièvre)

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)

Change History (52)

comment:1 Changed 2 years ago by gh-kliem

Authors: Jonathan Kliem
Branch: public/29936
Commit: 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:

23ed583fix doctest in hyperbolic_space/hyperbolic_point
5283dc4use 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: 5283dc49042f5bfd45e2cfa92c567e6843b22f6d7b244c03bb763bde8a3b95fa2ce5c129562c3260

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

7b244c0modify doctests to the extend that they hold with fuzz

comment:4 Changed 2 years ago by gh-kliem

Description: modified (diff)
Keywords: hyperbolic geodesic hyperbolic point added
Status: newneeds_review
Summary: Fix precicion errors in geometryFix 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
Description: modified (diff)
Summary: Fix wrong precicion claims in geometry doctestsFix wrong precision claims in geometry doctests

comment:6 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/29936public/29936-reb
Commit: 7b244c03bb763bde8a3b95fa2ce5c129562c32608292b04158fcc2f3e8fed87670026cbdfb418744

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:

da1c6bestart from a "random" random seed for doctesting
b7b836dmake random seed reproducible
eedbe5edocument random_seed
998b1b9default random seed 0 for now
1d7b00edash instead of underscore for command line options
b31e2d5Merge branch 'public/29962' of git://trac.sagemath.org/sage into public/29962-reb
2f30dd9small fixes
b62f781doctests do not start from a random seed by default yet
1d99129fix merge conflict
8292b04Merge branch 'public/29962-reb2' of git://trac.sagemath.org/sage into public/29936-reb

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

Replying to mjo:

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 ; Changed 2 years ago by Travis Scrimshaw

Replying to gh-kliem:

Replying to mjo:

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

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?

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

Replying to tscrim:

Replying to gh-kliem:

Replying to mjo:

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

Version 0, edited 2 years ago by Travis Scrimshaw (next)

comment:16 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

Replying to tscrim:

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 works as well, because The points are uniformly distributed over the rectangle [-10, 10] \times [0, 10i].

Last edited 2 years ago by gh-kliem (previous) (diff)

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

Replying to tscrim:

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-rebpublic/geometry/fix_hyperbolic_plane-29936
Commit: 8292b04158fcc2f3e8fed87670026cbdfb4187442028478e2e1ffd9294b7a5b0e96626479ce7f546
Status: needs_reviewneeds_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:

2028478Fixing perp bisector when coordinates out of order.

comment:22 Changed 2 years ago by gh-kliem

Authors: Jonathan Kliem
Description: modified (diff)
Keywords: point removed
Summary: Fix wrong precision claims in geometry doctestsSign 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
Reviewers: 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

Commit: 2028478e2e1ffd9294b7a5b0e96626479ce7f546ad8102e43ff5417ea8be2ce6ced027260829a29e

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

ad8102eAddressing 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.2sage-9.3

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

Milestone: sage-9.3sage-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

Commit: ad8102e43ff5417ea8be2ce6ced027260829a29e3ead57105bf098e57e1795391d05a0d86f24520c

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

3e6410bFixing perp bisector when coordinates out of order.
d0ee202Addressing reviewer comments.
3ead571Moebius transformations should not do conjugation.

comment:32 Changed 20 months ago by Travis Scrimshaw

Status: needs_workneeds_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

Commit: 3ead57105bf098e57e1795391d05a0d86f24520c8f00631626fa636d1979842a1e0b6761f1e2d399

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

8f00631Better check for trac #29936.

comment:36 Changed 20 months ago by Travis Scrimshaw

Status: needs_reviewneeds_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: 8f00631626fa636d1979842a1e0b6761f1e2d399767a045b913ac1858d8785f2f78a3aedcc0718fe

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

767a045Fixing issue with UHP.midpoint() being order dependent.

comment:39 Changed 20 months ago by Travis Scrimshaw

Status: needs_workneeds_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: 767a045b913ac1858d8785f2f78a3aedcc0718febd2d756c73c0f6d2f9677dfec9338752fcaeb6c8

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

bd2d756Doctest to check the distance in the hyperbolic metric.

comment:42 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

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.

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

Description: modified (diff)
Summary: Sign error in hyperbolic geodesicFix 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

Commit: bd2d756c73c0f6d2f9677dfec9338752fcaeb6c8234604b9faddb714d10b588325f258a0241304e2

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

234604bLast details of #29936.

comment:48 Changed 20 months ago by Travis Scrimshaw

Fixed.

comment:49 Changed 19 months ago by Travis Scrimshaw

Green patchbot.

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

Status: needs_reviewpositive_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-29936234604b9faddb714d10b588325f258a0241304e2
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.