HyperbolicGeodesicPD.plot can be wrong for CDF coordinates

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

This illustrates the wrong plotting:

sage: D = HyperbolicPlane().PD()
sage: a = 0
sage: b = 1e-15 + 0.7*I
sage: geodesics = [
....:     D.get_geodesic(a, b),
....:     D.get_geodesic(CC(a), CC(b)),
....:     D.get_geodesic(CDF(a), CDF(b))
....: ]
sage: show(graphics_array([geod.plot() for geod in geodesics]))

sage: print(geodesics[2].plot()[0])
Arc
  with ​center (-0.4210526315789469,1.0000000000000007)
 ​radii (0.42105263157894823,0.42105263157894823)
 ​angle 0.0
 ​inside the sector (-1.172273881128477,-0.6190660545826288)

The first two pictures show the expected output: a straight line from 0 to 0.7*I. The last picture shows the buggy output: a short arc between two points near 0.7*I. The printed output describes the arc in the buggy picture.

CC: @slel @sagetrac-jhonrubia6 @orlitzky

Component: geometry

Keywords: hyperbolic, geodesic, plot, precision

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

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

Output image from example code

[slel]

comment:1

Attachment: geodesic_bug.png

The current implementation of hyperbolic geodesics relies on a conversion to the half-plane model.

Any geodesic g has its corresponding half-plane model geodesic stored as g._cached_geodesic.

The example geodesic has endpoints I and 2.222222222222222e-14 + 5.666666666666666*I.

The corresponding complete geodesic is a euclidean semicircle. Due to rounding errors on its centre and radius, its endpoints are computed in CDF as (-2.375, 1400000000000002.0), whereas computations in CC give them much more correctly as (0.0, 1.4e15).

Converting back to the disc model, the effect is a disaster.

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

comment:2

Aha! I suspect it's possible to find the ideal endpoints of a disk model geodesic in a more numerically stable way. I can work on a test implementation.

If the test implementation works well, do you know what we'd have to do to incorporate it into hyperbolic_geodesic.py? Is it just a matter of overriding ideal_endpoints in HyperbolicGeodesicPD?

[slel]

Description changed:

--- 
+++ 
@@ -1,22 +1,25 @@
-**Input**
+This illustrates the wrong plotting:

 ```python
-D = HyperbolicPlane().PD()
-a = 0
-b = 1e-15 + 0.7*I
-geodesics = [
-  D.get_geodesic(a, b),
-  D.get_geodesic(CC(a), CC(b)),
-  D.get_geodesic(CDF(a), CDF(b))
-]
-show(graphics_array([geod.plot() for geod in geodesics]))
-print(geodesics[2].plot()[0])
+sage: D = HyperbolicPlane().PD()
+sage: a = 0
+sage: b = 1e-15 + 0.7*I
+sage: geodesics = [
+....:     D.get_geodesic(a, b),
+....:     D.get_geodesic(CC(a), CC(b)),
+....:     D.get_geodesic(CDF(a), CDF(b))
+....: ]
+sage: show(graphics_array([geod.plot() for geod in geodesics]))
+```
+![](https://github.com/sagemath/sage/files/ticket32362/b29d3ec9c6d948fa3a44f233c2c04338.png​)
+
+```
+sage: print(geodesics[2].plot()[0])
+Arc
+  with ​center (-0.4210526315789469,1.0000000000000007)
+ ​radii (0.42105263157894823,0.42105263157894823)
+ ​angle 0.0
+ ​inside the sector (-1.172273881128477,-0.6190660545826288)

-Result

-The first two pictures show the expected output: a straight line from 0 to 0.7*I. The last picture shows the buggy output: a short arc between two points near 0.7*I. The printed output describes the arc in the buggy picture:

- -Arc with center (-0.4210526315789469,1.0000000000000007) radii (0.42105263157894823,0.42105263157894823) angle 0.0 inside the sector (-1.172273881128477,-0.6190660545826288) - +The first two pictures show the expected output: a straight line from 0 to 0.7*I. The last picture shows the buggy output: a short arc between two points near 0.7*I. The printed output describes the arc in the buggy picture.

[slel]

comment:3

Replying to @Vectornaut:

Aha! I suspect it's possible to find the ideal endpoints of a disk model geodesic in a more numerically stable way. I can work on a test implementation.

Good!

If the test implementation works well, do you know what we'd have to do to incorporate it into hyperbolic_geodesic.py? Is it just a matter of overriding ideal_endpoints in HyperbolicGeodesicPD?

Yes.

[slel]

Changed keywords from none to hyperbolic, geodesic, plot, precision

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

Finding ideal endpoints of geodesics on the Poincaré disk: example implementation

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

Attachment: hyperbolic_geodesic_PD.sage.gz

Attachment: drawing-geodesics.pdf.gz

Finding ideal endpoints of geodesics on the Poincaré disk: method explanation

[f0d43079-b341-445a-b7cd-9fcda0a76ae7]

comment:4

I just attached an example implementation. I've done some basic tests, but no systematic challenges or comparison with the current implementation (via the upper half-plane).

[7e26d9b1-e225-43a9-ae35-f759ac5662bb]

comment:6

Maybe unrelated or not (If you think it's unrelated I can open a new ticket on this) I've just discovered working on ticket [ticket:23427] that ideal endpoints on PD geodesics can be mapped to imag()<0 points which is an error, and when plotted often draw a geodesic int the lower half plane. For example

Sage:    P = PD.get_point(-0.920704875684763 - 0.390259569889459*I)
Sage:    P
Boundary point in PD -0.920704875684763 - 0.390259569889459*I
Sage:    UHP(P)
Boundary point in UHP -0.662253938491478 - 1.66533453693773e-16*I

sage: g=HyperbolicPlane().PD().random_geodesic()
Sage:  g
Geodesic in PD from -0.920704875684763 - 0.390259569889459*I to 0.0632992206159446 + 0.997994593507106*I
sage: g.to_model('UHP')  
Geodesic in UHP from -0.662253938491478 - 5.55111512312578e-17*I to 31.5642842686728 - 8.34887714518118e-14*I

[hgranath]

Not being aware of this discussion I reported a similar issue here.

The one line change I mentioned there seems to me to fix the problem both in my use cases and in the case reported in the top post above. Maybe avoiding testing of exact equality of floting point numbers in the UHP code would be a good thing to do to improve the situation?

[hgranath]

I found a case with similar symptoms, but which I think reflects a separate bug of the same type:

r = exp((pi*I/2).n())
G = hyperbolic_triangle(0, r, r*(1 + I)/2, model='PD', fill=True, alpha=.3)
G._show_axes = False
show(G)

(The central angle should be pi/4.) I think this is the root cause

PD = HyperbolicPlane().PD()
UHP = HyperbolicPlane().UHP()
r = exp((pi*I/2).n())
p = PD.get_point(r)
UHP(p)

which yields Boundary point in UHP 3.26624787063907e16 - 1.00000000000000*I A fix might be to change this line to

    if abs(x - I) < EPSILON:

(and add from sage.geometry.hyperbolic_space.hyperbolic_constants import EPSILON to that file's preamble).