Revision history [back]

Hyperbolic polygons in the disc model

After noticing that there has been work towards this (though currently stalled), we explain how exploring the objects and their methods can suggest solutions.

Does Sage provide it already?

Regardless, let's explore a bit and build some solutions.

Exploration: an example

Let us explore a little starting from the code in the question (thanks for providing a starting point!).

sage: DD = HyperbolicPlane().PD()
sage: lines = [DD.get_geodesic(p, q) for p, q in [(1, I), (I, -1), (-1, 1)]]
sage: p = sum(line.plot() for line in lines)
sage: p.show(figsize=3)
Launched png viewer for Graphics object consisting of 6 graphics primitives

Notice that the boundary disc looks like it was plotted three times. It is indeed the case, as we shall see. A hint, in the Sage REPL (in the terminal) is that it mentions "6 graphics primitives", which we can guess are the three geodesics and three copies of the boundary circle.

To get started, let's give names to the geodesics.

sage: a, b, c = lines

Explore their methods using the TAB key: a.<TAB> reveals the method plot.

Give a name to the plots:

sage: pa, pb, pc = a.plot(), b.plot(), c.plot()

Explore their methods using the TAB key: pa.<TAB> reveals the method description.

sage: print(pa.description())
Arc with center (1.0,1.0) radii (1.0,1.0) angle 0.0
inside the sector (3.141592653589793,4.71238898038469)
Circle defined by (0.0,0.0) with r=1.0
sage: print(pb.description())
Arc with center (-1.0,1.0) radii (1.0,1.0) angle 0.0
inside the sector (-1.5707963267948966,0.0)
Circle defined by (0.0,0.0) with r=1.0
sage: print(pc.description())
Line defined by 2 points:   [(-1.0, 0.0), (1.0, 0.0)]
Circle defined by (0.0,0.0) with r=1.0

We see each description consists of two lines, one for an arc and a second one for the boundary circle. So indeed it was being plotted three times.

Could we get rid of the boundary circle? Exploring the plot method using a.plot? reveals the optional keyword boundary which defaults to True. Using False instead, we can get rid of the extra circles.

sage: pa_nb = a.plot(boundary=False)
sage: pb_nb = b.plot(boundary=False)
sage: pc_nb = c.plot(boundary=False)
sage: abc0 = pa_nb + pb_nb + pc_nb  # no boundary
sage: abc3 = pa + pb + pc  # boundary plotted 3 times
sage: abc1 = pa_nb + pb_nb + pc  # boundary plotted once
sage: abc_013 = graphics_array([abc0, abc3, abc1])
sage: abc_013.show(figsize=3)

Can we get the boundary by itself? By exploring DD with DD.<TAB>, we find a method get_background_graphic which gives it.

Besides the multiple plotting of the boundary, we observe that:

• geodesics through the origin are euclidean line segments; their description gives their two endpoints
• geodesics not through the origin are euclidean circle arcs; their description gives center and radius

This gives the idea of using region_plot (see the Sage graphics tour).

For that, we write inequalities to describe the region

• inside the unit circle:
• outside the circles for sides not through the origin:
• on the appropriate side of the lines for sides through the origin:

In our case:

sage: boundary = DD.get_background_graphic()
sage: p = sum(line.plot(boundary=False) for line in lines)
sage: x, y = SR.var('x, y')
sage: ieqs = [x^2 + y^2 < 1,
....:         (x - 1)^2 + (y - 1)^2 > 1,
....:         (x + 1)^2 + (y - 1)^2 > 1,
....:         y > 0]
sage: opt = {'incol': '#ffff7f', 'axes': False}
sage: inside = region_plot(ieqs, (x, -1, 1), (y, 0, 1), **opt)
sage: g = boundary + inside + p
sage: g.save('disc_triangle_filled_in.png', figsize=2)
sage: g.show(figsize=2)

Polygon plotting function using Sage's region plot function

Based on our exploration, we write the following function:

