Add support for directed graphs

README.md

@@ -7,6 +7,7 @@ partitions that are significant from the perspective of stochastic block model
 inference. It:
 
 - Works for single-layer and multi-layer networks
+- Supports both directed and undirected networks
 - Supports weighted networks via the `weight` edge attribute
 - Restricts focus to a fixed number of communities when desired
 - Contains a faster and more robust implementation of the CHAMP (Convex Hull of
@@ -19,6 +20,8 @@ information.
 Significantly more details can be found in the article's
 [Supplementary Information](https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-022-20142-6/MediaObjects/41598_2022_20142_MOESM1_ESM.pdf).
 
+Note: support for weighted and directed networks is an extension beyond the original publication.
+
 ## Installation
 
 This project is on [PyPI](https://pypi.org/project/modularitypruning/) and can

docs/index.rst

@@ -13,6 +13,7 @@ partitions that are significant from the perspective of stochastic block model
 inference. It:
 
 - Works for single-layer and multi-layer networks
+- Supports both directed and undirected networks
 - Supports weighted networks via the `weight` edge attribute
 - Restricts focus to a fixed number of communities when desired
 - Contains a faster and more robust implementation of the CHAMP (Convex Hull of
@@ -23,6 +24,8 @@ article's significantly more detailed `Supplementary Information`_.
 
 .. _Supplementary Information: https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-022-20142-6/MediaObjects/41598_2022_20142_MOESM1_ESM.pdf
 
+Note: support for weighted and directed networks is an extension beyond the original publication.
+
 Prior to version v1.3.0, we used the Louvain algorithm for modularity maximization instead of Leiden. The deprecated
 module ``modularitypruning.louvain_utilities`` now shims single-layer functions to their corresponding Leiden versions
 in ``modularitypruning.leiden_utilities`` (though it still contains the legacy multi-layer functions since they can be

experiments/lazega_law_firm/lazega_figures.py

@@ -1,4 +1,7 @@
 # Generates figures 5.7 through 5.10
+# TODO: We support directed parameter estimation now, but Roxana et al. only used the undirected theory.
+#   Small tweaks would be needed to this script to avoid changing the gamma estimates, probably by using a
+#   directed=False copy of the intralayer graph for the parameter estimation step.
 
 import os
 from modularitypruning.champ_utilities import CHAMP_3D

tests/shared_testing_functions.py

@@ -63,7 +63,7 @@ def generate_igraph_famous():
 
 
 def generate_connected_multilayer_intralayer_SBM(copying_probability, p_in, p_out, first_layer_membership, num_layers,
-                                                 weighted=False, low=0.5, high=5.0):
+                                                 weighted=False, low=0.5, high=5.0, directed=False):
     num_nodes_per_layer = len(first_layer_membership)
     community_labels_per_layer = [[0] * num_nodes_per_layer for _ in range(num_layers)]
     community_labels_per_layer[0] = list(first_layer_membership)
@@ -82,15 +82,22 @@ def generate_connected_multilayer_intralayer_SBM(copying_probability, p_in, p_ou
     combined_community_labels = sum(community_labels_per_layer, [])
     layer_membership = [i for i in range(num_layers) for _ in range(num_nodes_per_layer)]
 
-    for v in range(len(combined_community_labels)):
-        for u in range(v + 1, len(combined_community_labels)):
-            if layer_membership[v] == layer_membership[u]:
-                if combined_community_labels[v] == combined_community_labels[u]:
-                    if random() < p_in:
-                        intralayer_edges.append((u, v))
-                else:
-                    if random() < p_out:
-                        intralayer_edges.append((u, v))
+    if directed:
+        edges_per_layer = [(v, u) for v in range(len(combined_community_labels))
+                           for u in range(len(combined_community_labels))
+                           if layer_membership[v] == layer_membership[u] and v != u]
+    else:
+        edges_per_layer = [(v, u) for v in range(len(combined_community_labels))
+                           for u in range(v + 1, len(combined_community_labels))
+                           if layer_membership[v] == layer_membership[u]]
+
+    for v, u in edges_per_layer:
+        if combined_community_labels[v] == combined_community_labels[u]:
+            if random() < p_in:
+                intralayer_edges.append((v, u))
+        else:
+            if random() < p_out:
+                intralayer_edges.append((v, u))
 
     # force disconnected intralayer nodes to have at least one edge
     intralayer_degrees = [0] * len(combined_community_labels)
@@ -103,14 +110,72 @@ def generate_connected_multilayer_intralayer_SBM(copying_probability, p_in, p_ou
                           if layer_membership[u] == layer_membership[v] and u != v]
             intralayer_edges.append((v, same_layer[randint(0, len(same_layer) - 1)]))
 
-    G_intralayer = ig.Graph(n=len(combined_community_labels), edges=intralayer_edges, directed=False)
+    G_intralayer = ig.Graph(n=len(combined_community_labels), edges=intralayer_edges, directed=directed)
 
     if weighted:
         G_intralayer.es['weight'] = generate_random_values(G_intralayer.ecount(), low, high)
 
     return G_intralayer, layer_membership
 
 
+def generate_temporal_interlayer_edges(num_nodes_per_layer, num_layers):
+    """Directed interlayer temporal edges (each node connects to itself in the next layer)."""
+    return [(num_nodes_per_layer * layer + v, num_nodes_per_layer * (layer + 1) + v)
+            for layer in range(num_layers - 1) for v in range(num_nodes_per_layer)]
+
+
+def generate_multiplex_interlayer_edges(num_nodes_per_layer, num_layers):
+    """Directed interlayer multiplex edges (each node connects to itself in all other layers)."""
+    return [(num_nodes_per_layer * layer_start + v, num_nodes_per_layer * layer_end + v)
+            for layer_start in range(num_layers) for layer_end in range(num_layers) if layer_start != layer_end
+            for v in range(num_nodes_per_layer)]
+
+
+def generate_temporal_SBM(copying_probability, p_in, p_out, first_layer_membership, num_layers,
+                          weighted=False, directed=False):
+    """Generate an intralayer SBM with temporal interlayer edges.
+
+    As in Pamfil et al., the interlayer edges are directed in one direction."""
+    G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM(
+        copying_probability, p_in, p_out, first_layer_membership, num_layers, weighted=weighted,
+        directed=directed)
+    G_interlayer = ig.Graph(generate_temporal_interlayer_edges(len(first_layer_membership), num_layers),
+                            directed=True)
+    return G_intralayer, G_interlayer, layer_membership
+
+
+def generate_multiplex_SBM(copying_probability, p_in, p_out, first_layer_membership, num_layers,
+                           weighted=False, directed=False):
+    """Generate an intralayer SBM with multiplex interlayer edges.
+
+    As in Pamfil et al., the interlayer edges are directed in both directions."""
+    G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM(
+        copying_probability, p_in, p_out, first_layer_membership, num_layers, weighted=weighted,
+        directed=directed)
+    G_interlayer = ig.Graph(generate_multiplex_interlayer_edges(len(first_layer_membership), num_layers),
+                            directed=True)
+    return G_intralayer, G_interlayer, layer_membership
+
+
+def generate_SBM(q, directed=False, weighted=False, low=0.5, high=5.0):
+    """Generate a planted-partition SBM, optionally directed and/or weighted.
+
+    This is a setup analogous to FIG 1 of Newman's 'Equivalence between modularity optimization and maximum likelihood
+    methods for community detection'.
+    """
+    community_sizes = [250] * q
+    n = sum(community_sizes)
+    p_in = 16 * n / (q * 250 * 249)  # ~16 in-edges per node
+    p_out = 8 * n / (q * (q - 1) * 250 * 250)  # ~8 out-edges per node, spread over the other communities
+    pref_matrix = [[p_in if i == j else p_out for j in range(q)] for i in range(q)]
+    G = ig.Graph.SBM(pref_matrix, community_sizes, directed=directed)
+
+    if weighted:
+        G.es['weight'] = generate_random_values(G.ecount(), low, high)
+
+    return G, p_in, p_out
+
+
 def assert_almost_equal_or_both_none_or_nan(self, first, second, places=10):
     if first is None:
         self.assertIsNone(second)  # isnan only accepts floats

tests/test_champ_coefficients_2D.py

@@ -4,7 +4,7 @@
 from modularitypruning.champ_utilities import partition_coefficients_2D
 from modularitypruning.leiden_utilities import leiden_part_with_membership, repeated_leiden_from_gammas
 
-from .shared_testing_functions import generate_connected_ER, \
+from .shared_testing_functions import generate_connected_ER, generate_SBM, \
     generate_random_values, generate_random_partitions, generate_igraph_famous
 
 
@@ -87,6 +87,16 @@ def test_partition_coefficient_correctness_directed_varying_K_max(self):
                     self.assert_partition_coefficient_correctness_ER(directed=True, num_partitions=100, K_max=K_max,
                                                                      weighted=weighted)
 
+    def test_partition_coefficient_correctness_directed_SBM(self):
+        """Test on asymmetric directed SBMs, unweighted and weighted."""
+        for q in [2, 3, 5]:
+            for weighted in [False, True]:
+                with self.subTest(q=q, weighted=weighted):
+                    G, _, _ = generate_SBM(q, directed=True, weighted=weighted)
+                    partitions = generate_random_partitions(num_nodes=G.vcount(), num_partitions=100, K_max=q)
+                    coefficients = partition_coefficients_2D(G, partitions)
+                    self.assert_partition_coefficient_correctness(G, partitions, coefficients)
+
     def test_partition_coefficient_correctness_igraph_famous_leiden(self):
         """Test partition coefficient correctness on various famous graphs while obtaining partitions via Leiden.
 

tests/test_champ_coefficients_3D.py

@@ -4,7 +4,8 @@
 from modularitypruning.champ_utilities import partition_coefficients_3D
 from modularitypruning.leiden_utilities import multilayer_leiden_part_with_membership, leiden_part_with_membership
 
-from .shared_testing_functions import generate_connected_multilayer_ER, generate_random_partitions
+from .shared_testing_functions import generate_connected_multilayer_ER, generate_random_partitions, \
+    generate_temporal_SBM
 
 
 class TestCHAMPCoefficients3D(unittest.TestCase):
@@ -126,6 +127,27 @@ def test_partition_coefficient_correctness_directed_varying_K_max(self):
                 with self.subTest(K_max=K_max, weighted=weighted):
                     self.assert_partition_coefficient_correctness_ER(directed=True, K_max=K_max, weighted=weighted)
 
+    def assert_partition_coefficient_correctness_SBM(self, num_layers, weighted):
+        num_nodes_per_layer = 100
+        first_layer_membership = [i // (num_nodes_per_layer // 3) for i in range(num_nodes_per_layer)]
+
+        G_intralayer, G_interlayer, layer_membership = generate_temporal_SBM(
+            copying_probability=0.9, p_in=0.25, p_out=0.05, first_layer_membership=first_layer_membership,
+            num_layers=num_layers, weighted=weighted, directed=True)
+
+        partitions = generate_random_partitions(num_nodes=G_intralayer.vcount(), num_partitions=100, K_max=3)
+        coefficients = partition_coefficients_3D(G_intralayer, G_interlayer, layer_membership, partitions)
+
+        self.assert_partition_coefficient_correctness(G_intralayer, G_interlayer, layer_membership, partitions,
+                                                      coefficients)
+
+    def test_partition_coefficient_correctness_directed_SBM(self):
+        """Test on asymmetric directed SBMs, unweighted and weighted."""
+        for num_layers in [4, 5]:
+            for weighted in [False, True]:
+                with self.subTest(num_layers=num_layers, weighted=weighted):
+                    self.assert_partition_coefficient_correctness_SBM(num_layers=num_layers, weighted=weighted)
+
     def test_partition_coefficient_correctness_interleaved_directedness(self):
         """Test partition coefficient correctness when directedness of interlayer and intralayer edges do not match."""
         for weighted in [False, True]:

tests/test_monolayer_parameter_estimation.py

@@ -5,11 +5,11 @@
 import igraph as ig
 from modularitypruning.parameter_estimation import iterative_monolayer_resolution_parameter_estimation
 from modularitypruning.parameter_estimation_utilities import gamma_estimate
-from modularitypruning.partition_utilities import all_degrees, all_strengths
+from modularitypruning.partition_utilities import (all_strengths, in_strengths, out_strengths)
 from numpy import mean
 
-from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_igraph_famous, \
-    generate_random_partition, generate_random_values
+from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_SBM, \
+    generate_igraph_famous, generate_random_partition
 
 
 class TestMonolayerParameterEstimation(unittest.TestCase):
@@ -24,34 +24,16 @@ def test_newman_synthetic_networks(self):
                 with self.subTest(q=q, weighted=weighted):
                     self.assert_newman_synthetic_convergence(q, weighted=weighted)
 
-    @staticmethod
-    def generate_monolayer_SBM(q, weighted=False):
-        community_sizes = [250] * q  # q equally sized groups of 250 nodes
-        n = sum(community_sizes)
-        p_in = 16 * n / (q * 250 * 249)  # ~16 in-edges per node
-        p_out = 8 * n / (q * (q - 1) * 250 * 250)  # ~8 out-edges per node to each community
-        pref_matrix = [[p_in if i == j else p_out for j in range(q)] for i in range(q)]
-        G = ig.Graph.SBM(pref_matrix, community_sizes)
-
-        if weighted:
-            G.es['weight'] = generate_random_values(G.ecount(), 0.5, 5.0)
-
-        return G, p_in, p_out
-
     def assert_newman_synthetic_convergence(self, q, weighted=False):
-        G, p_in, p_out = self.generate_monolayer_SBM(q, weighted=weighted)
+        G, p_in, p_out = generate_SBM(q, weighted=weighted)
 
         # compute "ground truth" gamma estimate
-        if weighted:
-            k = mean(all_strengths(G))
-            m = sum(G.es['weight'])
-            mean_weight = m / G.ecount()
-            true_omega_in = mean_weight * p_in * (2 * m) / (k * k)
-            true_omega_out = mean_weight * p_out * (2 * m) / (k * k)
-        else:
-            k = mean(all_degrees(G))
-            true_omega_in = p_in * (2 * G.ecount()) / (k * k)
-            true_omega_out = p_out * (2 * G.ecount()) / (k * k)
+        # note that strengths are just degrees in unweighted graphs
+        m = sum(G.es['weight']) if weighted else G.ecount()
+        mean_weight = m / G.ecount() if weighted else 1.0
+        k = mean(all_strengths(G))
+        true_omega_in = mean_weight * p_in * (2 * m) / (k * k)
+        true_omega_out = mean_weight * p_out * (2 * m) / (k * k)
         true_gamma = (true_omega_in - true_omega_out) / (log(true_omega_in) - log(true_omega_out))
 
         gamma, _ = iterative_monolayer_resolution_parameter_estimation(G, gamma=1.0)
@@ -72,6 +54,52 @@ def test_directed_consistency_igraph_famous(self):
             with self.subTest(graph=name):
                 assert_almost_equal_or_both_none_or_nan(self, gamma_undirected, gamma_directed, places=10)
 
+    def assert_newman_directed_synthetic_convergence(self, q, weighted=False):
+        """Test directed analogue of Newman's synthetic test."""
+        G, p_in, p_out = generate_SBM(q, directed=True, weighted=weighted)
+
+        # compute "ground truth" gamma estimate for a directed graph
+        M = sum(G.es['weight']) if weighted else G.ecount()
+        mean_weight = M / G.ecount() if weighted else 1.0
+        k_in = mean(in_strengths(G))
+        k_out = mean(out_strengths(G))
+        true_omega_in = mean_weight * p_in * M / (k_in * k_out)
+        true_omega_out = mean_weight * p_out * M / (k_in * k_out)
+        true_gamma = (true_omega_in - true_omega_out) / (log(true_omega_in) - log(true_omega_out))
+
+        gamma, _ = iterative_monolayer_resolution_parameter_estimation(G, gamma=1.0)
+
+        # check we converged close to the ground truth "correct" value
+        self.assertLess(abs(true_gamma - gamma), 0.05)
+
+    def test_directed_newman_synthetic_networks(self):
+        """Directed analogue of test_newman_synthetic_networks."""
+
+        for weighted in [False, True]:
+            for q in range(3, 15):
+                with self.subTest(q=q, weighted=weighted):
+                    self.assert_newman_directed_synthetic_convergence(q, weighted=weighted)
+
+    def test_directed_difference(self):
+        """Simple test to confirm parameter estimate depends on edge direction.
+
+        The graph here has all edges pointing from the first community to the second."""
+        edges = [(0, 1), (2, 3), (0, 2), (0, 3), (1, 2), (1, 3)]
+        membership = (0, 0, 1, 1)
+
+        for i, weights in enumerate([[1.0] * 6, [1.0, 1.0, 3.0, 3.0, 3.0, 3.0]]):
+            with self.subTest(weight_set=i):
+                G = ig.Graph(4, edges, directed=True)
+                G.es['weight'] = weights
+                gamma_directed = gamma_estimate(G, membership)
+
+                # undirected version for comparison
+                G_undirected = ig.Graph(4, edges, directed=False)
+                G_undirected.es['weight'] = weights
+                gamma_undirected_treatment = gamma_estimate(G_undirected, membership)
+
+                self.assertGreater(abs(gamma_directed - gamma_undirected_treatment), 1e-2)
+
 
 if __name__ == "__main__":
     seed(0)