Couzin shoal with repulsion, alignment, and informed leaders
The Couzin shoal, in which a small informed minority leads a group that does not know who is informed.
- Name:
couzin_leadership - Field:
biology/collective_behaviour - Formalism: agent-based model
- References: Couzin et al. 2005; Mihailović et al. 2023
Factory
@register_model(
field='biology/collective_behaviour',
tags=('agents', 'schooling', 'leadership', 'collective decision', 'continuous space'),
references=('Couzin et al. 2005', 'Mihailović et al. 2023'),
defaults={'until': 300, 'dt': 1},
)
def couzin_leadership(
n_fish: int = 50,
*,
informed: float = 0.1,
omega: float = 0.5,
preferred: float = 0.0,
size: float = 100.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the Couzin shoal, in which a small informed minority leads a group that does not know who is informed.
The fish start in a disk of radius 5 at the center of a torus with random headings, and a fraction `informed` of
them prefers the direction at the angle `preferred` with the weight `omega`. Repulsion, attraction, and alignment
keep the shoal together, and the informed fish pull its mean heading toward their direction. The accuracy, the
closeness of the group's heading to the preferred direction, grows with the size of the minority, and a minority
of about 10% already steers a group of 50.
Args:
n_fish:
The number of fish.
informed:
The fraction of informed fish.
omega:
The weight of the preferred direction of the informed fish.
preferred:
The angle of the preferred direction.
size:
The side of the square torus.
random_state:
The seed or generator of the initial positions and headings.
Returns:
system:
The model, with the group `shoal`, the unit headings `vx` and `vy` of its fish, and the weight `omega` of
each, positive for the informed ones.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import couzin_leadership
>>> int((couzin_leadership(20, informed=0.2).elements[0].values['omega'] > 0).sum())
4
"""
rng = np.random.default_rng(random_state)
radius = 5.0 * np.sqrt(rng.uniform(0.0, 1.0, n_fish))
place = rng.uniform(0.0, 2 * math.pi, n_fish)
headings = rng.uniform(0.0, 2 * math.pi, n_fish)
weights = np.zeros(n_fish)
weights[: round(informed * n_fish)] = omega
locations = size / 2 + np.column_stack([radius * np.cos(place), radius * np.sin(place)])
shoal = GroupElement(
'shoal',
element=ShoalMemberElement,
size=n_fish,
values={'vx': np.cos(headings), 'vy': np.sin(headings), 'omega': weights, 'preferred': preferred, 'size': size},
location=[tuple(point) for point in locations],
)
return System(
elements=[shoal],
interactions=[
SchoolingInteraction('schooling', target='shoal', phase='turn'),
MapInteraction(
'deviation',
target='shoal',
state='deviation',
function=_draw_deviation,
form='agent',
schedule='member',
phase='turn',
),
SwimmingInteraction('swimming', target='shoal', phase='swim'),
],
space=ContinuousSpace(extent=(size, size), boundary='periodic'),
time=DiscreteTime(
schedule=PhasesScheduleSpec([PhaseSpec('turn', activation='simultaneous'), PhaseSpec('swim')]),
),
)
Components
class ShoalMemberElement(Element):
"""Fish of a Couzin shoal that avoids close neighbours, follows farther ones, and may know a preferred direction.
Every step `SchoolingInteraction` turns a fish toward its desired direction by at most `max_turn`. When neighbours
lie within the zone of repulsion `alpha`, it points away from them. Otherwise it points toward the neighbours
within the zone of interaction `rho` and along their headings and its own. An informed fish adds its preferred
direction with the weight `omega`, while a naive fish has `omega = 0`. It then draws a Gaussian angular noise
`deviation` of standard deviation `noise`, and `SwimmingInteraction` turns it by the noise and moves it at the
speed `speed`.
Args:
label:
The label of the fish or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
vx:
The horizontal component of the unit heading.
vy:
The vertical component of the unit heading.
omega:
The weight of the preferred direction, 0 for a naive fish.
deviation:
The angular noise of the last turn.
preferred:
The angle of the preferred direction of the informed fish.
alpha:
The radius of the zone of repulsion.
rho:
The radius of the zone of interaction.
speed:
The distance a fish swims per step.
max_turn:
The largest turning angle per step.
noise:
The standard deviation of the angular noise.
size:
The side of the square torus.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import ShoalMemberElement
>>> GroupElement('shoal', element=ShoalMemberElement, size=10, location=(5.0, 5.0)).size
10
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'omega'}),
*,
vx: Annotated[float, StateSpec(role='velocity')] = 1.0,
vy: Annotated[float, StateSpec(role='velocity')] = 0.0,
omega: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
deviation: Annotated[float, StateSpec()] = 0.0,
preferred: Annotated[float, ParameterSpec()] = 0.0,
alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
rho: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 6.0,
speed: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
max_turn: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
noise: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
size: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
) -> None:
super().__init__(label, visual=visual)
self.vx = vx
self.vy = vy
self.omega = omega
self.deviation = deviation
self.preferred = preferred
self.alpha = alpha
self.rho = rho
self.speed = speed
self.max_turn = max_turn
self.noise = noise
self.size = size
@classmethod
def _get_view(cls) -> type[ShoalMemberView]:
"""Return the view of a member of a shoal, `ShoalMemberView`."""
return ShoalMemberView
class SchoolingInteraction(CouplingInteraction):
"""Couzin's turn of every fish toward a desired direction set by zones around it, by at most `max_turn`.
When neighbours lie within the zone of repulsion `alpha`, the desired direction points away from them, and
repulsion overrides the other zones. Otherwise it points toward the neighbours within the zone of interaction `rho`
and along their headings and the fish's own, and an informed fish adds its preferred direction with the weight
`omega`. A fish with no neighbour keeps its heading. Offsets are the shortest ones on the torus of side `size`. It
is the agent form of coupling over two radii of different priorities, with a clamp of the turn.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group of fish.
phase:
The phase of the step the turn runs in, or `None` for the default phase.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import SchoolingInteraction
>>> [spec.priority for spec in SchoolingInteraction('schooling', target='shoal').radii]
[1, 0, 0]
"""
def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], phase: str | None = None) -> None:
super().__init__(
label,
target=target,
write=('vx', 'vy'),
form='agent',
combine=self._get_desired,
radii=(
RadiusSpec('alpha', self._get_repulsion, priority=1),
RadiusSpec('rho', self._get_attraction),
RadiusSpec('rho', self._get_alignment),
),
clamp=ClampSpec(max_angle='max_turn'),
phase=phase,
)
@staticmethod
def _unit(vector: Array) -> Array:
"""Return a vector scaled to unit length, or the vector itself when it is zero."""
norm = float(np.hypot(*vector))
return vector / norm if norm > 0 else vector
def _get_direction(self, agent: ShoalMemberView, other: ShoalMemberView) -> Array:
"""Return the unit direction from a fish to another along the shortest offset on the torus."""
size = agent.size
offset = np.asarray(other.location, dtype=float) - np.asarray(agent.location, dtype=float)
return self._unit((offset + size / 2) % size - size / 2)
def _get_repulsion(self, agent: ShoalMemberView, other: ShoalMemberView, model: ModelView) -> Array:
"""Return the direction away from a close neighbour."""
return -self._get_direction(agent, other)
def _get_attraction(self, agent: ShoalMemberView, other: ShoalMemberView, model: ModelView) -> Array:
"""Return the direction toward a neighbour."""
return self._get_direction(agent, other)
@staticmethod
def _get_alignment(agent: ShoalMemberView, other: ShoalMemberView, model: ModelView) -> Array:
"""Return the heading of a neighbour."""
return np.array([other.vx, other.vy])
def _get_desired(
self,
agent: AgentViewProtocol,
model: ModelViewProtocol,
reduction: Array | None,
priority: int | None,
) -> Array:
"""Return the desired direction from repulsion, or from attraction, alignment, and the preferred direction."""
fish = cast(ShoalMemberView, agent)
heading = np.array([fish.vx, fish.vy])
if reduction is None:
return heading
if priority == 1:
desired = self._unit(reduction)
else:
goal = np.array([math.cos(fish.preferred), math.sin(fish.preferred)])
desired = self._unit(self._unit(heading + reduction) + fish.omega * goal)
return desired if float(np.hypot(*desired)) > 0 else heading
class MapInteraction(Interaction):
"""Map of a state of discrete time, whose next value is a function, or the state plus it.
In the mode `'replace'` the next value of the state is `f`, and in the mode `'increment'` it is the state plus `f`.
Under the `'synchronous'` schedule every member of a bound group updates at once, from the states at the start of
the step. Under the `'sequential'` schedule the members update one after the other, in their order, and under the
`'random-sequential'` schedule in a new random order every step, each member reading the values the members before
it set. Under the `'member'` schedule of the agent form, every member updates in its own activation, together with
the rules of its element and the other member interactions, in the order the system lists them.
In the `'symbolic'` form, `function(namespace)` returns a symbolic expression, or a vector of them for a vector
state, where `namespace.target` holds the states and parameters of the target, `namespace.sources` the namespaces
of the sources, and `namespace.t` the time, and may hold draws such as `normal('noise')`. In the `'agent'` form,
the target is an agent group and `function(agent, model)` returns the value as a number, or as an array of the
shape of a vector state, so it can read the agent, its neighbours, and the history and storage of the model view.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the map sets.
state:
The name of the state of the target.
function:
The function `f(namespace)` in the symbolic form, or `f(agent, model)` in the agent form.
sources:
The paths of the elements, or members of groups such as `cells[0]`, that the function reads.
mode:
`'replace'` or `'increment'`.
schedule:
`'synchronous'`, `'sequential'`, `'random-sequential'`, or `'member'` for the agent form.
form:
`'symbolic'` for update laws of the states of the target, or `'agent'` for a rule over an agent group.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, the
function is not callable, or the agent form returns a value that does not match the state, raised by
`run` of the simulation.
Examples:
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import MapInteraction
>>> from skcomplex.core.spaces import DiscreteTime
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> from typing import Annotated, SupportsFloat
>>> class ValueElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.5):
... super().__init__(label)
... self.x = x
>>> logistic = MapInteraction(
... 'map', target='value', state='x', function=lambda s: 3 * s.target.x * (1 - s.target.x)
... )
>>> system = System(elements=[ValueElement('value')], interactions=[logistic], time=DiscreteTime())
>>> Simulation(system).check().prepare().run(until=2).trajectory_.values('value.x').tolist()
[0.5, 0.75, 0.5625]
>>> agents = MapInteraction('map', target='walkers', state='x', function=lambda a, m: a.x + 1, form='agent')
>>> type(agents.check().prepare().run(SimpleNamespace()).laws_[0]).__name__
'RuleLaw'
>>> from skcomplex.core.elements import GroupElement
>>> class WeightsElement(Element):
... def __init__(self, label, visual=None, *, w: Annotated[object, StateSpec(vector=True)] = 0.0):
... super().__init__(label, visual=visual)
... self.w = w
>>> group = GroupElement('weights', element=WeightsElement, size=1, values={'w': [1.0, 2.0]})
>>> doubling = MapInteraction(
... 'double', target='weights', state='w', function=lambda a, m: [2 * x for x in a.w], form='agent'
... )
>>> trajectory = Simulation(System(elements=[group], interactions=[doubling])).check().prepare().run(until=2)
>>> trajectory.trajectory_.agents('weights', 2).values['w'].tolist()
[[4.0, 8.0]]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str,
function: Callable[..., object],
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
mode: str = 'replace',
schedule: str = 'synchronous',
form: str = 'symbolic',
phase: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.function = function
self.sources = sources
self.mode = mode
self.schedule = schedule
self.form = form
self.phase = phase
def _check(self) -> None:
"""Raise an error for an unknown mode, schedule, or form, or a function that cannot be called.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, or the
function is not callable.
"""
super()._check()
self._check_variant('mode', self.mode, MODES, (BaseSpec,))
for name, value, variants in (
('schedule', self.schedule, SCHEDULES),
('form', self.form, FORMS),
):
if value not in variants:
error_msg = f'The {name} of map {self.label!r} is {value!r}. Use one of {", ".join(variants)}.'
raise ValidationError(error_msg)
if self.schedule == 'member' and self.form != 'agent':
error_msg = (
f'The map {self.label!r} updates every member in its activation, which only the agent form does.'
)
raise ValidationError(error_msg)
if not callable(self.function):
error_msg = f'The function of map {self.label!r} must be callable, but it is {self.function!r}.'
raise ValidationError(error_msg)
def _apply_mode(self, current: Value, value: Value) -> Value:
"""Return the next value of a state from its current value and the value of the function, by the mode."""
if self.mode == 'increment':
return current + value
if self.mode == 'replace':
return value
raise self._get_variant_error(self.mode, 'mode', '_apply_mode')
def _get_vector(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> Value:
"""Return the next value of a vector state of an agent."""
current = np.asarray(getattr(agent, self.state), dtype=np.float64)
value = np.asarray(self.function(agent, model), dtype=np.float64)
if value.shape != current.shape:
error_msg = (
f'The function of map {self.label!r} returns the shape {value.shape} for the state {self.state!r} of '
f'the shape {current.shape}. Return an array of the shape of the state.'
)
raise ValidationError(error_msg)
return self._apply_mode(current, value)
def _get_value(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> Value:
"""Return the next value of the state of an agent, a number or the array of a vector state."""
if np.ndim(getattr(agent, self.state)) > 0:
return self._get_vector(agent, model)
found = self.function(agent, model)
if not isinstance(found, SupportsFloat):
error_msg = f'The function of map {self.label!r} returns {found!r} for an agent. Return a number.'
raise ValidationError(error_msg)
value = float(found)
return self._apply_mode(float(getattr(agent, self.state)), value)
def _set_value(self, agent: AgentViewProtocol, value: Value) -> None:
"""Set the state of an agent to its next value."""
setattr(agent, self.state, value)
def _update_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Set the state of one member to its next value."""
if agent.alive:
self._set_value(agent, self._get_value(agent, model))
def update(self, model: ModelViewProtocol) -> None:
"""Set the state of every agent of the target to its next value, by the schedule.
Args:
model:
The model view of the interaction.
"""
self._update_agents(model, self.schedule, self._get_value, self._set_value)
def _get_order(self) -> object:
"""Return the order of the members of the symbolic form: none, their own order, or a uniform draw."""
if self.schedule == 'synchronous':
return None
if self.schedule == 'sequential':
return SR(0)
return uniform(re.sub(r'\W', '_', f'order_{self.label}'))
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the update of the state to its next value, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target` and `namespace.sources`.
Returns:
laws:
The update law of the state, or the rule that runs the map of the agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._update_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.update, phase=self.phase)]
value = getattr(self._get_member(namespace, 'target'), self.state)
mapped = cast(Value, self.function(namespace))
next_value = self._apply_mode(value, mapped)
return [UpdateLaw(value, next_value, increment=self.mode == 'increment', order=self._get_order())]
class SwimmingInteraction(MovementInteraction):
"""Swim of every fish by `speed` along its heading, turned first by its angular noise `deviation`.
The heading `(vx, vy)` stays a unit vector, and the fish moves on the torus of side `size`.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group of fish.
probability:
The probability that a fish swims in a step.
phase:
The phase of the step the swim runs in, or `None` for the default phase.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import SwimmingInteraction
>>> SwimmingInteraction('swimming', target='shoal', phase='swim').phase
'swim'
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
probability: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 1.0,
phase: str | None = None,
) -> None:
super().__init__(
label,
target=target,
step=self._swim,
probability=probability,
schedule='member',
phase=phase,
)
@staticmethod
def _swim(agent: AgentViewProtocol, model: ModelViewProtocol) -> tuple[float, ...]:
"""Turn the heading of a fish by its angular noise and return its position after one swim."""
fish = cast(ShoalMemberView, agent)
angle = math.atan2(fish.vy, fish.vx) + fish.deviation
fish.vx, fish.vy = math.cos(angle), math.sin(angle)
position = np.asarray(fish.location, dtype=float) + fish.speed * np.array([fish.vx, fish.vy])
return tuple(float(item) for item in np.mod(position, fish.size))