"""
Functions to identify states in coupled oscillators
"""
import warnings
import numpy as np
from numpy.linalg import norm
__all__ = [
"identify_state",
"identify_k_clusters",
"identify_q_twisted",
"order_parameter",
]
[docs]def identify_state(thetas, atol=1e-3):
"""
Identify the synchronization state.
Parameters
----------
thetas : array-like, shape (n_phases,)
Phases of the oscillators.
atol : float, optional
The absolute tolerance used for comparison with synchronization parameters.
Default is 1e-3.
Returns
-------
state : str
A string representing the identified synchronization state:
- "q-twisted" for q-twisted synchronization, where q is the winding number.
- "splay" for splay synchronization.
- "sync" for near-synchronised states with R1 ≈ 1.
- "2-cluster" for 2-cluster synchronization.
- "3-cluster" for 3-cluster synchronization.
- "other" for other or unsynchronized states.
Notes
-----
A perfectly synchronised state (all phases equal) is returned as "0-twisted",
not "sync". The "sync" label is reserved for noisy near-synchronised states
where the order parameter R1 ≈ 1 but the phase differences are not exactly
uniform.
A splay state where phases increase monotonically with node index is
indistinguishable from a 1-twisted state and will be returned as "1-twisted".
The "splay" label is only returned when phases are evenly distributed on the
circle but not ordered by node index.
See Also
--------
identify_q_twisted : Detect a q-twisted state and compute its winding number.
identify_k_clusters : Detect a k-cluster state and return cluster sizes.
order_parameter : Compute the Daido order parameter.
Examples
--------
>>> import numpy as np
>>> import hypersynchronization as hs
>>> hs.identify_state(np.zeros(10))
'0-twisted'
>>> theta = hs.generate_q_twisted_state(10, q=1, noise=0)
>>> hs.identify_state(theta)
'1-twisted'
>>> theta = np.array([0.0] * 5 + [np.pi] * 5)
>>> hs.identify_state(theta)
'2-cluster'
"""
R1 = order_parameter(thetas, order=1)
# R2 = order_parameter(thetas, order=2)
# R3 = order_parameter(thetas, order=3)
diff = np.diff(thetas, append=thetas[0]) % (2 * np.pi)
is_diff_zero = np.isclose(diff, 0, atol=atol) + np.isclose(
diff, 2 * np.pi, atol=atol
)
q, is_twisted = identify_q_twisted(thetas)
sorted_thetas = np.sort(thetas, axis=0) # sort along node axis
q_sorted, is_splay = identify_q_twisted(sorted_thetas)
try:
is_2clust, sizes2 = identify_k_clusters(thetas, k=2, atol=1e-2)
except Exception as err:
is_2clust = False
sizes2 = []
warnings.warn(str(err))
try:
is_3clust, sizes3 = identify_k_clusters(thetas, k=3, atol=1e-2)
except Exception as err:
is_3clust = False
sizes3 = []
warnings.warn(str(err))
if is_twisted:
return f"{q}-twisted"
elif is_splay and q_sorted == 1:
return "splay"
elif np.isclose(R1, 1, atol=atol) and np.all(is_diff_zero):
return "sync"
elif is_2clust:
return "2-cluster"
elif is_3clust:
return "3-cluster"
else:
return "other"
[docs]def identify_k_clusters(thetas, k, atol=1e-2):
"""
Check if k-cluster state.
A k-cluster state has k evenly spaced clusters on the circle.
Parameters
----------
thetas : array-like, shape (n_phases,)
Phases of the oscillators.
k : int
Number of clusters.
atol : float, optional
The absolute tolerance used for comparison with synchronization
parameters. Default is 1e-2.
Returns
-------
is_k_clusters : bool
True if the state is a k-cluster. False otherwise.
sizes : list of float
The relative sizes of each cluster.
See Also
--------
identify_state : Identify the overall synchronization state.
generate_k_clusters : Generate a k-cluster initial condition.
Examples
--------
>>> import numpy as np
>>> import hypersynchronization as hs
>>> theta = np.array([0.0] * 5 + [np.pi] * 5)
>>> is_2clust, sizes = hs.identify_k_clusters(theta, k=2)
>>> is_2clust, sizes
(True, [0.5, 0.5])
"""
n_clust = k
dist = 2 * np.pi / n_clust
N = len(thetas)
psi = thetas % (2 * np.pi)
psi = np.sort(psi)
diff = np.diff(psi)
# identify jumps larger than half the expected inter-cluster distance
idcs = np.where(diff >= 0.45 * dist)[0]
clusters = []
n_changes = len(idcs)
# partition phases in clusters separated by the jumps
for i in range(n_changes + 1):
start = idcs[i - 1] + 1 if i > 0 else None
end = idcs[i] + 1 if i < n_changes else None
clusters.append(psi[start:end])
if len(clusters) < k:
return False, []
is_k_clusters = True # changed below if False
sizes = [0] * k
# check if each cluster is compact enough
for i in range(n_changes + 1):
if np.mean(np.diff(clusters[i])) > atol: # cluster is not compact
is_k_clusters = False
# check if the clusters have the right distance between them
for i in range(n_changes):
dist_ij = abs(np.mean(clusters[i]) - np.mean(clusters[i + 1]))
if abs(dist_ij - dist) > atol:
is_k_clusters = False # clusters have wrong distance between them
if n_clust == len(clusters):
sizes = [len(cluster) / N for cluster in clusters]
elif n_clust == len(clusters) - 1:
sizes = [len(cluster) / N for cluster in clusters[:-1]]
sizes[0] += len(clusters[-1]) # 0th and last clusters are the same
else: # too many clusters found, can be a k-cluster state with larger k
is_k_clusters = False
# raise ValueError(f"k must be equal to or one unit below len(clusters)={len(clusters)}")
return is_k_clusters, sizes
[docs]def identify_q_twisted(thetas, atol=1e-1):
"""
Check if twisted state and identify its winding number.
The winding number indicates how many times the phase
angles wind around the unit circle.
Parameters
----------
thetas : array-like, shape (n_phases,)
Phases of the oscillators.
atol : float, optional
Tolerance for deciding whether phase differences are uniform enough
to be a twisted state. Default is 1e-1.
Returns
-------
q : int
The winding number, indicating how many times the phase
angles wind around the unit circle.
is_twisted_state : bool
True if the phase differences are close to their mean,
indicating a twisted state; False otherwise.
See Also
--------
identify_state : Identify the overall synchronization state.
generate_q_twisted_state : Generate a q-twisted initial condition.
Examples
--------
>>> import hypersynchronization as hs
>>> theta = hs.generate_q_twisted_state(10, q=2, noise=0)
>>> q, is_twisted = hs.identify_q_twisted(theta)
>>> q, is_twisted
(2, True)
"""
thetas = thetas % (2 * np.pi) # ensure it's mod 2 pi
diff = np.diff(thetas, prepend=thetas[-1])
# ensure phase diffs are in [-pi, pi]
diff = np.where(diff > np.pi, diff - 2 * np.pi, diff)
diff = np.where(diff < -np.pi, diff + 2 * np.pi, diff)
q = np.sum(diff)
q = round(q / (2 * np.pi)) # winding number
is_twisted_state = norm(diff - np.mean(diff)) < atol
return q, is_twisted_state
[docs]def order_parameter(thetas, order=1, complex=False, axis=0):
"""
Calculate the generalised Daido order parameter of order `order`.
The order parameter is a measure of how the phase angles are aligned.
It can be used to quantify the level of synchronization or coherence in a system.
Parameters
----------
thetas : array-like, shape (n_oscillators,) or (n_oscillators, n_times)
Phases of the oscillators over time.
order : int, optional
The order of the order parameter. Default is 1.
complex : bool, optional
If True, return the complex order parameter. If False (default),
return its magnitude.
axis : int, optional
The axis of length n_oscillators, along which the sum of the phases
is taken. Default is 0.
Returns
-------
result : float or numpy.ndarray
If `complex` is True, the complex order parameter.
If `complex` is False, the magnitude of the order parameter.
See Also
--------
identify_state : Identify the overall synchronization state.
Examples
--------
>>> import numpy as np
>>> import hypersynchronization as hs
>>> theta = np.zeros(10) # fully synchronized
>>> hs.order_parameter(theta)
1.0
>>> theta = np.linspace(0, 2 * np.pi, 10, endpoint=False) # splay state
>>> np.isclose(hs.order_parameter(theta), 0)
True
"""
N = len(thetas)
Z = np.sum(np.exp(1j * order * thetas), axis=axis) / N
if complex:
return Z
else:
return np.abs(Z)