Abelian sandpile model
You are encouraged to solve this task according to the task description, using any language you may know.
| 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
AArch64 Assembly
<lang AArch64 Assembly> /* ARM assembly AARCH64 Raspberry PI 3B or android 64 bits */ /* program abelian64.s */
/* run : abelian 256 12 12 */
/*******************************************/ /* Constantes file */ /*******************************************/ /* for this file see task include a file in language AArch64 assembly*/ .include "../includeConstantesARM64.inc" .equ MAXI, 25
/*********************************/ /* Initialized data */ /*********************************/ .data szMessValue: .asciz "@ " szMessErrParam: .asciz "error : command line = abelian size posx posy \n" szMessFin: .asciz "End display :\n" szCarriageReturn: .asciz "\n"
/*********************************/ /* UnInitialized data */ /*********************************/ .bss sZoneConv: .skip 24 iSandPile: .skip 8 * MAXI * MAXI /*********************************/ /* code section */ /*********************************/ .text .global main main: // entry of program
mov fp,sp ldr x4,[fp] // load number of parameters command line cmp x4,#3 // < 4 -> error ble 99f add x0,fp,#32 // load address param 4 = pos y ldr x0,[x0] bl conversionAtoD // conversion ascii -> numeric mov x3,x0 add x0,fp,#24 // load address param 3 = pos x ldr x0,[x0] bl conversionAtoD mov x2,x0 add x0,fp,#16 // load address param 2 = size begin pile ldr x0,[x0] bl conversionAtoD ldr x4,qAdriSandPile mov x5,#MAXI madd x5,x3,x5,x2 // compute offset = maxi * y + x str x0,[x4,x5,lsl #3] // and store size in pos x,y //mov x0,x4 // display start position //bl displaySandPile mov x0,x4 // sandpile address mov x1,x2 // pos x to start mov x2,x3 // pos y to start bl addSand ldr x0,qAdrszMessFin bl affichageMess mov x0,x4 bl displaySandPile b 100f
99: // line command error
ldr x0,qAdrszMessErrParam bl affichageMess
100: // standard end of the program
mov x0,0 // return code mov x8,EXIT // request to exit program svc 0 // perform the system call
qAdrszCarriageReturn: .quad szCarriageReturn qAdrsZoneConv: .quad sZoneConv qAdrszMessErrParam: .quad szMessErrParam qAdrszMessFin: .quad szMessFin qAdriSandPile: .quad iSandPile /***************************************************/ /* display sandpile */ /***************************************************/ // x0 contains address to sandpile displaySandPile:
stp x1,lr,[sp,-16]! // save registres stp x2,x3,[sp,-16]! // save registres stp x4,x5,[sp,-16]! // save registres stp x6,x7,[sp,-16]! // save registres mov x6,x0 mov x3,#0 // indice y mov x4,#MAXI
1:
mov x2,#0 // indice x
2:
madd x5,x3,x4,x2 // compute offset ldr x0,[x6,x5,lsl #3] // load value at pos x,y ldr x1,qAdrsZoneConv bl conversion10 // call decimal conversion add x1,x1,1 mov x7,#0 strb w7,[x1,x0] ldr x0,qAdrszMessValue ldr x1,qAdrsZoneConv // insert value conversion in message bl strInsertAtCharInc bl affichageMess add x2,x2,1 cmp x2,MAXI blt 2b ldr x0,qAdrszCarriageReturn bl affichageMess add x3,x3,1 cmp x3,MAXI blt 1b
100:
ldp x6,x7,[sp],16 // restaur des 2 registres ldp x4,x5,[sp],16 // restaur des 2 registres ldp x2,x3,[sp],16 // restaur des 2 registres ldp x1,lr,[sp],16 // restaur des 2 registres ret
qAdrszMessValue: .quad szMessValue /***************************************************/ /* display sandpile */ /***************************************************/ // x0 contains address to sanspile // x1 contains position x // x2 contains position y addSand:
stp x1,lr,[sp,-16]! // save registres stp x2,x3,[sp,-16]! // save registres stp x4,x5,[sp,-16]! // save registres mov x3,#MAXI madd x4,x3,x2,x1 // compute offset ldr x5,[x0,x4,lsl #3]
1:
cmp x5,#4 // 4 grains ? blt 100f sub x5,x5,4 // yes sustract str x5,[x0,x4,lsl #3] cmp x1,MAXI-1 // right position ok ? beq 2f add x1,x1,1 // yes bl add1Sand // add 1 grain bl addSand // and compute new pile sub x1,x1,1
2:
cmp x1,0 // left position ok ? beq 3f sub x1,x1,1 bl add1Sand bl addSand add x1,x1,1
3:
cmp x2,0 // higt position ok ? beq 4f sub x2,x2,1 bl add1Sand bl addSand add x2,x2,1
4:
cmp x2,MAXI-1 // low position ok ? beq 5f add x2,x2,1 bl add1Sand bl addSand sub x2,x2,1
5:
ldr x5,[x0,x4,lsl #3] // reload value b 1b // and loop
100:
ldp x4,x5,[sp],16 // restaur des 2 registres ldp x2,x3,[sp],16 // restaur des 2 registres ldp x1,lr,[sp],16 // restaur des 2 registres ret
/***************************************************/ /* add 1 grain of sand */ /***************************************************/ // x0 contains address to sanspile // x1 contains position x // x2 contains position y add1Sand:
stp x3,lr,[sp,-16]! // save registres stp x4,x5,[sp,-16]! // save registres mov x3,#MAXI madd x4,x3,x2,x1 // compute offset ldr x5,[x0,x4,lsl #3] // load value at pos x,y add x5,x5,1 str x5,[x0,x4,lsl #3] // and store
100:
ldp x4,x5,[sp],16 // restaur des 2 registres ldp x3,lr,[sp],16 // restaur des 2 registres ret
/********************************************************/ /* File Include fonctions */ /********************************************************/ /* for this file see task include a file in language AArch64 assembly */ .include "../includeARM64.inc" </lang>
- Output:
~/.../rosetta/asm1 $ abelian64 64 12 12 End display : 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 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 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 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 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 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 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 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 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 0 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 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 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 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 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 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 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 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 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 0 0 0 0
ARM Assembly
<lang ARM Assembly> /* ARM assembly Raspberry PI or android 32 bits */ /* program abelian.s */
/* run : abelian 256 12 12 */
/* REMARK 1 : this program use routines in a include file
see task Include a file language arm assembly for the routine affichageMess conversion10 see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */ /************************************/ /* Constantes */ /************************************/ .include "../constantes.inc" .equ MAXI, 25
/*********************************/ /* Initialized data */ /*********************************/ .data szMessValue: .asciz "@ " szMessErrParam: .asciz "error : command line = abelian size posx posy \n" szMessFin: .asciz "End display :\n" szCarriageReturn: .asciz "\n"
/*********************************/ /* UnInitialized data */ /*********************************/ .bss sZoneConv: .skip 24 iSandPile: .skip 4 * MAXI * MAXI /*********************************/ /* code section */ /*********************************/ .text .global main main: @ entry of program
mov fp,sp ldr r4,[fp] @ load number of parameters commend line cmp r4,#3 @ < 4 -> error ble 99f add r0,fp,#16 @ load address param 4 = pos y ldr r0,[r0] bl conversionAtoD @ conversion ascii -> numeric mov r3,r0 add r0,fp,#12 @ load address param 3 = pos x ldr r0,[r0] bl conversionAtoD mov r2,r0 add r0,fp,#8 @ load address param 2 = size begin pile ldr r0,[r0] bl conversionAtoD ldr r4,iAdriSandPile mov r5,#MAXI mul r5,r3,r5 @ compute offset = maxi * y add r5,r2 @ + x str r0,[r4,r5,lsl #2] @ and store size in pos x,y //mov r0,r4 @ display start position //bl displaySandPile mov r0,r4 @ sandpile address mov r1,r2 @ pos x to start mov r2,r3 @ pos y to start bl addSand ldr r0,iAdrszMessFin bl affichageMess mov r0,r4 bl displaySandPile b 100f
99: @ line command error
ldr r0,iAdrszMessErrParam bl affichageMess
100: @ standard end of the program
mov r0, #0 @ return code mov r7, #EXIT @ request to exit program svc #0 @ perform the system call
iAdrszCarriageReturn: .int szCarriageReturn iAdrsZoneConv: .int sZoneConv iAdrszMessErrParam: .int szMessErrParam iAdrszMessFin: .int szMessFin iAdriSandPile: .int iSandPile /***************************************************/ /* display sandpile */ /***************************************************/ // r0 contains address to sandpile displaySandPile:
push {r1-r6,lr} @ save registers
mov r6,r0
mov r3,#0 @ indice y
mov r4,#MAXI
1:
mov r2,#0 @ indice x
2:
mul r5,r3,r4 add r5,r2 @ compute offset ldr r0,[r6,r5,lsl #2] @ load value at pos x,y ldr r1,iAdrsZoneConv bl conversion10 @ call decimal conversion add r1,#1 mov r7,#0 strb r7,[r1,r0] ldr r0,iAdrszMessValue ldr r1,iAdrsZoneConv @ insert value conversion in message bl strInsertAtCharInc bl affichageMess add r2,#1 cmp r2,#MAXI blt 2b ldr r0,iAdrszCarriageReturn bl affichageMess add r3,#1 cmp r3,#MAXI blt 1b
100:
pop {r1-r6,lr} @ restaur registers
bx lr @ return
iAdrszMessValue: .int szMessValue /***************************************************/ /* display sandpile */ /***************************************************/ // r0 contains address to sanspile // r1 contains position x // r2 contains position y addSand:
push {r1-r5,lr} @ save registers
mov r3,#MAXI
mul r4,r3,r2
add r4,r1
ldr r5,[r0,r4,lsl #2]
1:
cmp r5,#4 @ 4 grains ? blt 100f sub r5,#4 @ yes sustract str r5,[r0,r4,lsl #2] cmp r1,#MAXI-1 @ right position ok ? beq 2f add r1,#1 @ yes bl add1Sand @ add 1 grain bl addSand @ and compute new pile sub r1,#1
2:
cmp r1,#0 @ left position ok ? beq 3f sub r1,#1 bl add1Sand bl addSand add r1,#1
3:
cmp r2,#0 @ higt position ok ? beq 4f sub r2,#1 bl add1Sand bl addSand add r2,#1
4:
cmp r2,#MAXI-1 @ low position ok ? beq 5f add r2,#1 bl add1Sand bl addSand sub r2,#1
5:
ldr r5,[r0,r4,lsl #2] @ reload value b 1b @ and loop
100:
pop {r1-r5,lr} @ restaur registers
bx lr @ return
/***************************************************/ /* add 1 grain of sand */ /***************************************************/ // r0 contains address to sanspile // r1 contains position x // r2 contains position y add1Sand:
push {r3-r5,lr} @ save registers
mov r3,#MAXI
mul r4,r3,r2
add r4,r1 @ compute offset
ldr r5,[r0,r4,lsl #2] @ load value at pos x,y
add r5,#1
str r5,[r0,r4,lsl #2] @ and store
100:
pop {r3-r5,lr} @ restaur registers
bx lr @ return
/***************************************************/ /* ROUTINES INCLUDE */ /***************************************************/ .include "../affichage.inc" </lang>
- Output:
pi@raspberrypi:~/asm32/rosetta32/ass10 $ abelian 512 12 12 End display : 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 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 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 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 0 2 2 2 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 2 3 2 3 2 3 2 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 3 0 2 0 3 0 2 0 3 3 2 0 0 0 0 0 0 0 0 0 0 0 2 3 0 3 3 3 3 3 3 3 3 3 0 3 2 0 0 0 0 0 0 0 0 0 1 1 3 3 0 3 1 3 3 3 1 3 0 3 3 1 1 0 0 0 0 0 0 0 0 3 2 0 3 3 2 2 0 3 0 2 2 3 3 0 2 3 0 0 0 0 0 0 0 1 0 3 2 3 1 2 2 2 3 2 2 2 1 3 2 3 0 1 0 0 0 0 0 0 1 2 2 0 3 3 0 2 0 3 0 2 0 3 3 0 2 2 1 0 0 0 0 0 0 1 2 3 3 3 3 3 3 3 0 3 3 3 3 3 3 3 2 1 0 0 0 0 0 0 1 2 2 0 3 3 0 2 0 3 0 2 0 3 3 0 2 2 1 0 0 0 0 0 0 1 0 3 2 3 1 2 2 2 3 2 2 2 1 3 2 3 0 1 0 0 0 0 0 0 0 3 2 0 3 3 2 2 0 3 0 2 2 3 3 0 2 3 0 0 0 0 0 0 0 0 1 1 3 3 0 3 1 3 3 3 1 3 0 3 3 1 1 0 0 0 0 0 0 0 0 0 2 3 0 3 3 3 3 3 3 3 3 3 0 3 2 0 0 0 0 0 0 0 0 0 0 0 2 3 3 0 2 0 3 0 2 0 3 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 2 3 2 3 2 3 2 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 0 2 2 2 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 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 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 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
C
Writes out the initial and final sand piles to the console and the final sand pile to a PPM file. <lang C>
- include<stdlib.h>
- include<string.h>
- include<stdio.h>
int main(int argc, char** argv) { int i,j,sandPileEdge, centerPileHeight, processAgain = 1,top,down,left,right; int** sandPile; char* fileName; static unsigned char colour[3];
if(argc!=3){
printf("Usage: %s <Sand pile side>
return 0; }
sandPileEdge = atoi(argv[1]); centerPileHeight = atoi(argv[2]);
if(sandPileEdge<=0 || centerPileHeight<=0){ printf("Sand pile and center pile dimensions must be positive integers."); return 0; }
sandPile = (int**)malloc(sandPileEdge * sizeof(int*));
for(i=0;i<sandPileEdge;i++){ sandPile[i] = (int*)calloc(sandPileEdge,sizeof(int)); }
sandPile[sandPileEdge/2][sandPileEdge/2] = centerPileHeight;
printf("Initial sand pile :\n\n");
for(i=0;i<sandPileEdge;i++){ for(j=0;j<sandPileEdge;j++){ printf("%3d",sandPile[i][j]); } printf("\n"); }
while(processAgain == 1){
processAgain = 0; top = 0; down = 0; left = 0; right = 0;
for(i=0;i<sandPileEdge;i++){ for(j=0;j<sandPileEdge;j++){ if(sandPile[i][j]>=4){ if(i-1>=0){ top = 1; sandPile[i-1][j]+=1; if(sandPile[i-1][j]>=4) processAgain = 1; } if(i+1<sandPileEdge){ down = 1; sandPile[i+1][j]+=1; if(sandPile[i+1][j]>=4) processAgain = 1; } if(j-1>=0){ left = 1; sandPile[i][j-1]+=1; if(sandPile[i][j-1]>=4) processAgain = 1; } if(j+1<sandPileEdge){ right = 1; sandPile[i][j+1]+=1; if(sandPile[i][j+1]>=4) processAgain = 1; } sandPile[i][j] -= (top + down + left + right); if(sandPile[i][j]>=4) processAgain = 1; } } } }
printf("Final sand pile : \n\n");
for(i=0;i<sandPileEdge;i++){ for(j=0;j<sandPileEdge;j++){ printf("%3d",sandPile[i][j]); } printf("\n"); }
fileName = (char*)malloc((strlen(argv[1]) + strlen(argv[2]) + 23)*sizeof(char));
strcpy(fileName,"Final_Sand_Pile_"); strcat(fileName,argv[1]); strcat(fileName,"_"); strcat(fileName,argv[2]); strcat(fileName,".ppm");
FILE *fp = fopen(fileName,"wb");
fprintf(fp,"P6\n%d %d\n255\n",sandPileEdge,sandPileEdge);
for(i=0;i<sandPileEdge;i++){ for(j=0;j<sandPileEdge;j++){ colour[0] = (sandPile[i][j] + i)%256; colour[1] = (sandPile[i][j] + j)%256; colour[2] = (sandPile[i][j] + i*j)%256; fwrite(colour,1,3,fp); } }
fclose(fp);
printf("\nImage file written to %s\n",fileName);
return 0; } </lang>
Console output :
abhishek_ghosh@Azure:~/doodles$ ./a.out 10 64 Initial sand pile : 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 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 64 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Final sand pile : 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 2 1 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 1 2 2 2 2 2 1 0 0 0 2 2 2 0 2 2 2 0 0 0 1 2 2 2 2 2 1 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 Image file written to Final_Sand_Pile_10_64.ppm
C++
<lang cpp>#include <iostream>
- include "xtensor/xarray.hpp"
- include "xtensor/xio.hpp"
- include "xtensor-io/ximage.hpp"
xt::xarray<int> init_grid (unsigned long x_dim, unsigned long y_dim) {
xt::xarray<int>::shape_type shape = { x_dim, y_dim };
xt::xarray<int> grid(shape);
grid(x_dim/2, y_dim/2) = 64000; return grid;
}
int print_grid(xt::xarray<int>& grid) {
// for output to the terminal uncomment next line // only makes sense for small grid < 32x32; // std::cout << grid << std::endl << std::endl;
// output result to an image
xt::dump_image("grid.jpg", grid);
return 0;
}
bool iterate_grid(xt::xarray<int>& grid, const unsigned long& x_dim, const unsigned long& y_dim) {
bool changed = false;
for (unsigned long i=0; i < x_dim; ++i)
{
for (unsigned long j=0; j < y_dim; ++j)
{
if ( grid(i, j) >= 4 )
{
grid(i, j) -= 4;
changed = true;
try
{
grid.at(i-1, j) += 1;
grid.at(i+1, j) += 1;
grid.at(i, j-1) += 1;
grid.at(i, j+1) += 1;
}
catch (const std::out_of_range& oor)
{
}
}
}
}
return changed;
}
int main(int argc, char* argv[]) {
const unsigned long x_dim { 200 };
const unsigned long y_dim { 200 };
xt::xarray<int> grid = init_grid(x_dim, y_dim);
bool changed { true };
iterate_grid(grid, x_dim, y_dim);
while (changed == true)
{
changed = iterate_grid(grid, x_dim, y_dim);
}
print_grid(grid);
return 0;
}</lang> Compile with following CMakeList.txt: <lang cmake>cmake_minimum_required(VERSION 3.1) project(abelian_sandpile)
find_package(xtl REQUIRED) find_package(xtensor REQUIRED)
- if xtensor was built with xsimd support:
- find_package(xsimd REQUIRED)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") include_directories(/usr/include/OpenImageIO) find_library(OIIO "OpenImageIO")
add_executable(abelian_sandpile src/abelian_sandpile.cpp)
target_compile_options(abelian_sandpile PRIVATE -march=native -std=c++14) target_link_libraries(abelian_sandpile xtensor ${OIIO})</lang>
Delphi
<lang Delphi> program Abelian_sandpile_model;
{$APPTYPE CONSOLE}
{$R *.res}
uses
System.SysUtils, Vcl.Graphics, System.Classes;
type
TGrid = array of array of Integer;
function Iterate(var Grid: TGrid): Boolean; var
changed: Boolean; i: Integer; j: Integer; val: Integer; Alength: Integer;
begin
Alength := length(Grid); changed := False;
for i := 0 to High(Grid) do
for j := 0 to High(Grid[0]) do
begin
val := Grid[i, j];
if val > 3 then
begin
Grid[i, j] := Grid[i, j] - 4;
if i > 0 then
Grid[i - 1, j] := Grid[i - 1, j] + 1;
if i < Alength - 1 then
Grid[i + 1, j] := Grid[i + 1, j] + 1;
if j > 0 then
Grid[i, j - 1] := Grid[i, j - 1] + 1;
if j < Alength - 1 then
Grid[i, j + 1] := Grid[i, j + 1] + 1;
changed := True;
end;
end;
Result := changed;
end;
procedure Simulate(var Grid: TGrid); var
changed: Boolean;
begin
while Iterate(Grid) do ;
end;
procedure Zeros(var Grid: TGrid; Size: Integer); var
i, j: Integer;
begin
SetLength(Grid, Size, Size);
for i := 0 to Size - 1 do
for j := 0 to Size - 1 do
Grid[i, j] := 0;
end;
procedure Println(Grid: TGrid); var
i, j: Integer;
begin
for i := 0 to High(Grid) do
begin
Writeln;
for j := 0 to High(Grid[0]) do
Write(Format('%3d', [Grid[i, j]]));
end;
Writeln;
end;
function Grid2Bmp(Grid: TGrid): TBitmap; const
Colors: array[0..2] of TColor = (clRed, clLime, clBlue);
var
Alength: Integer; i: Integer; j: Integer;
begin
Alength := Length(Grid);
Result := TBitmap.Create; Result.SetSize(Alength, Alength);
for i := 0 to Alength - 1 do
for j := 0 to Alength - 1 do
begin
Result.Canvas.Pixels[i, j] := Colors[Grid[i, j]];
end;
end;
procedure Grid2P6(Grid: TGrid; FileName: TFileName); var
f: text; i, j, Alength: Integer; ppm: TFileStream; Header: AnsiString;
const
COLORS: array[0..3] of array[0..2] of byte = // R, G, B
((0 , 0, 0),
(255 , 0, 0), (0 , 255, 0), (0 , 0, 255));
begin
Alength := Length(Grid);
ppm := TFileStream.Create(FileName, fmCreate);
Header := Format('P6'#10'%d %d'#10'255'#10, [Alength, Alength]);
writeln(Header);
ppm.Write(Tbytes(Header), Length(Header));
for i := 0 to Alength - 1 do
for j := 0 to Alength - 1 do
begin
ppm.Write(COLORS[Grid[i, j]], 3);
end;
ppm.Free;
end;
const
DIMENSION = 10;
var
Grid: TGrid; bmp: TBitmap;
begin
Zeros(Grid, DIMENSION);
Grid[4, 4] := 64;
Writeln('Before:');
Println(Grid);
Simulate(Grid);
Writeln(#10'After:'); Println(Grid);
// Output bmp
with Grid2Bmp(Grid) do
begin
SaveToFile('output.bmp');
free;
end;
// Output ppm Grid2P6(Grid, 'output.ppm');
Readln;
end. </lang>
Forth
<lang forth>#! /usr/bin/gforth -d 20M
\ Abelian Sandpile Model
0 assert-level !
\ command-line
- parse-number s>number? invert throw drop ;
- parse-size ." size : " next-arg parse-number dup . cr ;
- parse-height ." height: " next-arg parse-number dup . cr ;
- parse-args cr parse-size parse-height ;
parse-args constant HEIGHT constant SIZE
- allot-erase create here >r dup allot r> swap erase ;
- size^2 SIZE dup * cells ;
- 2cells [ 2 cells ] literal ;
- -2cells [ 2cells negate ] literal ;
size^2 allot-erase arr
\ array processing
- ix swap SIZE * + cells arr + ;
- center SIZE 2/ dup ;
- write-cell ix @ u. ;
- write-row SIZE 0 ?do dup i write-cell loop drop cr ;
- arr. SIZE 0 ?do i write-row loop ;
\ stack processing
- stack-empty? dup -1 = ;
- stack-full? stack-empty? invert ;
\ pgm-handling
- concat { a1 l1 a2 l2 } l1 l2 + allocate throw dup dup a1 swap l1 cmove a2 swap l1 + l2 cmove l1 l2 + ;
- write-pgm ." P2" cr SIZE u. SIZE u. cr ." 3" cr arr. ;
- u>s 0 <# #s #> ;
- filename s" sandpile-" SIZE u>s concat s" -" concat HEIGHT u>s concat s" .pgm" concat ;
- to-pgm filename w/o create-file throw ['] write-pgm over outfile-execute close-file throw ;
\ sandpile
- prep-arr HEIGHT center ix ! ;
- prep-stack -1 HEIGHT 4 u>= if center then ;
- prepare prep-arr prep-stack ;
- ensure if else 2drop 0 2rdrop exit then ;
- col>=0 dup 0>= ensure ;
- col<SIZE dup SIZE < ensure ;
- row>=0 over 0>= ensure ;
- row<SIZE over SIZE < ensure ;
- legal? col>=0 col<SIZE row>=0 row<SIZE 2drop true ;
- north 1. d- ;
- east 1+ ;
- south 1. d+ ;
- west 1- ;
- reduce 2dup ix dup -4 swap +! @ 4 < if 2drop then ;
- increase 2dup legal? if 2dup ix dup 1 swap +! @ 4 = if 2swap else 2drop then else 2drop then ;
- inc-north 2dup north increase ;
- inc-east 2dup east increase ;
- inc-south 2dup south increase ;
- inc-west 2dup west increase ;
- inc-all inc-north inc-east inc-south inc-west 2drop ;
- simulate prepare begin stack-full? while 2dup 2>r reduce 2r> inc-all repeat drop to-pgm ." written to " filename type cr ;
simulate bye</lang>
- Output:
sandpile with 5000 grains of sand:
./sandpile.fs 61 5000:
[1]
sandpile with 50000 grains of sand:
./sandpile.fs 201 50000:
[2]
sandpile with 500000 grains of sand:
./sandpile.fs 601 500000:
[3]
![[1]](http://commons.wikimedia.org/wiki/File:Sandpile-61-5000.png){kind=link}
![[2]](http://commons.wikimedia.org/wiki/File:Sandpile-201-50000.png){kind=link}
![[3]](http://commons.wikimedia.org/wiki/File:Sandpile-601-500000.png){kind=link}
Fortran
The Abelian sandpile operations are defined here. <lang fortran>module abelian_sandpile_m
implicit none
private public :: pile
type :: pile !! usage: !! 1) init !! 2) run
integer, allocatable :: grid(:,:) integer :: n(2)
contains procedure :: init procedure :: run
procedure, private :: process_node procedure, private :: inside end type
contains
logical function inside(this, i) class(pile), intent(in) :: this integer, intent(in) :: i(2)
inside = ((i(1) > 0) .and. (i(1) <= this%n(1)) .and. (i(2) > 0) .and. (i(2) <= this%n(2)) ) end function
recursive subroutine process_node(this, i) !! start process
class(pile), intent(inout) :: this
integer, intent(in) :: i(2)
!! node coordinates to process
integer :: i0(2,2), j(2), d, k
! if node has more than 4 grains -> redistribute
if (this%grid(i(1),i(2)) >= 4) then
! unit vectors: help shift only one dimension (see below)
i0 = reshape([1,0,0,1], [2,2])
! subtract 4 grains
this%grid(i(1),i(2)) = this%grid(i(1),i(2))-4
! add one grain to neighbor if not out of bound
do d = 1, 2 ! loop dimensions
do k = -1, 1, 2 ! loop +-1 step in direction d
j = i+k*i0(:,d) ! j = i, but one element is shifted by +-1
if (this%inside(j)) this%grid(j(1),j(2)) = this%grid(j(1),j(2)) + 1
end do
end do
! check neighbor nodes
do d = 1, 2 ! loop dimensions
do k = -1, 1, 2 ! loop +-1 step in direction d
j = i+k*i0(:,d) ! j = i, but one element is shifted by +-1
if (this%inside(j)) call this%process_node(j)
end do
end do
! check itself
call this%process_node(i)
end if
end subroutine
subroutine run(this) !! start process
class(pile), intent(inout) :: this
! only node that could be unstable is inital node call this%process_node(this%n/2) end subroutine
subroutine init(this, nx, ny, h)
class(pile), intent(out) :: this
integer, intent(in) :: nx, ny
!! grid dimensions
integer, intent(in) :: h
!! height of and grains in middle of grid
this%n = [nx, ny] allocate (this%grid(nx,ny), source=0) this%grid(nx/2, ny/2) = h end subroutine
end module</lang>
The main program calls the abelian_sandpile_m and creates an ppm bitmap file by loading rgbimage_m module, which is defined here.
<lang fortran>program main
use rgbimage_m use abelian_sandpile_m
implicit none
integer :: nx, ny, i, j
integer :: colors(0:3,3)
type(rgbimage) :: im type(pile) :: p
colors(0,:) = [255,255,255] colors(1,:) = [0,0,90] colors(2,:) = [0,0,170] colors(3,:) = [0,0,255]
nx = 200 ny = 100
call p%init(nx, ny, 2000) call p%run
call im%init(nx, ny)
do i = 1, nx
do j = 1, ny
call im%set_pixel(i, j, colors(p%grid(i,j),:))
end do
end do
call im%write('fig.ppm')
end program</lang>
Fōrmulæ
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website, However they run on execution servers. By default remote servers are used, but they are limited in memory and processing power, since they are intended for demonstration and casual use. A local server can be downloaded and installed, it has no limitations (it runs in your own computer). Because of that, example programs can be fully visualized and edited, but some of them will not run if they require a moderate or heavy computation/memory resources, and no local server is being used.
In this page you can see the program(s) related to this task and their results.
F#
<lang fsharp> // Abelian sandpile model. Nigel Galloway: July 20th., 2020 type Sandpile(x,y,N:int[])=
member private this.x=x
member private this.y=y
member private this.i=let rec topple n=match Array.tryFindIndex(fun n->n>3)n with
None->n
|Some g->let i=n.[g]/4
n.[g]<-n.[g]%4
match g%x,g/x with
(0,0)->n.[x]<-n.[x]+i;n.[1]<-n.[1]+i;topple n
|(α,0) when α=x-1->n.[g+x]<-n.[g+x]+i;n.[g-1]<-n.[g-1]+i;topple n
|(_,0)->n.[g-1]<-n.[g-1]+i;n.[g+1]<-n.[g+1]+i;n.[g+x]<-n.[g+x]+i;topple n
|(0,β) when β=y-1->n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;topple n
|(0,β)->n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;n.[g+x]<-n.[g+x]+i;topple n
|(α,β) when α=x-1 && β=y-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;topple n
|(α,_) when α=x-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+x]<-n.[g+x]+i;topple n
|(_,β) when β=y-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;topple n
|_->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+x]<-n.[g+x]+i;n.[g+1]<-n.[g+1]+i;topple n
topple N
static member (+) (n:Sandpile, g:Sandpile)=Sandpile(n.x,n.y,Array.map2(fun n g->n+g) n.i g.i)
member this.toS=sprintf "%A" (this.i|>Array.chunkBySize x|>array2D)
printfn "%s\n" (Sandpile(3,3,[|4;3;3;3;1;2;0;2;3|])).toS let e1=Array.zeroCreate<int> 25 in e1.[12]<-4; printfn "%s\n" (Sandpile(5,5,e1)).toS let e1=Array.zeroCreate<int> 25 in e1.[12]<-6; printfn "%s\n" (Sandpile(5,5,e1)).toS let e1=Array.zeroCreate<int> 25 in e1.[12]<-16; printfn "%s\n" (Sandpile(5,5,e1)).toS </lang>
- Output:
[[2; 1; 0] [0; 3; 3] [1; 2; 3]] [[0; 0; 0; 0; 0] [0; 0; 1; 0; 0] [0; 1; 0; 1; 0] [0; 0; 1; 0; 0] [0; 0; 0; 0; 0]] [[0; 0; 0; 0; 0] [0; 0; 1; 0; 0] [0; 1; 2; 1; 0] [0; 0; 1; 0; 0] [0; 0; 0; 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]]
Go
Stack management in Go is automatic, starting very small (2KB) for each goroutine and expanding as necessary until the maximum allowed size is reached.
<lang go>package main
import (
"fmt" "log" "os" "strings"
)
const dim = 16 // image size
func check(err error) {
if err != nil {
log.Fatal(err)
}
}
// Outputs the result to the terminal using UTF-8 block characters. func drawPile(pile [][]uint) {
chars:= []rune(" ░▓█")
for _, row := range pile {
line := make([]rune, len(row))
for i, elem := range row {
if elem > 3 { // only possible when algorithm not yet completed.
elem = 3
}
line[i] = chars[elem]
}
fmt.Println(string(line))
}
}
// Creates a .ppm file in the current directory, which contains // a colored image of the pile. func writePile(pile [][]uint) {
file, err := os.Create("output.ppm")
check(err)
defer file.Close()
// Write the signature, image dimensions and maximum color value to the file.
fmt.Fprintf(file, "P3\n%d %d\n255\n", dim, dim)
bcolors := []string{"125 0 25 ", "125 80 0 ", "186 118 0 ", "224 142 0 "}
var line strings.Builder
for _, row := range pile {
for _, elem := range row {
line.WriteString(bcolors[elem])
}
file.WriteString(line.String() + "\n")
line.Reset()
}
}
// Main part of the algorithm, a simple, recursive implementation of the model. func handlePile(x, y uint, pile [][]uint) {
if pile[y][x] >= 4 {
pile[y][x] -= 4
// Check each neighbor, whether they have enough "sand" to collapse and if they do,
// recursively call handlePile on them.
if y > 0 {
pile[y-1][x]++
if pile[y-1][x] >= 4 {
handlePile(x, y-1, pile)
}
}
if x > 0 {
pile[y][x-1]++
if pile[y][x-1] >= 4 {
handlePile(x-1, y, pile)
}
}
if y < dim-1 {
pile[y+1][x]++
if pile[y+1][x] >= 4 {
handlePile(x, y+1, pile)
}
}
if x < dim-1 {
pile[y][x+1]++
if pile[y][x+1] >= 4 {
handlePile(x+1, y, pile)
}
}
// Uncomment this line to show every iteration of the program.
// Not recommended with large input values.
// drawPile(pile)
// Finally call the function on the current cell again,
// in case it had more than 4 particles.
handlePile(x, y, pile)
}
}
func main() {
// Create 2D grid and set size using the 'dim' constant.
pile := make([][]uint, dim)
for i := 0; i < dim; i++ {
pile[i] = make([]uint, dim)
}
// Place some sand particles in the center of the grid and start the algorithm. hdim := uint(dim/2 - 1) pile[hdim][hdim] = 16 handlePile(hdim, hdim, pile) drawPile(pile)
// Uncomment this to save the final image to a file // after the recursive algorithm has ended. // writePile(pile)
}</lang>
- Output:
░
▓░▓
░░ ░░
▓░▓
░
Haskell
Using a custom monad to make the code cleaner.
<lang haskell>{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Rosetta.AbelianSandpileModel.ST
( simulate , test , toPGM ) where
import Control.Monad.Reader (asks, MonadReader (..), ReaderT, runReaderT) import Control.Monad.ST (runST, ST) import Control.Monad.State (evalStateT, forM_, lift, MonadState (..), StateT, modify, when) import Data.Array.ST (freeze, readArray, STUArray, thaw, writeArray) import Data.Array.Unboxed (array, assocs, bounds, UArray, (!)) import Data.Word (Word32) import System.IO (hPutStr, hPutStrLn, IOMode (WriteMode), withFile) import Text.Printf (printf)
type Point = (Int, Int) type ArrayST s = STUArray s Point Word32 type ArrayU = UArray Point Word32
newtype M s a = M (ReaderT (S s) (StateT [Point] (ST s)) a)
deriving (Functor, Applicative, Monad, MonadReader (S s), MonadState [Point])
data S s = S
{ bMin :: !Point
, bMax :: !Point
, arr :: !(ArrayST s)
}
runM :: M s a -> S s -> [Point]-> ST s a runM (M m) = evalStateT . runReaderT m
liftST :: ST s a -> M s a liftST = M . lift . lift
simulate :: ArrayU -> ArrayU simulate a = runST $ simulateST a
simulateST :: forall s. ArrayU -> ST s ArrayU simulateST a = do
let (p1, p2) = bounds a
s = [p | (p, c) <- assocs a, c >= 4]
b <- thaw a :: ST s (ArrayST s)
let st = S { bMin = p1
, bMax = p2
, arr = b
}
runM simulateM st s
simulateM :: forall s. M s ArrayU simulateM = do
ps <- get
case ps of
[] -> asks arr >>= liftST . freeze
p : ps' -> do
c <- changeArr p $ \x -> x - 4
when (c < 4) $ put ps'
forM_ [north, east, south, west] $ inc . ($ p)
simulateM
changeArr :: Point -> (Word32 -> Word32) -> M s Word32 changeArr p f = do
a <- asks arr oldC <- liftST $ readArray a p let newC = f oldC liftST $ writeArray a p newC return newC
inc :: Point -> M s () inc p = do
b <- inBounds p
when b $ do
c <- changeArr p succ
when (c == 4) $ modify $ (p :)
inBounds :: Point -> M s Bool inBounds p = do
st <- ask return $ p >= bMin st && p <= bMax st
north, east, south, west :: Point -> Point north (x, y) = (x, y + 1) east (x, y) = (x + 1, y) south (x, y) = (x, y - 1) west (x, y) = (x - 1, y)
toPGM :: ArrayU -> FilePath -> IO () toPGM a fp = withFile fp WriteMode $ \h -> do
let ((x1, y1), (x2, y2)) = bounds a
width = x2 - x1 + 1
height = y2 - y1 + 1
hPutStrLn h "P2"
hPutStrLn h $ show width ++ " " ++ show height
hPutStrLn h "3"
forM_ [y1 .. y2] $ \y -> do
forM_ [x1 .. x2] $ \x -> do
let c = min 3 $ a ! (x, y)
hPutStr h $ show c ++ " "
hPutStrLn h ""
initArray :: Int -> Word32 -> ArrayU initArray size height = array
((-size, -size), (size, size)) [((x, y), if x == 0 && y == 0 then height else 0) | x <- [-size .. size], y <- [-size .. size]]
test :: Int -> Word32 -> IO () test size height = do
printf "size = %d, height = %d\n" size height
let a = initArray size height
b = simulate a
fp = printf "sandpile_%d_%d.pgm" size height
toPGM b fp
putStrLn $ "wrote image to " ++ fp</lang>
- Output:
sandpile with 1000 grains of sand: test 15 1000: [4]
sandpile with 10000 grains of sand: test 40 10000: [5]
sandpile with 100000 grains of sand: test 150 100000: [6]
sandpile with 1000000 grains of sand: test 400 1000000: [7]
![[4]](http://commons.wikimedia.org/wiki/File:Sandpile_15_1000.jpg){kind=link}
![[5]](http://commons.wikimedia.org/wiki/File:Sandpile_40_10000.jpg){kind=link}
![[6]](http://commons.wikimedia.org/wiki/File:Sandpile_150_100000.jpg){kind=link}
![[7]](http://commons.wikimedia.org/wiki/File:Sandpile_400_1000000.jpg){kind=link}
J
<lang J>grid=: 4 : 'x (<<.-:2$y)} (2$y)$0' NB. y by y grid with x grains in middle ab=: - [: +/@(-"2 ((,-)=/~i.2)|.!.0]) 3&< NB. abelian sand pile for grid graph require 'viewmat' NB. viewmat utility viewmat ab ^: _ (1024 grid 25) NB. visual </lang>
Java
This is based on the JavaScript implementation linked to in the task description. <lang java>import java.awt.*; import java.awt.event.*; import javax.swing.*;
public class AbelianSandpile {
public static void main(String[] args) {
SwingUtilities.invokeLater(new Runnable() {
public void run() {
Frame frame = new Frame();
frame.setVisible(true);
}
});
}
private static class Frame extends JFrame {
private Frame() {
super("Abelian Sandpile Model");
setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
Container contentPane = getContentPane();
JPanel controlPanel = new JPanel(new FlowLayout(FlowLayout.LEFT));
JButton start = new JButton("Restart Simulation");
start.addActionListener(e -> restartSimulation());
JButton stop = new JButton("Stop Simulation");
stop.addActionListener(e -> stopSimulation());
controlPanel.add(start);
controlPanel.add(stop);
contentPane.add(controlPanel, BorderLayout.NORTH);
contentPane.add(canvas = new Canvas(), BorderLayout.CENTER);
timer = new Timer(100, e -> canvas.runAndDraw());
timer.start();
pack();
}
private void restartSimulation() {
timer.stop();
canvas.initGrid();
timer.start();
}
private void stopSimulation() {
timer.stop();
}
private Timer timer;
private Canvas canvas;
}
private static class Canvas extends JComponent {
private Canvas() {
setBorder(BorderFactory.createEtchedBorder());
setPreferredSize(new Dimension(600, 600));
}
public void paintComponent(Graphics g) {
int width = getWidth();
int height = getHeight();
g.setColor(Color.WHITE);
g.fillRect(0, 0, width, height);
int cellWidth = width/GRID_LENGTH;
int cellHeight = height/GRID_LENGTH;
for (int i = 0; i < GRID_LENGTH; ++i) {
for (int j = 0; j < GRID_LENGTH; ++j) {
if (grid[i][j] > 0) {
g.setColor(COLORS[grid[i][j]]);
g.fillRect(i * cellWidth, j * cellHeight, cellWidth, cellHeight);
}
}
}
}
private void initGrid() {
for (int i = 0; i < GRID_LENGTH; ++i) {
for (int j = 0; j < GRID_LENGTH; ++j) {
grid[i][j] = 0;
}
}
}
private void runAndDraw() {
for (int i = 0; i < 100; ++i)
addSand(GRID_LENGTH/2, GRID_LENGTH/2);
repaint();
}
private void addSand(int i, int j) {
int grains = grid[i][j];
if (grains < 3) {
grid[i][j]++;
}
else {
grid[i][j] = grains - 3;
if (i > 0)
addSand(i - 1, j);
if (i < GRID_LENGTH - 1)
addSand(i + 1, j);
if (j > 0)
addSand(i, j - 1);
if (j < GRID_LENGTH - 1)
addSand(i, j + 1);
}
}
private int[][] grid = new int[GRID_LENGTH][GRID_LENGTH]; }
private static final Color[] COLORS = {
Color.WHITE,
new Color(0x00, 0xbf, 0xff),
new Color(0xff, 0xd7, 0x00),
new Color(0xb0, 0x30, 0x60)
};
private static final int GRID_LENGTH = 300;
}</lang>
- Output:
See: abelian_sandpile.png (offsite PNG image)
Julia
Modified from code by Hayk Aleksanyan, viewable at github.com/hayk314/Sandpiles, license viewable there. <lang julia>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)
</lang>- Output:
Link to PNG output file for N=100000 ie. AbelSand.move(100000)
Link to PNG output file (run time >90 min) for N=1000000 (move(1000000))
Lua
<lang lua>local sandpile = {
init = function(self, dim, val)
self.cell, self.dim = {}, dim
for r = 1, dim do
self.cell[r] = {}
for c = 1, dim do
self.cell[r][c] = 0
end
end
self.cell[math.floor(dim/2)+1][math.floor(dim/2)+1] = val
end,
iter = function(self)
local dim, cel, more = self.dim, self.cell
repeat
more = false
for r = 1, dim do
for c = 1, dim do
if cel[r][c] >= 4 then
cel[r][c] = cel[r][c] - 4
if c > 1 then cel[r][c-1], more = cel[r][c-1]+1, more or cel[r][c-1]>=3 end
if c < dim then cel[r][c+1], more = cel[r][c+1]+1, more or cel[r][c+1]>=3 end
if r > 1 then cel[r-1][c], more = cel[r-1][c]+1, more or cel[r-1][c]>=3 end
if r < dim then cel[r+1][c], more = cel[r+1][c]+1, more or cel[r+1][c]>=3 end
end
more = more or cel[r][c] >= 4
end
end
until not more
end,
draw = function(self)
for r = 1, self.dim do
print(table.concat(self.cell[r]," "))
end
end,
} sandpile:init(15, 256) sandpile:iter() sandpile:draw()</lang>
- Output:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 1 0 0 0 0 0 0 0 0 0 2 1 0 2 0 1 2 0 0 0 0 0 0 0 2 3 3 3 2 3 3 3 2 0 0 0 0 0 2 3 2 2 2 3 2 2 2 3 2 0 0 0 1 1 3 2 2 3 0 3 2 2 3 1 1 0 0 2 0 3 2 3 2 3 2 3 2 3 0 2 0 0 2 2 2 3 0 3 0 3 0 3 2 2 2 0 0 2 0 3 2 3 2 3 2 3 2 3 0 2 0 0 1 1 3 2 2 3 0 3 2 2 3 1 1 0 0 0 2 3 2 2 2 3 2 2 2 3 2 0 0 0 0 0 2 3 3 3 2 3 3 3 2 0 0 0 0 0 0 0 2 1 0 2 0 1 2 0 0 0 0 0 0 0 0 0 1 2 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[sp] sp[s_List] + sp[n_Integer] ^:= sp[s] + sp[ConstantArray[n, Dimensions[s]]] sp[s_List] + sp[t_List] ^:= Module[{dim, r, tmp, neighbours}, dim = Dimensions[s];
r = s + t;
While[Max[r] > 3, r = ArrayPad[r, 1, 0];
tmp = Quotient[r, 4];
r -= 4 tmp;
r += RotateLeft[tmp, {0, 1}] + RotateLeft[tmp, {1, 0}] +
RotateLeft[tmp, {0, -1}] + RotateLeft[tmp, {-1, 0}];
r = ArrayPad[r, -1];];
sp[r]
]
u = sp[CenterArray[250, {15, 15}]]; u += sp[0]; StringRiffle[StringJoin /@ Map[ToString, u1, {2}], "\n"]</lang>
- Output:
000000000000000 000001222100000 000021020120000 000233323332000 002322222223200 011322232223110 020322030223020 022223323322220 020322030223020 011322232223110 002322222223200 000233323332000 000021020120000 000001222100000 000000000000000
Nim
Our program uses Rust algorithm (and also its colors 🙂) and the formula to compute grid size from number of particles comes from Pascal algorithm. Number of particles is an input from user. The program displays the values on the terminal if there are not too many and produce a PNG image. Code to produce a PPM image is also provided but not used. <lang Nim>
- Abelian sandpile.
from math import sqrt from nimPNG import savePNG24 from sequtils import repeat from strformat import fmt from strutils import strip, addSep, parseInt
- The grid represented as a sequence of sequences of int32.
type Grid = seq[seq[int32]]
- Colors to use for PPM and PNG files.
const Colors = [[byte 100, 40, 15],
[byte 117, 87, 30],
[byte 181, 134, 47],
[byte 245, 182, 66]]
- ---------------------------------------------------------------------------------------------------
func sideLength(initVal: int32): int32 =
# Return the grid side length needed for "initVal" particles. # We make sure that the returned value is odd. result = sqrt(initVal.toFloat / 1.75).int32 + 3 result += result and 1 xor 1
- ---------------------------------------------------------------------------------------------------
func doOneStep(grid: var Grid; boundary: var array[4, int]): bool =
## Compute one step.
result = false
for y in boundary[0]..boundary[2]:
for x in boundary[1]..boundary[3]:
if grid[y][x] >= 4:
let rem = grid[y][x] div 4
grid[y][x] = grid[y][x] mod 4
if y - 1 >= 0:
inc grid[y - 1][x], rem
if y == boundary[0]:
dec boundary[0]
if x - 1 >= 0:
inc grid[y][x - 1], rem
if x == boundary[1]:
dec boundary[1]
if y + 1 < grid.len:
inc grid[y + 1][x], rem
if y == boundary[2]:
inc boundary[2]
if x + 1 < grid.len:
inc grid[y][x + 1], rem
if x == boundary[3]:
inc boundary[3]
result = true
- ---------------------------------------------------------------------------------------------------
proc display(grid: Grid; initVal: int) =
## Display the grid as an array of values.
echo fmt"Starting with {initVal} particles."
echo ""
var line = newStringOfCap(2 * grid.len - 1)
for row in grid:
for value in row:
line.addSep(" ", 0)
line.add($value)
echo line
line.setLen(0)
echo ""
- ---------------------------------------------------------------------------------------------------
proc writePpmFile(grid: Grid; name: string) =
## Write a grid representation in a PPM file.
var file = open(name, fmWrite)
file.write(fmt"P6 {grid.len} {grid.len} 255 ")
for row in grid:
for value in row:
discard file.writeBytes(Colors[value], 0, 3)
file.close()
echo fmt"PPM image written in ""{name}""."
- ---------------------------------------------------------------------------------------------------
proc writePngFile(grid: Grid; name: string) =
## Write a grid representation in a PNG file.
var pixels = newSeq[byte](3 * grid.len * grid.len)
# Build pixel list.
var idx = 0
for row in grid:
for value in row:
pixels[idx..idx+2] = Colors[value]
inc idx, 3
discard savePNG24(name, pixels, grid.len, grid.len)
echo fmt"PNG image written in ""{name}""."
- ---------------------------------------------------------------------------------------------------
proc askInitVal(): int32 =
# Ask user for the number of particles.
while true:
stdout.write("Number of particles? ")
try:
let input = stdin.readLine().strip().parseInt()
if input in 4..int32.high:
return input.int32
echo fmt"Value not in expected range: 4..{int32.high}"
except ValueError:
echo "Invalid input"
except EOFError:
quit(QuitSuccess)
- ---------------------------------------------------------------------------------------------------
- Initialize the grid.
let initVal = askInitVal() let sideLen = sideLength(initVal) var grid = repeat(newSeq[int32](sideLen), sideLen) let origin = grid.len div 2 var boundaries: array[4, int] = [origin, origin, origin, origin] grid[origin][origin] = initVal
- Run the simulation.
while doOneStep(grid, boundaries):
discard
- Display grid.
if grid.len <= 20:
grid.display(initVal)
- grid.writePpmFile(fmt"grid_{initVal}.ppm")
grid.writePngFile(fmt"grid_{initVal}.png") </lang>
- Output:
Number of particles? 100 Starting with 100 particles. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 3 2 0 2 3 0 0 0 0 0 3 0 3 2 3 0 3 0 0 0 1 2 3 0 3 0 3 2 1 0 0 2 0 2 3 0 3 2 0 2 0 0 1 2 3 0 3 0 3 2 1 0 0 0 3 0 3 2 3 0 3 0 0 0 0 0 3 2 0 2 3 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 PNG image written in "grid_100.png".
Pascal
The main optimization was to spread the sand immediatly.mul := val DIV 4;//not only := val -4so that only (sand mod 4) stays in place.runtime for abelian(1e6) down to 1min 20 secs from 9 min
Memorizing the used colums of the rows has little effect when choosing the right size of the grid.Only 11 secs for abelian(1e6) -> 1min 9sec
Python shows 64 too.
<lang pascal>
program Abelian2;
{$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$CODEALIGN proc=16}{$ALIGN 16}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
SysUtils;
type
Tlimit = record
lmtLow,LmtHigh : LongWord;
end;
TRowlimits = array of Tlimit;
tOneRow = pLongWord;
tGrid = array of LongWord;
var
Grid: tGrid; Rowlimits:TRowlimits; s : AnsiString; maxval,maxCoor : NativeUint;
function CalcMaxCoor(maxVal : NativeUint):NativeUint; // maxVal = 10000;maxCoor = 77-2;// maxCoor*maxCoor *1,778; 0.009sec // maxVal = 100000;maxCoor = 236-2;// maxCoor*maxCoor *1.826; 0.825sec // maxVal = 1000000;maxCoor = 732-2;// maxCoor*maxCoor *1.877; 74 sec Begin
result := trunc(sqrt(maxval/1.75))+3;
end;
procedure clear; begin
setlength(Grid,0); setlength(Rowlimits,0); s := ;
end;
procedure InitGrid(var G:tGrid;InitVal:NativeUint); var
row,middle: nativeINt;
begin // setlength(Rowlimits,0); setlength(G,0);
MaxCoor := CalcMaxCoor(InitVal); setlength(G,sqr(maxCoor)); setlength(Rowlimits,maxCoor); fillchar(G[0],length(G)*SizeOf(G[0]),#0);
middle := (maxCoor) div 2;
Grid[middle*maxcoor+middle] := InitVal;
For row := 1 to maxCoor do
with Rowlimits[row] do
Begin
lmtLow := middle;
lmtHigh := middle;
end;
with Rowlimits[middle] do Begin lmtLow := middle; lmtHigh := middle; end;
end; procedure OutGridPPM(const G:tGrid;maxValue : NativeUint); const
color : array[0..3] of array[0..2] of Byte =
//R,G,B)
((0,0,0),
(255,0,0),
(0,255,0),
(0,0,255));
var
f :text; pActRow: tOneRow; col,row,sIdx,value : NativeInt;
Begin
Assignfile(f,'ppm/Grid_'+IntToStr(maxValue)+'.ppm');
rewrite(f);
write(f,Format('P6 %d %d %d ',[maxCoor-1,maxCoor-1,255]));
setlength(s,(maxCoor-1)*3);
pActRow :=@G[0];
For row := maxCoor-2 downto 0 do
Begin
inc(pActRow,maxCoor);
sIdx := 1;
For col := 1 to maxCoor-1 do
Begin
value := pActRow[col];
s[sIdx] := CHR(color[value,0]);
s[sIdx+1] := CHR(color[value,1]);
s[sIdx+2] := CHR(color[value,2]);
inc(sIdx,3);
end;
write(f,s);
end;
CloseFile(f);
end;
procedure OutGrid(const G:tGrid); //output of grid and test, if no sand is lost var
pActRow: tOneRow; col,row,sum,value : NativeUint;
Begin
setlength(s,maxcoor-1);
pActRow := @G[0];
sum := 0;
For row := maxCoor-1 downto 1 do
Begin
inc(pActRow,maxcoor);
For col := 1 to maxCoor-1 do
Begin
value := pActRow[col];
// IF value>=4 then writeln(row:5,col:5,value:13);
s[col] := chr(value+48);
inc(sum,value);
end;
if maxCoor <80 then
writeln(s);
end;
writeln('columns ',maxcoor-1,' checksum ',maxVal,' ?=? ',sum);
{
For row := 1 to maxCoor do
with Rowlimits[row] do
writeln(lmtLow:10,lmtHigh:10);
* }
end;
procedure Evolution(var G:tGrid); var
pActRow,pRowBefore,pRowAfter : tOneRow; col,row,mul,val,done : NativeUint;
begin
repeat
pRowBefore := @G[0];
pActRow := @G[maxcoor];
pRowAfter := @G[2*maxcoor];
done := 0;
For row := maxCoor-1 downto 1 do
Begin
with RowLimits[row] do
Begin
while (LmtLow >1) AND (pActRow[lmtLow]<> 0) do
dec(lmtLow);
while (lmtHigh < maxCoor) AND (pActRow[lmtHigh]<> 0) do
inc(lmtHigh);
For col := lmtLow to lmtHigh do
Begin
val := pActRow[col];
IF val >=4 then
Begin
mul := val DIV 4;
done := val;
inc(pRowBefore[col],mul);
inc(pActRow[col-1],mul);
pActRow[col] := val-4*Mul;
inc(pActRow[col+1],mul);
inc(pRowAfter[col],mul);
end;
end;
pRowBefore:= pActRow;
pActRow := pRowAfter;
inc(pRowAfter,maxcoor);
end;
end;
until done=0;
end;
procedure OneTurn(count:NativeUint); begin
Writeln(' Test abelian sandpile( ',count,' )');
MaxVal := count;
InitGrid(Grid,count);
Evolution(Grid);
OutGrid(Grid);
OutGridPPM(Grid,count);
clear;
end;
BEGIN
OneTurn(4); OneTurn(16); OneTurn(64); OneTurn(1000); OneTurn(10000); OneTurn(100000);
END. </lang>
- Output:
Test abelian sandpile( 4 ) 010 101 010 columns 3 checksum 4 ?=? 4 Test abelian sandpile( 16 ) 00100 02120 11011 02120 00100 columns 5 checksum 16 ?=? 16 Test abelian sandpile( 64 ) 00121000 02222200 12222210 22202220 12222210 02222200 00121000 00000000 columns 8 checksum 64 ?=? 64 Test abelian sandpile( 1000 ) 0000000001111111000000000 0000000130233320310000000 0000013223313133223100000 0000213222130312223120000 0002220123332333210222000 0011223233123213323221100 0033032313221223132303300 0122123203311133023212210 0322231023333333201322230 1032333332231322333332301 1231312332232322332131321 1313322133322233312233131 1330231131220221311320331 1313322133322233312233131 1231312332232322332131321 1032333332231322333332301 0322231023333333201322230 0122123203311133023212210 0033032313221223132303300 0011223233123213323221100 0002220123332333210222000 0000213222130312223120000 0000013223313133223100000 0000000130233320310000000 0000000001111111000000000 columns 25 checksum 1000 ?=? 1000 Test abelian sandpile( 10000 ) --shortened columns 77 checksum 10000 ?=? 10000 Test abelian sandpile( 100000 ) columns 241 checksum 100000 ?=? 100000 real 0m0,815s
OCaml
In Sandpile module (sandpile.ml) <lang OCaml> module Make =
functor (M : sig val m : int val n : int end) -> struct
let grid = Array.init M.m (fun _ -> Array.make M.n 0)
let print () =
for i = 0 to M.m - 1
do for j = 0 to M.n - 1
do Printf.printf "%d " grid.(i).(j)
done
; print_newline ()
done
let unstable = Hashtbl.create 10
let add_grain x y
= grid.(x).(y) <- grid.(x).(y) + 1
; if grid.(x).(y) >= 4 then
Hashtbl.replace unstable (x,y) () (* Use Hashtbl.replace for uniqueness *)
let topple x y
= grid.(x).(y) <- grid.(x).(y) - 4
; if grid.(x).(y) < 4
then Hashtbl.remove unstable (x,y)
; match (x,y) with
(* corners *)
| (0,0) -> add_grain 1 0
; add_grain 0 1
| (0,n) when n = M.n - 1
-> add_grain 1 n
; add_grain 0 (n-1)
| (m,0) when m = M.m - 1
-> add_grain m 1
; add_grain (m-1) 0
| (m,n) when m = M.m - 1 && n = M.n - 1
-> add_grain ( m ) (n-1)
; add_grain (m-1) ( n )
(* sides *)
| (0,y) -> add_grain 1 y
; add_grain 0 (y+1)
; add_grain 0 (y-1)
| (m,y) when m = M.m - 1
-> add_grain ( m ) (y-1)
; add_grain ( m ) (y+1)
; add_grain (m-1) ( y )
| (x,0) -> add_grain (x+1) 0
; add_grain (x-1) 0
; add_grain ( x ) 1
| (x,n) when n = M.n - 1
-> add_grain (x-1) ( n )
; add_grain (x+1) ( n )
; add_grain ( x ) (n-1)
(* else *)
| (x,y) -> add_grain ( x ) (y+1)
; add_grain ( x ) (y-1)
; add_grain (x+1) ( y )
; add_grain (x-1) ( y )
let add_sand n x y
= for i = 1 to n
do add_grain x y
done
let avalanche ()
= while Hashtbl.length unstable > 0
do
let unstable' = Hashtbl.fold (fun (x,y) () r -> (x,y) :: r) unstable []
in List.iter (fun (x,y) -> topple x y ) unstable'
done
end
(* testing *)
let ()
= let module S = Make (struct let m = 11 let n = 11 end) in S.add_sand 500 5 5 ; S.avalanche () ; S.print ()
</lang>
Perl
<lang Perl>#!/usr/bin/perl
use strict; # http://www.rosettacode.org/wiki/Abelian_sandpile_model use warnings;
my ($high, $wide) = split ' ', qx(stty size); my $mask = "\0" x $wide . ("\0" . "\177" x ($wide - 2) . "\0") x ($high - 5) .
"\0" x $wide;
my $pile = $mask =~ s/\177/ rand() < 0.02 ? chr 64 + rand 20 : "\0" /ger;
for ( 1 .. 1e6 )
{
print "\e[H", $pile =~ tr/\0-\177/ 1-~/r, "\n$_";
my $add = $pile =~ tr/\1-\177/\0\0\0\200/r; # set high bit for >=4
$add =~ /\200/ or last;
$pile =~ tr/\4-\177/\0-\173/; # subtract 4 if >=4
for ("\0$add", "\0" x $wide . $add, substr($add, 1), substr $add, $wide)
{
$pile |= $_;
$pile =~ tr/\200-\377/\1-\176/; # add one to each neighbor of >=4
$pile &= $mask;
}
select undef, undef, undef, 0.1; # comment out for full speed
}</lang>
Phix
Generates moving images similar to the julia output. The distributed version also has variable speed, additional display modes, and a random dropping toggle.
-- demo\rosetta\Abelian_sandpile_model.exw include pGUI.e Ihandle dlg, canvas cdCanvas cddbuffer sequence board = {{0,0,0}, {0,0,0}, {0,0,0}} procedure drop(integer y, x) sequence moves = {} while true do board[y,x] += 1 if board[y,x]>=4 then board[y,x] -= 4 moves &= {{y,x-1},{y,x+1},{y-1,x},{y+1,x}} end if -- extend board if rqd (maintain a border of zeroes) if x=1 then -- extend left for i=1 to length(board) do board[i] = prepend(board[i],0) end for for i=1 to length(moves) do moves[i][2] += 1 end for elsif x=length(board[1]) then -- extend right for i=1 to length(board) do board[i] = append(board[i],0) end for end if -- (copy the all-0 lines from the other end...) if y=1 then -- extend up board = prepend(board,board[$]) for i=1 to length(moves) do moves[i][1] += 1 end for elsif y=length(board) then -- extend down board = append(board,board[1]) end if if length(moves)=0 then exit end if {y,x} = moves[$] moves = moves[1..$-1] end while IupUpdate(canvas) end procedure function timer_cb(Ihandle /*ih*/) integer y = floor(length(board)/2)+1, x = floor(length(board[1])/2)+1 drop(y,x) return IUP_DEFAULT end function function redraw_cb(Ihandle ih, integer /*posx*/, integer /*posy*/) IupGLMakeCurrent(ih) cdCanvasActivate(cddbuffer) cdCanvasClear(cddbuffer) for y=1 to length(board) do for x=1 to length(board[1]) do integer c = board[y][x] if c!=0 then integer colour = {CD_VIOLET,CD_RED,CD_BLUE}[c] cdCanvasPixel(cddbuffer, x, y, colour) end if end for end for cdCanvasFlush(cddbuffer) return IUP_DEFAULT end function function map_cb(Ihandle ih) IupGLMakeCurrent(ih) atom res = IupGetDouble(NULL, "SCREENDPI")/25.4 cddbuffer = cdCreateCanvas(CD_GL, "300x100 %g", {res}) cdCanvasSetBackground(cddbuffer, CD_PARCHMENT) return IUP_DEFAULT end function procedure main() IupOpen() canvas = IupGLCanvas("RASTERSIZE=300x100") IupSetCallbacks({canvas}, {"ACTION", Icallback("redraw_cb"), "MAP_CB", Icallback("map_cb")}) dlg = IupDialog(canvas,"TITLE=\"Abelian sandpile model\"") IupShow(dlg) Ihandle timer = IupTimer(Icallback("timer_cb"), 10) IupMainLoop() IupClose() end procedure main()
Python
Python: Original, with output
<lang Python> import numpy as np import matplotlib.pyplot as plt
def iterate(grid):
changed = False
for ii, arr in enumerate(grid):
for jj, val in enumerate(arr):
if val > 3:
grid[ii, jj] -= 4
if ii > 0:
grid[ii - 1, jj] += 1
if ii < len(grid)-1:
grid[ii + 1, jj] += 1
if jj > 0:
grid[ii, jj - 1] += 1
if jj < len(grid)-1:
grid[ii, jj + 1] += 1
changed = True
return grid, changed
def simulate(grid):
while True:
grid, changed = iterate(grid)
if not changed:
return grid
if __name__ == '__main__':
start_grid = np.zeros((10, 10)) start_grid[4:5, 4:5] = 64 final_grid = simulate(start_grid.copy()) plt.figure() plt.gray() plt.imshow(start_grid) plt.figure() plt.gray() plt.imshow(final_grid)
</lang> Output: </n> Before: <lang Python> [[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. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0.64. 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. 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.]]
</lang> After: <lang Python> [[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 2. 1. 0. 0. 0. 0.] [0. 0. 2. 2. 2. 2. 2. 0. 0. 0.] [0. 1. 2. 2. 2. 2. 2. 1. 0. 0.] [0. 2. 2. 2. 0. 2. 2. 2. 0. 0.] [0. 1. 2. 2. 2. 2. 2. 1. 0. 0.] [0. 0. 2. 2. 2. 2. 2. 0. 0. 0.] [0. 0. 0. 1. 2. 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.]]
</lang>
An interactive variant to the above solution:
<lang python>
from os import system, name
from time import sleep
def clear(): if name == 'nt': _ = system('cls') else: _ = system('clear')
def exit(): import sys clear() sys.exit()
def make_area(x, y): area = [[0]*x for _ in range(y)] return area
def make_sandpile(area, loc, height): loc=list(n-1 for n in loc) x, y = loc
try: area[y][x]+=height except IndexError: pass
def run(area, by_frame=False): def run_frame(): for y_index, group in enumerate(area): y = y_index+1
for x_index, height in enumerate(group): x = x_index+1
if height < 4: continue
else: make_sandpile(area, (x+1, y), 1) make_sandpile(area, (x, y+1), 1)
if x_index-1 >= 0: make_sandpile(area, (x-1, y), 1) if y_index-1 >= 0: make_sandpile(area, (x, y-1), 1)
make_sandpile(area, (x, y), -4)
while any([any([pile>=4 for pile in group]) for group in area]): if by_frame: clear() run_frame() if by_frame: show_area(area); sleep(.05)
def show_area(area): display = [' '.join([str(item) if item else ' ' for item in group]) for group in area] [print(i) for i in display]
clear() if __name__ == '__main__': area = make_area(10, 10) print('\nBefore:') show_area(area) make_sandpile(area, (5, 5), 64) run(area) print('\nAfter:') show_area(area) </lang>
Output: <lang> Before: 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 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0
After: 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 1 2 2 2 2 2 1 0 0 0 2 2 2 0 2 2 2 0 0 0 1 2 2 2 2 2 1 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 1 2 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 </lang>
Python: using tkinter
<lang python> Python 3.6.5 code using Tkinter graphical user
interface (Canvas widget) to display final results.
from tkinter import *
- given a grid, display it on a tkinter Canvas:
class Sandpile:
def __init__(self, wn, grid):
self.window = wn
self.grid = grid
self.canvas = Canvas(wn, bg='lemon chiffon')
self.canvas.pack(fill=BOTH, expand=1)
colors = {0:'dodger blue',
1:'red',
2:'green',
3:'lemon chiffon'}
x = 10
y = 10
d = 5
for row in self.grid:
for value in row:
clr = colors[value]
self.canvas.create_rectangle(
x, y, x+d, y+d,
outline=clr,
fill = clr)
x += 5
x = 10
y += 5
class Grid:
def __init__(self, size, center):
self.size = size # rows/cols in (best if odd)
self.center = center # start value at center of grid
self.grid = [[0]*self.size for i in range(self.size)]
self.grid[self.size // 2][self.size // 2] = self.center
# print the grid:
def show(self, msg):
print(' ' + msg + ':')
for row in self.grid:
print(' '.join(str(x) for x in row))
print()
return
# dissipate piles of sand as required:
def abelian(self):
while True:
found = False
for r in range(self.size):
for c in range(self.size):
if self.grid[r][c] > 3:
self.distribute(self.grid[r][c], r, c)
found = True
if not found:
return
# distribute sand from a single pile to its neighbors:
def distribute(self, nbr, row, col):
qty, remain = divmod(nbr, 4)
self.grid[row][col] = remain
for r, c in [(row+1, col),
(row-1, col),
(row, col+1),
(row, col-1)]:
self.grid[r][c] += qty
return
# display the grid using tkinter:
def display(self):
root = Tk()
root.title('Sandpile')
root.geometry('700x700+100+50')
sp = Sandpile(root, self.grid)
root.mainloop()
- execute program for size, center value pair:
- just print results for a small grid
g = Grid(9,17) g.show('BEFORE') g.abelian() # scatter the sand g.show('AFTER')
- just show results in tkinter for a large grid
- I wish there was a way to attach a screen shot
- of the tkinter result here
g = Grid(131,25000) g.abelian() # scatter the sand g.display() # display result using tkinter
- OUTPUT:
- BEFORE:
- 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 0 0 0 0 0 0 0
- 0 0 0 0 17 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 0 0
- 0 0 0 0 0 0 0 0 0
- AFTER:
- 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
- 0 0 0 2 1 2 0 0 0
- 0 0 1 1 1 1 1 0 0
- 0 0 0 2 1 2 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
</lang>
Raku
(formerly Perl 6)
Terminal based
Defaults to a stack of 1000 and showing progress. Pass in a custom stack size if desired and -hide-progress to run without displaying progress (much faster.)
<lang perl6>sub cleanup { print "\e[0m\e[?25h\n"; exit(0) }
signal(SIGINT).tap: { cleanup(); exit(0) }
unit sub MAIN ($stack = 1000, :$hide-progress = False );
my @color = "\e[38;2;0;0;0m█",
"\e[38;2;255;0;0m█",
"\e[38;2;255;255;0m█",
"\e[38;2;0;0;255m█",
"\e[38;2;255;255;255m█"
;
my ($h, $w) = qx/stty size/.words».Int; my $buf = $w * $h; my @buffer = 0 xx $buf; my $done;
@buffer[$w * ($h div 2) + ($w div 2) - 1] = $stack;
print "\e[?25l\e[48;5;232m";
repeat {
$done = True;
loop (my int $row; $row < $h; $row = $row + 1) {
my int $rs = $row * $w; # row start
my int $re = $rs + $w; # row end
loop (my int $idx = $rs; $idx < $re; $idx = $idx + 1) {
if @buffer[$idx] >= 4 {
my $grains = @buffer[$idx] div 4;
@buffer[ $idx - $w ] += $grains if $row > 0;
@buffer[ $idx - 1 ] += $grains if $idx - 1 >= $rs;
@buffer[ $idx + $w ] += $grains if $row < $h - 1;
@buffer[ $idx + 1 ] += $grains if $idx + 1 < $re;
@buffer[ $idx ] %= 4;
$done = False;
}
}
}
unless $hide-progress {
print "\e[1;1H", @buffer.map( { @color[$_ min 4] }).join;
}
} until $done;
print "\e[1;1H", @buffer.map( { @color[$_ min 4] }).join;
cleanup;</lang>
Passing in 2048 as a stack size results in: Abelian-sandpile-model-perl6.png (offsite .png image)
SDL2 Animation
<lang perl6>use NativeCall; use SDL2::Raw;
my ($width, $height) = 900, 900;
unit sub MAIN ($stack = 10000);
my int ($w, $h) = 160, 160;
my $buf = $w * $h; my @buffer = 0 xx $buf;
@buffer[$w * ($h div 2) + ($w div 2) - 1] = $stack;
SDL_Init(VIDEO);
my SDL_Window $window = SDL_CreateWindow(
"Abelian sandpile - Raku", SDL_WINDOWPOS_CENTERED_MASK, SDL_WINDOWPOS_CENTERED_MASK, $width, $height, RESIZABLE
);
my SDL_Renderer $renderer = SDL_CreateRenderer( $window, -1, ACCELERATED +| TARGETTEXTURE );
my $asp_texture = SDL_CreateTexture($renderer, %PIXELFORMAT<RGB332>, STREAMING, $w, $h);
my $pixdatabuf = CArray[int64].new(0, $w, $h, $w);
my @color = 0x00, 0xDE, 0x14, 0xAA, 0xFF;
sub render {
my int $pitch; my int $cursor;
# work-around to pass the pointer-pointer. my $pixdata = nativecast(Pointer[int64], $pixdatabuf); SDL_LockTexture($asp_texture, SDL_Rect, $pixdata, $pitch);
$pixdata = nativecast(CArray[int8], Pointer.new($pixdatabuf[0]));
loop (my int $row; $row < $h; $row = $row + 1) {
my int $rs = $row * $w; # row start
my int $re = $rs + $w; # row end
loop (my int $idx = $rs; $idx < $re; $idx = $idx + 1) {
$pixdata[$idx] = @buffer[$idx] < 4 ?? @color[@buffer[$idx]] !! @color[4];
if @buffer[$idx] >= 4 {
my $grains = floor @buffer[$idx] / 4;
@buffer[ $idx - $w ] += $grains if $row > 0;
@buffer[ $idx - 1 ] += $grains if $idx - 1 >= $rs;
@buffer[ $idx + $w ] += $grains if $row < $h - 1;
@buffer[ $idx + 1 ] += $grains if $idx + 1 < $re;
@buffer[ $idx ] %= 4;
}
}
}
SDL_UnlockTexture($asp_texture);
SDL_RenderCopy($renderer, $asp_texture, SDL_Rect, SDL_Rect.new(:x(0), :y(0), :w($width), :h($height))); SDL_RenderPresent($renderer);
}
my $event = SDL_Event.new;
main: loop {
while SDL_PollEvent($event) {
my $casted_event = SDL_CastEvent($event);
given $casted_event {
when *.type == QUIT {
last main;
}
}
}
render(); print fps;
}
say ;
sub fps {
state $fps-frames = 0;
state $fps-now = now;
state $fps = ;
$fps-frames++;
if now - $fps-now >= 1 {
$fps = [~] "\b" x 40, ' ' x 20, "\b" x 20 ,
sprintf "FPS: %5.2f ", ($fps-frames / (now - $fps-now)).round(.01);
$fps-frames = 0;
$fps-now = now;
}
$fps
}</lang>
Passing in a stack size of 20000 results in: Abelian-sandpile-sdl2.png (offsite .png image)
Rust
<lang rust>// This is the main algorithm. // // It loops over the current state of the sandpile and updates it on-the-fly. fn advance(field: &mut Vec<Vec<usize>>, boundary: &mut [usize; 4]) -> bool {
// This variable is used to check whether we changed anything in the array. If no, the loop terminates. let mut done = false;
for y in boundary[0]..boundary[2]
{
for x in boundary[1]..boundary[3]
{
if field[y][x] >= 4
{
// This part was heavily inspired by the Pascal version. We subtract 4 as many times as we can
// and distribute it to the neighbors. Also, in case we have outgrown the current boundary, we
// update it to once again contain the entire sandpile.
// The amount that gets added to the neighbors is the amount here divided by four and (implicitly) floored.
// The remaining sand is just current modulo 4.
let rem: usize = field[y][x] / 4;
field[y][x] %= 4;
// The isize casts are necessary because usize can not go below 0.
// Also, the reason why x and y are compared to boundary[2]-1 and boundary[3]-1 is because for loops in
// Rust are upper bound exclusive. This means a loop like 0..5 will only go over 0,1,2,3 and 4.
if y as isize - 1 >= 0 {field[y-1][x] += rem; if y == boundary[0] {boundary[0]-=1;}}
if x as isize - 1 >= 0 {field[y][x-1] += rem; if x == boundary[1] {boundary[1]-=1;}}
if y+1 < field.len() {field[y+1][x] += rem; if x == boundary[2]-1 {boundary[2]+=1;}}
if x+1 < field.len() {field[y][x+1] += rem; if y == boundary[3]-1 {boundary[3]+=1;}}
done = true;
}
}
}
done
}
// This function can be used to display the sandpile in the console window. // // Each row is mapped onto chars and those characters are then collected into a string. // These are then printed to the console. // // Eg.: [0,1,1,2,3,0] -> [' ','░','░','▒','▓',' ']-> " ░░▒▓ " fn display(field: &Vec<Vec<usize>>) {
for row in field
{
let char_row = {
row.iter().map(|c| {match c {
0 => ' ',
1 => '░',
2 => '▒',
3 => '▓',
_ => '█'
}}).collect::<String>()
};
println!("{}", char_row);
}
}
// This function writes the end result to a file called "output.ppm". // // PPM is a very simple image format, however, it entirely uncompressed which leads to huge image sizes. // Even so, for demonstrative purposes it's perfectly good to use. For something more robust, look into PNG libraries. // // Read more about the format here: http://netpbm.sourceforge.net/doc/ppm.html fn write_pile(pile: &Vec<Vec<usize>>) {
use std::fs::File; use std::io::Write;
// We first create the file (or erase its contents if it already existed).
let mut file = File::create("./output.ppm").unwrap();
// Then we add the image signature, which is "P3 <newline>[width of image] [height of image]<newline>[maximum value of color]<newline>".
write!(file, "P3\n{} {}\n255\n", pile.len(), pile.len()).unwrap();
for row in pile {
// For each row, we create a new string which has more or less enough capacity to hold the entire row.
// This is for performance purposes, but shouldn't really matter much.
let mut line = String::with_capacity(row.len() * 14);
// We map each value in the field to a color.
// These are just simple RGB values, 0 being the background, the rest being the "sand" itself.
for elem in row {
line.push_str(match elem {
0 => "100 40 15 ",
1 => "117 87 30 ",
2 => "181 134 47 ",
3 => "245 182 66 ",
_ => unreachable!(),
});
}
// Finally we write this string into the file.
write!(file, "{}\n", line).unwrap();
}
}
fn main() {
// This is how big the final image will be. Currently the end result would be a 16x16 picture. let field_size = 16; let mut playfield = vec![vec![0; field_size]; field_size];
// We put the initial sand in the exact middle of the field. // This isn't necessary per se, but it ensures that sand can fully topple. // // The boundary is initially just the single tile which has the sand in it, however, as the algorithm // progresses, this will grow larger too. let mut boundary = [field_size/2-1, field_size/2-1, field_size/2, field_size/2]; playfield[field_size/2 - 1][field_size/2 - 1] = 16;
// This is the main loop. We update the field until it returns false, signalling that the pile reached its
// final state.
while advance(&mut playfield, &mut boundary) {};
// Once this happens, we simply display the result. Uncomment the line below to write it to a file. // Calling display with large field sizes is not recommended as it can easily become too large for the console. display(&playfield); //write_pile(&playfield);
}</lang>
Output:
░
▒░▒
░░ ░░
▒░▒
░
VBA
<lang VBA>Sub SetupPile(a As Integer, b As Integer) Application.ScreenUpdating = False For i = 1 To a For j = 1 To b Cells(i, j).value = "" Cells(i, j).Select
With Selection.Borders(xlEdgeLeft)
.LineStyle = xlContinuous .Weight = xlMedium
End With With Selection.Borders(xlEdgeTop)
.LineStyle = xlContinuous .Weight = xlMedium
End With With Selection.Borders(xlEdgeBottom)
.LineStyle = xlContinuous .Weight = xlMedium
End With With Selection.Borders(xlEdgeRight)
.LineStyle = xlContinuous .Weight = xlMedium
End With
With Selection
.HorizontalAlignment = xlCenter .VerticalAlignment = xlCenter
End With
Next j Next i Application.ScreenUpdating = True End Sub
Sub Abelian_Sandpile()
Dim PileWidth As Integer
Dim PileHeight As Integer
Dim FieldArray() As Integer
Debug.Print "Start:" & Now()
'Set Size of Playing Field PileWidth = 25 PileHeight = 25
ReDim FieldArray(PileWidth - 1, PileHeight - 1)
'Paint Basic Grid SetupPile PileWidth, PileHeight
'Drop sand amount into middle of playing field SandDropAmount = 1000 'Get around excel's incorrect rounding SandDropColumn = Round((PileWidth / 2) + 0.001, 0) SandDropRow = Round((PileHeight / 2) + 0.001, 0)
Cells(SandDropRow, SandDropColumn) = SandDropAmount FieldArray(SandDropRow - 1, SandDropColumn - 1) = SandDropAmount
Continue = False
'Check if Pile is already stabilized at the start For i = 1 To PileWidth 'Col For j = 1 To PileHeight 'Row If FieldArray(j - 1, i - 1) > 3 Then Continue = True Next j Next i
'While not stabilized While Continue For i = 1 To PileWidth For j = 1 To PileHeight
If FieldArray(j - 1, i - 1) > 3 Then 'Reduce by 4 FieldArray(j - 1, i - 1) = FieldArray(j - 1, i - 1) - 4 'Increase Neighbours 't If j >= 2 Then FieldArray(j - 2, i - 1) = FieldArray(j - 2, i - 1) + 1 'r If i < PileWidth Then FieldArray(j - 1, i) = FieldArray(j - 1, i) + 1 'b If j < PileHeight Then FieldArray(j, i - 1) = FieldArray(j, i - 1) + 1 'l If i >= 2 Then FieldArray(j - 1, i - 2) = FieldArray(j - 1, i - 2) + 1 'Next round GoTo Nextone End If
Next j Next i
Nextone:
'Check if now stabilized Continue = False For i = 1 To PileWidth For j = 1 To PileHeight 'Paint every step if needed 'Cells(j, i) = FieldArray(j - 1, i - 1)
If FieldArray(j - 1, i - 1) > 3 Then Continue = True Next j Next i
Wend
'Print out final step For i = 1 To PileWidth For j = 1 To PileHeight Cells(j, i) = FieldArray(j - 1, i - 1) Next j Next i
'Make field square and remove 0 Cells.Select Selection.ColumnWidth = 2 Selection.RowHeight = 13.5 Selection.Replace What:="0", Replacement:="", LookAt:=xlPart, SearchOrder:=xlByRows, MatchCase:=False, SearchFormat:=False, ReplaceFormat:=False Range("A1").Select
Range(Cells(1, 1), Cells(PileHeight, PileWidth)).Select
'Conditional Format Selection.FormatConditions.AddColorScale ColorScaleType:=3 Selection.FormatConditions(Selection.FormatConditions.Count).SetFirstPriority Selection.FormatConditions(1).ColorScaleCriteria(1).Type = xlConditionValueLowestValue With Selection.FormatConditions(1).ColorScaleCriteria(1).FormatColor
.Color = 8109667 .TintAndShade = 0
End With Selection.FormatConditions(1).ColorScaleCriteria(2).Type = xlConditionValuePercentile Selection.FormatConditions(1).ColorScaleCriteria(2).value = 50 With Selection.FormatConditions(1).ColorScaleCriteria(2).FormatColor
.Color = 8711167 .TintAndShade = 0
End With Selection.FormatConditions(1).ColorScaleCriteria(3).Type = xlConditionValueHighestValue With Selection.FormatConditions(1).ColorScaleCriteria(3).FormatColor
.Color = 7039480 .TintAndShade = 0
End With Range("A1").Select
Debug.Print "W,H,A:" & PileWidth & "," & PileHeight & "," & SandDropAmount Debug.Print "End:" & Now()
End Sub</lang> Output:
On Excel Page
Wren
<lang ecmascript>import "/fmt" for Fmt
class Sandpile {
// 'a' is a list of integers in row order
construct new(a) {
var count = a.count
_rows = count.sqrt.floor
if (_rows * _rows != count) Fiber.abort("The matrix of values must be square.")
_a = a
_neighbors = List.filled(count, 0)
for (i in 0...count) {
_neighbors[i] = []
if (i % _rows > 0) _neighbors[i].add(i-1)
if ((i + 1)%_rows > 0) _neighbors[i].add(i+1)
if (i - _rows >= 0) _neighbors[i].add(i-_rows)
if (i + _rows < count) _neighbors[i].add(i+_rows)
}
}
isStable { _a.all { |i| i <= 3 } }
// topples until stable
topple() {
while (!isStable) {
for (i in 0..._a.count) {
if (_a[i] > 3) {
_a[i] = _a[i] - 4
for (j in _neighbors[i]) _a[j] = _a[j] + 1
}
}
}
}
toString {
var s = ""
for (i in 0..._rows) {
for (j in 0..._rows) s = s + Fmt.swrite("$2d ", _a[_rows*i + j])
s = s + "\n"
}
return s
}
}
var printAcross = Fn.new { |str1, str2|
var r1 = str1.split("\n")
var r2 = str2.split("\n")
var rows = r1.count - 1
var cr = (rows/2).floor
for (i in 0...rows) {
var symbol = (i == cr) ? "->" : " "
Fmt.print("$s $s $s", r1[i], symbol, r2[i])
}
System.print()
}
var a1 = List.filled(25, 0) a1[12] = 4 var a2 = List.filled(25, 0) a2[12] = 6 var a3 = List.filled(25, 0) a3[12] = 16 var a4 = List.filled(100, 0) a4[55] = 64 for (a in [a1, a2, a3, a4]) {
var s = Sandpile.new(a) var str1 = s.toString s.topple() var str2 = s.toString printAcross.call(str1, str2)
}</lang>
- Output:
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 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 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 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 64 0 0 0 0 -> 0 0 2 2 2 0 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0