#!/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')