目标控制——只控制一部分节点需要多少驱动节点

# 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()