平均可控性与模态可控性（Gu 等 2015 的离散时间定义）及其与度的关系

# 离散时间模型 x(t+1) = A x(t) + B u(t)；按常见做法把对称加权连接矩阵归一化为 A / (1 + sigma_max(A))，保证稳定。
# 平均可控性：只驱动节点 i 时无限时域 Gramian 的迹 Tr(W_i)，用离散 Lyapunov 方程 W = A W A^T + e_i e_i^T 求解；
# 模态可控性：phi_i = sum_j (1 - lambda_j^2) v_ij^2（lambda_j、v_ij 为 A 的特征值与特征向量分量）。
# 再在“同度序列、权重随机打乱”的零模型上重复，检查与加权度的相关是否同样出现。
# 解读：两种指标分别与加权度强正相关、强负相关，且零模型也能复现——这正是 Tu 等（2018）提出的质疑之一。
import numpy as np
import networkx as nx
from scipy.linalg import solve_discrete_lyapunov

def ave_modal(W):
    A = W / (1 + np.linalg.svd(W, compute_uv=False)[0])
    lam, V = np.linalg.eigh(A)
    n = A.shape[0]
    ave = np.empty(n)
    for i in range(n):
        Bi = np.zeros((n, n)); Bi[i, i] = 1.0
        ave[i] = np.trace(solve_discrete_lyapunov(A, Bi))
    modal = ((1 - lam**2)[None, :] * V**2).sum(axis=1)
    return ave, modal

def rank_corr(x, y):
    rx, ry = np.argsort(np.argsort(x)), np.argsort(np.argsort(y))
    return np.corrcoef(rx, ry)[0, 1]

rng = np.random.default_rng(0)
G = nx.powerlaw_cluster_graph(120, 3, 0.3, seed=2)
Wt = np.zeros((120, 120))
for u, v in G.edges():
    Wt[u, v] = Wt[v, u] = rng.lognormal(0, 1)          # 模拟偏态分布的纤维计数权重
strength = Wt.sum(axis=1)
ave, modal = ave_modal(Wt)
print("真实（模拟）网络：Spearman(平均可控性, 加权度) = %.2f；Spearman(模态可控性, 加权度) = %.2f"
      % (rank_corr(ave, strength), rank_corr(modal, strength)))

H = nx.expected_degree_graph([d for _, d in G.degree()], seed=3, selfloops=False)
Wn = np.zeros((120, 120))
for u, v in H.edges():
    Wn[u, v] = Wn[v, u] = rng.lognormal(0, 1)
ave0, modal0 = ave_modal(Wn)
s0 = Wn.sum(axis=1)
print("随机零模型：      Spearman(平均可控性, 加权度) = %.2f；Spearman(模态可控性, 加权度) = %.2f"
      % (rank_corr(ave0, s0), rank_corr(modal0, s0)))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
for ax, y, y0, name in [(ax1, ave, ave0, "average controllability"), (ax2, modal, modal0, "modal controllability")]:
    ax.scatter(strength, y, s=20, alpha=0.7, label="network (rho = %.2f)" % rank_corr(y, strength))
    ax.scatter(s0, y0, s=20, alpha=0.7, marker="^", facecolors="none", edgecolors="tab:orange",
               label="null model (rho = %.2f)" % rank_corr(y0, s0))
    ax.set_xlabel("weighted degree (strength)"); ax.set_ylabel(name); ax.legend(fontsize=8)
ax1.set_title("Tracks strength ..."); ax2.set_title("... and the null model reproduces it")
plt.tight_layout(); plt.show()