#!/usr/bin/env python3
"""结直肠癌腹膜转移文献：2015–2019 vs 2020–2026 的主题词占比变化。

判定规则写在 ../claims.json，跑之前定死。零模型是年份标签随机打乱 1000 次。
"""
import json, collections
import numpy as np
from scipy.stats import fisher_exact

N_PERM, MIN_DOCS, SEED = 1000, 20, 20260827


def bh(p):
    p = np.asarray(p, float); n = len(p); o = np.argsort(p); adj = np.empty(n); prev = 1.0
    for r, i in enumerate(o[::-1]):
        prev = min(prev, p[i] * n / (n - r)); adj[i] = prev
    return adj


def screen(recs, labels, terms):
    """labels: True=后期。返回每个词的 (p, 前期占比, 后期占比)。"""
    n_late = int(labels.sum()); n_early = len(labels) - n_late
    out = []
    for t in terms:
        has = np.array([t in r["mesh"] for r in recs])
        a = int((has & labels).sum())          # 后期有
        b = int((has & ~labels).sum())         # 前期有
        _, p = fisher_exact([[a, n_late - a], [b, n_early - b]])
        out.append((t, float(p), b / n_early, a / n_late))
    return out


def main():
    d = json.load(open("corpus.json"))
    recs = [r for r in d["records"] if 2015 <= r["year"] <= 2026]
    labels = np.array([r["year"] >= 2020 for r in recs])
    print(f"[data] 文献 {len(recs)} 篇｜2015–2019 {int((~labels).sum())} 篇｜2020–2026 {int(labels.sum())} 篇")

    cnt = collections.Counter(m for r in recs for m in r["mesh"])
    terms = [t for t, c in cnt.items() if c >= MIN_DOCS]
    print(f"[data] 至少出现在 {MIN_DOCS} 篇里的 MeSH 词: {len(terms)} 个")

    obs = screen(recs, labels, terms)
    q = bh([x[1] for x in obs])
    sig = [(t, p, e, l, qq) for (t, p, e, l), qq in zip(obs, q) if qq < 0.05]
    sig.sort(key=lambda x: -(x[3] - x[2]))
    print(f"[实测] FDR<0.05 的主题词: {len(sig)} 个")

    rng = np.random.default_rng(SEED)
    null = []
    for _ in range(N_PERM):
        lab = labels.copy(); rng.shuffle(lab)
        pp = [x[1] for x in screen(recs, lab, terms)]
        null.append(int((bh(pp) < 0.05).sum()))
    null = np.array(null)
    p95 = float(np.percentile(null, 95))
    h1 = bool(len(sig) > p95 and len(sig) >= 5)
    print(f"[零模型] {N_PERM} 次打乱：均值 {null.mean():.2f}｜95分位 {p95:.1f}｜最大 {null.max()}")
    print(f"[H1] {len(sig)} > {p95:.1f} 且 >=5 → {'✅成立' if h1 else '❌被推翻'}")

    rising = [s for s in sig if s[3] > s[2]][:12]
    falling = sorted([s for s in sig if s[3] < s[2]], key=lambda x: x[3] - x[2])[:12]
    print("\n  升幅最大：")
    for t, p, e, l, qq in rising[:6]:
        print(f"    {t[:44]:<46} {e*100:5.1f}% → {l*100:5.1f}%  (FDR {qq:.1e})")
    print("  降幅最大：")
    for t, p, e, l, qq in falling[:6]:
        print(f"    {t[:44]:<46} {e*100:5.1f}% → {l*100:5.1f}%  (FDR {qq:.1e})")

    HOT = ("Immunotherapy", "Single-Cell", "Molecular Targeted", "Immune Checkpoint",
           "Antineoplastic Combined", "Tumor Microenvironment", "Biomarkers")
    hot_hits = [s[0] for s in rising if any(h.lower() in s[0].lower() for h in HOT)]
    h2 = len(hot_hits) > 0
    print(f"\n[H2] 升幅词里含免疫/单细胞/靶向类新兴主题 → {'✅成立' if h2 else '❌被推翻'}｜命中 {hot_hits}")

    json.dump({"seed": SEED, "n_docs": len(recs), "n_early": int((~labels).sum()),
               "n_late": int(labels.sum()), "n_terms": len(terms), "min_docs": MIN_DOCS,
               "n_significant": len(sig), "null_mean": float(null.mean()), "null_p95": p95,
               "null_max": int(null.max()),
               "empirical_p": float((null >= len(sig)).mean()),
               "H1_passes": h1, "H2_passes": h2, "hot_hits": hot_hits,
               "rising": [{"term": t, "early": e, "late": l, "fdr": qq} for t, p, e, l, qq in rising],
               "falling": [{"term": t, "early": e, "late": l, "fdr": qq} for t, p, e, l, qq in falling],
               "null_distribution": null.tolist()},
              open("results.json", "w"), ensure_ascii=False, indent=1)


if __name__ == "__main__":
    main()