Contact process

The contact process on a periodic lattice, which survives above a critical infection rate near 1.65.

• Name: contact_process
• Field: physics/statistical_mechanics
• Formalism: lattice dynamics
• References: Harris 1974; Marro and Dickman 1999; Pastor-Satorras et al. 2015

Factory

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'stochastic', 'absorbing state', 'phase transition', 'epidemics'),
    references=('Harris 1974', 'Marro and Dickman 1999', 'Pastor-Satorras et al. 2015'),
    defaults={'until': 200, 'dt': 1},
)
def contact_process(
    shape: Sequence[int] = (64, 64),
    *,
    rate: float = 3.0,
    tau: float = 0.1,
    infected: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the contact process on a periodic lattice, which survives above a critical infection rate near 1.65.

    Below the critical rate the infection dies out and the lattice falls into the absorbing healthy state.

    Args:
        shape:
            The number of rows and columns.

        rate:
            The infection rate.

        tau:
            The time slice of one step.

        infected:
            The initial fraction of infected sites.

        random_state:
            The seed or generator of the initial infections.

    Returns:
        system:
            The model, with the infections `lattice.infected`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import contact_process
        >>> contact_process((10, 10), rate=1.0).elements[0].values['rate']
        1.0
    """
    shape = tuple(shape)
    states = (np.random.default_rng(random_state).random(shape) < infected).astype(float)
    space = LatticeSpace(shape)
    return System(
        elements=[
            LatticeElement(
                'lattice',
                element=ContactSiteElement,
                space=space,
                values={'infected': states, 'rate': rate, 'tau': tau},
            ),
        ],
    )

Components

class ContactSiteElement(Element):
    """Site of the contact process, which infects a neighbour at the rate `rate / 4` and recovers at rate 1.

    One step is a time slice `tau`, in which an infected site recovers with probability `1 - exp(-tau)` and a
    healthy site is infected with probability `1 - exp(-rate tau f)`, where `f` is the fraction of its infected
    neighbours.

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

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

        infected:
            The initial state, 1 for infected and 0 for healthy.

        rate:
            The infection rate.

        tau:
            The time slice of one step.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ContactSiteElement
        >>> ContactSiteElement('site').rate
        3.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        infected: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
        rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
        tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    ) -> None:
        super().__init__(label, visual=visual)
        self.infected = infected
        self.rate = rate
        self.tau = tau

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the stochastic infection and recovery.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        return [
            LocalRuleLaw(s.infected, VonNeumannNeighborhood(center=True), ContactSiteElement._spread, stochastic=True),
        ]

    @staticmethod
    def _spread(infected: Array, *around: Array, rate: float | Array, tau: float | Array, random: Array) -> Array:
        """Return the next infections of the contact process over a time slice `tau`."""
        share = sum(around) / len(around)
        infects = random < -np.expm1(-rate * tau * share)
        recovers = random < -np.expm1(-tau)
        return np.where(infected == 1, np.where(recovers, 0.0, 1.0), np.where(infects, 1.0, 0.0))