Synchronization

Synchronization: coupled phase oscillators, pulse-coupled oscillators, swarmalators, chaos, and Kuramoto order.

AdaptiveCouplingInteraction(label, *, nodes, links, K=1.0, alpha=0.0, beta=0.0, epsilon=0.01)

Bases: NeighborInteraction

Adaptive Kuramoto-Sakaguchi coupling, whose weights relax towards a phase-lagged sine of the phase difference.

On every directed link from j to i, the phase of i turns by -K kappa_ij sin(theta_i - theta_j + alpha) and the weight relaxes as kappa_ij' = -epsilon (kappa_ij + sin(theta_i - theta_j + beta)), so it grows between oscillators whose phases differ by about -beta - pi / 2.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
nodes	Annotated[str, RoleSpec()]	The path of the group of phase oscillators on the nodes.	required
links	Annotated[str, RoleSpec()]	The path of the group of adaptive links on the edges.	required
K	Annotated[float, ParameterSpec()]	The coupling strength, 1 / N for the mean coupling of N oscillators.	1.0
alpha	Annotated[float, ParameterSpec()]	The phase lag of the coupling.	0.0
beta	Annotated[float, ParameterSpec()]	The phase lag of the adaptation.	0.0
epsilon	Annotated[float, ParameterSpec()]	The rate of the adaptation.	0.01

Examples:

>>> from skcomplex.systems.physics.synchronization import AdaptiveCouplingInteraction
>>> AdaptiveCouplingInteraction('adaptation', nodes='nodes', links='links').epsilon
0.01

Source code in src/skcomplex/systems/physics/synchronization/_adaptive.py

def __init__(
    self,
    label: str,
    *,
    nodes: Annotated[str, RoleSpec()],
    links: Annotated[str, RoleSpec()],
    K: Annotated[float, ParameterSpec()] = 1.0,
    alpha: Annotated[float, ParameterSpec()] = 0.0,
    beta: Annotated[float, ParameterSpec()] = 0.0,
    epsilon: Annotated[float, ParameterSpec()] = 0.01,
) -> None:
    super().__init__(label)
    self.nodes = nodes
    self.links = links
    self.K = K
    self.alpha = alpha
    self.beta = beta
    self.epsilon = epsilon

AdaptiveLinkElement(label, visual=None, *, kappa=0.0)

Bases: Element

Directed link of an adaptive network with the coupling weight kappa, a state that the phases of its ends adapt.

Parameters:

Name	Type	Description	Default
label	str	The label of the link or the group on the links.	required
visual	VisualSpec | None	How the link is drawn, or None for default_visual_.	None
kappa	Annotated[float, StateSpec()]	The initial coupling weight.	0.0

Attributes:

Name	Type	Description
default_visual_		The visual of the link when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import AdaptiveLinkElement
>>> AdaptiveLinkElement('link', kappa=0.5).kappa
0.5

Source code in src/skcomplex/systems/physics/synchronization/_adaptive.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    kappa: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.kappa = kappa

AggregateAveragingInteraction(label, *, source, target, r=4.0, kappa=0.4, kind='logistic', delay=0)

Bases: MapInteraction

Logistic or Bernoulli map of an aggregate of the value of a target and the value of a source, possibly delayed.

The map sets x_t(t + 1) = f(a) with the aggregate a = (1 - kappa) x_t(t) + kappa x_s(t - delay), and with the logistic map f(a) = r a (1 - a) or the Bernoulli shift f(a) = r a mod 1. A delayed source keeps delay lags.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
source	Annotated[str, RoleSpec()]	The path of the unit whose value the aggregate reads.	required
target	Annotated[str, RoleSpec()]	The path of the unit whose value the map sets.	required
r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The growth rate of the logistic map.	4.0
kappa	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The weight of the source in the aggregate.	0.4
kind	str	The map, 'logistic' or 'bernoulli'.	'logistic'
delay	int	The number of steps by which the source is delayed.	0

Raises:

Type	Description
ValidationError	The map is unknown, or the delay is not a nonnegative integer.

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.systems.physics.synchronization import AggregateAveragingInteraction
>>> averaging = AggregateAveragingInteraction('averaging', source='a', target='b').prepare()
>>> source, target = SimpleNamespace(x=var('u')), SimpleNamespace(x=var('v'))
>>> namespace = SimpleNamespace(source=source, target=target, r=4, kappa=var('kappa'))
>>> [law.value.factor() for law in averaging.run(namespace).laws_]
[-4*(kappa*u - kappa*v + v)*(kappa*u - kappa*v + v - 1)]

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

def __init__(
    self,
    label: str,
    *,
    source: Annotated[str, RoleSpec()],
    target: Annotated[str, RoleSpec()],
    r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.4,
    kind: str = 'logistic',
    delay: int = 0,
) -> None:
    super().__init__(label, target=target, state='x', function=self._average)
    self.source = source
    self.r = r
    self.kappa = kappa
    self.kind = kind
    self.delay = delay

DelayedCouplingInteraction(label, *, source, target, K=1.0, tau=1.0)

Bases: RateInteraction

Delayed sine coupling of one phase oscillator to another, theta_t' += K sin(theta_s(t - tau) - theta_t).

The target reads the phase of the source as it was a delay tau earlier.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
source	Annotated[str, RoleSpec()]	The path of the oscillator whose delayed phase the coupling reads.	required
target	Annotated[str, RoleSpec()]	The path of the oscillator whose phase the coupling turns.	required
K	Annotated[float, ParameterSpec()]	The coupling strength.	1.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The delay.	1.0

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.systems.physics.synchronization import DelayedCouplingInteraction
>>> coupling = DelayedCouplingInteraction('coupling', source='a', target='b').prepare()
>>> namespace = SimpleNamespace(source=SimpleNamespace(theta=var('a')), target=SimpleNamespace(theta=var('b')))
>>> [law.value for law in coupling.run(SimpleNamespace(**vars(namespace), K=var('K'), tau=1)).laws_]
[K*sin(-b + delayed(a, 1))]

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

def __init__(
    self,
    label: str,
    *,
    source: Annotated[str, RoleSpec()],
    target: Annotated[str, RoleSpec()],
    K: Annotated[float, ParameterSpec()] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, target=target, state='theta')
    self.source = source
    self.K = K
    self.tau = tau

DelayedUnitElement(label, visual=None, *, value=0.3, delay=1, variables=None)

Bases: DeclaredElement

Value x of a map that keeps its last delay values, which a delayed coupling reads.

The state x has delay lags, x_lag1 to x_lag<delay>, and before the start every lag holds the initial value.

Parameters:

Name	Type	Description	Default
label	str	The label of the unit.	required
visual	VisualSpec | None	How the unit is drawn, or None for default_visual_.	None
value	float	The initial value of x, also its value before the start.	0.3
delay	int	The number of past values the unit keeps, at least 1.	1
variables	Mapping[str, DeclaredVariableSpec] | None	The declared variables that replace the derived ones, or None for none.	None

Attributes:

Name	Type	Description
default_visual_		The visual of the unit when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import DelayedUnitElement
>>> DelayedUnitElement('unit', delay=2).get_value('x')
0.3

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    value: float = 0.3,
    delay: int = 1,
    variables: Mapping[str, DeclaredVariableSpec] | None = None,
) -> None:
    super().__init__(label, visual=visual, variables=variables)
    self.value = value
    self.delay = delay

DiffusiveCouplingInteraction(label, *, source, target, state='x2', k=0.1)

Bases: RateInteraction

Diffusive coupling of a state of one element to the same state of another, x_t' += k (x_s - x_t).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
source	Annotated[str, RoleSpec()]	The path of the element whose state the coupling reads.	required
target	Annotated[str, RoleSpec()]	The path of the element whose state the coupling changes.	required
state	str	The name of the coupled state.	'x2'
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The coupling strength.	0.1

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.systems.physics.synchronization import DiffusiveCouplingInteraction
>>> coupling = DiffusiveCouplingInteraction('coupling', source='a', target='b').prepare()
>>> source, target = SimpleNamespace(x2=var('u')), SimpleNamespace(x2=var('v'))
>>> [law.value for law in coupling.run(SimpleNamespace(source=source, target=target, k=var('k'))).laws_]
[k*(u - v)]

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

def __init__(
    self,
    label: str,
    *,
    source: Annotated[str, RoleSpec()],
    target: Annotated[str, RoleSpec()],
    state: str = 'x2',
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
    super().__init__(label, target=target, state=state)
    self.source = source
    self.k = k

EdgeCouplingInteraction(label, *, target, links, K=1.0, alpha=0.0)

Bases: KuramotoCouplingInteraction

Kuramoto coupling along the weighted edges of a network, theta_i' += K w_ij sin(theta_j - theta_i).

The coupling runs once per edge of the network of its oscillators and reads the weight w of the edge from the group on the edges.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group of oscillators on the nodes.	required
links	Annotated[str, RoleSpec()]	The path of the group on the edges, whose parameter w weighs every edge.	required
K	float	The coupling strength.	1.0
alpha	float	The phase lag.	0.0

Examples:

>>> from skcomplex.systems.physics.synchronization import EdgeCouplingInteraction
>>> EdgeCouplingInteraction('coupling', target='nodes', links='edges', K=0.5).weights
'w'

Source code in src/skcomplex/systems/physics/synchronization/_networks.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    links: Annotated[str, RoleSpec()],
    K: float = 1.0,
    alpha: float = 0.0,
) -> None:
    super().__init__(label, target=target, K=K, alpha=alpha, weights='w')
    self.links = links

GlobalCouplingInteraction(label, *, target, state='x', K=1.0)

Bases: CouplingInteraction

Diffusive coupling of one state of every member to the mean of the group, x_i' += K (mean_j x_j - x_i).

The mean runs over every member, the member itself included.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group.	required
state	str	The name of the coupled state.	'x'
K	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The coupling strength.	1.0

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.systems.physics.synchronization import GlobalCouplingInteraction
>>> coupling = GlobalCouplingInteraction('coupling', target='oscillators', state='y').prepare()
>>> member, other = SimpleNamespace(y=var('y')), SimpleNamespace(y=var('z'))
>>> coupling.function(SimpleNamespace(K=var('K'), member=member, other=other))
-K*(y - z)

Source code in src/skcomplex/systems/physics/synchronization/_amplitude_death.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    state: str = 'x',
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(
        label,
        target=target,
        read=state,
        write=state,
        function=self._pull,
        reduce='mean',
        mode='group',
        include_self=True,
    )
    self.state = state
    self.K = K

IntegrateAndFireOscillatorElement(label, visual=None, *, x=0.0, S=2.0, gamma=1.0, epsilon=0.1)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator or the group.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial charge.	0.0
S	Annotated[float, ParameterSpec(domain=open(0, oo))]	The charging input.	2.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The leak rate.	1.0
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of a pulse.	0.1

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import IntegrateAndFireOscillatorElement
>>> IntegrateAndFireOscillatorElement('oscillator').epsilon
0.1

Source code in src/skcomplex/systems/physics/synchronization/_pulses.py

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

LinkElement(label, visual=None, *, w=1.0)

Bases: Element

Edge of a network with a weight w, which couplings along the edge read.

Parameters:

Name	Type	Description	Default
label	str	The label of the edge or the group on the edges.	required
visual	VisualSpec | None	How the edge is drawn, or None for default_visual_.	None
w	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The weight of the edge.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the edge when it sets none.

Examples:

>>> from skcomplex.systems.physics.common.oscillators import LinkElement
>>> LinkElement('link', w=0.5).w
0.5

Source code in src/skcomplex/systems/physics/common/oscillators/_link.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    w: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.w = w

