Skip to content

Lenia and its expanded universe

Lenia, a continuous automaton whose ring-kernel convolutions grow gliding, lifelike solitons.

  • Name: lenia
  • Field: computer_science/automata
  • Formalism: discrete map
  • References: Chan 2019; Chan 2020

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('lenia')
entry = describe('lenia')
trajectory = create_simulation('lenia').prepare().run().trajectory_

Factory

@register_system(
    field='computer_science/automata',
    tags=('lattice', 'continuous automaton', 'convolution', 'soliton', 'artificial life', 'field'),
    references=('Chan 2019', 'Chan 2020'),
    defaults={'until': 200, 'dt': 1},
    observables=lambda system: [
        ObservableSpec(f'mass_{channel}', group='world', variable=channel, reduce=_mass) for channel in CHANNELS
    ],
)
def lenia(
    shape: Sequence[int] = (64, 64),
    *,
    initial: npt.ArrayLike | None = None,
    R: float = 13.0,
    T: float = 10.0,
    mu: float = 0.15,
    sigma: float = 0.015,
    beta: Sequence[float] = (1.0,),
    channels: int = 1,
    kernels: Sequence[Kernel] | None = None,
) -> System:
    """Return Lenia, a continuous automaton whose ring-kernel convolutions grow gliding, lifelike solitons.

    The world is a torus of sites with values in `[0, 1]`, and one step of `1 / T` adds the growth map of the
    convolution of the world with a ring kernel and clips the result, `A <- [A + G(K * A) / T]` in `[0, 1]`, with
    `G(u) = 2 exp(-(u - mu)^2 / (2 sigma^2)) - 1` and the kernel of radius `R`, whose concentric rings peak at `beta`.
    The expanded universe of Chan (2020) has several channels and several kernels
    `(source, destination, r, h, beta, mu, sigma)`, each of radius `r R` with its own growth map, which reads its
    source channel and grows its destination with the weight `h` over the sum of the weights into that destination,
    `A_j <- [A_j + sum_k (h_k / h_j) G_k(K_k * A_i) / T]`, on worlds of three or more dimensions. By default the world
    holds Orbium unicaudatus, which glides across the torus at a steady speed while its mass stays near 71. Every
    step of the system is one update, so a unit of Lenia's time takes `T` steps.

    Args:
        shape:
            The number of sites along every axis of the torus, two or more axes.

        initial:
            The initial values, an array of the shape shared by every channel, or one per channel, or `None` for
            Orbium in the middle of a two-dimensional world.

        R:
            The radius of the kernels in sites.

        T:
            The number of steps per unit time, the inverse of the increment of the update.

        mu:
            The centre of the growth map of the single kernel.

        sigma:
            The width of the growth map of the single kernel.

        beta:
            The peaks of the concentric rings of the single kernel.

        channels:
            The number of channels of the world, from one to three.

        kernels:
            The kernels `(source, destination, r, h, beta, mu, sigma)` of the expanded universe, with the channels as
            indices, or `None` for the single kernel of `mu`, `sigma`, and `beta` on the first channel.

    Returns:
        system:
            The system, with the values `world.A`, `world.B`, and `world.C` of the three channels of the sites, the
            last two at 0 for a world of one channel. Its standard observables `mass_A`, `mass_B`, and `mass_C` are the
            totals of the values of every channel.

    Raises:
        ParameterError:
            The world has fewer than two axes or not one to three channels, Orbium does not fit a default world, the
            initial values do not have the shape, or a kernel names a channel the world does not have.

    Examples:
        >>> from skcomplex.systems.computer_science.automata import lenia
        >>> lenia().elements[0].values['A'].shape
        (64, 64)
    """
    shape = tuple(shape)
    if len(shape) < MIN_AXES or not 1 <= channels <= len(CHANNELS):
        error_msg = f'Lenia needs two or more axes and one to three channels, but got {shape!r} and {channels!r}.'
        raise ParameterError(error_msg)
    values = _initial_values(shape, channels, initial)
    rings = ((0, 0, 1.0, 1.0, tuple(beta), mu, sigma),) if kernels is None else tuple(kernels)
    if any(not 0 <= index < channels for kernel in rings for index in kernel[:2]):
        error_msg = f'Every kernel needs a source and a destination among the {channels} channels.'
        raise ParameterError(error_msg)
    world = FieldElement(
        'world',
        element=LeniaSiteElement,
        space=LatticeSpace(shape),
        values={CHANNELS[index]: values[index] for index in range(channels)},
    )
    updates = [
        LeniaUpdateInteraction(
            f'update_{CHANNELS[index]}',
            target='world',
            channel=CHANNELS[index],
            kernels=[
                (CHANNELS[source], radius * R, weight, tuple(peaks), centre, width)
                for source, destination, radius, weight, peaks, centre, width in rings
                if destination == index
            ],
            T=T,
            n_dim=len(shape),
        )
        for index in range(channels)
        if any(kernel[1] == index for kernel in rings)
    ]
    return System(elements=[world], interactions=updates, time=DiscreteTime())

Components

Elements:

  • LeniaSiteElement, from Element

Interactions:

  • LeniaUpdateInteraction, from Interaction
