Recursive Abelian Sandpile algorithm won't return final value

I am trying to implement a simple abelian sandpile algorithm using a recursive call. The function takes a matrix of positive numbers as input. I structured my code as follows:

function[thematrix] = abelian(sandmatrix)
for c = 1:size(sandmatrix,2)
  for n = 1:size(sandmatrix,1)
      if sandmatrix(n,c) >= 4
          sandmatrix(n,c) = sandmatrix(n,c)-4
          sandmatrix(n+1,c) = sandmatrix(n+1,c) + 1
          sandmatrix(n,c+1) = sandmatrix(n,c+1) + 1
          sandmatrix(n-1,c) = sandmatrix(n-1,c) + 1
          sandmatrix(n,c-1) = sandmatrix(n,c-1) + 1
          if any(max(sandmatrix) >= 4)
              abelian(sandmatrix)
          else
              break
          end
      end
  end
end

The code seems to run well, and produces the sandpile behavior that I would expect. However, it isn't returning the final sandpile. Instead, it stores the one prior to the final one, (with a single 4 left in the matrix) and starts the recursion over again.

     0     0     0     0     0     0     0
     0     0     1     2     1     0     0
     0     1     0     3     0     1     0
     0     2     3     0     2     2     0
     0     1     0     2     4     0     0
     0     0     1     2     0     0     0
     0     0     0     0     0     0     0

Is my break in the wrong spot? Any help would be appreciated.

[jgg]

I'm not familiar with the algorithm, so I can't say what it should be doing, but this

 abelian(sandmatrix)

Looks wrong. You're not assigning the recursive call to anything here, so you never "pop" out of the recursion.

[Eric Faulk]

Shouldn't the if-else statement with the break be able to handle that?

[jgg]

Try walking through your code using the keyboard command and see.

[Eric Faulk]

I have, and it gets to the right answer, iterates through the rest of the elements in the matrix, then jumps up a level and reverts one step back. The sandmatrix is solved inside the first if statement, but it won't break out of the loop upon completion.

[jgg]

Do you have an example of a starting sandmatrix so I can try some things?

[Eric Faulk]

Sure. Here's one that finishes.

sandmat = zeros(7)
sandmat(4,4) = 27

And one that doesn't finish:

sandmat = zeros(7)
sandmat(4,4) = 28

[jgg]

Does this do what you intend for the second case:

 function[thematrix] = abelian(sandmatrix)

 global val;

 for c = 1:size(sandmatrix,2)
  for n = 1:size(sandmatrix,1)
      if sandmatrix(n,c) >= 4
          sandmatrix(n,c) = sandmatrix(n,c)-4;
          sandmatrix(n+1,c) = sandmatrix(n+1,c) + 1;
          sandmatrix(n,c+1) = sandmatrix(n,c+1) + 1;
          sandmatrix(n-1,c) = sandmatrix(n-1,c) + 1;
          sandmatrix(n,c-1) = sandmatrix(n,c-1) + 1;
          if any(max(sandmatrix) >= 4)
              abelian(sandmatrix);

          else

              val = sandmatrix;
              break;

          end
      end
  end

thematrix = val;
end

There must be a more efficient way to do this though. Hopefully someone else has a suggestion.

Answer 1

[Geoff Hayes]

Eric - since this algorithm is an example of a cellular automaton, it may be that for each iteration (or generation) you will want to take a "snapshot" of all those cells that have four or more grains (like Walter indicated in his answer) and then perform the grain distribution to each of its neighbouring cells.

On subsequent iterations, we would execute the same query and so "capture" all those cells who, when receiving a grain of sand, now have four grains and so are required (by the rules of the algorithm) to distribute their four grains of sand appropriately. In this manner, you can avoid the recursion (which may be problematic when a cell has thousands of grains of sand). Something like the following will work

 function [sandmatrix] = abelian(sandmatrix)
    close all;
    GRAIN_THRESHOLD = 4;
    % query for all cells that meet or exceed the threshold
    [rowIdcs,colIdcs] = find(sandmatrix>=GRAIN_THRESHOLD);
    while ~isempty(rowIdcs)
        for k=1:length(rowIdcs)
            % remove the GRAIN_THRESHOLD grains
            sandmatrix(rowIdcs(k),colIdcs(k)) = sandmatrix(rowIdcs(k),colIdcs(k)) - GRAIN_THRESHOLD;
            % distribute the grains to neighbouring cells
            sandmatrix = addGrainToCell(sandmatrix,rowIdcs(k)+1,colIdcs(k));
            sandmatrix = addGrainToCell(sandmatrix,rowIdcs(k),colIdcs(k)+1);
            sandmatrix = addGrainToCell(sandmatrix,rowIdcs(k)-1,colIdcs(k));
            sandmatrix = addGrainToCell(sandmatrix,rowIdcs(k),colIdcs(k)-1);
        end
        % display the current state of the model
        image(sandmatrix);
        pause(0.1);
        % query for all cells that meet or exceed the threshold
        [rowIdcs,colIdcs] = find(sandmatrix>=GRAIN_THRESHOLD);
    end  
 end

 function [sandmatrix] = addGrainToCell(sandmatrix,r,c)
    if r < 1 || r > size(sandmatrix,1) || c < 1 || c > size(sandmatrix,2)
        return;
    end
    sandmatrix(r,c) = sandmatrix(r,c) + 1;
 end

[Eric Faulk]

Geoff--your code works beautifully, thank you. The graphic is a nice touch. I suppose the while loop is really a more efficient recursive mechanism than a recursive function call.

Answer 2

[Walter Roberson]

function [sandmatrix] = abelian(sandmatrix)
for c = 1:size(sandmatrix,2)
  for n = 1:size(sandmatrix,1)
      if sandmatrix(n,c) >= 4
          sandmatrix(n,c) = sandmatrix(n,c)-4
          sandmatrix(n+1,c) = sandmatrix(n+1,c) + 1
          sandmatrix(n,c+1) = sandmatrix(n,c+1) + 1
          sandmatrix(n-1,c) = sandmatrix(n-1,c) + 1
          sandmatrix(n,c-1) = sandmatrix(n,c-1) + 1
          if any(max(sandmatrix) >= 4)
              sandmatrix = abelian(sandmatrix);
          else
              break
          end
      end
  end
end

I am not convinced this should be looping over c and n though. I suspect it should be searching for a point that is >= 4 rather than looping. I could be wrong, though, as I am not familiar with this algorithm.