Find the identity sandpile

Question

This question is about abelian sandpiles. Read this previous challenge and watch this numberphile video to learn more.

---

An abelian sandpile of size n by n is a grid containing the number 0, 1, 2 and 3 (representing the number of grains of sand). Adding two sandpiles works by first adding element by element, and then toppling any element that goes above 3. The order in which you topple doesn't matter, the end result is the same. When a cell topples its number decreases by 4, and each of its direct neighbors increases by 1. This can cause a chain reaction. If a cell is on the edge of the grid, any grains that fall off the grid while toppling disappear.

For example, I'm adding two 3 by 3 sandpiles (giving a rather extreme chain reaction):

3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2
3 3 3 + 2 1 2 = 5 4 5 -> 6 0 6 -> 2 4 2 -> 3 0 3 -> 5 0 5 -> 1 4 1 -> 2 0 2 -> 4 0 4 -> 0 4 0 -> 1 0 1
3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2

In this challenge we're interested in a subset of all possible n by n sandpiles. This subset contains any sandpile that you can get by adding an arbitrary sandpile to the all-3s n by n sandpile. For example, just above we saw that 212 | 101 | 212 is in the subset, because we got it by adding something to the all-3 sandpile.

Now this subset has an interesting element: the identity element. If you take this element and add it to any other element in the subset, the sum is unchanged. In other words, this sandpile acts like a zero of this subset. It just so happens that 212 | 101 | 212 is the zero element for the subset of 3 by 3. For example:

2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2
2 2 2 + 1 0 1 = 3 2 3 -> 5 2 5 -> 1 6 1 -> 2 2 2
2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2

Now this is your challenge: given n, find the identity element of the subset of the n by n grid. Output it by assigning a unique color with sufficient contrast of your choice to each of 0, 1, 2, 3 and outputting an n by n image. Your code must be able to produce the 50 by 50 case in under a minute on a reasonable modern PC.

---

For example, the 500 by 500 identity element:

Here is blue = 3, green = 2, red = 1, white = 0. But you don't have to use this color scheme in your answer.

Answer 1

Octave, 120 113 bytes

function a=W(a);while nnz(b=a>3);a+=conv2(b,[t=[0 1 0];!t;t],'same')-b*4;end;end;@(n)imagesc(W((x=ones(n)*6)-W(x)))

Thanks to JungHwan Min for providing a link to the reference paper in his Mathematica answer.
Thanks to Stewie Griffin saved me 7 bytes [any(any(x)) -> nnz(x)]

Here two function is used:

1.f: for stabilization of a matrix
2. An anonymous function that takes n as input and shows the identity matrix.

Try It On rextester! for generation of a 50 * 50 matrix

Elapsed time for computation of the matrix: 0.0844409 seconds.

Explanation:

Consider a function f that stabilizes a matrix the task of finding the identity is simply

