Add algorithms to compute rational pleating rays

In farey.py we now have an implementation of Newton's algorithm
and other paraphernalia to compute points on pleating rays.

There is also a new example, parabolic_slice_pleating_rays.py.

This is the first half of the features listed in #18.

CHANGELOG

@@ -1,6 +1,7 @@
 v0.1.1
- - Add accidental parabolic example and an example of (1;2)-compression body groups.
  - Add continued fractions algorithms
+ - Add algorithms to compute rational pleating rays, see examples/parabolic_slice_pleating_rays.py.
+ - Add accidental parabolic example and an example of (1;2)-compression body groups.
  - Compute Lyndon-Ullman bounds in the Riley slice explorer
 
 v0.1.0

bella/farey.py

@@ -13,6 +13,13 @@ class FractionOutOfRangeException(Exception):
     pass
 
 
+##############################################
+###
+### Farey words
+###
+##############################################
+
+
 @functools.cache
 def farey_word(r,s):
     """ Compute the Farey word of slope r/s using the cutting sequence definition.
@@ -132,6 +139,14 @@ def peripheral_splittings(w, include_conjugates = True):
         splittings.append((u, v))
     return splittings
 
+
+
+##############################################
+###
+### Farey sequences and misc functions on fractions
+###
+##############################################
+
 @functools.cache
 def next_neighbour(p,q):
     """ Return the larger Farey neighbour of p/q in the Farey sequence of denominator q.
@@ -202,6 +217,14 @@ def walk_tree_bfs(end = None):
                 yield (r,s)
 
 
+
+
+##############################################
+###
+### Farey and Riley polynomials
+###
+##############################################
+
 def numpy_to_mpmath_polynomial(np_poly):
     """ Convert a numpy polynomial to a polynomial that mpmath will accept. """
     return list(reversed(np_poly.coef))
@@ -298,6 +321,15 @@ def riley_polynomial(r,s):
     return p
 
 
+
+
+##############################################
+###
+### Continued fraction approximations
+###
+##############################################
+
+
 def euclidean_algorithm(a,b):
     """ Run the Euclidean algorithm for a/b.
 
@@ -405,3 +437,56 @@ def collapse_continued_fraction(expansion):
     expansion = expansion + [1]
     return collapse_continued_fraction(expansion)
 
+
+
+
+##############################################
+###
+### Methods to compute rational pleating rays
+###
+##############################################
+
+def newtons_method(f, z0, df = None, tol_re = 1e-10, tol_im = 1e-20):
+    """ Run Newton's method starting at z0 until reaching the given tolerance.
+
+        Real and imaginary tolerances can be given separately.
+    """
+
+    if df == None:
+        df = f.deriv()
+
+    w0 = f(z0)
+    while mp.fabs(w0.real) > tol_re or mp.fabs(w0.imag) > tol_im:
+        z0 = z0 - w0/df(z0)
+        w0 = f(z0)
+
+    return z0
+
+def real_point_on_circle(f, R, angle, tol=None):
+    """ Find the point on the circle |z| = R where Im f(z) = 0 with arg closest to angle.
+    """
+    actual_function = lambda theta: ( f(R*mp.exp(1j * theta)) ).imag
+    return R*mp.exp(1j*mp.findroot(actual_function, angle, tol=tol))
+
+def approximate_pleating_ray(r, s, p, q, R = 20, N = 10):
+    """ Return points on the r/s pleating ray of the (p,q)-Riley slice.
+
+        We will return N+1 points [z0,...,zN] where z0 is a point on the r/s pleating ray
+        with |z0| = R and where zN is the r/s cusp point.
+
+        Arguments:
+          r,s -- coprime integers representing the slope
+          p,q -- orders of the group generators, so we are working in R^{p,q}
+          R -- radius (in CC) to start approximating at
+          N -- number of points to compute (minus 1)
+    """
+    f = farey_polynomial_classic(r,s,p,q)
+    df = f.deriv()
+    z = [real_point_on_circle(f, R, mp.pi*r/s)]
+    L = [f(z[0])]
+    epsilon = mp.fabs(L[0] + 2)/N
+    for i in range(1,N+1):
+        L.append( (1 - i/N) * L[-1] - 2*(i/N) )
+        z.append(newtons_method(f - L[-1], z[-1], df))
+
+    return z

examples/README.md

@@ -3,7 +3,8 @@
 ## Basic examples
 Run these with plain vanilla `python [filename]`. The utility bash script `generate_zoo.sh` will run all of the examples that just produce PNG files.
  - [cayley_graph_speed.py](cayley_graph_speed.py) -- demonstrate that walking the Cayley group of a non-free group incurs a massive speed overhead.
- - [parabolic_slice_hidef.py](parabolic_slice_hidef.py) -- draw a PNG file containing a high-res picture of the parabolic Riley slice to high definition.
+ - [parabolic_slice_hidef.py](parabolic_slice_hidef.py) -- draw a PNG file containing a high-res picture of the parabolic Riley slice.
+ - [parabolic_slice_pleating_rays.py](parabolic_slice_pleating_rays.py) -- plot rational pleating rays around the parabolic Riley slice.
  - [apollonian_gasket.py](apollonian_gasket.py) -- draw a PNG file containing a picture of the Apollonian Gasket.
  - [padic.py](padic.py) -- demonstrate that it is possible to use other number types (e.g. $` p `$-adic numbers) as long as you are sufficiently masochistic.
  - [farey.py](farey.py) -- produce TeX tables of Farey and Riley words and polynomials.

examples/generate_zoo.sh

@@ -1,6 +1,6 @@
 #!/bin/sh
 
-image_scripts="parabolic_slice_hidef.py apollonian_gasket.py padic.py jorgensen_marden.py geometrically_infinite.py modular_group.py connected_component_bound.py apanasov.py web.py elementary.py beads.py zarrow.py accidental_parabolic.py atom.py riley_slice_cusps.py"
+image_scripts="parabolic_slice_hidef.py parabolic_slice_pleating_rays.py apollonian_gasket.py padic.py jorgensen_marden.py geometrically_infinite.py modular_group.py connected_component_bound.py apanasov.py web.py elementary.py beads.py zarrow.py accidental_parabolic.py atom.py riley_slice_cusps.py"
 
 echo "generate_zoo.sh: generating all basic examples that just produce images."
 for i in $image_scripts

examples/parabolic_slice_pleating_rays.py

@@ -0,0 +1,23 @@
+""" Example: plotting the parabolic Riley slice at high definition, as well as various pleating rays.
+"""
+
+from bella import slices, farey
+from mpmath import mp
+import holoviews as hv
+hv.extension('bokeh')
+
+# Just compute the parabolic slice and save it to an image file.
+depth = 50
+df = slices.parabolic_exterior_from_farey(depth)
+
+rays = []
+for r,s in farey.walk_tree_bfs(6):
+  rays.append(farey.approximate_pleating_ray(r,s,mp.inf,mp.inf, R=20, N=100))
+
+roots = hv.Scatter(df, kdims=['x'],vdims=['y']).opts(marker = "dot", size = 4, frame_width=1600, frame_height=800, data_aspect=1, color='gray')\
+          .redim(x=hv.Dimension('x', range=(-5,5)),y=hv.Dimension('y', range=(-2.5, 2.5)))
+
+for ray in rays:
+    roots *= hv.Path([(float(z.real), float(z.imag)) for z in ray]).opts(color='black')
+
+hv.save(roots, 'parabolic_slice_pleating_rays.png', fmt='png')