Speed up limit set computations

Replace `coloured_limit_set_mc` with `coloured_limit_set_fast` that
doesn't walk the Cayley graph but just moves a single point around.

CHANGELOG

@@ -2,6 +2,7 @@ v0.0.99
  - a myriad of new examples
  - new method in cayley to construct groups by taking orientation-preserving half of group generated by circle inversions
  - improved test coverage
+ - faster limit set calculation
 
 v0.0.98.post1
  - delete ununused code in riley.py that depends on `deprecated` package not in dependencies

bella/cayley.py

@@ -5,7 +5,6 @@
 """
 
 from mpmath import mp
-from numpy.linalg import inv, eig
 import itertools
 import functools
 import random
@@ -215,11 +214,11 @@ def cayley_graph_mc(self, depth, count, yield_shorter=True):
                 if yield_shorter or n == depth:
                     yield word
 
-    def coloured_limit_set_mc(self, depth, count, seed = 0):
+    def coloured_limit_set_fast(self, count, seed=0):
         """ Monte-carlo search for points in the limit set.
 
-            Produce `depth`*count` translates of the element `seed`, thus approximating the limit set,
-            by computing the Cayley graph as returned by `cayley_graph_mc(depth, count)`.
+            Produce `count` translates of the element `seed`, thus approximating the limit set,
+            by doing a random walk.
 
             Generates: a dataframe with columns [ x, y, colour ] where x+yi is a point in the limit set
             and colour is the index of the first element in the word indexing that limit set.
@@ -228,15 +227,16 @@ def coloured_limit_set_mc(self, depth, count, seed = 0):
             base = self._underlying_matrix_t([[1],[0]])
         else:
             base = self._underlying_matrix_t([[seed],[1]])
+        def _internal_generator(base):
+            last = next(self.free_cayley_graph_mc(1,1))
+            for _ in range(count):
+                base = self[last] @ base
+                if base[1] != 0:
+                    cpx = base[0]/base[1]
+                    yield (float(cpx.real), float(cpx.imag), last[0])
+                last = self.free_random_walk_locally(last)[:1]
 
-        def _internal_generator():
-            for w in self.free_cayley_graph_mc(depth,count):
-                point = self[w] @ base
-                if point[1] != 0:
-                    cpx = point[0]/point[1]
-                    yield (float(cpx.real), float(cpx.imag), w[0])
-
-        return pd.DataFrame(_internal_generator(), columns=['x','y','colour'])
+        return pd.DataFrame(_internal_generator(base), columns=['x','y','colour'])
 
     def isometric_circle(self, word):
         """ Return the isometric circle of the word.
@@ -253,7 +253,7 @@ def isometric_circle(self, word):
     def coloured_isometric_circles_mc(self, depth, count):
         """ Monte-carlo search for isometric circles in the limit set.
 
-            Produce `depth`*count` isometric circles, thus approximating the limit set,
+            Produce `depth`*`count` isometric circles, thus approximating the limit set,
             by computing the Cayley graph as returned by `cayley_graph_mc(depth, count)`.
 
             Generates: a dataframe with columns [ x, y, radius, colour ] where (x,y) is the centre
@@ -263,7 +263,7 @@ def coloured_isometric_circles_mc(self, depth, count):
 
         def _internal_generator():
             for w in self.cayley_graph_mc(depth,count):
-                centre, radius = isometric_circle(w)
+                centre, radius = self.isometric_circle(w)
                 if centre == mp.inf:
                     continue
                 yield [float(centre.real), float(centre.imag), float(radius), w[0]]

examples/animate_limits_as_slices_vary.py

@@ -9,11 +9,11 @@
 import os
 
 # Output one frame of the animation. See below for explanation of the parameters.
-def one_frame(kk,θ,depth,logpoints,first,scale):
+def one_frame(kk,θ,num_points,first,scale):
     print(f"#{kk}")
     η = (kk/scale) * mp.pi
     G = riley.RileyGroup(θ,η,2j) # <- we fix the value of μ.