LogisticCouplingInteraction(label, *, target, state='x')

Bases: CouplingInteraction

Coupled logistic maps, x_i(t + 1) = (1 - epsilon) f(x_i) + epsilon mean_j f(x_j) with f(x) = r x (1 - x).

The map of every node mixes with the mean of the maps of its neighbours j on its lattice or network, with the parameters r and epsilon of the node. A node with no neighbours keeps (1 - epsilon) f(x_i). It is the map form of coupling, with the mean of the maps of the neighbours combined with the map of the node.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group of maps, on the nodes of a lattice or network.	required
state	str	The name of the state that the maps iterate.	'x'

Examples:

>>> import networkx as nx
>>> from skcomplex.core.components.elements import LatticeElement
>>> from skcomplex.core.components.spaces import DiscreteTime, NetworkSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.physics.synchronization import LogisticCouplingInteraction, LogisticNodeElement
>>> maps = LatticeElement(
...     'maps',
...     element=LogisticNodeElement,
...     space=NetworkSpace(nx.path_graph(2)),
...     values={'x': [0.25, 0.5], 'epsilon': 0.5},
... )
>>> coupling = LogisticCouplingInteraction('coupling', target='maps')
>>> system = System(elements=[maps], interactions=[coupling], time=DiscreteTime())
>>> Simulation(system).prepare().run(until=1).trajectory_.values('maps.x')[-1].tolist()
[0.875, 0.875]

Source code in src/skcomplex/systems/physics/synchronization/_networks.py

def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], state: str = 'x') -> None:
    super().__init__(
        label,
        target=target,
        read=state,
        write=state,
        function=self._get_map,
        reduce='mean',
        form='map',
        combine=self._get_mix,
    )
    self.state = state

LogisticNodeElement(label, visual=None, *, x=0.5, r=4.0, epsilon=0.5)

Bases: Element

Logistic map on a node of a network, with its value, growth rate, and weight of the neighbours.

LogisticCouplingInteraction iterates x_i(t + 1) = (1 - epsilon) f(x_i) + epsilon mean_j f(x_j), where the map f(x) = r x (1 - x) of the node mixes with the mean of the maps of its neighbours j on the network.

Parameters:

Name	Type	Description	Default
label	str	The label of the node or the group on the nodes.	required
visual	VisualSpec | None	How the node is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial value.	0.5
r	Annotated[float, ParameterSpec()]	The growth rate of the logistic map.	4.0
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The weight of the neighbours.	0.5

Attributes:

Name	Type	Description
default_visual_		The visual of the node when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import LogisticNodeElement
>>> LogisticNodeElement('node').epsilon
0.5

Source code in src/skcomplex/systems/physics/synchronization/_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.5,
    r: Annotated[float, ParameterSpec()] = 4.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.r = r
    self.epsilon = epsilon

LogisticUnitElement(label, visual=None, *, x=0.3)

Bases: Element

Value x in [0, 1] of a logistic map, which maps of the coupling iterate.

Parameters:

Name	Type	Description	Default
label	str	The label of the unit.	required
visual	VisualSpec | None	How the unit is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial value.	0.3

Attributes:

Name	Type	Description
default_visual_		The visual of the unit when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import LogisticUnitElement
>>> LogisticUnitElement('unit').x
0.3

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x

RosslerOscillatorElement(label, visual=None, *, x1=1.0, x2=0.0, x3=0.0, omega=1.0, a=0.2, b=0.2, c=5.7, p=1.0)

Bases: Element

Rössler oscillator x1' = -omega x2 - x3, x2' = omega x1 + a x2, x3' = b + x3 (x1 - c).

It is chaotic for a = b = 0.2 and c = 5.7. The time scale p divides every rate, p x' = f(x), so p < 1 speeds the oscillator up.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x1	Annotated[float, StateSpec()]	The initial first variable.	1.0
x2	Annotated[float, StateSpec()]	The initial second variable.	0.0
x3	Annotated[float, StateSpec()]	The initial third variable.	0.0
omega	Annotated[float, ParameterSpec()]	The natural frequency.	1.0
a	Annotated[float, ParameterSpec()]	The first Rössler parameter.	0.2
b	Annotated[float, ParameterSpec()]	The second Rössler parameter.	0.2
c	Annotated[float, ParameterSpec()]	The third Rössler parameter.	5.7
p	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time scale, which divides every rate.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.systems.physics.synchronization import RosslerOscillatorElement
>>> RosslerOscillatorElement('oscillator').c
5.7

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x1: Annotated[float, StateSpec()] = 1.0,
    x2: Annotated[float, StateSpec()] = 0.0,
    x3: Annotated[float, StateSpec()] = 0.0,
    omega: Annotated[float, ParameterSpec()] = 1.0,
    a: Annotated[float, ParameterSpec()] = 0.2,
    b: Annotated[float, ParameterSpec()] = 0.2,
    c: Annotated[float, ParameterSpec()] = 5.7,
    p: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x1 = x1
    self.x2 = x2
    self.x3 = x3
    self.omega = omega
    self.a = a
    self.b = b
    self.c = c
    self.p = p

SimplexCouplingInteraction(label, *, target, K1=0.0, K2=1.0)

Bases: HyperedgeInteraction

Kuramoto coupling through the simplices of a hypergraph, pairs and triangles.

On a link, the focal oscillator i gets K1 sin(theta_j - theta_i), and on a triangle K2 sin(theta_j + theta_k - 2 theta_i) from its two other members j and k.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group of oscillators on the nodes.	required
K1	Annotated[float, ParameterSpec()]	The coupling of every link.	0.0
K2	Annotated[float, ParameterSpec()]	The coupling of every triangle.	1.0

Examples:

>>> from skcomplex.systems.physics.synchronization import SimplexCouplingInteraction
>>> SimplexCouplingInteraction('coupling', target='oscillators', K2=0.5).K2
0.5

Source code in src/skcomplex/systems/physics/synchronization/_hypergraphs.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    K1: Annotated[float, ParameterSpec()] = 0.0,
    K2: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.K1 = K1
    self.K2 = K2

SwarmalatorCouplingInteraction(label, first, second, cutoff=1000000.0, position='x', *, A=1.0, B=1.0, J=1.0, K=0.0, period=None, n_units=1)

Bases: PairCouplingInteraction

Coupling of two swarmalators through their separation r and their phase difference dtheta.

Each swarmalator moves toward the other at the speed (A + J cos(dtheta)) / N, away from it at B / (N |r|), and turns its phase toward the other's at the rate K sin(dtheta) / (N |r|), as in O'Keeffe, Ceron, and Petersen. It is a pair coupling of the positions and the phases jointly.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
first	Annotated[str, RoleSpec()]	The path of the group of the first members.	required
second	Annotated[str, RoleSpec()]	The path of the group of the second members.	required
cutoff	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance beyond which pairs do not interact, large for the all-to-all system.	1000000.0
position	str	The name of the vector state that holds the positions of the members.	'x'
A	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the attraction.	1.0
B	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the repulsion.	1.0
J	Annotated[float, ParameterSpec()]	The influence of the phase similarity on the attraction.	1.0
K	Annotated[float, ParameterSpec()]	The phase coupling, synchronizing when positive.	0.0
period	float | None	The side of a periodic box, for minimum image distances, or None for an open space.	None
n_units	int	The number of swarmalators N that normalizes the coupling.	1

Examples:

>>> from skcomplex.systems.physics.synchronization import SwarmalatorCouplingInteraction
>>> SwarmalatorCouplingInteraction('coupling', first='units', second='units', J=0.5).J
0.5

Source code in src/skcomplex/systems/physics/synchronization/_swarmalators.py

def __init__(
    self,
    label: str,
    first: Annotated[str, RoleSpec()],
    second: Annotated[str, RoleSpec()],
    cutoff: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1000000.0,
    position: str = 'x',
    *,
    A: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    B: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    J: Annotated[float, ParameterSpec()] = 1.0,
    K: Annotated[float, ParameterSpec()] = 0.0,
    period: float | None = None,
    n_units: int = 1,
) -> None:
    super().__init__(
        label,
        first,
        second,
        cutoff=cutoff,
        position=position,
        write=(position, 'theta'),
        function=self._get_terms,
        period=period,
    )
    self.A = A
    self.B = B
    self.J = J
    self.K = K
    self.n_units = n_units

SwarmalatorElement(label, visual=None, *, theta=0.0, omega=1.0, x=0.0)

Bases: SitedOscillatorElement

Phase oscillator whose place x is a position that moves by the rates of interactions.

Parameters:

Name	Type	Description	Default
label	str	The label of the swarmalator or the group.	required
visual	VisualSpec | None	How the swarmalator is drawn, or None for default_visual_.	None
theta	float	The initial phase.	0.0
omega	float	The natural frequency.	1.0
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position. A number applies to every component.	0.0

Examples:

>>> from skcomplex.systems.physics.synchronization import SwarmalatorElement
>>> SwarmalatorElement('unit', x=[0.0, 1.0]).x
[0.0, 1.0]

Source code in src/skcomplex/systems/physics/synchronization/_swarmalators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    theta: float = 0.0,
    omega: float = 1.0,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
) -> None:
    super().__init__(label, visual=visual, theta=theta, omega=omega, x=x)

WinfreeMeanFieldInteraction(label, *, target, K=1.0)

Bases: CouplingInteraction

Winfree coupling theta_i' += K mean_j P(theta_j) R(theta_i) over every member j of the group.

The influence P(theta) = 1 + cos(theta) that an oscillator sends and the sensitivity R(theta) = -sin(theta) with which it responds make every oscillator respond to the mean influence of the group.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group of oscillators.	required
K	Annotated[float, ParameterSpec()]	The coupling strength.	1.0

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.systems.physics.synchronization import WinfreeMeanFieldInteraction
>>> coupling = WinfreeMeanFieldInteraction('coupling', target='oscillators')
>>> member, other = SimpleNamespace(theta=var('a')), SimpleNamespace(theta=var('b'))
>>> coupling.function(SimpleNamespace(K=var('K'), member=member, other=other))
-K*(cos(b) + 1)*sin(a)

Source code in src/skcomplex/systems/physics/synchronization/_pulses.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    K: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(
        label,
        target=target,
        read='theta',
        write='theta',
        function=self._response,
        reduce='mean',
        mode='group',
        include_self=True,
    )
    self.K = K

adaptive_kuramoto(n_oscillators=50, *, alpha=0.3 * math.pi, beta=0.23 * math.pi, epsilon=0.01, omega=0.0, seed=0)

Return globally coupled phase oscillators whose coupling weights adapt to their phase differences.

The oscillators follow the Kuramoto-Sakaguchi model theta_i' = omega - (1 / N) Σ_j kappa_ij sin(theta_i - theta_j + alpha), and every weight of the directed link from j to i adapts as kappa_ij' = -epsilon (kappa_ij + sin(theta_i - theta_j + beta)), Eqs. 28-29 of Berner et al. The weights are states of a group on the edges of the complete directed graph. Even identical oscillators split into frequency clusters: groups whose members share one mean frequency, strongly coupled within and weakly coupled across, with the weights inside a cluster near -sin(theta_i - theta_j + beta). The phases and the weights start uniformly at random in [0, 2 pi) and [-1, 1].

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	50
alpha	float	The phase lag of the coupling.	0.3 * pi
beta	float	The phase lag of the adaptation.	0.23 * pi
epsilon	float	The rate of the adaptation, slow for small values.	0.01
omega	float	The natural frequency of every oscillator.	0.0
seed	int	The seed of the random initial phases and weights.	0

Returns:

