Drossel-Schwabl forest fire

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

• Name: forest_fire
• Field: physics/criticality
• Formalism: lattice dynamics
• References: Drossel and Schwabl 1992; Charbonneau 2017

Factory

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

Components

class ForestElement(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.

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

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

        cover:
            The initial cover, 0, 1, or 2.

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

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

    Examples:
        >>> from skcomplex.models.physics.criticality import ForestElement
        >>> ForestElement('site').growth
        0.05
    """

    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

class FireInteraction(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.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the lattice or network group of sites.

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

    Examples:
        >>> from skcomplex.models.physics.criticality import FireInteraction
        >>> FireInteraction('fire', target='forest').stochastic
        True
    """

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

    @staticmethod
    def _burn(cover: Array, *around: Array, growth: float | Array, lightning: float | Array, random: Array) -> Array:
        """Return the covers after trees catch fire, fires burn out, and trees grow."""
        burning = np.any(np.stack(around) == BURNING, axis=0) if around else np.zeros(np.shape(cover), dtype=bool)
        ignites = burning | (random < lightning)
        tree = np.where(ignites, float(BURNING), 1.0)
        empty = np.where(random < growth, 1.0, 0.0)
        return np.where(cover == BURNING, 0.0, np.where(cover == 1, tree, empty))