f(ones(n)*6 - f(ones(n)*6).

that ones(n)*6 means a n*n matrix of 6.

so for n=3:

M = [6 6 6
     6 6 6
     6 6 6];

The result will be f(M-f(M))

For stabilization function 2D convolution used to speed the task; In each iteration we make a binary matrix b with the same size of the input matrix and set it to 1 if corresponding element of the input matrix is >3. Then we apply a 2D convolution of the binary matrix with the following mask

0 1 0
1 0 1
0 1 0

representing four direct neighbors .
Result of convolution is added to the matrix and 4 times the binary matrix subtracted from it.

The loop continued until all element of the matrix are <= 3

Ungolfed version:

function a=stabilize(a)
    mask = [t=[0 1 0];!t;t];
    while any(any(b=a>3))
        a+=conv2(b,mask,'same')-b*4;
    end
end
n= 50;
M = ones(n)*6;
result = stabilize(M-stabilize(M));
imagesc(result);

Answer 2

Mathematica, 177 157 135 133 bytes

Colorize[f=BlockMap[⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&,#~ArrayPad~1,{3,3},1]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&

Takes a number n. The output is the identity sandpile. 0 is black, 1 is light gray, 2 is magenta, and 3 is blue-gray.

Sadly, Mathematica does not have a builtin for this...

Uses the algorithm stated in Scott Corry and David Perkinson's paper.

Takes 91.7 seconds on my 5-year-old laptop to calculate the 50x50 identity sandpile. I am confident that a reasonable modern desktop computer is more than 50% faster. (I also have a way faster code at the end).

Explanation

f= ...

Define function f (the input is a sandpile matrix): a function that ...

BlockMap[ ... ]~FixedPoint~#&

... repeats the BlockMap operation until the output does not change. BlockMap operation: ...

---

#~ArrayPad~1

... pad the input array with one layer of 0 ...

{3,3},1

... partition it into 3x3 matrices, with offset 1 ...

⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&

... and for each partition, add the number of sand grains toppled onto the center cell and the center cell value mod 4.

i.e. the output of f is the stabilized version of the input.

---

k=Table[6,#,#]

Define k as an n by n array of 6s.

f[k-f@k]]

Compute f(k - f(k)).

Colorize[ ... ]

Apply colors to the result.

Faster version (142 bytes)

Colorize[r=RotateLeft;a=ArrayPad;f=a[b=#~a~1;b+r[g=⌊b/4⌋,s={0,1}]+g~r~-s+r[g,1-s]+r[g,s-1]-4g,-1]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&

Same code, but uses built-in list rotation instead of BlockMap. Computes n=50 in 4.0 seconds on my laptop.

Answer 3

Python 3 + Numpy + PIL, 385 370 364 bytes

import numpy as q,PIL.Image as w
n=int(input())
z=n,n
def r(p):
 while len(p[p>3]):
  for x,y in q.ndindex(z):
   if p[x,y]>3:
    p[x,y]-=4;p[x-1,y]+=x>0;p[x,y-1]+=y>0
    if~-n>x:p[x+1,y]+=1
    if~-n>y:p[x,y+1]+=1
s=q.full(z,6)
t=s.copy()
r(t)
i=s-t
r(i)
w.fromarray(q.uint8(q.array(q.vectorize(lambda x:[x//1*65]*3,otypes=[object])(i).tolist()))).save('i.png')

Takes input on STDIN. Outputs the image as greyscale to i.png. Black corresponds to 0, dark grey to 1, light grey to 2, and white to 0.

Uses the formula I = R(S - R(S)), where I is the identity element, S is the matrix filled with sixes, and R is the reduction function.

I could probably save some bytes by switching to Python 2 and doing from numpy import*, but (1) I don't have Numpy installed on Python 2 and (2) the program wasn't terminating with from numpy import*.

Ungolfed:

import numpy as np
from PIL import Image

# Compute the identity element

n = int(input('Size of the sandpile: '))

def reduce_pile(sandpile):
  while any(element >= 4 for element in np.nditer(sandpile)):
    for x, y in np.ndindex((n, n)):
      if sandpile[x, y] >= 4:
        sandpile[x, y] -= 4
        if x > 0: sandpile[x - 1, y] += 1
        if y > 0: sandpile[x, y - 1] += 1
        if x < n - 1: sandpile[x + 1, y] += 1
        if y < n - 1: sandpile[x, y + 1] += 1

s = np.full((n, n), 6, dtype=np.int32)
s_prime = np.copy(s)

reduce_pile(s_prime)

identity = s - s_prime
reduce_pile(identity)

# Output it to identity.png as an image

colours = [[255, 255, 255], [255, 0, 0], [0, 255, 0], [0, 0, 255]]
img_array = np.vectorize(lambda x: colours[x], otypes=[object])(identity)
img_array = np.array(img_array.tolist(), dtype=np.uint8)

img = Image.fromarray(img_array)
img.save('identity.png')