# Self-Organized Criticality Simulation (SocSim)

Faculty of Physics. University of Warsaw

Documentation Status Build Status codecov

Project is created as part of subject: Team student projects Faculty of Physics

Programs in Python that simulate dynamical systems that have a critical point as an attractor. So called self-organized criticality(SOC)

### 1. Theoretical problem description

Self-organized criticality wiki:

> In physics, self-organized criticality (SOC) is a property of dynamical systems that have a critical point as an attractor. Their macroscopic behavior thus displays the spatial or temporal scale-invariance characteristic of the critical point of a phase transition, but without the need to tune control parameters to a precise value, because the system, effectively, tunes itself as it evolves towards criticality.

> The concept was put forward by Per Bak, Chao Tang and Kurt Wiesenfeld (“BTW”) in a paper published in 1987 in Physical Review Letters, and is considered to be one of the mechanisms by which complexity arises in nature. Its concepts have been enthusiastically applied across fields as diverse as geophysics, physical cosmology, evolutionary biology and ecology, bio-inspired computing and optimization (mathematics), economics, quantum gravity, sociology, solar physics, plasma physics, neurobiology and others.

> SOC is typically observed in slowly driven non-equilibrium systems with a large number of degrees of freedom and strongly nonlinear dynamics. Many individual examples have been identified since BTW’s original paper, but to date there is no known set of general characteristics that guarantee a system will display SOC.

### 2. Program structure, installation and use cases

#### 2.1 Project folder structure

Project folder structure is inspired by these sources: sources1 source2 and Kwant project.

##### socsim:

- docsrc - holds Sphinx scripts used for documentation generation.
- docs - GitHub configuration folder, which holds web-page of project.
- resource - Non executable files.
- results - folder used for holding results of simulation, Jupiter notebooks and different use cases.
- SOC - main project folder, which holds all source code.
  - models - contains different SOC models, like: Abelian sandpile model, forest-fire model, etc..
  - common - common code between all models
  - tests - unit tests of code

#### 2.2 Installation and dependencies

##### Dependencies

Mostly numerical libraries, visualsation, web page generation etc. For whole list take a look at requirements.txt

#### 2.3 Use cases

##### Running program

Program is designed in next way:

- Framework part - placed under `SOC` folder.
- Research part - consists of jupyter notebooks(which can be easily deployed to web-page) and is placed under `research folder`

##### Developing the program

use `python setup.py develop` to install a basic set of dependencies and link the package to be importable in your current Python environment.

##### Running test cases

To make folder SOC an import package, run only once:

    python setup.py develop

After that, simply use `pytest SOC` to automatically find and execute all existing test cases.

##### Web-page generation

Web page is generated using Sphinx library.

Under terminal enter into `docsrc` folder and type:

    make html

Web-page will be generated into `./docsrc/build/html` folder. If you want to update web-page, copy generated web-page into `/docs` folder.

### Results

Click the model’s name for more examples:

#### BTW

power-law

#### Manna model(Abelian/non-Abelian)

mannaexample

#### OFC

ofcexample

#### Forest fire

### Summary of intial goals:

- Make wide coverage of all self-organized criticality models.
- Figuring out the best algorithms.
  - Easy scaling on many processors machine(multithreading, GPGPU, CUDA, numba).
- Using some best practice of programming:
  - Coding conventions
  - Unit Tests.
  - Creation of common modules.
  - Automatic documentation generation.
  - Readability of code and easy of use(between clarity and speed, we should choose clarity).

#### 0.2 Commitments

Here are described code formatting style and other conventions, to make code more uniform. Also this section is for newcomers and contributors.

##### Unit Tests

SocSim uses the lovely PyTest for its unit testing needs. Tests are run automatically on every commit using TravisCI.

##### Documentation

Most popular documentation generator for Python - Sphinx. Good tutorial about using Sphinx here. Here is example of good Google style docstring standardized by PEP-484.

    pip install -U sphinx
    pip install sphinx_rtd_theme
    pip install nbsphinx

pandoc

Dependencies of sphinx: `recommonmark`.

### What’s next?

