# Kuramoto in Stan (PyStan)

tl;dr: Project on github: https://github.com/laszukdawid/pystan-kuramoto

Stan is a programming language focused on probabilistic computations. Although it’s a rather recent language it’s been nicely received in data science/Bayesian community for its focus on designing model, rather than programming and getting stuck with computational details. Depending what is your background you might have heard about it either from Facebook and its Prophet project or as a native implementation for Hamiltonian Monte Carlo (HMC) and its optimised variation – No-U-Turn Sampler for HMC (NUTS).

For ease of including models in other programmes there are some interfaces/wrappers available, including RStan and PyStan.

Stan is not the easiest language to go through. Currently there are about 600 pages of documentation and then separate “documentations” for wrappers, which for PyStan isn’t very helpful. Obviously there’s no need for reading all of it, but it took me a while to actually understand what goes where an why. The reward, however, is very satisfying.

Since I’ve written a bit about Kuramoto models on this blog, it’s consistent if I share its implementation in Stan as well. Pystan-kuramoto project uses PyStan, but the actual Stan code is platform independent.

Currently there are two implementations (couldn’t come up with better names):

• All-to-all, where Kuramoto model fit is performed to phase vector $\vec{\Phi}$ with distinct oscillators, i.e. $\vec{\Phi}_{N}(t) = \{\phi_1(t), \phi_2(t), \dots, \phi_N(t)\}$.
• All-to-one, where the model fits superposition (sum) of oscillators to phase time series $\Phi_{N}(t) = \sum_{n=1}^{N} \phi_n(t)$.

In all honesty, this seems to be rather efficient. Optimisation is performed using penalized maximum likelihood estimation with optimization (L-BFGS). Before using it I wasn’t very knowledgeable in the algorithm, but now I’m simply amazed with its speed and accuracy. Thumbs up.

# Spellchecker in VIM

Only recently I’ve started using VIM as my default editor and… I’m amazed.

I’m not advocating here for its superiority over any other text editor. You might be into design or simply wish for a graphical interface. I’m into not-using-mouse too much, so it actually fits quite well my purposes. Having said that, it took me a while to discover it with more commands than saving and exiting.

One of the features that from the outsider point of view, i.e. myself a month ago, that it couldn’t possibly ever have is spellchecker. Not that I actively thought that vim doesn’t have a spellchecker, but rather that passively ignoring the idea of suggesting corrections in a terminal. I was wrong. Vim has spellchecker and it’s pretty good.

The magical command is:
:setlocal spell spelllang=XX
where for English XX is one of {en, en_au, en_ca, en_gb, en_nz, en_us}. If you have a problem with setting up British, like I had, try command:
:setlocal spell! spelllang=en_gb
i.e. with extra exclamation mark.

Actually I suggest mapping command to some keyboard key, so when you press it turn on/off highlighting mode. To do this edit your .vimrc file (Ubuntu: ~/.vimrc) and add following line:
:map :setlocal spell! spelllang=en_gb
which means that whenever you press key F7 it will toggle command.

