Mirollo-Strogatz pulse-coupled oscillators
Integrate-and-fire oscillators coupled by pulses, which end up firing in unison.
- Name:
mirollo_strogatz - Field:
physics/synchronization - Formalism: ordinary differential equations
- References: Mirollo and Strogatz 1990; Arenas et al. 2008
Factory
@register_model(
field='physics/synchronization',
tags=('synchronization', 'events', 'pulse coupling', 'integrate and fire'),
references=('Mirollo and Strogatz 1990', 'Arenas et al. 2008'),
defaults={'until': 20.0, 'dt': 0.1},
)
def mirollo_strogatz(
n_oscillators: int = 5,
*,
S: float = 2.0,
gamma: float = 1.0,
epsilon: float = 0.1,
seed: int = 0,
) -> System:
"""Return integrate-and-fire oscillators coupled by pulses, which end up firing in unison.
Every oscillator of the group charges as `x' = S - gamma x`, a concave rise, and fires when it reaches 1. Firing is
an event whose guard is the group maximum of the charges, and it resets the most charged oscillators to 0 and lifts
every other one by `epsilon`. An oscillator lifted to 1 fires at the same instant and is absorbed into the group
that fired, so the groups merge until all oscillators fire together.
Args:
n_oscillators:
The number of oscillators.
S:
The charging input.
gamma:
The leak rate, below `S` so that the oscillators fire.
epsilon:
The strength of a pulse.
seed:
The seed of the initial charges.
Returns:
system:
The model, with the charges `oscillators.x` of the group members and the firings as `fire` events.
Examples:
>>> from skcomplex.models.physics.synchronization import mirollo_strogatz
>>> mirollo_strogatz(3).elements[0].check().prepare().size_
3
"""
charges = np.random.default_rng(seed).uniform(0.0, 1.0, n_oscillators)
oscillators = GroupElement(
'oscillators',
element=IntegrateAndFireOscillatorElement,
size=n_oscillators,
values={'x': charges, 'S': S, 'gamma': gamma, 'epsilon': epsilon},
)
fire = CrossingInteraction(
'fire',
target='oscillators',
guard=lambda s: maximum(s.target.x) - 1,
reset={
'x': lambda s: where(
s.target.x >= maximum(s.target.x),
0,
where(s.target.x + s.target.epsilon >= 1, 0, s.target.x + s.target.epsilon),
),
},
direction=1,
)
return System(elements=[oscillators], interactions=[fire])
Components
class IntegrateAndFireOscillatorElement(Element):
"""Integrate-and-fire oscillator that charges as `x' = S - gamma x`, and sends pulses of strength `epsilon`.
The oscillator charges towards `S / gamma`. Its firing at 1 is a crossing of the group it belongs to, which resets
the most charged oscillators to 0 and lifts every other one by `epsilon`.
Args:
label:
The label of the oscillator or the group.
visual:
How the oscillator is drawn, or `None` for `default_visual_`.
x:
The initial charge.
S:
The charging input.
gamma:
The leak rate.
epsilon:
The strength of a pulse.
Attributes:
default_visual_:
The visual of the oscillator when it sets none.
Examples:
>>> from skcomplex.models.physics.synchronization import IntegrateAndFireOscillatorElement
>>> IntegrateAndFireOscillatorElement('oscillator').epsilon
0.1
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec()] = 0.0,
S: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.S = S
self.gamma = gamma
self.epsilon = epsilon
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the charging of the oscillator.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rate of the charge.
"""
s = namespace
return [RateLaw(s.x, s.S - s.gamma * s.x)]
class CrossingInteraction(Interaction):
"""Event at the crossing of zero by a guard, which sets states of the target to new values or ends the run.
The guard and the new values are functions of the namespace, which holds `namespace.target`, the states and
parameters of the target, and the time `namespace.t`. A guard that reads the members of a group gives every member
its own event, which sets the states of that member alone. Solvers locate the crossing within their tolerance, and
the trajectory records every crossing among its events, under the label of the interaction and with the member.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the guard reads and the event sets.
guard:
The function `guard(namespace)` of the expression whose crossing of zero triggers the event.
reset:
A mapping from the name of a state of the target to the function `f(namespace)` of its new value, applied
all at once, or `None` to keep the states.
direction:
`1` for crossings from negative to positive, `-1` for the opposite, and `0` for both.
terminal:
Whether the crossing ends the run.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import CrossingInteraction
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class ClockElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.x = x
...
... def _get_laws(self, s):
... return [RateLaw(s.x, 1.0)]
>>> wrap = CrossingInteraction(
... 'wrap', target='clock', guard=lambda s: s.target.x - 1, reset={'x': lambda s: 0}, direction=1
... )
>>> system = System(elements=[ClockElement('clock')], interactions=[wrap])
>>> trajectory = Simulation(system, dt=0.5).check().prepare().run(until=2.2).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
guard: Callable[[SimpleNamespace], object],
reset: Mapping[str, Callable[[SimpleNamespace], object]] | None = None,
direction: int = 0,
terminal: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.guard = guard
self.reset = reset
self.direction = direction
self.terminal = terminal
def _check(self) -> None:
"""Raise an error for a direction that is not -1, 0, or 1, or a guard or a new value that is not a function.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
"""
super()._check()
if self.direction not in DIRECTIONS:
error_msg = f'The direction of crossing {self.label!r} is {self.direction!r}. Use -1, 0, or 1.'
raise ValidationError(error_msg)
if not callable(self.guard) or not all(callable(function) for function in (self.reset or {}).values()):
error_msg = f'The guard and the new values of crossing {self.label!r} must be functions of the namespace.'
raise ValidationError(error_msg)
@staticmethod
def _get_law_value(value: object) -> object:
"""Return what a function of the namespace returned: a function as it is, or a symbolic expression."""
return value if BaseLaw._is_function(value) else SR(value)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing law of the guard with the new values of the states.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The crossing law.
"""
reset = {
getattr(self._get_member(namespace, 'target'), name): self._get_law_value(f(namespace))
for name, f in (self.reset or {}).items()
}
return [
CrossingLaw(
self._get_law_value(self.guard(namespace)),
reset or None,
self.direction,
self.terminal,
label=self.label,
),
]
Equations
oscillators.x[i]: diff(x_oscillators(t), t) == -gamma_oscillators*x_oscillators(t) + S_oscillators