#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
穩定度結果的診斷與分層
======================
    python _pipeline\\穩定度診斷.py

只讀 40_案件歷程/合議庭穩定度_逐案.csv，不重掃終結案件資料，數十秒可跑完。

要處理的三個問題：

D1 分母不對稱
   逐案指標的分母是「掃描時保留下來的案件」。被聲請案件的法官全部是目標法官，
   其全部案件都被保留，分母完整；對照案件若全由非目標法官組成，只有落在
   「被聲請案件之法院×年度」範圍內的案件被保留，分母被截斷，穩定度會被低估。
   本診斷把對照組切成「含目標法官」與「不含目標法官」兩層分別比較，
   前者與被聲請組的分母條件相同，是乾淨的對照。

D2 法院異質性
   最高法院與高等法院的庭運作方式差異極大，合併比較會互相稀釋。
   本診斷逐法院分層重跑。

D3 合議庭人數
   僅列一名法官者無合議可言，panel_n 與 pair_mean 失去意義，應排除。

輸出於 40_案件歷程/：
    穩定度診斷.csv
並更新 00_索引/穩定度分析報告.md（附加「五、診斷與分層」）
"""

import csv, sys, random, math, collections
from pathlib import Path

csv.field_size_limit(10 ** 9)
random.seed(20260727)

ROOT = Path(__file__).resolve().parent.parent
D = ROOT / "40_案件歷程"
M = ["panel_n", "pair_mean", "hhi_mean"]


try:
    import numpy as np
except ImportError:
    np = None


def med(v):
    v = sorted(v)
    if not v:
        return None
    m = len(v) // 2
    return v[m] if len(v) % 2 else (v[m - 1] + v[m]) / 2


def mwu(a, b):
    """Mann-Whitney U 檢定（常態近似，含結值校正）。
    回傳 (p 雙尾, rank-biserial 效果量)。
    效果量 r 介於 -1 與 1，正值表示 A 組傾向大於 B 組，0.1／0.3／0.5 約為小／中／大。
    對偏態且樣本數懸殊的資料，秩檢定比中位數置換穩健，且不需重抽樣。"""
    n1, n2 = len(a), len(b)
    if n1 < 3 or n2 < 3:
        return None, None
    allv = a + b
    order = sorted(range(len(allv)), key=lambda i: allv[i])
    ranks = [0.0] * len(allv)
    i = 0
    ties = 0.0
    while i < len(order):
        j = i
        while j + 1 < len(order) and allv[order[j + 1]] == allv[order[i]]:
            j += 1
        r = (i + j) / 2.0 + 1
        t = j - i + 1
        if t > 1:
            ties += t ** 3 - t
        for k in range(i, j + 1):
            ranks[order[k]] = r
        i = j + 1
    r1 = sum(ranks[:n1])
    u1 = r1 - n1 * (n1 + 1) / 2.0
    mu = n1 * n2 / 2.0
    n = n1 + n2
    sd2 = n1 * n2 / 12.0 * ((n + 1) - ties / (n * (n - 1))) if n > 1 else 0
    if sd2 <= 0:
        return None, None
    z = (u1 - mu) / (sd2 ** 0.5)
    p = 2 * (1 - 0.5 * (1 + math.erf(abs(z) / (2 ** 0.5))))
    rb = 2 * u1 / (n1 * n2) - 1
    return round(max(p, 1e-16), 5), round(rb, 3)


def hl_shift(a, b, cap=400):
    """Hodges-Lehmann 位移估計（兩組差的中位數）。樣本過大時隨機抽樣以控制成本。"""
    if len(a) < 3 or len(b) < 3:
        return None
    sa = a if len(a) <= cap else random.sample(a, cap)
    sb = b if len(b) <= cap else random.sample(b, cap)
    if np is not None:
        d = (np.array(sa)[:, None] - np.array(sb)[None, :]).ravel()
        return round(float(np.median(d)), 3)
    d = [x - y for x in sa for y in sb]
    return round(med(d), 3)


def van_elteren(strata):
    """分層 Wilcoxon（van Elteren 檢定），權重 1/(N_i+1)。
    每一層是一個 (法院, 年度, 案由) 格，層內比較被聲請與對照，再跨層合成。
    這是處理「各層基準值不同、且兩組在各層佔比不同」時的正確做法；
    直接合併比較會產生組成效應（Simpson 悖論），把層內的無差異變成整體的假差異。
    回傳 (p 雙尾, z, 有效層數, 加權層內 HL 位移中位數)。"""
    num = 0.0
    den = 0.0
    used = 0
    shifts = []
    for a, b in strata:
        n1, n2 = len(a), len(b)
        if n1 < 1 or n2 < 1:
            continue
        N = n1 + n2
        allv = a + b
        order = sorted(range(N), key=lambda i: allv[i])
        ranks = [0.0] * N
        i = 0
        ties = 0.0
        while i < N:
            j = i
            while j + 1 < N and allv[order[j + 1]] == allv[order[i]]:
                j += 1
            r = (i + j) / 2.0 + 1
            t = j - i + 1
            if t > 1:
                ties += t ** 3 - t
            for k in range(i, j + 1):
                ranks[order[k]] = r
            i = j + 1
        W = sum(ranks[:n1])
        E = n1 * (N + 1) / 2.0
        V = n1 * n2 / 12.0 * ((N + 1) - ties / (N * (N - 1))) if N > 1 else 0.0
        if V <= 0:
            continue
        w = 1.0 / (N + 1)
        num += w * (W - E)
        den += w * w * V
        used += 1
        shifts.append((med(a) - med(b), N))
    if used < 2 or den <= 0:
        return None, None, used, None
    z = num / (den ** 0.5)
    p = 2 * (1 - 0.5 * (1 + math.erf(abs(z) / (2 ** 0.5))))
    tot = sum(n for _, n in shifts)
    wsh = sum(d * n for d, n in shifts) / tot if tot else None
    return round(max(p, 1e-16), 5), round(z, 3), used, (round(wsh, 3) if wsh is not None else None)


def main():
    try:
        sys.stdout.reconfigure(encoding="utf-8")
    except Exception:
        pass
    src = D / "合議庭穩定度_逐案.csv"
    if not src.exists():
        print("找不到", src, "\n請先執行 _pipeline\\合議庭穩定度.py")
        sys.exit(1)

    targets = {r["法官"].strip() for r in csv.DictReader(open(D / "法官參與.csv", encoding="utf-8-sig"))
               if r.get("法官", "").strip()}
    outcome = {}
    if (D / "釋憲結果.csv").exists():
        for r in csv.DictReader(open(D / "釋憲結果.csv", encoding="utf-8-sig")):
            outcome[r["釋憲案"]] = r["R碼"]

    rows = []
    n = 0
    for r in csv.DictReader(open(src, encoding="utf-8-sig")):
        n += 1
        try:
            js = [x for x in r["法官"].split("、") if x]
            rows.append({"法院": r["法院"], "年": r["年"], "案由": r["案由"],
                         "pet": r["被聲請"] == "1", "nj": int(r["n_judge"]),
                         "tgt": any(j in targets for j in js),
                         "釋憲案": r["釋憲案"],
                         "panel_n": float(r["panel_n"]), "pair_mean": float(r["pair_mean"]),
                         "hhi_mean": float(r["hhi_mean"])})
        except (KeyError, ValueError):
            continue
        if n % 500000 == 0:
            print(f"      讀入 {n}")
    print(f"讀入 {len(rows)} 列")

    P = [r for r in rows if r["pet"]]
    cells = {(r["法院"], r["年"], r["案由"]) for r in P}
    C = [r for r in rows if not r["pet"] and (r["法院"], r["年"], r["案由"]) in cells]
    print(f"被聲請 {len(P)}；同格對照 {len(C)}")

    out = []

    def add(tag, a, b, na, nb):
        row = {"比較": tag, "A組": na, "A件數": len(a), "B組": nb, "B件數": len(b)}
        for m in M:
            va = [x[m] for x in a]
            vb = [x[m] for x in b]
            row[f"{m}_A中位"] = med(va)
            row[f"{m}_B中位"] = med(vb)
            row[f"{m}_HL位移"] = hl_shift(va, vb)
            p, rb = mwu(va, vb)
            row[f"{m}_效果量r"] = rb
            row[f"{m}_p"] = p
        out.append(row)

    # 全體（重現原報告）
    add("原比較", P, C, "被聲請", "同格對照")

    # D3 僅合議（三人以上）
    P3 = [r for r in P if r["nj"] >= 3]
    C3 = [r for r in C if r["nj"] >= 3]
    add("D3 僅三人以上合議", P3, C3, "被聲請", "同格對照")

    # D1 對照組依是否含目標法官分層
    Ct = [r for r in C3 if r["tgt"]]
    Cn = [r for r in C3 if not r["tgt"]]
    add("D1 對照僅含目標法官（分母條件相同）", P3, Ct, "被聲請", "對照(含目標法官)")
    add("D1 對照不含目標法官（分母被截斷）", P3, Cn, "被聲請", "對照(無目標法官)")
    add("D1 偏誤大小：兩類對照互比", Ct, Cn, "對照(含目標)", "對照(無目標)")

    # D2 逐法院
    bycourt = collections.Counter(r["法院"] for r in P3)
    for court, cnt in bycourt.most_common(6):
        if cnt < 8:
            continue
        a = [r for r in P3 if r["法院"] == court]
        b = [r for r in Ct if r["法院"] == court]
        if len(b) < 8:
            b = [r for r in C3 if r["法院"] == court]
            add(f"D2 {court}（對照未分層）", a, b, "被聲請", "同格對照")
        else:
            add(f"D2 {court}（對照含目標法官）", a, b, "被聲請", "對照(含目標法官)")

    # 釋憲結果分組（僅三人以上合議）
    VIO, CON = {"R1", "R2", "R3", "R6"}, {"R4", "R5"}
    gv, gc = [], []
    for r in P3:
        codes = {outcome[k] for k in r["釋憲案"].split("、") if k in outcome}
        if codes & VIO:
            gv.append(r)
        elif codes & CON:
            gc.append(r)
    add("釋憲結果 違憲類 vs 合憲類（僅合議）", gv, gc, "違憲類", "合憲類")

    # D4 分層檢定：以 (法院,年,案由) 為層，處理組成效應
    ve = []
    for label, ctrl in (("對照含目標法官", Ct), ("全部同格對照", C3)):
        bycell_p = collections.defaultdict(list)
        bycell_c = collections.defaultdict(list)
        for r in P3:
            bycell_p[(r["法院"], r["年"], r["案由"])].append(r)
        for r in ctrl:
            bycell_c[(r["法院"], r["年"], r["案由"])].append(r)
        keys = [k for k in bycell_p if k in bycell_c]
        for m in M:
            strata = [([x[m] for x in bycell_p[k]], [x[m] for x in bycell_c[k]]) for k in keys]
            p, z, used, wsh = van_elteren(strata)
            ve.append({"對照": label, "指標": m, "有效層數": used,
                       "加權層內位移": wsh, "z": z, "p": p})
    with open(D / "穩定度分層檢定.csv", "w", encoding="utf-8-sig", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=["對照", "指標", "有效層數", "加權層內位移", "z", "p"])
        w.writeheader()
        for r in ve:
            w.writerow(r)
    print("  --- D4 分層（van Elteren）---")
    for r in ve:
        print(f"  {r['對照']:14s} {r['指標']:10s} 層數={r['有效層數']:3d} 層內位移={r['加權層內位移']} z={r['z']} p={r['p']}")

    fields = ["比較", "A組", "A件數", "B組", "B件數"]
    for m in M:
        fields += [f"{m}_A中位", f"{m}_B中位", f"{m}_HL位移", f"{m}_效果量r", f"{m}_p"]
    with open(D / "穩定度診斷.csv", "w", encoding="utf-8-sig", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=fields)
        w.writeheader()
        for r in out:
            w.writerow(r)

    rep = ["", "---", "", "## 五、診斷與分層", "",
           "指標分母為掃描時保留之案件。被聲請案件的法官全屬目標法官，分母完整；",
           "對照案件若全由非目標法官組成，分母被截斷，穩定度會被低估。",
           "下表把對照組切成「含目標法官」與「不含目標法官」兩層，前者與被聲請組分母條件相同。",
           "另排除僅列一名法官之案件，並逐法院分層。", "",
           "| 比較 | A組 | A件數 | B組 | B件數 | panel_n HL位移 | r | p | hhi_mean HL位移 | r | p |",
           "|---|---|---|---|---|---|---|---|---|---|---|"]
    for r in out:
        rep.append("| {比較} | {A組} | {A件數} | {B組} | {B件數} | {a} | {b} | {c} | {d} | {e} | {f} |".format(
            **r, a=r["panel_n_HL位移"], b=r["panel_n_效果量r"], c=r["panel_n_p"],
            d=r["hhi_mean_HL位移"], e=r["hhi_mean_效果量r"], f=r["hhi_mean_p"]))
    rep += ["", "### D4 分層檢定（van Elteren）", "",
            "以 (法院, 年度, 案由) 為層，層內比較後跨層合成。合併比較會受組成效應影響："
            "兩組在各層的佔比不同，而各層基準值又不同，就會把層內的無差異變成整體的假差異。", "",
            "| 對照 | 指標 | 有效層數 | 加權層內位移 | z | p |", "|---|---|---|---|---|---|"]
    for r in ve:
        rep.append(f"| {r['對照']} | {r['指標']} | {r['有效層數']} | {r['加權層內位移']} | {r['z']} | {r['p']} |")
    rep += ["", "完整三項指標見 `40_案件歷程/穩定度診斷.csv` 與 `穩定度分層檢定.csv`。", "",
            "註：HL 位移為 Hodges-Lehmann 兩組差之中位數（樣本超過 400 時隨機抽樣估計）；",
            "r 為 rank-biserial 效果量，介於 -1 與 1，正值表示 A 組傾向較大，0.1／0.3／0.5 約為小／中／大；",
            "p 為 Mann-Whitney U 檢定（常態近似，含結值校正，雙尾）。",
            "樣本數懸殊時 p 值容易因對照組龐大而顯著，請以效果量 r 為主要判準。"]
    f = ROOT / "00_索引" / "穩定度分析報告.md"
    f.write_text(f.read_text(encoding="utf-8") + "\n".join(rep) + "\n", encoding="utf-8")

    print("=" * 60)
    for r in out:
        print(f"  {r['比較'][:32]:34s} panel_n HL={r['panel_n_HL位移']} r={r['panel_n_效果量r']} p={r['panel_n_p']}   hhi r={r['hhi_mean_效果量r']} p={r['hhi_mean_p']}")
    print("報告已附加第五節：", f)
    print("=" * 60)


if __name__ == "__main__":
    main()