Just a couple of days ago I posted about the Numberphile video about Frog Jumping (see here) where one of the forum members suggested I should check out the video about sandpiles.

If you haven't seen it already, you should. I must admit it is a fairly long video that starts off boring, but it is pure gold at the end. Anyhow, how can we do such calculations in the Wolfram Language I van Veen was wondering. So here we go.

We represent our sandpile as a simple matrix (following the example around 8:37 in the video):

a = sp[{{3, 3, 3}, {3, 3, 3}, {3, 3, 3}}]
b = sp[{{0, 0, 0}, {0, 1, 0}, {0, 0, 0}}]
c = a + b

If we execute this now, nothing really happens, that's because we have not yet told the computer how to add up sp objects. So let's add a definition:

ClearAll[sp]
sp[s_List] + sp[t_List] ^:= Module[{dim, r, tmp, neighbours},
  dim = Dimensions[s];
  r = s + t;
  While[Max[r] > 3,
   r = ArrayPad[r, 1, 0];
   tmp = Quotient[r, 4];
   r -= 4 tmp;
   r += RotateLeft[tmp, {0, 1}] + RotateLeft[tmp, {1, 0}] + 
     RotateLeft[tmp, {0, -1}] + RotateLeft[tmp, {-1, 0}];
   r = ArrayPad[r, -1];
   ];
  sp[r]
  ]

What is done here, is the following:

1. Add up the raw matrices
2. While there is a number still above 3 loop:
3. Pad the entire matrix with 0s around it
4. Calculate Floor[values/4], this is done using quotient, and gives a matrix of the number that can be toppled, and how many times it can be toppled.
5. Subtract 4 times this matrix.
6. Add shifted version of this matrix to redistribute the sand to the neighbours (using RotateLeft).
7. Remove the padding again using Arraypad[...,-1]

Note that I have to use :^= rather than := because Plus is protected by default. Moreover, I don't want to touch the definition of Plus, so by using ^:= the definition is added to sp rather than Plus. I also tried out several other different methods, I found this was one of the fastest ways. Other options included using Position, finding the neighbours and subtracting manually, convolutions, or using BlockMap. Perhaps one could even coerce CellularAutomaton to do it for you.

Now that we have addition figures out, we can re-execute our above code:

a = sp[{{3, 3, 3}, {3, 3, 3}, {3, 3, 3}}];
b = sp[{{0, 0, 0}, {0, 1, 0}, {0, 0, 0}}];
c = a + b

sp[{{1, 3, 1}, {3, 0, 3}, {1, 3, 1}}]

So it is working! Now make some convenient functions to make it look nicer:

colorrules = {0 -> Black, 1 -> Yellow, 2 -> Blue, 3 -> Red, _Integer -> Pink};
Format[sp[t_]] := ArrayPlot[t /. colorrules, ImageSize -> 200, Frame -> False]

Using Format we can tell the Wolfram Language how to show an sp object. Re-executing now gives us an image, rather than a dull matrix:

However we can not get the data from such an image any more, so it is a good idea to define a Normal and some other functions as well, here the total code:

ClearAll[sp]
sp[x_SparseArray]:=sp[Normal[x]]
Normal[sp[s_List]]^:=s
sp[s_List]+sp[n_Integer]^:=sp[s]+sp[ConstantArray[n,Dimensions[s]]]
sp[s_List]+sp[t_List]^:=Module[{dim,r,tmp,neighbours},
    dim=Dimensions[s];
    r=s+t;
    While[Max[r]>3,
        r=ArrayPad[r,1,0];
        tmp=Quotient[r,4];
        r-=4tmp;
        r+=RotateLeft[tmp,{0,1}]+RotateLeft[tmp,{1,0}]+RotateLeft[tmp,{0,-1}]+RotateLeft[tmp,{-1,0}];
        r=ArrayPad[r,-1];
    ];
    sp[r]
]
colorrules={0->Black,1->Yellow,2->Blue,3->Red,_Integer->Pink};
Format[sp[t_]]:=ArrayPlot[t/.colorrules,ImageSize->200,Frame->False]

Now we can type:

Normal[c]

to get the values again:

{{1, 3, 1}, {3, 0, 3}, {1, 3, 1}}

Or add a constant sandpile:

sp[{{3, 3, 3}, {3, 3, 3}, {3, 3, 3}}] + sp[1]

So let's create the iconic image the Numberphile video showed:

n = 100;
u = sp[ConstantArray[6, {n, n}]];
start = AbsoluteTime[];
u = u + sp[0];  (* the addition of 0 will trigger the toppling *)
out = u = sp[ConstantArray[6, {n, n}] - Normal[u]] + sp[0]
{n, AbsoluteTime[] - start}

Letting it run takes a few seconds:

Or run it for large n = 512:

The authors also shows what happens when you put a large pile of sand on a single cell and see how it spreads, we can do that too:

u = sp[CenterArray[50000, {201, 201}]];
AbsoluteTiming[out = u + sp[0];]
out

Note that for larger numbers of sand it will take longer and longer to converge; this is because sand will take a long time to reach the edge of the matrix. The grid increases as n^2, but the border only as n. So it becomes harder and harder to lose sand. I'm not completely sure how to speed up the code without going through all the calculations...

Hope you enjoyed this, now that you have the basic code, one can try non-rectangular patterns, different neighbourhoods, et cetera. Also check out the other Numberphile inspired posts:

• [Numberphile] - Frog Jumping - Solving the puzzle
• [Numberphile] - The Illumination Problem