共激活模式（CAP）：种子阈值选帧 + k-means

# 思路来自 Liu & Duyn（2013, PNAS）：只取种子区信号超过阈值（此处取最高 15%）的时间帧，
# 对这些帧的全脑空间模式做 k-means，每类求平均得到 CAP。
# 模拟：20 个脑区；种子区（0 号）高活动时，伴随两种不同的空间模式之一（模式 A：1–6 号；模式 B：7–12 号）。
# 输出解读：全部帧平均（相当于种子相关图）把两种模式混在一起；k = 2 的 CAP 能把它们分开。
import numpy as np
from scipy.cluster.vq import kmeans2

rng = np.random.default_rng(4)
N, T = 20, 3000
X = rng.standard_normal((N, T))
events = rng.choice(T, 300, replace=False)
kind = rng.integers(0, 2, events.size)               # 每个事件随机属于模式 A 或 B
for e, k in zip(events, kind):
    X[0, e] += 3.0
    members = np.arange(1, 7) if k == 0 else np.arange(7, 13)
    X[members, e] += 2.5
X = (X - X.mean(1, keepdims=True)) / X.std(1, keepdims=True)

thr = np.quantile(X[0], 0.85)                         # 种子信号最高 15% 的帧
frames = np.where(X[0] > thr)[0]
F = X[:, frames].T                                    # (帧数, N)
_, lab = kmeans2(F, 2, seed=0, minit="++")
seed_map = F.mean(0)
print(f"选中帧数 = {frames.size}（占 {frames.size / T:.0%}）")
print("全部选中帧平均（类似种子相关图）：区 1–6 均值 %.2f，区 7–12 均值 %.2f"
      % (seed_map[1:7].mean(), seed_map[7:13].mean()))
for c in range(2):
    cap = F[lab == c].mean(0)
    print(f"CAP {c}：帧数 {np.sum(lab == c):3d}，区 1–6 均值 {cap[1:7].mean():.2f}，"
          f"区 7–12 均值 {cap[7:13].mean():.2f}")

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(10, 3.6), sharey=True)
maps = [("Mean of all selected frames", seed_map)] + \
       [(f"CAP {c} ({np.sum(lab == c)} frames)", F[lab == c].mean(0)) for c in range(2)]
for ax, (ttl, mp) in zip(axes, maps):
    ax.bar(np.arange(N), mp, color=np.where(mp >= 0, "#CC6677", "#4477AA"))
    ax.axvspan(0.5, 6.5, color="#DDCC77", alpha=0.25)
    ax.axvspan(6.5, 12.5, color="#88CCEE", alpha=0.25)
    ax.set(title=ttl, xlabel="Region (0 = seed)")
axes[0].set_ylabel("Mean z-score")
plt.tight_layout(); plt.show()