Criticality

Criticality: sandpiles, forest fires, absorbing-state transitions, learners at the edge, and avalanches.

ActiveFieldElement(label, visual=None, *, rho=0.5, nu=1.0, a=-1.0, b=1.0, c=1.0, d=1.0, e=0.2, f=0.0, g=0.0)

Bases: Element

Site of the Langevin equation of directed percolation with a conserved field, an activity rho and a field nu.

The activity grows at the rate a + b nu, saturates, diffuses, and fluctuates with the multiplicative noise e sqrt(rho) of an absorbing state, rho_t = (a + b nu) rho - c rho^2 + d laplacian(rho) + e sqrt(rho) eta. The field diffuses, gains f, and is consumed by the activity, nu_t = laplacian(nu) + f - g rho. Both stay nonnegative.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
rho	Annotated[float, StateSpec()]	The initial activity.	0.5
nu	Annotated[float, StateSpec()]	The initial field.	1.0
a	Annotated[float, ParameterSpec()]	The linear growth rate of the activity.	-1.0
b	Annotated[float, ParameterSpec()]	The growth of the activity per unit of the field.	1.0
c	Annotated[float, ParameterSpec()]	The saturation of the activity.	1.0
d	Annotated[float, ParameterSpec()]	The diffusion coefficient of the activity.	1.0
e	Annotated[float, ParameterSpec()]	The amplitude of the noise.	0.2
f	Annotated[float, ParameterSpec()]	The input of the field.	0.0
g	Annotated[float, ParameterSpec()]	The consumption of the field by the activity.	0.0

Attributes:

Name	Type	Description
default_visual_		The visual of the site when it sets none.

Examples:

>>> from skcomplex.systems.physics.criticality import ActiveFieldElement
>>> ActiveFieldElement('site').e
0.2

Source code in src/skcomplex/systems/physics/criticality/_absorbing.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    rho: Annotated[float, StateSpec()] = 0.5,
    nu: Annotated[float, StateSpec()] = 1.0,
    a: Annotated[float, ParameterSpec()] = -1.0,
    b: Annotated[float, ParameterSpec()] = 1.0,
    c: Annotated[float, ParameterSpec()] = 1.0,
    d: Annotated[float, ParameterSpec()] = 1.0,
    e: Annotated[float, ParameterSpec()] = 0.2,
    f: Annotated[float, ParameterSpec()] = 0.0,
    g: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.rho = rho
    self.nu = nu
    self.a = a
    self.b = b
    self.c = c
    self.d = d
    self.e = e
    self.f = f
    self.g = g

AvalancheInteraction(label, *, target, state='height', threshold=4, mode='redistribute')

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction and of its avalanche events.	required
target	Annotated[str, RoleSpec()]	The path of the lattice group of sites.	required
state	str	The name of the state that holds the grains.	'height'
threshold	int	The number of grains at which a site topples.	4
mode	str	'redistribute' or 'random'.	'redistribute'

Examples:

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

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

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,
    )

BranchingTreeElement(label, visual=None, *, p=0.3, z=1, g=0)

Bases: Element

Self-organized branching process: avalanches on a binary tree whose branching probability the avalanches tune.

The probability p that an active site relaxes into two active sites of the next generation, the number z of active sites of the current generation of the avalanche, and its generation g. Map interactions advance the avalanche by one generation per step.

Parameters:

Name	Type	Description	Default
label	str	The label of the tree.	required
visual	VisualSpec | None	How the tree is drawn, or None for default_visual_.	None
p	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial branching probability.	0.3
z	Annotated[float, StateSpec(domain=NN)]	The initial number of active sites.	1
g	Annotated[float, StateSpec(domain=NN)]	The initial generation.	0

Attributes:

Name	Type	Description
default_visual_		The visual of the tree when it sets none.

Examples:

>>> from skcomplex.systems.physics.criticality import BranchingTreeElement
>>> BranchingTreeElement('tree').p
0.3

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    p: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.3,
    z: Annotated[float, StateSpec(domain=NN)] = 1,
    g: Annotated[float, StateSpec(domain=NN)] = 0,
) -> None:
    super().__init__(label, visual=visual)
    self.p = p
    self.z = z
    self.g = g

CracklingNoiseAnalysis(variables=None, *, threshold=0.0, durations=(), min_count=10, min_duration=1, size_exponent=None, duration_exponent=None)

Bases: EventAnalysis

Crackling-noise scaling of the avalanches of an activity series or of the relaxations of a lattice.

It finds the avalanches as the event analysis of kind 'excursions', then fits the exponent gamma of the mean size against the duration, <S>(T) ~ T^gamma, predicts it as (tau_T - 1) / (tau_S - 1) from the exponents of the size and duration laws, such as those FitAnalysis(law='power') fits, and collapses the mean shapes of the chosen durations, which agree after scaling at a critical point.

Parameters:

Name	Type	Description	Default
variables	Sequence[str] | None	The one variable of a trajectory to read, or None for the only scalar series.	None
threshold	float	The activity an avalanche exceeds at every step, subtracted from its profile.	0.0
durations	Sequence[int]	The durations whose mean shapes are compared.	()
min_count	int	The fewest avalanches of a duration used in the fit of the mean size.	10
min_duration	int	The shortest duration used in the fit of the mean size.	1
size_exponent	float | None	The exponent tau_S of the power law of the sizes, or None to leave out the prediction.	None
duration_exponent	float | None	The exponent tau_T of the power law of the durations, or None to leave out the prediction.	None

Attributes:

Name	Type	Description
exponent_		The exponent gamma of the mean size against the duration.
predicted_		The exponent (tau_T - 1) / (tau_S - 1), or nan without both exponents.
mean_shapes_		The mean shape of every chosen duration, scaled to unit mean height, on a grid of 0 to 1.
shape_spread_		The mean over the inner 80 percent of the scaled time of the relative difference between the largest and the smallest mean shape of the chosen durations, nan with fewer than 2 of them.

Examples:

>>> from skcomplex.systems.physics.criticality import CracklingNoiseAnalysis
>>> activity = [0, 1, 0, 1, 2, 1, 0, 1, 0, 2, 3, 2, 0]
>>> crackling = CracklingNoiseAnalysis(min_count=1).prepare(activity).run()
>>> round(crackling.exponent_, 2)
1.55

Source code in src/skcomplex/systems/physics/criticality/_crackling_noise.py

def __init__(
    self,
    variables: Sequence[str] | None = None,
    *,
    threshold: float = 0.0,
    durations: Sequence[int] = (),
    min_count: int = 10,
    min_duration: int = 1,
    size_exponent: float | None = None,
    duration_exponent: float | None = None,
) -> None:
    super().__init__('excursions', variables, threshold=threshold)
    self.durations = durations
    self.min_count = min_count
    self.min_duration = min_duration
    self.size_exponent = size_exponent
    self.duration_exponent = duration_exponent

EnergySiteElement(label, visual=None, *, energy=0.0)

Bases: Element

Site of the Zhang sandpile, which stores a continuous energy.

A cascade loads the energy in random packets and resets the sites above the threshold.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
energy	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial energy.	0.0

Attributes:

Name	Type	Description
default_visual_		The visual of the site when it sets none.

Examples:

>>> from skcomplex.systems.physics.criticality import EnergySiteElement
>>> EnergySiteElement('site').energy
0.0

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    energy: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.energy = energy

FireInteraction(label, *, target, neighborhood=NEAREST)

Bases: LocalRuleInteraction

Stochastic rule of the Drossel-Schwabl forest fire on a cover of 0 for empty, 1 for a tree, and 2 for burning.

Every step a tree catches fire from a burning neighbour or from lightning with the probability lightning of the site, a burning site becomes empty, and a tree grows on an empty site with the probability growth of the site.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the lattice or network group of sites.	required
neighborhood	BaseNeighborhood	The neighbours a fire spreads to, the von Neumann neighbourhood of a square lattice, the six neighbours of a triangular lattice, or the neighbours of a node, each with the site itself.	NEAREST

Examples:

>>> from skcomplex.systems.physics.criticality import FireInteraction
>>> FireInteraction('fire', target='forest').stochastic
True

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    neighborhood: BaseNeighborhood = NEAREST,
) -> None:
    super().__init__(
        label,
        target=target,
        state='cover',
        rule=self._burn,
        neighborhood=neighborhood,
        stochastic=True,
    )

ForestElement(label, visual=None, *, cover=0, growth=0.05, lightning=0.0001)

Bases: Element

Site of the Drossel-Schwabl forest fire: empty at 0, a tree at 1, or burning at 2.

A FireInteraction grows trees on empty sites with probability growth, sets them on fire from a burning neighbour or from lightning with probability lightning, and empties burning sites.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
cover	Annotated[float, StateSpec(domain=NN)]	The initial cover, 0, 1, or 2.	0
growth	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability that a tree grows on an empty site per step.	0.05
lightning	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability that lightning strikes a tree per step.	0.0001

Examples:

>>> from skcomplex.systems.physics.criticality import ForestElement
>>> ForestElement('site').growth
0.05

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    cover: Annotated[float, StateSpec(domain=NN)] = 0,
    growth: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    lightning: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0001,
) -> None:
    super().__init__(label, visual=visual)
    self.cover = cover
    self.growth = growth
    self.lightning = lightning

LearnerElement(label, visual=None, *, incentive=0.0, spin=1.0, temperature=1.0, fitness=0.0, reward=0.0, memory=0.75, max_temperature=2.0, couplings=None, evolving=True, settling=2)

Bases: Element

Agent of the evolutionary Sherrington-Kirkpatrick game, who learns a binary action and may be replaced.

The agent keeps an incentive Q, the value of action +1, and plays the action s = ±1 by the softmax policy at its temperature T. Every step the incentive becomes (1 - alpha) Q + alpha h for the local field h = Σ_j J_ij s_j, and the fitness is s Q. When the agent evolves and its fitness is the least of the group, it is replaced by a newcomer with no incentive and a temperature drawn uniformly up to max_temperature, who learns from the field at once. The reward s h of the last actions is kept.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
incentive	Annotated[float, StateSpec()]	The initial incentive Q.	0.0
spin	Annotated[float, StateSpec(domain=ZZ)]	The initial action, -1 or 1.	1.0
temperature	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial temperature T of the softmax policy.	1.0
fitness	Annotated[float, StateSpec()]	The initial fitness s Q.	0.0
reward	Annotated[float, StateSpec()]	The initial reward s h.	0.0
memory	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The weight alpha of the latest field in the incentive, so small values keep a long memory.	0.75
max_temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The largest temperature of a newcomer.	2.0
couplings	ArrayLike | None	The coupling matrix J, one row and one column per agent.	None
evolving	bool	Whether the least fit agent is replaced.	True
settling	int	The number of steps of learning alone after every replacement, so replacements happen when t is a multiple of settling + 1.	2

Examples:

>>> from skcomplex.systems.physics.criticality import LearnerElement
>>> LearnerElement('agent').memory
0.75

Source code in src/skcomplex/systems/physics/criticality/_evolutionary_learning.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    incentive: Annotated[float, StateSpec()] = 0.0,
    spin: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    temperature: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    fitness: Annotated[float, StateSpec()] = 0.0,
    reward: Annotated[float, StateSpec()] = 0.0,
    memory: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
    max_temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    couplings: npt.ArrayLike | None = None,
    evolving: bool = True,
    settling: int = 2,
) -> None:
    super().__init__(label, visual=visual)
    self.incentive = incentive
    self.spin = spin
    self.temperature = temperature
    self.fitness = fitness
    self.reward = reward
    self.memory = memory
    self.max_temperature = max_temperature
    self.couplings = couplings
    self.evolving = evolving
    self.settling = settling

LearningAgentElement(label, visual=None, *, sigma=1.0, Q=0.0, alpha=0.5, temperature=0.5, couplings=None)

Bases: Element

Agent of the SK game that plays an action of ±1 from a softmax policy over its incentive to play +1.

The incentive follows the reinforcement rule Q(t + 1) = (1 - alpha) Q(t) + alpha sum_j J_ij sigma_j(t), the discounted reward of playing +1, and the agent plays sigma = ±1 with probability 1 / (1 + exp(-+2 Q / T)).

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
sigma	Annotated[float, StateSpec(domain=ZZ)]	The initial action, -1 or 1.	1.0
Q	Annotated[float, StateSpec()]	The initial incentive.	0.0
alpha	Annotated[float, ParameterSpec(domain=open(0, oo))]	The memory parameter, the weight of the latest reward.	0.5
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The exploration temperature T of the softmax policy.	0.5
couplings	ArrayLike | None	The coupling matrix J, one row and one column per agent.	None

Examples:

>>> from skcomplex.systems.physics.criticality import LearningAgentElement
>>> LearningAgentElement('agent').alpha
0.5

Source code in src/skcomplex/systems/physics/criticality/_learning.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    sigma: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    Q: Annotated[float, StateSpec()] = 0.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    couplings: npt.ArrayLike | None = None,
) -> None:
    super().__init__(label, visual=visual)
    self.sigma = sigma
    self.Q = Q
    self.alpha = alpha
    self.temperature = temperature
    self.couplings = couplings

SandpileCellElement(label, visual=None, *, height=0)

Bases: Element

Site of a sandpile, which holds a number of grains.

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

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
height	Annotated[float, StateSpec(domain=NN)]	The initial number of grains.	0

Attributes:

Name	Type	Description
default_visual_		The visual of the site when it sets none.

Examples:

>>> from skcomplex.systems.physics.common import SandpileCellElement
>>> SandpileCellElement('site').height
0

Source code in src/skcomplex/systems/physics/common/_sandpile_cell.py

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

ToppleInteraction(label, *, target, threshold=4)

Bases: LocalRuleInteraction

Parallel toppling of a sandpile without drive, where every site at or above the threshold topples once a step.

Every step, every site that holds at least threshold grains loses threshold grains, and every site gains one grain from each toppling neighbour, all at once. On a periodic lattice no grain is lost, so the number of grains is conserved, and the activity either dies out in an absorbing configuration or lasts forever.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the lattice group of sites.	required
threshold	int	The number of grains at which a site topples, the number of its neighbours.	4

Examples:

>>> from skcomplex.systems.physics.criticality import ToppleInteraction
>>> ToppleInteraction('topple', target='pile').threshold
4

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], threshold: int = 4) -> None:
    super().__init__(
        label,
        target=target,
        state='height',
        rule=self._topple,
        neighborhood=VonNeumannNeighborhood(center=True),
    )
    self.threshold = threshold

directed_percolation_langevin(shape=(64, 64), *, a=-1.0, b=1.0, c=1.0, d=1.0, e=0.2, f=0.0, g=0.0, energy=1.2, density=0.5, random_state=0)

Return the Langevin field theory of an absorbing-state transition, an active density fed by an energy field.

The active density rho grows at the local rate a + b nu, saturates through -c rho^2, spreads with the diffusivity d, and fluctuates with the multiplicative noise e sqrt(rho), which vanishes in the absorbing state rho = 0. The energy field nu diffuses, gains the drive f, and loses the dissipation g rho, so rho_t = (a + b nu) rho - c rho^2 + d laplacian(rho) + e sqrt(rho) xi and nu_t = laplacian(nu) + f - g rho (Tian et al., Eqs. 3-4). Both fields stay nonnegative. Without drive and dissipation the energy is a fixed control parameter: below the threshold nu_c = -a / b activity dies out, and above it the density settles near (a + b nu) / c.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along every axis of the periodic lattice.	(64, 64)
a	float	The bare growth rate of the activity.	-1.0
b	float	The gain of the growth rate per unit of energy.	1.0
c	float	The saturation of the activity.	1.0
d	float	The diffusivity of the activity.	1.0
e	float	The amplitude of the multiplicative noise.	0.2
f	float	The drive of the energy.	0.0
g	float	The dissipation of the energy per unit of activity.	0.0
energy	float	The initial energy at every site.	1.2
density	float	The mean initial activity, drawn uniformly between 0 and twice this value.	0.5
random_state	int | Generator | None	The seed or generator of the initial activity.	0

Returns:

Name	Type	Description
system	System	The system, with the activity neurons.rho and the energy neurons.nu.

Examples:

>>> from skcomplex.systems.physics.criticality import directed_percolation_langevin
>>> float(directed_percolation_langevin((8, 8), energy=0.5).elements[0].values['nu'].mean())
0.5

Source code in src/skcomplex/systems/physics/criticality/_absorbing.py

@register_system(
    field='physics/criticality',
    tags=('field', 'absorbing state', 'directed percolation', 'multiplicative noise', 'spde'),
    references=('Tian et al. 2022', 'Henkel, Hinrichsen, and Lübeck 2008', 'Bonachela and Muñoz 2009'),
    defaults={'until': 100.0, 'dt': 1.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def directed_percolation_langevin(
    shape: Sequence[int] = (64, 64),
    *,
    a: float = -1.0,
    b: float = 1.0,
    c: float = 1.0,
    d: float = 1.0,
    e: float = 0.2,
    f: float = 0.0,
    g: float = 0.0,
    energy: float = 1.2,
    density: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Langevin field theory of an absorbing-state transition, an active density fed by an energy field.

    The active density `rho` grows at the local rate `a + b nu`, saturates through `-c rho^2`, spreads with the
    diffusivity `d`, and fluctuates with the multiplicative noise `e sqrt(rho)`, which vanishes in the absorbing state
    `rho = 0`. The energy field `nu` diffuses, gains the drive `f`, and loses the dissipation `g rho`, so
    `rho_t = (a + b nu) rho - c rho^2 + d laplacian(rho) + e sqrt(rho) xi` and `nu_t = laplacian(nu) + f - g rho`
    (Tian et al., Eqs. 3-4). Both fields stay nonnegative. Without drive and dissipation the energy is a fixed control
    parameter: below the threshold `nu_c = -a / b` activity dies out, and above it the density settles near
    `(a + b nu) / c`.

    Args:
        shape:
            The number of sites along every axis of the periodic lattice.

        a:
            The bare growth rate of the activity.

        b:
            The gain of the growth rate per unit of energy.

        c:
            The saturation of the activity.

        d:
            The diffusivity of the activity.

        e:
            The amplitude of the multiplicative noise.

        f:
            The drive of the energy.

        g:
            The dissipation of the energy per unit of activity.

        energy:
            The initial energy at every site.

        density:
            The mean initial activity, drawn uniformly between 0 and twice this value.

        random_state:
            The seed or generator of the initial activity.

    Returns:
        system:
            The system, with the activity `neurons.rho` and the energy `neurons.nu`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import directed_percolation_langevin
        >>> float(directed_percolation_langevin((8, 8), energy=0.5).elements[0].values['nu'].mean())
        0.5
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    values = {
        'rho': rng.uniform(0.0, 2 * density, shape),
        'nu': np.full(shape, energy),
        'a': a,
        'b': b,
        'c': c,
        'd': d,
        'e': e,
        'f': f,
        'g': g,
    }
    return System(
        elements=[FieldElement('neurons', element=ActiveFieldElement, space=LatticeSpace(shape), values=values)],
    )

evolutionary_sk_game(n_agents=256, *, memory=0.75, temperature=None, max_temperature=2.0, random_state=0)

Return the evolutionary Sherrington-Kirkpatrick game, which self-organises to the edge of ergodicity breaking.

After Lesmana et al. (2026, arXiv 2604.15669), N agents play ±1 on a Sherrington-Kirkpatrick landscape with symmetric Gaussian couplings of variance 1 / N. Each learns the incentive of its actions with memory alpha and plays by a softmax policy at its own temperature, and the least fit agent is replaced, as in the Bak-Sneppen system, by a newcomer whose temperature is drawn uniformly on [0, max_temperature]. The temperatures evolve to the edge between ergodic and frozen learning, the replacement cascades have sizes distributed as S^(-3/2), and the mean reward exceeds that of the game at any fixed temperature. With a fixed temperature, no agent is replaced and the system is the Sherrington-Kirkpatrick game of Galla and Farmer.

Parameters:

Name	Type	Description	Default
n_agents	int	The number of agents N.	256
memory	float	The weight alpha of the latest field in the incentive.	0.75
temperature	ArrayLike | None	The temperature of every agent, or one per agent, of the game without replacement, or None for the evolving game.	None
max_temperature	float	The largest temperature of a newcomer, and of the initial temperatures drawn uniformly.	2.0
random_state	int | Generator | None	The seed or generator of the couplings, the initial actions, and the initial temperatures.	0

Returns:

Name	Type	Description
system	System	The system, with the group agents of states incentive, spin, temperature, fitness, and reward, and the observables mean_reward, mean_temperature, and least_fitness.

Examples:

>>> from skcomplex.systems.physics.criticality import evolutionary_sk_game
>>> agents = evolutionary_sk_game(16).elements[0]
>>> couplings = agents.values['couplings']
>>> agents.size, bool((couplings == couplings.T).all()), float(couplings.trace())
(16, True, 0.0)

Source code in src/skcomplex/systems/physics/criticality/_evolutionary_learning.py

@register_system(
    field='physics/criticality',
    tags=('self-organized criticality', 'reinforcement learning', 'spin glass', 'extremal dynamics', 'discrete map'),
    references=('Lesmana, Feng, Chen, and Lai 2026', 'Bak and Sneppen 1993', 'Sherrington and Kirkpatrick 1975'),
    defaults={'until': 3000, 'dt': 1},
    observables=_observables,
)
def evolutionary_sk_game(
    n_agents: int = 256,
    *,
    memory: float = 0.75,
    temperature: npt.ArrayLike | None = None,
    max_temperature: float = 2.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the evolutionary Sherrington-Kirkpatrick game, which self-organises to the edge of ergodicity breaking.

    After Lesmana et al. (2026, arXiv 2604.15669), `N` agents play ±1 on a Sherrington-Kirkpatrick landscape with
    symmetric Gaussian couplings of variance `1 / N`. Each learns the incentive of its actions with memory `alpha` and
    plays by a softmax policy at its own temperature, and the least fit agent is replaced, as in the Bak-Sneppen
    system, by a newcomer whose temperature is drawn uniformly on `[0, max_temperature]`. The temperatures evolve to the
    edge between ergodic and frozen learning, the replacement cascades have sizes distributed as `S^(-3/2)`, and the
    mean reward exceeds that of the game at any fixed temperature. With a fixed temperature, no agent is replaced and
    the system is the Sherrington-Kirkpatrick game of Galla and Farmer.

    Args:
        n_agents:
            The number of agents `N`.

        memory:
            The weight `alpha` of the latest field in the incentive.

        temperature:
            The temperature of every agent, or one per agent, of the game without replacement, or `None` for the
            evolving game.

        max_temperature:
            The largest temperature of a newcomer, and of the initial temperatures drawn uniformly.

        random_state:
            The seed or generator of the couplings, the initial actions, and the initial temperatures.

    Returns:
        system:
            The system, with the group `agents` of states `incentive`, `spin`, `temperature`, `fitness`, and `reward`,
            and the observables `mean_reward`, `mean_temperature`, and `least_fitness`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import evolutionary_sk_game
        >>> agents = evolutionary_sk_game(16).elements[0]
        >>> couplings = agents.values['couplings']
        >>> agents.size, bool((couplings == couplings.T).all()), float(couplings.trace())
        (16, True, 0.0)
    """
    rng = np.random.default_rng(random_state)
    couplings = _sk_couplings(n_agents, rng)
    spins = rng.choice([-1.0, 1.0], size=n_agents)
    evolving = temperature is None
    if temperature is None:
        temperatures = rng.uniform(0.0, max_temperature, n_agents)
    else:
        temperatures = np.broadcast_to(np.asarray(temperature, dtype=float), (n_agents,)).copy()
    agents = GroupElement(
        'agents',
        element=LearnerElement,
        size=n_agents,
        values={
            'spin': spins,
            'temperature': temperatures,
            'memory': memory,
            'max_temperature': max_temperature,
            'couplings': couplings,
            'evolving': evolving,
        },
        visual=TEMPERATURE_VISUAL,
    )
    return System(elements=[agents])

fixed_energy_sandpile(shape=(32, 32), *, density=2.3, random_state=0)

Return the fixed-energy sandpile, a Bak-Tang-Wiesenfeld pile with no drive and no loss on a periodic lattice.

Every step, every site with at least four grains gives one grain to each of its four neighbours, all at once. The lattice is periodic, so the density of grains density is conserved, and it controls an absorbing phase transition. Below the critical density, about 2.125 on the square lattice, the topplings stop in an absorbing configuration, and above it the pile stays active forever. The density of active sites, the observable activity, is the order parameter. The grains start scattered uniformly at random over the sites.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(32, 32)
density	float	The mean number of grains per site.	2.3
random_state	int | Generator | None	The seed or generator of the initial grains.	0

Returns:

Name	Type	Description
system	System	The system, with the heights pile.height and the observable activity, the density of active sites.

Examples:

>>> from skcomplex.systems.physics.criticality import fixed_energy_sandpile
>>> fixed_energy_sandpile((8, 8), density=2.0).elements[0].prepare().size_
64

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

@register_system(
    field='physics/criticality',
    tags=('lattice', 'absorbing states', 'phase transition', 'conservation', 'sandpile'),
    references=('Dickman, Vespignani, and Zapperi 1998', 'Vespignani et al. 2000', 'Gros 2015'),
    defaults={'until': 500, 'dt': 1},
    observables=lambda system: [ObservableSpec('activity', _active, group='pile')],
)
def fixed_energy_sandpile(
    shape: Sequence[int] = (32, 32),
    *,
    density: float = 2.3,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the fixed-energy sandpile, a Bak-Tang-Wiesenfeld pile with no drive and no loss on a periodic lattice.

    Every step, every site with at least four grains gives one grain to each of its four neighbours, all at once.
    The lattice is periodic, so the density of grains `density` is conserved, and it controls an absorbing phase
    transition. Below the critical density, about 2.125 on the square lattice, the topplings stop in an absorbing
    configuration, and above it the pile stays active forever. The density of active sites, the observable `activity`,
    is the order parameter. The grains start scattered uniformly at random over the sites.

    Args:
        shape:
            The number of rows and columns.

        density:
            The mean number of grains per site.

        random_state:
            The seed or generator of the initial grains.

    Returns:
        system:
            The system, with the heights `pile.height` and the observable `activity`, the density of active sites.

    Examples:
        >>> from skcomplex.systems.physics.criticality import fixed_energy_sandpile
        >>> fixed_energy_sandpile((8, 8), density=2.0).elements[0].prepare().size_
        64
    """
    size = int(np.prod(shape))
    rng = np.random.default_rng(random_state)
    grains = rng.integers(size, size=round(density * size))
    heights = np.bincount(grains, minlength=size).reshape(tuple(shape)).astype(float)
    space = LatticeSpace(tuple(shape), boundary='periodic')
    return System(
        elements=[LatticeElement('pile', element=SandpileCellElement, space=space, values={'height': heights})],
        interactions=[ToppleInteraction('topple', target='pile', threshold=THRESHOLD)],
    )

forest_fire(shape=(100, 100), *, growth=0.05, lightning=0.0001, density=0.5, random_state=0)

Return the Drossel-Schwabl forest fire, where slow growth and rare lightning give fires of every size.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(100, 100)
growth	float	The probability that a tree grows on an empty site per step.	0.05
lightning	float	The probability that lightning strikes a tree per step.	0.0001
density	float	The initial fraction of sites with a tree.	0.5
random_state	int | Generator | None	The seed or generator of the initial forest.	0

Returns:

Name	Type	Description
system	System	The system, with the covers forest.cover.

Examples:

>>> from skcomplex.systems.physics.criticality import forest_fire
>>> forest_fire((10, 10)).elements[0].values['growth']
0.05

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

@register_system(
    field='physics/criticality',
    tags=('lattice', 'fire', 'stochastic', 'power law'),
    references=('Drossel and Schwabl 1992', 'Charbonneau 2017'),
    defaults={'until': 500, 'dt': 1},
)
def forest_fire(
    shape: Sequence[int] = (100, 100),
    *,
    growth: float = 0.05,
    lightning: float = 1e-4,
    density: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Drossel-Schwabl forest fire, where slow growth and rare lightning give fires of every size.

    Args:
        shape:
            The number of rows and columns.

        growth:
            The probability that a tree grows on an empty site per step.

        lightning:
            The probability that lightning strikes a tree per step.

        density:
            The initial fraction of sites with a tree.

        random_state:
            The seed or generator of the initial forest.

    Returns:
        system:
            The system, with the covers `forest.cover`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import forest_fire
        >>> forest_fire((10, 10)).elements[0].values['growth']
        0.05
    """
    cover = (np.random.default_rng(random_state).random(tuple(shape)) < density).astype(float)
    space = LatticeSpace(tuple(shape))
    return System(
        elements=[
            LatticeElement(
                'forest',
                element=ForestElement,
                space=space,
                values={'cover': cover, 'growth': growth, 'lightning': lightning},
            ),
        ],
        interactions=[FireInteraction('fire', target='forest')],
    )

hexagonal_forest_fire(shape=(80, 80), *, graph=None, growth=0.02, lightning=1e-05, density=0.5, cover=None, random_state=0)

Return the Drossel-Schwabl forest fire on a triangular lattice of hexagonal cells, or on the nodes of a graph.

Every site of the periodic triangular lattice has six neighbours, and every node of a graph has the neighbours along its edges. Slow growth and rare lightning give fires of every size.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns of the triangular lattice, unused on a graph.	(80, 80)
graph	Graph | None	The graph whose nodes are the sites, or None for the triangular lattice.	None
growth	float	The probability that a tree grows on an empty site per step.	0.02
lightning	float	The probability that lightning strikes a tree per step.	1e-05
density	float	The initial fraction of sites with a tree.	0.5
cover	Sequence[float] | None	The initial cover of every site, 0 for empty, 1 for a tree, and 2 for burning, in lattice shape or node order, or None for random trees at the given density.	None
random_state	int | Generator | None	The seed or generator of the initial forest.	0

Returns:

Name	Type	Description
system	System	The system, with the covers forest.cover.

Examples:

>>> import networkx as nx
>>> from skcomplex.systems.physics.criticality import hexagonal_forest_fire
>>> hexagonal_forest_fire((10, 10)).elements[0].prepare().size_
100
>>> hexagonal_forest_fire(graph=nx.cycle_graph(12)).elements[0].prepare().size_
12

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

@register_system(
    field='physics/criticality',
    tags=('lattice', 'network', 'fire', 'stochastic', 'power law'),
    references=('Drossel and Schwabl 1992', 'Argun et al. 2021'),
    defaults={'until': 400, 'dt': 2},
)
def hexagonal_forest_fire(
    shape: Sequence[int] = (80, 80),
    *,
    graph: nx.Graph | None = None,
    growth: float = 0.02,
    lightning: float = 1e-5,
    density: float = 0.5,
    cover: Sequence[float] | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Drossel-Schwabl forest fire on a triangular lattice of hexagonal cells, or on the nodes of a graph.

    Every site of the periodic triangular lattice has six neighbours, and every node of a graph has the neighbours
    along its edges. Slow growth and rare lightning give fires of every size.

    Args:
        shape:
            The number of rows and columns of the triangular lattice, unused on a graph.

        graph:
            The graph whose nodes are the sites, or `None` for the triangular lattice.

        growth:
            The probability that a tree grows on an empty site per step.

        lightning:
            The probability that lightning strikes a tree per step.

        density:
            The initial fraction of sites with a tree.

        cover:
            The initial cover of every site, 0 for empty, 1 for a tree, and 2 for burning, in lattice shape or node
            order, or `None` for random trees at the given density.

        random_state:
            The seed or generator of the initial forest.

    Returns:
        system:
            The system, with the covers `forest.cover`.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.systems.physics.criticality import hexagonal_forest_fire
        >>> hexagonal_forest_fire((10, 10)).elements[0].prepare().size_
        100
        >>> hexagonal_forest_fire(graph=nx.cycle_graph(12)).elements[0].prepare().size_
        12
    """
    space = (LatticeSpace(tuple(shape), geometry='triangular') if graph is None else NetworkSpace(graph)).prepare()
    size = int(space.size_ or 0)
    initial = (
        (np.random.default_rng(random_state).random(size) < density).astype(float)
        if cover is None
        else np.ravel(np.asarray(cover, dtype=float))
    )
    forest = LatticeElement(
        'forest',
        element=ForestElement,
        space=space,
        values={'cover': initial, 'growth': growth, 'lightning': lightning},
    )
    neighborhood = SITE_NEIGHBORHOODS['triangular' if graph is None else 'graph']
    return System(elements=[forest], interactions=[FireInteraction('fire', target='forest', neighborhood=neighborhood)])

manna_sandpile(shape=(32, 32), *, random_state=0)

Return the Manna sandpile, a stochastic sandpile of another universality class than Bak-Tang-Wiesenfeld.

The pile starts from random stable heights of 0 or 1 grain, and its avalanches are recorded as avalanche events with their sizes, durations, and areas.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(32, 32)
random_state	int | Generator | None	The seed or generator of the initial heights.	0

Returns:

Name	Type	Description
system	System	The system, with the heights pile.height.

Examples:

>>> from skcomplex.systems.physics.criticality import manna_sandpile
>>> manna_sandpile((10, 10)).elements[0].prepare().size_
100

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

@register_system(
    field='physics/criticality',
    tags=('lattice', 'avalanches', 'power law', 'stochastic'),
    references=('Manna 1991', 'Jensen 2023'),
    defaults={'until': 10000, 'dt': 50},
)
def manna_sandpile(shape: Sequence[int] = (32, 32), *, random_state: int | np.random.Generator | None = 0) -> System:
    """Return the Manna sandpile, a stochastic sandpile of another universality class than Bak-Tang-Wiesenfeld.

    The pile starts from random stable heights of 0 or 1 grain, and its avalanches are recorded as `avalanche`
    events with their sizes, durations, and areas.

    Args:
        shape:
            The number of rows and columns.

        random_state:
            The seed or generator of the initial heights.

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

    Examples:
        >>> from skcomplex.systems.physics.criticality import manna_sandpile
        >>> manna_sandpile((10, 10)).elements[0].prepare().size_
        100
    """
    heights = np.random.default_rng(random_state).integers(0, MANNA_THRESHOLD, 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', threshold=MANNA_THRESHOLD, mode='random')],
    )

sandpile(shape=(50, 50), *, random_state=0)

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.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(50, 50)
random_state	int | Generator | None	The seed or generator of the initial heights.	0

Returns:

Name	Type	Description
system	System	The system, with the heights pile.height.

Examples:

>>> from skcomplex.systems.physics.criticality import sandpile
>>> sandpile((10, 10)).elements[0].prepare().size_
100

Source code in src/skcomplex/systems/physics/criticality/_criticality.py

@register_system(
    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 system, with the heights `pile.height`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import sandpile
        >>> sandpile((10, 10)).elements[0].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')],
    )

self_organized_branching(p=0.3, *, generations=10)

Return the self-organized branching process of Zapperi, Lauritsen, and Stanley, a mean-field sandpile.

An avalanche starts from one active site at the top of a binary tree of generations generations. Every active site relaxes with the probability p into two active sites of the next generation, and stays still otherwise, so every generation has 2 Binomial(z, p) active sites. The sigma active sites that reach the last generation leave the system as dissipation. After every avalanche one unit of energy has entered and sigma have left, so the probability changes by (1 - sigma) / N with N = 2^(generations + 1) - 1 sites. A supercritical p makes large dissipating avalanches that lower it, and a subcritical one lets energy accumulate, so p organizes itself to the critical value p_c = 1 / 2, where the branching ratio 2 p is 1. The process advances one generation per step, so a step is a time of the avalanche and avalanches follow each other without pause.

Parameters:

Name	Type	Description	Default
p	float	The initial branching probability.	0.3
generations	int	The number of generations of the tree, at least 1.	10

Returns:

Name	Type	Description
system	System	The system, with the probability tree.p, the active sites tree.z, and the generation tree.g.

Examples:

>>> from skcomplex.systems.physics.criticality import self_organized_branching
>>> self_organized_branching(0.4).elements[0].p
0.4

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

@register_system(
    field='physics/criticality',
    tags=('branching process', 'avalanches', 'mean field', 'self-organized criticality', 'stochastic'),
    references=('Zapperi, Lauritsen, and Stanley 1995', 'Jensen 2023'),
    defaults={'until': 30000, 'dt': 1},
)
def self_organized_branching(p: float = 0.3, *, generations: int = 10) -> System:
    """Return the self-organized branching process of Zapperi, Lauritsen, and Stanley, a mean-field sandpile.

    An avalanche starts from one active site at the top of a binary tree of `generations` generations. Every active
    site relaxes with the probability `p` into two active sites of the next generation, and stays still otherwise, so
    every generation has `2 Binomial(z, p)` active sites. The `sigma` active sites that reach the last generation leave
    the system as dissipation. After every avalanche one unit of energy has entered and `sigma` have left, so the
    probability changes by `(1 - sigma) / N` with `N = 2^(generations + 1) - 1` sites. A supercritical `p` makes large
    dissipating avalanches that lower it, and a subcritical one lets energy accumulate, so `p` organizes itself to the
    critical value `p_c = 1 / 2`, where the branching ratio `2 p` is 1. The process advances one generation per step,
    so a step is a time of the avalanche and avalanches follow each other without pause.

    Args:
        p:
            The initial branching probability.

        generations:
            The number of generations of the tree, at least 1.

    Returns:
        system:
            The system, with the probability `tree.p`, the active sites `tree.z`, and the generation `tree.g`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import self_organized_branching
        >>> self_organized_branching(0.4).elements[0].p
        0.4
    """
    return System(
        elements=[BranchingTreeElement('tree', p=p)],
        interactions=[
            MapInteraction('branching', target='tree', state='z', function=partial(_branch, generations)),
            MapInteraction('generation', target='tree', state='g', function=partial(_advance, generations)),
            MapInteraction('tuning', target='tree', state='p', function=partial(_tune, generations)),
        ],
        time=DiscreteTime(),
    )

sk_game(n_agents=256, *, alpha=0.5, temperature=0.5, random_state=0)

Return the SK game, reinforcement learners whose rewards come from quenched Sherrington-Kirkpatrick couplings.

Every agent plays sigma_i = ±1 and earns sigma_i sum_j J_ij sigma_j, with symmetric Gaussian couplings of mean 0, variance 1 / N, and zero diagonal. It discounts its rewards with the memory alpha into the incentive Q_i and plays from a softmax policy at the temperature T (Lesmana et al., Eqs. 1-3). All agents start with zero incentives, so their first actions are random, and update in parallel. At low temperature the agents freeze into one of many ordered states with a high collective reward, and at high temperature they explore noisily with a low one.

Parameters:

Name	Type	Description	Default
n_agents	int	The number of agents N.	256
alpha	float	The memory parameter.	0.5
temperature	float	The exploration temperature of every agent.	0.5
random_state	int | Generator | None	The seed or generator of the couplings and of the initial actions.	0

Returns:

Name	Type	Description
system	System	The system, with the actions agents.sigma, the incentives agents.Q, and the coupling matrix in the values of the group. Its standard observable reward is the mean reward per agent.

Examples:

>>> from skcomplex.systems.physics.criticality import sk_game
>>> couplings = sk_game(10).elements[0].values['couplings']
>>> bool((couplings == couplings.T).all())
True

Source code in src/skcomplex/systems/physics/criticality/_learning.py

@register_system(
    field='physics/criticality',
    tags=('spin glass', 'reinforcement learning', 'game', 'disorder', 'discrete map'),
    references=('Lesmana et al. 2026', 'Sherrington and Kirkpatrick 1975', 'Watkins and Dayan 1992'),
    defaults={'until': 200, 'dt': 1},
    observables=_observables,
)
def sk_game(
    n_agents: int = 256,
    *,
    alpha: float = 0.5,
    temperature: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the SK game, reinforcement learners whose rewards come from quenched Sherrington-Kirkpatrick couplings.

    Every agent plays `sigma_i = ±1` and earns `sigma_i sum_j J_ij sigma_j`, with symmetric Gaussian couplings of mean
    0, variance `1 / N`, and zero diagonal. It discounts its rewards with the memory `alpha` into the incentive `Q_i`
    and plays from a softmax policy at the temperature `T` (Lesmana et al., Eqs. 1-3). All agents start with zero
    incentives, so their first actions are random, and update in parallel. At low temperature the agents freeze into
    one of many ordered states with a high collective reward, and at high temperature they explore noisily with a low
    one.

    Args:
        n_agents:
            The number of agents `N`.

        alpha:
            The memory parameter.

        temperature:
            The exploration temperature of every agent.

        random_state:
            The seed or generator of the couplings and of the initial actions.

    Returns:
        system:
            The system, with the actions `agents.sigma`, the incentives `agents.Q`, and the coupling matrix in the
            values of the group. Its standard observable `reward` is the mean reward per agent.

    Examples:
        >>> from skcomplex.systems.physics.criticality import sk_game
        >>> couplings = sk_game(10).elements[0].values['couplings']
        >>> bool((couplings == couplings.T).all())
        True
    """
    rng = np.random.default_rng(random_state)
    upper = np.triu(rng.normal(0.0, 1.0 / np.sqrt(n_agents), (n_agents, n_agents)), 1)
    agents = GroupElement(
        'agents',
        element=LearningAgentElement,
        size=n_agents,
        values={
            'sigma': rng.choice([-1.0, 1.0], size=n_agents),
            'Q': 0.0,
            'alpha': alpha,
            'temperature': temperature,
            'couplings': upper + upper.T,
        },
    )
    return System(elements=[agents])

zhang_sandpile(shape=(32, 32), *, threshold=1.0, random_state=0)

Return the Zhang sandpile, a continuous sandpile whose unstable sites give away all their energy.

Every step loads a site chosen at random with a packet of energy of uniform size in [0, threshold]. A site whose energy reaches the threshold relaxes: its energy falls to 0 and each of its four neighbours gains a quarter of it, all at once, until no site is unstable. Energy that would go beyond the open edges is lost. The pile organizes itself into a critical state with avalanches of every size, recorded as avalanche events with their sizes, durations, and areas, and the energies of the sites cluster near a few values. The energies start uniformly at random below the threshold.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(32, 32)
threshold	float	The energy at which a site relaxes, which is also the largest packet.	1.0
random_state	int | Generator | None	The seed or generator of the initial energies.	0

Returns:

Name	Type	Description
system	System	The system, with the energies pile.energy.

Examples:

>>> from skcomplex.systems.physics.criticality import zhang_sandpile
>>> zhang_sandpile((8, 8)).interactions[0].mode
'reset'

Source code in src/skcomplex/systems/physics/criticality/_sandpiles.py

@register_system(
    field='physics/criticality',
    tags=('lattice', 'avalanches', 'power law', 'continuous energy', 'sandpile'),
    references=('Zhang 1989', 'Jensen 2023'),
    defaults={'until': 10000, 'dt': 100},
)
def zhang_sandpile(
    shape: Sequence[int] = (32, 32),
    *,
    threshold: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Zhang sandpile, a continuous sandpile whose unstable sites give away all their energy.

    Every step loads a site chosen at random with a packet of energy of uniform size in `[0, threshold]`. A site whose
    energy reaches the threshold relaxes: its energy falls to 0 and each of its four neighbours gains a quarter of it,
    all at once, until no site is unstable. Energy that would go beyond the open edges is lost. The pile organizes
    itself into a critical state with avalanches of every size, recorded as `avalanche` events with their sizes,
    durations, and areas, and the energies of the sites cluster near a few values. The energies start uniformly at
    random below the threshold.

    Args:
        shape:
            The number of rows and columns.

        threshold:
            The energy at which a site relaxes, which is also the largest packet.

        random_state:
            The seed or generator of the initial energies.

    Returns:
        system:
            The system, with the energies `pile.energy`.

    Examples:
        >>> from skcomplex.systems.physics.criticality import zhang_sandpile
        >>> zhang_sandpile((8, 8)).interactions[0].mode
        'reset'
    """
    energies = np.random.default_rng(random_state).uniform(0.0, threshold, size=tuple(shape))
    space = LatticeSpace(tuple(shape), boundary='fixed')
    avalanche = CascadeInteraction(
        'avalanche',
        target='pile',
        state='energy',
        threshold=threshold,
        share=1 / 4,
        mode='reset',
        drive=partial(_load, threshold),
    )
    return System(
        elements=[LatticeElement('pile', element=EnergySiteElement, space=space, values={'energy': energies})],
        interactions=[avalanche],
    )