# 依赖：networkx>=3, numpy, scipy, matplotlib
# 数据：有向网络用边列表 [(u, v), ...] 表示有向边 u -> v（节点编号 0..N-1）；系统矩阵约定 A[v, u] != 0
#       表示 u 影响 v。示例网络由 networkx 生成（有向 ER 随机图、BA 网络、小世界网络等）或用静态模型生成有向无标度网络；
#       边权用随机数生成，仅用于演示方法。非线性与布尔网络示例使用小型人工网络。

# %% 示例：用 Kalman 秩判据检验可控性
# 系统 dx/dt = A x + B u，A[i, j] != 0 表示有向边 j -> i；B 的每一列是一个独立输入，作用在一个驱动节点上。
# 可控性矩阵 C = [B, AB, A^2 B, ..., A^{N-1} B]，满秩（rank = N）则可控。
# 检查三种情形：(a) 链 0->1->2->3 只驱动节点 0；(b) 星 0->1, 0->2, 0->3 只驱动节点 0；
# (c) 星形再加驱动节点 2、3。解读：(a) 可控，(b) 不可控（叶子只能收到同一信号的不同倍数），(c) 可控。
import numpy as np

def system_matrix(edges, n, seed=0):
    rng = np.random.default_rng(seed)
    A = np.zeros((n, n))
    for u, v in edges:
        A[v, u] = rng.uniform(0.5, 1.5)        # 随机非零权重
    return A

def input_matrix(drivers, n):
    B = np.zeros((n, len(drivers)))
    for j, d in enumerate(drivers):
        B[d, j] = 1.0
    return B

def ctrb_rank(A, B):
    n = A.shape[0]
    C = np.hstack([np.linalg.matrix_power(A, k) @ B for k in range(n)])
    return np.linalg.matrix_rank(C)

chain = [(0, 1), (1, 2), (2, 3)]
star = [(0, 1), (0, 2), (0, 3)]
cases = []
for name, edges, drv in [("链, 驱动{0}", chain, [0]), ("星, 驱动{0}", star, [0]), ("星, 驱动{0,2,3}", star, [0, 2, 3])]:
    r = ctrb_rank(system_matrix(edges, 4), input_matrix(drv, 4))
    cases.append((edges, drv, r))
    print("%-16s rank = %d / 4  -> %s" % (name, r, "可控" if r == 4 else "不可控"))

import networkx as nx
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(10, 3.4))
titles = ["chain, drive {0}", "star, drive {0}", "star, drive {0,2,3}"]
for ax, (edges, drv, r), title in zip(axes, cases, titles):
    G = nx.DiGraph(edges)
    pos = {0: (0, 0), 1: (1, 0), 2: (2, 0), 3: (3, 0)} if edges is chain else {0: (0, 0), 1: (1, 1), 2: (1, 0), 3: (1, -1)}
    nx.draw_networkx_edges(G, pos, ax=ax, edge_color="0.5", arrows=True, arrowsize=14, node_size=400)
    nx.draw_networkx_nodes(G, pos, ax=ax, node_size=400, node_color=["tab:red" if v in drv else "0.85" for v in G])
    nx.draw_networkx_labels(G, pos, ax=ax, font_size=9)
    ax.set_title("%s\nrank = %d / 4 (%s)" % (title, r, "controllable" if r == 4 else "not controllable"), fontsize=10)
    ax.axis("off"); ax.margins(0.2)
plt.tight_layout(); plt.show()

# %% 示例：用最大匹配找驱动节点，并用 Kalman 判据验证
# 把有向图拆成二分图：每个节点 u 有出端 ('o', u) 和入端 ('i', u)，边 u->v 连接 ('o',u)-('i',v)。
# 用 Hopcroft–Karp 求最大匹配后，入端未被匹配的节点就是驱动节点；若全部匹配，则任选 1 个节点。
# 再给随机有向网络的边赋随机权重，只在这些驱动节点上施加输入，检查 Kalman 秩是否等于 N。
# 解读：秩 = N 说明最大匹配给出的驱动节点集足以控制（几乎所有权重下）；N_D 是最小驱动节点数，集合本身不唯一。
import numpy as np
import networkx as nx

def driver_nodes(G):
    Bg = nx.Graph()
    top = [("o", u) for u in G]
    Bg.add_nodes_from(top)
    Bg.add_nodes_from([("i", v) for v in G])
    Bg.add_edges_from((("o", u), ("i", v)) for u, v in G.edges())
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=top)
    D = [v for v in G if ("i", v) not in M]
    return D if D else [next(iter(G))]

G = nx.gnm_random_graph(12, 18, directed=True, seed=4)
D = driver_nodes(G)
print("驱动节点:", D, " N_D =", len(D))

n = G.number_of_nodes()
rng = np.random.default_rng(0)
A = np.zeros((n, n))
for u, v in G.edges():
    A[v, u] = rng.uniform(0.5, 1.5)
