结构截断与“伪异配”——简单图约束本身就会产生负的 r

# 用配置模型生成幂律度序列（gamma=2.2，N=5000）。允许重边与自环时（多重图），网络在期望上是中性的；
# 把重边合并、自环删除得到简单图后，r 变为负值。
# 同时比较结构截断 k_s ~ sqrt(<k> N) 与最大度 k_max：k_max 远大于 k_s 时，枢纽之间“应有”的多条边
# 只能保留一条，枢纽被迫连向低度节点。
# 解读：重尾度分布的简单图出现负 r，不能直接当作“生物学上的枢纽回避”。
import numpy as np
import networkx as nx

rng = np.random.default_rng(7)
N, gamma, kmin = 5000, 2.2, 2
# 从离散幂律 p_k ~ k^(-gamma)（k >= kmin）抽样，用连续分布逆变换后取整
u = rng.random(N)
deg = np.floor(kmin * (1 - u) ** (-1 / (gamma - 1))).astype(int)
deg = np.minimum(deg, N - 1)
if deg.sum() % 2:
    deg[0] += 1

Gm = nx.configuration_model(deg.tolist(), seed=11)          # MultiGraph，含重边与自环
r_multi = nx.degree_assortativity_coefficient(Gm)
n_multi = Gm.number_of_edges() - nx.Graph(Gm).number_of_edges()
n_loops = nx.number_of_selfloops(Gm)

Gs = nx.Graph(Gm)                                           # 合并重边
Gs.remove_edges_from(nx.selfloop_edges(Gs))                 # 删除自环
r_simple = nx.degree_assortativity_coefficient(Gs)

k_mean = deg.mean()
k_s = np.sqrt(k_mean * N)
print(f"N = {N}, <k> = {k_mean:.2f}, 最大度 k_max = {deg.max()}, 结构截断 k_s = sqrt(<k>N) ≈ {k_s:.0f}")
print(f"度超过 k_s 的节点数: {(deg > k_s).sum()}")
print(f"多重图（保留重边/自环）: r = {r_multi:+.3f}；重边 {n_multi} 条，自环 {n_loops} 个")
print(f"简单图（合并重边、删自环）: r = {r_simple:+.3f}")

# 对照：把度截断在 k_s 以下再生成，伪异配基本消失
deg_cut = np.minimum(deg, int(k_s))
if deg_cut.sum() % 2:
    deg_cut[0] += 1
Gc = nx.Graph(nx.configuration_model(deg_cut.tolist(), seed=12))
Gc.remove_edges_from(nx.selfloop_edges(Gc))
print(f"度序列截断在 k_s 以下后的简单图: r = {nx.degree_assortativity_coefficient(Gc):+.3f}")

# 对照：泊松度分布（ER 图）没有结构截断问题
Ge = nx.gnm_random_graph(N, int(k_mean * N / 2), seed=13)
print(f"同样平均度的 ER 随机图: r = {nx.degree_assortativity_coefficient(Ge):+.3f}, 最大度 = {max(d for _, d in Ge.degree())}")

# 作图：度分布的互补累积分布与结构截断 k_s；四种构造下的 r
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ks = np.sort(deg)
ax[0].loglog(ks, 1 - np.arange(len(ks)) / len(ks), ".", ms=3, color="C0")
ax[0].axvline(k_s, color="C3", ls="--", label=f"structural cutoff k_s = {k_s:.0f}")
ax[0].set_xlabel("degree k"); ax[0].set_ylabel("P(K >= k)"); ax[0].legend()
ax[0].set_title(f"Power-law degrees (gamma = {gamma}), k_max = {deg.max()}")
r_vals = [r_multi, r_simple, nx.degree_assortativity_coefficient(Gc), nx.degree_assortativity_coefficient(Ge)]
ax[1].bar(["multigraph", "simple graph", "simple, k<=k_s", "ER"], r_vals, color=["0.6", "C3", "C0", "0.6"])
ax[1].axhline(0, color="k", lw=0.8); ax[1].set_ylabel("degree assortativity r")
ax[1].set_title("Removing multi-edges alone makes r negative")
plt.tight_layout(); plt.show()