Fix moebius_transform, midpoint and perpendicular_bisector
Closed this issue · 71 comments
We fix several related issues in hyperbolic geometry:
moebius_transformwrong in the orientation-reversing
(negative determinant) casemidpointandperpendicular_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
Branch: public/29936
Author: Jonathan Kliem
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 |
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
Branch pushed to git repo; I updated commit sha1. New commits:
7b244c0 | modify doctests to the extend that they hold with fuzz |
I tested 10 times and not a single test failed. I hope that is good enough.
Description changed:
---
+++
@@ -2,4 +2,7 @@
Also there is an equality in `geometry/hyperbolic_space/hyperbolic_point.py` that doesn't always hold.
-We fix those tests to also hold at random state.
+We modify precision to an extend that those tests do pass usually.
+
+Follow up:
+ #29939 the problem.Changed keywords from none to hyperbolic geodesic, hyperbolic point
Description changed:
---
+++
@@ -1,8 +1,8 @@
-There are some claims in doctests of `geometry/hyperbolic_space/hyperbolic_geodesic.py` about precicion that simply don't hold unless you get lucky.
+Some precision claims in doctests in `geometry/hyperbolic_space/hyperbolic_geodesic.py` only hold in lucky cases.
-Also there is an equality in `geometry/hyperbolic_space/hyperbolic_point.py` that doesn't always hold.
+Likewise, an equality in `geometry/hyperbolic_space/hyperbolic_point.py` doesn't always hold.
-We modify precision to an extend that those tests do pass usually.
+We add tolerance so those tests pass usually.
Follow up:
- #29939 the problem.
+- #29939 fix the precision problemI'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?
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 |
Changed branch from public/29936 to public/29936-reb
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 tolto 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).
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 tolto 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.
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?
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 tolto 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
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.
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 = zAt 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.
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.
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.
Replying to @tscrim:
Note that removing the
conjugate()makes the testRP*gP == gPinHyperbolicGeodesic.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].
So the original doctests are fine up to a sign issue.
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) +InfinitySo I think the issue comes from this in the
moebius_transformfunction:if a * d - b * c < 0: w = z.conjugate() # Reverses orientation else: w = zAt least, I only seem to notice it when the determinant of
Ahas a negative real value. So if I remove theconjugate(), 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, ooin 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.
Changed branch from public/29936-reb to public/geometry/fix_hyperbolic_plane-29936
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. |
Description changed:
---
+++
@@ -1,8 +1,13 @@
-Some precision claims in doctests in `geometry/hyperbolic_space/hyperbolic_geodesic.py` only hold in lucky cases.
+Fuzzing doctests for `geometry/hyperbolic_space/hyperbolic_geodesic.py` revealed a sign error.
-Likewise, an equality in `geometry/hyperbolic_space/hyperbolic_point.py` doesn't always hold.
-
-We add tolerance so those tests pass usually.
-
-Follow up:
-- #29939 fix the precision problem
+```
+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)
+```Changed keywords from hyperbolic geodesic, hyperbolic point to hyperbolic geodesic
Changed author from Jonathan Kliem to none
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)
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).
Reviewer: 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)
Author: Travis Scrimshaw
Branch pushed to git repo; I updated commit sha1. New commits:
ad8102e | Addressing reviewer comments. |
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...
Setting new milestone based on a cursory review of ticket status, priority, and last modification date.
The patchbot's pyflakes plugin complains that moebius_transform is undefined.
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.
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.
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
Branch pushed to git repo; I updated commit sha1. New commits:
8f00631 | Better check for trac #29936. |
Good point, which uncovered that there is still an error in the second test...
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
Branch pushed to git repo; I updated commit sha1. New commits:
767a045 | Fixing issue with UHP.midpoint() being order dependent. |
This now fixes the issue with midpoint(), which was the same issue as in the perp bisector.
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.
Branch pushed to git repo; I updated commit sha1. New commits:
bd2d756 | Doctest to check the distance in the hyperbolic metric. |
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.
The gap between iterations is no problem. Thank you for taking the time to review this.
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.
Description changed:
---
+++
@@ -1,4 +1,13 @@
-Fuzzing doctests for `geometry/hyperbolic_space/hyperbolic_geodesic.py` revealed a sign error.
+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()
@@ -11,3 +20,4 @@
(Point in UHP -4.92591958004554 + 7.81075974592370*I,
Point in UHP -8.44416583993011 + 3.87025183089797*I)
```
+I will make that change tomorrow.
Ah. right, I forgot to take care of that pyflakes issue...
Branch pushed to git repo; I updated commit sha1. New commits:
234604b | Last details of #29936. |
Fixed.
Green patchbot.
Good to see these issues solved.
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...
Changed branch from public/geometry/fix_hyperbolic_plane-29936 to 234604b