# Simulating Watts and Strogatz's Small-World Model¶

In [287]:
import numpy as np
import math
import matplotlib.pyplot as plt
import networkx as nx
import random
from matplotlib import rc
rc('text', usetex=True)
import matplotlib.patches as mpatches

In [288]:
newimage = 3

In [289]:
P = [i/10 for i in range(1,10)]

In [290]:
def circle(N,ms,i,p,graph=None):
plt.figure(figsize=(10,10))
plt.axis('off')
theta = [i*(2*np.pi)/N for i in range(1,N+1)]
X = np.cos(theta)
Y = np.sin(theta)
plt.scatter(X,Y,color='k',s=ms)

if graph:

edges = graph.edges()
for edge in edges:
a,b = edge
thetA = (a*2*np.pi)/N
thetB = (b*2*np.pi)/N
xa,ya = np.cos(thetA),np.sin(thetA)
xb,yb = np.cos(thetB),np.sin(thetB)
plt.plot((xa,xb),(ya,yb))
plt.savefig('Figs/N'+str(N)+'/Graphp'+str(int(p*10))+'/graph'+str(i+1)+'.png',bbox_inches='tight')


In [291]:
class WattsStrogatz:
def __init__(self,N,k,p):
'''
initialize graph with n nodes and average degree k
'''
self.N = N
self.k = k
self.p = p
self.node_set = {n for n in range(N)}
self.graph = nx.Graph()

self.begin_edges()

def begin_edges(self):
'''
build all initial edges on graph
'''
for node in self.graph.nodes():
edge_tups = [(node,(node+i)%len(self.node_set)) for i in range(1,int(self.k/2)+1)]

def rewire(self,change,plist,spot,newedges,n):
for i,edgelist in enumerate(change):
for j,edge in enumerate(edgelist):
if not plist[i][j]:
a,b = edge
else:
#this is where we rewire
#print('rewiring!')
if spot == 1:  #changing the FIRST number in the tuple [0]
#print('changing the first half')
a,b = edge
E = (c,b)
while E in newedges.edges():
E = (c,b)
else:  #changing the SECOND number in the tuple [1]
#print('changing the second half')
a,b = edge
E = (a,c)
while E in newedges.edges():
E = (a,c)
#for GIF-ing
#'''
if self.p in P:
if (n+1)%(self.N*(self.k/2)/20) == 0:
circle(self.N,2,n,self.p,graph=newedges)
#print(newedges.graph)
n+=1
#'''
return n , newedges

def run(self):
'''
Do the shit (run through each node and rewire all links to new parters with probability p)
'''
edges = list(self.graph.edges())
#print(len(edges))
rows = []
for i in range(0,len(edges),int(self.k/2)):
rows.append(edges[i:i+(int(self.k/2))])
#print(len(rows))
changefirst = [rows[i] for i in range(len(rows)) if i%2 == 0]
changesecond= [rows[i] for i in range(len(rows)) if i%2 == 1]
pfirst = [np.random.rand(int(k/2)) < self.p for i in changefirst]
#print(changefirst)
#print(pfirst)
psecond = [np.random.rand(int(k/2)) < self.p for i in changesecond]
newedges = nx.Graph()
n=0
n2=0
#print('I will now change the first set of edges = '+str(len(changefirst)*5))
n2,newedges = self.rewire(changefirst,pfirst,1,newedges,n)
#print('I will now change the second set of edges = '+str(len(changefirst)*5))
nend,newedges.graph = self.rewire(changesecond,psecond,2,newedges,n2)
#print(self.graph)
return newedges.graph



## To Recreate the Watts+Strogatz Graph Goodness¶

In [370]:
N = 1000
k = 10
probs = [0]+list(np.geomspace(0.0001,1,16))
Cs = {}
Ls = {}
for key in probs:
Cs[key]=[]
Ls[key]=[]
for t in range(15):
#C = []
#L = []
#for i in range(0,10000):
#    p = float(i)/10000
for p in probs:
a = WattsStrogatz(N,k,p)
newedges = a.run()
clusters = nx.clustering(newedges.graph)
ckeys = [k for k in clusters]
cvals = [clusters[k] for k in ckeys]
lvals = nx.average_shortest_path_length(newedges.graph)
if p == 0:
C0 = np.mean(cvals)
L0 = lvals
#print(np.mean(cvals)/C0)
#print(lvals/L0)
#C.append(np.mean(cvals)/C0)
#L.append(lvals/L0)
Cs[p].append(np.mean(cvals)/C0)
Ls[p].append(lvals/L0)
#if i%100 == 0:
#print('p = '+str(p)+' is done!')
print('trial '+str(t+1)+' of 15 is done!')
#ps = [i/10000 for i in range(0,10000)]
C=[]
L=[]
for key in Cs.keys():
C.append(np.mean(Cs[key]))
L.append(np.mean(Ls[key]))