-    df = G.coloured_limit_set_mc(depth,10**logpoints)
+    df = G.coloured_limit_set_fast(num_points)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 1,  color = 'colour', width=1000, height=1000, data_aspect=1, cmap='Category10')\
                 .redim(x=hv.Dimension('x', range=(-1,1)),y=hv.Dimension('y', range=(-1, 1)))
@@ -23,10 +23,9 @@ def one_frame(kk,θ,depth,logpoints,first,scale):
 # We use this __main__ pattern since we have to use multiprocessing, see
 # the documentation: https://docs.python.org/dev/library/multiprocessing.html#multiprocessing-programming
 if __name__=='__main__':
-    # θ is fixed. Depth and logpoints are passed directly into GroupCache.coloured_limit_set_mc().
+    # θ is fixed. num_points is passed directly into GroupCache.coloured_limit_set_fast().
     θ = 0
-    depth = 20
-    logpoints = 4
+    num_points = 20*10**4
 
     # We will compute η as (kk/scale)*π, where kk runs from first to last. Hence to adjust the number of frames
     # adjust scale and modify last so that last/scale is the endpoint value of the anumation that you want.
@@ -45,4 +44,4 @@ def one_frame(kk,θ,depth,logpoints,first,scale):
     # module, the only way we get out of running out of memory is to start a new process each time.) As a nice side-effect, it's fast.
     multiprocessing.set_start_method('spawn') # Using "fork" also leaks memory.
     with multiprocessing.Pool(4, maxtasksperchild=1) as pool:
-        pool.starmap(one_frame, [ (kk,θ,depth,logpoints,first,scale) for kk in range(first, last) ], chunksize=1)
+        pool.starmap(one_frame, [ (kk,θ,num_points,first,scale) for kk in range(first, last) ], chunksize=1)

examples/apanasov.py

@@ -29,11 +29,10 @@ def __init__(self):
         print(f"det Y = {cayley.simple_det(Y)}")
         super().__init__([X,Y])
 
-depth = 30
-logpoints = 5
+num_points = 30*10**5
 G = ApanasovGroup()
 seed = G.fixed_points((0,1))[0]
-df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+df = G.coloured_limit_set_fast(num_points, seed=seed)
 scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
             .opts(marker = "dot", size = 1,  color = 'colour', width=1000, height=1000, data_aspect=1, cmap='Set1')
 

examples/apollonian_gasket.py

@@ -10,7 +10,7 @@
 logpoints = 4
 G = riley.RileyGroup(0, 0, 2j)
 seed = G.fixed_points((0,1))[0]
-df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+df = G.coloured_limit_set_fast(depth*10**logpoints, seed=seed)
 scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
             .opts(marker = "dot", size = 0.1,  color = 'colour', width=1600, height=1600, data_aspect=1, cmap='Set1')\
               .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))

examples/beads.py

@@ -53,8 +53,7 @@ def __init__(self, path):
 
         super().__init__(cayley.generators_from_circle_inversions(self.circles, []))
 
-depth = 50
-logpoints = 4
+num_points = 10**5
 
 def write_limit_set(G,filename):
     # We put fixed points and beads onto the picture, here is the calculation.
@@ -63,7 +62,7 @@ def write_limit_set(G,filename):
     fixed_points = [hv.Points([ [float(p.real), float(p.imag)] for p in G.fixed_points((n,))]).opts(marker = "dot", size = 20) for n in range(0,len(G))]
 
     # compute the limit set
-    df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+    df = G.coloured_limit_set_fast(num_points, seed=seed)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 0.1,  color = 'colour', width=1600, height=1600, data_aspect=1, cmap='Set1')
 

examples/connected_component_bound.py

@@ -18,11 +18,10 @@ def __init__(self, N):
             gens.append(mp.matrix([[1+2j*m, 4*m**2],[1,1-2j*m]]))
         super().__init__(gens)
 
-depth = 20
-logpoints = 4
+num_points = 10**5
 G = BoundGroup(8)
 seed = G.fixed_points((0,1))[0]
