Skip to content

Vicsek model

The Vicsek model: 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

Factory

@register_model(
    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,
) -> System:
    """Return the Vicsek model: 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.

    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.

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

    Examples:
        >>> from skcomplex.models.physics.active_matter import vicsek
        >>> vicsek(10).elements[0].size
        10
    """
    alignment = HeadingAlignmentInteraction('alignment', target='flock')
    return _flock(VicsekParticleElement, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)

Components

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.models.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.models.physics.active_matter import PropulsionStepInteraction
        >>> step = PropulsionStepInteraction('propulsion', target='flock').check().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.

    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.

    Examples:
        >>> from skcomplex.models.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) -> 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),
        )
        self.nematic = nematic

    @staticmethod
    def _get_heading(namespace: SimpleNamespace) -> Expression:
        """Return the heading of the neighbour."""
        return namespace.other.theta

    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)