Sandpiles

Implementing the sandpile in from Numberphile Video

We start by definding a type definition

In [1]:
type SandPile
    M::Array{Int,2}
    n::Tuple{Int,Int}
    SandPile(a::Array{Int,2})=new(a,size(a))
end

We create the avalanche function to propagate the overflowing in the Sandpile

In [2]:
function avalanche(A::Array{Int,2},limit=3)
    while any(A.>limit) 
        R = CartesianRange(size(A))
        B=zeros(A)
        for i in R
            if A[i]>limit
                B[i]-=(limit+1)
                modPts=map(c->i+CartesianIndex(c),vcat([(z,0) for z=(-1,1)],[(0,z) for z=(-1,1)]))
                for j in modPts
                    (j in R) && (B[j]+=1)
                end
            end
        end
        A=A+B
    end
    A
end
Out[2]:
avalanche (generic function with 2 methods)

We extend the addition operation to handle sandpiles

In [3]:
import Base.+
function +(a::SandPile,b::SandPile) 
    a.n==b.n || error("Sandpiles must have the same dimenstion")
    c=a.M+b.M
    SandPile(avalanche(c))
end
Out[3]:
+ (generic function with 164 methods)

We write down sandpile zero for the 3x3 sandpile

In [4]:
sp_zero=SandPile([2 1 2
                  1 0 1
                  2 1 2])
Out[4]:
SandPile([2 1 2; 1 0 1; 2 1 2],(3,3))

We test that zero+zero=zero

In [5]:
sp_zero+sp_zero
Out[5]:
SandPile([2 1 2; 1 0 1; 2 1 2],(3,3))

We do another test for a number in $\mathcal{S}$

In [6]:
SandPile([2 2 2
          2 2 2
          2 2 2])+sp_zero
Out[6]:
SandPile([2 2 2; 2 2 2; 2 2 2],(3,3))

It works!

Now we test adding and number to its Negative

In [7]:
s3=SandPile(ones(Int,3,3)*3)
ns3=ones(Int,3,3)*3
ns3[2,2]=1
s̄3=SandPile(ns3)
s3+s̄3
Out[7]:
SandPile([2 1 2; 1 0 1; 2 1 2],(3,3))

It works again!

Now create some sanpile in $\mathcal{S}$ and try again.

In [8]:
@show sp_rnd=rand(0:3,3,3)
new_sp=SandPile(avalanche(ones(Int,3,3)*3+sp_rnd))
sp_rnd = rand(0:3,3,3) = [3 1 3; 3 0 2; 3 3 0]
Out[8]:
SandPile([3 3 3; 2 2 0; 0 2 1],(3,3))
In [9]:
new_sp+sp_zero
Out[9]:
SandPile([3 3 3; 2 2 0; 0 2 1],(3,3))

$a+0=a$... we are done!