A /= max(1.0, np.abs(np.linalg.eigvals(A)).max())   # 缩放，避免高次幂数值溢出
B = np.zeros((n, len(D)))
for j, d in enumerate(D):
    B[d, j] = 1.0
C = np.hstack([np.linalg.matrix_power(A, k) @ B for k in range(n)])
print("Kalman 秩 =", np.linalg.matrix_rank(C), "/", n)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
pos = nx.kamada_kawai_layout(G.to_undirected())
nx.draw_networkx_edges(G, pos, ax=ax, edge_color="0.5", alpha=0.6, arrows=True, arrowsize=10, node_size=250)
nx.draw_networkx_nodes(G, pos, ax=ax, node_size=250, node_color=["tab:red" if v in D else "0.85" for v in G])
nx.draw_networkx_labels(G, pos, ax=ax, font_size=8)
ax.set_title("Driver nodes from maximum matching (red): N_D = %d, rank = %d / %d"
             % (len(D), np.linalg.matrix_rank(C), n), fontsize=10)
ax.axis("off")
plt.tight_layout(); plt.show()

# %% 示例：驱动节点比例与平均度、度异质性
# 比较 ER 有向网络与无标度有向网络（N=1000，静态模型，度指数 gamma=2.5）在不同平均度下的
# 驱动节点比例 n_D = N_D/N，并检查驱动节点的平均度是否低于全网平均度（驱动节点避开枢纽）。
# 解读：平均度越高 n_D 越低；同等平均度下无标度网络 n_D 更高；驱动节点的平均度低于全网平均度。
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

def driver_nodes(G):
    Bg = nx.Graph()
    top = [("o", u) for u in G]
    Bg.add_nodes_from(top); Bg.add_nodes_from([("i", v) for v in G])
    Bg.add_edges_from((("o", u), ("i", v)) for u, v in G.edges())
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=top)
    D = [v for v in G if ("i", v) not in M]
    return D if D else [next(iter(G))]

