"""
=============================================================================
CS230 · 第 11 讲工具模块：一个合成的"类胸片"医学影像世界（代码 05 用）
=============================================================================

对应讲义：lectures/lecture11_full_cycle_dl_project.md §9（AI + 医疗）

-----------------------------------------------------------------------------
【先说清楚】
-----------------------------------------------------------------------------
这**不是**真实的医学影像，也不模拟任何真实疾病的影像学表现。它只是一个 16×16 的玩具：
"身体"背景上有两块暗的"肺野"，"病灶"是肺野里一块亮斑。我们用它来演示医疗 AI 项目里三个和数据有关的难题：
专家标签有噪声、"超过医生"这类结论怎样读、换一家医院性能会怎样变。数字只说明机制，不代表任何真实系统。

-----------------------------------------------------------------------------
【世界是怎么生成的】
-----------------------------------------------------------------------------
每个"病人"有一个看不见的严重程度 s ∈ [0, 1]（s = U^1.7，U 在 [0, 1] 上均匀；约 1/3 的病人 s > 0.5）。
**真相**：s > 0.5 记为"阳性"。现实中没有人知道真相，这里我们知道，所以能检查每一种"参考标准"离真相有多远。

    图像：身体背景（亮度 0.50–0.60）+ 两块暗的肺野（亮度 0.25–0.35，位置随机抖动 ±1 像素）
          + 一个病灶：肺野里随机位置的高斯亮斑，幅度 0.30·s，半径 1 + s（s 越大越亮、越大；s 很小时几乎看不见）
          + 像素噪声 + 扫描仪（增益、偏置、噪声）
    读片的医生：医生 r 看到的严重程度 = s + 偏差 bias_r + 噪声 N(0, σ_r²)，超过 0.5 就判阳性。
          5 位医生的 (bias, σ) 不同：有人偏保守、有人偏积极、有人更稳定。s 在 0.5 附近的病例，医生之间会有分歧。
    医院：hospital="A"（训练用）和 "B"（外部测试）的扫描仪不同（增益、偏置、噪声），病人的阳性比例不同；
          另外 A 院的重症病人常用**床旁便携机**拍，片子左上角会有一个亮的标记（2×2 像素）——
          这和病情相关，却不是病本身（一种"捷径"，Zech et al. 2018 报告过真实的类似现象）。B 院所有病人都有一半概率带这个标记。

【形状约定】make_patients(...) 返回 dict：
    "X"      : (m, 16, 16) float32
    "s"      : (m,) 真实严重程度；"y_true"：(m,) int64，s > 0.5
    "reads"  : (m, 5) int64，5 位医生各自的判断（0/1）
    "marker" : (m,) bool，片子上有没有便携机标记

【注释标签】【思路】【参数】【返回】【边界】【DL 联系】【numpy 小知识】

直接运行本文件会打印几张示意图和基本统计：
    python3 clinic_world.py          # < 1 秒
"""

import numpy as np

IMG = 16
YS, XS = np.mgrid[0:IMG, 0:IMG].astype(np.float64)

# 5 位"医生"：(偏差, 噪声标准差)。偏差 > 0 表示更容易判阳性
READERS = [(-0.04, 0.07), (0.00, 0.09), (0.03, 0.10), (-0.02, 0.12), (0.05, 0.14)]

# 两家医院：扫描仪（增益、偏置、噪声）、阳性病人的大致比例（通过 s 的分布控制）、便携机标记的规则
HOSPITALS = {
    "A": {"gain": 1.0, "offset": 0.0, "noise": 0.04, "power": 1.7, "marker_if_pos": 0.7, "marker_if_neg": 0.1},
    "B": {"gain": 0.85, "offset": 0.06, "noise": 0.05, "power": 2.6, "marker_if_pos": 0.5, "marker_if_neg": 0.5},
}


def make_patients(m, rng, hospital="A"):
    """
    生成 m 个病人的片子、真实严重程度和 5 位医生的读片结果。

    【思路】先抽 s，再画图（背景 → 肺野 → 病灶 → 便携机标记 → 噪声 → 扫描仪），最后模拟 5 位医生读片。
    【参数】m：病人数；rng：np.random.Generator；hospital："A" 或 "B"
    【返回】dict，见文件开头
    【边界】power 越大，s 越集中在 0 附近，阳性比例越低（A 院约 1/3，B 院约 1/4）。
    【numpy 小知识】(XS - cx[:, None, None]) ** 2：(16,16) 和 (m,1,1) 广播成 (m,16,16)，一次算出 m 张图上每个像素到中心的距离。
    """
    h = HOSPITALS[hospital]
    s = rng.uniform(size=m) ** h["power"]
    y_true = (s > 0.5).astype(np.int64)
    body = rng.uniform(0.50, 0.60, size=m)
    lung = rng.uniform(0.25, 0.35, size=m)
    jx, jy = rng.uniform(-1, 1, size=m), rng.uniform(-1, 1, size=m)
    img = np.broadcast_to(body[:, None, None], (m, IMG, IMG)).copy()
    for cx0 in (4.5, 10.5):                                       # 左肺、右肺：竖着的椭圆
        inside = ((XS - (cx0 + jx[:, None, None])) / 2.6) ** 2 + ((YS - (7.5 + jy[:, None, None])) / 5.5) ** 2 < 1
        img = np.where(inside, lung[:, None, None], img)
    # 病灶：随机选左肺或右肺，在肺野里随机放一个亮斑
    side = np.where(rng.uniform(size=m) < 0.5, 4.5, 10.5)
    cx = side + jx + rng.uniform(-1.2, 1.2, size=m)
    cy = 7.5 + jy + rng.uniform(-3.5, 3.5, size=m)
    r = 1.0 + s
    blob = np.exp(-((XS - cx[:, None, None]) ** 2 + (YS - cy[:, None, None]) ** 2) / (2 * r[:, None, None] ** 2))
    img = img + (0.30 * s)[:, None, None] * blob
    # 便携机标记：左上角 2×2 的亮块
    p_marker = np.where(y_true == 1, h["marker_if_pos"], h["marker_if_neg"])
    marker = rng.uniform(size=m) < p_marker
    img[:, 0:2, 0:2] = np.where(marker[:, None, None], 0.95, img[:, 0:2, 0:2])
    img = img + rng.normal(0.0, h["noise"], size=(m, IMG, IMG))
    img = np.clip(h["gain"] * img + h["offset"], 0.0, 1.0)
    # 5 位医生读片
    reads = np.stack([(s + b + rng.normal(0.0, sd, size=m) > 0.5) for b, sd in READERS], axis=1).astype(np.int64)
    return {"X": img.astype(np.float32), "s": s, "y_true": y_true, "reads": reads, "marker": marker}


def majority(reads):
    """
    多数投票：reads 形状 (m, k)，返回 (m,) 的 0/1。k 为偶数时平票记为阳性（这里只用奇数个医生，不会平票）。
    """
    return (reads.mean(axis=1) >= 0.5).astype(np.int64)


def ascii_xray(img, chars=" .:-=+*#%@"):
    """固定亮度范围 [0.2, 0.9] 的字符画，返回 16 个字符串。"""
    n = len(chars)
    idx = np.clip(((np.asarray(img) - 0.2) / 0.7 * n).astype(int), 0, n - 1)
    return ["".join(chars[v] for v in row) for row in idx]


if __name__ == "__main__":
    d = make_patients(4000, np.random.default_rng(0), "A")
    dB = make_patients(4000, np.random.default_rng(1), "B")
    print(f"A 院：阳性比例 {d['y_true'].mean():.3f}，带便携机标记 {d['marker'].mean():.3f}"
          f"（阳性中 {d['marker'][d['y_true'] == 1].mean():.2f}，阴性中 {d['marker'][d['y_true'] == 0].mean():.2f}）")
    print(f"B 院：阳性比例 {dB['y_true'].mean():.3f}，带便携机标记 {dB['marker'].mean():.3f}")
    acc = [(d["reads"][:, k] == d["y_true"]).mean() for k in range(len(READERS))]
    print("5 位医生相对真相的准确率：" + "、".join(f"{a:.3f}" for a in acc))
    idx = [int(np.flatnonzero(d["s"] < 0.1)[0]), int(np.flatnonzero((d["s"] > 0.45) & (d["s"] < 0.55))[0]),
           int(np.flatnonzero(d["s"] > 0.8)[0])]
    panels = [ascii_xray(d["X"][i]) for i in idx]
    print("\n三张示意图（s ≈ " + "、".join(f"{d['s'][i]:.2f}" for i in idx) + "）：")
    for row in range(IMG):
        print("   ".join(p[row] for p in panels))