Chaos Ferns

James Gleick’s book Chaos: Making a New Science is full of fascinating examples of chaos theory. One such example is Michael Barnsley’s “chaos game”. In this “game”, an image of a fractal fern is created by iteratively coloring points in an initially blank page. The points are colored with only four rules, or transformations, selected at random based on predefined probabilities. The starting point \(\mathbf{x}_0\) is \((0, 0)\), and, from there, the next points are calculated as:

\[ \mathbf{x}_{n+1} = f(\mathbf{x}_{n}), \tag{1}\]

where the transformation \(f\) is given by [1]: \[ \begin{align} f_1(x,y) &= \begin{bmatrix} 0.00 & 0.00 \\ 0.00 & 0.16 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \ \text{ , probability 0.01} \\ f_2(x,y) &= \begin{bmatrix} 0.85 & 0.04 \\ -0.04 & 0.85 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} 0.00 \\ 1.60 \end{bmatrix} \ \text{ , probability 0.85} \\ f_3(x,y) &= \begin{bmatrix} 0.20 & -0.26 \\ 0.23 & 0.22 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} 0.00 \\ 1.60 \end{bmatrix} \ \text{ , probability 0.07} \\ f_4(x,y) &= \begin{bmatrix} -0.15 & 0.28 \\ 0.26 & 0.24 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} 0.00 \\ 0.44 \end{bmatrix} \ \text{ , probability 0.07} \end{align} \tag{2}\]

Let’s make some imports:

Code

import random

import colorcet
import datashader as ds
import numpy as np
import pandas as pd
from numba import jit

Now, let’s create a function that encodes the rules in Equation 2:

Code

@jit(nopython=True)
def get_next_x_y(x, y):
    random_number = random.random()

    if random_number < 0.01:
        xn = 0
        yn = 0.16 * y
    elif random_number < 0.86:
        xn = 0.85 * x + 0.04 * y
        yn = -0.04 * x + 0.85 * y + 1.6
    elif random_number < 0.93:
        xn = 0.2 * x - 0.26 * y
        yn = 0.23 * x + 0.22 * y + 1.6
    else:
        xn = -0.15 * x + 0.28 * y
        yn = 0.26 * x + 0.24 * y + 0.44

    return xn, yn

We are using Numba’s @jit decorator to optimize the function with a just-in-time compiler. This makes the code run faster because it translates the Python code into machine code, taking advantage of the specialized features of the target architecture.

Now, all we have to do is call the get_next_x_y function millions of times and plot the \(x, y\) pairs we get. Let’s first generate ten million points:

Code

n_points = 10_000_000
x_arr = np.zeros(n_points)
y_arr = np.zeros(n_points)

x, y = 0, 0

for i in range(n_points):
    x_arr[i] = x
    y_arr[i] = y
    x, y = get_next_x_y(x, y)

Now, we’re going to plot the results using Datashader instead of the more commonly used Matplotlib library due to the sheer number of points. Datashader efficiently transforms the data behind the scenes into a visually interpretable representation that avoids many of the pitfalls that can occur when plotting large datasets.

Code

points = pd.DataFrame({"x": x_arr, "y": y_arr})

cvs = ds.Canvas(plot_width=800, plot_height=800)
agg = cvs.points(points, "x", "y")
ds.tf.set_background(ds.tf.shade(agg, cmap=colorcet.kgy), "black")

Figure 1: Ten million-point Barnsley fern.

Nice fern! It is fascinating to see how such a complex shape emerges from only four transformations. Equation 2 shows that this particular fern shape emerges from a set of transformations defined by the coefficients in some matrices and vectors. What happens if we change these coefficients? Would we still get ferns, some other shape, or purely random blobs of points?

References

[1]

Barnsley, M., “Fractals Everywhere,” Morgan Kaufmann, Oxford, England, 2000.