

 This page uses content from Wikipedia. The original article was at Abelian sandpile model. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



Implement the Abelian sandpile model also known as Bak–Tang–Wiesenfeld model. Its history, mathematical definition and properties can be found under its wikipedia article.

The task requires the creation of a 2D grid of arbitrary size on which "piles of sand" can be placed. Any "pile" that has 4 or more sand particles on it collapses, resulting in four particles being subtracted from the pile and distributed among its neighbors.

It is recommended to display the output in some kind of image format, as terminal emulators are usually too small to display images larger than a few dozen characters tall. As an example of how to accomplish this, see the Bitmap/Write a PPM file task.

Examples up to 2^30, wow!
javascript running on web
Examples:

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 4 0 0 -> 0 1 0 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 6 0 0 -> 0 1 2 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0  0 0  0  0    0 0 1 0 0
0  0 0  0  0    0 2 1 2 0
0  0 16 0  0 -> 1 1 0 1 1
0  0 0  0  0    0 2 1 2 0
0  0 0  0  0    0 0 1 0 0

module AbelSand
 
# supports output functionality for the results of the sandpile simulations
# outputs the final grid in CSV format, as well as an image file
 
using CSV, DataFrames, Images
 
function TrimZeros(A)
    # given an array A trims any zero rows/columns from its borders
    # returns a 4 tuple of integers, i1, i2, j1, j2, where the trimmed array corresponds to A[i1:i2, j1:j2]
    # A can be either numeric or a boolean array
 
    i1, j1 = 1, 1
    i2, j2 = size(A)
 
    zz = typeof(A[1, 1])(0)    # comparison of a value takes into account the type as well
 
    # i1 is the first row which has non zero element
    for i = 1:size(A, 1)
        q = false
        for k = 1:size(A, 2)
            if A[i, k] != zz
                q = true
                i1 = i
                break
            end
        end
 
        if q == true
            break
        end
    end
 
    # i2 is the first from below row with non zero element
    for i in size(A, 1):-1:1
        q = false
        for k = 1:size(A, 2)
            if A[i, k] != zz
                q = true
                i2 = i
                break
            end
        end
 
        if q == true
            break
        end
    end
 
    # j1 is the first column with non zero element
 
    for j = 1:size(A, 2)
        q = false
        for k = 1:size(A, 1)
            if A[k, j] != zz
                j1 = j
                q = true
                break
            end
        end
 
        if q == true
            break
        end
    end
 
    # j2 is the last column with non zero element
 
    for j in size(A, 2):-1:1
        q=false
        for k=1:size(A,1)
            if A[k, j] != zz
                j2 = j
                q=true
                break
            end
        end
 
        if q==true
            break
        end
    end
 
    return i1, i2, j1, j2
end
 
function addLayerofZeros(A, extraLayer)
    # adds layer of zeros from all corners to the given array A
 
    if extraLayer <= 0
        return A
    end
 
    N, M = size(A)
 
 
    Z = zeros( typeof(A[1,1]), N + 2*extraLayer, M + 2*extraLayer)
    Z[(extraLayer+1):(N + extraLayer ), (extraLayer+1):(M+extraLayer)] = A
 
    return Z
 
end
 
function printIntoFile(A, extraLayer, strFileName, TrimSmallValues = false)
    # exports a 2d matrix A into a csv file
    # @extraLayer is an integers adding layer of 0-s sorrounding the output matrix
 
    # trimming off very small values; tiny values affect the performance of CSV export
    if TrimSmallValues == true
        A = map(x -> if (abs(x - floor(x)) < 0.01) floor(x) else x end, A) 
    end
 
    i1, i2, j1, j2  = TrimZeros( A )
    A = A[i1:i2, j1:j2]
 
    A = addLayerofZeros(A, extraLayer)
 
    CSV.write(string(strFileName,".csv"), DataFrame(A), writeheader = false)
 
    return A
 
end
 
function Array_magnifier(A, cell_mag, border_mag)
    # A is the main array; @cell_mag is the magnifying size of the cell,
    # @border_mag is the magnifying size of the border between lattice cells
 
    # creates a new array where each cell of the original array A appears magnified by size = cell_mag
 
 
    total_factor = cell_mag + border_mag
 
    A1 = zeros(typeof(A[1, 1]), total_factor*size(A, 1), total_factor*size(A, 2))
 
    for i = 1:size(A,1), j = 1:size(A,2), u = ((i-1)*total_factor+1):(i*total_factor),
                                          v = ((j-1)*total_factor+1):(j*total_factor)
        if(( u - (i - 1) * total_factor <= cell_mag) && (v - (j - 1) * total_factor <= cell_mag))
            A1[u, v] = A[i, j] 
        end
    end
 
    return A1
 
end
 