Once you have this set and someone introduces mistakes into your text/code, you can simply navigate between these using ]s and [s for next and previous, respectively. Obviously you can also just navigate using keys or hjkl, but who doesn’t want to feel like hacker?

Assuming someone has played with your text and make it less perfect, you have two options:

z=   acceptance of the mistake and listening for looking up some suggestions, or
zw   teach your dictionary this new fabulous word, or
zg   denial, blissful ignorance, snooze button.

If you change your mind with your last command you can undo it using zuw or zug. I suggest quick read about these commands in Quick Start from spell documentation or :help spell.

Now! The default setting for spellchecking performance are against large files/projects, like PhD thesis. If you notice that after certain line spellchecker doesn’t inform you about the obvious mistake then you should use this StackOverflow advice and increase your syn sync parameters. In short, but I suggest reading the SO question and answers (and comments), you can increase minlines and maxlines in (Ubuntu) ~/.vim/after/syntax/tex.vim with:

syn sync maxlines=20000
syn sync minlines=500


# Easy clipboard in bash

Those who know me, know that I’m not a fan of using mouse. Whenever possible I try to avoid using it. Although, after years of experience, I am fluent in keyboarding there are still some tasks that I need mouse. Well, until quite recent.

The number one of annoying tasks is copying things from terminal to clipboard. Depending what exactly I needed to do with it I’d either select something with mouse and then copy, or stream output to file and then select it within editor. Some people found it weird, but yes, often copying through editor is much faster.

Revolution came with xclip. I’m not yet fluent in doing everything with it, but even with limited experience that I have it feels like a superpower. This program allows to copy into/from X clipboard.

Most unix distro have it installed by default or at least in package manager. In case of Ubuntu one can install it with:
$sudo apt-get update$ sudo apt-get install xclip 

Examples:
Copy current directory into clipboard:
$pwd | xclip Paste whatever you have in clipboard: $ xclip -o

Copy to global clipboard so you can use in any other program:
$echo$PATH | xclip -selection clipboard

If this is too long to type, one can always use aliases, i.e. setting shorter name for long command. Either type directly into terminal (but that is only for current terminal), or update your ~/.bashrc file with:
alias "c=xclip" alias "v=xclip -o" alias "cs=xclip -selection clipboard" alias "vs=xclip -o -selection clipboard" 

Then copying content of a file to a clipboard:
$cat file.txt | cs … and CTRL+V wherever one wishes. Yes, this is awesome! # Kuramoto in Python Code for Kuramoto in Python is available here or from code subpage. Explanation on how to use it is on the bottom of this post. ## Tiny introduction Kuramoto[1, 2] is probably one of the most popular and successful models for coupled oscillators. There is plenty of information about it, but in brief summary it models oscillators’ phases to be dependent on scaled phase differences for all pairs of oscillator. A while ago I have actually wrote on a specific method that I used to find optimal parameters for a Kuramoto model given some observation. (See Bayesian Dynamic Inference here.) Mathematically speaking this model is defined by a set of coupled ordinary differential equations (ODEs). Given N oscillators, dynamics for each oscillator’s phase $\phi_i$ is defined as $i\in N$ is defined by $\dot\phi_i = \omega_i + \sum_{j=0}^N k_{ij} \sin(\phi_i-\phi_j)$, where the summation is over all others oscillators. Note that one can leave coupling with itself, because in such case $\sin(\phi_i - \phi_i) = 0$, thus it doesn’t input anything into the final solution. Quick note that Kuramoto coupling model can be used to reference model with coupling function that can be represented as a sum of harmonic functions. For example, second order means that we are including both $k \sin(\Delta\phi)$ and $k^{(2)}\sin(2\Delta\phi)$. In general Kuramoto model of order $M$ would describe $\dot\phi_i = \omega_i + \sum_{m=0}^{M} \sum_{j=0}^N k_{ij}^{(m)} \sin(m (\phi_i-\phi_j))$. ## Examples Below are two examples. Both use the same core values with the exception that for the second system the coupling function defined by two harmonic terms. This naturally changes dynamics, but it shouldn’t change much average value, which is close to intrinsic frequency. Exact values for the first experiment are presented in table below. Each row is a different oscillator with initial phase Φ0 and intrinsic frequency Ω. The following columns denote coupling strength between respective pairs of oscillators. Φ0 ω k.1 k.2 k.3 1 0 28 0.2 1.1 2 π 19 0.5 -0.7 3 0 11 0.3 0.9 Phase dynamics, i.e. time derivative of obtained phases, are presented in Fig. 1. One can see that all plots are centred on intrinsic frequency with some modulations. Fig. 1. Phase dynamics, i.e. time derivative, in simple Kuramoto model. Table below shows values used in the second simulation. Extra 3 columns denote scaling values used in second harmonic. Φ0 ω k(1).1 k(1).2 k(1).3 k(2).1 k(2).2 k(2).3 1 0 28 0.2 1.1 -0.5 0.2 2 π 19 0.5 -0.7 -0.4 1.0 3 0 11 0.3 0.9 0.8 0.8 Again phase dynamics have been presented in a form of plot (Fig. 2). I think the difference is clear. Despite having similar mean values (approximately equal to intrinsic frequency), their modulations have change. Not only the frequency content of these modulations has changed, but also their amplitude. Fig. 2. Phase dynamics, i.e. time derivative, in Kuramoto model with included second harmonic. ## Using code Except for downloading code from either github or code subpage, one is expected to have SciPy module. Kuramoto uses it to solve differential equations. Script below shows an example of how one can use the Kuramoto module. When you run this, make sure that kuramoto.py is either in your path. Note also that most of the code is to just defining initial parameters and plotting the results. Actual execution of the module are two lines. Fig. 1 is the expected output. import numpy as np import pylab as plt from kuramoto import Kuramoto # Defining time array t0, t1, dt = 0, 40, 0.01 T = np.arange(t0, t1, dt) # Y0, W, K are initial phase, intrinsic freq and # coupling K matrix respectively Y0 = np.array([0, np.pi, 0]) W = np.array([28, 19, 11]) K1 = np.array([[ 0, 0.2, 1.1], [0.5, 0, -0.7], [0.3, 0.9, 0]]) K2 = np.array([[ 0, -0.5, 0.2], [-0.4, 0, 1.0], [ 0.8, 0.8, 0]]) K = np.dstack((K1, K2)).T # Passing parameters as a dictionary init_params = {'W':W, 'K':K, 'Y0':Y0} # Running Kuramoto model kuramoto = Kuramoto(init_params) odePhi = kuramoto.solve(T) # Computing phase dynamics phaseDynamics = np.diff(odePhi)/dt # Plotting response nOsc = len(W) for osc in range(nOsc): plt.subplot(nOsc, 1, 1+osc) plt.plot(T[:-1], phaseDynamics[osc]) plt.ylabel("$\dot\phi_{%i}\$" %(osc+1))
plt.show()


[1] Y. Kuramoto, “Chemical Oscillations, Waves, and Turbulence,” 1984, doi: 10.1007/978-3-642-69689-3.
[2] Steven H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” 2000, doi: 10.1016/S0167-2789(00)00094-4.

# Update: Particle Swarm Optimisation in Python

It came to my attention that my PSO for Python is actually quite popular. Only after few people contacted me I’ve noticed that the public version was not the same that I’ve been using. It wasn’t bad, but definitely not as good. Thus, obviously, I’ve updated the version and cleaned it a bit.

Update programme is available from my github or from Code subpage.

What’s the difference? There are few.
– Initial values, unless provided, are psuedo-random generated using Halton sequence. This prevents from artificial clustering.
– Perturbation in form of a Gaussian noise should mitigate false local minima by forcing particles to search surrounding area.
– Added max_repetition threshold, which states the maximum number of obtaining the same optimum value. Upon reaching threshold program finishes.
– General improvement in performance.
– Improved usage documentation within the file.
– Program is now compatible with Python3.

Feel free to request any features.

There is an idea of adding progressive save, which would quit, resume and modify parameters at any point of computation.

# Common dtype for NumPy arrays

Recent challenge I had was to convert two given numpy arrays such it’ll be possible to have them in common type without losing information. One can follow this link for stackoverflow entry.

Initially I thought about comparing dtypes, because in Python 2 it is allowed to do something like:

>>> np.float32 < np.float64
True
>>> np.int32 < np.float32
True
>>> np.int64 > np.float16
False


… which kind of makes sense(?). Well, except that int64 vs float16 which looks suspicious. It turns out that these are type objects and Python is simply comparing their locations in memory[Citation needed] and that is obviously not reliable. Actually, in Python 3 such comparison is forbidden and it fails.

As the answer to my stackoverflow question suggests one could try to use dtype.kind and dtype.itemsize to create own ordering. This is not difficult, but it should contain all types, such as (unsigned) ints, floats, bools, complex…

Fortunately, for my purposes, there is NumPy’s method which does what I want, i.e. numpy.find_common_type. It determines common type following standard coercion rules. With a help of this function my common conversion looks like:

import numpy as np

def common_dtype(x, y):
dtype = np.find_common_type([x.dtype, y.dtype], [])
if x.dtype != dtype: x = x.astype(dtype)
if y.dtype != dtype: y = y.astype(dtype)

return x, y


What to expect from such function?
float32 for float16 and float32
float32 for int16 and float32
int64 for int32 and uint16
float64 for int32 and float16

Behaves as it should.

# Halton sequence in Python

Sometimes when we ask for random we don’t actually mean random by just random. Yes, pseudo-random.

Consider unitary distribution with ranges 0 and 1. Say you want to draw 5 samples. Selecting them at random would mean that we might end up with set of {0, 0.1, 0.02, 0.09, 0.01} or {0.11, 0.99, 0.09, 0.91, 0.01}. Yes, these values don’t seem very random, but that’s the thing about randomness, that it randomly can seem to not be random.

Depending on the purpose of our selection, these values might be just OK. After all, they came from that distribution. However, if our goal is to reconstruct the distribution, or extract information about with limited number of samples, it is often better to draw those samples in pseudo-random way. For example, in accordance to van der Corput sequences for 1D distributions or its generalized version Halton sequence.

The best practice for sampling N dimensional distribution is to use different prime numbers for each dimension. For example, when I need to sample a 5 dimensional unitary distribution, or search space, I will use bases of (5, 7, 11, 13, 17). This is to prevent periodic visits of the same position.

In case you are wondering what’s the difference between actual random and pseudo-random, here is a gist:

Both are good, but the actual random can produce many empty holes. What we like to have is a fair representation of all areas of our search space.

Thus, without further ado, here are some code snippets.

This is a definition of my prime generating generator:

def next_prime():
def is_prime(num):
"Checks if num is a prime value"
for i in range(2,int(num**0.5)+1):
if(num % i)==0: return False
return True

prime = 3
while(1):
if is_prime(prime):
yield prime
prime += 2


As for Halton sequence, as mentioned before it uses van der Corput sequence. Again, here is the definition:

def vdc(n, base=2):
vdc, denom = 0, 1
while n:
denom *= base
n, remainder = divmod(n, base)
vdc += remainder/float(denom)
return vdc


And finally, definition for the Halton sequence:

def halton_sequence(size, dim):
seq = []
primeGen = next_prime()
next(primeGen)
for d in range(dim):
base = next(primeGen)
seq.append([vdc(i, base) for i in range(size)])
return seq


To use all of this simply call halton_sequence(size, dim). These variables refer to the number of size of sample poll and the dimension of your problem. So if one wants to sample 3 dimensional space with 10 samples each it would be called as below. (Notice: first dimension has prime value 5, then it’s 7, 11, and following prime values.)

>>> halton_sequence(10, 3)
[
[0, 0.2, 0.4, 0.6, 0.8, 0.04, 0.24000000000000002, 0.44, 0.64, 0.8400000000000001],
[0, 0.14285714285714285, 0.2857142857142857, 0.42857142857142855, 0.5714285714285714, 0.7142857142857143, 0.8571428571428571, 0.02040816326530612, 0.16326530612244897, 0.30612244897959184],
[0, 0.09090909090909091, 0.18181818181818182, 0.2727272727272727, 0.36363636363636365, 0.45454545454545453, 0.5454545454545454, 0.6363636363636364, 0.7272727272727273, 0.8181818181818182]
]