Directed modularity and tests

1 files changed, 97 insertions, 4 deletions

diff --git a/networkx/linalg/modularitymatrix.py b/networkx/linalg/modularitymatrix.py
index b0b7faac..11a26d76 100644
--- a/networkx/linalg/modularitymatrix.py
+++ b/networkx/linalg/modularitymatrix.py
@@ -13,8 +13,9 @@ from networkx.utils import not_implemented_for
 __author__ = "\n".join(['Aric Hagberg <aric.hagberg@gmail.com>',
                         'Pieter Swart (swart@lanl.gov)',
                         'Dan Schult (dschult@colgate.edu)',
-                        'Jean-Gabriel Young <jean.gabriel.young@gmail.com>'])
-__all__ = ['modularity_matrix']
+                        'Jean-Gabriel Young (Jean.gabriel.young@gmail.com)'])
+
+__all__ = ['modularity_matrix', 'directed_modularity_matrix']
 
 
 @not_implemented_for('directed')
@@ -24,11 +25,16 @@ def modularity_matrix(G, nodelist=None):
 
     The modularity matrix is the matrix B = A - <A>, where A is the adjacency
     matrix and <A> is the average adjacency matrix, assuming that the graph
-    is descibred by the configuration model.
+    is described by the configuration model.
+
+    More specifically, the element B_ij of B is defined as
+        A_ij - k_i k_j/m
+    where k_i(in) is the degree of node i, and were m is the number of edges
+    in the graph.
 
     Parameters
     ----------
-    G : graph
+    G : Graph
        A NetworkX graph
 
     nodelist : list, optional
@@ -40,6 +46,20 @@ def modularity_matrix(G, nodelist=None):
     B : Numpy matrix
       The modularity matrix of G.
 
+    Examples
+    --------
+    >>> import networkx as nx
+    >>> k =[3, 2, 2, 1, 0]
+    >>> G = nx.havel_hakimi_graph(k)
+    >>> B = nx.modularity_matrix(G)
+    >>> print B
+    [[-1.125  0.25   0.25   0.625  0.   ]
+     [ 0.25  -0.5    0.5   -0.25   0.   ]
+     [ 0.25   0.5   -0.5   -0.25   0.   ]
+     [ 0.625 -0.25  -0.25  -0.125  0.   ]
+     [ 0.     0.     0.     0.     0.   ]]
+
+
     See Also
     --------
     to_numpy_matrix
@@ -61,6 +81,79 @@ def modularity_matrix(G, nodelist=None):
     return A - X
 
 
+@not_implemented_for('undirected')
+@not_implemented_for('multigraph')
+def directed_modularity_matrix(G, nodelist=None):
+    """Return the directed modularity matrix of G.
+
+    The modularity matrix is the matrix B = A - <A>, where A is the adjacency
+    matrix and <A> is the expected adjacency matrix, assuming that the graph
+    is described by the configuration model.
+
+    More specifically, the element B_ij of B is defined as
+        B_ij = A_ij - k_i(out) k_j(in)/m
+    where k_i(in) is the in degree of node i, and k_j(out) is the out degree
+    of node j, with m the number of edges in the graph.
+
+    Parameters
+    ----------
+    G : DiGraph
+       A NetworkX DiGraph
+
+    nodelist : list, optional
+       The rows and columns are ordered according to the nodes in nodelist.
+       If nodelist is None, then the ordering is produced by G.nodes().
+
+    Returns
+    -------
+    B : Numpy matrix
+      The modularity matrix of G.
+
+    Examples
+    --------
+    >>> import networkx as nx
+    >>> G = nx.DiGraph()
+    >>> G.add_edges_from(((1,2), (1,3), (3,1), (3,2), (3,5), (4,5), (4,6),
+    ...                   (5,4), (5,6), (6,4)))
+    >>> B = nx.directed_modularity_matrix(G)
+    >>> print B
+    [[-0.2  0.6  0.8 -0.4 -0.4 -0.4]
+     [ 0.   0.   0.   0.   0.   0. ]
+     [ 0.7  0.4 -0.3 -0.6  0.4 -0.6]
+     [-0.2 -0.4 -0.2 -0.4  0.6  0.6]
+     [-0.2 -0.4 -0.2  0.6 -0.4  0.6]
+     [-0.1 -0.2 -0.1  0.8 -0.2 -0.2]]
+
+
+    Notes
+    _____
+    NetworkX defines the element A_ij of the adjacency matrix as 1 if there
+    is a link going from node i to node j. Leicht and Newman use the opposite
+    definition. This explains the different expression for B_ij.
+
+    See Also
+    --------
+    to_numpy_matrix
+    adjacency_matrix
+    laplacian_matrix
+
+    References
+    ----------
+    .. [1] E. A. Leicht, M. E. J. Newman, 
+       "Community structure in directed networks",
+        Phys. Rev Lett., vol. 100, no. 11, p. 118703, 2008.
+    """
+    if nodelist is None:
+        nodelist = G.nodes()
+    A = nx.to_scipy_sparse_matrix(G, nodelist=nodelist, format='csr')
+    k_in = A.sum(axis=0)
+    k_out = A.sum(axis=1)
+    m = G.number_of_edges()
+    # Expected adjacency matrix
+    X = k_out * k_in / m
+    return A - X
+
+
 # fixture for nose tests
 def setup_module(module):
     from nose import SkipTest