From 94a0f84eacd8ce3f145eedda7067989d7207915e Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 19:34:32 -0400 Subject: [PATCH 01/13] Add directed SBM parameter estimation --- utilities/parameter_estimation_utilities.py | 75 ++++++++++++++++----- 1 file changed, 58 insertions(+), 17 deletions(-) diff --git a/utilities/parameter_estimation_utilities.py b/utilities/parameter_estimation_utilities.py index 1814cd5..7ac87c6 100644 --- a/utilities/parameter_estimation_utilities.py +++ b/utilities/parameter_estimation_utilities.py @@ -24,6 +24,9 @@ def estimate_singlelayer_SBM_parameters( See https://doi.org/10.1103/PhysRevE.94.052315 for more details. + For directed graphs, the convention of Leicht and Newman is used: each directed edge is counted once, and the total + edge weight is :math:`M = \sum_i k_i^{in}`. Self-loops count as internal edges. + :param G: graph of interest :type G: igraph.Graph :param partition: partition of interest @@ -44,15 +47,31 @@ def estimate_singlelayer_SBM_parameters( community: list[int] = partition.membership m_in: float = sum(e['weight'] * (community[e.source] == community[e.target]) for e in G.es) - kappa_r_list: list[float] = [0] * len(partition) - for e in G.es: - kappa_r_list[community[e.source]] += e['weight'] - kappa_r_list[community[e.target]] += e['weight'] - sum_kappa_sqr: float = sum(x ** 2 for x in kappa_r_list) - omega_in: float = (2 * m_in) / (sum_kappa_sqr / (2 * m)) - # guard for div by zero with single community partition - omega_out: float = (2 * m - 2 * m_in) / (2 * m - sum_kappa_sqr / (2 * m)) if len(partition) > 1 else 0 + if G.is_directed(): + # directed modularity of Leicht and Newman + # null model term sum_c k_c^{in} * k_c^{out} / M + # internal edge weights are counted once per directed edge + kappa_in_list: list[float] = [0] * len(partition) + kappa_out_list: list[float] = [0] * len(partition) + for e in G.es: + kappa_out_list[community[e.source]] += e['weight'] + kappa_in_list[community[e.target]] += e['weight'] + sum_kappa_in_kappa_out: float = sum(k_in * k_out for k_in, k_out in zip(kappa_in_list, kappa_out_list)) + + omega_in: float = m_in / (sum_kappa_in_kappa_out / m) + # guard for div by zero with single community partition + omega_out: float = (m - m_in) / (m - sum_kappa_in_kappa_out / m) if len(partition) > 1 else 0 + else: + kappa_r_list: list[float] = [0] * len(partition) + for e in G.es: + kappa_r_list[community[e.source]] += e['weight'] + kappa_r_list[community[e.target]] += e['weight'] + sum_kappa_sqr: float = sum(x ** 2 for x in kappa_r_list) + + omega_in = (2 * m_in) / (sum_kappa_sqr / (2 * m)) + # guard for div by zero with single community partition + omega_out = (2 * m - 2 * m_in) / (2 * m - sum_kappa_sqr / (2 * m)) if len(partition) > 1 else 0 # return estimates for omega_in, omega_out return omega_in, omega_out @@ -117,17 +136,39 @@ def estimate_multilayer_SBM_parameters( m_t_in[layer_vec[e.source]] += e['weight'] kappa_t_r_list: list[list[int]] = [[0] * K for _ in range(T)] - for e in G_intralayer.es: - layer = layer_vec[e.source] - kappa_t_r_list[layer][community[e.source]] += e['weight'] - kappa_t_r_list[layer][community[e.target]] += e['weight'] - sum_kappa_t_sqr: list[float] = [sum(x ** 2 for x in kappa_t_r_list[t]) for t in range(T)] - theta_in: float = sum(2 * m_t_in[t] for t in range(T)) / sum(sum_kappa_t_sqr[t] / (2 * m_t[t]) for t in range(T)) + if G_intralayer.is_directed(): + # directed modularity per layer as in estimate_singlelayer_SBM_parameters + kappa_t_in_list: list[list[float]] = [[0] * K for _ in range(T)] + kappa_t_out_list: list[list[float]] = [[0] * K for _ in range(T)] + for e in G_intralayer.es: + layer = layer_vec[e.source] + kappa_t_out_list[layer][community[e.source]] += e['weight'] + kappa_t_in_list[layer][community[e.target]] += e['weight'] + sum_kappa_in_kappa_out_t: list[float] = [ + sum(k_in * k_out for k_in, k_out in zip(kappa_t_in_list[t], kappa_t_out_list[t])) + for t in range(T)] + + theta_in: float = sum(m_t_in[t] for t in range(T)) / sum(sum_kappa_in_kappa_out_t[t] / m_t[t] + for t in range(T)) + + # guard for div by zero with e.g. a single community partition + theta_out_numerator: float = sum(m_t[t] - m_t_in[t] for t in range(T)) + theta_out_denominator: float = sum(m_t[t] - sum_kappa_in_kappa_out_t[t] / m_t[t] for t in range(T)) + else: + for e in G_intralayer.es: + layer = layer_vec[e.source] + kappa_t_r_list[layer][community[e.source]] += e['weight'] + kappa_t_r_list[layer][community[e.target]] += e['weight'] + sum_kappa_t_sqr: list[float] = [sum(x ** 2 for x in kappa_t_r_list[t]) for t in range(T)] + + theta_in = sum(2 * m_t_in[t] for t in range(T)) / sum(sum_kappa_t_sqr[t] / (2 * m_t[t]) + for t in range(T)) + + # guard for div by zero with e.g. a single community partition + theta_out_numerator = sum(2 * m_t[t] - 2 * m_t_in[t] for t in range(T)) + theta_out_denominator = sum(2 * m_t[t] - sum_kappa_t_sqr[t] / (2 * m_t[t]) for t in range(T)) - # guard for div by zero with e.g. a single community partition - theta_out_numerator: float = sum(2 * m_t[t] - 2 * m_t_in[t] for t in range(T)) - theta_out_denominator: float = sum(2 * m_t[t] - sum_kappa_t_sqr[t] / (2 * m_t[t]) for t in range(T)) if theta_out_denominator == 0: theta_out: float = 0 else: From 1adc9a6055900d7e249c72ff46b610c7eb0b049d Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 19:58:09 -0400 Subject: [PATCH 02/13] Add directed SBM generation for tests --- tests/shared_testing_functions.py | 48 ++++++++++++++++++++++++------- 1 file changed, 37 insertions(+), 11 deletions(-) diff --git a/tests/shared_testing_functions.py b/tests/shared_testing_functions.py index 766a62b..ebc1da4 100644 --- a/tests/shared_testing_functions.py +++ b/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,7 +110,7 @@ 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) @@ -111,6 +118,25 @@ def generate_connected_multilayer_intralayer_SBM(copying_probability, p_in, p_ou return G_intralayer, 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 From e91bd1e7ecf9a23d7072e79754543409ac48fcb6 Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 20:11:07 -0400 Subject: [PATCH 03/13] Add more temporal and multiplex helpers for tests --- tests/shared_testing_functions.py | 39 +++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/tests/shared_testing_functions.py b/tests/shared_testing_functions.py index ebc1da4..9c74993 100644 --- a/tests/shared_testing_functions.py +++ b/tests/shared_testing_functions.py @@ -118,6 +118,45 @@ def generate_connected_multilayer_intralayer_SBM(copying_probability, p_in, p_ou 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. From 6e490cad74eb3f90515d8c47c644e9e730d4003a Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 20:15:23 -0400 Subject: [PATCH 04/13] Add directed monolayer estimation tests --- tests/test_monolayer_parameter_estimation.py | 68 ++++++++++++++------ 1 file changed, 50 insertions(+), 18 deletions(-) diff --git a/tests/test_monolayer_parameter_estimation.py b/tests/test_monolayer_parameter_estimation.py index ded178e..c00bc67 100644 --- a/tests/test_monolayer_parameter_estimation.py +++ b/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_degrees, 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,22 +24,8 @@ 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: @@ -72,6 +58,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) From 8d0969fe602dc28b71f81f9687e04be8312e8e8c Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 20:47:55 -0400 Subject: [PATCH 05/13] Add directed temporal/multilevel parameter estimation tests --- tests/test_monolayer_parameter_estimation.py | 1 + ...emporal_multilevel_parameter_estimation.py | 84 +++++++++++++------ 2 files changed, 58 insertions(+), 27 deletions(-) diff --git a/tests/test_monolayer_parameter_estimation.py b/tests/test_monolayer_parameter_estimation.py index c00bc67..25d3713 100644 --- a/tests/test_monolayer_parameter_estimation.py +++ b/tests/test_monolayer_parameter_estimation.py @@ -28,6 +28,7 @@ def assert_newman_synthetic_convergence(self, q, weighted=False): G, p_in, p_out = generate_SBM(q, weighted=weighted) # compute "ground truth" gamma estimate + # TODO: Why mean weight scaling for weighted graphs? Confirm this. if weighted: k = mean(all_strengths(G)) m = sum(G.es['weight']) diff --git a/tests/test_temporal_multilevel_parameter_estimation.py b/tests/test_temporal_multilevel_parameter_estimation.py index b2867e3..e98f57d 100644 --- a/tests/test_temporal_multilevel_parameter_estimation.py +++ b/tests/test_temporal_multilevel_parameter_estimation.py @@ -6,39 +6,34 @@ from modularitypruning.leiden_utilities import repeated_leiden_from_gammas_omegas from modularitypruning.parameter_estimation import iterative_multilayer_resolution_parameter_estimation from modularitypruning.parameter_estimation_utilities import gamma_omega_estimate -from modularitypruning.partition_utilities import num_communities, all_degrees +from modularitypruning.partition_utilities import num_communities, all_degrees, in_degrees, out_degrees from numpy import mean from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_random_partition, \ - generate_connected_multilayer_intralayer_SBM + generate_temporal_SBM class TestTemporalAndMultilevelParameterEstimation(unittest.TestCase): - def generate_temporal_SBM(self, copying_probability, p_in, p_out, first_layer_membership, num_layers, - weighted=False): - G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM( - copying_probability, p_in, p_out, first_layer_membership, num_layers, weighted=weighted - ) - - # connect each node to itself in the next layer - num_nodes_per_layer = len(first_layer_membership) - interlayer_edges = [(num_nodes_per_layer * layer + v, num_nodes_per_layer * layer + v + num_nodes_per_layer) - for layer in range(num_layers - 1) for v in range(num_nodes_per_layer)] - G_interlayer = ig.Graph(interlayer_edges, directed=True) - - return G_intralayer, G_interlayer, layer_membership - def assert_temporal_SBM_correct_convergence(self, first_layer_membership, copying_probability=0.75, num_layers=12, - p_in=0.25, p_out=0.05, tolerance=0.1, weighted=False): + p_in=0.25, p_out=0.05, tolerance=0.1, weighted=False, + directed=False): K = num_communities(first_layer_membership) - G_intralayer, G_interlayer, layer_membership = self.generate_temporal_SBM(copying_probability, p_in, p_out, - first_layer_membership, - num_layers, weighted=weighted) + G_intralayer, G_interlayer, layer_membership = generate_temporal_SBM(copying_probability, p_in, p_out, + first_layer_membership, + num_layers, weighted=weighted, + directed=directed) # compute ground truth gamma - k = mean(all_degrees(G_intralayer)) - true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers - true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers + # TODO: Why not strengths for weighted graphs? Confirm this. + if directed: + k_in = mean(in_degrees(G_intralayer)) + k_out = mean(out_degrees(G_intralayer)) + true_theta_in = p_in * G_intralayer.ecount() / (k_in * k_out) / num_layers + true_theta_out = p_out * G_intralayer.ecount() / (k_in * k_out) / num_layers + else: + k = mean(all_degrees(G_intralayer)) + true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers + true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers true_gamma = (true_theta_in - true_theta_out) / (log(true_theta_in) - log(true_theta_out)) # compute ground truth omega. For some reason, Pamfil et al. scale this by 1/2 (perhaps due to the directedness @@ -140,13 +135,48 @@ def test_temporal_SBM_correct_convergence_varying_num_layers(self): ) ) + def test_directed_temporal_SBM_correct_convergence(self): + """Directed analogue of the temporal/multilevel convergence tests.""" + for weighted in [False, True]: + with self.subTest(weighted=weighted): + self.retry_test_across_graph_sizes( + test_fn=lambda n, w=weighted: self.assert_temporal_SBM_correct_convergence( + copying_probability=0.75, num_layers=8, + first_layer_membership=generate_random_partition(num_nodes=n, K=3), weighted=w, + directed=True) + ) + + def test_directed_asymmetry_temporal(self): + """Simple test to confirm multilayer parameter estimates depend on edge direction. + + The intralayer graph has all edges pointing from the first community to the second.""" + one_way_edges = [(0, 1), (2, 3), (0, 2), (0, 3), (1, 2), (1, 3)] + edge_list = one_way_edges + [(u + 4, v + 4) for (u, v) in one_way_edges] + weights = [1.0, 2.0, 3.0, 3.0, 3.0, 3.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0] + membership = [0, 0, 1, 1, 0, 0, 1, 1] + layer_membership = [0, 0, 0, 0, 1, 1, 1, 1] + G_interlayer = ig.Graph([(0, 4), (1, 5), (2, 6), (3, 7)], directed=True) + + for model in ['temporal', 'multilevel']: + with self.subTest(model=model): + G = ig.Graph(8, edge_list, directed=True) + G.es['weight'] = weights + gamma_directed, _ = gamma_omega_estimate(G, G_interlayer, layer_membership, membership, model=model) + + # undirected version for comparison + G_undirected = ig.Graph(8, edge_list, directed=False) + G_undirected.es['weight'] = weights + gamma_undirected_treatment, _ = gamma_omega_estimate(G_undirected, G_interlayer, layer_membership, + membership, model=model) + self.assertGreater(abs(gamma_directed - gamma_undirected_treatment), 1e-2) + def test_directed_consistency_temporal_SBM_leiden(self): """Test parameter estimate consistency on a temporal SBM when the intralayer edges are directed.""" membership = [0] * 25 + [1] * 25 + [2] * 25 - G_intralayer, G_interlayer, layer_membership = self.generate_temporal_SBM(copying_probability=0.9, - p_in=0.25, p_out=0.05, - first_layer_membership=membership, - num_layers=25) + G_intralayer, G_interlayer, layer_membership = generate_temporal_SBM(copying_probability=0.9, + p_in=0.25, p_out=0.05, + first_layer_membership=membership, + num_layers=25) partitions = repeated_leiden_from_gammas_omegas(G_intralayer, G_interlayer, layer_membership, gammas=[0.5, 1.0, 1.5], omegas=[0.5, 1.0, 1.5]) From 6d3eedf3c41c46b6c8bc841fd62a1142fe070a5b Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 21:01:24 -0400 Subject: [PATCH 06/13] Add directed multiplex parameter estimation tests --- tests/test_multiplex_parameter_estimation.py | 85 +++++++++++++------- 1 file changed, 54 insertions(+), 31 deletions(-) diff --git a/tests/test_multiplex_parameter_estimation.py b/tests/test_multiplex_parameter_estimation.py index 77cb040..1002b01 100644 --- a/tests/test_multiplex_parameter_estimation.py +++ b/tests/test_multiplex_parameter_estimation.py @@ -6,41 +6,32 @@ from modularitypruning.leiden_utilities import repeated_leiden_from_gammas_omegas from modularitypruning.parameter_estimation import iterative_multilayer_resolution_parameter_estimation from modularitypruning.parameter_estimation_utilities import gamma_omega_estimate -from modularitypruning.partition_utilities import num_communities, all_degrees +from modularitypruning.partition_utilities import num_communities, all_degrees, in_degrees, out_degrees from numpy import mean from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_random_partition, \ - generate_connected_multilayer_intralayer_SBM + generate_multiplex_SBM class TestMultiplexParameterEstimation(unittest.TestCase): - def generate_multiplex_SBM(self, copying_probability, p_in, p_out, first_layer_membership, num_layers, - weighted=False): - G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM( - copying_probability, p_in, p_out, first_layer_membership, num_layers, weighted=weighted - ) - - # connect each node to itself in all other layers - num_nodes_per_layer = len(first_layer_membership) - interlayer_edges = [(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)] - G_interlayer = ig.Graph(interlayer_edges, directed=True) - - return G_intralayer, G_interlayer, layer_membership - def assert_multiplex_SBM_correct_convergence(self, first_layer_membership, copying_probability=0.75, num_layers=10, - p_in=0.25, p_out=0.05, weighted=False): + p_in=0.25, p_out=0.05, weighted=False, directed=False): K = num_communities(first_layer_membership) - G_intralayer, G_interlayer, layer_membership = self.generate_multiplex_SBM(copying_probability, p_in, p_out, - first_layer_membership, - num_layers, weighted=weighted) + G_intralayer, G_interlayer, layer_membership = generate_multiplex_SBM(copying_probability, p_in, p_out, + first_layer_membership, + num_layers, weighted=weighted, + directed=directed) # compute ground truth gamma - k = mean(all_degrees(G_intralayer)) - true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers - true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers + if directed: + k_in = mean(in_degrees(G_intralayer)) + k_out = mean(out_degrees(G_intralayer)) + true_theta_in = p_in * G_intralayer.ecount() / (k_in * k_out) / num_layers + true_theta_out = p_out * G_intralayer.ecount() / (k_in * k_out) / num_layers + else: + k = mean(all_degrees(G_intralayer)) + true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers + true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers true_gamma = (true_theta_in - true_theta_out) / (log(true_theta_in) - log(true_theta_out)) # compute ground truth omega @@ -100,19 +91,51 @@ def test_multiplex_SBM_correct_convergence_varying_num_layers(self): self.assert_multiplex_SBM_correct_convergence(first_layer_membership=membership, num_layers=num_layers, weighted=weighted) + def test_directed_multiplex_SBM_correct_convergence(self): + """Directed analogue of the multiplex convergence tests.""" + for weighted in [False, True]: + with self.subTest(weighted=weighted): + membership = generate_random_partition(num_nodes=300, K=3) + self.assert_multiplex_SBM_correct_convergence(copying_probability=0.75, num_layers=8, + first_layer_membership=membership, weighted=weighted, + directed=True) + + def test_directed_asymmetry_multiplex(self): + """Simple test to confirm multilayer parameter estimates depend on edge direction. + + The intralayer graph has all edges pointing from the first community to the second.""" + one_way_edges = [(0, 1), (2, 3), (0, 2), (0, 3), (1, 2), (1, 3)] + edge_list = one_way_edges + [(u + 4, v + 4) for (u, v) in one_way_edges] + weights = [1.0, 2.0, 3.0, 3.0, 3.0, 3.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0] + membership = [0, 0, 1, 1, 0, 0, 1, 1] + layer_membership = [0, 0, 0, 0, 1, 1, 1, 1] + G_interlayer = ig.Graph([(0, 4), (4, 0), (1, 5), (5, 1), (2, 6), (6, 2), (3, 7), (7, 3)], directed=True) + + with self.subTest(model='multiplex'): + G = ig.Graph(8, edge_list, directed=True) + G.es['weight'] = weights + gamma_directed, _ = gamma_omega_estimate(G, G_interlayer, layer_membership, membership, model='multiplex') + + # undirected version for comparison + G_undirected = ig.Graph(8, edge_list, directed=False) + G_undirected.es['weight'] = weights + gamma_undirected_treatment, _ = gamma_omega_estimate(G_undirected, G_interlayer, layer_membership, + membership, model='multiplex') + self.assertGreater(abs(gamma_directed - gamma_undirected_treatment), 1e-2) + def test_directed_consistency_multiplex_SBM_leiden(self): """Test parameter estimate consistency on a multiplex SBM when the intralayer edges are directed.""" membership = [0] * 25 + [1] * 25 + [2] * 25 - G_intralayer, G_interlayer, layer_membership = self.generate_multiplex_SBM(copying_probability=0.9, - p_in=0.25, p_out=0.05, - first_layer_membership=membership, - num_layers=25) + G_intralayer, G_interlayer, layer_membership = generate_multiplex_SBM(copying_probability=0.9, + p_in=0.25, p_out=0.05, + first_layer_membership=membership, + num_layers=25) partitions = repeated_leiden_from_gammas_omegas(G_intralayer, G_interlayer, layer_membership, gammas=[0.5, 1.0, 1.5], omegas=[0.5, 1.0, 1.5]) - for partition in partitions: - with self.subTest(partition=partition): + for i, partition in enumerate(partitions): + with self.subTest(partition=i): # here, undirected/directed refers to the intralayer edges only # in Pamfil et al.'s multiplex networks, interlayer edges are taken to be directed gamma_undirected, omega_undirected = gamma_omega_estimate(G_intralayer, G_interlayer, layer_membership, From c8a86b5645489e01a09d341f380a71adbbe8ba39 Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 21:04:42 -0400 Subject: [PATCH 07/13] Add directed SBM CHAMP tests --- tests/test_champ_coefficients_2D.py | 12 +++++++++++- tests/test_champ_coefficients_3D.py | 24 +++++++++++++++++++++++- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/tests/test_champ_coefficients_2D.py b/tests/test_champ_coefficients_2D.py index e45a176..fa0227e 100644 --- a/tests/test_champ_coefficients_2D.py +++ b/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. diff --git a/tests/test_champ_coefficients_3D.py b/tests/test_champ_coefficients_3D.py index 4c33b1e..a47e9e2 100644 --- a/tests/test_champ_coefficients_3D.py +++ b/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]: From 9c0c782704384f2cf92c148058156051874487c4 Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 21:06:31 -0400 Subject: [PATCH 08/13] Slightly refactor performance tests to new helper functions --- tests/test_parallel_leiden_performance.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/tests/test_parallel_leiden_performance.py b/tests/test_parallel_leiden_performance.py index 7adcc77..a9f3529 100644 --- a/tests/test_parallel_leiden_performance.py +++ b/tests/test_parallel_leiden_performance.py @@ -13,7 +13,8 @@ repeated_leiden_from_gammas_omegas, repeated_parallel_leiden_from_gammas_omegas) -from .shared_testing_functions import generate_connected_ER, generate_connected_multilayer_intralayer_SBM +from .shared_testing_functions import generate_connected_ER, generate_connected_multilayer_intralayer_SBM, \ + generate_temporal_interlayer_edges # this set of tests ensures that we achieve >= 75% parallel performance compared to perfect scaling of # single-threaded jobs to multiple cores (with no memory contention). This threshold will be decreased in @@ -118,9 +119,8 @@ def test_tiny_multilayer_graph_many_runs(self): G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM( copying_probability=0.9, p_in=0.8, p_out=0.2, first_layer_membership=[0] * 5 + [1] * 5, num_layers=5 ) - interlayer_edges = [(10 * layer + v, 10 * layer + v + 10) - for layer in range(5 - 1) for v in range(10)] - G_interlayer = ig.Graph(interlayer_edges, directed=True) + G_interlayer = ig.Graph(generate_temporal_interlayer_edges(num_nodes_per_layer=10, num_layers=5), + directed=True) gammas = np.linspace(0.0, 2.0, 100) omegas = np.linspace(0.0, 2.0, 100) @@ -133,9 +133,8 @@ def test_larger_multilayer_graph_few_runs(self): G_intralayer, layer_membership = generate_connected_multilayer_intralayer_SBM( copying_probability=0.9, p_in=0.15, p_out=0.05, first_layer_membership=[0] * 50 + [1] * 50, num_layers=25 ) - interlayer_edges = [(100 * layer + v, 100 * layer + v + 100) - for layer in range(25 - 1) for v in range(100)] - G_interlayer = ig.Graph(interlayer_edges, directed=True) + G_interlayer = ig.Graph(generate_temporal_interlayer_edges(num_nodes_per_layer=100, num_layers=25), + directed=True) gammas = np.linspace(0.0, 2.0, 7) omegas = np.linspace(0.0, 2.0, 7) From 7f3ea557ffd35795a4b1fc3189e970ef629b6fcc Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 21:47:55 -0400 Subject: [PATCH 09/13] Improve consistency of ground truth parameter test calculations --- tests/test_monolayer_parameter_estimation.py | 19 +++++--------- tests/test_multiplex_parameter_estimation.py | 25 +++++++++++------- ...emporal_multilevel_parameter_estimation.py | 26 ++++++++++++------- 3 files changed, 39 insertions(+), 31 deletions(-) diff --git a/tests/test_monolayer_parameter_estimation.py b/tests/test_monolayer_parameter_estimation.py index 25d3713..a39af6d 100644 --- a/tests/test_monolayer_parameter_estimation.py +++ b/tests/test_monolayer_parameter_estimation.py @@ -5,7 +5,7 @@ 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, in_strengths, out_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_SBM, \ @@ -28,17 +28,12 @@ def assert_newman_synthetic_convergence(self, q, weighted=False): G, p_in, p_out = generate_SBM(q, weighted=weighted) # compute "ground truth" gamma estimate - # TODO: Why mean weight scaling for weighted graphs? Confirm this. - 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) diff --git a/tests/test_multiplex_parameter_estimation.py b/tests/test_multiplex_parameter_estimation.py index 1002b01..1e2750d 100644 --- a/tests/test_multiplex_parameter_estimation.py +++ b/tests/test_multiplex_parameter_estimation.py @@ -6,7 +6,7 @@ from modularitypruning.leiden_utilities import repeated_leiden_from_gammas_omegas from modularitypruning.parameter_estimation import iterative_multilayer_resolution_parameter_estimation from modularitypruning.parameter_estimation_utilities import gamma_omega_estimate -from modularitypruning.partition_utilities import num_communities, all_degrees, in_degrees, out_degrees +from modularitypruning.partition_utilities import num_communities, all_strengths, in_strengths, out_strengths from numpy import mean from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_random_partition, \ @@ -22,16 +22,23 @@ def assert_multiplex_SBM_correct_convergence(self, first_layer_membership, copyi num_layers, weighted=weighted, directed=directed) - # compute ground truth gamma + # compute "ground truth" gamma estimate + # note that strengths are just degrees in unweighted graphs + if weighted: + total_weight = sum(G_intralayer.es['weight']) + mean_weight = total_weight / G_intralayer.ecount() + else: + total_weight = G_intralayer.ecount() + mean_weight = 1.0 + twom = total_weight if directed else 2 * total_weight + if directed: - k_in = mean(in_degrees(G_intralayer)) - k_out = mean(out_degrees(G_intralayer)) - true_theta_in = p_in * G_intralayer.ecount() / (k_in * k_out) / num_layers - true_theta_out = p_out * G_intralayer.ecount() / (k_in * k_out) / num_layers + k_in = mean(in_strengths(G_intralayer)) + k_out = mean(out_strengths(G_intralayer)) else: - k = mean(all_degrees(G_intralayer)) - true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers - true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers + k_in = k_out = mean(all_strengths(G_intralayer)) + true_theta_in = mean_weight * p_in * twom / (k_in * k_out) / num_layers + true_theta_out = mean_weight * p_out * twom / (k_in * k_out) / num_layers true_gamma = (true_theta_in - true_theta_out) / (log(true_theta_in) - log(true_theta_out)) # compute ground truth omega diff --git a/tests/test_temporal_multilevel_parameter_estimation.py b/tests/test_temporal_multilevel_parameter_estimation.py index e98f57d..1b6a640 100644 --- a/tests/test_temporal_multilevel_parameter_estimation.py +++ b/tests/test_temporal_multilevel_parameter_estimation.py @@ -6,7 +6,7 @@ from modularitypruning.leiden_utilities import repeated_leiden_from_gammas_omegas from modularitypruning.parameter_estimation import iterative_multilayer_resolution_parameter_estimation from modularitypruning.parameter_estimation_utilities import gamma_omega_estimate -from modularitypruning.partition_utilities import num_communities, all_degrees, in_degrees, out_degrees +from modularitypruning.partition_utilities import num_communities, all_strengths, in_strengths, out_strengths from numpy import mean from .shared_testing_functions import assert_almost_equal_or_both_none_or_nan, generate_random_partition, \ @@ -23,17 +23,23 @@ def assert_temporal_SBM_correct_convergence(self, first_layer_membership, copyin num_layers, weighted=weighted, directed=directed) - # compute ground truth gamma - # TODO: Why not strengths for weighted graphs? Confirm this. + # compute "ground truth" gamma estimate + # note that strengths are just degrees in unweighted graphs + if weighted: + total_weight = sum(G_intralayer.es['weight']) + mean_weight = total_weight / G_intralayer.ecount() + else: + total_weight = G_intralayer.ecount() + mean_weight = 1.0 + twom = total_weight if directed else 2 * total_weight + if directed: - k_in = mean(in_degrees(G_intralayer)) - k_out = mean(out_degrees(G_intralayer)) - true_theta_in = p_in * G_intralayer.ecount() / (k_in * k_out) / num_layers - true_theta_out = p_out * G_intralayer.ecount() / (k_in * k_out) / num_layers + k_in = mean(in_strengths(G_intralayer)) + k_out = mean(out_strengths(G_intralayer)) else: - k = mean(all_degrees(G_intralayer)) - true_theta_in = p_in * (2 * G_intralayer.ecount()) / (k * k) / num_layers - true_theta_out = p_out * (2 * G_intralayer.ecount()) / (k * k) / num_layers + k_in = k_out = mean(all_strengths(G_intralayer)) + true_theta_in = mean_weight * p_in * twom / (k_in * k_out) / num_layers + true_theta_out = mean_weight * p_out * twom / (k_in * k_out) / num_layers true_gamma = (true_theta_in - true_theta_out) / (log(true_theta_in) - log(true_theta_out)) # compute ground truth omega. For some reason, Pamfil et al. scale this by 1/2 (perhaps due to the directedness From 715eac5a8a97b4c3ad33dc79512a5ef64aa68f83 Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 22:14:04 -0400 Subject: [PATCH 10/13] Add documentation notes on new directed+weighted support --- README.md | 3 +++ docs/index.rst | 3 +++ utilities/champ_utilities.py | 4 ++++ utilities/leiden_utilities.py | 4 ++++ utilities/parameter_estimation_utilities.py | 11 +++++++++-- 5 files changed, 23 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index f9d680b..a2eefb7 100644 --- a/README.md +++ b/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 diff --git a/docs/index.rst b/docs/index.rst index bd9c2cd..2738fc9 100644 --- a/docs/index.rst +++ b/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 diff --git a/utilities/champ_utilities.py b/utilities/champ_utilities.py index a135bef..d558469 100644 --- a/utilities/champ_utilities.py +++ b/utilities/champ_utilities.py @@ -107,6 +107,8 @@ def CHAMP_2D( See https://doi.org/10.3390/a10030093 for more details. + Directed, undirected, weighted, and unweighted graphs are supported. + :param G: graph of interest :type G: igraph.Graph :param all_parts: partitions to prune @@ -176,6 +178,8 @@ def CHAMP_3D( See https://doi.org/10.3390/a10030093 for more details. + Both the intralayer and interlayer graphs may be directed or undirected, weighted or unweighted. + :param G_intralayer: intralayer graph of interest :type G_intralayer: igraph.Graph :param G_interlayer: interlayer graph of interest diff --git a/utilities/leiden_utilities.py b/utilities/leiden_utilities.py index 751863e..a3dd3be 100644 --- a/utilities/leiden_utilities.py +++ b/utilities/leiden_utilities.py @@ -51,6 +51,8 @@ def singlelayer_leiden( ) -> SingleLayerLeidenResult: r"""Run the Leiden modularity maximization algorithm at a single :math:`\gamma` value. + Directed, undirected, weighted, and unweighted graphs are supported. + :param G: graph of interest :type G: igraph.Graph :param gamma: gamma (resolution parameter) to run Leiden at @@ -157,6 +159,8 @@ def multilayer_leiden( ) -> MultilayerLeidenResult: r"""Run the Leiden modularity maximization algorithm at a single (:math:`\gamma, \omega`) value. + Both the intralayer and interlayer graphs may be directed or undirected, weighted or unweighted. + WARNING: Optimization can be EXTREMELY slow for large numbers of layers! Leidenalg does not properly implement multilayer optimization. diff --git a/utilities/parameter_estimation_utilities.py b/utilities/parameter_estimation_utilities.py index 7ac87c6..f5685bd 100644 --- a/utilities/parameter_estimation_utilities.py +++ b/utilities/parameter_estimation_utilities.py @@ -85,6 +85,9 @@ def estimate_multilayer_SBM_parameters( ) -> MultilayerSBMParameters: """Estimates multilayer SBM parameters from a graph and a partition + For directed intralayer graphs, the convention of Leicht and Newman is used per layer (as in + :meth:`~modularitypruning.parameter_estimation_utilities.estimate_singlelayer_SBM_parameters`). + :param G_intralayer: intralayer graph of interest :type G_intralayer: igraph.Graph :param G_interlayer: interlayer graph of interest @@ -203,6 +206,9 @@ def gamma_estimate(G: ig.Graph, partition: LeidenPartition | Membership) -> Gamm See https://doi.org/10.1103/PhysRevE.94.052315 for more details. + For directed graphs, the convention of Leicht and Newman is used (see + :meth:`~modularitypruning.parameter_estimation_utilities.estimate_singlelayer_SBM_parameters`). + :param G: graph of interest :type G: igraph.Graph :param partition: partition of interest @@ -597,7 +603,7 @@ def prune_to_stable_partitions( Runs our full pruning pipeline on a singlelayer network. Returns the pruned list of stable partitions whose parameter estimates are within the provided ``gamma_start`` and ``gamma_end`` bounds. - Weighted graphs are supported via the ``weight`` edge attribute. Edge weights must be non-negative. + Directed, undirected, weighted, and unweighted graphs are supported. Edge weights must be non-negative. See https://doi.org/10.1038/s41598-022-20142-6 for more details. @@ -661,7 +667,8 @@ def prune_to_multilayer_stable_partitions( parameter estimates are within the provided ``gamma_start``, ``gamma_end``, ``omega_start``, and ``omega_end`` bounds. - Weighted graphs are supported via the ``weight`` edge attribute. Edge weights must be non-negative. + Directed, undirected, weighted, and unweighted graphs are supported. Interlayer graphs must be directed. Edge + weights must be non-negative. Interlayer weights are not fully supported by the underlying multilayer SBM theory (uniform interlayer coupling). As such, they are used in the CHAMP step, but parameter estimation treats interlayer edges as unweighted. The code From 8ec9269d8f74cf0bfad1b8ef8259f53765b04f5d Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 22:18:48 -0400 Subject: [PATCH 11/13] Make one warning in prune_to_multilayer_stable_partitions more precise --- utilities/parameter_estimation_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/parameter_estimation_utilities.py b/utilities/parameter_estimation_utilities.py index f5685bd..ebe5369 100644 --- a/utilities/parameter_estimation_utilities.py +++ b/utilities/parameter_estimation_utilities.py @@ -739,7 +739,7 @@ def prune_to_multilayer_stable_partitions( if G_interlayer.is_weighted() and any(x != 1.0 for x in G_interlayer.es['weight']): warnings.warn("Parameter estimation treats interlayer edges as unweighted and will proceed as though all " - "interlayer weights are 1.0.") + "interlayer weights are 1.0, but CHAMP will still use the actual interlayer weights!") if isinstance(parts, leidenalg.RBConfigurationVertexPartition): # convert to (canonically represented) membership vectors if necessary From d53984b5f0f4804ca0c42aabb15a90ce1e37905d Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sat, 17 Oct 2026 22:28:03 -0400 Subject: [PATCH 12/13] Add new directedness note to original Lazega experiment --- experiments/lazega_law_firm/lazega_figures.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/experiments/lazega_law_firm/lazega_figures.py b/experiments/lazega_law_firm/lazega_figures.py index 1964bef..8819fcf 100644 --- a/experiments/lazega_law_firm/lazega_figures.py +++ b/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 From ea9a4407e20a3848dba3cd4307fefa5c1f6a611a Mon Sep 17 00:00:00 2001 From: Ryan Gibson Date: Sun, 18 Oct 2026 00:01:00 -0400 Subject: [PATCH 13/13] Improve robustness of degenerate partition guards --- utilities/parameter_estimation_utilities.py | 40 +++++++++++++-------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/utilities/parameter_estimation_utilities.py b/utilities/parameter_estimation_utilities.py index ebe5369..063a1de 100644 --- a/utilities/parameter_estimation_utilities.py +++ b/utilities/parameter_estimation_utilities.py @@ -58,20 +58,23 @@ def estimate_singlelayer_SBM_parameters( kappa_out_list[community[e.source]] += e['weight'] kappa_in_list[community[e.target]] += e['weight'] sum_kappa_in_kappa_out: float = sum(k_in * k_out for k_in, k_out in zip(kappa_in_list, kappa_out_list)) + denominator: float = sum_kappa_in_kappa_out / m - omega_in: float = m_in / (sum_kappa_in_kappa_out / m) - # guard for div by zero with single community partition - omega_out: float = (m - m_in) / (m - sum_kappa_in_kappa_out / m) if len(partition) > 1 else 0 + # guard for div by zero with degenerate partitions (e.g., a single community, communities with only incoming + # or only outgoing edges) + omega_in: float = m_in / denominator if denominator > 0 else 0 + omega_out: float = (m - m_in) / (m - denominator) if denominator < m else 0 else: kappa_r_list: list[float] = [0] * len(partition) for e in G.es: kappa_r_list[community[e.source]] += e['weight'] kappa_r_list[community[e.target]] += e['weight'] sum_kappa_sqr: float = sum(x ** 2 for x in kappa_r_list) + denominator = sum_kappa_sqr / (2 * m) - omega_in = (2 * m_in) / (sum_kappa_sqr / (2 * m)) - # guard for div by zero with single community partition - omega_out = (2 * m - 2 * m_in) / (2 * m - sum_kappa_sqr / (2 * m)) if len(partition) > 1 else 0 + # guard for div by zero with a degenerate partition (e.g., a single community) + omega_in = (2 * m_in) / denominator if denominator > 0 else 0 + omega_out = (2 * m - 2 * m_in) / (2 * m - denominator) if denominator < 2 * m else 0 # return estimates for omega_in, omega_out return omega_in, omega_out @@ -138,8 +141,6 @@ def estimate_multilayer_SBM_parameters( if community[e.source] == community[e.target] and layer_vec[e.source] == layer_vec[e.target]: m_t_in[layer_vec[e.source]] += e['weight'] - kappa_t_r_list: list[list[int]] = [[0] * K for _ in range(T)] - if G_intralayer.is_directed(): # directed modularity per layer as in estimate_singlelayer_SBM_parameters kappa_t_in_list: list[list[float]] = [[0] * K for _ in range(T)] @@ -152,26 +153,27 @@ def estimate_multilayer_SBM_parameters( sum(k_in * k_out for k_in, k_out in zip(kappa_t_in_list[t], kappa_t_out_list[t])) for t in range(T)] - theta_in: float = sum(m_t_in[t] for t in range(T)) / sum(sum_kappa_in_kappa_out_t[t] / m_t[t] - for t in range(T)) + # guard for div by zero with degenerate partitions (e.g., communities with only incoming + # or only outgoing edges) + theta_in_denominator: float = sum(sum_kappa_in_kappa_out_t[t] / m_t[t] for t in range(T)) + theta_in: float = sum(m_t_in[t] for t in range(T)) / theta_in_denominator if theta_in_denominator > 0 else 0 - # guard for div by zero with e.g. a single community partition theta_out_numerator: float = sum(m_t[t] - m_t_in[t] for t in range(T)) theta_out_denominator: float = sum(m_t[t] - sum_kappa_in_kappa_out_t[t] / m_t[t] for t in range(T)) else: + kappa_t_r_list: list[list[int]] = [[0] * K for _ in range(T)] for e in G_intralayer.es: layer = layer_vec[e.source] kappa_t_r_list[layer][community[e.source]] += e['weight'] kappa_t_r_list[layer][community[e.target]] += e['weight'] sum_kappa_t_sqr: list[float] = [sum(x ** 2 for x in kappa_t_r_list[t]) for t in range(T)] + theta_in_denominator = sum(sum_kappa_t_sqr[t] / (2 * m_t[t]) for t in range(T)) + theta_in = sum(2 * m_t_in[t] for t in range(T)) / theta_in_denominator if theta_in_denominator > 0 else 0 - theta_in = sum(2 * m_t_in[t] for t in range(T)) / sum(sum_kappa_t_sqr[t] / (2 * m_t[t]) - for t in range(T)) - - # guard for div by zero with e.g. a single community partition theta_out_numerator = sum(2 * m_t[t] - 2 * m_t_in[t] for t in range(T)) theta_out_denominator = sum(2 * m_t[t] - sum_kappa_t_sqr[t] / (2 * m_t[t]) for t in range(T)) + # guard for div by zero with a degenerate partition (e.g., a single community) if theta_out_denominator == 0: theta_out: float = 0 else: @@ -271,6 +273,10 @@ def multiplex_omega_estimate_from_parameters( if p == 0: return 0 + # guard for degenerate partition with no within-community edges + if theta_in == 0: + return 0 + # if p is 1, the optimal omega is infinite (here, omega_max) if p >= 1.0 or theta_in == 1.0: return omega_max @@ -303,6 +309,10 @@ def temporal_multilevel_omega_estimate_from_parameters( if p == 0: return 0 + # guard for degenerate partition with no within-community edges + if theta_in == 0: + return 0 + if theta_out == 0: return log(1 + p * K / (1 - p)) / (2 * log(theta_in)) if p < 1.0 else omega_max # if p is 1, the optimal omega is infinite (here, omega_max)