deterministic kuramoto models

kuramoto.py

import numpy as np
from types import *
sin = np.sin

#    Deterministic Kuramoto model like functions
#    Copyright (C) 2013  \0/\o/\0/\o/\0/

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU Affero General Public License as
#    published by the Free Software Foundation, either version 3 of the
#    License, or (at your option) any later version.

#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU Affero General Public License for more details.

#    You should have received a copy of the GNU Affero General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

def kuramoto(theta, t, freqs, A, K = False, Normed = False):

    ''' This is a general sinusoidal coupled phase system with strength times coupling K adjacency matrix A , it's also uses  the diagonal terms A_ii to store the speed field as it account for the components.
        In the case that K argument is given (not false), than it will try to evaluate  the function  according to K type. If K is a number than it will be assigned as coupling constant K*A. If K is a function it will try to calculate the Ks for the chain case. TODO-If it has a tuple with length > 1 parameters in position 1 of the tuple, than it will try to use the last element of theta and the parameters to update the K value.
        Normed parameter may be used for mean-field studies (i.e. in the limit of N -> infnity number of oscillators with as many as possible connections), it will scale the coupling with (1/N).
    
        Topologies examples:

        Unidimensional chain with uncoupled border: A = np.diag(np.ones(N-1),1)+ np.diag(np.ones(N-1),-1);

        Unidimensional Ring: A = np.diag(np.ones(N-1),1)+ np.diag(np.ones(N-1),-1); A[0,N-1] = 1; A[N-1, 0] = 1;

    '''         

    N = len(A);
    freqs =  np.diag(freqs);
    if K != False:
        tipo = type(K);
        if (tipo == int) | (tipo == float):
            A = K*A;
        if tipo == tuple:
            l = len(K);
            #if l > 1:
            # TODO - caso com acoplamento plastico 
        if type(K) is FunctionType:
            Ks = K(freqs);           
            A = Ks*A;

    if Normed == True:
        A = (1/N)*(A - np.diag(np.diag(A))) + freqs
    else:
        A = (A - np.diag(np.diag(A))) + freqs;
    A[0,0] =  A[0,0] + (A[0,1:]*np.sin(theta[1:]-theta[0])).sum();
 
    for i in range(1,N):
        A[i,i] =  A[i,i] + (A[i,0:i]*np.sin(theta[0:i]-theta[i])).sum(); 
        A[i,i] =  A[i,i] + (A[i,i+1:N]*np.sin(theta[i+1:]-theta[i])).sum();
    return np.diag(A);

def coherence_r(freqs,A):
    ''' Computes r(t) e^{1j Psi(t)} = (1/N)*\sum_{i}{e^{1j \theta_{i}(t)}} '''

    N = len(A);
    theta = np.random.uniform(- np.pi, np.pi,N);
    t = np.linspace(0.0,40.0,800);
    r = np.zeros(800);
    psi = np.zeros(800);
    y = integ.odeint(kuramoto,theta,t,(freqs,A));    
    for i in range(800):
        a = np.exp(1j*y[i,:]).sum()/float(N);      
        psi[i] = np.arctan(a.imag/a.real);
        a = a*np.exp(-1j*psi[i]);
        r[i] = abs(a);
    return (r,psi);

def r_k(freqs, A):
    l = len(freqs);
    r_m = np.zeros(50);
    k = np.zeros(50);
    c = 0;
    psi_m = np.zeros(50);
    ran = np.linspace(0,6.0,50);
    for i in ran:
        B = i*A;
        (r,psi) = coherence_r(freqs,B);
        r_m[c] = np.mean(r[400:801]);
        psi_m[c] = np.mean(psi[400:801]);
        k[c] = i;
        c = c + 1;
    return (k,r_m,psi_m);