Flocking with active and passive agents
The active and passive agents of Charbonneau (2017, section 10.4) in a periodic box.
- Name:
active_passive_flocking - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Charbonneau 2017; Silverberg et al. 2013
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('active_passive_flocking')
entry = describe('active_passive_flocking')
trajectory = create_simulation('active_passive_flocking').prepare().run().trajectory_
Factory
@register_system(
field='physics/active_matter',
tags=('active matter', 'sde', 'flocking', 'phase separation', 'pair interaction'),
references=('Charbonneau 2017', 'Silverberg et al. 2013'),
defaults={'until': 50.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=TIME_STEP), 'backend': NumbaBackend()},
observables=_observables,
)
def active_passive_flocking(
n_agents: int = 342,
*,
n_active: int = 45,
size: float = 1.0,
radius: float = 0.025,
epsilon: float = 25.0,
exponent: float = 1.5,
flocking_radius: float = 0.1,
alpha: float = 1.0,
v0: float = 0.02,
mu: float = 10.0,
eta: float = 0.1,
time_step: float = TIME_STEP,
initial_speed: float = 1.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the active and passive agents of Charbonneau (2017, section 10.4) in a periodic box.
All agents repel each other within `2 r_0` with the force `epsilon (1 - r / (2 r_0))^exponent`, and all feel the
force `mu (v_t - |v|) v / |v|`, with the target speed `v_t = v0` for an active agent and 0 for a passive one. An
active agent also feels the flocking force `alpha V / |V|`, with `V` the sum of the velocities of the agents of both
kinds within the flocking radius, itself excluded, and a random force. Agents start at uniform random positions
with velocity components uniform up to `initial_speed`. The defaults are those of figure 10.7 at the compactness
`C = pi N r_0^2 / size^2 = 0.672`. The values `exponent=2.5`, `epsilon=100`, `mu=0.05`, `radius=1`,
`flocking_radius=4`, and `v0=1` are those of Silverberg et al. (2013).
Args:
n_agents:
The number of agents.
n_active:
The number of active agents.
size:
The side of the periodic box.
radius:
The radius `r_0` of an agent, so agents overlap below the distance `2 r_0`.
epsilon:
The strength of the repulsion.
exponent:
The exponent of the repulsion.
flocking_radius:
The radius within which an active agent reads the velocities of the others.
alpha:
The size of the flocking force on an active agent, 0 for no flocking.
v0:
The target speed of the active agents.
mu:
The strength of the propulsion.
eta:
The half-width of the uniform random force per step of the book. The noise is a Wiener process with the
same variance per step, of amplitude `eta sqrt(time_step / 3)`.
time_step:
The step of the solver that `eta` refers to.
initial_speed:
The largest size of a component of an initial velocity.
random_state:
The seed or generator of the initial positions and velocities.
Returns:
system:
The system, with the group `agents` that holds the active agents first, the positions `agents.x`, the
velocities `agents.x_d1`, the marks `agents.active`, and the observables `polarization`, `active_distance`,
and `passive_distance`.
Raises:
ValidationError:
The number of active agents is not between 1 and `n_agents`.
Examples:
>>> from skcomplex.systems.physics.active_matter import active_passive_flocking
>>> [element.label for element in active_passive_flocking(114).elements]
['agents']
"""
if not 1 <= n_active <= n_agents:
error_msg = f'The system has {n_active} active agents among {n_agents}. Give at least one, and at most all.'
raise ValidationError(error_msg)
rng = np.random.default_rng(random_state)
agents = GroupElement(
'agents',
element=FlockingAgentElement,
size=n_agents,
values={
'x': rng.uniform(0.0, size, (n_agents, 2)),
'x_d1': rng.uniform(-initial_speed, initial_speed, (n_agents, 2)),
'active': np.where(np.arange(n_agents) < n_active, ACTIVE_AGENT, PASSIVE_AGENT),
},
visual=VisualSpec(
'points',
{'position': 'x', 'radius': radius, 'color_by': 'active', 'color_range': (0.0, 1.0)},
),
)
energy = functools.partial(_repulsion_energy, epsilon=epsilon, radius=radius, exponent=exponent)
return System(
elements=[agents],
interactions=[
PairCouplingInteraction('repulsion', 'agents', 'agents', cutoff=2 * radius, potential=energy, period=size),
RateInteraction(
'propulsion',
target='agents',
rate=functools.partial(_propulsion, mu=mu, v0=v0),
noise=functools.partial(_random_force, amplitude=eta * math.sqrt(time_step / 3)),
),
RateInteraction(
'flocking',
target='agents',
rate=functools.partial(_flocking, alpha=alpha, reach=flocking_radius),
),
],
space=ContinuousSpace(extent=(size, size), boundary='periodic'),
)
Components
Elements:
FlockingAgentElement, fromPointElement
Interactions:
PairCouplingInteractionRateInteraction
class FlockingAgentElement(PointElement):
"""Agent of a crowd, a point with a mark that tells whether it is active or passive.
Args:
label:
The label of the agent or the group.
x:
The initial position.
x_d1:
The initial velocity.
visual:
How the agent is drawn, or `None` for `default_visual_`.
active:
The mark, 1 for an active agent or 0 for a passive one.
Examples:
>>> from skcomplex.systems.physics.active_matter import FlockingAgentElement
>>> FlockingAgentElement('agent', active=0.0).active
0.0
"""
def __init__(
self,
label: str,
x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2, role='position')] = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
*,
active: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([PASSIVE_AGENT, ACTIVE_AGENT]))] = ACTIVE_AGENT,
) -> None:
super().__init__(label, x=x, x_d1=x_d1, visual=visual)
self.active = active
class PairCouplingInteraction(PairInteraction):
"""Coupling of every pair of members closer than a cutoff, through a function of the pair or a pair potential.
Every member of a pair adds a value to the rate of its state `write`, which gives the highest derivative of a state
of any order. With `function`, the value is `function(namespace)`, where `namespace.first` is the member that adds
it and `namespace.second` the other member. With `potential`, a function `potential(namespace, r)` of the distance
`r` between the values of the vector state `position` of the pair, the value is minus the gradient of the potential
with respect to the member's own `position`. A `divisor` names a state or parameter of the members that divides
the value. With several states `write`, the function returns one value per state, in their order, so one pair
couples several states jointly.
Args:
label:
The label of the interaction, unique among the parts of its system.
first:
The path of the group whose members are the first of every pair.
second:
The path of the group whose members are the second of every pair, which may be the first group.
cutoff:
The largest distance between the members of a pair.
position:
The name of the vector state whose distances select the pairs, and which the potential is a function of.
write:
The name of the state whose rate the coupling adds to, or a sequence of names, one per value of the
function.
function:
The value a member adds, or a sequence of values, one per state `write`, a callable of the namespace of the
pair, or `None` when the potential is given.
potential:
The pair potential, a callable of the namespace of the pair and the distance, or `None` when the function
is given.
divisor:
The name of the state or parameter of the members that divides the value, or `None` for none.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
Raises:
ValidationError:
Not exactly one of the function and the potential is a callable, or the potential has several states to
write, raised by `prepare`.
Examples:
>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var, vector
>>> from skcomplex.core.components.interactions import PairCouplingInteraction
>>> first = SimpleNamespace(x=vector([var('x1')]), m=var('m1'))
>>> second = SimpleNamespace(x=vector([var('x2')]), m=var('m2'))
>>> spring = PairCouplingInteraction('spring', 'a', 'b', potential=lambda s, r: r**2 / 2, divisor='m')
>>> [law.value for law in spring.prepare().run(SimpleNamespace(first=first, second=second)).laws_]
[(-(x1 - x2)/m1), ((x1 - x2)/m2)]
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw,)
def __init__(
self,
label: str,
first: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('position', 'write', 'divisor'))],
second: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('position', 'write', 'divisor'))],
cutoff: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
position: str = 'x',
*,
write: str | Sequence[str] = 'x',
function: NamespaceFunction | None = None,
potential: PairPotential | None = None,
divisor: str | None = None,
period: float | None = None,
) -> None:
super().__init__(label, first, second, cutoff=cutoff, position=position)
self.write = write
self.function = function
self.potential = potential
self.divisor = divisor
self.period = period
def _check(self) -> None:
"""Raise an error unless exactly one of the function and the potential is a callable.
Raises:
ValidationError:
The function and the potential are both given or both missing, the given one is not callable, or the
potential has several states to write.
"""
super()._check()
given = [item for item in (self.function, self.potential) if item is not None]
if len(given) != 1 or not callable(given[0]):
error_msg = (
f'{type(self).__name__} {self.label!r} has the function {self.function!r} and the potential '
f'{self.potential!r}. Give exactly one of them as a callable.'
)
raise ValidationError(error_msg)
if self.potential is not None and not isinstance(self.write, str):
error_msg = f'{type(self).__name__} {self.label!r} writes {self.write!r}, but a potential writes one state.'
raise ValidationError(error_msg)
def _get_gradient_value(self, namespace: SimpleNamespace) -> Target:
"""Return minus the gradient of the pair potential with respect to the position of the first member."""
assert self.potential is not None
distance = SR.var('r', domain='positive')
first, second = self._get_member(namespace, 'first'), self._get_member(namespace, 'second')
separation = getattr(first, self.position) - getattr(second, self.position)
r = sqrt(separation.dot_product(separation))
slope = SR(self.potential(namespace, distance)).diff(distance).subs({distance: r})
return -slope * separation / r
def _get_values(self, namespace: SimpleNamespace) -> list[Target]:
"""Return the values the first member of a pair adds to the rates of its states, in the order of `write`."""
if self.function is None:
values = [self._get_gradient_value(namespace)]
elif isinstance(self.write, str):
values = [self.function(namespace)]
else:
values = list(cast(Sequence[Target], self.function(namespace)))
if self.divisor is None:
return values
return [value / getattr(self._get_member(namespace, 'first'), self.divisor) for value in values]
def _get_names(self) -> list[str]:
"""Return the names of the states the coupling writes."""
return [self.write] if isinstance(self.write, str) else list(self.write)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rate of the written state of both members of a pair.
Args:
namespace:
The namespace of the symbols the laws can use, with `first` and `second` for the pair.
Returns:
laws:
The rate laws of the written states of the two members.
"""
first, second = self._get_member(namespace, 'first'), self._get_member(namespace, 'second')
swapped = SimpleNamespace(
**{**vars(namespace), self._get_role('first'): second, self._get_role('second'): first},
)
names = self._get_names()
laws: list[BaseLaw] = []
for member, pair in ((first, namespace), (second, swapped)):
values = self._get_values(pair)
laws += [RateLaw(getattr(member, name), value) for name, value in zip(names, values, strict=True)]
return laws
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
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 rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.components.elements import CompartmentElement
>>> from skcomplex.core.components.interactions import RateInteraction
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw, NoiseLaw, FractionalNoiseLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'divisor'))],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws