Lecture 13: Self-organized criticality
======================================

Before this class you should:

.. include:: prep13.txt

Before next class you should:

.. include:: prep14.txt

Note taker: Ken Phanthavong


Critical Systems
----------------

A **critical point** is a point in a system where it is in-between states.
This in-between state is often called the **critical state**.
Systems under a critical state are called **critical systems**.

Most systems require a parameter that changes the state to be tuned precisely to reach a critical state.
When a system has **self-organized criticality**, the system will tend toward a critical state on its own
without the parameter being tuned.


Sand Pile Model
---------------

The first example of a model that showed the behavior of self organized criticality was the **sand pile model**
by Bak, Tang, and Weisenfeld. The sand pile model is not meant to realistically simulate the behavior of sand piles,
but to show the general behaviors of systems with large numbers of interacting elements and critical systems.

The sand pile model is a 2D cellular automata model where each cell is the slope at that location.
When a cell exceeds some threshold slope, the sand topples and transfers one unit of slope
to each of the neighbouring cells. Single grains of sand are dropped at a time to perturb the system.
When a perturbation causes one or more cells to topple, it is called an **avalanche**.

The following implementation of the sand pile has the constructor set all the cells to an equal level.

.. code-block:: python

    class SandPile(Cell2D):

        kernel = np.array([[0, 1, 0],
                           [1,-4, 1],
                           [0, 1, 0]])

        def __init__(self, n, m, level=9):
            self.array = np.ones((n, m)) * level

|

The :func:`step` method topples all cells above the threshold and distributes its value to its neighbours.
The :func:`run` method calls the :func:`step` method continuously until there are no cells left that can be toppled (when it reaches equilibrium).

.. code-block:: python

    def step(self, K=3):
        toppling = self.array > K
        num_toppled = np.sum(toppling)
        c = correlate2d(toppling, self.kernel, mode='same')
        self.array += c
        return num_toppled

    def run(self):
        total = 0
        for i in itertools.count(1):
            num_toppled = self.step()
            total += num_toppled
            if num_toppled == 0:
                return i, total

|

The :func:`drop` method adds a value of 1 to a random cell.

.. code-block:: python

    def drop(self):
        a = self.array
        n, m = a.shape
        index = np.random.randint(n), np.random.randint(m)
        a[index] += 1

|

Here is an example of initializing and running a sand pile until it reaches equilibrium.
The level set to every cell is 10 in a 20x20 cell grid.

.. code-block:: python

    pile = SandPile(n=20, level=10)
    pile.run()

.. figure:: images/lecture13/sandpile_run.svg

    The sand pile in equilibrium.


Heavy-Tailed Distribution
-------------------------

Critical systems have heavy tailed distributions in some of their properties.
Sometimes these follow the power law, which is a type of heavy tailed distribution that is scale invariant.

The distribution of a sand pile model can be shown to have a heavy-tailed distribution in its avalanche size (S) and duration (T).
This can be measured by performing many random grain drops and calculating the probabilities of avalanche sizes and durations.

.. code-block:: python

    pile2 = SandPile(n=50, level=30)
    pile2.run()
    iters = 100000
    res = [pile2.drop_and_run() for _ in range(iters)]
    T, S = np.transpose(res)
    T = T[T>1]
    S = S[S>0]
    pmfT = Pmf(T)
    pmfS = Pmf(S)

.. figure:: images/lecture13/sandpile_heavy_tailed1.svg

    Probability mass functions of sand pile avalanche size and duration.

.. figure:: images/lecture13/sandpile_heavy_tailed2.svg

    Probability mass functions of sand pile avalanche size and duration on log-log scale to show power law similarity.


Fractals
--------

Another property of critical systems is that they often result in fractal geometry forming.
This can be shown in an initialized sand pile under equilibrium. The CA looks like it has fractal patterns.

.. code-block:: python

    pile3 = SandPile(n=131, level=22)
    pile3.run()

.. figure:: images/lecture13/sandpile_fractal1.svg

    Sand pile in equilibrium showing fractal-like patterns.


Under equilibrium, there are only 4 levels in the sand pile CA: levels 0-3.
The fractal patterns and dimensions can be found for cells at each level.

.. code-block:: python

    def draw_four(viewer, levels=range(4)):
        thinkplot.preplot(rows=2, cols=2)
        a = viewer.viewee.array

        for i, level in enumerate(levels):
            thinkplot.subplot(i+1)
            viewer.draw_array(a==level, vmax=1)

.. figure:: images/lecture13/sandpile_fractal2.svg

    Levels 0-3 of the sand pile from the top left to the bottom right.


The bounding box method is used to measure the amount of cells in a growing box on the sand pile.
The fractal dimension can be estimated with the number of cells for each box size.

.. code-block:: python

    def count_cells(a):
        n, m = a.shape
        end = min(n, m)

        res = []
        for i in range(1, end, 2):
            top = (n-i) // 2
            left = (m-i) // 2
            box = a[top:top+i, left:left+i]
            total = np.sum(box)
            res.append((i, i**2, total))

        return np.transpose(res)

    res = count_cells(pile.array==level)
    steps, steps2, cells = res

.. figure:: images/lecture13/sandpile_fractal3.svg

    Fractal dimensions of levels 0-3 of the sand pile from the top left to the bottom right.


Pink Noise
----------

:math:`1/f` noise, also known as "pink noise" is a type of behavior present in signals of critical systems.
**Signals** are properties or quantities that change with time.
These signals can be represented as the combination of sinusoidal waves/functions
with different magnitudes, essentially breaking them down into different frequencies with different strengths.
These strengths or magnitudes are also called power, and a **power spectrum** is a function that shows the power at each
frequency. The power of pink noise is equal to :math:`1/f`, :math:`f` being the frequency.

The Welch method implemented in the SciPy function calculates the power spectrum given a sequence of values
and the time between each value, ``nperseg``. The power spectrum is plotted for the number of toppled cells of the earlier sand pile.

.. code-block:: python

    from scipy.signal import welch

    signal = pile2.toppled_seq
    nperseg = 2048
    freqs, spectrum = welch(signal, nperseg=nperseg, fs=nperseg)

.. figure:: images/lecture13/sandpile_noise.svg

    Power spectrum of the number of toppled cells from grain drops on a log-log scale showing relation to pink noise.


Reductionism and Holism
-----------------------

Holistic Model

* Focuses on similar behaviors in different systems

* Contains only the elements necessary to describe those systems


Reductionist Model

* Focuses on realism

* Describes the fundamental parts of a system and their interactions


The sand pile is meant to simulate a behavior that appears in many systems that show self-organized criticality.
Self-organized criticality is a holistic model of naturally occurring critical systems.


Causation and Prediction
------------------------

Self-organized critical systems may not have a specific cause for rare and unexpected events since they
contain heavy-tailed distributions. Instead, it could be naturally occurring as a property of the system
as a whole rather than a specific interaction. This means that they are fundamentally unpredictable and unexplainable.