def sf_directed(N, k_mean, gamma=2.5, seed=0):
    """静态模型：节点权重 w_i ∝ i^{-1/(gamma-1)}，按权重分别抽取起点和终点。"""
    rng = np.random.default_rng(seed)
    w = np.arange(1, N + 1) ** (-1.0 / (gamma - 1)); w /= w.sum()
    L = int(k_mean * N)
    codes = np.array([], dtype=np.int64)
    while codes.size < L:                       # 批量抽样，去掉自环和重边
        u = rng.choice(N, 4 * L, p=w); v = rng.choice(N, 4 * L, p=w)
        codes = np.concatenate([codes, (u * N + v)[u != v]])
        _, idx = np.unique(codes, return_index=True)
        codes = codes[np.sort(idx)]
    G = nx.DiGraph(); G.add_nodes_from(range(N))
    G.add_edges_from((int(c // N), int(c % N)) for c in codes[:L])
    return G

N = 1000
ks = [1, 2, 3, 4, 6]
nd_er, nd_sf = [], []
for k in ks:
    G_er = nx.gnm_random_graph(N, k * N, directed=True, seed=k)
    G_sf = sf_directed(N, k, seed=k)
    nd_er.append(len(driver_nodes(G_er)) / N)
    nd_sf.append(len(driver_nodes(G_sf)) / N)
    print("<k>=%d  n_D(ER) = %.3f  n_D(SF) = %.3f" % (k, nd_er[-1], nd_sf[-1]))

plt.plot(ks, nd_er, "o-", label="ER"); plt.plot(ks, nd_sf, "s-", label="SF (gamma=2.5)")
plt.xlabel("<k> (out-degree)"); plt.ylabel("n_D"); plt.legend(); plt.title("Driver node fraction")
plt.show()

G = sf_directed(N, 3, seed=9)
D = set(driver_nodes(G))
deg = dict(G.degree())
print("SF 全网平均度 = %.2f，驱动节点平均度 = %.2f" % (np.mean(list(deg.values())), np.mean([deg[v] for v in D])))

# %% 示例：PBH 判据与精确可控性——同一拓扑、不同权重需要的输入数不同
# PBH 判据：(A,B) 可控 <=> 对 A 的每个特征值 lambda，rank[lambda I - A, B] = N。
# Yuan 等（2013）的精确可控性：给定 A 时最少独立输入数 N_D = max_lambda [N - rank(lambda I - A)]（最大几何重数）。
# 比较 (a) 全 1 权重的完全有向图（N_D = N-1）、(b) 同一拓扑的随机权重（N_D = 1，与最大匹配一致）、
# (c) 稀疏 ER 网络：用 N - rank(A) 近似（Yuan 2013 式 (5)），并与最大匹配给出的 N_D 对比。
# 解读：结构可控性给的是“几乎所有权重”下的答案；权重相等这类约束会把系统推进例外集合。
import numpy as np
import networkx as nx

def exact_ND(A, tol=1e-8):
    """最大几何重数：对每个（数值上不同的）特征值计算 N - rank(lambda I - A)。"""
    N = A.shape[0]
    lams = np.linalg.eigvals(A)
    best, seen = 1, []
    for lam in lams:
        if any(abs(lam - s) < 1e-6 for s in seen):
            continue
        seen.append(lam)
        r = np.linalg.matrix_rank(lam * np.eye(N) - A, tol=tol * max(1.0, np.abs(A).max()) * N)
        best = max(best, N - r)
    return best

def matching_ND(G):
    Bg = nx.Graph()
    top = [("o", u) for u in G]
    Bg.add_nodes_from(top); Bg.add_nodes_from([("i", v) for v in G])
    Bg.add_edges_from((("o", u), ("i", v)) for u, v in G.edges())
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=top)
    return max(1, sum(1 for v in G if ("i", v) not in M))

N = 6
A_unit = np.ones((N, N)) - np.eye(N)                      # 完全图，权重全为 1
rng = np.random.default_rng(1)
A_rand = A_unit * rng.uniform(0.5, 1.5, size=(N, N))      # 同一拓扑，独立随机权重
print("完全图 N=%d：全 1 权重 N_D = %d；随机权重 N_D = %d；最大匹配 N_D = %d"
      % (N, exact_ND(A_unit), exact_ND(A_rand), matching_ND(nx.complete_graph(N, create_using=nx.DiGraph))))

# PBH 检查：全 1 权重时，只驱动 1 个节点在 lambda = -1 处秩不足
B1 = np.zeros((N, 1)); B1[0, 0] = 1
print("全 1 权重、驱动 1 个节点：rank[-I-A, B] = %d / %d"
      % (np.linalg.matrix_rank(np.hstack([-np.eye(N) - A_unit, B1])), N))

G = nx.gnm_random_graph(300, 450, directed=True, seed=3)
A = nx.to_numpy_array(G).T                                 # A[v,u] = 1 表示 u -> v
print("稀疏有向 ER（N=300, L=450，单位权重）：N - rank(A) = %d；最大几何重数 = %d；最大匹配 N_D = %d"
      % (300 - np.linalg.matrix_rank(A), exact_ND(A), matching_ND(G)))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ev_u, ev_r = np.linalg.eigvals(A_unit), np.linalg.eigvals(A_rand)
ax1.plot(ev_u.real, ev_u.imag, "o", ms=12, mfc="none", mew=2, label="unit weights (N_D = %d)" % exact_ND(A_unit))
ax1.plot(ev_r.real, ev_r.imag, "x", ms=8, mew=2, label="random weights (N_D = %d)" % exact_ND(A_rand))
ax1.annotate("lambda = -1, multiplicity %d" % (N - 1), (-1, 0), xytext=(0.2, 0.3), arrowprops=dict(arrowstyle="->"))
ax1.set_ylim(-0.6, 0.6)
ax1.set_xlabel("Re lambda"); ax1.set_ylabel("Im lambda"); ax1.legend(fontsize=8)
ax1.set_title("Complete digraph N = %d: degenerate spectrum" % N)
ev = np.linalg.eigvals(A)
ax2.plot(ev.real, ev.imag, ".", c="tab:blue", alpha=0.6)
ax2.set_xlabel("Re lambda"); ax2.set_ylabel("Im lambda")
ax2.set_title("Sparse ER (unit weights): %d of %d eigenvalues near 0" % ((np.abs(ev) < 1e-6).sum(), len(ev)), fontsize=10)
plt.tight_layout(); plt.show()

# %% 示例：控制中心性——单个节点能控制的子空间维数
# Liu、Slotine 与 Barabási（2012）把“只驱动节点 i 时可控子空间的通用维数”定义为节点 i 的控制中心性。
# 可控子空间就是 Krylov 子空间 span{b, Ab, A^2 b, ...}；这里给非零边赋随机权重（对应“几乎所有权重”），
# 用逐步正交化求其维数，避免直接计算 A^{N-1} 带来的数值病态。
# 对无环有向图（DAG），控制中心性等于节点在层次结构中的层号（Hosoe 定理的特例），本例同时打印两者对照。
# 解读：控制中心性高的节点未必是度最高的节点；它反映的是从该节点出发能“展开”多长的不相交路径与环。
import numpy as np
import networkx as nx

def krylov_dim(A, b, tol=1e-9):
    Q = []
    v = b / np.linalg.norm(b)
    for _ in range(A.shape[0]):
        for q in Q:                              # 两次 Gram–Schmidt，保证正交性
            v = v - (q @ v) * q
        for q in Q:
            v = v - (q @ v) * q
        nv = np.linalg.norm(v)
        if nv < tol:
            break
        v = v / nv
        Q.append(v)
        v = A @ v
    return len(Q)

def control_centrality(G, seed=0):
    nodes = list(G)
    idx = {u: k for k, u in enumerate(nodes)}
    rng = np.random.default_rng(seed)
    A = np.zeros((len(nodes), len(nodes)))
    for u, v in G.edges():
        A[idx[v], idx[u]] = rng.uniform(0.5, 1.5)
    out = {}
    for u in nodes:
        b = np.zeros(len(nodes)); b[idx[u]] = 1.0
        out[u] = krylov_dim(A, b)
    return out

# 一个小的有向无环网络（层次结构）
G = nx.DiGraph([(0, 1), (0, 2), (1, 3), (2, 3), (3, 4), (4, 5), (1, 6), (6, 5), (7, 6)])
cc = control_centrality(G)
# 层号：从节点出发的最长有向路径所含的节点数（DAG 情形，对应 Liu 2012 的层次标记）
longest = {u: 1 + max([0] + [len(p) - 1 for w in nx.descendants(G, u) for p in nx.all_simple_paths(G, u, w)]) for u in G}
for u in sorted(G):
    print("节点 %d：出度 %d，控制中心性 %d，最长下游路径节点数 %d" % (u, G.out_degree(u), cc[u], longest[u]))

# 含环网络：随机有向网络中控制中心性的分布
H = nx.gnm_random_graph(60, 90, directed=True, seed=2)
ccH = control_centrality(H)
vals = np.array(list(ccH.values()))
outdeg = np.array([H.out_degree(u) for u in H])
print("随机有向网络 N=60：控制中心性 最大 %d，中位数 %.0f；与出度的 Spearman 相关 = %.2f"
      % (vals.max(), np.median(vals),
         np.corrcoef(np.argsort(np.argsort(vals)), np.argsort(np.argsort(outdeg)))[0, 1]))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
layers = {}
for u in sorted(G):
    layers.setdefault(longest[u], []).append(u)
pos = {u: (-L, j - (len(us) - 1) / 2) for L, us in layers.items() for j, u in enumerate(us)}
nx.draw_networkx_edges(G, pos, ax=ax1, edge_color="0.5", arrows=True, arrowsize=12, node_size=350)
nodes_ = nx.draw_networkx_nodes(G, pos, ax=ax1, node_size=350, node_color=[cc[u] for u in G], cmap="YlGnBu", vmin=0, vmax=7, edgecolors="0.3")
nx.draw_networkx_labels(G, pos, ax=ax1, font_size=8)
plt.colorbar(nodes_, ax=ax1, label="control centrality"); ax1.axis("off")
ax1.set_title("DAG: control centrality = layer index")
jit = np.random.default_rng(0).uniform(-0.15, 0.15, (2, len(vals)))
ax2.scatter(outdeg + jit[0], vals + jit[1], s=30, alpha=0.6)
ax2.set_xlabel("out-degree"); ax2.set_ylabel("control centrality (dim of controllable subspace)")
ax2.set_title("Random digraph N = 60: not determined by degree")
plt.tight_layout(); plt.show()

# %% 示例：目标控制——只控制一部分节点需要多少驱动节点
# Gao 等（2014）的贪心算法：右侧放当前目标集合，左侧放指向它们的上游节点，求二分图最大匹配；
# 未被匹配的目标成为驱动节点，被匹配到的上游节点成为下一轮的“目标”，迭代直到集合为空或重复出现。
# 最后用随机权重的目标可控性秩判据 rank(C_T [B, AB, ..., A^{N-1}B]) = |T| 检验结果。
# 解读：控制一小部分目标所需的驱动节点通常远少于全控制的 N_D；结果是近似（上界），不保证最小。
import numpy as np
import networkx as nx

def matching(left_of, targets, n):
    # 左侧节点 u 编号为 u，右侧目标 t 编号为 n + t（整数编号使结果可复现）
    targets = sorted(targets)
    L = sorted(set().union(*[left_of[t] for t in targets])) if targets else []
    if not L:
        return targets, set()
    Bg = nx.Graph()
    Bg.add_nodes_from(L); Bg.add_nodes_from(n + t for t in targets)
    Bg.add_edges_from((u, n + t) for t in targets for u in sorted(left_of[t]))
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=L)
    unmatched = [t for t in targets if n + t not in M]
    matched_up = {M[n + t] for t in targets if n + t in M}
    return unmatched, matched_up

def target_drivers(G, T):
    """G 的节点须为 0..n-1 的整数。"""
    n = G.number_of_nodes()
    left_of = {v: set(G.predecessors(v)) for v in G}
    drivers, current, seen = set(), set(T), []
    while current and current not in seen:
        seen.append(set(current))
        unmatched, current = matching(left_of, current, n)
        drivers |= set(unmatched)
    if not drivers:                              # 目标全部落在环上：至少需要 1 个输入
        drivers = {min(T)}
    return drivers

def target_rank(G, drivers, T, seed=0):
    nodes = list(G); idx = {u: k for k, u in enumerate(nodes)}; n = len(nodes)
    rng = np.random.default_rng(seed)
    A = np.zeros((n, n))
    for u, v in G.edges():
        A[idx[v], idx[u]] = rng.uniform(0.5, 1.5)
    A /= max(1.0, np.abs(np.linalg.eigvals(A)).max())
    B = np.zeros((n, len(drivers)))
    for j, d in enumerate(drivers):
        B[idx[d], j] = 1.0
    blocks, X = [], B.copy()
    for _ in range(n):
        blocks.append(X); X = A @ X
    C = np.hstack(blocks)[[idx[t] for t in T], :]
    s = np.linalg.svd(C, compute_uv=False)
    return int((s > 1e-10 * s[0]).sum())

G = nx.gnm_random_graph(40, 60, directed=True, seed=7)
full = target_drivers(G, set(G))
print("全控制：驱动节点 %d 个（应等于最大匹配的 N_D）" % len(full))
rng = np.random.default_rng(0)
bars = []
for frac in [0.1, 0.25, 0.5]:
    T = set(rng.choice(40, int(frac * 40), replace=False).tolist())
    D = target_drivers(G, T)
    bars.append((len(T), len(D)))
    print("目标比例 %.2f：|T| = %d，驱动节点 %d 个，目标可控性秩 = %d / %d"
          % (frac, len(T), len(D), target_rank(G, D, sorted(T)), len(T)))

import matplotlib.pyplot as plt
bars.append((40, len(full)))
fig, ax = plt.subplots(figsize=(6.4, 4.2))
x = np.arange(len(bars)); labels = ["10%", "25%", "50%", "100% (full)"]
ax.bar(x - 0.18, [b[0] for b in bars], 0.36, color="0.7", label="targets |T|")
ax.bar(x + 0.18, [b[1] for b in bars], 0.36, color="tab:blue", label="driver nodes")
ax.set_xticks(x, labels); ax.set_xlabel("target fraction"); ax.set_ylabel("number of nodes"); ax.legend()
ax.set_title("Target control needs far fewer drivers than full control")
plt.tight_layout(); plt.show()

# %% 示例：控制能量、Gramian 条件数与“理论可控、数值不可控”
# 无向网络的稳定线性系统 A = -(L + 0.1 I)（L 为拉普拉斯矩阵），随机选 N_D 个驱动节点，
# 用 Lyapunov 方程 A W + W A^T + B B^T = 0 求无限时域可控性 Gramian W。
# 最小能量的上下界分别是 1/lambda_min(W) 与 1/lambda_max(W)（Yan 等 2012）；W 的条件数决定数值上能否求逆。
# 解读：驱动节点减少时，最难方向的能量与条件数近似指数增长（Yan 等 2015 的 E_max ~ e^{N/N_D}）；
# 只驱动 1–2 个节点时 W 的最小特征值已低于双精度可分辨的量级（约 1e-16 倍最大特征值），算出负值：
# 此时数值上已无法判断是否可控，更谈不上实际操控。
import numpy as np
import networkx as nx
from scipy.linalg import solve_continuous_lyapunov

N = 40
G = nx.connected_watts_strogatz_graph(N, 4, 0.2, seed=1)
A = -(nx.laplacian_matrix(G).toarray().astype(float) + 0.1 * np.eye(N))
rng = np.random.default_rng(0)
spectra, cond_med = {}, []
for ND in [40, 20, 10, 5, 2, 1]:
    lmins, conds = [], []
    for rep in range(10):
        drv = rng.choice(N, ND, replace=False)
        B = np.zeros((N, ND)); B[drv, np.arange(ND)] = 1.0
        W = solve_continuous_lyapunov(A, -B @ B.T)
        ev = np.linalg.eigvalsh((W + W.T) / 2)
        lmins.append(ev[0]); conds.append(ev[-1] / max(abs(ev[0]), 1e-300))
    lm = np.median(lmins)
    emax = "%8.1e" % (1 / lm) if lm > 0 else "无法求逆"
    print("N_D = %2d：lambda_min(W) 中位数 = %9.2e（E_max ~ %s），条件数中位数 = %8.1e，算出非正特征值 %d/10 次"
          % (ND, lm, emax, np.median(conds), sum(l <= 0 for l in lmins)))
    spectra[ND] = ev; cond_med.append(np.median(conds))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
for ND, ls in zip(spectra, ["-", "--", "-.", ":", "-", "--"]):
    ev = np.sort(np.abs(spectra[ND]))[::-1]
    ax1.semilogy(np.arange(1, N + 1), ev, ls=ls, marker=".", label="N_D = %d" % ND)
ax1.axhline(1e-16 * max(np.abs(spectra[1])), c="tab:red", lw=0.8, ls=":")
ax1.set_xlabel("index"); ax1.set_ylabel("|eigenvalue of W| (sorted)"); ax1.legend(fontsize=8, ncol=2)
ax1.set_title("Gramian spectrum (one draw per N_D)")
ax2.semilogy(list(spectra), cond_med, "ko-")
ax2.set_xscale("log"); ax2.set_xlabel("number of driver nodes N_D"); ax2.set_ylabel("median condition number of W")
ax2.set_title("Fewer drivers: numerically uncontrollable")
plt.tight_layout(); plt.show()

# %% 示例：平均可控性与模态可控性（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()

# %% 示例：可观测性——根强连通分支与最少观测节点
# Liu、Slotine 与 Barabási（2013）：画推断图（若 x_j 出现在 x_i 的方程中，则连 x_i -> x_j），
# 分解为强连通分支（SCC）；没有入边的“根 SCC”每个至少要放一个观测节点——这是任何（非线性）系统可观测的必要条件。
# 线性系统可借助对偶性：(A, C) 可观测 <=> (A^T, C^T) 可控，所以最少观测节点数 = 反向网络的最少驱动节点数。
# 解读：根 SCC 个数是下界；线性系统的充分观测集合通常更大（还要消除“扩张”）。
import numpy as np
import networkx as nx

def matching_drivers(G):
    Bg = nx.Graph()
    top = [("o", u) for u in G]
    Bg.add_nodes_from(top); Bg.add_nodes_from([("i", v) for v in G])
    Bg.add_edges_from((("o", u), ("i", v)) for u, v in G.edges())
    M = nx.bipartite.hopcroft_karp_matching(Bg, top_nodes=top)
    D = [v for v in G if ("i", v) not in M]
    return D if D else [next(iter(G))]

# 系统有向图 S：u -> v 表示 u 影响 v（u 出现在 v 的方程中）；推断图 = S 的反向图
S = nx.gnm_random_graph(200, 300, directed=True, seed=5)
inference = S.reverse(copy=True)
cond = nx.condensation(inference)
root_sccs = [c for c in cond if cond.in_degree(c) == 0]
sensors_lower = [min(cond.nodes[c]["members"]) for c in root_sccs]     # 每个根 SCC 任取一个节点
sensors_linear = matching_drivers(inference)                            # 线性情形：反向图上的驱动节点
print("根 SCC 个数（任意动力学的必要观测节点数）= %d" % len(root_sccs))
print("线性动力学的最少观测节点数（对偶：反向图最大匹配）= %d" % len(sensors_linear))

# 用 Kalman 可观测性矩阵在随机权重下检验线性观测集合
n = S.number_of_nodes()
rng = np.random.default_rng(0)
A = np.zeros((n, n))
for u, v in S.edges():
    A[v, u] = rng.uniform(0.5, 1.5)
A /= max(1.0, np.abs(np.linalg.eigvals(A)).max())
C = np.zeros((len(sensors_linear), n)); C[np.arange(len(sensors_linear)), sensors_linear] = 1.0
blocks, X = [], C.copy()
for _ in range(n):
    blocks.append(X); X = X @ A
O = np.vstack(blocks)
s = np.linalg.svd(O, compute_uv=False)
print("线性观测集合的可观测性秩 = %d / %d" % (int((s > 1e-10 * s[0]).sum()), n))
C2 = np.zeros((len(sensors_lower), n)); C2[np.arange(len(sensors_lower)), sensors_lower] = 1.0
blocks, X = [], C2.copy()
for _ in range(n):
    blocks.append(X); X = X @ A
s2 = np.linalg.svd(np.vstack(blocks), compute_uv=False)
print("只放根 SCC 观测节点时的线性可观测性秩 = %d / %d（不足说明线性系统需要更多观测）"
      % (int((s2 > 1e-10 * s2[0]).sum()), n))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.semilogy(np.arange(1, len(s) + 1), s / s[0], "-", lw=2, label="linear sensor set (%d sensors)" % len(sensors_linear))
ax.semilogy(np.arange(1, len(s2) + 1), np.maximum(s2 / s2[0], 1e-20), "--", lw=2,
            label="root-SCC sensors only (%d)" % len(sensors_lower))
ax.axhline(1e-10, c="tab:red", lw=0.8, ls=":", label="rank tolerance")
ax.set_ylim(1e-20, 2); ax.set_xlabel("singular value index"); ax.set_ylabel("normalized singular value")
ax.legend(fontsize=8); ax.set_title("Observability matrix spectrum: root SCCs are not enough")
plt.tight_layout(); plt.show()

# %% 示例：反馈顶点集控制——锁定少数节点把非线性网络推向指定吸引子
# 四基因网络：X 与 Y 互相抑制（双稳态开关），X 激活 Z，Y 抑制 W，Z 与 W 互相激活形成第二个环。
# 反馈顶点集（FVS）是删去后网络无环的节点集合；Fiedler 与 Mochizuki 等（2013）证明：在温和条件下，
# 把 FVS 节点的轨迹锁定为目标吸引子上的值，其余节点会收敛到该吸引子（Zañudo 等 2017 在此基础上加入源节点）。
# 本例先求两个稳态，再从稳态 1 出发锁定最小 FVS，检查其余节点是否收敛到稳态 2。
# 解读：锁定 FVS（2 个节点）即可在不依赖具体参数的意义上完成状态切换；只锁定 Z（不是 FVS）失败。
# 在这组具体参数下只锁定 X 也成功——FVS 是“对所有满足条件的非线性函数都有效”的充分集合，具体模型可能需要更少。
import itertools
import numpy as np
import networkx as nx
from scipy.integrate import solve_ivp

edges = {("X", "Y"): -1, ("Y", "X"): -1, ("X", "Z"): +1, ("Y", "W"): -1, ("Z", "W"): +1, ("W", "Z"): +1}
nodes = ["X", "Y", "Z", "W"]
G = nx.DiGraph(list(edges))
hill = lambda x, s: x**4 / (0.5**4 + x**4) if s > 0 else 0.5**4 / (0.5**4 + x**4)

def rhs(t, x, clamp):
    d = dict(zip(nodes, x)); out = []
    for v in nodes:
        if v in clamp:
            out.append(0.0); continue
        regs = [hill(d[u], s) for (u, w), s in edges.items() if w == v]
        out.append(np.mean(regs) - d[v])
    return out

def run(x0, clamp=()):
    return solve_ivp(rhs, (0, 60), x0, args=(set(clamp),), rtol=1e-8).y[:, -1]

s1 = run([1.0, 0.0, 1.0, 1.0]); s2 = run([0.0, 1.0, 0.0, 0.0])
print("稳态 1:", np.round(s1, 3), " 稳态 2:", np.round(s2, 3))

def is_fvs(S):
    H = G.copy(); H.remove_nodes_from(S)
    return nx.is_directed_acyclic_graph(H)
fvs = next(set(S) for k in range(1, 5) for S in itertools.combinations(nodes, k) if is_fvs(S))
print("最小反馈顶点集:", sorted(fvs))

for S in [sorted(fvs), ["Z"], ["X"]]:
    x0 = s1.copy()
    for v in S:
        x0[nodes.index(v)] = s2[nodes.index(v)]      # 把锁定节点设为目标吸引子上的值并保持不变
    xf = run(x0, clamp=S)
    print("锁定 %-10s -> 终态 %s，与稳态 2 的距离 %.3f" % (S, np.round(xf, 3), np.linalg.norm(xf - s2)))

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for ax, S in zip(axes, [sorted(fvs), ["Z"]]):
    x0 = s1.copy()
    for v in S:
        x0[nodes.index(v)] = s2[nodes.index(v)]
    sol = solve_ivp(rhs, (0, 30), x0, args=(set(S),), t_eval=np.linspace(0, 30, 300))
    for i, (v, ls) in enumerate(zip(nodes, ["-", "--", "-.", ":"])):
        l, = ax.plot(sol.t, sol.y[i], ls=ls, lw=2, label=v + (" (clamped)" if v in S else ""))
        ax.axhline(s2[i], c=l.get_color(), lw=0.6, alpha=0.5)
    ax.set_xlabel("time"); ax.set_title("Clamp %s: %s" % ("+".join(S), "reaches state 2" if np.linalg.norm(sol.y[:, -1] - s2) < 0.05 else "fails"))
    ax.legend(fontsize=8)
axes[0].set_ylabel("expression level (thin lines = state 2)")
plt.tight_layout(); plt.show()

# %% 示例：布尔网络的吸引子与最小钉扎控制集（穷举）
# 6 个节点的同步更新布尔网络，先穷举 2^6 个状态找出全部吸引子，再穷举节点子集：
# 把子集中的节点固定为目标吸引子上的取值（钉扎），检查是否从所有初始状态都到达目标吸引子。
# 布尔网络的控制一般是 NP 难问题（Akutsu 等 2007），穷举只适用于很小的网络；稳定基序等方法利用结构降低搜索量。
# 解读：最小钉扎集通常比“全部节点”小得多，并与网络中的反馈环位置有关。
import itertools

nodes = ["A", "B", "C", "D", "E", "F"]
rules = {
    "A": lambda s: not s["B"],               # A、B 互相抑制：双稳态开关
    "B": lambda s: not s["A"],
    "C": lambda s: s["A"] or s["D"],         # C、D 互相激活：正反馈记忆环
    "D": lambda s: s["C"] and not s["F"],
    "E": lambda s: s["C"],
    "F": lambda s: s["B"] and not s["E"],
}

def step(state, pin=None):
    s = dict(zip(nodes, state))
    new = tuple(int(rules[v](s)) for v in nodes)
    if pin:
        new = tuple(pin.get(v, x) for v, x in zip(nodes, new))
    return new

def attractor_from(state, pin=None):
    seen = {}
    while state not in seen:
        seen[state] = len(seen); state = step(state, pin)
    start = seen[state]
    cyc = [s for s, k in seen.items() if k >= start]
    return frozenset(cyc)

all_states = list(itertools.product([0, 1], repeat=6))
attractors = {attractor_from(s) for s in all_states}
fixed = sorted([sorted(a)[0] for a in attractors if len(a) == 1])
print("吸引子个数 %d，其中不动点：%s" % (len(attractors), fixed))

target = max(fixed, key=lambda s: s[nodes.index("C")])     # 选一个 C=1 的不动点作为目标
print("目标不动点（A..F）：", target)
for k in range(1, 7):
    sols = []
    for S in itertools.combinations(nodes, k):
        pin = {v: target[nodes.index(v)] for v in S}
        if all(attractor_from(s, pin) == frozenset([target]) for s in all_states):
            sols.append(S)
    if sols:
        print("最小钉扎集大小 %d，共 %d 种：%s" % (k, len(sols), sols))
        break

# 作图：64 个状态的状态转移图，按所属吸引子（吸引域）着色
import networkx as nx
import matplotlib.pyplot as plt
att_list = sorted(attractors, key=lambda a: sorted(a)[0])
basin = {s: att_list.index(attractor_from(s)) for s in all_states}
STG = nx.DiGraph([(s, step(s)) for s in all_states])
fig, ax = plt.subplots(figsize=(6.4, 4.6))
pos = {}
for b, comp in enumerate(sorted(nx.weakly_connected_components(STG), key=lambda c: basin[next(iter(c))])):
    sub = nx.kamada_kawai_layout(STG.subgraph(comp).to_undirected()) if len(comp) > 1 else {next(iter(comp)): (0, 0)}
    pos.update({s: (xy[0] + 2.6 * (b % 3), xy[1] - 2.6 * (b // 3)) for s, xy in sub.items()})
nx.draw_networkx_edges(STG, pos, ax=ax, edge_color="0.6", alpha=0.5, arrowsize=6, node_size=40)
cmap = plt.get_cmap("tab10")
for b, a in enumerate(att_list):
    members = [s for s in all_states if basin[s] == b]
    lab = "attractor %d (%s, %d states)" % (b, "fixed point" if len(a) == 1 else "cycle %d" % len(a), len(members))
    nx.draw_networkx_nodes(STG, pos, nodelist=members, ax=ax, node_size=[120 if s in a else 40 for s in members],
                           node_color=[cmap(b)], label=lab)
ax.legend(fontsize=7, loc="upper center", bbox_to_anchor=(0.5, -0.02), ncol=2); ax.axis("off")
ax.set_title("State transition graph (large dots = attractor states)", fontsize=10)
plt.tight_layout(); plt.show()

# %% 示例：随机失效与蓄意攻击下的结构鲁棒性
# 对 BA 网络和同样边数的 ER 网络，逐步删除比例 f 的节点（随机删除 / 按初始度从高到低删除），
# 记录巨连通分支占原网络的比例；再用 Molloy–Reed 判据 f_c = 1 - 1/(kappa0 - 1)，kappa0 = <k^2>/<k>，
# 估计随机失效的临界删除比例。
# 解读：BA 网络 kappa0 更大、f_c 更接近 1，对随机失效最鲁棒，但对按度攻击最脆弱。
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

N = 1000
G_ba = nx.barabasi_albert_graph(N, 2, seed=1)
G_er = nx.gnm_random_graph(N, G_ba.number_of_edges(), seed=1)
fs = np.linspace(0, 0.9, 19)

def giant_curve(G, order):
    out = []
    for f in fs:
        H = G.copy()
        H.remove_nodes_from(order[: int(f * N)])
        out.append(max((len(c) for c in nx.connected_components(H)), default=0) / N)
    return out

rng = np.random.default_rng(0)
for name, G, style in [("ER", G_er, "-"), ("BA", G_ba, "--")]:
    rand_order = list(rng.permutation(N))
    attack_order = [v for v, _ in sorted(G.degree(), key=lambda x: -x[1])]
    rand_curve = giant_curve(G, rand_order)
    attack_curve = giant_curve(G, attack_order)
    plt.plot(fs, rand_curve, style, label=name + " random")
    plt.plot(fs, attack_curve, style, label=name + " attack")
    k = np.array([d for _, d in G.degree()], dtype=float)
    kappa0 = (k**2).mean() / k.mean()
    fc = 1 - 1 / (kappa0 - 1)
    i20 = int(np.argmin(np.abs(fs - 0.2)))
    print("%s: kappa0 = %.2f, f_c 估计 = %.3f；删除 20%% 后巨分支：随机 %.2f，攻击 %.2f"
          % (name, kappa0, fc, rand_curve[i20], attack_curve[i20]))
plt.xlabel("fraction removed f"); plt.ylabel("giant component / N"); plt.legend(); plt.title("Robustness")
plt.show()