Source code for netrd.dynamics.kuramoto

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'`. For more information on the Kuramoto model, see the review essay 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. """ 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) # adjust phases TS = TS % (2 * np.pi) self.results["internal_frequencies"] = omega self.results["ground_truth"] = G self.results["TS"] = TS return TS