-df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+df = G.coloured_limit_set_fast(num_points, seed=seed)
 scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
             .opts(marker = "dot", size = 1,  color = 'colour', width=400, height=1500, data_aspect=1, cmap='Set1')\
               .redim(x=hv.Dimension('x', range=(-2,2)),y=hv.Dimension('y', range=(-1, 14)))

examples/elementary.py

@@ -19,13 +19,12 @@ def __init__(self, λ,μ):
         self.lines = [ (0,λ), (0,μ), (λ,λ+μ), (μ,λ+μ) ]
         super().__init__(cayley.generators_from_circle_inversions([], self.lines))
 
-depth = 30
-logpoints = 3
+num_points = 10**4
 
 # Roots of unity
 G = IncompleteGroup(1, 1+1j) #<- if you modify this so that the parallelogram has angles not submultiples of unity, the resulting group is no longer discrete.
 seed = 0
-df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+df = G.coloured_limit_set_fast(num_points, seed=seed)
 scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
             .opts(marker = "dot", size = 4,  color = 'colour', width=800, height=800, data_aspect=1, cmap='Set1')\
               .redim(x=hv.Dimension('x', range=(-8,8)),y=hv.Dimension('y', range=(-8, 8)))

examples/generate_zoo.sh

@@ -1,6 +1,7 @@
 #!/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 atom.py zarrow.py"
+image_scripts="jorgensen_marden.py geometrically_infinite.py modular_group.py connected_component_bound.py apanasov.py web.py elementary.py beads.py atom.py zarrow.py"
+# 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 atom.py zarrow.py"
 
 echo "generate_zoo.sh: generating all basic examples that just produce images."
 for i in $image_scripts

examples/geometrically_infinite.py

@@ -23,11 +23,10 @@ def __init__(self, m):
         X = mp.matrix([[-self.λ*self.x, -(1+self.x**2)],[1,self.λ**-1*self.x]])
         super().__init__([T,X])
 
-depth = 80
-numpoints = 20000
+num_points = 10**6
 G = GeomInfGroup(3)
 seed = G.fixed_points((0,1))[0]
-df = G.coloured_limit_set_mc(depth,numpoints, seed=seed)
+df = G.coloured_limit_set_fast(num_points, seed=seed)
 scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
             .opts(marker = "dot", size = 0.1,  color = 'colour', width=2000, height=2000, data_aspect=1, cmap='Set1')\
               .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))

examples/grandma.py

@@ -26,12 +26,12 @@ def __init__(self, a,b):
 #       `colour` attribute returned from GroupCache.coloured_limit_set,
 #       i.e. it refers to the first letter in the word indexing that limit point); and
 #   (b) three points, which give the fixed points of the generators X and Y (the fourth fixed point is at infinity).
-def limit_set_points(are=1,aim=0,bre=0, bim=1, depth=15, logpoints=3):
+def limit_set_points(are=1,aim=0,bre=0, bim=1, logpoints=3):
     G = GrandmaGroup(are+1j*aim,bre+1j*bim)
     fixed_points_X = [ [float(p.real), float(p.imag)] for p in G.fixed_points((0,))]
     fixed_points_Y = [ [float(p.real), float(p.imag)] for p in G.fixed_points((1,))]
     seed = G.fixed_points((0,1))[0]
-    df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+    df = G.coloured_limit_set_fast(10**logpoints, seed=seed)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 0.1,  color = 'colour', width=800, height=800, data_aspect=1, cmap='Category10')\
                   .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))
@@ -43,6 +43,5 @@ def limit_set_points(are=1,aim=0,bre=0, bim=1, depth=15, logpoints=3):
                                               hv.Dimension('aim', label='Im(t_a)', range=(-4.0,4.0), step=.01, default=0),
                                               hv.Dimension('bre', label='Re(t_b)', range=(-4.0,4.0), step=.01, default=2),
                                               hv.Dimension('bim', label='Im(t_b)', range=(-4.0,4.0), step=.01, default=0),
