Infomap 的目标函数——两层映射方程（map equation）

# 对无向无权网络，随机游走的稳态访问概率 p_α = k_α / 2m；模块 i 的退出概率 q_i = (i 的外连边数)/2m。
# 两层映射方程（Rosvall 与 Bergstrom 2008）给出描述一步游走所需的平均码长（bit）：
#   L(M) = plogp(Σq_i) - 2 Σ plogp(q_i) - Σ_α plogp(p_α) + Σ_i plogp(q_i + Σ_{α∈i} p_α)，plogp(x)=x·log2 x。
# (a) 比较空手道网络的几种划分（单一模块、真实两派、Louvain 划分）的 L；
# (b) 在“K5 团组成的环”上比较“每个团一个模块”与“相邻两团合并”两种划分的 Q 与 L。
# 输出：(a) L 越小代表划分越能压缩游走的描述；(b) 团足够多时 Q 偏好合并（分辨率极限），
#       而映射方程在这一构造下仍偏好把每个团分开。
import networkx as nx
import numpy as np


def plogp(x):
    x = np.asarray(x, float)
    return np.where(x > 0, x * np.log2(np.where(x > 0, x, 1)), 0.0)


def map_equation(G, communities):
    m2 = 2 * G.number_of_edges()
    label = {v: i for i, c in enumerate(communities) for v in c}
    p = np.array([G.degree(v) / m2 for v in G])
    q = np.zeros(len(communities))
    pin = np.zeros(len(communities))
    for v in G:
        pin[label[v]] += G.degree(v) / m2
    for u, v in G.edges():
        if label[u] != label[v]:
            q[label[u]] += 1 / m2
            q[label[v]] += 1 / m2
    return float(plogp(q.sum()) - 2 * plogp(q).sum() - plogp(p).sum() + plogp(q + pin).sum())


K = nx.karate_club_graph()
G = nx.Graph(K.edges())
club = [{v for v in K if K.nodes[v]["club"] == "Mr. Hi"}, {v for v in K if K.nodes[v]["club"] != "Mr. Hi"}]
louv = nx.community.louvain_communities(G, seed=1)
for name, part in [("单一模块", [set(G)]), ("真实两派", club), ("Louvain(%d 个)" % len(louv), louv)]:
    print("(a) %-12s Q = %.3f   L = %.3f bit" % (name, nx.community.modularity(G, part), map_equation(G, part)))

ring = []                                           # 存下结果用于作图
for n_c in [10, 30, 60]:
    R = nx.ring_of_cliques(n_c, 5)                  # 第 j 个团为节点 5j..5j+4
    single = [set(range(5 * j, 5 * j + 5)) for j in range(n_c)]
    pairs = [single[j] | single[j + 1] for j in range(0, n_c, 2)]
    print("(b) %2d 个 K5: 每团一模块 Q=%.4f L=%.3f | 两两合并 Q=%.4f L=%.3f"
          % (n_c, nx.community.modularity(R, single), map_equation(R, single),
             nx.community.modularity(R, pairs), map_equation(R, pairs)))
    ring.append((n_c, nx.community.modularity(R, single), nx.community.modularity(R, pairs),
                 map_equation(R, single), map_equation(R, pairs)))

# 作图：团环上“每团一模块”与“两两合并”的 Q（越大越好）和 L（越小越好）
import matplotlib.pyplot as plt

r = np.array(ring)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
for ax, j, lab in [(ax1, 1, "modularity Q (higher is better)"), (ax2, 3, "map equation L, bits (lower is better)")]:
    ax.plot(r[:, 0], r[:, j], "o-", label="one module per clique")
    ax.plot(r[:, 0], r[:, j + 1], "s--", label="adjacent cliques merged")
    ax.set_xlabel("number of K5 cliques in the ring"); ax.set_ylabel(lab); ax.legend(fontsize=8)
ax1.set_title("Q prefers merging for many cliques (resolution limit)", fontsize=10)
ax2.set_title("Map equation keeps cliques separate", fontsize=10)
plt.tight_layout(); plt.show()