NetDynamics.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat May 11 08:56:41 2019

@author: nolte

D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019)
"""

# https://www.python-course.eu/networkx.php
# https://networkx.github.io/documentation/stable/tutorial.html
# https://networkx.github.io/documentation/stable/reference/functions.html

import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt
import networkx as nx
from UserFunction import linfit
import time

tstart = time.time()

plt.close('all')



Nfac = 25   # 25
N = 100      # 50
width = 0.2

# model_case 1 = Complete Graph
# model_case 2 = Barabasi
# model_case 3 = Power Law
# model_case 4 = D-dimensional Hypercube
# model_case 5 = Erdos Renyi
# model_case 6 = Random Regular
# model_case 7 = Strogatz
# model_case 8 = Hexagonal lattice
# model_case 9 = Tree
# model_case 10 == 2D square lattice

model_case = int(input('Input Model Case (1-10)'))
if model_case == 1: # Complete Graph
    facoef = 0.2
    nodecouple = nx.complete_graph(N)
elif model_case == 2: # Barabasi
    facoef = 2
    k = 3
    nodecouple = nx.barabasi_albert_graph(N, k, seed=None)
elif model_case == 3: # Power law
    facoef = 3
    k = 3
    triangle_prob = 0.3
    nodecouple = nx.powerlaw_cluster_graph(N, k, triangle_prob)        
elif model_case == 4:
    Dim = 6
    facoef = 3
    nodecouple = nx.hypercube_graph(Dim)
    N = nodecouple.number_of_nodes()
elif model_case == 5: # Erdos-Renyi
    facoef = 5
    prob = 0.1
    nodecouple = nx.erdos_renyi_graph(N, prob, seed=None, directed=False)
elif model_case == 6: # Random
    facoef = 5
    nodecouple = nx.random_regular_graph(3, N, seed=None)
elif model_case == 7: # Watts
    facoef = 7
    k = 4;              # nearest neighbors
    rewire_prob = 0.2   # rewiring probability
    nodecouple = nx.watts_strogatz_graph(N, k, rewire_prob, seed=None)
elif model_case == 8:
    facoef = 8
    rows = 4
    colm = 8
    nodecouple = nx.hexagonal_lattice_graph(rows, colm, periodic=True, with_positions=False)
    N = nodecouple.number_of_nodes()
elif model_case == 9: # k-fold tree
    facoef = 16
    k = 3  # degree
    h = 3  # height
    sm = 0
    for lp in range(h+1):
        sm = sm + k**lp
    N = sm
    nodecouple = nx.balanced_tree(k, h)
elif model_case == 10:  # square lattice
    facoef = 3
    m = 6
    n = 6
    nodecouple = nx.grid_2d_graph(m, n, periodic=True)
    N = nodecouple.number_of_nodes()

plt.figure(1)
nx.draw(nodecouple)
#nx.draw_circular(nodecouple)
#nx.draw_spring(nodecouple)
#nx.draw_spectral(nodecouple)
print('Number of nodes = ',nodecouple.number_of_nodes())
print('Number of edges = ',nodecouple.number_of_edges())
#print('Average degree = ',nx.k_nearest_neighbors(nodecouple))
        
# function: omegout, yout = coupleN(G)
def coupleN(G):

    # function: yd = flow_deriv(x_y)
    def flow_deriv(y,t0):
                
        yp = np.zeros(shape=(N,))
        ind = -1
        for omloop in G.node:
            ind = ind + 1
            temp = omega[ind]
            linksz = G.node[omloop]['numlink']
            for cloop in range(linksz):
                cindex = G.node[omloop]['link'][cloop]
                indx = G.node[cindex]['index']
                g = G.node[omloop]['coupling'][cloop]
                

                temp = temp + g*np.sin(y[indx]-y[ind])
            
            yp[ind] = temp
        
        yd = np.zeros(shape=(N,))
        for omloop in range(N):
            yd[omloop] = yp[omloop]
        
        return yd
    # end of function flow_deriv(x_y)

    mnomega = 1.0
    
    ind = -1
    for nodeloop in G.node:
        ind = ind + 1
        omega[ind] = G.node[nodeloop]['element']
    
    x_y_z = omega    
    
    # Settle-down Solve for the trajectories
    tsettle = 100
    t = np.linspace(0, tsettle, tsettle)
    x_t = integrate.odeint(flow_deriv, x_y_z, t)
    x0 = x_t[tsettle-1,0:N]
    
    t = np.linspace(0,1000,1000)
    y = integrate.odeint(flow_deriv, x0, t)
    siztmp = np.shape(y)
    sy = siztmp[0]
        
    # Fit the frequency
    m = np.zeros(shape = (N,))
    w = np.zeros(shape = (N,))
    mtmp = np.zeros(shape=(4,))
    btmp = np.zeros(shape=(4,))
    for omloop in range(N):
        
        if np.remainder(sy,4) == 0:
            mtmp[0],btmp[0] = linfit(t[0:sy//2],y[0:sy//2,omloop]);
            mtmp[1],btmp[1] = linfit(t[sy//2+1:sy],y[sy//2+1:sy,omloop]);
            mtmp[2],btmp[2] = linfit(t[sy//4+1:3*sy//4],y[sy//4+1:3*sy//4,omloop]);
            mtmp[3],btmp[3] = linfit(t,y[:,omloop]);
        else:
            sytmp = 4*np.floor(sy/4);
            mtmp[0],btmp[0] = linfit(t[0:sytmp//2],y[0:sytmp//2,omloop]);
            mtmp[1],btmp[1] = linfit(t[sytmp//2+1:sytmp],y[sytmp//2+1:sytmp,omloop]);
            mtmp[2],btmp[2] = linfit(t[sytmp//4+1:3*sytmp/4],y[sytmp//4+1:3*sytmp//4,omloop]);
            mtmp[3],btmp[3] = linfit(t[0:sytmp],y[0:sytmp,omloop]);

        
        #m[omloop] = np.median(mtmp)
        m[omloop] = np.mean(mtmp)
        
        w[omloop] = mnomega + m[omloop]
     
    omegout = m
    yout = y
    
    return omegout, yout
    # end of function: omegout, yout = coupleN(G)



Nlink = N*(N-1)//2      
omega = np.zeros(shape=(N,))
omegatemp = width*(np.random.rand(N)-1)
meanomega = np.mean(omegatemp)
omega = omegatemp - meanomega
sto = np.std(omega)



lnk = np.zeros(shape = (N,), dtype=int)
ind = -1
for loop in nodecouple.node:
    ind = ind + 1
    nodecouple.node[loop]['index'] = ind
    nodecouple.node[loop]['element'] = omega[ind]
    nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop))
    nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop)))
    lnk[ind] = len(list(nx.neighbors(nodecouple,loop)))


avgdegree = np.mean(lnk)
mnomega = 1

facval = np.zeros(shape = (Nfac,))
yy = np.zeros(shape=(Nfac,N))
xx = np.zeros(shape=(Nfac,))
for facloop in range(Nfac):
    print(facloop)

    fac = facoef*(16*facloop/(Nfac))*(1/(N-1))*sto/mnomega
    ind = -1
    for nodeloop in nodecouple.node:
        ind = ind + 1
        nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],))
        for linkloop in range (lnk[ind]):
            nodecouple.node[nodeloop]['coupling'][linkloop] = fac

    facval[facloop] = fac*avgdegree
    
    omegout, yout = coupleN(nodecouple)                           # Here is the subfunction call for the flow

    for omloop in range(N):
        yy[facloop,omloop] = omegout[omloop]

    xx[facloop] = facval[facloop]

plt.figure(2)
lines = plt.plot(xx,yy)
plt.setp(lines, linewidth=0.5)
plt.show()

elapsed_time = time.time() - tstart
print('elapsed time = ',format(elapsed_time,'.2f'),'secs')