An a posteriori measure of network modularity -- code associated to the F1000R paper

qr.py

import networkx as nx
import numpy as np
import scipy as sp
from scipy import stats
import matplotlib.pyplot as plt
import progressbar as pb # Get it from http://code.google.com/p/python-progressbar/
import community # Get it from http://perso.crans.org/aynaud/communities/

def nullModel(graph):
	"""
	This function returns as randomized graph.
	The distribution of in and out degrees is the same as in the original graph.
	"""
	N = nx.DiGraph()
	N.add_nodes_from(graph.node)
	S = N.number_of_nodes()
	for n1 in N:
		for n2 in N:
			Pout = graph.out_degree()[n1]/float(S)
			Pin  =  graph.in_degree()[n2]/float(S)
			Pint = (Pout+Pin)/float(2)
			if np.random.binomial(1,Pint) == 1:
				N.add_edge(n1,n2)
	return N

def Qr(G):
	"""
	Measures the realized modularity, expressed in the 0-1 range.
	"""
	E = G.number_of_edges()
	W = 0
	for e in G.edges():
		if G.node[e[0]]['m'] == G.node[e[1]]['m']:
			W += 1
	return 2*(W/float(E))-1

def Mod(G,usebest=True,l=1):
	"""
	A wrapper arround the functions in community.
	Returns
		the modularity
		the number of modules
		the realized modularity
		the number of edges
		the number of nodes
	"""
	D = G.to_undirected()
	dendo = community.generate_dendogram(D, None)
	if usebest:
		level = len(dendo)-1
	else:
		level = l
	partition = community.partition_at_level(dendo,level)
	mod = community.modularity(partition, D)
	for n in G:
		G.node[n]['m'] = partition[n]
	nmod = len(list(set(partition.values())))
	return [mod, nmod, Qr(G), G.number_of_edges(), G.number_of_nodes()]

####### EXAMPLE APPLICATION

filename = 'network.txt'

# Reading the network
G = nx.read_edgelist(filename, create_using=nx.DiGraph())

# Modularity of the empirical network
m = Mod(G,True)

# Null model
mQ = []
mQr = []

widgets = ['Randomization ', pb.Percentage(), ' ', pb.Bar(),' ', pb.ETA()]
progress = pb.ProgressBar(widgets=widgets)
for re in progress(xrange(1000)):
	nG = nullModel(G)
	try:
		nullm = Mod(nG,True)
		mQr.append(nullm[2])
		mQ.append(nullm[0])
	except:
		pass

# Mean and confidence interval
nQ = stats.bayes_mvs(mQ,0.95)[0]
nQr = stats.bayes_mvs(mQr,0.95)[0]

print nQ
print nQr

# Quick plot of modularity
plt.hist(mQ,bins=20,histtype='step')
plt.title("Modularity")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.axvline(x=m[0])
plt.show()

# Quick plot of realized modularity
plt.hist(mQr,bins=20,histtype='step')
plt.title("Modularity")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.axvline(x=m[2])
plt.show()