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.graph.add_nodes_from(self.node_set)
        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)]
            self.graph.add_edges_from(edge_tups)
            
    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
                    newedges.add_edge(a,b)
                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
                        bad_c = {a}
                        c = random.sample(self.node_set - bad_c,1)[0]
                        E = (c,b)
                        while E in newedges.edges():
                            bad_c.add(c)
                            c = random.sample(self.node_set - bad_c,1)[0]
                            E = (c,b)
                        newedges.add_edge(c,b)
                    else:  #changing the SECOND number in the tuple [1]
                            #print('changing the second half')
                        a,b = edge
                        bad_c = {b}
                        c = random.sample(self.node_set - bad_c,1)[0]
                        E = (a,c)
                        while E in newedges.edges():
                            bad_c.add(c)
                            c = random.sample(self.node_set - bad_c,1)[0]
                            E = (a,c)
                        newedges.add_edge(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))
    ax = fig.add_subplot(111,projection='polar')
    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()