""" ============================================================================= CS230 · 第 18 讲工具模块：二维"圆环上的 8 个高斯"世界（GAN 和扩散模型共用） ============================================================================= 对应讲义：lectures/lecture18_adversarial_gan_diffusion.md §5–§8（代码 02 的 GAN、代码 03 的 DDPM 都用它） ----------------------------------------------------------------------------- 【为什么用二维的玩具数据】 ----------------------------------------------------------------------------- 生成模型要学的是一个**分布** p(x)，而不是一个 x → y 的函数。真实图像的分布有几万维，我们既画不出来， 也说不清楚"模型漏掉了哪些东西"。二维平面上的分布就不一样了： · 真相是我们自己写下的公式，**完全已知**； · 能直接在终端里画出"样本落在哪里"（字符画密度图）； · 能精确地数出模型覆盖了几个"模式"（mode）、有多少样本是"好样本"。 "圆环上均匀摆放的 8 个高斯"是 GAN 文献里常用的玩具数据（例如 Metz et al. 2017 的 Unrolled GAN）， 专门用来观察 **mode collapse（模式坍缩）**：一个只会生成其中一两团的生成器，一眼就能看出来。 ----------------------------------------------------------------------------- 【这个世界长什么样】 ----------------------------------------------------------------------------- K = 8 个模式，中心均匀地摆在半径 R = 2 的圆上：c_k = (2 cos(2πk/8), 2 sin(2πk/8))，k = 0..7 每个样本：先等概率地选一个模式 k，再从 N(c_k, σ² I) 里抽一个点，σ = 0.05 所以每一"团"都很小（半径约 0.15），团与团之间相距约 1.53 —— 中间是大片的"空白"。 （参数 R、σ 是我在写实验之前定下的：σ 小一点，"样本落在两团之间"这种错误才容易看出来。之后没有为了结果改过。） ----------------------------------------------------------------------------- 【怎样给一个生成器打分】（讲义 §6.5 讨论为什么这件事很难） ----------------------------------------------------------------------------- 因为真相已知，我们可以用几个"知道答案"的指标（真实图像上做不到，那时只能用 FID 之类的近似指标）： · 高质量样本比例 hq_frac：离最近的模式中心不超过 3σ 的样本占多少。 真实数据本身的这个比例是 1 − e^{−4.5} ≈ 0.989（二维高斯落在 3σ 圆内的概率）。 · 覆盖的模式数 covered：一个模式至少分到 MIN_SHARE = 1% 的样本（而且是高质量样本），才算"被覆盖"。 （公平的份额是 1/8 = 12.5%；1% 的门槛只是用来区分"有"和"几乎没有"。） · 最少的模式占公平份额的比例 min_share：8 个模式里分到样本最少的那个，相对 12.5% 有多少。1 表示完全均衡。 · 二维的 Fréchet 距离（FID 的公式，但直接用在 (x, y) 坐标上）：只比较均值和协方差 —— 讲义 §6.5 用它说明 "只看前两阶统计量的指标会被骗"。 ----------------------------------------------------------------------------- 【形状约定】PyTorch 的习惯：每一行一个样本。样本数组的形状一律是 (n, 2)。 【注释标签】【思路】【参数】【返回】【边界】【DL 联系】【numpy 小知识】【Python 小知识】 直接运行本文件做自检（真实数据的各项指标、几种"假生成器"的指标、一张字符画；约 1 秒）： python3 ring_world.py """ import math import os import sys import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) for _d in ("lecture01", "lecture02"): sys.path.insert(0, os.path.join(HERE, "..", _d)) from console import pad, print_header # noqa: E402,F401 （第 1 讲：对齐表格、分隔标题；脚本从这里导入） from text_plot import SERIES, figure_path, get_pyplot, style_axes # noqa: E402,F401 （第 2 讲：画图设置） K_MODES = 8 # 模式个数 RADIUS = 2.0 # 圆环半径 SIGMA = 0.05 # 每个模式的标准差 HQ_RADIUS = 3 * SIGMA # "高质量样本"的门槛：离最近中心 ≤ 3σ MIN_SHARE = 0.01 # 一个模式至少分到 1% 的样本，才算被覆盖 # 8 个中心，形状 (8, 2)。列表推导式：对 k = 0..7 各算一个 (x, y)。 CENTERS = np.array([[RADIUS * math.cos(2 * math.pi * k / K_MODES), RADIUS * math.sin(2 * math.pi * k / K_MODES)] for k in range(K_MODES)]) # ============================================================================= # 1. 数据 # ============================================================================= def sample_ring(n, rng): """ 从"圆环上的 8 个高斯"里抽 n 个样本。 【思路】两步抽样（ancestral sampling，混合分布的标准抽法）：先抽"属于哪个模式" k ~ 均匀{0..7}， 再抽 x ~ N(c_k, σ² I) = c_k + σ · (标准正态)。 【参数】n：样本数；rng：np.random.Generator（每个演示用自己的 rng，保证可复现） 【返回】(n, 2) 的 float32 数组 【numpy 小知识】CENTERS[k]：k 是长度 n 的整数数组时，"花式索引"一次取出 n 行，形状 (n, 2)。 """ k = rng.integers(0, K_MODES, size=n) return (CENTERS[k] + SIGMA * rng.standard_normal((n, 2))).astype(np.float32) # ============================================================================= # 2. 指标：模式覆盖与样本质量 # ============================================================================= def mode_report(x): """ 给一批生成的样本打分：覆盖了几个模式、高质量样本占多少、各模式分到多少。 【思路】对每个样本算它到 8 个中心的距离，取最近的那个中心作为"它属于的模式"； 距离 ≤ 3σ 的叫高质量样本；只数高质量样本来判断一个模式有没有被覆盖。 【参数】x：(n, 2) 数组（numpy 或能转成 numpy 的东西） 【返回】dict： "covered" 覆盖的模式数（0..8） "hq_frac" 高质量样本比例（0..1） "counts" 长度 8 的整数数组：每个模式分到的高质量样本数 "min_share" 分到最少的那个模式 / 公平份额（n / 8 的高质量样本）。0 表示至少丢了一个模式 "mean_dist" 样本到最近中心的平均距离（衡量"离数据有多远"，不需要门槛） 【边界】x 里有 nan / inf（训练发散）时，距离也是 nan：这里把这种样本当成"离所有中心无穷远"， 于是它们不算高质量，平均距离返回 inf。 【numpy 小知识】x[:, None, :] - CENTERS[None, :, :]：(n, 1, 2) − (1, 8, 2) → 广播成 (n, 8, 2)， 一次算出每个样本到每个中心的差；np.linalg.norm(..., axis=2) 沿最后一维求长度 → (n, 8)。 np.bincount(a, minlength=8)：数 a 里 0, 1, ..., 7 各出现几次。 """ x = np.asarray(x, dtype=np.float64) finite = np.all(np.isfinite(x), axis=1) d = np.full((len(x), K_MODES), np.inf) d[finite] = np.linalg.norm(x[finite][:, None, :] - CENTERS[None, :, :], axis=2) nearest = d.argmin(axis=1) dmin = d.min(axis=1) hq = dmin <= HQ_RADIUS counts = np.bincount(nearest[hq], minlength=K_MODES) fair = len(x) / K_MODES return { "covered": int(np.sum(counts >= MIN_SHARE * len(x))), "hq_frac": float(hq.mean()), "counts": counts, "min_share": float(counts.min() / fair), "mean_dist": float(dmin.mean()) if np.all(finite) else float("inf"), } def frechet_2d(a, b): """ 两组二维样本之间的 Fréchet 距离：把每组样本看成一个高斯（只用均值和协方差），再算两个高斯之间的距离。 【公式】FD = ‖μ_a − μ_b‖² + Tr(Σ_a + Σ_b − 2 (Σ_a Σ_b)^{1/2}) 这正是 FID（Fréchet Inception Distance，Heusel et al. 2017）的公式；FID 把它用在 Inception 网络第 2048 维特征上，这里直接用在 (x, y) 坐标上。 【思路】Tr((Σ_a Σ_b)^{1/2}) = Σ_i √λ_i，λ_i 是矩阵 Σ_a Σ_b 的特征值（两个半正定矩阵的乘积，特征值是非负实数）。 2×2 的矩阵不需要 scipy 的 sqrtm，求特征值就够了。 【参数】a、b：(n, 2) 和 (m, 2) 的样本 【返回】float，≥ 0；两组样本来自同一个分布时接近 0（有限样本会有一点点正的偏差） 【边界】它**只看前两阶统计量**：两个均值、协方差相同但形状完全不同的分布，FD = 0（讲义 §6.5、代码 02 演示 3）。 【numpy 小知识】np.cov(a, rowvar=False)：每一行是一个样本时要写 rowvar=False，得到 (2, 2) 的协方差矩阵。 """ a = np.asarray(a, dtype=np.float64) b = np.asarray(b, dtype=np.float64) mu_a, mu_b = a.mean(axis=0), b.mean(axis=0) cov_a, cov_b = np.cov(a, rowvar=False), np.cov(b, rowvar=False) eig = np.linalg.eigvals(cov_a @ cov_b) tr_sqrt = float(np.sum(np.sqrt(np.clip(eig.real, 0.0, None)))) # 数值误差可能给出极小的负数或虚部，裁掉 return float(np.sum((mu_a - mu_b) ** 2) + np.trace(cov_a) + np.trace(cov_b) - 2.0 * tr_sqrt) # ============================================================================= # 3. 字符画：二维密度图 # ============================================================================= def ascii_density(x, lim=3.0, width=37, height=19, chars=" .:-=+*#%@"): """ 把一批二维样本画成字符密度图：每个格子数一数落进来的样本，越多字符越"浓"。 【思路】把 [−lim, lim]² 切成 height 行 × width 列的格子（终端里一个字符大约是 1:2 的竖长方形， 所以列数取行数的两倍左右，画出来的圆才不会被压扁）。计数用平方根缩放， 让"样本少的格子"也能看见（否则最密的那一团会把其他格子都压成空白）。 落在画框外的样本不画，但会在最后一行报告有多少个。 【参数】x：(n, 2)；lim：画框的半边长；width、height：格子数；chars：从空到满的字符表 【返回】list[str]，height + 1 行（最后一行是说明） 【边界】x 为空或全是 nan 时，画一张空图。 【numpy 小知识】np.histogram2d(y, x, bins=..., range=...) 返回二维计数矩阵；第一维对应第一个参数。 这里把 y 放在第一维、并把行倒过来（[::-1]），这样 y 大的样本画在上面，和数学里的坐标系一致。 """ x = np.asarray(x, dtype=np.float64) x = x[np.all(np.isfinite(x), axis=1)] inside = np.all(np.abs(x) <= lim, axis=1) H, _, _ = np.histogram2d(x[inside, 1], x[inside, 0], bins=(height, width), range=[[-lim, lim], [-lim, lim]]) H = H[::-1] peak = max(float(H.max()), 1.0) idx = np.clip(np.ceil(np.sqrt(H / peak) * (len(chars) - 1)).astype(int), 0, len(chars) - 1) lines = ["".join(chars[v] for v in row) for row in idx] lines.append(f"框外 {int(np.sum(~inside))} 个") return lines def side_by_side(blocks, titles, gap=3): """ 把几块字符画（每块是 list[str]）并排拼起来，每块上面一行标题。 【参数】blocks：list[list[str]]（每块行数相同）；titles：list[str]；gap：块之间的空格数 【返回】list[str]，可以逐行 print 【Python 小知识】zip(*blocks)：把"每块的第 r 行"凑成一组。* 把列表"拆开"成多个参数传给 zip。 """ widths = [max(len(line) for line in b) for b in blocks] sep = " " * gap lines = [sep.join(pad(t, w) for t, w in zip(titles, widths))] for rows in zip(*blocks): lines.append(sep.join(pad(r, w) for r, w in zip(rows, widths))) return lines def format_report(r): """把 mode_report 的结果写成一行短文字，例如 "8 个模式，高质量 0.95，最少模式 0.81"。""" return f"{r['covered']} 个模式，高质量 {r['hq_frac']:.3f}，最少模式 {r['min_share']:.2f}" def save_scatter_grid(panels, titles, path, lim=3.0, suptitle=None): """ （可选）把几组二维样本画成一排散点图，存成 PNG。没装 matplotlib 时什么也不做，返回 False。 【参数】panels：list of (n, 2) 数组；titles：每张小图的标题；path：保存路径；lim：坐标范围；suptitle：总标题 【返回】是否保存成功 """ plt = get_pyplot() if plt is None: return False fig, axes = plt.subplots(1, len(panels), figsize=(2.6 * len(panels), 2.9)) axes = np.atleast_1d(axes) for ax, x, t in zip(axes, panels, titles): style_axes(ax) x = np.asarray(x) ax.scatter(x[:, 0], x[:, 1], s=2, alpha=0.4, color=SERIES[0], linewidths=0) ax.scatter(CENTERS[:, 0], CENTERS[:, 1], s=14, marker="x", color=SERIES[1], linewidths=1) ax.set_xlim(-lim, lim) ax.set_ylim(-lim, lim) ax.set_aspect("equal") ax.set_title(t, fontsize=9) ax.set_xticks([-2, 0, 2]) ax.set_yticks([-2, 0, 2]) if suptitle: fig.suptitle(suptitle, fontsize=10) fig.tight_layout() fig.savefig(path, dpi=130) plt.close(fig) return True # ============================================================================= # 自检 # ============================================================================= if __name__ == "__main__": print_header("ring_world 自检：真实数据和几个'假生成器'的指标") rng = np.random.default_rng(0) real = sample_ring(2000, rng) real2 = sample_ring(2000, rng) print(f"8 个中心（前 3 个）：{np.round(CENTERS[:3], 3).tolist()}；相邻两个中心相距 " f"{np.linalg.norm(CENTERS[0] - CENTERS[1]):.3f}，每团的 3σ 半径 {HQ_RADIUS:.2f}") print(f"理论上真实数据的高质量比例 = 1 − e^(−4.5) = {1 - math.exp(-4.5):.4f}") rows = [("真实数据（另一批 2000 个）", real2), ("只生成第 0 团（完全坍缩）", (CENTERS[0] + SIGMA * rng.standard_normal((2000, 2)))), ("全部落在中心 (0, 0)", np.zeros((2000, 2)))] print(pad("样本", 30) + pad("覆盖", 6, ">") + pad("高质量", 9, ">") + pad("最少模式", 10, ">") + pad("FD", 9, ">")) for name, x in rows: r = mode_report(x) print(pad(name, 30) + f"{r['covered']:>6}{r['hq_frac']:>9.3f}{r['min_share']:>10.2f}{frechet_2d(real, x):>9.3f}") print() for line in side_by_side([ascii_density(real)], ["真实数据 2000 个点"]): print(line)