-                                              hv.Dimension('depth', label='maximal length of word', range=(5,50), step=1, default=10),
-                                              hv.Dimension('logpoints', label='log10(number of words)', range=(2,6), default=3)])
+                                              hv.Dimension('logpoints', label='log10(number of points)', range=(2,8), default=4)])
 pn.panel(plot).servable(title="Grandma's Recipe groups")

examples/indras_necklace.py

@@ -29,12 +29,12 @@ def __init__(self, y, v):
 #       `colour` attribute returned from GroupCache.coloured_limit_set,
 #       i.e. it refers to the first letter in the word indexing that limit point); and
 #   (b) three points, which give the fixed points of the generators X and Y (the fourth fixed point is at infinity).
-def limit_set_points(y,v, depth=15, logpoints=3):
+def limit_set_points(y,v, logpoints=3):
     G = NecklaceGroup(y,v)
     fixed_points_X = [ [float(p.real), float(p.imag)] for p in G.fixed_points((0,))]
     fixed_points_Y = [ [float(p.real), float(p.imag)] for p in G.fixed_points((1,))]
     seed = G.fixed_points((0,1))[0]
-    df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+    df = G.coloured_limit_set_fast(10**logpoints, seed=seed)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 1,  color = 'colour', width=800, height=800, data_aspect=1, cmap='Category10')\
                   .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))
@@ -44,6 +44,5 @@ def limit_set_points(y,v, depth=15, logpoints=3):
 # Now we make a DynamicMap so that the user can modify the parameters.
 plot = hv.DynamicMap(limit_set_points, kdims=[hv.Dimension('y', label='y', range=(0.01,4.0), step=.01, default=1),
                                               hv.Dimension('v', label='v', range=(0.01,4.0), step=.01, default=1),
-                                              hv.Dimension('depth', label='maximal length of word', range=(5,50), step=1, default=10),
-                                              hv.Dimension('logpoints', label='log10(number of words)', range=(2,6), default=3)])
+                                              hv.Dimension('logpoints', label='log10(number of points)', range=(2,8), default=4)])
 pn.panel(plot).servable(title="Indra's Necklace groups")

examples/jorgensen_marden.py

@@ -28,11 +28,10 @@ def __init__(self):
         Y = (1/b)*mp.matrix([[2+1j,1j], [-1,-1j]])
         super().__init__([X,Y])
 
-depth = 50
-logpoints = 4
+num_points = 10**5
 def write_limit_set(G,filename):
     seed = G.fixed_points((0,1))[0]
-    df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+    df = G.coloured_limit_set_fast(num_points, seed=seed)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 0.1,  color = 'colour', width=1600, height=1600, data_aspect=1, cmap='Set1')\
                   .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))

examples/maskit.py

@@ -22,12 +22,12 @@ def __init__(self, μ):
 #       `colour` attribute returned from GroupCache.coloured_limit_set,
 #       i.e. it refers to the first letter in the word indexing that limit point); and
 #   (b) three points, which give the fixed points of the generators X and Y (the fourth fixed point is at infinity).
-def limit_set_points(μre=1,μim=0, depth=15, logpoints=3):
+def limit_set_points(μre=1,μim=0, logpoints=3):
     G = MaskitGroup(μre+1j*μim)
     fixed_points_X = [ [float(p.real), float(p.imag)] for p in G.fixed_points((0,))]
     fixed_points_Y = [ [float(p.real), float(p.imag)] for p in G.fixed_points((1,))]
     seed = G.fixed_points((0,1))[0]
-    df = G.coloured_limit_set_mc(depth,10**logpoints, seed=seed)
+    df = G.coloured_limit_set_fast(10**logpoints, seed=seed)
     scatter = hv.Scatter(df, kdims = ['x'], vdims = ['y','colour'])\
                 .opts(marker = "dot", size = 0.1,  color = 'colour', width=800, height=800, data_aspect=1, cmap='Category10')\
                   .redim(x=hv.Dimension('x', range=(-4,4)),y=hv.Dimension('y', range=(-4, 4)))
