10.6. Fractals¶
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 22
.
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 level
and False
otherwise. 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)
The parameter, 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 n
and m
.
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: i
, i**2
, and total
.
Here’s how we use count_cells
:
res = count_cells(pile.array==level)
steps, steps2, cells = res
The first line creates a boolean array that contains True
where the array equals level
, calls count_cells
, and gets an array with three rows.
The second line unpacks the rows and assigns them to steps
, steps2
, and 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[0]
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.