- How we can apply Keras? Predictions, finding hidden parameters, etc.
- More tests for the batch processes running (Dask)
- Convenient selection of different boundary conditions of a system (lattice)
- Parametrization of the earthquake model for the coverage of wider selection of submodels (by including eg.: drive with random loading, toppling to the neighbours in a specific state (crack model), delay of the fracture initiation, threshold for the fracture propagation) like in Lomnitz-Adler (1993)
- Database of the simulations

### 3. Links/References

- Bak, P., Tang, C. and Wiesenfeld, K. (1987). “Self-organized criticality: an explanation of 1/f noise”. Physical Review Letters. 59 (4): 381–384. Bibcode:1987PhRvL..59..381B. doi:10.1103/PhysRevLett.59.381. PMID 10035754. Papercore summary: http://papercore.org/Bak1987.
- Abelian sandpile model
- Forest-fire model
- Theoretical Models of Self-Organized Criticality (SOC) Systems
- Pink noise
- Introduction to Self-Organized Criticality & Earthquakes
- 25 Years of Self-Organized Criticality: Solar and Astrophysics
- SOC computer simulations
  - Studies in self-organized criticality
- Theoretical Models of SOC Systems

## Manna model

The Manna model is similar in concept to the BTW model. However, where BTW dissipates its “sand grains” deterministically, the Manna model introduces some randomness. Let’s take a look:

    [3]:

    from SOC import Manna
    model = Manna(L=3, save_every=1)
    model.critical_value

    [3]:

    1

This means that the model begins toppling its “sandpiles” once we put two grains somewhere. Let’s try that:

    [20]:

    model = Manna(L=3, save_every=1)
    model.values[2,2] = 2
    model.plot_state(with_boundaries=True);
    model.AvalancheLoop()
    model.plot_state(with_boundaries=True);

Why these two target locations in particular? It actually is random! Let’s rerun that:

    [21]:

    model = Manna(L=3, save_every=1)
    model.values[2,2] = 2
    model.plot_state(with_boundaries=True);
    model.AvalancheLoop()
    model.plot_state(with_boundaries=True);

    [22]:

    model = Manna(L=3, save_every=1)
    model.values[2,2] = 2
    model.plot_state(with_boundaries=True);
    model.AvalancheLoop()
    model.plot_state(with_boundaries=True);

Oh, that’s a bit weird, isn’t it? These seem to have moved awfully far. The trick is that the two grains that fall from the toppling location pick their location at random independently, and here they both picked (1, 1) at first.

Let’s run it for some more time:

    [25]:

    model = Manna(L=5, save_every=1)
    model.run(1000)
    model.animate_states(notebook=True)

    Waiting for wait_for_n_iters=10 iterations before collecting data. This should let the system thermalize.

Your browser does not support the video tag.

Let’s run a larger simulation instead:

    [37]:

    model = Manna(L=10, save_every=1)
    model.run(1000)
    model.animate_states(notebook=True)

    Waiting for wait_for_n_iters=10 iterations before collecting data. This should let the system thermalize.

Your browser does not support the video tag.

If you look closely, you’ll see that the system begins to exhibit very large avalanches very soon:

    [38]:

    model.data_df.AvalancheSize.plot()

    [38]:

    <matplotlib.axes._subplots.AxesSubplot at 0x7f2828870050>

Let’s take a look at how modifying the critical value affects the simulation. We’ll do some more iterations, so the system has the opportunity to “fill up” better. We’ll also skip some animation frames.

    [47]:

    model = Manna(L=10, critical_value=4, save_every=10)
    model.run(5000, wait_for_n_iters = 1000)
    model.animate_states(notebook=True)

    Waiting for wait_for_n_iters=1000 iterations before collecting data. This should let the system thermalize.

Your browser does not support the video tag.

    [48]:

    model.data_df.AvalancheSize.plot()

    [48]:

    <matplotlib.axes._subplots.AxesSubplot at 0x7f281fc7fb90>

