基于graph-tool的图模块化实现结果与论文数值不符的问题排查求助

我现在正试着自己实现一种图的模块化定义，用Python的graph-tool库构建单层图结构，还照着一篇属性图聚类相关论文里的例子（图4）搭建了几个测试图。但碰到了个棘手的问题：我算出的模块化数值和论文里报告的对不上——我得到的是0.42，可论文里写的是0.28。我翻来覆去检查了代码，还是找不到问题出在哪，想请大家帮忙排查一下！

我的单图构建方法 #

下面是我用来构建单层图的函数，支持有向/无向图和边颜色设置：

def singlelayer_graph(edges: Edges, directed: bool = True, edge_colors: dict[tuple[int, int]] = None) -> gt.Graph:
    """
    Initializes a single-layer directed graph with optional edge colors.

    Args:
        edges: A list of tuples of the form (source_id, target_id, weight).
        directed: Optional. A boolean flag that defines the directionality of edges (directed by default)
        edge_colors: Optional. A dictionary where keys are (source, target) tuples and values
                     are RGB color lists for the edges.

    Returns:
        A gt.Graph object with the following properties:
        - original_id (vertex property): The original node ID mapping.
        - weight (edge property): The weight of each edge.
        - edge_color (edge property, if provided): The color of each edge.
    """
    G = gt.Graph(directed=directed)

    # Define graph properties
    G.gp["num_edges"] = G.new_graph_property("int")
    G.gp["num_nodes"] = G.new_graph_property("int")

    G.gp["num_edges"] = len(edges)
    num_nodes = len(set(node for edge in edges for node in edge[:2]))
    G.gp["num_nodes"] = num_nodes

    # Define vertex properties
    G.vp["original_id"] = G.new_vertex_property("int")  # Store original node ID mapping.
    G.ep["weight"] = G.new_edge_property("int")

    # Add edge color property if edge_colors is provided
    if edge_colors:
        G.ep["edge_color"] = G.new_edge_property("vector<double>")

    # Create a mapping from original node IDs to graph vertices
    node_map = {}

    for source, target, weight in edges:
        # Ensure both source and target nodes exist in the graph
        if source not in node_map:
            v_source = G.add_vertex()
            G.vp["original_id"][v_source] = source
            node_map[source] = v_source

        if target not in node_map:
            v_target = G.add_vertex()
            G.vp["original_id"][v_target] = target
            node_map[target] = v_target

        # Add edge using mapped vertices
        edge = G.add_edge(node_map[source], node_map[target])
        G.ep["weight"][edge] = weight  # Assign weight

        # Assign edge color if edge_colors is provided
        if edge_colors and (source, target) in edge_colors:
            G.ep["edge_color"][edge] = edge_colors[(source, target)]  # Assign color

    return G

测试图的构建代码 #

这是我用来复现论文中例子的具体图构建函数，还预设了社区划分：

def bothorel_graph_3():

    es = [(1, 2, 1),
          (1, 3, 1),
          (2, 3, 1),
          (3, 4, 1),
          (2, 4, 1),
          (4, 5, 1),
          (5, 6, 1),
          (4, 7, 1),
          (6, 7, 1),
          (7, 8, 1)]

    # Create graph
    G = singlelayer_graph(edges=es, directed=False)

    # Assign communities by (original) node ID
    communities = {
        1: 0, 2: 0, 3: 0,
        4: 1, 5: 1, 6: 1, 7: 1, 8: 1
    }

    # Create a new vertex property for communities
    G.vp["community"] = G.new_vertex_property("int")

    # Assign communities based on original IDs
    for v in G.vertices():
        G.vp["community"][v] = communities[G.vp["original_id"][v]]

    return G

我实现的模块化计算函数 #

下面是核心的模块化计算逻辑，支持有向和无向图：

from decimal import Decimal, ROUND_HALF_UP

def graph_modularity(graph: gt.Graph) -> float:
    """
    Computes the modularity of the input graph (supports both directed and undirected graphs).

    Args:
        graph (gt.Graph): The graph-tool Graph object.

    Returns:
        float: The modularity score of the graph.
    """
    # Ensure the graph has community assignments
    if "community" not in graph.vp:
        raise ValueError("Community assignment missing. Ensure 'community' vertex property is set.")

    # Check if the graph is directed
    is_directed = graph.is_directed()

    print(f"Graph is directed: {is_directed}")

    # Set m to the sum of weights (total number of edges in an undirected graph)
    # Convert edge weight sum to Decimal
    m = Decimal(sum(Decimal(graph.ep["weight"][e]) for e in graph.edges()))

    print(f"Total edge weight (m): {m}")

    modularity = Decimal(0)

    # List of all vertices and their corresponding communities
    all_vertices = list(graph.vertices())

    # Compute modularity contribution for all pairs (i, j)
    for i in range(len(all_vertices)):
        for j in range(len(all_vertices)):  # Consider all pairs for directed graphs

            if i == j:
                continue  # Skip self-loops

            vi, vj = all_vertices[i], all_vertices[j]
            ci, cj = graph.vp["community"][vi], graph.vp["community"][vj]

            # Check if there is an edge between vi and vj
            edge = graph.edge(vi, vj)
            A_ij = Decimal(graph.ep["weight"][edge]) if edge is not None else Decimal(0)  # Edge weight or 0 if no edge

            # Compute degree based on directed or undirected graph
            if is_directed:
                k_i_out = sum(Decimal(graph.ep["weight"][e]) for e in vi.out_edges())  # Out-degree
                k_j_in = sum(Decimal(graph.ep["weight"][e]) for e in vj.in_edges())    # In-degree
                expected_weight = (k_i_out * k_j_in) / m
            else:
                k_i = sum(Decimal(graph.ep["weight"][e]) for e in vi.all_edges())
                k_j = sum(Decimal(graph.ep["weight"][e]) for e in vj.all_edges())
                expected_weight = (k_i * k_j) / (2 * m)

            if ci == cj:
                contribution = A_ij - expected_weight
            else:
                contribution = Decimal(0)

            modularity += contribution

            # Debugging prints for each pair (optional)
            # Print detailed debug info
            print(f"Pair ({graph.vp['original_id'][vi]}, {graph.vp['original_id'][vj]}):")
            print(f"  Community of {graph.vp['original_id'][vi]}: {ci}, Community of {graph.vp['original_id'][vj]}: {cj}")
            print(f"  A_ij: {A_ij}, k_i: {k_i}, k_j: {k_j}")
            print(f"  Expected Weight: {expected_weight:.6f}")
            print(f"  Contribution: {contribution:.6f}")
            print(f"  Modularity So Far: {modularity:.6f}\n")
    if is_directed:
        # Normalize and return modularity with precision
        modularity = (modularity / m).quantize(Decimal("0.000001"), rounding=ROUND_HALF_UP)
    else:
        # Normalize and return modularity with precision
        modularity = (modularity / (2 * m)).quantize(Decimal("0.000001"), rounding=ROUND_HALF_UP)

    return modularity

调试输出的部分内容 #

我加了调试打印，下面是部分输出结果，能看到每对节点的贡献计算过程：

Graph is directed: False
Total edge weight (m): 10
Pair (1, 2):
  Community of 1: 0, Community of 2: 0
  A_ij: 1, k_i: 2, k_j: 3
  Expected Weight: 0.300000
  Contribution: 0.700000
  Modularity So Far: 0.700000

Pair (1, 3):
  Community of 1: 0, Community of 3: 0
  A_ij: 1, k_i: 2, k_j: 3
  Expected Weight: 0.300000
  Contribution: 0.700000
  Modularity So Far: 1.400000

Pair (1, 4):
  Community of 1: 0, Community of 4: 1
  A_ij: 0, k_i: 2, k_j: 4
  Expected Weight: 0.400000
  Contribution: 0.000000
  Modularity So Far: 1.400000

Pair (1, 5):
  Community of 1: 0, Community of 5: 1
  A_ij: 0, k_i: 2, k_j: 2
  Expected Weight: 0.200000
  Contribution: 0.000000
  Modularity So Far: 1.400000

Pair (1, 6):
  Community of 1: 0, Community of 6: 1
  A_ij: 0, k_i: 2, k_j: 2
  Expected Weight: 0.200000
  Contribution: 0.000000
  Modularity So Far: 1.400000

Pair (1, 7):
  Community of 1: 0, Community of 7: 1
  A_ij: 0, k_i: 2, k_j: 3
  Expected Weight: 0.300000
  Contribution: 0.000000
  Modularity So Far: 1.400000

Pair (1, 8):
  Community of 1: 0, Community of 8: 1
  A_ij: 0, k_i: 2, k_j: 1
  Expected Weight: 0.100000
  Contribution: 0.000000
  Modularity So Far: 1.400000

Pair (2, 1):
  Community of 2: 0, Community of 1: 0
  A_ij: 1, k_i: 3, k_j: 2
  Expected Weight: 0.300000
  Contribution: 0.700000
  Modularity So Far: 2.100000

Pair (2, 3):
  Community of 2: 0, Community of 3: 0
  A_ij: 1, k_i: 3, k_j: 3
  Expected Weight: 0.450000
  Contribution: 0.550000
  Modularity So Far: 2.650000

Pair (2, 4):
  Community of 2: 0, Community of 4: 1
  A_ij: 1, k_i: 3, k_j: 4
  Expected Weight: 0.600000
  Contribution: 0.000000
  Modularity So Far: 2.650000

...（剩余调试输出略）

备注：内容来源于stack exchange，提问作者Anne Flank