Abelian Sandpile Model
We call this Mathematica function to retrieve some background information from Wiki page: WikipediaData["Abelian sandpile model"]
The Abelian sandpile model, also known as the Bak[Dash]Tang[Dash]Wiesenfeld model, was the first discovered example of a dynamical system displaying self-organized criticality. It was introduced by Per Bak, Chao Tang and Kurt Wiesenfeld in a 1987 paper.The model is a cellular automaton.
Also in NKS there is discussion about self-organized criticality and some simple demo codes in Mathematica. You may also find this Numberhile video on this topic interesting, typically it talks about where "Abelian" comes from based on group theory point of view.
I found a GitHub repo with a nice implementation of the model in C++ and it works quite well for large number ( ~ 1 million initial sand particles at origin). I think this can be a good example to show how to use LTemplate
package to call this C++ package from Mathematica if you do not want to start from scratch using the basic LibraryLink C API's.
Implementation
Follow this link to install LTemplate
package and refer to the document page coming with this package to find the basic example and some type conversion information (well, not like Mathematica you need to specify the type by yourself). In particular we need mma::IntMatrixRef
as return type and to transfer a 2D array back to Mathematica. Meanwhile, we also need ra = mma::makeMatrix<mint>
(ra
stands for raw array) to store the grid information during computation. This ra
also replaces the double pointer in the original code to handle the array for plot later. Here is the code
#include <LTemplate.h>
#include <numeric>
#include <iomanip>
class sandpile {
public:
const int dx[4] = { 1,0,-1,0 };
const int dy[4] = { 0,1,0,-1 };
const int N_size = 781;
const int N_size2 = 390;
struct L_coord
{
unsigned int x;
unsigned int y;
};
mma::IntMatrixRef CreateSandpile(int n ) {
auto ra = mma::makeMatrix<mint>(N_size,N_size);
std::fill(ra.begin(),ra.end(),0);
bool** V_sites = new bool*[N_size];
bool** Move_ = new bool*[N_size];
unsigned int** Odometer = new unsigned int*[N_size];
for (int k = 0; k < N_size; k++)
{
V_sites[k] = new bool[N_size];
Odometer[k] = new unsigned int[N_size];
Move_[k] = new bool[N_size];
for (int i = 0; i<N_size; i++)
{
Move_[k][i] = false; V_sites[k][i] = false;
}
}
L_coord* walking_ = new L_coord[N_size*N_size];
V_sites[N_size2][N_size2] = true;
Move_[N_size2][N_size2] = true;
ra(N_size2,N_size2) = n;
walking_[0].x = N_size2;
walking_[0].y = N_size2;
unsigned int x = 0, y = 0;
int top_ = 0;
unsigned int lx = 0, ly = 0;
unsigned int d = 0;
unsigned int max_of_top = 0;
unsigned long N_of_Moves = 0;
while (top_ >= 0)
{
N_of_Moves += 1;
if (max_of_top < top_) { max_of_top = top_; }
x = walking_[top_].x; y = walking_[top_].y;
top_--; Move_[x][y] = false;
d = (ra(x,y) >> 2);
ra(x,y) = ra(x,y) - (d << 2);
for (int k = 0; k < 4; k++)
{
lx = x + dx[k];
ly = y + dy[k];
V_sites[lx][ly] = true;
ra(lx,ly) += d;
if (Move_[lx][ly] == false && ra(lx,ly) >= 4)
{
walking_[++top_].x = lx;
walking_[top_].y = ly;
Move_[lx][ly] = true;
}
}
}
for (int k = 0; k<N_size; k++)
{
delete[] V_sites[k]; delete[] Odometer[k]; delete[] Move_[k];
}
delete[] V_sites; delete[] Odometer; delete[] Move_; delete[] walking_;
return ra;
}
};
I keep most part intact, only made several essential modification according to the document of LTemplate:
- Include
LTemplate
head file
- Create a class called
sandpile
which wraps around the CreateSandpile
function to be imported
- Return
mma::IntMatrixRef
object instead of Clock
object used to print computation time cost
Compile and Load
In a Mathematica notebook, I run the following lines. First I need to load the package
<< LTemplate`
SetDirectory[NotebookDirectory[]] (*put sandpile.h in the same directory as the working nb*)
we need to tell Mathematica which class to load and to pick the right member function at the same time
template =LClass["sandpile",{
LFun["CreateSandpile", {Integer}, {Integer, 2}]
}];
If our compiler is set up correctly, with verbose log, we should compile the class to dll just like in VS developer command or xcode or some clang/gcc environment on Linux.
CompileTemplate[template,
"ShellCommandFunction" -> Print, "ShellOutputFunction"->Print
]
The last line is the dll in the $UserBaseDirectory
. Then we can load the dll with
LoadTemplate[template]
Import the class
sp = Make["sandpile"]
and call the member function:
sp@"CreateSandpile"[100]
I omitted one ArrayTrim function in C++ code to simplify the process. We actually need it to trim margins with all zeros in a row or column. I choose to implement it quickly in Mathematica with much shorter lines than the c++ version, even though we might transferred many null points back to kernel wasted. Nevertheless,
cleanRow[arr_] := DeleteCases[arr, _?(Total@# == 0 &)]
For small set, the
With[{arr = cleanRow@Transpose@cleanRow@(sp@"CreateSandpile"[5000])},
ArrayPlot[arr, ColorRules -> {0 -> Black, 1 -> Magenta, 2 -> Red, 3 -> Blue}]]
Growing complexity
Style[Grid[Partition[Map[With[{arr=cleanRow@Transpose@cleanRow@(sp@"CreateSandpile"[#])},
ArrayPlot[arr,ColorRules->{0->Black,1->Magenta,2-> Red,3-> Blue}]]&,
{10,20,50,80,100,200,500,800,1000,5000,10000,50000,100000,200000,500000}],5]],
ImageSizeMultipliers->0.4]
I tested the code on 10 million initial points and no problem at all to display it nicely with ArrayPlot
. The animation at the beginning is only shows a quarter of the whole graph due to symmetry. The color scheme is the same as the original post.
1 Million particles at center at t-zero
With[{arr = cleanRow@Transpose@cleanRow@(sp2@"CreateSandpile"[1000000])},
ArrayPlot[arr,ColorRules -> {0 -> Black, 1 -> Magenta, 2 -> Red, 3 -> Blue}]]
10 Million particles at center at t-zero (only show quarter) and the result is attached in csv format.
ImageCrop@
Colorize[sp1@"CreateSandpile"[10000000],
ColorRules -> {0 -> White, 1 -> Green, 2 -> Purple,
3 -> RGBColor["#D4AF37"]}]
Summary
- Use
LibraryLink
based package LTemplate
to call a C++ function with one integer input and 2D integer output.
- Demonstrate how to pass 2D array in C++ back to Mathematica
- Use
ArrayPlot
to visualize tensor computed in C++ environment
- Safely wrap
Manipulate
or ListAnimate
around external library function
- Mathematica implementation can also be found here
- Analytic research on this topic's complexity
Attachments: