Grégoire-Chaté cohesive flocks
The Grégoire-Chaté model: aligning particles that also attract and repel, forming cohesive flocks.
- Name:
gregoire_chate - Field:
physics/active_matter - Formalism: discrete map
- References: Grégoire and Chaté 2004; Vicsek and Zafeiris 2012
Factory
@register_model(
field='physics/active_matter',
tags=('active matter', 'flocking', 'cohesion', 'discrete time', 'neighborhood'),
references=('Grégoire and Chaté 2004', 'Vicsek and Zafeiris 2012'),
defaults={'until': 300.0, 'dt': 1.0},
observables=lambda system: _orders('flock'),
)
def gregoire_chate(
n_particles: int = 200,
*,
size: float = 5.0,
v0: float = 0.05,
eta: float = 0.1,
alpha: float = 1.0,
beta: float = 0.5,
seed: int = 0,
) -> System:
"""Return the Grégoire-Chaté model: aligning particles that also attract and repel, forming cohesive flocks.
At every step particle `i` takes the direction of `alpha sum_j exp(i theta_j) + beta sum_j f(r_ij) e_ij +
eta n_i xi`, over the `n_i` particles within the unit distance, where `e_ij` points from `i` to `j`, `xi` is a
random unit vector, and the force `f` is a strong repulsion below `r_c = 0.2`, zero at the preferred distance
`r_e = 0.5`, and attractive up to `r_a = 0.8`. The particle then moves `v0` along the new heading. At low noise and
moderate cohesion the particles gather into cohesive, aligned flocks, and the polar order parameter grows towards
1. A strong cohesion locks them in place, and the order is lost.
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 strength of the vectorial noise.
alpha:
The strength of the alignment.
beta:
The strength of the cohesion.
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 gregoire_chate
>>> gregoire_chate(10, beta=2.0).elements[0].values['beta']
2.0
"""
return _flock(CohesiveParticleElement, n_particles, size, seed, v0=v0, eta=eta, alpha=alpha, beta=beta)
Components
class CohesiveParticleElement(Element):
"""Particle of a Grégoire-Chaté flock, which aligns, repels close neighbours, and attracts distant ones.
The new heading is the direction of `alpha sum_j exp(i theta_j) + beta sum_j f(r_ij) e_ij + eta n_i xi` over the
`n_i` particles within the unit distance, where `e_ij` points from `i` to `j`, `xi` is a random unit vector, and the
force `f` is a strong repulsion below `r_c`, zero at `r_e`, and attractive up to `r_a`. 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 strength of the vectorial noise.
alpha:
The strength of the alignment.
beta:
The strength of the cohesion.
r_c:
The distance below which the neighbours repel strongly.
r_e:
The preferred distance.
r_a:
The distance of the strongest attraction.
Attributes:
default_visual_:
The visual of the particle when it sets none.
Examples:
>>> from skcomplex.models.physics.active_matter import CohesiveParticleElement
>>> CohesiveParticleElement('bird').r_e
0.5
"""
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,
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
r_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = CLOSE,
r_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
r_a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.8,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.theta = theta
self.eta = eta
self.alpha = alpha
self.beta = beta
self.r_c = r_c
self.r_e = r_e
self.r_a = r_a
def _attraction(self, distance: Expression, namespace: SimpleNamespace) -> Expression:
"""Return the cohesion force, repulsive below `r_c` and towards `r_e` beyond it."""
s = namespace
return where(distance < s.r_c, -REPULSION, (distance - s.r_e) / (4 * (s.r_a - s.r_e)))
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the alignment, the cohesion, and the vectorial noise of the heading.
Args:
namespace:
The namespace of the symbols the laws can use, with the neighbour `namespace.other`.
Returns:
laws:
The update of the heading.
"""
s = namespace
near, others = RadiusNeighborhood(1.0, center=True), RadiusNeighborhood(1.0)
dx, dy = s.other.x[0] - s.x[0], s.other.x[1] - s.x[1]
distance = sqrt(dx**2 + dy**2)
force = self._attraction(distance, s)
count = neighbors(1, over=near)
angle = uniform('noise', -math.pi, math.pi)
along_x = (
s.alpha * neighbors(cos(s.other.theta), over=near)
+ s.beta * neighbors(force * dx / distance, over=others)
+ s.eta * count * cos(angle)
)
along_y = (
s.alpha * neighbors(sin(s.other.theta), over=near)
+ s.beta * neighbors(force * dy / distance, over=others)
+ s.eta * count * sin(angle)
)
return [UpdateLaw(s.theta, atan2(along_y, along_x))]
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)])
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) == arctan2(_neighbours_2*eta_flock*sin(6.283185307179586*uniform_noise_flock - 3.141592653589793) + _neighbours_3*alpha_flock + _neighbours_4*beta_flock, _neighbours_2*eta_flock*cos(6.283185307179586*uniform_noise_flock - 3.141592653589793) + _neighbours_0*alpha_flock + _neighbours_1*beta_flock)