传播溯源——Jordan 中心与基于观测点的估计

# SI 传播从未知源头出发，只观测到 (a) 某一时刻的感染快照，或 (b) 少数观测点被感染的时刻。
# (a) Jordan 中心：在感染子图中，到其他感染节点最大距离最小的节点；
# (b) 观测点法：对每个候选源 s，若它是源头，则 t_o - d(s,o) 对所有观测点应近似相等，取方差最小者。
# 观测点的选择影响很大：这里比较随机选取与按介数最高选取的观测点（各占 5%）。
# 解读：报告估计源与真实源的平均跳数误差与命中率。介数最高的节点集中在网络中部、彼此信息冗余，
#       反而不如随机分布的观测点；“选谁当观测点”本身是一个需要优化的问题。
import numpy as np
import networkx as nx

G = nx.connected_watts_strogatz_graph(400, 6, 0.1, seed=3)
dist = dict(nx.all_pairs_shortest_path_length(G))
nodes = list(G)
bet = nx.betweenness_centrality(G)

def si_times(src, rng, beta=0.5):
    t_inf = {src: 0}; t = 0; frontier = {src}
    while len(t_inf) < len(G):
        t += 1; new = set()
        for u in list(t_inf):
            for v in G.neighbors(u):
                if v not in t_inf and rng.random() < beta:
                    new.add(v)
        for v in new:
            t_inf[v] = t
    return t_inf

def jordan(t_inf, T):
    inf = [v for v, t in t_inf.items() if t <= T]
    return min(inf, key=lambda s: max(dist[s][v] for v in inf))

def observer_est(t_inf, obs):
    return min(nodes, key=lambda s: np.var([t_inf[o] - dist[s][o] for o in obs]))

rng = np.random.default_rng(0)
n_obs = 20
obs_rand = list(rng.choice(nodes, n_obs, replace=False))
obs_bet = sorted(nodes, key=lambda v: -bet[v])[:n_obs]
err = {"Jordan 中心(快照)": [], "观测点-随机": [], "观测点-介数": []}
for _ in range(100):
    src = int(rng.integers(400)); t_inf = si_times(src, rng)
    err["Jordan 中心(快照)"].append(dist[src][jordan(t_inf, 6)])
    err["观测点-随机"].append(dist[src][observer_est(t_inf, obs_rand)])
    err["观测点-介数"].append(dist[src][observer_est(t_inf, obs_bet)])
for k_, v in err.items():
    print("%-14s 平均误差 = %.2f 跳，命中率 = %.2f" % (k_, np.mean(v), np.mean(np.array(v) == 0)))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
names = ["Jordan center (snapshot)", "observers: random", "observers: top betweenness"]
mx = max(max(v) for v in err.values())
for (k_, v), name, ls in zip(err.items(), names, ["-", "--", ":"]):
    h = np.bincount(v, minlength=mx + 1) / len(v)
    ax.plot(np.arange(mx + 1), h, ls=ls, marker="o", ms=4, label="%s (mean %.2f)" % (name, np.mean(v)))
ax.set_xlabel("hop distance from true source"); ax.set_ylabel("fraction of trials"); ax.legend(fontsize=8)
ax.set_title("Source localization error (100 SI outbreaks)")
plt.tight_layout(); plt.show()