Another property of critical systems is fractal geometry. The initial configuration in Figure 10.1 (left) resembles a fractal, but you can’t always tell by looking. A more reliable way to identify a fractal is to estimate its fractal dimension, as we saw in Section 9.6 and Section 9.7.
We’ll start by making a bigger sand pile, with
n=131 and initial level
pile3 = SandPile(n=131, level=22) pile3.run()
It takes 28,379 steps for this pile to reach equilibrium, with more than 200 million cells toppled.
To see the resulting pattern more clearly, we can select cells with levels 0, 1, 2, and 3, and plot them separately:
def draw_four(viewer, levels=range(4)): thinkplot.preplot(rows=2, cols=2) a = viewer.viewee.array for i, level in enumerate(levels): thinkplot.subplot(i+1) viewer.draw_array(a==level, vmax=1)
draw_four takes a
SandPileViewer object, which is defined in
Sand.py in the repository for this book. The parameter
levels is the list of levels we want to plot; the default is the range 0 through 3. If the sand pile has run until equilibrium, these are the only levels that should exist.
Inside the loop, it uses
a==level to make a boolean array that’s
True where the array is
draw_array treats these booleans as 1s and 0s.
Figure 10.4 shows the results for
pile3. Visually, these patterns resemble fractals, but looks can be deceiving. To be more confident, we can estimate the fractal dimension for each pattern using box-counting, as we saw in Section 9.6.
We’ll count the number of cells in a small box at the center of the pile, then see how the number of cells increases as the box gets bigger. Here’s an implementation:
def count_cells(a): n, m = a.shape end = min(n, m) res =  for i in range(1, end, 2): top = (n-i) // 2 left = (m-i) // 2 box = a[top:top+i, left:left+i] total = np.sum(box) res.append((i, i**2, total)) return np.transpose(res)
a, is a boolean array. The size of the box is initially 1. Each time through the loop, it increases by 2 until it reaches
end, which is the smaller of
Each time through the loop,
box is a set of cells with width and height
i, centered in the array.
total is the number of “on” cells in the box.
The result is a list of tuples, where each tuple contains
i, \(i^2\), and the number of cells in the box.
When we pass this result to
transpose, NumPy converts it to an array with three columns, and then transposes it; that is, it makes the columns into rows and the rows into columns. The result is an array with 3 rows:
Here’s how we use
res = count_cells(pile.array==level) steps, steps2, cells = res
The first line creates a boolean array that contains
True where the array equals
count_cells, and gets an array with three rows.
The second line unpacks the rows and assigns them to
cells, which we can plot like this:
thinkplot.plot(steps, steps2, linestyle='dashed') thinkplot.plot(steps, cells) thinkplot.plot(steps, steps, linestyle='dashed')
Figure 10.5 shows the results. On a log-log scale, the cell counts form nearly straight lines, which indicates that we are measuring fractal dimension over a valid range of box sizes.
To estimate the slopes of these lines, we can use the SciPy function
linregress, which fits a line to the data by linear regression.
from scipy.stats import linregress params = linregress(np.log(steps), np.log(cells)) slope = params
The estimated fractal dimensions are:
0 1.871 1 3.502 2 1.781 3 2.084
The fractal dimension for levels 0, 1, and 2 seems to be clearly non-integer, which indicates that the image is fractal.
The estimate for level 3 is indistinguishable from 2, but given the results for the other values, the apparent curvature of the line, and the appearance of the pattern, it seems likely that it is also fractal.
One of the exercises in the notebook for this chapter asks you to run this analysis again with different values of
n and the initial
level to see if the estimated dimensions are consistent.