I'm trying to do the Abelian sandpile model using a simple numpy array.
When a 'pile' is 4 >=, then it collapse among its neighbors.
I understand how the "gravity" thing works, but I can't think of a way of making it.
Here's the code to make my array :
import numpy as np
spile = np.zeros((5, 5), dtype=np.uint32)
spile[2, 2] = 16
Which gives me the following :
array([[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 16, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0]], dtype=uint32)
Now, I need the "gravity" code that does these steps of calculation :
array([[ 0, 0, 0, 0, 0],
[ 0, 0, 4, 0, 0],
[ 0, 4, 0, 4, 0],
[ 0, 0, 4, 0, 0],
[ 0, 0, 0, 0, 0]], dtype=uint32)
array([[ 0, 0, 1, 0, 0],
[ 0, 2, 1, 2, 0],
[ 1, 1, 0, 1, 1],
[ 0, 2, 1, 2, 0],
[ 0, 0, 1, 0, 0]], dtype=uint32)
The last array is the final result I'm trying to get.
I'm not trying to make you guys code for me, I just need some ideas as I've never ever did such a thing (but feel free to provide a code if you're that kind :p ).
Use np.divmod to identify where the cells tumble and how much tumbles. Then use array slicing to shift the amounts tumbled and add back into the sandpile.
import numpy as np
spile = np.zeros((5, 5), dtype=np.uint32)
spile[2, 2] = 16
def do_add( spile, tumbled ):
""" Updates spile in place """
spile[ :-1, :] += tumbled[ 1:, :] # Shift N and add
spile[ 1:, :] += tumbled[ :-1, :] # Shift S
spile[ :, :-1] += tumbled[ :, 1:] # Shift W
spile[ :, 1:] += tumbled[ :, :-1] # Shift E
def tumble( spile ):
while ( spile > 3 ).any():
tumbled, spile = np.divmod( spile, 4 )
do_add( spile, tumbled )
# print( spile, '\n' ) # Uncomment to print steps
return spile
print( tumble( spile ) )
# or tumble( spile ); print( spile )
# [[0 0 1 0 0]
# [0 2 1 2 0]
# [1 1 0 1 1]
# [0 2 1 2 0]
# [0 0 1 0 0]]
Uncommented print statement prints these results
[[0 0 0 0 0]
[0 0 4 0 0]
[0 4 0 4 0]
[0 0 4 0 0]
[0 0 0 0 0]]
[[0 0 1 0 0]
[0 2 0 2 0]
[1 0 4 0 1]
[0 2 0 2 0]
[0 0 1 0 0]]
[[0 0 1 0 0]
[0 2 1 2 0]
[1 1 0 1 1]
[0 2 1 2 0]
[0 0 1 0 0]]