Note how the avalanche size grows until a certain period, and then starts to fluctuate randomly at pretty large values. We can try to investigate the histogram of those avalanche sizes. We’ll also fit a line to the linear segment (picked purely subjectively, visually and heuristically).

    [53]:

    model.get_exponent(low=5, high=20)

    y = 3117.481 exp(-1.5522 x)

    [53]:

    {'exponent': -1.552219845297347, 'intercept': 3.4938038036304393}

One thing for sure, there is a region where the scaling in log-log scale is linear. The line fits rather well.

Let’s run a larger simulation and try to estimate the scaling exponent. We’ll wait for a good while so that the system can thermalize well:

    [61]:

    model = Manna(L=40, save_every=100)
    model.run(100000, wait_for_n_iters=50000)

    Waiting for wait_for_n_iters=50000 iterations before collecting data. This should let the system thermalize.

    [62]:

    model.animate_states(notebook=True)

Your browser does not support the video tag.

    [63]:

    model.get_exponent(low = 10, high=100)

    y = 39083.385 exp(-1.3356 x)

    [63]:

    {'exponent': -1.3356064470598588, 'intercept': 4.591992174363189}

Let’s see how reproducible this is:

    [65]:

    model = Manna(L=40, save_every=10000)
    model.run(100000, wait_for_n_iters=50000)
    model.get_exponent(low = 10, high=100)

    Waiting for wait_for_n_iters=50000 iterations before collecting data. This should let the system thermalize.

    y = 38482.710 exp(-1.3329 x)

    [65]:

    {'exponent': -1.3328829112845815, 'intercept': 4.585265642606212}

    [66]:

    model = Manna(L=40, save_every=10000)
    model.run(100000, wait_for_n_iters=50000)
    model.get_exponent(low = 10, high=100)

    Waiting for wait_for_n_iters=50000 iterations before collecting data. This should let the system thermalize.

    y = 36309.665 exp(-1.3169 x)

    [66]:

    {'exponent': -1.316892681945754, 'intercept': 4.5600222427865065}

    [67]:

    model = Manna(L=40, save_every=10000)
    model.run(100000, wait_for_n_iters=50000)
    model.get_exponent(low = 10, high=100)

    Waiting for wait_for_n_iters=50000 iterations before collecting data. This should let the system thermalize.

    y = 35625.884 exp(-1.3077 x)

    [67]:

    {'exponent': -1.3077133473189786, 'intercept': 4.551765652488615}

And the exponent for this model was reported by others (see Table 6.1 in Pruessner) to be around 1.25-1.30.

Which is pretty darn close!

There is another quantity we could calculate here, the number of iterations for an avalanche to finish:

    [68]:

    model.get_exponent("number_of_iterations", low = 10, high=100)

    y = 55075.511 exp(-1.4428 x)

    [68]:

    {'exponent': -1.4428414522987587, 'intercept': 4.7409585345785565}

## Olami–Feder–Christensen (OFC) Earthquake Model

Olami Z., Feder H., Christensen K., Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes, Phys. Rev. Lett. 68, 1992, https://doi.org/10.1103/PhysRevLett.68.1244

modified as in eqns (1) from: Grassberger P., 1994. Efficient large-scale simulations of a uniformly driven system, Phys. Rev. E, 49, 2436–2444, https://doi.org/10.1103/PhysRevE.49.2436

    [2]:

    from SOC.models import OFC

    [3]:

    sim0 = OFC(conservation_lvl=0.2, L=30, save_every = 1)
    sim0.run(1000, wait_for_n_iters=1000)

    Waiting for wait_for_n_iters=1000 iterations before collecting data. This should let the system thermalize.

    [4]:

    sim0.animate_states(notebook=True)

