Skip to content

Vicsek system with metric, topological, vision-cone, and delayed alignment

The Vicsek system: particles that move at a constant speed along the mean heading of their neighbours.

  • Name: vicsek
  • Field: physics/active_matter
  • Formalism: discrete map
  • References: Vicsek et al. 1995; Argun et al. 2021; Jensen 2023; Vicsek and Zafeiris 2012

Select this system

from skcomplex.systems import create_simulation, create_system, describe

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

Factory

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'alignment', 'discrete time', 'neighborhood'),
    references=('Vicsek et al. 1995', 'Argun et al. 2021', 'Jensen 2023', 'Vicsek and Zafeiris 2012'),
    defaults={'until': 200.0, 'dt': 1.0},
    observables=lambda system: _orders('flock'),
)
def vicsek(
    n_particles: int = 400,
    *,
    size: float = 10.0,
    v0: float = 0.5,
    eta: float = 0.1,
    seed: int = 0,
    neighbors: int | None = None,
    vision: float | None = None,
    delay: int = 0,
    memory: int = 5,
) -> System:
    """Return the Vicsek system: particles that move at a constant speed along the mean heading of their neighbours.

    At every step particle `i` takes the direction of the sum of the headings of the particles within the unit
    distance, itself included, plus a uniform angle in `[-eta pi, eta pi]`, and moves `v0` along it in a periodic box.
    At low noise and high density the particles align into a flock, and the polar order parameter
    `|<exp(i theta)>|` grows from about `1 / sqrt(N)` to near 1. Raising the noise or lowering the density destroys
    the order.

    Two variants of Argun et al. (2021, Section 8.3) make the interaction non-metric or non-reciprocal. With
    `neighbors`, a particle aligns with its `neighbors` nearest particles and itself at any distance, a topological
    interaction that need not be reciprocal, and the clusters grow with `neighbors`. With `vision`, a particle sees only
    the particles within the unit distance whose direction from it makes an angle of at most `vision` with its heading,
    a vision cone that breaks reciprocity. The two variants do not combine, since the nearest neighbourhood of the core
    has no cone.

    With a positive `delay` `h`, the sensory delay of Section 8.2, a particle aligns with the headings its neighbours
    had `h` steps ago, `theta_j(n + 1) = <theta_k(n - h)> + noise`. Argun et al. report that a short delay stabilizes
    the clusters, and a long delay hinders them. A negative delay is a prediction: the headings of the neighbours `|h|`
    steps ahead are extrapolated by a least-squares line through their last `memory` headings, and it hinders
    clustering. Every particle then keeps its last 25 headings, so the delay lies in `[-25, 25]`, and before the start
    they equal the initial heading.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box, in units of the interaction radius.

        v0:
            The distance a particle moves in one step.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

        seed:
            The seed of the uniform initial positions and headings.

        neighbors:
            The number of nearest particles a particle aligns with, or `None` for the particles within the unit
            distance.

        vision:
            The half-angle of the vision cone in `(0, pi]`, or `None` for the whole disk.

        delay:
            The delay of the alignment in steps, negative for a prediction.

        memory:
            The number of past headings a prediction fits, from 2 to 26.

    Returns:
        system:
            The system, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Raises:
        ParameterError:
            Both `neighbors` and `vision` are given, or the delay or the memory exceeds the 25 kept headings.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import vicsek
        >>> vicsek(10).elements[0].size
        10
        >>> type(vicsek(10, neighbors=4).interactions[1].neighborhood).__name__
        'NearestNeighborhood'
    """
    if neighbors is not None and vision is not None:
        error_msg = 'A Vicsek flock aligns with its nearest neighbours or within a vision cone, not both.'
        raise ParameterError(error_msg)
    if abs(delay) > MEMORY or not FIT_POINTS <= memory <= MEMORY + 1:
        error_msg = (
            f'The delay {delay!r} and the memory {memory!r} need |delay| <= {MEMORY} and 2 <= memory <= {MEMORY + 1}.'
        )
        raise ParameterError(error_msg)
    neighborhood: BaseNeighborhood | None = None
    if neighbors is not None:
        neighborhood = NearestNeighborhood(neighbors, center=True)
    if vision is not None:
        neighborhood = ConeNeighborhood(1.0, vision, direction='theta', center=True)
    alignment = HeadingAlignmentInteraction(
        'alignment',
        target='flock',
        neighborhood=neighborhood,
        latency=delay,
        memory=memory,
    )
    element = VicsekParticleElement if delay == 0 else DelayedVicsekParticleElement
    return _flock(element, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)

Components

Elements:

  • VicsekParticleElement, from Element

Interactions:

  • PropulsionStepInteraction, from MapInteraction
  • HeadingAlignmentInteraction, from CouplingInteraction
class VicsekParticleElement(Element):
    """Particle of a Vicsek flock, whose heading becomes the mean heading of its neighbours plus noise every step.

    `HeadingAlignmentInteraction` sets the new heading to the circular mean of the headings within the unit distance,
    itself included, plus a uniform angle in `[-eta pi, eta pi]`. A propulsion step moves the particle along its
    heading.

    Args:
        label:
            The label of the particle or the group.

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

        x:
            The initial position.

        theta:
            The initial heading.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

    Attributes:
        default_visual_:
            The visual of the particle when it sets none.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import VicsekParticleElement
        >>> VicsekParticleElement('bird').eta
        0.1
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
        theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
        eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.theta = theta
        self.eta = eta

class PropulsionStepInteraction(MapInteraction):
    """Step `x(t + 1) = x + speed (cos theta, sin theta)` of a position in the plane along a heading, in discrete time.

    The step reads the heading at the start of the step, as every map of the step does.

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

        target:
            The path of the particle or the group whose positions move.

        speed:
            The distance moved in one step.

        state:
            The name of the vector state that holds the position.

        heading:
            The name of the state that holds the heading.

    Examples:
        >>> from types import SimpleNamespace
        >>> from sage.all__sagemath_symbolics import var, vector
        >>> from skcomplex.systems.physics.common.particles import PropulsionStepInteraction
        >>> step = PropulsionStepInteraction('propulsion', target='flock').prepare()
        >>> target = SimpleNamespace(x=vector([var('x1'), var('x2')]), theta=var('theta'))
        >>> [law.value for law in step.run(SimpleNamespace(target=target, speed=var('v'))).laws_]
        [(v*cos(theta) + x1, v*sin(theta) + x2)]
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
        state: str = 'x',
        heading: str = 'theta',
    ) -> None:
        super().__init__(label, target=target, state=state, function=self._step, mode='increment')
        self.speed = speed
        self.heading = heading

    def _step(self, namespace: SimpleNamespace) -> FreeModuleElement:
        """Return the step `speed (cos theta, sin theta)`."""
        heading = getattr(namespace.target, self.heading)
        return namespace.speed * vector(SR, [cos(heading), sin(heading)])

class HeadingAlignmentInteraction(CouplingInteraction):
    """Vicsek alignment, the turn of every particle to the circular mean heading of its neighbours plus noise.

    Every step the heading `theta` of every particle becomes the circular mean of the headings within the unit
    distance, itself included, plus a uniform angle in `[-eta pi, eta pi]` with the parameter `eta` of the particle.
    With `nematic`, every heading is flipped to the side nearer the particle's own heading first, as for particles
    with an axis rather than a direction. It is the map form of coupling with the circular mean over a radius. Another
    neighbourhood, such as the `k` nearest particles or a vision cone, replaces the unit disk.

    With a positive `latency` `h`, every particle reads the headings of its neighbours `h` steps ago. With a negative
    latency it reads their headings predicted `|h|` steps ahead, by a least-squares line through their last `memory`
    headings, unwrapped by taking every change between two steps as the smallest angle. A latent alignment needs
    particles that keep enough past headings, such as `DelayedVicsekParticleElement`.

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

        target:
            The path of the group of particles.

        nematic:
            Whether every heading is flipped to the side nearer the particle's own heading.

        neighborhood:
            The neighbourhood of every particle, itself included, or `None` for the unit disk.

        latency:
            The number of steps by which the headings of the neighbours are read late, or predicted ahead when
            negative.

        memory:
            The number of past headings a prediction fits, at least 2.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import HeadingAlignmentInteraction
        >>> HeadingAlignmentInteraction('alignment', target='flock', nematic=True).reduce
        'circular-mean'
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        nematic: bool = False,
        neighborhood: BaseNeighborhood | None = None,
        latency: int = 0,
        memory: int = 5,
    ) -> None:
        super().__init__(
            label,
            target=target,
            read='theta',
            write='theta',
            function=self._get_heading,
            reduce='circular-mean',
            weights=self._get_side,
            form='map',
            combine=self._get_noisy,
            neighborhood=RadiusNeighborhood(1.0, center=True) if neighborhood is None else neighborhood,
        )
        self.nematic = nematic
        self.latency = latency
        self.memory = memory

    def _get_heading(self, namespace: SimpleNamespace) -> Expression:
        """Return the heading of the neighbour, delayed or predicted."""
        heading = namespace.other.theta
        if self.latency >= 0:
            return lag(heading, self.latency)
        unwrapped = [heading]
        for step in range(1, self.memory):
            change = lag(heading, step - 1) - lag(heading, step)
            unwrapped.append(unwrapped[-1] - atan2(sin(change), cos(change)))
        times = [-step for step in range(self.memory)]
        mean_time = sum(times) / self.memory
        mean_heading = sum(unwrapped) / self.memory
        spread = sum((time - mean_time) ** 2 for time in times)
        slope = sum((time - mean_time) * (value - mean_heading) for time, value in zip(times, unwrapped, strict=True))
        return mean_heading + slope / spread * (-self.latency - mean_time)

    def _get_side(self, namespace: SimpleNamespace) -> Expression:
        """Return -1 for a neighbour heading away from the particle when nematic, and 1 otherwise."""
        if not self.nematic:
            return SR(1)
        return where(cos(namespace.other.theta - namespace.member.theta) >= 0, 1, -1)

    @staticmethod
    def _get_noisy(namespace: SimpleNamespace, reduction: Target) -> Target:
        """Return the mean heading plus a uniform angle in `[-eta pi, eta pi]`."""
        eta = namespace.member.eta
        return reduction + uniform('noise', -eta * math.pi, eta * math.pi)

Equations

flock.x[i, 0]: x_flock_0(t + 1) == speed_propulsion*cos(theta_flock(t)) + x_flock_0(t)
flock.x[i, 1]: x_flock_1(t + 1) == speed_propulsion*sin(theta_flock(t)) + x_flock_1(t)
flock.theta[i]: theta_flock(t + 1) == 6.283185307179586*eta_flock*uniform_noise_alignment - 3.141592653589793*eta_flock + arctan2(_neighbours_1, _neighbours_0)