search

fabridamicelli / kuramoto

Python implementation of the Kuramoto model
GNU General Public License v3.0
123 stars 25 forks source link

readme

Downloads

Like the package? Don't forget to give it a GitHub ⭐ to help others find and trust it!

kuramoto

Python implementation of the Kuramoto model.

Install

pip install kuramoto

Usage

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import seaborn as sns

from kuramoto import Kuramoto, plot_phase_coherence, plot_activity

sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.6)

# Interactions are represented as an adjacency matrix _A_, a 2D numpy ndarray.
# Instantiate a random graph and transform into an adjacency matrix
graph_nx = nx.erdos_renyi_graph(n=100, p=1) # p=1 -> all-to-all connectivity
adj_mat = nx.to_numpy_array(graph_nx)

# Instantiate model with parameters
model = Kuramoto(coupling=3, dt=0.01, T=10, n_nodes=len(adj_mat))

# Run simulation - output is time series for all nodes (node vs time)
activity = model.run(adj_mat=adj_mat)

# Plot all the time series
plot_activity(activity)

png

# Plot evolution of global order parameter R_t
plot_phase_coherence(act_mat)

png

# Plot oscillators in complex plane at times t = 0, 250, 500
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(15, 5),
                         subplot_kw={
                             "ylim": (-1.1, 1.1),
                             "xlim": (-1.1, 1.1),
                             "xlabel": r'$\cos(\theta)$',
                             "ylabel": r'$\sin(\theta)$',
                         })

times = [0, 200, 500]
for ax, time in zip(axes, times):
    ax.plot(np.cos(act_mat[:, time]),
            np.sin(act_mat[:, time]),
            'o',
            markersize=10)
    ax.set_title(f'Time = {time}')

png

As a sanity check, let's look at the phase transition of the global order parameter (Rt) as a function of coupling (K) (find critical coupling Kc) and compare with numerical results already published by English, 2008 (see Ref.) – we will match those model parameters.

# Instantiate a random graph and transform into an adjacency matrix
n_nodes = 500
graph_nx = nx.erdos_renyi_graph(n=n_nodes, p=1) # p=1 -> all-to-all connectivity
graph = nx.to_numpy_array(graph_nx)

# Run model with different coupling (K) parameters
coupling_vals = np.linspace(0, 0.6, 100)
runs = []
for coupling in coupling_vals:
    model = Kuramoto(coupling=coupling, dt=0.1, T=500, n_nodes=n_nodes)
    model.natfreqs = np.random.normal(1, 0.1, size=n_nodes)  # reset natural frequencies
    act_mat = model.run(adj_mat=graph)
    runs.append(act_mat)

# Check that natural frequencies are correct (we need them for prediction of Kc)
plt.figure()
plt.hist(model.natfreqs)
plt.xlabel('natural frequency')
plt.ylabel('count')

png

# Plot all time series for all coupling values (color coded)
runs_array = np.array(runs)

plt.figure()
for i, coupling in enumerate(coupling_vals):
    plt.plot(
        [model.phase_coherence(vec)
         for vec in runs_array[i, ::].T],
        c=str(1-coupling),  # higher -> darker   
    )
plt.ylabel(r'order parameter ($R_t$)')
plt.xlabel('time')

png

# Plot final Rt for each coupling value
plt.figure()
for i, coupling in enumerate(coupling_vals):
    r_mean = np.mean([model.phase_coherence(vec)
                      for vec in runs_array[i, :, -1000:].T]) # mean over last 1000 steps
    plt.scatter(coupling, r_mean, c='steelblue', s=20, alpha=0.7)

# Predicted Kc – analytical result (from paper)
Kc = np.sqrt(8 / np.pi) * np.std(model.natfreqs) # analytical result (from paper)
plt.vlines(Kc, 0, 1, linestyles='--', color='orange', label='analytical prediction')

plt.legend()
plt.grid(linestyle='--', alpha=0.8)
plt.ylabel('order parameter (R)')
plt.xlabel('coupling (K)')
sns.despine()

png

Kuramoto model 101

jpg

where K is the coupling parameter and Mi is the number of oscillators interacting with oscillator i. A is the adjacency matrix enconding the interactions - typically binary and undirected (symmetric), such that if node i interacts with node j, Aij = 1, otherwise 0. The basic idea is that, given two oscillators, the one running ahead is encouraged to slow down while the one running behind to accelerate.

In particular, the classical set up has M = N, since the interactions are all-to-all (i.e., a complete graph). Otherwise, Mi is the degree of node i.

Kuramoto model 201

A couple of facts in order to gain intuition about the model's behaviour:

For more and better details, this talk by the great Steven Strogatz is a nice primer.

Requirements

Tests

Run tests with

make test

Citing

If you find this package useful for a publication, then please use the following BibTeX to cite it:

@misc{kuramoto,
  author = {Damicelli, Fabrizio},
  title = {Python implementation of the Kuramoto model},
  year = {2019},
  publisher = {GitHub},
  journal = {GitHub repository},
  howpublished = {\url{https://github.com/fabridamicelli/kuramoto}},
}

References & links