Name	Type	Description
system	System	The system, with the phases nodes.theta and the weights links.kappa, one per directed link.

Examples:

>>> from skcomplex.systems.physics.synchronization import adaptive_kuramoto
>>> adaptive_kuramoto(5).elements[1].prepare().size_
20

Source code in src/skcomplex/systems/physics/synchronization/_adaptive.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'adaptive network', 'frequency clusters', 'phase', 'ode'),
    references=('Berner et al. 2023', 'Berner et al. 2019', 'Kasatkin et al. 2017'),
    defaults={'until': 2000.0, 'dt': 1.0},
)
def adaptive_kuramoto(
    n_oscillators: int = 50,
    *,
    alpha: float = 0.3 * math.pi,
    beta: float = 0.23 * math.pi,
    epsilon: float = 0.01,
    omega: float = 0.0,
    seed: int = 0,
) -> System:
    """Return globally coupled phase oscillators whose coupling weights adapt to their phase differences.

    The oscillators follow the Kuramoto-Sakaguchi model `theta_i' = omega - (1 / N) Σ_j kappa_ij sin(theta_i - theta_j +
    alpha)`, and every weight of the directed link from `j` to `i` adapts as
    `kappa_ij' = -epsilon (kappa_ij + sin(theta_i - theta_j + beta))`, Eqs. 28-29 of Berner et al. The weights are
    states of a group on the edges of the complete directed graph. Even identical oscillators split into frequency
    clusters: groups whose members share one mean frequency, strongly coupled within and weakly coupled across, with
    the weights inside a cluster near `-sin(theta_i - theta_j + beta)`. The phases and the weights start uniformly at
    random in `[0, 2 pi)` and `[-1, 1]`.

    Args:
        n_oscillators:
            The number of oscillators.

        alpha:
            The phase lag of the coupling.

        beta:
            The phase lag of the adaptation.

        epsilon:
            The rate of the adaptation, slow for small values.

        omega:
            The natural frequency of every oscillator.

        seed:
            The seed of the random initial phases and weights.

    Returns:
        system:
            The system, with the phases `nodes.theta` and the weights `links.kappa`, one per directed link.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import adaptive_kuramoto
        >>> adaptive_kuramoto(5).elements[1].prepare().size_
        20
    """
    rng = np.random.default_rng(seed)
    graph = nx.complete_graph(n_oscillators, create_using=nx.DiGraph)
    space = NetworkSpace(graph)
    nodes = GroupElement(
        'nodes',
        element=PhaseOscillatorElement,
        space=space,
        on='sites',
        values={'theta': rng.uniform(0.0, 2 * math.pi, n_oscillators), 'omega': omega},
        visual=VisualSpec('sites', {'states': ['theta'], 'low': 0.0, 'high': 2 * math.pi}),
    )
    links = GroupElement(
        'links',
        element=AdaptiveLinkElement,
        space=space,
        on='edges',
        values={'kappa': rng.uniform(-1.0, 1.0, graph.number_of_edges())},
    )
    adaptation = AdaptiveCouplingInteraction(
        'adaptation',
        nodes='nodes',
        links='links',
        K=1 / n_oscillators,
        alpha=alpha,
        beta=beta,
        epsilon=epsilon,
    )
    return System(elements=[nodes, links], interactions=[adaptation])

aggregate_averaging_synchronization(x1=0.3, x2=0.6, *, r=4.0, kappa=0.4, kind='logistic', delay=0)

Return two chaotic maps coupled by aggregate averaging, which synchronize their chaos for strong coupling.

Each map iterates on the aggregate of its own value and the other's, x_i(t + 1) = f(a_i) with a_1 = (1 - kappa) x_1(t) + kappa x_2(t - T) and a_2 = (1 - kappa) x_2(t) + kappa x_1(t - T), for the delay T. The logistic map is f(a) = r a (1 - a). Without delay, its synchronized orbit is stable for kappa > (r - 1) / (2 r), which is 3 / 8 for r = 4, and the maps agree in one step for kappa = 1 / 2. The Bernoulli shift f(a) = r a mod 1 has the constant slope r, so a difference d of the maps follows d(t + 1) = r ((1 - kappa) d(t) - kappa d(t - T)) exactly, and the maps synchronize when every root of z^(T + 1) - r (1 - kappa) z^T + r kappa lies in the unit disk, which for T = 0 is |r (1 - 2 kappa)| < 1. The delay shrinks the window of synchronizing couplings: for r = 1.5 it is 0.17 < kappa < 0.83 without delay, 0.17 < kappa < 0.66 for T = 1, and 0.17 < kappa < 0.44 for T = 2, and it closes for T >= 5. Before the start, the delayed values equal the initial values. An integer slope such as 2 loses one binary digit of the value per step in floating point and falls onto 0, so a slope that is not an integer, such as 1.5, keeps the chaos.

Parameters:

Name	Type	Description	Default
x1	float	The initial value of the first map.	0.3
x2	float	The initial value of the second map.	0.6
r	float	The growth rate of the logistic map.	4.0
kappa	float	The coupling strength.	0.4
kind	str	The map, 'logistic' or 'bernoulli'.	'logistic'
delay	int	The delay T of the coupling, in steps.	0

Returns:

Name	Type	Description
system	System	The system, with the values first.x and second.x, and their lags first.x_lag<j> and second.x_lag<j> for a delay.

Examples:

>>> from skcomplex.systems.physics.synchronization import aggregate_averaging_synchronization
>>> aggregate_averaging_synchronization(kappa=0.1).interactions[0].kappa
0.1

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'chaos', 'map', 'logistic map', 'bernoulli map', 'delay'),
    references=('Gros 2015',),
    defaults={'until': 50, 'dt': 1},
)
def aggregate_averaging_synchronization(
    x1: float = 0.3,
    x2: float = 0.6,
    *,
    r: float = 4.0,
    kappa: float = 0.4,
    kind: str = 'logistic',
    delay: int = 0,
) -> System:
    """Return two chaotic maps coupled by aggregate averaging, which synchronize their chaos for strong coupling.

    Each map iterates on the aggregate of its own value and the other's, `x_i(t + 1) = f(a_i)` with `a_1 = (1 - kappa)
    x_1(t) + kappa x_2(t - T)` and `a_2 = (1 - kappa) x_2(t) + kappa x_1(t - T)`, for the delay `T`. The logistic map is
    `f(a) = r a (1 - a)`. Without delay, its synchronized orbit is stable for `kappa > (r - 1) / (2 r)`, which is `3 /
    8` for `r = 4`, and the maps agree in one step for `kappa = 1 / 2`. The Bernoulli shift `f(a) = r a mod 1` has the
    constant slope `r`, so a difference `d` of the maps follows `d(t + 1) = r ((1 - kappa) d(t) - kappa d(t - T))`
    exactly, and the maps synchronize when every root of `z^(T + 1) - r (1 - kappa) z^T + r kappa` lies in the unit
    disk, which for `T = 0` is `|r (1 - 2 kappa)| < 1`. The delay shrinks the window of synchronizing couplings: for `r
    = 1.5` it is `0.17 < kappa < 0.83` without delay, `0.17 < kappa < 0.66` for `T = 1`, and `0.17 < kappa < 0.44` for
    `T = 2`, and it closes for `T >= 5`. Before the start, the delayed values equal the initial values. An integer slope
    such as 2 loses one binary digit of the value per step in floating point and falls onto 0, so a slope that is not an
    integer, such as 1.5, keeps the chaos.

    Args:
        x1:
            The initial value of the first map.

        x2:
            The initial value of the second map.

        r:
            The growth rate of the logistic map.

        kappa:
            The coupling strength.

        kind:
            The map, `'logistic'` or `'bernoulli'`.

        delay:
            The delay `T` of the coupling, in steps.

    Returns:
        system:
            The system, with the values `first.x` and `second.x`, and their lags `first.x_lag<j>` and
            `second.x_lag<j>` for a delay.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import aggregate_averaging_synchronization
        >>> aggregate_averaging_synchronization(kappa=0.1).interactions[0].kappa
        0.1
    """
    if delay:
        units: list[Element] = [
            DelayedUnitElement('first', value=x1, delay=delay),
            DelayedUnitElement('second', value=x2, delay=delay),
        ]
    else:
        units = [LogisticUnitElement('first', x=x1), LogisticUnitElement('second', x=x2)]
    averaging = [
        AggregateAveragingInteraction(
            f'averaging_{target}',
            source=source,
            target=target,
            r=r,
            kappa=kappa,
            kind=kind,
            delay=delay,
        )
        for source, target in (('second', 'first'), ('first', 'second'))
    ]
    return System(elements=units, interactions=averaging, time=DiscreteTime())

coupled_map_lattices_on_networks(n_nodes=100, *, network='erdos_renyi', degree=20, rewiring=0.1, r=4.0, epsilon=0.9, seed=0)

Return chaotic logistic maps on the nodes of a network, each pulled towards the mean of its neighbours.

Node i follows x_i(t + 1) = (1 - epsilon) f(x_i) + (epsilon / k_i) sum_j A_ij f(x_j) with the logistic map f(x) = r x (1 - x), a mean over the neighbours of the network. For a coupling epsilon large enough against the Lyapunov exponent of the map and the spectral gap of the network, the maps synchronize on one chaotic orbit, and the variance of the states over the nodes falls to zero. Below it they stay incoherent.

Parameters:

Name	Type	Description	Default
n_nodes	int	The number of maps, the nodes of the network.	100
network	str	The generator, 'watts_strogatz', 'erdos_renyi', or 'barabasi_albert'.	'erdos_renyi'
degree	int	The mean degree of the network.	20
rewiring	float	The probability of rewiring every edge of the Watts-Strogatz ring.	0.1
r	float	The parameter of the logistic map, chaotic at 4.	4.0
epsilon	float	The coupling between 0 and 1.	0.9
seed	int	The seed of the network and the initial states.	0

Returns:

Name	Type	Description
system	System	The system, with the states maps.x and their variance over the nodes, spread.

Raises:

Type	Description
ValueError	The network is unknown.

Examples:

>>> from skcomplex.systems.physics.synchronization import coupled_map_lattices_on_networks
>>> coupled_map_lattices_on_networks(20).elements[0].prepare().size_
20

Source code in src/skcomplex/systems/physics/synchronization/_networks.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'map', 'chaos', 'network'),
    references=('Arenas et al. 2008', 'Kaneko 1989'),
    defaults={'until': 100.0, 'dt': 1.0},
    observables=lambda system: [ObservableSpec('spread', group='maps', variable='x', reduce='var')],
)
def coupled_map_lattices_on_networks(
    n_nodes: int = 100,
    *,
    network: str = 'erdos_renyi',
    degree: int = 20,
    rewiring: float = 0.1,
    r: float = 4.0,
    epsilon: float = 0.9,
    seed: int = 0,
) -> System:
    """Return chaotic logistic maps on the nodes of a network, each pulled towards the mean of its neighbours.

    Node `i` follows `x_i(t + 1) = (1 - epsilon) f(x_i) + (epsilon / k_i) sum_j A_ij f(x_j)` with the logistic map
    `f(x) = r x (1 - x)`, a mean over the neighbours of the network. For a coupling `epsilon` large enough against the
    Lyapunov exponent of the map and the spectral gap of the network, the maps synchronize on one chaotic orbit, and
    the variance of the states over the nodes falls to zero. Below it they stay incoherent.

    Args:
        n_nodes:
            The number of maps, the nodes of the network.

        network:
            The generator, `'watts_strogatz'`, `'erdos_renyi'`, or `'barabasi_albert'`.

        degree:
            The mean degree of the network.

        rewiring:
            The probability of rewiring every edge of the Watts-Strogatz ring.

        r:
            The parameter of the logistic map, chaotic at 4.

        epsilon:
            The coupling between 0 and 1.

        seed:
            The seed of the network and the initial states.

    Returns:
        system:
            The system, with the states `maps.x` and their variance over the nodes, `spread`.

    Raises:
        ValueError:
            The network is unknown.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import coupled_map_lattices_on_networks
        >>> coupled_map_lattices_on_networks(20).elements[0].prepare().size_
        20
    """
    space = _network(network, n_nodes, degree, rewiring, seed)
    rng = np.random.default_rng(seed + 1)
    maps = LatticeElement(
        'maps',
        element=LogisticNodeElement,
        space=space,
        values={'x': rng.uniform(0.0, 1.0, n_nodes), 'r': r, 'epsilon': epsilon},
        visual=VisualSpec('sites', {'states': ['x'], 'low': 0.0, 'high': 1.0}),
    )
    coupling = LogisticCouplingInteraction('coupling', target='maps')
    return System(elements=[maps], interactions=[coupling], time=DiscreteTime())