def discgon(points, color='teal', fill_color='plum',
            plot_points=65, thickness=1.5, **opt):
    r"""
    Return a hyperbolic polygon in the disc model as a graphics object.

    INPUT:

    - ``points`` -- a list of complex numbers in the unit disc,
      consecutive vertices of a convex hyperbolic polygon
    - ``color`` -- (default: 'teal') line thickness for sides
    - ``fill_color`` -- (default: 'plum') fill color for polygon
    - ``plot_points`` -- (default: 257) number of points on the arcs
    - ``thickness`` -- (default: 1.5) line thickness for sides
    - ``opt`` -- optional extra plotting options

    EXAMPLES::

        sage: a = [1, I, -1]
        sage: b = [0, 1/2, 1/2 + 1/2*I, 1/2*I]
        sage: c = [0, 1/2, 1/2 + 1/2*I, I, -1/2 + 1/2*I]
        sage: d = [(3/4 if k % 2 else 1)*E(8)^k for k in range(8)]
        sage: abcd = graphics_array([discgon(p) for p in [a, b, c, d]])
        sage: abcd.show(figsize=5)
        Launched png viewer for Graphics Array of size 1 x 4
        sage: fivegons = [discgon([r/4*E(5)^k for k in range(5)])
        ....:             for r in range(1, 5)]
        sage: fives = graphics_array(fivegons)
        sage: fives.show(figsize=5)
        Launched png viewer for Graphics Array of size 1 x 4
    """
    line_opt = {'color': color, 'thickness': thickness}
    fill_opt = {'incol': fill_color, 'plot_points': plot_points}
    n = len(points)
    if n < 3:
        raise ValueError(f"need >= 3 points for polygon, got {n}")
    xy = [(real(z), imag(z)) for z in points]
    xx, yy = zip(*(list(z) for z in xy))
    xab = min(xx), max(xx)
    yab = min(yy), max(yy)
    x, y = SR.var('x, y')
    DD = HyperbolicPlane().PD()
    relations = [x^2 + y^2 < RDF.one()]
    side_plots = []
    for k, q in enumerate(points, start=-1):
        p = points[k]
        pq = DD.get_geodesic(p, q).plot(boundary=False, **line_opt, **opt)
        side_plots.append(pq)
        d = pq.description()
        if d.startswith('Arc'):
            xc, yc, r, _ = (RDF(z) for s in d.split(')')[:-2]
                            for z in s.split('(')[1].split(','))
            relations.append((x - xc)^2 + (y - yc)^2 > r^2)
        elif d.startswith('Line'):
            xa, ya, xb, yb = (RDF(z) for s in d.split(')')[:-1]
                              for z in s.split('(')[1].split(','))
            a = yb - ya
            b = xa - xb
            lin_form = a * x + b * y
            aff_fc = lin_form - lin_form.subs({x: p.real(), y: p.imag()})
            a, r, j = max((abs(ft), ft, j)
                          for j, t in enumerate(xy)
                          for ft in [aff_fc.subs({x: t[0], y: t[1]})])
            relations.append(sign(r)*aff_fc > 0)
        else:
            raise ValueError("unexpected side type")
    inside = region_plot(relations, xab, yab, **fill_opt, **opt)
    return sum(side_plots, DD.get_background_graphic(**opt) + inside)

We test our function:

sage: a = [1, I, -1]
sage: b = [0, 1/2, 1/2 + 1/2*I, 1/2*I]
sage: c = [0, 1/2, 1/2 + 1/2*I, I, -1/2 + 1/2*I]
sage: d = [(3/4 if k % 2 else 1)*E(8)^k for k in range(8)]
sage: abcd = graphics_array([discgon(p) for p in [a, b, c, d]])
sage: abcd.show(figsize=5)
Launched png viewer for Graphics Array of size 1 x 4

sage: fives = [[r/4*E(5)^k for k in range(5)] for r in range(1, 5)]
sage: fivegons = graphics_array([discgon(f) for f in fives])
sage: fivegons.show(figsize=5)
Launched png viewer for Graphics Array of size 1 x 4

Improvement ideas

The approach above is somewhat unsatisfactory for several reasons.

First, the region plot poorly fills near acute vertices, especially if we sample too few points.

Second, sampling is wasteful in itself, since we already know the outline of the region to plot.

This suggests approximating that outline by a euclidean polygon and using Sage's 2D polygon plotting function, which can fill.

Replacing sides by polygonal lines is easy: sides through the origin are already segments, and those not through the origin are circle arcs which we can approximate by polygonal lines.

Polygon plotting function using Sage's (euclidean) 2D polygon function

Let us warm up with the initial example. The side descriptions not only include the center and radius of the corresponding euclidean circle, but also the two extreme angles for the arc.

So we can write something like:

sage: nsteps = 64
sage: s, t = (3.141592653589793, 4.71238898038469)
sage: step = (t - s) / nsteps
sage: stop = s - step / 2
sage: aa = [(1 + cos(x), 1 + sin(x)) for x in srange(t, stop, -step)]
sage: s, t = (-1.5707963267948966, 0.0)
sage: step = (t - s) / nsteps
sage: stop = s - step / 2
sage: bb = [(-1 + cos(x), 1 + sin(x)) for x in srange(t, stop, -step)]
sage: cc = [(-1.0, 0.0), (1.0, 0.0)]
sage: abc = aa + bb + cc
sage: abc_gon = polygon2d(abc, color='plum', edgecolor='teal', axes=False)
sage: abc_gon.show(figsize=2)

Turning this into a function, we get:

def disc_ngon(points, color='teal', fill_color='plum',
              thickness=1.5, nsteps=256, **opt):
    r"""
    Return a hyperbolic polygon in the disc model as a graphics object.

    INPUT:

    - ``points`` -- list of complex numbers in the unit disc for vertices
      of a convex hyperbolic polygon listed counterclockwise
    - ``color`` -- (default: 'teal') the line color for the sides
    - ``fill_color`` -- (default: 'plum') the fill color for the polygon
    - ``thickness`` -- (default: 1.5) the line thickness for the sides
    - ``nsteps`` -- (default: 256) number of steps for circular sides
    - ``opt`` -- optional extra options

    Circular sides are turned into polygonal lines with ``nsteps`` segments
    for filling the polygon.

    EXAMPLES::

        sage: a = [1, I, -1]
        sage: b = [0, 1/2, 1/2 + 1/2*I, 1/2*I]
        sage: c = [0, 1/2, 1/2 + 1/2*I, I, -1/2 + 1/2*I]
        sage: d = [(3/4 if k % 2 else 1)*E(8)^k for k in range(8)]
        sage: abcd = graphics_array([discgon(p) for p in [a, b, c, d]])
        sage: abcd.show(figsize=5)
        Launched png viewer for Graphics Array of size 1 x 4
        sage: fivegons = [discgon([r/4*E(5)^k for k in range(5)])
        ....:             for r in range(1, 5)]
        sage: fives = graphics_array(fivegons)
        sage: fives.show(figsize=5)
        Launched png viewer for Graphics Array of size 1 x 4
    """
    line_opt = {'color': color, 'thickness': thickness}
    fill_opt = {'color': fill_color, 'thickness': 0}
    n = len(points)
    if n < 3:
        raise ValueError(f"need >= 3 points for polygon, got {n}")
    xy = [(real(z), imag(z)) for z in points]
    DD = HyperbolicPlane().PD()
    side_plots = []
    vertices = []
    for k, q in enumerate(points, start=-1):
        p = points[k]
        vertices.append(xy[k])
        pq = DD.get_geodesic(p, q).plot(boundary=False, **line_opt, **opt)
        side_plots.append(pq)
        d = pq.description()
        if d.startswith('Arc'):
            xc, yc, r, _, s, t = (RDF(z) for s in d.split(')')[:-1]
                                  for z in s.split('(')[1].split(','))
            step = (t - s) / nsteps
            vertices.extend([(xc + r*cos(x), yc + r*sin(x))
                             for x in srange(t + step/2, s, -step)])
    inside = polygon2d(vertices, **fill_opt, **opt)
    return sum(side_plots, DD.get_background_graphic(**opt) + inside)

We test our function:

sage: a = [1, I, -1]
sage: b = [0, 1/2, 1/2 + 1/2*I, 1/2*I]
sage: c = [0, 1/2, 1/2 + 1/2*I, I, -1/2 + 1/2*I]
sage: d = [(3/4 if k % 2 else 1)*E(8)^k for k in range(8)]
sage: abcd = graphics_array([disc_ngon(p) for p in [a, b, c, d]])
sage: abcd.show(figsize=5)
Launched png viewer for Graphics Array of size 1 x 4

sage: fives = [[r/4*E(5)^k for k in range(5)] for r in range(1, 5)]
sage: fivegons = graphics_array([disc_ngon(f) for f in fives])
sage: fivegons.show(figsize=5)
Launched png viewer for Graphics Array of size 1 x 4