@@ -37,6 +37,5 @@ def limit_set_points(μre=1,μim=0, depth=15, logpoints=3):
 # Now we make a DynamicMap so that the user can modify the parameters.
 plot = hv.DynamicMap(limit_set_points, kdims=[hv.Dimension('μre', label='Re(μ)', range=(-6.0,6.0), step=.01, default=0),
                                               hv.Dimension('μim', label='Im(μ)', range=(-6.0,6.0), step=.01, default=2),
-                                              hv.Dimension('depth', label='maximal length of word', range=(5,50), step=1, default=10),
-                                              hv.Dimension('logpoints', label='log10(number of words)', range=(2,6), default=3)])
+                                              hv.Dimension('logpoints', label='log10(number of points)', range=(2,8), default=4)])
 pn.panel(plot).servable(title="Maskit groups")

examples/peripheral_subgroups.py

@@ -12,13 +12,13 @@
 def rotate(g):
     return  g[-1:] + g[:-1]
 
-def limit_set_points(r=1,s=3,mure=2, muim=2, depth=15, logpoints=3):
+def limit_set_points(r=1,s=3,mure=2, muim=2, logpoints=3):
     xbounds=(-1.5,1.5)
     ybounds=(-1.5,1.5)
 
     numpts = (10**logpoints)
     G = riley.ClassicalRileyGroup(mp.inf,mp.inf, mure+muim*1j)
-    big_limit_set = G.coloured_limit_set_mc(depth, numpts//3, (G.fixed_points((1,1)))[0])
+    big_limit_set = G.coloured_limit_set_fast(numpts//3, (G.fixed_points((1,1)))[0])
     big_scatter = hv.Scatter(big_limit_set, 'x','y')\
         .opts(marker = "dot", size = 1,  color = "gray", alpha=0.3, width=800, height=800, data_aspect=1)\
             .redim(x=hv.Dimension('x', range=xbounds),y=hv.Dimension('y', range=ybounds))
@@ -34,12 +34,12 @@ def limit_set_points(r=1,s=3,mure=2, muim=2, depth=15, logpoints=3):
 
     # Compute the two peripheral subgroups using GroupCache.subgroup().
     P = G.subgroup([G[left], G[left_middle]])
-    small_limit_set1 = P.coloured_limit_set_mc(depth, numpts//2, (P.fixed_points((1,1)))[0])
+    small_limit_set1 = P.coloured_limit_set_fast(numpts//2, (P.fixed_points((1,1)))[0])
     small_scatter1 = hv.Scatter(small_limit_set1, 'x','y')\
         .opts(marker = "dot", size = 1,  color = "red", width=800, height=800, data_aspect=1)\
             .redim(x=hv.Dimension('x', range=xbounds),y=hv.Dimension('y', range=ybounds))
     P = G.subgroup([G[right_middle], G[right]])
-    small_limit_set2 = P.coloured_limit_set_mc(depth, numpts//2, (P.fixed_points((1,1)))[0])
+    small_limit_set2 = P.coloured_limit_set_fast(numpts//2, (P.fixed_points((1,1)))[0])
     small_scatter2 = hv.Scatter(small_limit_set2, 'x','y')\
         .opts(marker = "dot", size = 1,  color = "blue", width=800, height=800, data_aspect=1)\
             .redim(x=hv.Dimension('x', range=xbounds),y=hv.Dimension('y', range=ybounds))
@@ -49,7 +49,6 @@ def limit_set_points(r=1,s=3,mure=2, muim=2, depth=15, logpoints=3):
                                               hv.Dimension('s', range=(0,5), default=3),
                                               hv.Dimension('mure', label='Re(μ)', range=(-8.0,8.0), step=.01, default=1.55),
                                               hv.Dimension('muim', label='Im(μ)', range=(-8.0,8.0), step=.01, default=1.41),
-                                              hv.Dimension('depth', label='maximal length of word', range=(5,50), step=1, default=15),
-                                              hv.Dimension('logpoints', label='log10(number of words)', range=(2,6), default=3)])
+                                              hv.Dimension('logpoints', label='log10(number of points)', range=(2,8), default=4)])
 
 pn.panel(plot).servable(title='Peripheral subgroups')