Bak-Tang-Wiesenfeld sandpile

The Bak-Tang-Wiesenfeld sandpile on a lattice with open edges, started from random stable heights.

• Name: sandpile
• Field: physics/criticality
• Formalism: lattice dynamics
• References: Bak, Tang, and Wiesenfeld 1987; Charbonneau 2017

Factory

@register_model(
    field='physics/criticality',
    tags=('lattice', 'avalanches', 'power law'),
    references=('Bak, Tang, and Wiesenfeld 1987', 'Charbonneau 2017'),
    defaults={'until': 20000, 'dt': 100},
)
def sandpile(shape: Sequence[int] = (50, 50), *, random_state: int | np.random.Generator | None = 0) -> System:
    """Return the Bak-Tang-Wiesenfeld sandpile on a lattice with open edges, started from random stable heights.

    Its avalanches are recorded as `avalanche` events with their sizes, the number of topplings, their durations,
    and their areas.

    Args:
        shape:
            The number of rows and columns.

        random_state:
            The seed or generator of the initial heights.

    Returns:
        system:
            The model, with the heights `pile.height`.

    Examples:
        >>> from skcomplex.models.physics.criticality import sandpile
        >>> sandpile((10, 10)).elements[0].check().prepare().size_
        100
    """
    heights = np.random.default_rng(random_state).integers(0, 4, size=tuple(shape)).astype(float)
    space = LatticeSpace(tuple(shape), boundary='fixed')
    return System(
        elements=[LatticeElement('pile', element=SandpileCellElement, space=space, values={'height': heights})],
        interactions=[AvalancheInteraction('avalanche', target='pile')],
    )

Components

class SandpileCellElement(Element):
    """Site of a sandpile, which holds a number of grains.

    An `AvalancheInteraction` drops the grains and topples the sites that hold too many.

    Args:
        label:
            The label of the site or the lattice group.

        visual:
            How the site is drawn, or `None` for `default_visual_`.

        height:
            The initial number of grains.

    Attributes:
        default_visual_:
            The visual of the site when it sets none.

    Examples:
        >>> from skcomplex.models.physics.criticality import SandpileCellElement
        >>> SandpileCellElement('site').height
        0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        height: Annotated[float, StateSpec(domain=NN)] = 0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.height = height

class AvalancheInteraction(CascadeInteraction):
    """Slowly driven sandpile: every step drops one grain on a random site, and the full sites topple in avalanches.

    A site that holds `threshold` grains topples. In the mode `'redistribute'` it sends one grain to each of its
    neighbours, as in the Bak-Tang-Wiesenfeld sandpile, and in the mode `'random'` every grain goes to a neighbour drawn
    at random, as in the Manna sandpile. Grains that topple over a fixed edge leave the pile, and every avalanche is
    recorded as an event under the label of the interaction.

    Args:
        label:
            The label of the interaction and of its avalanche events.

        target:
            The path of the lattice group of sites.

        state:
            The name of the state that holds the grains.

        threshold:
            The number of grains at which a site topples.

        mode:
            `'redistribute'` or `'random'`.

    Examples:
        >>> from skcomplex.models.physics.criticality import AvalancheInteraction
        >>> AvalancheInteraction('avalanche', target='pile', threshold=2, mode='random').share
        1
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        state: str = 'height',
        threshold: int = 4,
        mode: str = 'redistribute',
    ) -> None:
        super().__init__(
            label,
            target=target,
            state=state,
            threshold=threshold,
            share=1,
            mode=mode,
            drive=self._drop_grain,
        )

    @staticmethod
    def _drop_grain(heights: Array, *, random: np.random.Generator) -> Array:
        """Return the heights with one grain added on a random site."""
        following = heights.copy()
        following[tuple(random.integers(size) for size in heights.shape)] += 1
        return following