Your browser does not support the video tag.

    [17]:

    sim1 = OFC(L=30)
    sim1.run(100000, wait_for_n_iters = 10000)

    Waiting for wait_for_n_iters=10000 iterations before collecting data. This should let the system thermalize.

    [22]:

    sim1.get_exponent(low=20, high=300)

    y = 42477.635 exp(-1.2325 x)

    [22]:

    {'exponent': -1.232457941524679, 'intercept': 4.628160328647776}

    [23]:

    sim1.get_exponent(column='NumberOfReleases', low=2, high=300)

    y = 28415.887 exp(-1.2413 x)

    [23]:

    {'exponent': -1.2412823768665933, 'intercept': 4.453561223266139}

    [8]:

    sim2 = OFC(conservation_lvl=0.1, L=30)
    sim2.run(500000, wait_for_n_iters = 10000)

    Waiting for wait_for_n_iters=10000 iterations before collecting data. This should let the system thermalize.

    [9]:

    sim2.get_exponent(low=10, high=40)

    y = 14940969.613 exp(-3.4342 x)

    [9]:

    {'exponent': -3.4341725493293325, 'intercept': 7.1743787824846805}

    [16]:

    sim2.get_exponent(column='NumberOfReleases', low=2, high=30)

    y = 251979.007 exp(-2.7298 x)

    [16]:

    {'exponent': -2.7297970986317637, 'intercept': 5.401364360728629}

    [ ]:

    [ ]:

## The BAK–TANG–WIESENFELD Sandpile

source page 85

### General features

- First published by Bak, Tang, and Wiesenfeld (1987).
- Motivated by avalanching behaviour of a real sandpile.
- In one dimension rules represent downward movement of sand grains.
- Defined in any dimension, exactly solved (trivial) in one.
- Stochastic (bulk) drive, deterministic relaxation.
- Non-Abelian in its original definition.
- Many results actually refer to Dhar’s (1990a) Abelian sandpile, Sec. 4.2.
- Simple scaling behaviour disputed, multiscaling proposed.
- Exponents listed in Table 4.1, p. 92, are for the Abelian BTW Model.

### Rules

- d dimensional (usually) hyper-cubic lattice and q the coordination number (on cubic lattices q = 2d).
- Choose (arbitrary) critical slope z^c = q − 1.
- Each site n ∈ {1,…, L}^d has slope z_n.
- Initialisation: irrelevant, model studied in the stationary state.
- Driving: add a grain at n0 chosen at random and update all uphill nearest neighbours n‘0 of n0: z_n0 →z_n0 + q/2 z_n0 →z_n‘0 − 1.
- Toppling: for each site n with z_n > z^c distribute q grains among its nearest neighbours n’ : z_n →z_n − q ∀n’.nn.n z_n →z_n + 1. In one dimension site n = L relaxes according to z_L → z_L − 1 z_L−1 → z_L−1 + 1.
- Dissipation: grains are lost at open boundaries.
- Parallel update: discrete microscopic time, sites exceeding zc at time t topple at t + 1 (updates in sweeps).
- Separation of time scales: drive only once all sites are stable, i.e. z_n ≤ z^c (quiescence).
- Key observables (see Sec. 1.3): avalanche sizes, the total number of topplings until quiescence; avalanche duration T , the total number of parallel updates until quiescence

    [20]:

    from SOC.models import BTW

### Empty model

    [21]:

    b = BTW(L = 50, save_every = 50)

### Running model

    [22]:

    b.run(200000, wait_for_n_iters = 100)

    Waiting for wait_for_n_iters=100 iterations before collecting data. This should let the system thermalize.

    [23]:

    b.data_df.describe()

    [23]:

| AvalancheSize | NumberOfReleases | number_of_iterations
--- | --- | --- | ---
count | 200000.000000 | 200000.000000 | 200000.000000
mean | 87.584035 | 92.235385 | 13.013605
std | 251.254297 | 344.621080 | 32.636649
min | 0.000000 | 0.000000 | 0.000000
25% | 0.000000 | 0.000000 | 0.000000
50% | 0.000000 | 0.000000 | 0.000000
75% | 27.000000 | 12.000000 | 7.000000
max | 2337.000000 | 10060.000000 | 519.000000

    [24]:

    b.get_exponent(low = 2, high = 40)

    y = 17528.587 exp(-1.0186 x)

    [24]:

    {'exponent': -1.018594541688584, 'intercept': 4.243746902353518}

    [25]:

    b.plot_state();

    [26]:

    %matplotlib inline
    a = BTW(30, save_every = 1)
    a.run(10000)

    Waiting for wait_for_n_iters=10 iterations before collecting data. This should let the system thermalize.

    [27]:

    a.animate_states(notebook = True)

Your browser does not support the video tag.

