diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..496ee2c --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.DS_Store \ No newline at end of file diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..0a8d966 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,7 @@ +Copyright (c) 2017 Florian Roscheck + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..0faf991 --- /dev/null +++ b/README.md @@ -0,0 +1,172 @@ +Visual Bak-Tang-Wiesenfeld Sandpile Model for Matlab +==================================================== +This Matlab library helps to visualize the avalanche and power law characteristics of the popular sandpile model. + +The library is meant for visualizing the Bak-Tang-Wiesenfeld sandpile model and plotting the sandpile's statistics +on the fly. The visualization has been optimized for performance. + +General Description of the Sandpile Model +----------------------------------------- +The Bak-Tang-Wiesenfeld sandpile model (also called *Abelian Sandpile Model*) represents a sandpile to which sand +grains are added one by one. It is described in an intuitive way in Per Bak's famous book +[How nature works: the science of self-organized criticality](https://books.google.com/books/about/How_nature_works.html?id=e5XuAAAAMAAJ). +The model gives an easily understandable example of the concept of self-organized criticality and shows how large-scale, +complex behavior can emerge from local, simple patterns in nature. + +As observed in the real world, parts of a sandpile might start "sliding" off its slope when enough grains have +accumulated at a location. Once a sand grain has started its sliding movement, it might also bump into other sand grains +and cause them to slide. This event is known as "avalanche". In his research, Per Bak has observed that the number of +grains that take part in an avalanche (as a measure of the avalanche size) follows a power-law distribution. + +A power-law distribution means that some quantity *D* can be expressed as some power of another quantity *s*, for +example as *D(s) = k*s^-τ*, where *τ* and *k* are constants. When plotted on a log-log plot, a power-law distribution +shows through the plotted data resembling a straight line. + +![Log-log plot showing avalanche size distribution](sandpile_log-log_plot.png "Avalanche Sizes Follow Power Law") + +The plot shown above has been created from data gathered by this sandpile model. Up to an avalanche size of 100 grains, +the data supports the claim of a power-law distribution: The data resembles a straight line, approximated by the graph +of a linear fit function shown in green. + +The power-law characteristics of the sandpile system are very important for Bak's line of thought. He remarks: + +> The power law indicates that the stationary state is critical. +> We conclude that he pile has self-organized into a critical state. +> -- Bak, P. *How nature works : the science of self-organized criticality*. New York: Springer, 1996, 57. + +In this quote, the terms *critical* and *self-organized* jump out. The *critical* state, in the context of the sandpile +model, relates to how individual sand grains behave *as a group*. + +A pile of sand that consists of many randomly arranged sand grains might not have a critical state. In this kind of +sandpile, any grain that is added, at any position, will probably not have global effects on the sandpile. A newly added +grain might cause a few grains to topple in its neighborhood, but there is no remarkable pattern of topples or +avalanches this grain causes. + +However, a sandpile that is in a critical state is different. In such a sandpile, a grain that is added at a random +position will eventually cause global effects in the sandpile. The added grain might start an avalanche of toppling +processes that change the entire geometry of the sandpile, reaching far beyond the local domain. For a sandpile that +has a typical conical form, the critical state is reached when its slope has reached a certain angle. Keeping all +parameters the same and adding more sand grains, sand sliding off its slope will prevent the sandpile's slope angle +from becoming any steeper. + +To put it in Bak's words: In a non-critical state, individual grains follow their own local dynamics, while in a +critical state, the emergent dynamics are global. In the critical state, a pattern emerges from the observation of the +sandpile: Avalanche sizes (and other measurable quantities) show a power-law distribution. + +The sandpile does organize itself without external influences into the critical state. It is entirely *self-organized*. +This means that even if grains are initially distributed randomly and not in a specific way triggering a critical state, +the sandpile will show critical behavior after adding a certain number of sand grains to arbitrary positions. + +How This Sandpile Model Works +----------------------------- +In this sandpile model, the sandpile is built onto a square grid of a specified `pile_width` The total number of +positions in the square grid is `pile width * pile width`, so for a sandpile with a `pile_width` of 12, 144 grid +positions are available. + +On every grid position, an arbitrary amount of sand grains can be stacked. However, if the number of sand grains on a +grid position is greater than or equal to 4, a toppling event takes place, originating at that grid position. + +In a toppling event, four grains get taken away from the originating grid position and are distributed equally among the +neighboring grid positions north, east, south, and west of the originating grid position. As a result, the originating +grid position possesses four grains less, while the number of grains on all its four neighboring positions has increased +by one. This behavior eventually leads to a chain reaction, in which a further toppling event on one of the neighboring +positions is caused. Such a chain reaction is called "avalanche". If, through a toppling event, a sand grain is to be +added to a position beyond the boundary of the pile, it is ignored and no longer considered in the model. This can be +imagined as a sand grain "falling off" the surface on which the sandpile is created. + +In this model, new sand grains are added to random grid positions. The sand grains are added one by one, except during +toppling or avalanche events. In this case, no grains are added until these events are over. The model is initialized +with a random number of up to 3 grains at every grid position. (As noted [before](#general-description-of-the-sandpile-model), +this random initialization is not necessarily critical, the critical behavior evolves over time.) + +Matlab Plots for Visualization +------------------------------ +During execution of this script, Matlab visualizes the sandpile by means of a sandpile plot and a plot showing the +avalanche size distribution. The sandpile plot can be turned off to increase performance. + +### Sandpile Plot +A sample of the sandpile plot is shown in the picture below. The plot shows a "top view" of the entire pile. The +pile has a width of 22 positions, so 22*22=484 positions in total. A square is shown for every grid position. The +lighter the square in color, the more grains are stacked on the grid position. On the color scale, there is no +discrimination between positions that have four and more grains, all of them are white. A black square indicates that +there are no grains placed on that position. + +The peach-colored square overlaying a grid position shows the location where the most recent sand grain is dropped. The +square is not shown during an avalanche. + +The sandpile plot is updated automatically as the model is running. + +![Matlab sandpile plot](sandpile_grid.png "Matlab plot showing sandpile") + +### Plot of Avalanche Sizes +The plot of avalanche sizes is a log-log plot that is created and updated as the model runs. It shows the quantity of +avalanches and their sizes. In this model, the size of an avalanche is defined as the number of direct and indirect +toppling events that were caused by the addition of a single sand grain to the pile. + +Through the log-log scale the plot enables an assessment of whether the avalanche distribution follows a power law. As +there are note many observations for the biggest avalanches, the plot usually shows some clutter for the biggest +avalanches. + +![Matlab plot of avalanche sizes](sandpile_size_plot.png "Matlab plot showing avalanche sizes") + +Usage +----- +The code is thoroughly documented, the documentation can be accessed through the `doc simulateSandpile` command. To run +the model, call the `simulateSandpile` function, the only user-facing function in the script. + +Configurable **input** parameters of `simulateSandpile` are: + +* `pile_width`: Side length of the square pile + +* `no_of_grains`: Number of grains that are added one by one to the sandpile as the model runs. The simulation stops +after the specified number of grains has been added. + +* `draw_speed`: Speed of animation of the plots. A speed of 0 skips the sandpile animation entirely, and will only plot +a chart of avalanche sizes once the simulation has finished running. A value of 0.5 yields a relatively slow animation, smaller +values yield a faster one. + +The following code is a simple usage example: +```matlab +pile_width = 40; % Pile represented through 40x40 grid +no_of_grains = 4000; % Run the simulation until 4000 new sand grains have been added +draw_speed = 0.2 % Animate the pile relatively slow +avalanche_output = simulateSandpile(pile_width, no_of_grains, draw_speed) +``` + +The `simulateSandpile` function returns a single **output** variable, `avalanche_output`. This variable is a matrix with +the shape *max. avalanche size x 2*. The matrix contains information about how many avalanches with a specific size +occurred during the run. ([The avalanche plot description](#plot-of-avalanche-sizes) shows how the avalanche size is +defined in this model.) + +Here is an annotated example of a typical `avalanche_output`: +```matlab +% Content of avalanche_output: +[1 50; % 50 avalanches causing 1 toppling event + 2 32; % 32 avalanches causing 2 toppling events + 3 11; % 11 avalanche causing 3 toppling events + ...]; +``` + +Download +-------- +An archive with the library can be downloaded from the [releases page](https://github.com/flrs/visual_sandpile/releases). + +Installation +------------ +To install the library, extract all files into a folder and add them to your Matlab path. + +Contribution +------------ +I am happy about any contribution or feedback. Please let me know about your comments via the Issues tab on +[GitHub](https://github.com/flrs/visual_sandpile/issues). + +Acknowledgements +---------------- +The main source of inspiration for writing this library is Per Bak's book +[How nature works: the science of self-organized criticality](https://books.google.com/books/about/How_nature_works.html?id=e5XuAAAAMAAJ). + +The function `rldCumsumDiff` that repeats copies of array elements is taken from [a post byDivakar's on stackoverflow.com](http://stackoverflow.com/a/29079288/2778484). + +License +------- +The Visual Bak-Tang-Wiesenfeld Sandpile Model for Matlab is distributed under the [MIT License (MIT)](https://github.com/flrs/visual_sandpile/blob/master/LICENSE.md). \ No newline at end of file diff --git a/file_exchange_thumb.png b/file_exchange_thumb.png new file mode 100644 index 0000000..89fed04 Binary files /dev/null and b/file_exchange_thumb.png differ diff --git a/sandpile/plotPile.m b/sandpile/plotPile.m new file mode 100644 index 0000000..75994de --- /dev/null +++ b/sandpile/plotPile.m @@ -0,0 +1,90 @@ +function plotPile(pile_store, pile_store_add,pile_img, pointer_patch,... + no_grains, avalanche_plt, avalanche_ct_plot, avalanche_desc_text,... + draw_speed) +%plotPile - Plot the sandpile and a chart counting avalanche sizes +%The function plots all events happening in a single avalanche. The +%function has been optimized for performance. +% +% Syntax: plotPile(pile_store, pile_store_add, no_grains, avalanche_plt,... +% pointer_patch, pile_img, avalanche_ct_plot, avalanche_desc_text,... +% draw_speed) +% +% Inputs: +% pile_store - Pile history, matrix of shape (pile width, pile width, +% no. of history time steps) +% pile_store_add - Vector indicating locations of added grains in +% history +% pile_img - Handle to the image containing the sandpile +% pointer_patch - Handle of the patch that shows where new grains have +% been dropped +% no_grains - Total count of grains added to sandpile so far +% avalanche_plt - Vector containing counts of avalanche sizes +% avalanche_ct_plot - Handle of avalanche counter plot +% avalanche_desc_text - Handle of descriptive text on avalanche counter +% plot +% draw_speed - Speed of animation +% +% Outputs: +% none +% +% Other m-files required: none +% Subfunctions: none +% MAT-files required: none +% +% See also: setupPlots +% +% Author: Florian Roscheck +% Website: http://github.com/flrs/visual_sandpile +% January 2017; Last revision: 27-January-2017 + +%------------- BEGIN CODE -------------- + +% plot pointer for new grain and sandpile +if draw_speed>0 + % pointer for new grain + pile_width = size(pile_store(:, :, 1), 1); + + % fast ind2sub (see http://tipstrickshowtos.blogspot.com/2011/09/fast-r + % eplacement-for-ind2sub.html, checked on 2017-01-26) + add_row = rem(pile_store_add-1, pile_width)+1; + add_col = (pile_store_add-add_row)/pile_width + 1; + + set(pointer_patch, 'XData',... + [add_col-0.25 add_col-0.25 add_col+0.25 add_col+0.25],... + 'YData', [add_row+0.25 add_row-0.25 add_row-0.25 add_row+0.25]); + + drawnow + + pause(2*draw_speed) % pause for visualization + + set(pointer_patch, 'XData', 0, 'YData', 0); % hide new grain pointer + + % sandpile + ct_goal = size(pile_store, 3); % show all avalanche time steps + + for ct = 1:ct_goal + set(pile_img, 'cdata', pile_store(:, :, ct)+1); + + drawnow + + pause(draw_speed) % pause for visualization + end +end + +% plot avalanche chart +set(avalanche_ct_plot, 'XData', 1:numel(avalanche_plt),... + 'YData', avalanche_plt); + +% update text in avalanche chart +figure(get(get(avalanche_ct_plot, 'Parent'), 'Parent')); + +cur_xlim=xlim; +cur_ylim=ylim; +cur_string = get(avalanche_desc_text, 'String'); +cur_string{1} = [num2str(no_grains) ' sand grains']; + +set(avalanche_desc_text, 'position', [cur_xlim(2) cur_ylim(2) 0]); +set(avalanche_desc_text, 'String', cur_string); + +drawnow; +%------------- END CODE -------------- diff --git a/sandpile/resolvePeaks.m b/sandpile/resolvePeaks.m new file mode 100644 index 0000000..8f2f5a0 --- /dev/null +++ b/sandpile/resolvePeaks.m @@ -0,0 +1,93 @@ +function [pile, intermediate_piles] = resolvePeaks(pile, peak_pos) +%resolvePeaks - Resolve all peaks in a pile +% +% Syntax: [pile, intermediate_piles] = resolvePeaks(pile, peak_pos) +% +% Inputs: +% pile - Matrix of shape (pile width, pile width, +% no. of history time steps), with integer values from 0 to 4 +% peak_pos - Vector containing positions of all peaks +% +% Outputs: +% pile - Matrix of shape (pile width, pile width), with integer values +% from 0 to 4, with peaks in initial pile resolved (might now contain +% peaks resulting from resolving the initial peaks) +% intermediate_piles - Matrix of shape (pile width, pile width, no. of +% intermediate time steps), with integer values from 0 to 4, containing +% all intermediate steps taken in resolving the peaks +% +% Example: +% [pile, intermediate_piles] = resolvePeaks([4 1;3 2], 1) +% +% Other m-files required: none +% Subfunctions: none +% MAT-files required: none +% +% See also: scanPileForPeaks +% +% Author: Florian Roscheck +% Website: http://github.com/flrs/visual_sandpile +% January 2017; Last revision: 27-January-2017 + +%------------- BEGIN CODE -------------- +%% initialize +pile_width = size(pile,1); +prealloc_size = round(pile_width^1.3); % preallocate empty array depending + % on pile side length, it is + % expected that no. of piles + % increases exponentially with pile + % side length, exponent 1.3 is + % arbitrary and a tradeoff between + % unneccessarily slow intialization + % and unneccessary overhead when + % expanding the array later in the + % code +intermediate_piles = zeros(pile_width,pile_width,prealloc_size); +intermediate_pile_ct = 1; + +peak_pattern = [0 1 0;1 -4 1;0 1 0;]; % pattern for resolving peaks, + % characteristic of Abelian sandpile + +pile_frame = zeros(pile_width+2); % construct frame around pile to catch + % falling off the grid +pile_frame(2:end-1,2:end-1) = pile; % insert pile into frame + +%% process peaks +% resolve peaks one by one +for peak = 1:numel(peak_pos) + % fast ind2sub (see http://tipstrickshowtos.blogspot.com/2011/09/fast-r + % eplacement-for-ind2sub.html, checked on 2017-01-26) + peakY = rem(peak_pos(peak)-1, pile_width)+1; + peakX = (peak_pos(peak)-peakY)/pile_width + 1; + + % resolve peaks + pile_frame(peakY:peakY+2, peakX:peakX+2) = ... + pile_frame(peakY:peakY+2, peakX:peakX+2)+peak_pattern; + + % extract new pile from frame + pile = pile_frame(2:end-1, 2:end-1); + + % expand intermediate pile array when it has reached its size limit + if intermediate_pile_ct>size(intermediate_piles, 3) + intermediate_piles = ... + cat(3, intermediate_piles, ... + zeros(pile_width, pile_width, prealloc_size)); + end + + % append new pile to intermediate pile array + intermediate_piles(:, :, intermediate_pile_ct) = pile; + + intermediate_pile_ct = intermediate_pile_ct+1; +end + +if intermediate_pile_ct>1 + % eliminate unused, preallocated entries from intermediate pile array + intermediate_piles = ... + intermediate_piles(:, :, 1:intermediate_pile_ct-1); +else + % no piles resolved + intermediate_piles = []; +end + +end + diff --git a/sandpile/rldCumsumDiff.m b/sandpile/rldCumsumDiff.m new file mode 100644 index 0000000..b2dc482 --- /dev/null +++ b/sandpile/rldCumsumDiff.m @@ -0,0 +1,8 @@ +function out = rldCumsumDiff(vals,runlens) +% Divakar's solution to run-length decoding +% (http://stackoverflow.com/a/29079288/2778484) +clens = cumsum(runlens); +idx(clens(end))=0; +idx([1 clens(1:end-1)+1]) = diff([0 vals]); +out = cumsum(idx); +return \ No newline at end of file diff --git a/sandpile/scanPileForPeaks.m b/sandpile/scanPileForPeaks.m new file mode 100644 index 0000000..b709c85 --- /dev/null +++ b/sandpile/scanPileForPeaks.m @@ -0,0 +1,29 @@ +function peak_pos = scanPileForPeaks(pile) +%scanPileForPeaks - Find all peaks in a sandpile +% +% Syntax: peak_pos = scanPileForPeaks(pile) +% +% Inputs: +% pile - Matrix of shape (pile width, pile width) with integer values +% from 0 to 4 +% +% Outputs: +% peak_pos - Vector containing positions of all peaks +% +% Example: +% peak_pos = scanPileForPeaks([0 1 3; 4 1 2; 2 1 0]); +% +% Other m-files required: none +% Subfunctions: none +% MAT-files required: none +% +% See also: simulateSandpile +% +% Author: Florian Roscheck +% Website: http://github.com/flrs/visual_sandpile +% January 2017; Last revision: 27-January-2017 + +%------------- BEGIN CODE -------------- +peak_pos = find(pile>=4); +%------------- END CODE -------------- + diff --git a/sandpile/setupPlots.m b/sandpile/setupPlots.m new file mode 100644 index 0000000..280b703 --- /dev/null +++ b/sandpile/setupPlots.m @@ -0,0 +1,77 @@ +function [pointer_patch, pile_img, avalanche_ct_plot,... + avalanche_desc_text] = setupPlots(pile_width, draw_speed) +%setupPlots - This function sets up the sandpile and the avalanche plot +% +% Syntax: [pointer_patch, pile_img avalanche_ct_plot,... +% avalanche_desc_text] = setupPlots(pile_width, draw_speed) +% +% Inputs: +% pile_width - Side length of the square pile +% draw_speed - Speed of animation +% +% Outputs: +% pointer_patch - Handle of the patch that shows where new grains have +% been dropped +% pile_img - Handle to the image containing the sandpile +% avalanche_ct_plot - Handle of avalanche counter plot +% avalanche_desc_text - Handle of descriptive text on avalanche counter +% plot +% +% Example: +% [pointer_patch, pile_img avalanche_ct_plot, avalanche_desc_text] = ... +% setupPlots(20, 0.25) +% +% Other m-files required: none +% Subfunctions: none +% MAT-files required: none +% +% See also: plotPile +% +% Author: Florian Roscheck +% Website: http://github.com/flrs/visual_sandpile +% January 2017; Last revision: 27-January-2017 + +%------------- BEGIN CODE -------------- +%% Set up avalanche size plot +figure('position', [750 200 700 500], 'Color', [1 1 1]); +avalanche_ct_plot = loglog(0, 0, '.-k', 'LineWidth', 1.5, 'MarkerSize', 10); + +title('Avalanche Sizes Follow Power Law'); +xlabel('Avalanche size D(s)'); +ylabel('No. of observed avalanches s'); + +grid on +set(gca, 'TickDir', 'out') +box off + +% set up descriptive text +avalanche_desc_text = text(1, 1,... + {'0 sand grains'; + [num2str(pile_width) 'x' num2str(pile_width) ' pile size']},... + 'HorizontalAlignment','right',... + 'BackgroundColor','w'); + +%% Set up sandpile plot +if draw_speed + figure('position', [200 200 500 500]); + set(gcf, 'Units', 'normal') + set(gca, 'Position', [0 0 1 1]) + set(gcf, 'Color', [1 1 1]) + + pile_img = image(zeros(pile_width)); + + % set image-specific properties + colormap(gray(5)); + xlim([0.5 pile_width+0.5]); + ylim([0.5 pile_width+0.5]); + set(gca, 'xtick', [],'ytick', []); + + % initialize patch for new grains + pointer_patch = patch([0 0 0 0], [0 0 0 0], 0,... + 'EdgeColor', 'none', 'FaceColor', [244 165 130]/255); +else + % draw_speed == 0, so we are not even plotting the sandpile + pile_img = []; + pointer_patch = []; +end +%------------- END CODE -------------- diff --git a/sandpile/simulateSandpile.m b/sandpile/simulateSandpile.m new file mode 100644 index 0000000..e368605 --- /dev/null +++ b/sandpile/simulateSandpile.m @@ -0,0 +1,135 @@ +function avalanche_output = simulateSandpile(pile_width, no_of_grains, ... + draw_speed) +%simulateSandpile - Simulate a square, planar Abelian sandpile +%The Abelian sandpile is also known as the Bak-Tang-Wiesenfeld model. +% +%In this function, it is assumed that there is a square grid with the side +%length `pile_width`. When initiating the model, up to three sand grains +%are added randomly to every position in the grid. Then, a total of +%`no_of_grains` sand grains are added to the grid sequentially at random +%positions. If a sand grain causes the number of grains present at a grid +%position to exceed 3, an avalanche is started from the grid position. In +%an avalanche, four grains are removed from the starting grid position. One +%grain of these four grains is added to the grid position toward the top, +%right, bottom and left, respectively. If this causes any of these grid +%positions to exceed three grains, an avalanche is started from the +%respective positions. The avalanche process continues until all positions +%in the grid carry not more than three grains. Then, the next sand grain is +%added at a random position, as described above. +%If a pile is resolved that is at the boundary of the grid, some grains +%will fall "off the grid" and no longer be considered in the model. +% +% Syntax: avalanche_output = simulateSandpile(pile_width, no_of_grains, ... +% draw_speed) +% +% Inputs: +% pile_width - Side length of the square pile +% no_of_grains - No. of grains that should be added to the sandpile +% draw_speed - Speed of animation. A speed of 0 skips the animation +% entirely, and will only plot a chart of avalanche sizes once the +% simulation has finished running. A value of 0.5 yields a relatively +% slow animation, smaller values yield a faster one. +% +% Outputs: +% avalanche_output - Matrix containing count of observed avalanche +% sizes. The second column shows the no. of avalanches +% observed, the first column equals the no. of observations in such +% configuration. For example. a vector like [1 8; 2 3; 3 1] would +% mean that 8 avalanches with only 1 grid point exceeding 3 grains +% have be observed, 3 avalanches with 2 grid points exceeding 3 +% grains, and 1 avalanche with 3 grid points exceeding 3 grains. +% +% Example: +% avalanche_output = simulateSandpile(15, 5000, 0.25) +% +% Other m-files required: setupPlots.m, scanPileForPeaks.m, resolvePeaks.m, +% plotPile.m +% Subfunctions: none +% MAT-files required: none +% +% See also: none +% +% Author: Florian Roscheck +% Website: http://github.com/flrs/visual_sandpile +% January 2017; Last revision: 27-January-2017 + +%------------- BEGIN CODE -------------- +%% check inputs +assert(nargin >= 3, 'The function %s needs at least 3 inputs.', mfilename); + +check_vars = {pile_width, 'pile_width'; + no_of_grains, 'no_of_grains'; + draw_speed, 'draw_speed'}; + +for ct = 1:size(check_vars,1) + assert(isfloat(check_vars{ct,1}), ['The input variable "%s" needs to'... + ' be a floating point number.'], check_vars{ct,2}); + assert(numel(check_vars{ct,1})==1, ['The input variable "%s" needs to'... + ' be a scalar.'], check_vars{ct,2}); +end + +for ct = 1:2 + assert(check_vars{ct,1}>0, ['The value of the input variable "%s" '... + 'needs to be greater than zero.'], check_vars{ct,2}); +end + +%% initialize +% initialize variables +pile = zeros(pile_width); +pile_store = zeros([pile_width, pile_width, 1]); +pile_store_add = 0; +avalanche_plt = 0; + +% initialize pile +genRand = @(inp) round(rand()*3); +pile = arrayfun(genRand, pile); + +% initialize plots +[pointer_patch, pile_img, avalanche_ct_plot, avalanche_desc_text] = ... + setupPlots(pile_width, draw_speed); + +%% run model +% generate pile +for ct = 1:no_of_grains + % add grain to pile + fprintf('Adding grain %.0f of %.0f...\n', ct, no_of_grains); + add_pos = round(1+rand()*((pile_width^2)-1)); + pile(add_pos) = pile(add_pos) + 1; + pile_store = pile; + pile_store_add = add_pos; + avalanche_size = 0; + + % resolve peaks + peaks = scanPileForPeaks(pile); + intermediate_piles = []; + while numel(peaks) ~= 0 + [pile, intermediate_piles] = resolvePeaks(pile, peaks); + if ~isempty(intermediate_piles) + if draw_speed + pile_store = cat(3, pile_store, intermediate_piles); + end + avalanche_size = avalanche_size + size(intermediate_piles, 3); + end + peaks = scanPileForPeaks(pile); + end + + % update avalanche counter + if avalanche_size > 0 + if numel(avalanche_plt)>=avalanche_size + avalanche_plt(avalanche_size) = ... + avalanche_plt(avalanche_size)+1; + else + avalanche_plt(avalanche_size) = 1; + end + fprintf(['Captured avalanche with duration of %.0f time ',... + 'steps.\n'], avalanche_size); + end + + % call plot function + plotPile(pile_store, pile_store_add, pile_img, pointer_patch,... + ct, avalanche_plt, avalanche_ct_plot,... + avalanche_desc_text, draw_speed); +end + +avalanche_output = [1:numel(avalanche_plt);avalanche_plt]'; +%------------- END CODE -------------- diff --git a/sandpile_grid.png b/sandpile_grid.png new file mode 100644 index 0000000..5cff0bf Binary files /dev/null and b/sandpile_grid.png differ diff --git a/sandpile_log-log_plot.png b/sandpile_log-log_plot.png new file mode 100644 index 0000000..de0d5e4 Binary files /dev/null and b/sandpile_log-log_plot.png differ diff --git a/sandpile_size_plot.png b/sandpile_size_plot.png new file mode 100644 index 0000000..654913d Binary files /dev/null and b/sandpile_size_plot.png differ