# Source code for netrd.dynamics.kuramoto

"""
kuramoto.py
-----------
Kuramoto model of oscillators.

author: Harrison Hartle
"""

from .base import BaseDynamics
import networkx as nx
import numpy as np
import scipy.integrate as it
from ..utilities import unweighted

[docs]class Kuramoto(BaseDynamics):
"""Kuramoto model of oscillators."""

[docs]    @unweighted
def simulate(self, G, L, dt=0.01, strength=1, phases=None, freqs=None):
r"""Simulate Kuramoto model on a ground truth network.

Kuramoto oscillators model synchronization processes. At each time
step, each node adjusts its phase :math:\theta_i according to the
equation

.. math::
\theta_i = \omega_i + \frac{\lambda}{N}\sum_{j=1}^{N}\sin\left(\theta_j - \theta_i\right),

where :math:\lambda, is a coupling strength parameter and each node
has an internal frequency :math:\omega_i; the freqs function
parameter provides the option to initialize these frequencies with
user-defined values (or leave as None to randomly initialize). Each
node's initial phase :math:\theta_{i0} can be randomly initialized
(the default behavior) or set by specifying the phases parameter.

The results dictionary also stores the ground truth network as
'ground_truth' and the internal frequencies of the process as
'internal_frequencies'.

included below.

Parameters
----------

G (nx.Graph)
the input (ground-truth) graph with :math:N nodes.

L (int)
the length of the desired time series.

dt (float)
size of timestep for numerical integration.

strength (float)
coupling strength (prefactor for interaction terms).

phases (np.ndarray)
an :math:N \times 1 array of initial phases.

freqs (np.ndarray)
an :math:N \times 1 array of internal frequencies.

Returns
-------

TS (np.ndarray)
an :math:N \times L array of synthetic time series data.

Examples
--------

.. code:: python

G = nx.ring_of_cliques(4,16)
N = G.number_of_nodes()
L = int(1e4)
omega = np.random.uniform(0.95, 1.05, N)
dynamics = Kuramoto()
TS = dynamics.simulate(G, L, dt=0.01, strength=0.3, freqs=omega)

References
----------
.. [1] F. Rodrigues, T. Peron, P. Ji, J. Kurths.
The Kuramoto model in complex networks.
https://arxiv.org/abs/1511.07139

"""
A = nx.to_numpy_array(G)
N = G.number_of_nodes()

try:
if phases is not None:
assert len(phases) == N
theta_0 = phases
else:
theta_0 = 2 * np.pi * np.random.rand(N)

if freqs is not None:
assert len(freqs) == N
omega = freqs
else:
omega = np.random.uniform(0.9, 1.1, N)

except AssertionError:
raise ValueError("Initial conditions must be None or lists of length N.")

t = np.linspace(dt, L * dt, L)  # time-vector
one = np.ones(N)  # define a rate of change function

def ddt_theta(theta, t, g, strength, A):
prefactor = strength / N
first = np.outer(one, theta)
second = np.outer(theta, one)

return g + prefactor * (A * np.sin(first - second)).dot(one)

# integrate the equations of motion numerically
args = (omega, strength, A)
TS_T = it.odeint(ddt_theta, theta_0, t, args=args)

# odeint returns LxN result
# transposing yields reversed-order nodes => apply flipud.
TS = np.flipud(TS_T.T)