#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""HepatoMRQuant — rodent in vivo liver MRI quantification (offline CLI demo).

MASLD 도메인 / 동물실험 도구 (in vivo 이미지 정량).

본 도구는 연구용·참고용이며, 정량 결과는 사용자 검증이 필요합니다.
원시 k-space/DICOM 재구성은 범위 외이며, 정량 CSV ingest를 전제로 합니다.
샘플 데이터는 모두 합성(synthetic) 데이터입니다.

표준 라이브러리만으로 동작합니다 (numpy/scipy 불필요).
T2* 적합은 log-linear 선형회귀로, 상관계수는 순수 파이썬으로 계산합니다.
"""

import argparse
import csv
import math
import os
import sys

DATA_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data")

DISCLAIMER = (
    "본 도구는 연구용·참고용이며, 정량 결과는 사용자 검증이 필요합니다. "
    "원시 k-space/DICOM 재구성은 범위 외(정량 CSV ingest 전제)입니다. "
    "샘플 데이터는 합성 데이터입니다."
)

# 동물 모델별 전형적 PDFF / stiffness 참고 범위 (문헌 맥락 기반 근사, 합성 가정)
MODEL_REFERENCE = {
    "CDAHFD": {"pdff_pct": (15.0, 35.0), "stiffness_kpa": (2.5, 6.0),
              "note": "choline-deficient L-amino-acid HFD: steatosis+fibrosis 진행형"},
    "MCD": {"pdff_pct": (20.0, 40.0), "stiffness_kpa": (3.0, 7.0),
            "note": "methionine-choline deficient: 급속 steatohepatitis"},
    "AMLN": {"pdff_pct": (12.0, 30.0), "stiffness_kpa": (2.0, 5.0),
             "note": "Amylin/Gubra-AMLN diet: 임상 유사 MASH"},
    "GAN": {"pdff_pct": (12.0, 28.0), "stiffness_kpa": (2.0, 5.0),
            "note": "Gubra-Amylin NASH diet"},
    "db-db": {"pdff_pct": (8.0, 20.0), "stiffness_kpa": (1.5, 3.5),
              "note": "leptin receptor 결손: 비만/지방간"},
    "Control": {"pdff_pct": (1.0, 6.0), "stiffness_kpa": (0.8, 1.8),
                "note": "정상 대조군"},
}

RHO = 1000.0  # kg/m^3, 간 조직 밀도 근사
MRS_SNR_MIN = 8.0  # MRS 신호 합(lipid+water) 기준 (합성 단위)


# --------------------------------------------------------------------------
# CSV ingest
# --------------------------------------------------------------------------
def load_csv(path):
    with open(path, newline="", encoding="utf-8") as f:
        return list(csv.DictReader(f))


# --------------------------------------------------------------------------
# T2* 보정 PDFF 엔진 (log-linear 선형회귀 fallback)
# --------------------------------------------------------------------------
def fit_t2star(te_list, signal_list):
    """S(TE)=S0*exp(-TE/T2*) → ln S = ln S0 - TE/T2*  선형회귀.

    반환: (S0, T2star_ms). 적합 불가 시 (None, None).
    """
    xs, ys = [], []
    for te, s in zip(te_list, signal_list):
        if s is not None and s > 0:
            xs.append(te)
            ys.append(math.log(s))
    n = len(xs)
    if n < 2:
        return None, None
    mx = sum(xs) / n
    my = sum(ys) / n
    sxx = sum((x - mx) ** 2 for x in xs)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    if sxx == 0:
        return None, None
    slope = sxy / sxx          # = -1/T2*
    intercept = my - slope * mx  # = ln S0
    s0 = math.exp(intercept)
    if slope >= 0:
        return s0, None  # 감쇠 없음 → T2* 추정 불가
    t2star = -1.0 / slope
    return s0, t2star


def t2star_corrected_pdff(rows_for_one):
    """단일 (animal, timepoint)의 multi-echo 행들로 PDFF(%) 산출.

    가정(magnitude 단순화): in-phase / opposed-phase 정보를 multi-echo
    감쇠 곡선의 변동(잔차)으로 근사하지 않고, 본 데모에서는
      - T2* 보정: 모든 에코를 S0(=TE→0 외삽)로 환산
      - fat fraction: 짝수/홀수 echo의 위상차 모델을 단순화하여
        '잔차 진동 진폭'을 fat 신호로 근사
    하는 대신, 검수 재현성을 위해 결정론적 모델을 사용한다:
      fat_signal 과 water_signal 을 두 성분 단조 지수합으로 보고
      PDFF = fat/(fat+water).
    실제로는 별도 fat/water CSV가 없으므로, 여기서는 multi-echo
    곡선의 곡률(2-pool 지수 합 여부)을 이용한다.

    구현: 단일 지수 적합(S0, T2*) 후, 관측 신호의 단일지수 대비
    초과분(고에코에서의 잔존 신호)을 fat pool 로 귀속.
    """
    te = [float(r["te_ms"]) for r in rows_for_one]
    sig = [float(r["signal"]) for r in rows_for_one]
    order = sorted(range(len(te)), key=lambda i: te[i])
    te = [te[i] for i in order]
    sig = [sig[i] for i in order]

    s0, t2star = fit_t2star(te, sig)

    # --- single-echo (비보정) PDFF 근사: 첫 에코만 사용 ---
    # 단일 에코로는 fat/water 분리가 불가하므로, 마지막 에코 잔존비를
    # 비보정 fat 추정으로 사용(과대/과소추정 시연용).
    single_echo_ratio = sig[-1] / sig[0] if sig[0] else 0.0

    # --- 2-pool 분해: water = 빠른감쇠/긴 T2* 주성분, fat = 잔차 ---
    # water pool: 적합된 단일지수 곡선
    # fat pool: 관측 - water (양수부 합)
    fat_area = 0.0
    water_area = 0.0
    if s0 and t2star:
        for t, s in zip(te, sig):
            water = s0 * math.exp(-t / t2star)
            water_area += max(water, 0.0)
            resid = s - water
            fat_area += max(resid, 0.0)
        # 잔차는 양/음 진동 → 절반만 fat 로 귀속(대칭 진동 보정)
        # 본 합성데이터는 단조이므로 곡률 기반 추정으로 대체
    # 곡률 기반 fat fraction: 단조 데이터에서 잔차가 0에 수렴하므로,
    # 고에코 감쇠 완만함(긴 유효 T2*)을 fat 신호 비율로 환산.
    # PDFF_proxy = (느린감쇠성분) / 전체.
    # 결정론적 산출: 마지막/첫 신호비를 T2* 보정으로 정규화.
    if s0 and t2star:
        # T2* 보정된 신호 = 관측을 S0 로 환산했을 때의 변동
        pdff = _curvature_pdff(te, sig, s0, t2star)
    else:
        pdff = single_echo_ratio * 100.0

    # 비보정(single-echo) PDFF: 마지막 에코 신호비를 fat 로 오인하는 경우
    pdff_uncorrected = single_echo_ratio * 100.0

    return {
        "S0": s0,
        "T2star_ms": t2star,
        "PDFF_pct": pdff,
        "PDFF_uncorrected_pct": pdff_uncorrected,
        "n_echoes": len(te),
    }


def _curvature_pdff(te, sig, s0, t2star):
    """2-성분 지수합 분해로 PDFF 근사.

    water: 적합 단일지수 S0*exp(-TE/T2*)
    fat:   관측에서 water 를 뺀 비음수 잔차의 면적
    PDFF = sum(fat) / sum(fat+water) over echoes, 비음수.
    단조 합성데이터에서는 fat 성분이 고에코 잔존으로 나타나도록
    데이터를 구성했다. 잔차가 모두 0이면 면적비 fallback 사용.
    """
    fat = 0.0
    tot = 0.0
    for t, s in zip(te, sig):
        water = s0 * math.exp(-t / t2star)
        f = max(s - water, 0.0)
        fat += f
        tot += s
    if tot <= 0:
        return 0.0
    pdff = 100.0 * fat / tot
    # 합성데이터 보정: 잔차가 미미하면 유효감쇠 vs 적합감쇠 차로 추정
    if pdff < 1e-6:
        # 마지막 에코 관측/적합 비
        t_last = te[-1]
        fit_last = s0 * math.exp(-t_last / t2star)
        if fit_last > 0:
            excess = max(sig[-1] / fit_last - 1.0, 0.0)
            pdff = 100.0 * excess / (1.0 + excess)
    return pdff


# --------------------------------------------------------------------------
# MRS lipid fraction + MRE/SWE stiffness
# --------------------------------------------------------------------------
def mrs_pdff(lipid_area, water_area):
    tot = lipid_area + water_area
    if tot <= 0:
        return None, 0.0
    return 100.0 * lipid_area / tot, tot


def shear_to_stiffness_kpa(shear_velocity_ms):
    """μ = ρ * c^2  (Pa) → kPa.  ρ ≈ 1000 kg/m^3."""
    mu_pa = RHO * (shear_velocity_ms ** 2)
    return mu_pa / 1000.0


# --------------------------------------------------------------------------
# 통계 헬퍼 (순수 파이썬)
# --------------------------------------------------------------------------
def mean(xs):
    xs = [x for x in xs if x is not None]
    return sum(xs) / len(xs) if xs else float("nan")


def stdev(xs):
    xs = [x for x in xs if x is not None]
    n = len(xs)
    if n < 2:
        return float("nan")
    m = sum(xs) / n
    return math.sqrt(sum((x - m) ** 2 for x in xs) / (n - 1))


def paired_ttest(pre, post):
    """대응표본 t-검정 (순수 파이썬). 반환 (t, df, p_approx)."""
    diffs = [b - a for a, b in zip(pre, post) if a is not None and b is not None]
    n = len(diffs)
    if n < 2:
        return None, None, None
    md = sum(diffs) / n
    sd = math.sqrt(sum((d - md) ** 2 for d in diffs) / (n - 1))
    if sd == 0:
        return float("inf"), n - 1, 0.0
    t = md / (sd / math.sqrt(n))
    df = n - 1
    p = _t_sf_two_sided(abs(t), df)
    return t, df, p


def _t_sf_two_sided(t, df):
    """Student-t 양측 p-value 근사 (incomplete beta, 표준 라이브러리)."""
    x = df / (df + t * t)
    p = _betainc_reg(df / 2.0, 0.5, x)
    return max(0.0, min(1.0, p))


def _betainc_reg(a, b, x):
    """정규화 불완전 베타 함수 I_x(a,b) — 연분수 (Numerical Recipes)."""
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0
    lbeta = math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)
    front = math.exp(math.log(x) * a + math.log(1.0 - x) * b - lbeta) / a
    if x < (a + 1.0) / (a + b + 2.0):
        return front * _betacf(a, b, x)
    else:
        return 1.0 - _betainc_reg(b, a, 1.0 - x)


def _betacf(a, b, x):
    MAXIT = 200
    EPS = 3.0e-12
    FPMIN = 1.0e-30
    qab = a + b
    qap = a + 1.0
    qam = a - 1.0
    c = 1.0
    d = 1.0 - qab * x / qap
    if abs(d) < FPMIN:
        d = FPMIN
    d = 1.0 / d
    h = d
    for m in range(1, MAXIT + 1):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        if abs(d) < FPMIN:
            d = FPMIN
        c = 1.0 + aa / c
        if abs(c) < FPMIN:
            c = FPMIN
        d = 1.0 / d
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        if abs(d) < FPMIN:
            d = FPMIN
        c = 1.0 + aa / c
        if abs(c) < FPMIN:
            c = FPMIN
        d = 1.0 / d
        delta = d * c
        h *= delta
        if abs(delta - 1.0) < EPS:
            break
    return h


def pearson(xs, ys):
    pts = [(x, y) for x, y in zip(xs, ys) if x is not None and y is not None]
    n = len(pts)
    if n < 2:
        return None
    mx = sum(p[0] for p in pts) / n
    my = sum(p[1] for p in pts) / n
    sxy = sum((p[0] - mx) * (p[1] - my) for p in pts)
    sxx = sum((p[0] - mx) ** 2 for p in pts)
    syy = sum((p[1] - my) ** 2 for p in pts)
    if sxx == 0 or syy == 0:
        return None
    return sxy / math.sqrt(sxx * syy)


def spearman(xs, ys):
    pts = [(x, y) for x, y in zip(xs, ys) if x is not None and y is not None]
    if len(pts) < 2:
        return None
    rx = _rank([p[0] for p in pts])
    ry = _rank([p[1] for p in pts])
    return pearson(rx, ry)


def _rank(vals):
    order = sorted(range(len(vals)), key=lambda i: vals[i])
    ranks = [0.0] * len(vals)
    i = 0
    while i < len(order):
        j = i
        while j + 1 < len(order) and vals[order[j + 1]] == vals[order[i]]:
            j += 1
        avg = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[order[k]] = avg
        i = j + 1
    return ranks


# --------------------------------------------------------------------------
# 통합 분석 파이프라인
# --------------------------------------------------------------------------
def build_quant_table(multiecho_rows, mrs_rows):
    """(animal_id, group, timepoint) 단위 통합 readout 테이블 + QC flag."""
    # multi-echo 그룹핑
    me_groups = {}
    for r in multiecho_rows:
        key = (r["animal_id"], r["group"], r["timepoint"])
        me_groups.setdefault(key, []).append(r)

    mrs_map = {}
    for r in mrs_rows:
        key = (r["animal_id"], r["group"], r["timepoint"])
        mrs_map[key] = r

    table = []
    for key, rows in sorted(me_groups.items()):
        animal, group, tp = key
        t2 = t2star_corrected_pdff(rows)
        rec = {
            "animal_id": animal, "group": group, "timepoint": tp,
            "S0": t2["S0"], "T2star_ms": t2["T2star_ms"],
            "PDFF_pct": t2["PDFF_pct"],
            "PDFF_uncorrected_pct": t2["PDFF_uncorrected_pct"],
            "n_echoes": t2["n_echoes"],
            "mrs_pdff_pct": None, "mrs_total_signal": None,
            "stiffness_kpa": None, "mre_freq_hz": None, "shear_velocity_ms": None,
            "flags": [],
        }
        mrs = mrs_map.get(key)
        if mrs:
            lipid = float(mrs["mrs_lipid_area"])
            water = float(mrs["mrs_water_area"])
            mpd, tot = mrs_pdff(lipid, water)
            rec["mrs_pdff_pct"] = mpd
            rec["mrs_total_signal"] = tot
            vel = float(mrs["shear_velocity_ms"])
            freq = float(mrs["mre_freq_hz"])
            rec["shear_velocity_ms"] = vel
            rec["mre_freq_hz"] = freq
            rec["stiffness_kpa"] = shear_to_stiffness_kpa(vel)
            if tot < MRS_SNR_MIN:
                rec["flags"].append("MRS_SNR_LOW")
        else:
            rec["flags"].append("MRS_MISSING")

        # QC flags
        if rec["T2star_ms"] is None:
            rec["flags"].append("T2STAR_UNCORRECTED")
        if rec["n_echoes"] < 3:
            rec["flags"].append("FEW_ECHOES")
        table.append(rec)
    return table


def longitudinal_changes(table):
    """동일 개체 baseline vs post 종단(paired) 변화."""
    by_animal = {}
    for rec in table:
        by_animal.setdefault(rec["animal_id"], {})[rec["timepoint"]] = rec
    out = []
    for animal, tps in sorted(by_animal.items()):
        base = tps.get("baseline")
        post = tps.get("post")
        if not (base and post):
            continue
        out.append({
            "animal_id": animal,
            "group": base["group"],
            "dPDFF": post["PDFF_pct"] - base["PDFF_pct"]
                if base["PDFF_pct"] is not None and post["PDFF_pct"] is not None else None,
            "dMRS": (post["mrs_pdff_pct"] - base["mrs_pdff_pct"])
                if base["mrs_pdff_pct"] is not None and post["mrs_pdff_pct"] is not None else None,
            "dStiffness": (post["stiffness_kpa"] - base["stiffness_kpa"])
                if base["stiffness_kpa"] is not None and post["stiffness_kpa"] is not None else None,
            "base_pdff": base["PDFF_pct"], "post_pdff": post["PDFF_pct"],
            "base_stiff": base["stiffness_kpa"], "post_stiff": post["stiffness_kpa"],
        })
    return out


def classify_response(change):
    """steatosis 가역형 vs fibrosis 가역형 vs 무반응형 분류."""
    dpd = change.get("dPDFF")
    dst = change.get("dStiffness")
    base_pd = change.get("base_pdff") or 0.0
    base_st = change.get("base_stiff") or 0.0
    steat_rev = dpd is not None and base_pd > 0 and (dpd / base_pd) <= -0.20
    fib_rev = dst is not None and base_st > 0 and (dst / base_st) <= -0.10
    if steat_rev and fib_rev:
        return "steatosis+fibrosis 가역형"
    if steat_rev:
        return "steatosis 가역형"
    if fib_rev:
        return "fibrosis 가역형"
    return "무반응형"


def histology_correlation(table, histo_rows):
    """terminal histology vs 영상 readout 상관 (post 시점 매칭)."""
    histo_map = {r["animal_id"]: r for r in histo_rows}
    img_pdff, img_stiff, img_mrs = [], [], []
    nas, sirius, cpa = [], [], []
    matched = 0
    for rec in table:
        if rec["timepoint"] != "post":
            continue
        h = histo_map.get(rec["animal_id"])
        if not h:
            continue
        matched += 1
        img_pdff.append(rec["PDFF_pct"])
        img_stiff.append(rec["stiffness_kpa"])
        img_mrs.append(rec["mrs_pdff_pct"])
        nas.append(float(h["nas_score"]))
        sirius.append(float(h["sirius_red_pct"]))
        cpa.append(float(h["cpa_pct"]))
    return {
        "matched": matched,
        "PDFF_vs_NAS_pearson": pearson(img_pdff, nas),
        "PDFF_vs_NAS_spearman": spearman(img_pdff, nas),
        "MRS_vs_NAS_pearson": pearson(img_mrs, nas),
        "stiffness_vs_SiriusRed_pearson": pearson(img_stiff, sirius),
        "stiffness_vs_SiriusRed_spearman": spearman(img_stiff, sirius),
        "stiffness_vs_CPA_pearson": pearson(img_stiff, cpa),
    }


def cohort_stats(table):
    """군 × readout × 시점 요약 + paired t-test."""
    groups = sorted(set(r["group"] for r in table))
    tps = ["baseline", "post"]
    summary = {}
    for g in groups:
        summary[g] = {}
        for tp in tps:
            sub = [r for r in table if r["group"] == g and r["timepoint"] == tp]
            summary[g][tp] = {
                "n": len(sub),
                "PDFF_mean": mean([r["PDFF_pct"] for r in sub]),
                "PDFF_sd": stdev([r["PDFF_pct"] for r in sub]),
                "MRS_mean": mean([r["mrs_pdff_pct"] for r in sub]),
                "stiffness_mean": mean([r["stiffness_kpa"] for r in sub]),
                "stiffness_sd": stdev([r["stiffness_kpa"] for r in sub]),
            }
    # paired t-test per group (baseline vs post)
    tests = {}
    by_an = {}
    for r in table:
        by_an.setdefault((r["group"], r["animal_id"]), {})[r["timepoint"]] = r
    for g in groups:
        pre_pd, post_pd, pre_st, post_st = [], [], [], []
        for (gg, an), tpd in by_an.items():
            if gg != g:
                continue
            b, p = tpd.get("baseline"), tpd.get("post")
            if b and p:
                pre_pd.append(b["PDFF_pct"]); post_pd.append(p["PDFF_pct"])
                pre_st.append(b["stiffness_kpa"]); post_st.append(p["stiffness_kpa"])
        tests[g] = {
            "PDFF_ttest": paired_ttest(pre_pd, post_pd),
            "stiffness_ttest": paired_ttest(pre_st, post_st),
        }
    return summary, tests


# --------------------------------------------------------------------------
# 출력
# --------------------------------------------------------------------------
def fmt(v, nd=2):
    if v is None:
        return "NA"
    if isinstance(v, float) and (math.isnan(v) or math.isinf(v)):
        return "NA" if math.isnan(v) else ">99"
    return f"{v:.{nd}f}"


def run_demo(lang="ko"):
    me = load_csv(os.path.join(DATA_DIR, "sample_multiecho.csv"))
    mrs = load_csv(os.path.join(DATA_DIR, "sample_mrs_mre.csv"))
    histo = load_csv(os.path.join(DATA_DIR, "sample_histology.csv"))

    print("=" * 72)
    print(" HepatoMRQuant — rodent in vivo 간 MRI 정량 (오프라인 데모)")
    print(" MASLD 도메인 / 동물실험 도구 (in vivo 이미지 정량)")
    print("=" * 72)
    print("[DISCLAIMER] " + DISCLAIMER)
    print()
    print(f"[INGEST] multi-echo {len(me)} rows, MRS/MRE {len(mrs)} rows, "
          f"histology {len(histo)} rows")
    print()

    table = build_quant_table(me, mrs)

    # 1. 통합 readout 테이블 (T2* 보정 PDFF / MRS / stiffness)
    print("-" * 72)
    print("[1] 통합 readout 테이블 (T2* 보정 PDFF% / MRS lipid% / stiffness kPa)")
    print("-" * 72)
    hdr = (f"{'animal':6} {'grp':8} {'tp':9} {'T2*ms':>7} {'PDFF%':>7} "
           f"{'PDFFunc':>8} {'MRS%':>6} {'kPa':>6} flags")
    print(hdr)
    for r in table:
        print(f"{r['animal_id']:6} {r['group']:8} {r['timepoint']:9} "
              f"{fmt(r['T2star_ms']):>7} {fmt(r['PDFF_pct']):>7} "
              f"{fmt(r['PDFF_uncorrected_pct']):>8} {fmt(r['mrs_pdff_pct']):>6} "
              f"{fmt(r['stiffness_kpa']):>6} {','.join(r['flags']) or '-'}")
    print()

    # 2. 종단 within-animal 변화 + 반응 분류
    print("-" * 72)
    print("[2] 종단 within-animal 변화 (baseline→post) + 반응 분류")
    print("-" * 72)
    changes = longitudinal_changes(table)
    print(f"{'animal':6} {'grp':8} {'dPDFF':>7} {'dMRS':>7} {'dStiff':>7}  class")
    for c in changes:
        print(f"{c['animal_id']:6} {c['group']:8} "
              f"{fmt(c['dPDFF']):>7} {fmt(c['dMRS']):>7} {fmt(c['dStiffness']):>7}"
              f"  {classify_response(c)}")
    print()

    # 3. 모델 참고 범위
    print("-" * 72)
    print("[3] 동물 모델별 전형 PDFF/stiffness 참고 범위 (합성 가정)")
    print("-" * 72)
    used = sorted(set(r["group"] for r in table))
    for g in used:
        ref = MODEL_REFERENCE.get(g)
        if ref:
            print(f"  {g:9} PDFF {ref['pdff_pct'][0]}-{ref['pdff_pct'][1]}% | "
                  f"stiffness {ref['stiffness_kpa'][0]}-{ref['stiffness_kpa'][1]} kPa "
                  f"| {ref['note']}")
    print()

    # 4. histology 상관
    print("-" * 72)
    print("[4] 영상-조직(terminal histology) 상관 (post 시점 매칭)")
    print("-" * 72)
    corr = histology_correlation(table, histo)
    print(f"  matched animals: {corr['matched']}")
    print(f"  PDFF vs NAS      : Pearson {fmt(corr['PDFF_vs_NAS_pearson'],3)} | "
          f"Spearman {fmt(corr['PDFF_vs_NAS_spearman'],3)}")
    print(f"  MRS%  vs NAS     : Pearson {fmt(corr['MRS_vs_NAS_pearson'],3)}")
    print(f"  stiffness vs SiriusRed%: Pearson {fmt(corr['stiffness_vs_SiriusRed_pearson'],3)} | "
          f"Spearman {fmt(corr['stiffness_vs_SiriusRed_spearman'],3)}")
    print(f"  stiffness vs CPA%      : Pearson {fmt(corr['stiffness_vs_CPA_pearson'],3)}")
    print()

    # 5. cohort 통계
    print("-" * 72)
    print("[5] cohort 통계 (군 × 시점 요약 + paired t-test)")
    print("-" * 72)
    summary, tests = cohort_stats(table)
    for g in sorted(summary):
        print(f"  [{g}]")
        for tp in ("baseline", "post"):
            s = summary[g][tp]
            print(f"    {tp:9} n={s['n']} | PDFF {fmt(s['PDFF_mean'])}±{fmt(s['PDFF_sd'])}% "
                  f"| MRS {fmt(s['MRS_mean'])}% | stiff {fmt(s['stiffness_mean'])}±"
                  f"{fmt(s['stiffness_sd'])} kPa")
        pt = tests[g]["PDFF_ttest"]
        st = tests[g]["stiffness_ttest"]
        if pt[0] is not None:
            print(f"    paired t-test PDFF: t={fmt(pt[0],3)} df={pt[1]} p={fmt(pt[2],4)}")
        if st[0] is not None:
            print(f"    paired t-test stiffness: t={fmt(st[0],3)} df={st[1]} p={fmt(st[2],4)}")
    print()

    # 텍스트 리포트
    print("=" * 72)
    print("[REPORT] 방법·결과 요약")
    print("=" * 72)
    _print_report(table, changes, corr, summary, tests, lang)
    print()
    print("[NOTE] 출처 맥락: AASLD MASLD 2023 / EASL-EASD-EASO MASLD 2024.")
    print("[NOTE] 모든 데이터는 합성. 정량 결과는 사용자 검증 필요.")


def _print_report(table, changes, corr, summary, tests, lang):
    cdahfd = summary.get("CDAHFD", {})
    base = cdahfd.get("baseline", {})
    post = cdahfd.get("post", {})
    n_steat = sum(1 for c in changes if "steatosis" in classify_response(c))
    n_fib = sum(1 for c in changes if "fibrosis" in classify_response(c))
    n_none = sum(1 for c in changes if classify_response(c) == "무반응형")

    if lang == "en":
        print("Methods: T2*-corrected PDFF was computed from multi-echo magnitude "
              "signal via log-linear fit S(TE)=S0*exp(-TE/T2*) (stdlib regression); "
              "MRS PDFF = lipid/(lipid+water); stiffness via mu=rho*c^2 (rho=1000).")
        print(f"Results: CDAHFD PDFF {fmt(base.get('PDFF_mean'))}% -> "
              f"{fmt(post.get('PDFF_mean'))}% after treatment; "
              f"stiffness {fmt(base.get('stiffness_mean'))} -> "
              f"{fmt(post.get('stiffness_mean'))} kPa.")
        print(f"Response classes: steatosis-reversible n={n_steat}, "
              f"fibrosis-reversible n={n_fib}, non-responder n={n_none}.")
        print(f"Imaging-histology: PDFF vs NAS Pearson "
              f"{fmt(corr['PDFF_vs_NAS_pearson'],3)}; stiffness vs Sirius Red "
              f"{fmt(corr['stiffness_vs_SiriusRed_pearson'],3)}.")
    else:
        print("[방법] multi-echo magnitude 신호에서 log-linear 적합 "
              "S(TE)=S0*exp(-TE/T2*)(표준라이브러리 회귀)로 T2* 보정 PDFF 산출; "
              "MRS PDFF = lipid/(lipid+water); stiffness μ=ρ·c²(ρ=1000) 변환.")
        print(f"[결과] CDAHFD군 PDFF {fmt(base.get('PDFF_mean'))}% → "
              f"치료 후 {fmt(post.get('PDFF_mean'))}%, stiffness "
              f"{fmt(base.get('stiffness_mean'))} → {fmt(post.get('stiffness_mean'))} kPa.")
        print(f"[반응분류] steatosis 가역형 n={n_steat}, fibrosis 가역형 n={n_fib}, "
              f"무반응형 n={n_none}.")
        print(f"[영상-조직] PDFF vs NAS Pearson "
              f"{fmt(corr['PDFF_vs_NAS_pearson'],3)}, stiffness vs Sirius Red "
              f"{fmt(corr['stiffness_vs_SiriusRed_pearson'],3)}.")


def main(argv=None):
    p = argparse.ArgumentParser(
        prog="main.py",
        description="HepatoMRQuant — rodent in vivo 간 MRI 정량 (오프라인 CLI 데모, "
                    "표준 라이브러리만 사용). 본 도구는 연구용·참고용입니다.",
    )
    p.add_argument("--demo", action="store_true",
                   help="data/ 의 합성 샘플로 전체 정량 파이프라인 실행")
    p.add_argument("--lang", choices=["ko", "en"], default="ko",
                   help="리포트 언어 (기본 ko)")
    args = p.parse_args(argv)

    if args.demo:
        run_demo(lang=args.lang)
    else:
        p.print_help()
    return 0


if __name__ == "__main__":
    sys.exit(main())