plt.close('all')
Cc = mpatches.Patch(color='c', label=r'$\frac{C(p)}{C(0)}$')
Lg = mpatches.Patch(color='g', label=r'$\frac{L(p)}{L(0)}$')
#ps = [float(i)/10000 for i in range(0,10000) ]
plt.figure(figsize=(12,10))
plt.scatter(np.log10(probs),C,color='c')
plt.scatter(np.log10(probs),L,color='g')
#plt.xscale('log')
plt.xlabel(r'$log_{10}(p)$',Fontsize=14)
plt.ylabel(r'$\frac{L(p)}{L(0)}$ or $\frac{C(p)}{C(0)}$',Fontsize=14)
plt.title(r'Characteristic Path Length L(p) and Clustering Coefficient C(p) for the family of graphs rewired with probability p',Fontsize=16)
plt.legend(handles=[Cc,Lg])
plt.grid(linestyle='dotted',linewidth=1)
plt.savefig('Figs/LCvpCurveN-'+str(N),bbox_inches='tight')
plt.show()
#newimage += 1

trial 1 of 15 is done!

/Applications/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:45: RuntimeWarning: divide by zero encountered in log10
/Applications/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:46: RuntimeWarning: divide by zero encountered in log10


## For the Purposes of GIF'ing¶

In [286]:
N = 1000
k = 10
probs = [i/10 for i in range(1,10)]

for p in probs:
a = WattsStrogatz(N,k,p)
newedges = a.run()
print('p = '+str(p)+' is done!')

/Applications/anaconda/lib/python3.6/site-packages/matplotlib/pyplot.py:524: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (matplotlib.pyplot.figure) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam figure.max_open_warning).
max_open_warning, RuntimeWarning)

p = 0.1 is done!
p = 0.2 is done!
p = 0.3 is done!
p = 0.4 is done!
p = 0.5 is done!
p = 0.6 is done!
p = 0.7 is done!
p = 0.8 is done!
p = 0.9 is done!


### test for 1 value of p¶

In [209]:
N = 1000
k = 10
p = 0.01
a = WattsStrogatz(N,k,p)
newedges = a.run()
clusters = nx.clustering(newedges.graph)
ckeys = [k for k in clusters]
cvals = [clusters[k] for k in ckeys]
lvals = nx.average_shortest_path_length(newedges.graph)
print(np.mean(cvals)/C0)
print(lvals/L0)

0.974863551167
0.2003500196463654


# More Peculiar Power Law Tails¶

In [402]:
import numpy as np
import math
import matplotlib.pyplot as plt
import networkx as nx
import random
from matplotlib import rc
rc('text', usetex=True)
import matplotlib.patches as mpatches
from scipy import stats

In [337]:
n = 1000
Ns = [10**(i) for i in range(1,7)]
Max = {}
for key in Ns:
Max[key] = []

for N in Ns:
for i in range(n):
Z = np.random.rand(N)
z = max(Z)
x = (1-z)**(-2/3)
Max[N].append(x)

In [409]:
ns = [i for i in range(n)]
Maxs = [np.mean(Max[N]) for N in Ns]
for N in Ns:
plt.figure(figsize = (13,5))
plt.plot(ns,Max[N],color='k',linewidth=0.5)
plt.xlabel('n = trial number out of '+str(n),Fontsize=16)
plt.ylabel(r'k_{max}',Fontsize=16)
plt.title(r'Max Values $(k_{max})$ Taken from Distribution $P_k = \frac{3}{2}k^{-5/2}$ for '+str(n)+' Iterations of '+str(N)+' Samples',Fontsize=20)
plt.savefig('Figs/kMaxPlotN'+str(N),bbox_inches='tight')
plt.show()
i = [i for i in range(1,7)]
plt.figure(figsize = (13,10))
plt.scatter(np.log10(Ns),np.log10(Maxs),s=100,marker='*',color='k')
#plt.plot(i,[j*(2/3)+0.4 for j in i],linestyle='dashed',linewidth=0.5,color='olive')
plt.xlabel(r'$log_{10}(N)$',Fontsize=16)
plt.ylabel(r'$log_{10}(k_{max})$',Fontsize=16)
plt.grid(linestyle='dotted',linewidth=0.5)
plt.title(r'Average $k_{max}$ for Sample Size, N out of '+str(n)+' Trials',Fontsize=20)
slope, intercept, r_value, p_value, std_err = stats.linregress(np.log10(Ns),np.log10(Maxs))
f = lambda x: slope*x + intercept
x = np.array([0.5,6])
plt.plot(x,f(x),lw = 0.5, color = 'olive',label = r'Regression, $m=$'+str(slope))
plt.legend(loc=2,fontsize=14)
#plt.text('loc=2,fontsize=14)
plt.savefig('Figs/AveragekMaxvsN',bbox_inches='tight')
plt.show()


### Polar coordinate sillyness¶

In [64]:
def toPolar(x,y):
r = np.sqrt(x**2+y**2)
theta = np.arctan(y/x)
return r,theta

def polarcircle(N, WS = None):
fig = plt.figure(figsize=(12,12))
total = 0
for i in range(N):
ax.scatter(total,2,color='k',s=2)
total += 2*np.pi/N
plt.xticks([])
plt.yticks([])
plt.axis('off')
'''
if WS:
edges = WS.graph.edges()
for edge in edges:
x,y = edge
if x == 0:
x = 10**(-10)
r, theta = toPolar(x,y)
ax.plot(theta,r)
'''
plt.show()