coupled_rossler(n_oscillators=2, *, k=0.04, mismatch=0.02, a=0.2, b=0.2, c=5.7, seed=0)

Return chaotic Rössler oscillators on a ring, diffusively coupled through their second variable.

Oscillator i follows x1' = -ω_i x2 - x3, x2' = ω_i x1 + a x2 + k sum_j (x2_j - x2_i), and x3' = b + x3 (x1 - c), summed over its neighbours on the ring. Two oscillators with frequencies 1 ± mismatch pass from no synchronization to chaotic phase synchronization as k grows. Identical oscillators, with no mismatch, synchronize completely when every nonzero eigenvalue of the ring's Laplacian times k falls where the master stability function is negative.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators on the ring.	2
k	float	The coupling strength.	0.04
mismatch	float	The half width of the interval of natural frequencies around 1.	0.02
a	float	The first Rössler parameter.	0.2
b	float	The second Rössler parameter.	0.2
c	float	The third Rössler parameter.	5.7
seed	int	The seed of the initial states.	0

Returns:

Name	Type	Description
system	System	The system, with the states oscillator_<i>.x1, oscillator_<i>.x2, and oscillator_<i>.x3.

Examples:

>>> from skcomplex.systems.physics.synchronization import coupled_rossler
>>> [element.omega for element in coupled_rossler().elements]
[0.98, 1.02]

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'chaos', 'ode', 'master stability function', 'networks'),
    references=('Rössler 1976', 'Pecora and Carroll 1998', 'Datseris and Parlitz 2022', 'Arenas et al. 2008'),
    defaults={'until': 100.0, 'dt': 0.1},
)
def coupled_rossler(
    n_oscillators: int = 2,
    *,
    k: float = 0.04,
    mismatch: float = 0.02,
    a: float = 0.2,
    b: float = 0.2,
    c: float = 5.7,
    seed: int = 0,
) -> System:
    """Return chaotic Rössler oscillators on a ring, diffusively coupled through their second variable.

    Oscillator `i` follows `x1' = -ω_i x2 - x3`, `x2' = ω_i x1 + a x2 + k sum_j (x2_j - x2_i)`, and
    `x3' = b + x3 (x1 - c)`, summed over its neighbours on the ring. Two oscillators with frequencies `1 ± mismatch`
    pass from no synchronization to chaotic phase synchronization as `k` grows. Identical oscillators, with no mismatch,
    synchronize completely when every nonzero eigenvalue of the ring's Laplacian times `k` falls where the master
    stability function is negative.

    Args:
        n_oscillators:
            The number of oscillators on the ring.

        k:
            The coupling strength.

        mismatch:
            The half width of the interval of natural frequencies around 1.

        a:
            The first Rössler parameter.

        b:
            The second Rössler parameter.

        c:
            The third Rössler parameter.

        seed:
            The seed of the initial states.

    Returns:
        system:
            The system, with the states `oscillator_<i>.x1`, `oscillator_<i>.x2`, and `oscillator_<i>.x3`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import coupled_rossler
        >>> [element.omega for element in coupled_rossler().elements]
        [0.98, 1.02]
    """
    omegas = np.linspace(1 - mismatch, 1 + mismatch, n_oscillators)
    starts = np.random.default_rng(seed).uniform(-1.0, 1.0, (n_oscillators, 2))
    oscillators = [
        RosslerOscillatorElement(
            f'oscillator_{index}',
            x1=float(starts[index, 0]),
            x2=float(starts[index, 1]),
            x3=0.0,
            omega=round(float(omegas[index]), 12),
            a=a,
            b=b,
            c=c,
        )
        for index in range(n_oscillators)
    ]
    links = [(index, index + 1) for index in range(n_oscillators - 1)]
    if n_oscillators > PAIR:
        links.append((n_oscillators - 1, 0))
    couplings = [
        DiffusiveCouplingInteraction(
            f'coupling_{source}_{target}',
            source=f'oscillator_{source}',
            target=f'oscillator_{target}',
            k=k,
        )
        for first, second in links
        for source, target in ((first, second), (second, first))
    ]
    return System(elements=oscillators, interactions=couplings)

delayed_kuramoto(n_oscillators=5, *, K=2.0, tau=1.0, spread=0.2, seed=0)

Return Kuramoto phase oscillators coupled all to all through signals that arrive after a delay.

Oscillator i follows theta_i' = omega_i + (K / N) sum_j sin(theta_j(t - tau) - theta_i), with natural frequencies spread evenly around 1 and initial phases spread around the circle.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	5
K	float	The coupling strength.	2.0
tau	float	The delay of the coupling.	1.0
spread	float	The width of the interval of natural frequencies around 1.	0.2
seed	int	The seed of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases oscillator_<i>.theta.

Examples:

>>> from skcomplex.systems.physics.synchronization import delayed_kuramoto
>>> len(delayed_kuramoto(n_oscillators=3).interactions)
6

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'delay', 'oscillators'),
    references=('Yeung and Strogatz 1999', 'Kuramoto 1984'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def delayed_kuramoto(
    n_oscillators: int = 5,
    *,
    K: float = 2.0,
    tau: float = 1.0,
    spread: float = 0.2,
    seed: int = 0,
) -> System:
    """Return Kuramoto phase oscillators coupled all to all through signals that arrive after a delay.

    Oscillator `i` follows `theta_i' = omega_i + (K / N) sum_j sin(theta_j(t - tau) - theta_i)`, with natural
    frequencies spread evenly around 1 and initial phases spread around the circle.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        tau:
            The delay of the coupling.

        spread:
            The width of the interval of natural frequencies around 1.

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The system, with the phases `oscillator_<i>.theta`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import delayed_kuramoto
        >>> len(delayed_kuramoto(n_oscillators=3).interactions)
        6
    """
    phases = [2 * math.pi * ((seed * 0.618 + index / n_oscillators) % 1.0) for index in range(n_oscillators)]
    omegas = [1.0 + spread * (index / max(n_oscillators - 1, 1) - 0.5) for index in range(n_oscillators)]
    elements = [
        PhaseOscillatorElement(f'oscillator_{index}', theta=phases[index], omega=omegas[index])
        for index in range(n_oscillators)
    ]
    interactions = [
        DelayedCouplingInteraction(
            f'coupling_{source}_{target}',
            source=f'oscillator_{source}',
            target=f'oscillator_{target}',
            K=K / n_oscillators,
            tau=tau,
        )
        for source in range(n_oscillators)
        for target in range(n_oscillators)
        if source != target
    ]
    return System(elements=elements, interactions=interactions)

heterogeneous_global_coupling(n_systems=100, *, k=2.0, beta=0.1, gamma=2.0, heterogeneity=0.8, random_state=0)

Return globally coupled oscillators with heterogeneous time scales, whose fixed point the coupling restabilizes.

System j follows x_j' = τ_j g(x_j, y_j) + k (X - x_j) and y_j' = τ_j h(x_j, y_j) + k (Y - y_j), where X and Y are the population means and g + ih = (beta + i gamma) z - |z|² z is the Hopf normal form of z = x + iy. Its fixed point 0 is an unstable focus with eigenvalues beta ± i gamma. The time scales τ_j are spread evenly over [1 - heterogeneity, 1 + heterogeneity], with mean 1 and variance sigma². The order parameters X, Y, V = <δτ δx>, and W = <δτ δy> close on four equations. The global fixed point is stable when their eigenvalues nu - k / 2 ± √(k² / 4 + sigma² nu²), for nu = beta ± i gamma, have negative real parts, which holds roughly for 2 beta (1 + sigma) < k < (gamma² - beta²) sigma² / beta.

Parameters:

Name	Type	Description	Default
n_systems	int	The number of coupled systems.	100
k	float	The coupling constant.	2.0
beta	float	The growth rate β of the unstable focus.	0.1
gamma	float	The rotation rate gamma of the unstable focus.	2.0
heterogeneity	float	The half width of the interval of time scales, below 1.	0.8
random_state	int | Generator | None	The seed or generator of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the states systems.x and systems.y and the observables amplitude, the mean amplitude, and the order parameters X, Y, V, and W.

Examples:

>>> from skcomplex.systems.physics.synchronization import heterogeneous_global_coupling
>>> heterogeneous_global_coupling(4).elements[0].values['tau']
array([0.4, 0.8, 1.2, 1.6])

Source code in src/skcomplex/systems/physics/synchronization/_amplitude_death.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'amplitude death', 'mean field', 'heterogeneity', 'ode'),
    references=("De Monte, d'Ovidio, and Mosekilde 2003", 'Bertin 2021'),
    defaults={'until': 40.0, 'dt': 0.2},
    observables=lambda system: [
        ObservableSpec('amplitude', _amplitude, group='systems'),
        ObservableSpec('X', group='systems', variable='x'),
        ObservableSpec('Y', group='systems', variable='y'),
        ObservableSpec('V', _time_scale_deviation_x, group='systems'),
        ObservableSpec('W', _time_scale_deviation_y, group='systems'),
    ],
)
def heterogeneous_global_coupling(
    n_systems: int = 100,
    *,
    k: float = 2.0,
    beta: float = 0.1,
    gamma: float = 2.0,
    heterogeneity: float = 0.8,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return globally coupled oscillators with heterogeneous time scales, whose fixed point the coupling restabilizes.

    System `j` follows `x_j' = τ_j g(x_j, y_j) + k (X - x_j)` and `y_j' = τ_j h(x_j, y_j) + k (Y - y_j)`, where `X`
    and `Y` are the population means and `g + ih = (beta + i gamma) z - |z|² z` is the Hopf normal form of
    `z = x + iy`. Its fixed point 0 is an unstable focus with eigenvalues `beta ± i gamma`. The time scales `τ_j` are
    spread evenly over `[1 - heterogeneity, 1 + heterogeneity]`, with mean 1 and variance `sigma²`. The order
    parameters `X`, `Y`, `V = <δτ δx>`, and `W = <δτ δy>` close on four equations. The global fixed point is stable
    when their eigenvalues `nu - k / 2 ± √(k² / 4 + sigma² nu²)`, for `nu = beta ± i gamma`, have negative real parts,
    which holds roughly for `2 beta (1 + sigma) < k < (gamma² - beta²) sigma² / beta`.

    Args:
        n_systems:
            The number of coupled systems.

        k:
            The coupling constant.

        beta:
            The growth rate `β` of the unstable focus.

        gamma:
            The rotation rate `gamma` of the unstable focus.

        heterogeneity:
            The half width of the interval of time scales, below 1.

        random_state:
            The seed or generator of the initial phases.

    Returns:
        system:
            The system, with the states `systems.x` and `systems.y` and the observables `amplitude`, the mean amplitude,
            and the order parameters `X`, `Y`, `V`, and `W`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import heterogeneous_global_coupling
        >>> heterogeneous_global_coupling(4).elements[0].values['tau']
        array([0.4, 0.8, 1.2, 1.6])
    """
    taus = 1 + heterogeneity * (2 * (np.arange(n_systems) + 1 / 2) / n_systems - 1)
    x, y = _spread_on_circle(n_systems, random_state)
    radius = math.sqrt(abs(beta))
    systems = GroupElement(
        'systems',
        element=StuartLandauElement,
        size=n_systems,
        values={'x': radius * x, 'y': radius * y, 'mu': beta, 'omega': gamma, 'tau': taus},
    )
    return System(
        elements=[systems],
        interactions=[
            GlobalCouplingInteraction(f'coupling_{state}', target='systems', state=state, K=k) for state in 'xy'
        ],
    )

hypergraph_kuramoto(n_oscillators=60, *, K1=0.0, K2=8.0, width=0.5, links=1.0, triangles=1.0, phases=None, seed=0)

Return Kuramoto oscillators coupled in pairs and by three-body interactions on the triangles of a hypergraph.

Oscillator i follows theta_i' = omega_i + (K1 / <k1>) sum_j sin(theta_j - theta_i) + (K2 / <k2>) sum_jk sin(theta_j + theta_k - 2 theta_i) over the links ij and the triangles ijk of a random hypergraph, where <k1> and <k2> are the mean numbers of links and of triangles of a node. All-to-all, the triangle term is the three-body coupling of Skardal and Arenas, and the Ott-Antonsen reduction of Lorentzian frequencies of half width width gives r' = -width r + K1 (r - r^3) / 2 + K2 (r^3 - r^5) / 2. Without links, the incoherent state stays stable for every K2, and for K2 > 8 width a synchronized state of r^2 = (1 + sqrt(1 - 8 width / K2)) / 2 coexists with it, so the order parameter depends on the initial phases and sweeps of the coupling show hysteresis.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators, the nodes of the hypergraph.	60
K1	float	The pairwise coupling, divided by the mean number of links of a node.	0.0
K2	float	The three-body coupling, divided by the mean number of triangles of a node.	8.0
width	float	The half width of the Lorentzian distribution of the natural frequencies.	0.5
links	float	The probability that a pair of oscillators is a link.	1.0
triangles	float	The probability that a triple of oscillators is a triangle.	1.0
phases	Sequence[float] | None	The initial phases, or None for phases drawn uniformly from the circle.	None
seed	int	The seed of the hypergraph, the frequencies, and the phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases oscillators.theta and the order parameter order.

Examples:

>>> from skcomplex.systems.physics.synchronization import hypergraph_kuramoto
>>> system = hypergraph_kuramoto(6)
>>> len(system.elements[0].space.prepare().hyperedges_)
35

Source code in src/skcomplex/systems/physics/synchronization/_hypergraphs.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'hypergraph'),
    references=('Skardal and Arenas 2019', 'Majhi et al. 2022'),
    defaults={'until': 40.0, 'dt': 0.1},
    observables=lambda system: order_parameter('order', 'oscillators'),
)
def hypergraph_kuramoto(
    n_oscillators: int = 60,
    *,
    K1: float = 0.0,
    K2: float = 8.0,
    width: float = 0.5,
    links: float = 1.0,
    triangles: float = 1.0,
    phases: Sequence[float] | None = None,
    seed: int = 0,
) -> System:
    """Return Kuramoto oscillators coupled in pairs and by three-body interactions on the triangles of a hypergraph.

    Oscillator `i` follows `theta_i' = omega_i + (K1 / <k1>) sum_j sin(theta_j - theta_i) + (K2 / <k2>) sum_jk
    sin(theta_j + theta_k - 2 theta_i)` over the links `ij` and the triangles `ijk` of a random hypergraph, where `<k1>`
    and `<k2>` are the mean numbers of links and of triangles of a node. All-to-all, the triangle term is the
    three-body coupling of Skardal and Arenas, and the Ott-Antonsen reduction of Lorentzian frequencies of half width
    `width` gives `r' = -width r + K1 (r - r^3) / 2 + K2 (r^3 - r^5) / 2`. Without links, the incoherent state stays
    stable for every `K2`, and for `K2 > 8 width` a synchronized state of `r^2 = (1 + sqrt(1 - 8 width / K2)) / 2`
    coexists with it, so the order parameter depends on the initial phases and sweeps of the coupling show hysteresis.

    Args:
        n_oscillators:
            The number of oscillators, the nodes of the hypergraph.

        K1:
            The pairwise coupling, divided by the mean number of links of a node.

        K2:
            The three-body coupling, divided by the mean number of triangles of a node.

        width:
            The half width of the Lorentzian distribution of the natural frequencies.

        links:
            The probability that a pair of oscillators is a link.

        triangles:
            The probability that a triple of oscillators is a triangle.

        phases:
            The initial phases, or `None` for phases drawn uniformly from the circle.

        seed:
            The seed of the hypergraph, the frequencies, and the phases.

    Returns:
        system:
            The system, with the phases `oscillators.theta` and the order parameter `order`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import hypergraph_kuramoto
        >>> system = hypergraph_kuramoto(6)
        >>> len(system.elements[0].space.prepare().hyperedges_)
        35
    """
    hypergraph = _simplices(n_oscillators, links, triangles, seed)
    space = NetworkSpace(nx.complete_graph(n_oscillators), hyperedges=hypergraph)
    sizes = [len(members) for members in hypergraph.edges.members(dtype=list)]
    mean_links = max(PAIR * sizes.count(PAIR) / n_oscillators, 1.0)
    mean_triangles = max(3 * sizes.count(3) / n_oscillators, 1.0)
    rng = np.random.default_rng(seed + 1)
    omega = width * stats.cauchy.ppf(rng.uniform(size=n_oscillators))
    theta = rng.uniform(0, 2 * math.pi, n_oscillators) if phases is None else np.asarray(phases, dtype=float)
    oscillators = LatticeElement(
        'oscillators',
        element=PhaseOscillatorElement,
        space=space,
        values={'theta': theta, 'omega': omega},
        visual=VisualSpec('sites', {'states': ['theta'], 'low': 0.0, 'high': 2 * math.pi}),
    )
    coupling = SimplexCouplingInteraction(
        'coupling',
        target='oscillators',
        K1=K1 / mean_links,
        K2=K2 / mean_triangles,
    )
    return System(elements=[oscillators], interactions=[coupling])

kuramoto(n_oscillators=50, *, K=4.0, sigma=1.0, seed=0)

Return Kuramoto oscillators coupled all to all, theta_i' = omega_i + (K / N) sum_j sin(theta_j - theta_i).

Every oscillator reads the group means of cos theta and sin theta, so the coupling costs time linear in the number of oscillators. The natural frequencies are drawn from a normal distribution of width sigma, and the phases uniformly from the circle. Above the critical coupling K_c = 2 sigma sqrt(2 / pi) a synchronized cluster forms, and the order parameter r, the modulus of the mean of exp(i theta), grows from about 1 / sqrt(N) towards 1.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	50
K	float	The coupling strength.	4.0
sigma	float	The standard deviation of the natural frequencies.	1.0
seed	int	The seed of the natural frequencies and the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases oscillators.theta and the order parameter order.

Examples:

>>> from skcomplex.systems.physics.synchronization import kuramoto
>>> system = kuramoto(10)
>>> system.elements[0].size, system.space
(10, None)

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'mean field'),
    references=('Kuramoto 1984', 'Datseris and Parlitz 2022', 'Gros 2015', 'Jensen 2023'),
    defaults={'until': 50.0, 'dt': 0.5},
    observables=lambda system: order_parameter('order', 'oscillators'),
)
def kuramoto(n_oscillators: int = 50, *, K: float = 4.0, sigma: float = 1.0, seed: int = 0) -> System:
    """Return Kuramoto oscillators coupled all to all, `theta_i' = omega_i + (K / N) sum_j sin(theta_j - theta_i)`.

    Every oscillator reads the group means of `cos theta` and `sin theta`, so the coupling costs time linear in the
    number of oscillators. The natural frequencies are drawn from a normal distribution of width `sigma`, and the
    phases uniformly from the circle. Above the critical coupling `K_c = 2 sigma sqrt(2 / pi)` a synchronized cluster
    forms, and the order parameter `r`, the modulus of the mean of `exp(i theta)`, grows from about `1 / sqrt(N)`
    towards 1.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        sigma:
            The standard deviation of the natural frequencies.

        seed:
            The seed of the natural frequencies and the initial phases.

    Returns:
        system:
            The system, with the phases `oscillators.theta` and the order parameter `order`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import kuramoto
        >>> system = kuramoto(10)
        >>> system.elements[0].size, system.space
        (10, None)
    """
    rng = np.random.default_rng(seed)
    oscillators = GroupElement(
        'oscillators',
        element=PhaseOscillatorElement,
        size=n_oscillators,
        values={'theta': _spread_phases(n_oscillators, seed), 'omega': rng.normal(0.0, sigma, n_oscillators)},
    )
    coupling = KuramotoCouplingInteraction(
        'coupling',
        target='oscillators',
        K=K,
        mode='group',
        reduce='mean',
        include_self=True,
    )
    return System(
        elements=[oscillators],
        interactions=[coupling],
    )

kuramoto_on_networks(n_oscillators=200, *, network='watts_strogatz', degree=6, rewiring=0.1, K=1.0, sigma=0.5, seed=0)

Return Kuramoto oscillators coupled along the edges of a network, theta_i' = omega_i + K sum_j A_ij sin(...).

Oscillator i follows theta_i' = omega_i + (K / <k>) sum_j w_ij sin(theta_j - theta_i) over its neighbours j, where the weights w_ij = 1 live on a group on the edges and <k> is the mean degree. The coupling runs once per edge. The natural frequencies are drawn from a normal distribution of width sigma, and the phases uniformly from the circle. Above a critical coupling of the order of the all-to-all value 2 sigma sqrt(2 / pi), the order parameter r = |<exp(i theta)>| grows from about 1 / sqrt(N) towards 1, and small-world shortcuts lower it.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators, the nodes of the network.	200
network	str	The generator, 'watts_strogatz', 'erdos_renyi', or 'barabasi_albert'.	'watts_strogatz'
degree	int	The mean degree of the network, even for the Watts-Strogatz ring.	6
rewiring	float	The probability of rewiring every edge of the Watts-Strogatz ring.	0.1
K	float	The coupling strength, divided by the mean degree.	1.0
sigma	float	The standard deviation of the natural frequencies.	0.5
seed	int	The seed of the network, the frequencies, and the phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases oscillators.theta, the weights links.w, and the order parameter order.

Raises:

Type	Description
ValueError	The network is unknown.

Examples:

>>> from skcomplex.systems.physics.synchronization import kuramoto_on_networks
>>> system = kuramoto_on_networks(20, degree=4)
>>> system.elements[1].prepare().size_
40

Source code in src/skcomplex/systems/physics/synchronization/_networks.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'network'),
    references=('Arenas et al. 2008', 'Jensen 2023'),
    defaults={'until': 40.0, 'dt': 0.2},
    observables=lambda system: order_parameter('order', 'oscillators'),
)
def kuramoto_on_networks(
    n_oscillators: int = 200,
    *,
    network: str = 'watts_strogatz',
    degree: int = 6,
    rewiring: float = 0.1,
    K: float = 1.0,
    sigma: float = 0.5,
    seed: int = 0,
) -> System:
    """Return Kuramoto oscillators coupled along the edges of a network, `theta_i' = omega_i + K sum_j A_ij sin(...)`.

    Oscillator `i` follows `theta_i' = omega_i + (K / <k>) sum_j w_ij sin(theta_j - theta_i)` over its neighbours `j`,
    where the weights `w_ij = 1` live on a group on the edges and `<k>` is the mean degree. The coupling runs once per
    edge. The natural frequencies are drawn from a normal distribution of width `sigma`, and the phases uniformly from
    the circle. Above a critical coupling of the order of the all-to-all value `2 sigma sqrt(2 / pi)`, the order
    parameter `r = |<exp(i theta)>|` grows from about `1 / sqrt(N)` towards 1, and small-world shortcuts lower it.

    Args:
        n_oscillators:
            The number of oscillators, the nodes of the network.

        network:
            The generator, `'watts_strogatz'`, `'erdos_renyi'`, or `'barabasi_albert'`.

        degree:
            The mean degree of the network, even for the Watts-Strogatz ring.

        rewiring:
            The probability of rewiring every edge of the Watts-Strogatz ring.

        K:
            The coupling strength, divided by the mean degree.

        sigma:
            The standard deviation of the natural frequencies.

        seed:
            The seed of the network, the frequencies, and the phases.

    Returns:
        system:
            The system, with the phases `oscillators.theta`, the weights `links.w`, and the order parameter `order`.

    Raises:
        ValueError:
            The network is unknown.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import kuramoto_on_networks
        >>> system = kuramoto_on_networks(20, degree=4)
        >>> system.elements[1].prepare().size_
        40
    """
    space = _network(network, n_oscillators, degree, rewiring, seed)
    rng = np.random.default_rng(seed + 1)
    mean_degree = 2 * space.graph.number_of_edges() / n_oscillators
    oscillators = LatticeElement(
        'oscillators',
        element=PhaseOscillatorElement,
        space=space,
        values={'theta': rng.uniform(0, 2 * math.pi, n_oscillators), 'omega': rng.normal(0.0, sigma, n_oscillators)},
        visual=VisualSpec('sites', {'states': ['theta'], 'low': 0.0, 'high': 2 * math.pi}),
    )
    links = GroupElement('links', element=LinkElement, space=space, on='edges')
    coupling = EdgeCouplingInteraction('coupling', target='oscillators', links='links', K=K / mean_degree)
    return System(elements=[oscillators, links], interactions=[coupling])

kuramoto_sakaguchi_chimera(n_oscillators=32, *, A=0.2, beta=0.1, spread=0.1, seed=0)

Return two populations of identical phase oscillators with phase lag, where one synchronizes and one drifts.

Oscillator i of population p follows theta_i' = sum_q (K_pq / N) sum_j sin(theta_j - theta_i - alpha), with the stronger coupling mu = (1 + A) / 2 inside a population, the weaker nu = (1 - A) / 2 between them, and the lag alpha = pi / 2 - beta. Every population reads the means of cos theta and sin theta of both populations. The first population starts nearly synchronized and the second spread over the circle. The populations settle in a chimera, where one is fully synchronized and the other keeps drifting with a partial order parameter, which breathes for A = 0.2. Which population ends up coherent depends on the initial phases.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators of every population.	32
A	float	The disparity of the couplings inside and between the populations.	0.2
beta	float	The distance of the phase lag from pi / 2.	0.1
spread	float	The width of the initial phases of the first population.	0.1
seed	int	The seed of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases first.theta and second.theta and the order parameters first_order and second_order.

Examples:

>>> from skcomplex.systems.physics.synchronization import kuramoto_sakaguchi_chimera
>>> len(kuramoto_sakaguchi_chimera(4).interactions)
4

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'chimera'),
    references=('Abrams et al. 2008', 'Jensen 2023'),
    defaults={'until': 100.0, 'dt': 1.0},
    observables=lambda system: [*order_parameter('first_order', 'first'), *order_parameter('second_order', 'second')],
)
def kuramoto_sakaguchi_chimera(
    n_oscillators: int = 32,
    *,
    A: float = 0.2,
    beta: float = 0.1,
    spread: float = 0.1,
    seed: int = 0,
) -> System:
    """Return two populations of identical phase oscillators with phase lag, where one synchronizes and one drifts.

    Oscillator `i` of population `p` follows `theta_i' = sum_q (K_pq / N) sum_j sin(theta_j - theta_i - alpha)`, with
    the stronger coupling `mu = (1 + A) / 2` inside a population, the weaker `nu = (1 - A) / 2` between them, and the
    lag `alpha = pi / 2 - beta`. Every population reads the means of `cos theta` and `sin theta` of both populations.
    The first population starts nearly synchronized and the second spread over the circle. The populations settle in a
    chimera, where one is fully synchronized and the other keeps drifting with a partial order parameter, which
    breathes for `A = 0.2`. Which population ends up coherent depends on the initial phases.

    Args:
        n_oscillators:
            The number of oscillators of every population.

        A:
            The disparity of the couplings inside and between the populations.

        beta:
            The distance of the phase lag from `pi / 2`.

        spread:
            The width of the initial phases of the first population.

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The system, with the phases `first.theta` and `second.theta` and the order parameters `first_order` and
            `second_order`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import kuramoto_sakaguchi_chimera
        >>> len(kuramoto_sakaguchi_chimera(4).interactions)
        4
    """
    rng = np.random.default_rng(seed)
    alpha, mu, nu = math.pi / 2 - beta, (1 + A) / 2, (1 - A) / 2
    populations = [
        GroupElement(
            'first',
            element=PhaseOscillatorElement,
            size=n_oscillators,
            values={'theta': rng.normal(0.0, spread, n_oscillators), 'omega': 0.0},
            visual=VisualSpec('phase_ring', {'phase': 'theta'}),
        ),
        GroupElement(
            'second',
            element=PhaseOscillatorElement,
            size=n_oscillators,
            values={'theta': _spread_phases(n_oscillators, seed + 1), 'omega': 0.0},
            visual=VisualSpec('phase_ring', {'phase': 'theta'}),
        ),
    ]
    interactions: list[BaseInteraction] = [
        KuramotoCouplingInteraction(
            f'{target}_from_{source}',
            target=target,
            source=source,
            K=mu if target == source else nu,
            alpha=alpha,
            reduce='mean',
            mode='group',
            include_self=True,
        )
        for target in ('first', 'second')
        for source in ('first', 'second')
    ]
    return System(
        elements=populations,
        interactions=interactions,
    )

mirollo_strogatz(n_oscillators=5, *, S=2.0, gamma=1.0, epsilon=0.1, seed=0)

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.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	5
S	float	The charging input.	2.0
gamma	float	The leak rate, below S so that the oscillators fire.	1.0
epsilon	float	The strength of a pulse.	0.1
seed	int	The seed of the initial charges.	0

Returns:

Name	Type	Description
system	System	The system, with the charges oscillators.x of the group members and the firings as fire events.

Examples:

>>> from skcomplex.systems.physics.synchronization import mirollo_strogatz
>>> mirollo_strogatz(3).elements[0].prepare().size_
3

Source code in src/skcomplex/systems/physics/synchronization/_pulses.py

@register_system(
    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 system, with the charges `oscillators.x` of the group members and the firings as `fire` events.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import mirollo_strogatz
        >>> mirollo_strogatz(3).elements[0].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])

nonlocal_ring_chimera(n_oscillators=100, *, R=14, alpha=1.46, omega=0.0, seed=0)

Return identical phase oscillators on a ring, each coupled to its R nearest neighbours on both sides.

Oscillator i follows phi_i' = omega - (1 / 2R) sum_{j = i - R}^{i + R} sin(phi_i - phi_j + alpha), where the sum includes the oscillator itself. Near alpha = pi / 2, coherent and incoherent domains coexist on the ring as a chimera state, which is a long transient on small rings.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators on the ring.	100
R	int	The range of the coupling.	14
alpha	float	The phase lag.	1.46
omega	float	The natural frequency.	0.0
seed	int	The seed of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases ring.theta and the places of the oscillators ring.x.

Examples:

>>> from skcomplex.systems.physics.synchronization import nonlocal_ring_chimera
>>> nonlocal_ring_chimera(20, R=3).interactions[0].neighborhood.radius
3.5

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'chimera', 'nonlocal coupling'),
    references=('Kuramoto and Battogtokh 2002', 'Datseris and Parlitz 2022'),
    defaults={'until': 100.0, 'dt': 1.0},
    observables=lambda system: order_parameter('order', 'ring'),
)
def nonlocal_ring_chimera(
    n_oscillators: int = 100,
    *,
    R: int = 14,
    alpha: float = 1.46,
    omega: float = 0.0,
    seed: int = 0,
) -> System:
    """Return identical phase oscillators on a ring, each coupled to its `R` nearest neighbours on both sides.

    Oscillator `i` follows `phi_i' = omega - (1 / 2R) sum_{j = i - R}^{i + R} sin(phi_i - phi_j + alpha)`, where the sum
    includes the oscillator itself. Near `alpha = pi / 2`, coherent and incoherent domains coexist on the ring as a
    chimera state, which is a long transient on small rings.

    Args:
        n_oscillators:
            The number of oscillators on the ring.

        R:
            The range of the coupling.

        alpha:
            The phase lag.

        omega:
            The natural frequency.

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The system, with the phases `ring.theta` and the places of the oscillators `ring.x`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import nonlocal_ring_chimera
        >>> nonlocal_ring_chimera(20, R=3).interactions[0].neighborhood.radius
        3.5
    """
    weight = 1 / (2 * R)
    ring = GroupElement(
        'ring',
        element=SitedOscillatorElement,
        size=n_oscillators,
        values={
            'theta': _spread_phases(n_oscillators, seed),
            'omega': omega,
            'x': np.arange(n_oscillators, dtype=float).reshape(-1, 1),
        },
        visual=VisualSpec('phase_ring', {'phase': 'theta', 'position': 'x'}),
    )
    return System(
        elements=[ring],
        interactions=[
            KuramotoCouplingInteraction(
                'coupling',
                target='ring',
                K=weight,
                alpha=alpha,
                include_self=True,
                neighborhood=RadiusNeighborhood(R + 0.5),
            ),
        ],
        space=ContinuousSpace((float(n_oscillators),), boundary='periodic'),
    )

oscillator_death(n_oscillators=100, *, K=3.0, width=1.5, distribution='lorentzian', random_state=0)

Return globally coupled Stuart-Landau oscillators, whose oscillations all die for strong coupling and spread.

Oscillator j follows z_j' = (1 + iω_j) z_j - |z_j|² z_j + (K / N) Σ_k (z_k - z_j) for z_j = x_j + i y_j. Alone it settles on the unit circle. The frequencies are the midpoint quantiles of a density of width ω0, the Lorentzian ω0 / (π (ω² + ω0²)) or the quartic g0 / (ω⁴ + ω0⁴). For many oscillators, the common fixed point z = 0 is stable for K > 1 and ω0 above a critical width, 1 for the Lorentzian and (1 + √(2K - 1)) / √2 for the quartic density, and every amplitude then decays to zero. For the Lorentzian, the slowest perturbation decays at the rate ω0 - 1.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	100
K	float	The coupling strength.	3.0
width	float	The width ω0 of the density of frequencies.	1.5
distribution	str	The density of frequencies, 'lorentzian' or 'quartic'.	'lorentzian'
random_state	int | Generator | None	The seed or generator of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the states oscillators.x and oscillators.y and the observables amplitude, the mean amplitude, X and Y, the mean field, and order, its modulus.

Raises:

Type	Description
ValidationError	The distribution is neither 'lorentzian' nor 'quartic'.

Examples:

>>> from skcomplex.systems.physics.synchronization import oscillator_death
>>> oscillator_death(4, distribution='quartic').elements[0].prepare().size_
4

Source code in src/skcomplex/systems/physics/synchronization/_amplitude_death.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'amplitude death', 'mean field', 'ode'),
    references=('Ermentrout 1990', 'Mirollo and Strogatz 1990', 'Bertin 2021'),
    defaults={'until': 20.0, 'dt': 0.2},
    observables=lambda system: [
        ObservableSpec('amplitude', _amplitude, group='oscillators'),
        ObservableSpec('X', group='oscillators', variable='x'),
        ObservableSpec('Y', group='oscillators', variable='y'),
        ObservableSpec('order', _order),
    ],
)
def oscillator_death(
    n_oscillators: int = 100,
    *,
    K: float = 3.0,
    width: float = 1.5,
    distribution: str = 'lorentzian',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return globally coupled Stuart-Landau oscillators, whose oscillations all die for strong coupling and spread.

    Oscillator `j` follows `z_j' = (1 + iω_j) z_j - |z_j|² z_j + (K / N) Σ_k (z_k - z_j)` for `z_j = x_j + i y_j`.
    Alone it settles on the unit circle. The frequencies are the midpoint quantiles of a density of width `ω0`, the
    Lorentzian `ω0 / (π (ω² + ω0²))` or the quartic `g0 / (ω⁴ + ω0⁴)`. For many oscillators, the common fixed point
    `z = 0` is stable for `K > 1` and `ω0` above a critical width, 1 for the Lorentzian and
    `(1 + √(2K - 1)) / √2` for the quartic density, and every amplitude then decays to zero. For the Lorentzian, the
    slowest perturbation decays at the rate `ω0 - 1`.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        width:
            The width `ω0` of the density of frequencies.

        distribution:
            The density of frequencies, `'lorentzian'` or `'quartic'`.

        random_state:
            The seed or generator of the initial phases.

    Returns:
        system:
            The system, with the states `oscillators.x` and `oscillators.y` and the observables `amplitude`, the mean
            amplitude, `X` and `Y`, the mean field, and `order`, its modulus.

    Raises:
        ValidationError:
            The distribution is neither `'lorentzian'` nor `'quartic'`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import oscillator_death
        >>> oscillator_death(4, distribution='quartic').elements[0].prepare().size_
        4
    """
    omegas = _frequency_quantiles(n_oscillators, width, distribution)
    x, y = _spread_on_circle(n_oscillators, random_state)
    oscillators = GroupElement(
        'oscillators',
        element=StuartLandauElement,
        size=n_oscillators,
        values={'x': x, 'y': y, 'mu': 1.0, 'omega': omegas, 'tau': 1.0},
    )
    return System(
        elements=[oscillators],
        interactions=[
            GlobalCouplingInteraction(f'coupling_{state}', target='oscillators', state=state, K=K) for state in 'xy'
        ],
    )

ott_antonsen(r=0.1, *, K=2.0, Delta=0.5, omega=0.0)

Return the Ott-Antonsen equation for the order parameter of infinitely many Kuramoto oscillators.

For natural frequencies with a Lorentzian distribution of center omega and half width Delta, the complex order parameter z = x + i y follows exactly z' = (i omega - Delta) z + (K / 2)(1 - |z|²) z, a Stuart-Landau oscillator with the growth rate K / 2 - Delta and the saturation K / 2. Its modulus settles at sqrt(1 - 2 Delta / K) above the critical coupling K_c = 2 Delta, and at 0 below it.

Parameters:

Name	Type	Description	Default
r	float	The initial modulus of the order parameter.	0.1
K	float	The coupling strength.	2.0
Delta	float	The half width of the distribution of natural frequencies.	0.5
omega	float	The center of the distribution of natural frequencies.	0.0

Returns:

Name	Type	Description
system	System	The system, with the real part field.x and the imaginary part field.y of the order parameter.

Examples:

>>> from skcomplex.systems.physics.synchronization import ott_antonsen
>>> ott_antonsen(K=3.0).elements[0].kappa
1.5

Source code in src/skcomplex/systems/physics/synchronization/_kuramoto.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'ode', 'mean field', 'dimension reduction'),
    references=('Ott and Antonsen 2008', 'Bick et al. 2020'),
    defaults={'until': 30.0, 'dt': 0.1},
)
def ott_antonsen(
    r: float = 0.1,
    *,
    K: float = 2.0,
    Delta: float = 0.5,
    omega: float = 0.0,
) -> System:
    """Return the Ott-Antonsen equation for the order parameter of infinitely many Kuramoto oscillators.

    For natural frequencies with a Lorentzian distribution of center `omega` and half width `Delta`, the complex order
    parameter `z = x + i y` follows exactly `z' = (i omega - Delta) z + (K / 2)(1 - |z|²) z`, a Stuart-Landau
    oscillator with the growth rate `K / 2 - Delta` and the saturation `K / 2`. Its modulus settles at
    `sqrt(1 - 2 Delta / K)` above the critical coupling `K_c = 2 Delta`, and at 0 below it.

    Args:
        r:
            The initial modulus of the order parameter.

        K:
            The coupling strength.

        Delta:
            The half width of the distribution of natural frequencies.

        omega:
            The center of the distribution of natural frequencies.

    Returns:
        system:
            The system, with the real part `field.x` and the imaginary part `field.y` of the order parameter.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import ott_antonsen
        >>> ott_antonsen(K=3.0).elements[0].kappa
        1.5
    """
    return System(elements=[StuartLandauElement('field', x=r, y=0.0, mu=K / 2 - Delta, omega=omega, kappa=K / 2)])

