含负权的模块度与功能角色划分（Guimerà–Amaral 功能制图）

# 对含正负权的相关矩阵，用 Rubinov & Sporns (2011) 的非对称模块度
#   Q* = Q+ − [v−/(v+ + v−)] Q−
# 评价给定划分；再对正权网络计算模块内度 z 分数与参与系数 P，按 Guimerà & Amaral (2005) 的阈值
# （z ≥ 2.5 为枢纽；P 分界 0.30/0.75 用于枢纽）给节点分角色。
# 输出解读：真实系统划分的 Q* 明显高于随机打乱标签的 Q*；模拟中的跨系统节点 9 应有最高的 P。
# 小网络里很少有节点达到 z ≥ 2.5，这是该阈值源自大型网络的缘故，阈值不宜机械套用。
import numpy as np

def simulate_fc(n_mod=4, per=10, T=300, seed=1):
    rng = np.random.default_rng(seed)        # 模拟数据，仅用于演示方法
    labels = np.repeat(np.arange(n_mod), per)
    sys_sig = rng.standard_normal((n_mod, T))
    sys_sig[1] -= 0.6 * sys_sig[0]           # 系统 0 与 1 反相关，产生负权
    load = rng.uniform(0.4, 1.0, size=n_mod * per)
    X = load[:, None] * sys_sig[labels] + rng.standard_normal((n_mod * per, T))
    X[9] += 0.9 * sys_sig[2] + 0.9 * sys_sig[3]   # 节点 9：跨系统连接
    return np.corrcoef(X), labels

def modularity_signed(W, labels):
    """Rubinov & Sporns (2011) 非对称加权模块度 Q*。"""
    Wp, Wn = np.clip(W, 0, None), np.clip(-W, 0, None)
    same = labels[:, None] == labels[None, :]
    def q(A):
        v = A.sum()
        if v == 0:
            return 0.0, 0.0
        k = A.sum(1)
        return ((A - np.outer(k, k) / v) * same).sum() / v, v
    qp, vp = q(Wp)
    qn, vn = q(Wn)
    return qp - vn / (vp + vn) * qn

def roles(W, labels):
    """正权网络的模块内强度 z 分数与参与系数。"""
    A = np.clip(W, 0, None)
    k = A.sum(1)
    M = labels.max() + 1
    kis = np.column_stack([A[:, labels == m].sum(1) for m in range(M)])
    P = 1 - ((kis / k[:, None]) ** 2).sum(1)
    k_in = kis[np.arange(len(labels)), labels]
    z = np.zeros(len(labels))
    for m in range(M):
        idx = labels == m
        z[idx] = (k_in[idx] - k_in[idx].mean()) / (k_in[idx].std() + 1e-12)
    return z, P

R, labels = simulate_fc()
np.fill_diagonal(R, 0)
rng = np.random.default_rng(0)
q_true = modularity_signed(R, labels)
q_perm = [modularity_signed(R, rng.permutation(labels)) for _ in range(200)]
print(f"真实系统划分 Q* = {q_true:.3f}；随机标签 Q* 均值 = {np.mean(q_perm):.3f}，最大 = {np.max(q_perm):.3f}")
print("负权占全部权重绝对值的比例 =", round(np.clip(-R, 0, None).sum() / np.abs(R).sum(), 3))

z, P = roles(R, labels)
print("参与系数最高的 3 个节点：", np.argsort(P)[::-1][:3].tolist(), "P =", np.round(np.sort(P)[::-1][:3], 2).tolist())
print("模块内 z 最高的 3 个节点：", np.argsort(z)[::-1][:3].tolist(), "z =", np.round(np.sort(z)[::-1][:3], 2).tolist())
hub = z >= 2.5
print("z ≥ 2.5 的枢纽数 =", int(hub.sum()), "（40 个节点的小网络中通常为 0）")
for i in np.where(hub)[0]:
    role = "省内枢纽" if P[i] < 0.30 else ("连接枢纽" if P[i] < 0.75 else "交汇枢纽")
    print(f"  节点 {i}: z = {z[i]:.2f}, P = {P[i]:.2f} → {role}")

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
vmax = np.abs(R).max()
im = axes[0].imshow(R, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
for b in np.where(np.diff(labels))[0] + 0.5:
    axes[0].axhline(b, color="k", lw=0.6); axes[0].axvline(b, color="k", lw=0.6)
fig.colorbar(im, ax=axes[0], shrink=0.8, label="r")
axes[0].set(title=f"Signed FC: Q* = {q_true:.2f} (null max {np.max(q_perm):.2f})", xlabel="Region", ylabel="Region")
for m, mk in zip(range(4), "osD^"):
    idx = labels == m
    axes[1].scatter(P[idx], z[idx], marker=mk, s=50, alpha=0.8, label=f"system {m}")
axes[1].annotate("node 9", (P[9], z[9]), xytext=(6, 4), textcoords="offset points", color="crimson")
axes[1].axhline(2.5, color="gray", ls="--", lw=1)
for x0 in (0.30, 0.75):
    axes[1].axvline(x0, color="gray", ls=":", lw=1)
axes[1].set(xlabel="Participation coefficient P", ylabel="Within-module z", title="Functional roles (Guimera-Amaral plane)")
axes[1].legend(fontsize=8)
plt.tight_layout(); plt.show()