function saveAsGrayImage(A, fileName, cell_mag, border_mag, TrimSmallValues = false)
    # given a 2d matrix A, we save it as a gray image after magnifying by the given factors
    A1 = Array_magnifier(A, cell_mag, border_mag)
    A1 = A1/maximum(maximum(A1))
 
    # trimming very small values from A1 to improve performance
    if TrimSmallValues == true
        A1 = map(x -> if ( x < 0.01) 0.0 else round(x, digits = 2) end, A1) 
    end
 
    save(string(fileName, ".png") , colorview(Gray, A1)) 
end
 
function saveAsRGBImage(A, fileName, color_codes, cell_mag, border_mag)
    # color_codes is a dictionary, where key is a value in A and value is an RGB triplet
    # given a 2d array A, and color codes (mapping from values in A to RGB triples), save A
    # into fileName as png image after applying the magnifying factors
 
    A1 = Array_magnifier(A, cell_mag, border_mag)
    color_mat = zeros(UInt8, (3, size(A1, 1), size(A1, 2)))
 
    for i = 1:size(A1,1)
        for j = 1:size(A1,2)
            color_mat[:, i, j]  = get(color_codes, A1[i, j] , [0, 0, 0])
        end
    end
 
    save(string(fileName, ".png") , colorview(RGB, color_mat/255)) 
end
 
const N_size = 700       # the radius of the lattice Z^2, the actual size becomes (2*N+1)x(2*N+1)
const dx = [1, 0, -1, 0] # for a given (x,y) in Z^2, (x + dx, y + dy) for all (dx,dy) covers the neighborhood of (x,y)
const dy = [0, 1, 0, -1]
 
struct L_coord
    # represents a lattice coordinate
    x::Int
    y::Int
end
 
function FindCoordinate(Z::Array{L_coord,1}, a::Int, b::Int)
    # in the given array Z of coordinates finds the (first) index of the tuple (a,b)
    # if no match, returns -1
 
    for i=1:length(Z)
        if (Z[i].x == a) && (Z[i].y == b)
            return i
        end
    end
 
    return -1
end
 
function move(N)
    # the main function moving the pile sand grains of size N at the origin of Z^2 until the sandpile becomes stable
 
    Z_lat = zeros(UInt8, 2 * N_size + 1, 2 * N_size + 1)     # models the integer lattice Z^2, we will have at most 4 sands on each vertex
    V_sites = falses(2 * N_size + 1, 2 * N_size + 1)         # all sites which are visited by the sandpile process, are being marked here
    Odometer = zeros(UInt64, 2 * N_size + 1, 2 * N_size + 1) # stores the values of the odometer function
 
 
    walking = L_coord[]    # the coordinates of sites which need to move
 
    V_sites[N_size + 1, N_size + 1] = true
 
    # i1, ... j2  -> show the boundaries of the box which is visited by the sandpile process
    i1, i2, j1, j2 = N_size + 1, N_size + 1, N_size + 1, N_size + 1 
    n = N
 
    t1 = time_ns()
 
    while n > 0
        n -= 1
 
        Z_lat[N_size + 1, N_size + 1] += 1
        if (Z_lat[N_size + 1, N_size + 1] >= 4)
            push!(walking, L_coord(N_size + 1, N_size + 1))
        end
 
        while(length(walking) > 0)
            w = pop!(walking)
            x = w.x
            y = w.y
 
            Z_lat[x, y] -= 4
            Odometer[x, y] += 4
 
            for k = 1:4
                Z_lat[x + dx[k], y + dy[k]] += 1
                V_sites[x + dx[k], y + dy[k]] = true
                if Z_lat[x + dx[k], y + dy[k]] >= 4
                    if FindCoordinate(walking, x + dx[k] , y + dy[k]) == -1
                        push!(walking, L_coord( x + dx[k], y + dy[k]))
                    end
                end
            end
 
            i1 = min(i1, x - 1)
            i2 = max(i2, x + 1)
            j1 = min(j1, y - 1)
            j2 = max(j2, y + 1)
        end
 
 
    end #end of the main while
    t2 = time_ns()
 
    println("The final boundaries are:: ", (i2 - i1 + 1),"x",(j2 - j1 + 1), "\n")
    print("time elapsed: " , (t2 - t1) / 1.0e9, "\n")
 
    Z_lat = printIntoFile(Z_lat, 0, string("Abel_Z_", N))
    Odometer = printIntoFile(Odometer, 1, string("Abel_OD_", N))
 
    saveAsGrayImage(Z_lat, string("Abel_Z_", N), 20, 0)
    color_code = Dict(1=>[255, 128, 255], 2=>[255, 0, 0],3 => [0, 128, 255])
    saveAsRGBImage(Z_lat, string("Abel_Z_color_", N), color_code, 20, 0)
 
    # for the total elapsed time, it's better to use the @time macros on the main call
 
    return Z_lat, Odometer # these are trimmed in output module
 
end # end of function move
 
 
end # module
 
 
using .AbelSand
 
Z_lat, Odometer = AbelSand.move(100000)
 