### Uniform model load

    [28]:

    c = BTW(100)
    c.values[1:-1,1:-1] = 10

    [29]:

    c.plot_state();

    [30]:

    c.AvalancheLoop()

    [30]:

    {'AvalancheSize': 10000,
     'NumberOfReleases': 28191892,
     'number_of_iterations': 6920}

    [31]:

    c.plot_state();

#### Does pattern depend on height of load?

    [32]:

    for i in [4, 5, 6, 10, 20]:
        mdl = BTW(100)
        mdl.values[1:-1,1:-1] = i
        mdl.AvalancheLoop()
        mdl.plot_state();

As we can see patter is changeing depending of pile height. Is there some maximum load?

    [33]:

    mdl = BTW(100)
    mdl.values[1:-1,1:-1] = 40
    mdl.AvalancheLoop()
    mdl.plot_state();

    [34]:

    mdl = BTW(100)
    mdl.values[1:-1,1:-1] = 50
    mdl.AvalancheLoop()
    mdl.plot_state();

Pattern depends on load height.

#### Loading model in the single point at centr.

    [35]:

    for i in [10**2, 10**3, 10**4, 10**5]:
        mdl = BTW(101)
        mdl.values[50,50] = i
        mdl.AvalancheLoop()
        mdl.plot_state();

### How exponent dependce on size of system?

    [36]:

    zip([10, 50, 100, 500],[1,1,1,1])

    [36]:

    <zip at 0x7fb61fdc2320>

    [37]:

    for l, w in zip([10, 50, 100, 200], [3*100**2, 3*100**2, 3*100**2, 5*100**2]):
        mdl = BTW(L = l, save_every = 100)
        mdl.run(2*w, wait_for_n_iters = w)
        mdl.get_exponent(low = 2, high = 1.3*l)

    Waiting for wait_for_n_iters=30000 iterations before collecting data. This should let the system thermalize.

    y = 6951.140 exp(-0.9537 x)
    Waiting for wait_for_n_iters=30000 iterations before collecting data. This should let the system thermalize.

    y = 5898.736 exp(-1.0567 x)
    Waiting for wait_for_n_iters=30000 iterations before collecting data. This should let the system thermalize.

    y = 5673.058 exp(-1.0858 x)
    Waiting for wait_for_n_iters=50000 iterations before collecting data. This should let the system thermalize.

    y = 9486.563 exp(-1.2028 x)

## Forest fire model

    [1]:

    from SOC.models import Forest

    /progs/miniconda3/lib/python3.7/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
      import pandas.util.testing as tm

    [2]:

    T chance_of_thunder = 0.001
    model = Forest(p=0.04, f=chance_of_thunder, L=50, save_every = 1)

    model.plot_state(False)

    [2]:

    [3]:

    model.run(1000, wait_for_n_iters=1000)

    Waiting for wait_for_n_iters=1000 iterations before collecting data. This should let the system thermalize.

- 0 - ash
- 1 - tree
- 2 - burning

    [4]:

    model.animate_states(notebook=True, interval=100)

Your browser does not support the video tag.

    [5]:

    model.data_df

    [5]:

| AvalancheSize | NumberOfReleases | number_of_iterations
--- | --- | --- | ---
0 | 0 | 0 | 75
1 | 0 | 0 | 74
2 | 0 | 0 | 67
3 | 0 | 0 | 74
4 | 0 | 0 | 77
... | ... | ... | ...
995 | 0 | 0 | 68
996 | 0 | 0 | 60
997 | 0 | 0 | 48
998 | 0 | 0 | 47
999 | 0 | 0 | 47

1000 rows × 3 columns

    [14]:

    model.data_df.number_of_iterations.plot()

    [14]:

    <matplotlib.axes._subplots.AxesSubplot at 0x7f9ea2102f90>

## Input/Output

### Saving state of simulation

implementaion is based on zarr library

    [1]:

    from SOC.models import BTW

    a = BTW(15)
    a.run(100)
    a.plot_state()
    root = a.save()
    root.tree()

    100%|██████████| 100/100 [00:03<00:00, 26.37it/s]

    [1]: