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.models.physics.criticality import ActiveFieldElement
>>> ActiveFieldElement('site').e
0.2

Source code in src/skcomplex/models/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.models.physics.criticality import AvalancheInteraction
>>> AvalancheInteraction('avalanche', target='pile', threshold=2, mode='random').share
1

Source code in src/skcomplex/models/physics/criticality/_avalanche.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,
    )

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.models.physics.criticality import CracklingNoiseAnalysis
>>> activity = [0, 1, 0, 1, 2, 1, 0, 1, 0, 2, 3, 2, 0]
>>> crackling = CracklingNoiseAnalysis(min_count=1).check().prepare(activity).run()
>>> round(crackling.exponent_, 2)
1.55

Source code in src/skcomplex/models/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

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.models.physics.criticality import FireInteraction
>>> FireInteraction('fire', target='forest').stochastic
True

Source code in src/skcomplex/models/physics/criticality/_fire.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.models.physics.criticality import ForestElement
>>> ForestElement('site').growth
0.05

Source code in src/skcomplex/models/physics/criticality/_forest.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.models.physics.criticality import LearnerElement
>>> LearnerElement('agent').memory
0.75

Source code in src/skcomplex/models/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.models.physics.criticality import LearningAgentElement
>>> LearningAgentElement('agent').alpha
0.5

Source code in src/skcomplex/models/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.models.physics.criticality import SandpileCellElement
>>> SandpileCellElement('site').height
0

Source code in src/skcomplex/models/physics/criticality/_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

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 model, with the activity neurons.rho and the energy neurons.nu.

Examples:

>>> from skcomplex.models.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/models/physics/criticality/_absorbing.py

@register_model(
    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 model, with the activity `neurons.rho` and the energy `neurons.nu`.

    Examples:
        >>> from skcomplex.models.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 model, 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 model 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 model, 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.models.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/models/physics/criticality/_evolutionary_learning.py

@register_model(
    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
    model, 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 model 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 model, 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.models.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])

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 model, with the covers forest.cover.

Examples:

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

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

@register_model(
    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 model, with the covers `forest.cover`.

    Examples:
        >>> from skcomplex.models.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 model, with the covers forest.cover.

Examples:

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

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

@register_model(
    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 model, with the covers `forest.cover`.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.physics.criticality import hexagonal_forest_fire
        >>> hexagonal_forest_fire((10, 10)).elements[0].check().prepare().size_
        100
        >>> hexagonal_forest_fire(graph=nx.cycle_graph(12)).elements[0].check().prepare().size_
        12
    """
    space = (
        (LatticeSpace(tuple(shape), geometry='triangular') if graph is None else NetworkSpace(graph)).check().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 model, with the heights pile.height.

Examples:

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

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

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

    Examples:
        >>> from skcomplex.models.physics.criticality import manna_sandpile
        >>> manna_sandpile((10, 10)).elements[0].check().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 model, with the heights pile.height.

Examples:

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

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

@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')],
    )

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 model, 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.models.physics.criticality import sk_game
>>> couplings = sk_game(10).elements[0].values['couplings']
>>> bool((couplings == couplings.T).all())
True

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

@register_model(
    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 model, 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.models.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])