class LeniaSiteElement(Element):
    """Site of a Lenia world, with a value in `[0, 1]` for each of up to three channels.

    The updates of the world change the channels, and a world of one channel leaves `B` and `C` at 0.

    Args:
        label:
            The label of the site or the field.

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

        A:
            The initial value of the first channel.

        B:
            The initial value of the second channel.

        C:
            The initial value of the third channel.

    Examples:
        >>> from skcomplex.systems.computer_science.automata import LeniaSiteElement
        >>> LeniaSiteElement('site', A=0.5).A
        0.5
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        A: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
        B: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
        C: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.A = A
        self.B = B
        self.C = C

class LeniaUpdateInteraction(Interaction):
    """Update of a channel of a Lenia world by the growth maps of the convolutions of channels with ring kernels.

    Every kernel `(source, radius, h, beta, mu, sigma)` reads the state of its source channel, and its concentric
    rings peak at `beta`. At the distance `q` in units of its radius, with `B` rings, it is `beta_i K_C(B q - i)` on
    the ring `i = floor(B q)`, with the bump `K_C(x) = exp(4 - 1 / (x (1 - x)))`, normalised to a sum of 1 over the
    sites. Its growth map is `G(u) = 2 exp(-(u - mu)^2 / (2 sigma^2)) - 1`, and one step adds the growths weighted by
    `h` over the sum of the weights and clips the channel, `A <- [A + sum_k (h_k / h) G_k(K_k * A_k) / T]` in
    `[0, 1]`.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the world.

        channel:
            The name of the state of the channel that the update changes.

        kernels:
            The kernels `(source, radius, h, beta, mu, sigma)`, with the name of the state of the source channel and
            the radius in sites.

        T:
            The number of steps per unit time.

        n_dim:
            The number of dimensions of the world.

    Examples:
        >>> from skcomplex.systems.computer_science.automata import LeniaUpdateInteraction
        >>> LeniaUpdateInteraction('update', target='world').T
        10.0
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        channel: str = 'A',
        kernels: Sequence[RingKernel] = (('A', 13.0, 1.0, (1.0,), 0.15, 0.015),),
        T: float = 10.0,
        n_dim: int = 2,
    ) -> None:
        super().__init__(label)
        self.target = target
        self.channel = channel
        self.kernels = kernels
        self.T = T
        self.n_dim = n_dim

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the update of the channel by the kernels.

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

        Returns:
            laws:
                The laws of the component.
        """
        distance = sqrt(sum(displacement(axis) ** 2 for axis in range(self.n_dim)))
        total = sum(kernel[2] for kernel in self.kernels)
        growth: Expression = 0
        for source, radius, weight, peaks, mu, sigma in self.kernels:
            kernel = self._shell(distance / radius, peaks) / self._norm(radius, peaks)
            potential = convolve(getattr(s.target, source), kernel, math.ceil(radius))
            growth += weight / total * (2 * exp(-((potential - mu) ** 2) / (2 * sigma**2)) - 1)
        value = getattr(s.target, self.channel)
        return [UpdateLaw(value, clip(value + growth / self.T, 0, 1))]

    @staticmethod
    def _shell(q: Expression, peaks: Sequence[float]) -> Expression:
        """Return the rings of a kernel at the distance `q` in units of its radius, zero outside the unit ball."""
        rings = len(peaks)
        shell: Expression = 0
        for index, peak in enumerate(peaks):
            x = rings * q - index
            inside = both(x > 0, x < 1)
            shell += where(inside, peak * exp(4 - 1 / where(inside, x * (1 - x), 1)), 0)
        return shell

    def _norm(self, radius: float, peaks: Sequence[float]) -> float:
        """Return the sum of the rings of a kernel over the sites within its radius."""
        steps = np.arange(-math.ceil(radius), math.ceil(radius) + 1)
        grid = np.meshgrid(*[steps] * self.n_dim, indexing='ij')
        q: Array = np.sqrt(sum(axis**2 for axis in grid)) / radius
        rings = len(peaks)
        values = np.zeros_like(q)
        for index, peak in enumerate(peaks):
            x = rings * q - index
            inside = (x > 0) & (x < 1)
            values += np.where(inside, peak * np.exp(4 - 1 / np.where(inside, x * (1 - x), 1)), 0)
        return float(values.sum())

Equations

world.A[i]: A_world(t + 1) == where(A_world(t) + 0.2*e^(-2222.222222222222*(convolve(A_world(t), 0.00490726701172641*where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 1.0*e^(-1/where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2)*(-0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) + 1), 1) + 4), 0), 13) - 0.15)^2) - 0.1 < 0, 0, where(A_world(t) + 0.2*e^(-2222.222222222222*(convolve(A_world(t), 0.00490726701172641*where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 1.0*e^(-1/where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2)*(-0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) + 1), 1) + 4), 0), 13) - 0.15)^2) - 0.1 > 1, 1, A_world(t) + 0.2*e^(-2222.222222222222*(convolve(A_world(t), 0.00490726701172641*where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 1.0*e^(-1/where(both(0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) > 0, 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) < 1), 0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2)*(-0.07692307692307693*sqrt(displacement_0^2 + displacement_1^2) + 1), 1) + 4), 0), 13) - 0.15)^2) - 0.1))
world.B[i]: B_world(t + 1) == B_world(t)
world.C[i]: C_world(t + 1) == C_world(t)