rossler_driven_lorenz(k=8.0, *, auxiliary=True, p=1 / 6)

Return a chaotic Rössler system that drives a Lorenz-63 system through its second variable.

The drive follows p x1' = -x2 - x3, p x2' = x1 + 0.2 x2, p x3' = 0.2 + x3 (x1 - 5.7), with the time scale p = 1 / 6 that brings it close to the Lorenz system, and the response follows y1' = 10 (y2 - y1), y2' = 28 y1 - y2 - y1 y3 + k x2, y3' = y1 y2 - 2.666 y3. The two systems have nothing in common, yet for a coupling above k ≈ 6.66 the response forgets its initial state and becomes a function of the state of the drive, generalized synchronization. With auxiliary, a second copy of the response from another initial state receives the same drive, and the two copies converge to each other exactly when generalized synchronization holds, the auxiliary system method.

Parameters:

Name	Type	Description	Default
k	float	The strength of the coupling.	8.0
auxiliary	bool	Whether a second copy of the response, auxiliary, receives the drive too.	True
p	float	The time scale of the drive.	1 / 6

Returns:

Name	Type	Description
system	System	The system, with the drive drive.x1, drive.x2, drive.x3, the response response.x, response.y, response.z, and with auxiliary the copy auxiliary.x, auxiliary.y, auxiliary.z.

Examples:

>>> from skcomplex.systems.physics.synchronization import rossler_driven_lorenz
>>> [element.label for element in rossler_driven_lorenz(auxiliary=False).elements]
['drive', 'response']

Source code in src/skcomplex/systems/physics/synchronization/_chaos.py

@register_system(
    field='physics/synchronization',
    tags=('chaos', 'generalized synchronization', 'drive-response', 'ode'),
    references=('Rulkov et al. 1995', 'Abarbanel, Rulkov, and Sushchik 1996', 'Datseris and Parlitz 2022'),
    defaults={'until': 100.0, 'dt': 0.01},
)
def rossler_driven_lorenz(k: float = 8.0, *, auxiliary: bool = True, p: float = 1 / 6) -> System:
    """Return a chaotic Rössler system that drives a Lorenz-63 system through its second variable.

    The drive follows `p x1' = -x2 - x3`, `p x2' = x1 + 0.2 x2`, `p x3' = 0.2 + x3 (x1 - 5.7)`, with the time scale
    `p = 1 / 6` that brings it close to the Lorenz system, and the response follows `y1' = 10 (y2 - y1)`,
    `y2' = 28 y1 - y2 - y1 y3 + k x2`, `y3' = y1 y2 - 2.666 y3`. The two systems have nothing in common, yet for a
    coupling above `k ≈ 6.66` the response forgets its initial state and becomes a function of the state of the drive,
    generalized synchronization. With `auxiliary`, a second copy of the response from another initial state receives
    the same drive, and the two copies converge to each other exactly when generalized synchronization holds, the
    auxiliary system method.

    Args:
        k:
            The strength of the coupling.

        auxiliary:
            Whether a second copy of the response, `auxiliary`, receives the drive too.

        p:
            The time scale of the drive.

    Returns:
        system:
            The system, with the drive `drive.x1`, `drive.x2`, `drive.x3`, the response `response.x`, `response.y`,
            `response.z`, and with `auxiliary` the copy `auxiliary.x`, `auxiliary.y`, `auxiliary.z`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import rossler_driven_lorenz
        >>> [element.label for element in rossler_driven_lorenz(auxiliary=False).elements]
        ['drive', 'response']
    """
    drive = RosslerOscillatorElement('drive', x1=1.0, x2=0.0, x3=0.0, p=p)
    responses = [LorenzElement('response', x=1.0, y=1.0, z=1.0, beta=2.666)]
    if auxiliary:
        responses.append(LorenzElement('auxiliary', x=-5.0, y=3.0, z=20.0, beta=2.666))
    couplings = [
        RateInteraction(
            f'coupling_{response.label}',
            target=response.label,
            state='y',
            sources=['drive'],
            rate=lambda s: k * s.sources[0].x2,
        )
        for response in responses
    ]
    return System(elements=[drive, *responses], interactions=couplings)

swarmalators(n_units=100, *, J=1.0, K=0.0, A=1.0, B=1.0, random_state=0)

Return swarmalators, identical oscillators whose phases set their attraction and whose positions their coupling.

Every swarmalator i moves by (1 / N) sum_j [(x_j - x_i) / |x_j - x_i| (A + J cos(theta_j - theta_i)) - B (x_j - x_i) / |x_j - x_i|^2] and turns by (K / N) sum_j sin(theta_j - theta_i) / |x_j - x_i|. With J > 0 like phases attract. With K = 0 and J = 1 the swarm settles into a static phase wave, a disk whose phase winds once around its center, so one of the order parameters S± = |mean(exp(i (phi ± theta)))| of the spatial angle phi and the phase approaches 1. With K > 0 the phases synchronize into a uniform disk, and with K < 0 and small J they stay asynchronous.

Parameters:

Name	Type	Description	Default
n_units	int	The number of swarmalators N.	100
J	float	The influence of the phase similarity on the attraction.	1.0
K	float	The phase coupling.	0.0
A	float	The strength of the attraction.	1.0
B	float	The strength of the repulsion.	1.0
random_state	int | Generator | None	The seed or generator of the initial positions in the square of side 2 and of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the positions units.x and the phases units.theta.

Examples:

>>> from skcomplex.systems.physics.synchronization import swarmalators
>>> swarmalators(10, K=-0.1).interactions[0].K
-0.1

Source code in src/skcomplex/systems/physics/synchronization/_swarmalators.py

@register_system(
    field='physics/synchronization',
    tags=('swarmalators', 'synchronization', 'swarming', 'pair interaction', 'ode'),
    references=("O'Keeffe, Ceron, and Petersen 2017", 'te Vrugt and Wittkowski 2024'),
    defaults={'until': 100.0, 'dt': 1.0, 'solver': RK4Solver(dt=0.1)},
)
def swarmalators(
    n_units: int = 100,
    *,
    J: float = 1.0,
    K: float = 0.0,
    A: float = 1.0,
    B: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return swarmalators, identical oscillators whose phases set their attraction and whose positions their coupling.

    Every swarmalator `i` moves by `(1 / N) sum_j [(x_j - x_i) / |x_j - x_i| (A + J cos(theta_j - theta_i)) -
    B (x_j - x_i) / |x_j - x_i|^2]` and turns by `(K / N) sum_j sin(theta_j - theta_i) / |x_j - x_i|`. With `J > 0`
    like phases attract. With `K = 0` and `J = 1` the swarm settles into a static phase wave, a disk whose phase winds
    once around its center, so one of the order parameters `S± = |mean(exp(i (phi ± theta)))|` of the spatial angle
    `phi` and the phase approaches 1. With `K > 0` the phases synchronize into a uniform disk, and with `K < 0` and
    small `J` they stay asynchronous.

    Args:
        n_units:
            The number of swarmalators `N`.

        J:
            The influence of the phase similarity on the attraction.

        K:
            The phase coupling.

        A:
            The strength of the attraction.

        B:
            The strength of the repulsion.

        random_state:
            The seed or generator of the initial positions in the square of side 2 and of the initial phases.

    Returns:
        system:
            The system, with the positions `units.x` and the phases `units.theta`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import swarmalators
        >>> swarmalators(10, K=-0.1).interactions[0].K
        -0.1
    """
    rng = np.random.default_rng(random_state)
    units = GroupElement(
        'units',
        element=SwarmalatorElement,
        size=n_units,
        values={'x': rng.uniform(-1.0, 1.0, (n_units, 2)), 'theta': rng.uniform(-np.pi, np.pi, n_units), 'omega': 0.0},
        visual=VisualSpec('points', {'position': 'x', 'color_by': 'theta', 'radius': 0.03}),
    )
    coupling = SwarmalatorCouplingInteraction(
        'coupling',
        first='units',
        second='units',
        A=A,
        B=B,
        J=J,
        K=K,
        n_units=n_units,
    )
    return System(elements=[units], interactions=[coupling], space=ContinuousSpace(n_dim=2))

topological_kuramoto(n_nodes=30, *, density=0.4, coupling=2.0, width=1.0, seed=0, model='NL')

Return phases on the nodes and links of a clique complex, coupled through its incidence by order parameters.

Node phases follow theta' = omega + K R_links codifferential(sin(d theta)). Link phases follow phi' = omega + K R_nodes d(sin(codifferential(phi))) - K codifferential(sin(d phi)). With the signs of codifferential, these are the Nodes-Links system of Ghorbanchian et al. (2021), Eq. 6 of the review of Majhi, Perc, and Ghosh (2022, arXiv 2203.06601 §III), which is the source at hand. Here R_nodes is the order parameter of the node phases and R_links that of the projection of the link phases on the nodes. The review takes frequencies of mean Omega and standard deviation 1 / tau, with Omega = 2 and tau = 1 for its Fig. 2, on complexes of 500 nodes. The mean 0, the size, and the random graph are choices of this system.

The Nodes-Links-Triangles system, model='NLT', Eq. 8 of the review, scales the coupling of the link phases through the nodes by R_nodes R_up and through the triangles by R_links, where R_up is the order parameter of the projection of the link phases on the triangles. Every projection then couples to another, so all three order parameters jump together at one explosive transition, while in the Nodes-Links system R_up moves on its own.

Parameters:

Name	Type	Description	Default
n_nodes	int	The number of nodes of the random graph.	30
density	float	The probability of every edge of the random graph.	0.4
coupling	float	The coupling K.	2.0
width	float	The standard deviation of the natural frequencies of the nodes and the links.	1.0
seed	int	The seed of the graph, the frequencies, and the initial phases.	0
model	str	The coupling, 'NL' for the Nodes-Links system or 'NLT' for the Nodes-Links-Triangles system.	'NL'

Returns:

Name	Type	Description
system	System	The system, with a group triangles of unit weights on the faces for the Nodes-Links-Triangles system, which only gives the order parameter of the projection on the triangles.

Raises:

Type	Description
ValidationError	The model is neither 'NL' nor 'NLT'.

Examples:

>>> from skcomplex.systems.physics.synchronization import topological_kuramoto
>>> [element.size_ for element in topological_kuramoto(6, density=0.8).prepare().elements]
[6, 13]

Source code in src/skcomplex/systems/physics/synchronization/_topological.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'simplicial complex', 'cell complex', 'explosive'),
    references=('Majhi, Perc, and Ghosh 2022', 'Ghorbanchian et al. 2021'),
    defaults={'until': 40.0, 'dt': 0.2},
    observables=_topological_observables,
)
def topological_kuramoto(
    n_nodes: int = 30,
    *,
    density: float = 0.4,
    coupling: float = 2.0,
    width: float = 1.0,
    seed: int = 0,
    model: str = 'NL',
) -> System:
    """Return phases on the nodes and links of a clique complex, coupled through its incidence by order parameters.

    Node phases follow `theta' = omega + K R_links codifferential(sin(d theta))`. Link phases follow `phi' = omega + K
    R_nodes d(sin(codifferential(phi))) - K codifferential(sin(d phi))`. With the signs of `codifferential`, these are
    the Nodes-Links system of Ghorbanchian et al. (2021), Eq. 6 of the review of Majhi, Perc, and Ghosh (2022,
    arXiv 2203.06601 §III), which is the source at hand. Here `R_nodes` is the order parameter of the node phases and
    `R_links` that of the projection of the link phases on the nodes. The review takes frequencies of mean `Omega` and
    standard deviation `1 / tau`, with `Omega = 2` and `tau = 1` for its Fig. 2, on complexes of 500 nodes. The mean 0,
    the size, and the random graph are choices of this system.

    The Nodes-Links-Triangles system, `model='NLT'`, Eq. 8 of the review, scales the coupling of the link phases through
    the nodes by `R_nodes R_up` and through the triangles by `R_links`, where `R_up` is the order parameter of the
    projection of the link phases on the triangles. Every projection then couples to another, so all three order
    parameters jump together at one explosive transition, while in the Nodes-Links system `R_up` moves on its own.

    Args:
        n_nodes:
            The number of nodes of the random graph.

        density:
            The probability of every edge of the random graph.

        coupling:
            The coupling `K`.

        width:
            The standard deviation of the natural frequencies of the nodes and the links.

        seed:
            The seed of the graph, the frequencies, and the initial phases.

        model:
            The coupling, `'NL'` for the Nodes-Links system or `'NLT'` for the Nodes-Links-Triangles system.

    Returns:
        system:
            The system, with a group `triangles` of unit weights on the faces for the Nodes-Links-Triangles system,
            which only gives the order parameter of the projection on the triangles.

    Raises:
        ValidationError:
            The model is neither `'NL'` nor `'NLT'`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import topological_kuramoto
        >>> [element.size_ for element in topological_kuramoto(6, density=0.8).prepare().elements]
        [6, 13]
    """
    if model not in MODELS:
        error_msg = f'Unknown model {model!r}. Use one of {", ".join(MODELS)}.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(seed)
    graph = nx.gnp_random_graph(n_nodes, density, seed=seed)
    space = CellComplexSpace(NetworkSpace(graph), faces='cliques').prepare()
    n_links = space.incidences_[0].shape[0]
    nodes = GroupElement(
        'nodes',
        element=PhaseOscillatorElement,
        space=space,
        on='sites',
        values={'theta': rng.uniform(-np.pi, np.pi, n_nodes), 'omega': rng.normal(0.0, width, n_nodes)},
    )
    links = GroupElement(
        'links',
        element=PhaseOscillatorElement,
        space=space,
        on='edges',
        values={'theta': rng.uniform(-np.pi, np.pi, n_links), 'omega': rng.normal(0.0, width, n_links)},
        visual=VisualSpec('cell', {'state': 'theta'}),
    )
    interactions = [
        RateInteraction(
            'node_coupling',
            target='nodes',
            state='theta',
            sources=['links'],
            rate=partial(_node_rate, coupling),
        ),
    ]
    if model == 'NL':
        link_coupling = RateInteraction(
            'link_coupling',
            target='links',
            state='theta',
            sources=['nodes'],
            rate=partial(_link_rate, coupling),
        )
        return System(elements=[nodes, links], interactions=[*interactions, link_coupling])
    triangles = GroupElement('triangles', element=LinkElement, space=space, on='hyperedges')
    link_coupling = RateInteraction(
        'link_coupling',
        target='links',
        state='theta',
        sources=['nodes', 'triangles'],
        rate=partial(_triangle_link_rate, coupling),
    )
    return System(elements=[nodes, links, triangles], interactions=[*interactions, link_coupling])

winfree(n_oscillators=50, *, K=0.5, spread=0.1, seed=0)

Return the Winfree system, theta_i' = omega_i + (K / N) sum_j P(theta_j) R(theta_i), with the sum over all j.

The influence P(theta) = 1 + cos(theta) and the sensitivity R(theta) = -sin(theta) make every oscillator respond to the mean influence of the population, which every oscillator reads as a group mean. The natural frequencies are spread evenly over [1 - spread, 1 + spread]. Weak coupling leaves the oscillators incoherent, moderate coupling locks them to one frequency, and strong coupling stops them, which is oscillator death.

Parameters:

Name	Type	Description	Default
n_oscillators	int	The number of oscillators.	50
K	float	The coupling strength.	0.5
spread	float	The half width of the interval of natural frequencies around 1.	0.1
seed	int	The seed of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases oscillators.theta and the order parameter order.

Examples:

>>> from skcomplex.systems.physics.synchronization import winfree
>>> len(winfree(5).interactions)
1

Source code in src/skcomplex/systems/physics/synchronization/_pulses.py

@register_system(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'oscillator death'),
    references=('Winfree 1967', 'Ariaratnam and Strogatz 2001', 'Bick et al. 2020'),
    defaults={'until': 50.0, 'dt': 0.5},
    observables=lambda system: order_parameter('order', 'oscillators'),
)
def winfree(n_oscillators: int = 50, *, K: float = 0.5, spread: float = 0.1, seed: int = 0) -> System:
    """Return the Winfree system, `theta_i' = omega_i + (K / N) sum_j P(theta_j) R(theta_i)`, with the sum over all `j`.

    The influence `P(theta) = 1 + cos(theta)` and the sensitivity `R(theta) = -sin(theta)` make every oscillator respond
    to the mean influence of the population, which every oscillator reads as a group mean. The natural frequencies are
    spread evenly over `[1 - spread, 1 + spread]`. Weak coupling leaves the oscillators incoherent, moderate coupling
    locks them to one frequency, and strong coupling stops them, which is oscillator death.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        spread:
            The half width of the interval of natural frequencies around 1.

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The system, with the phases `oscillators.theta` and the order parameter `order`.

    Examples:
        >>> from skcomplex.systems.physics.synchronization import winfree
        >>> len(winfree(5).interactions)
        1
    """
    oscillators = GroupElement(
        'oscillators',
        element=PhaseOscillatorElement,
        size=n_oscillators,
        values={
            'theta': np.random.default_rng(seed).uniform(0, 2 * np.pi, n_oscillators),
            'omega': np.linspace(1 - spread, 1 + spread, n_oscillators),
        },
    )
    return System(
        elements=[oscillators],
        interactions=[WinfreeMeanFieldInteraction('coupling', target='oscillators', K=K)],
    )