"""
core.py — MASLDNITSchedule-Kor 핵심 계산 모듈

MASLD/MASH 종단 관찰연구의 비침습검사(NIT) 측정 스케줄 설계를 위한 순수 계산 함수.
의존성: numpy, scipy 만 사용 (statsmodels 등 금지).

통계 공식 출처 주석:
 - MDC(최소검출변화, Minimal Detectable Change):
     MDC = z * sqrt(2) * SEM,  여기서 SEM = 측정-재측정 SD (한 측정의 표준오차)
     일반적 형태 MDC95 = 1.96 * sqrt(2) * SEM
     (출처: 측정학 표준 공식. de Vet 등 'Measurement in Medicine'; Weir 2005 J Strength Cond Res)
 - 평균(반복측정 평균)의 MDC: n회 반복 시 평균의 SE = SEM/sqrt(n)
     MDC_mean = 1.96 * sqrt(2) * SEM / sqrt(n)
 - 종단 기울기(slope) 검정력: 선형혼합모형 기울기 분산 근사
     Var(slope) ~ sigma_e^2 / (n_subj * Sxx),  Sxx = sum((t_i - tbar)^2)
     (출처: 종단자료 표본수 공식. Diggle 등 'Analysis of Longitudinal Data';
      Fitzmaurice 등 'Applied Longitudinal Analysis')

면책: 본 도구는 참고용·연구용이며, 실제 연구 설계는 생물통계학자 검토가 필수다.
NIT 변동성 값은 문헌 기반 추정치다.
"""

import json
import os
import math

try:
    import numpy as np
except ImportError as e:  # pragma: no cover
    raise ImportError(
        "numpy가 필요합니다. README의 설치 안내를 참고하세요. (pip install --user numpy)"
    ) from e

try:
    from scipy import stats as _scipy_stats
    _HAS_SCIPY = True
except ImportError:  # pragma: no cover
    _HAS_SCIPY = False


DISCLAIMER = (
    "[면책] 본 도구(MASLDNITSchedule-Kor)는 참고용·연구용 설계 보조 도구입니다. "
    "실제 관찰연구 설계는 생물통계학자 검토가 필수이며, NIT 변동성·진행률·단가 값은 "
    "발표 문헌 기반 추정치입니다. 정밀 값이 아닙니다."
)

DATA_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "nit_variability.json")


# ---------------------------------------------------------------------------
# 정규분포 보조 (scipy 없을 때 대비)
# ---------------------------------------------------------------------------
def _z_from_alpha(alpha_two_sided=0.05):
    """양측 검정 임계 z 값."""
    if _HAS_SCIPY:
        return float(_scipy_stats.norm.ppf(1.0 - alpha_two_sided / 2.0))
    # 근사: 흔한 값만
    table = {0.05: 1.959963985, 0.10: 1.644853627, 0.01: 2.575829304}
    return table.get(round(alpha_two_sided, 3), 1.959963985)


def _z_from_power(power=0.80):
    """검정력에 대응하는 z 값 (단측)."""
    if _HAS_SCIPY:
        return float(_scipy_stats.norm.ppf(power))
    table = {0.80: 0.841621234, 0.90: 1.281551566, 0.95: 1.644853627}
    return table.get(round(power, 2), 0.841621234)


def _norm_cdf(x):
    if _HAS_SCIPY:
        return float(_scipy_stats.norm.cdf(x))
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))


# ---------------------------------------------------------------------------
# 데이터 로드
# ---------------------------------------------------------------------------
def load_nit_library(path=None):
    """NIT 변동성 라이브러리 JSON 로드."""
    path = path or DATA_PATH
    with open(path, "r", encoding="utf-8") as f:
        return json.load(f)


def list_nits(library=None):
    """사용 가능한 NIT 키 목록."""
    library = library or load_nit_library()
    return list(library["nits"].keys())


def get_nit(library, nit_key):
    """NIT 설정 반환 (키 검증 포함)."""
    if nit_key not in library["nits"]:
        raise KeyError(f"알 수 없는 NIT: {nit_key}. 가능: {list(library['nits'].keys())}")
    return library["nits"][nit_key]


def resolve_sem(nit_cfg):
    """
    NIT 설정에서 측정-재측정 표준오차(SEM, 절대단위)를 산출.
    test_retest_sd_abs 우선, 없으면 CV(%) x typical_baseline_value 로 환산.
    반환: (sem_abs, unit, baseline_value)
    """
    v = nit_cfg["variability"]
    baseline = v.get("typical_baseline_value")
    sd_abs = v.get("test_retest_sd_abs")
    cv = v.get("test_retest_cv_percent")
    if sd_abs is not None:
        sem = float(sd_abs)
    elif cv is not None and baseline is not None:
        sem = float(cv) / 100.0 * float(baseline)
    else:
        raise ValueError("변동성 정보 부족: SD 또는 (CV + baseline) 필요")
    return sem, nit_cfg.get("unit", ""), baseline


# ---------------------------------------------------------------------------
# 기능 1 보조: 스케줄 빌더
# ---------------------------------------------------------------------------
def build_schedule(baseline_month=0, n_followups=2, interval_months=12):
    """
    측정 visit 시점(개월) 리스트 생성.
    baseline + n_followups 개의 추적 visit, 각 interval_months 간격.
    반환: dict { 'visit_months': [...], 'total_months': , 'n_visits': }
    """
    if n_followups < 1:
        raise ValueError("추적 visit는 최소 1회 필요")
    if interval_months <= 0:
        raise ValueError("visit 간격은 양수")
    months = [float(baseline_month)]
    for i in range(1, n_followups + 1):
        months.append(float(baseline_month + i * interval_months))
    return {
        "visit_months": months,
        "total_months": months[-1] - months[0],
        "n_visits": len(months),
    }


# ---------------------------------------------------------------------------
# 기능 2: MDC (최소검출변화) 산출
# ---------------------------------------------------------------------------
def compute_mdc(sem, n_repeat=1, confidence=0.95):
    """
    두 시점 간 변화의 최소검출변화(MDC).
      MDC = z * sqrt(2) * SE_mean,  SE_mean = SEM / sqrt(n_repeat)
    n_repeat: 각 visit에서의 반복 측정 횟수(예: VCTE를 같은 날 2회).
    (출처: 측정학 표준 MDC 공식, de Vet 'Measurement in Medicine')
    """
    if sem <= 0:
        raise ValueError("SEM은 양수")
    if n_repeat < 1:
        raise ValueError("반복 횟수는 1 이상")
    alpha = 1.0 - confidence
    z = _z_from_alpha(alpha)
    se_mean = sem / math.sqrt(n_repeat)
    mdc = z * math.sqrt(2.0) * se_mean
    return {
        "mdc": mdc,
        "sem": sem,
        "se_mean": se_mean,
        "n_repeat": n_repeat,
        "confidence": confidence,
        "z": z,
    }


def mdc_table_for_nit(library, nit_key, repeats=(1, 2, 3)):
    """NIT별, 반복측정 횟수별 MDC 표."""
    nit = get_nit(library, nit_key)
    sem, unit, baseline = resolve_sem(nit)
    rows = []
    for n in repeats:
        r = compute_mdc(sem, n_repeat=n)
        rows.append({
            "n_repeat": n,
            "mdc": r["mdc"],
            "unit": unit,
        })
    return {
        "nit": nit_key,
        "label": nit.get("label", nit_key),
        "unit": unit,
        "sem": sem,
        "baseline": baseline,
        "rows": rows,
    }


def mdc_warning(mdc, meaningful_change, unit=""):
    """
    MDC가 임상적으로 의미 있는 변화보다 크면 경고 문구 생성.
    반환: (is_warning: bool, message: str)
    """
    if meaningful_change is None:
        return False, "임상적 의미 변화 기준 미설정"
    if mdc > meaningful_change:
        return True, (
            f"경고: 이 스케줄의 MDC({mdc:.2f} {unit})가 임상적 의미 변화 기준"
            f"({meaningful_change:.2f} {unit})보다 큽니다. "
            f"{meaningful_change:.2f} {unit} 미만 변화는 측정오차와 구분 불가합니다."
        )
    return False, (
        f"양호: MDC({mdc:.2f} {unit}) <= 의미 변화 기준({meaningful_change:.2f} {unit}). "
        f"의미 있는 변화를 검출할 수 있습니다."
    )


# ---------------------------------------------------------------------------
# 기능 3: 연구기간 vs 진행률 tradeoff
# ---------------------------------------------------------------------------
def time_to_detectable_change(mdc, progression_per_year):
    """
    자연 진행률(연간 변화량) 기준, MDC만큼의 변화에 도달하는 데 걸리는 개월 수.
      months = MDC / |progression_per_year| * 12
    반환: 개월 수 (float). 진행률 0이면 inf.
    """
    rate = abs(progression_per_year)
    if rate <= 0:
        return float("inf")
    years = mdc / rate
    return years * 12.0


def tradeoff_curve(library, nit_key, subgroup=None, n_repeat=1,
                   max_months=60, step_months=3):
    """
    연구기간(개월) 대비 검출 가능 여부 곡선.
    각 시점에서 누적 기대 변화 vs MDC 비교.
    반환: dict { months: [...], expected_change: [...], mdc: , detectable_from_month: }
    """
    nit = get_nit(library, nit_key)
    sem, unit, baseline = resolve_sem(nit)
    mdc = compute_mdc(sem, n_repeat=n_repeat)["mdc"]

    prog = nit["progression"]
    rate = prog["estimate_per_year"]
    if subgroup and subgroup in prog.get("subgroups", {}):
        rate = prog["subgroups"][subgroup]

    months = list(range(0, max_months + 1, step_months))
    expected = [abs(rate) * (m / 12.0) for m in months]
    detect_month = None
    for m, e in zip(months, expected):
        if e >= mdc:
            detect_month = m
            break

    return {
        "nit": nit_key,
        "unit": unit,
        "months": months,
        "expected_change": expected,
        "mdc": mdc,
        "progression_per_year": rate,
        "subgroup": subgroup,
        "detectable_from_month": detect_month,
        "time_to_mdc_months": time_to_detectable_change(mdc, rate),
    }


# ---------------------------------------------------------------------------
# 기능 4: visit window 최적화 + 결측 영향 Monte-Carlo
# ---------------------------------------------------------------------------
def slope_power_analytic(visit_months, n_subjects, sem,
                         true_slope_per_year, alpha=0.05,
                         subject_slope_cv=0.6):
    """
    선형혼합모형(임의기울기) 기울기 검정의 해석적(근사) 검정력.

    피험자별 기울기 추정량 분산은 두 성분의 합:
      Var(slope_i) = tau^2 + sigma_e^2 / Sxx
        - tau^2     : 피험자 간 진짜 기울기 분산 (집단 내 이질성)
        - sigma_e^2 : 측정오차 분산 (= sem^2)
        - Sxx       : sum((t_i - tbar)^2),  t 단위는 '년'
    집단 평균 기울기의 SE = sqrt(Var(slope_i) / n_subj).
      검정통계량 z = mean_slope / SE,  비중심모수 ncp = |slope| / SE
      power = P(|Z| > z_alpha/2 | ncp)
    (출처: 종단자료 검정력 근사. Fitzmaurice 'Applied Longitudinal Analysis';
     Diggle 'Analysis of Longitudinal Data')

    subject_slope_cv: 피험자 간 기울기 표준편차를 |true_slope| 대비 비율로 가정
      (기본 0.6 — 진행률은 개인차가 크다는 보수적 가정. 0이면 측정오차만 반영).
    """
    t_years = np.array(visit_months, dtype=float) / 12.0
    tbar = t_years.mean()
    sxx = float(np.sum((t_years - tbar) ** 2))
    if sxx <= 0 or n_subjects < 2:
        return 0.0
    tau2 = (subject_slope_cv * true_slope_per_year) ** 2
    var_slope_i = tau2 + (sem ** 2) / sxx
    se_mean = math.sqrt(var_slope_i / n_subjects)
    if se_mean <= 0:
        return 1.0
    ncp = abs(true_slope_per_year) / se_mean  # 비중심 모수
    z_alpha = _z_from_alpha(alpha)
    power = (1.0 - _norm_cdf(z_alpha - ncp)) + _norm_cdf(-z_alpha - ncp)
    return float(min(max(power, 0.0), 1.0))


def _ideal_detect_rate(t_plan_years, sem, tau, true_slope, mdc):
    """
    결측·window 지터 없는 이상적 조건에서 피험자 단위 변화 검출률(해석적).
    baseline~종료 관찰변화 Δ = slope_i*span + (e_end - e_base)
      slope_i ~ N(slope, tau^2)  →  Δ ~ N(slope*span, tau^2*span^2 + 2*sem^2)
    P(|Δ| > MDC) 를 정규분포로 계산.
    """
    span = float(t_plan_years[-1] - t_plan_years[0])
    mean_delta = true_slope * span
    var_delta = (tau ** 2) * (span ** 2) + 2.0 * (sem ** 2)
    if var_delta <= 0:
        return 1.0 if abs(mean_delta) > mdc else 0.0
    sd = math.sqrt(var_delta)
    p_above = 1.0 - _norm_cdf((mdc - mean_delta) / sd)
    p_below = _norm_cdf((-mdc - mean_delta) / sd)
    return float(min(max(p_above + p_below, 0.0), 1.0))


def simulate_missingness(visit_months, n_subjects, sem, true_slope_per_year,
                         window_weeks=4, missing_rate=0.15,
                         n_sim=400, alpha=0.05, seed=42,
                         subject_slope_cv=0.6):
    """
    visit window 폭과 결측률을 반영한 종단 기울기 검정력 Monte-Carlo 시뮬레이션.

    모형: 각 피험자 i, visit j 에서
        y_ij = b0_i + slope_i * t_ij + e_ij
      - b0_i    ~ N(0, sem^2)                    (피험자 임의절편)
      - slope_i ~ N(true_slope, tau^2)           (피험자 임의기울기; 진행률 개인차)
                  tau = subject_slope_cv * |true_slope|
      - e_ij    ~ N(0, sem^2)                    (측정오차)
      - t_ij = 계획시점 + 균등(-window, +window) 지터  (window 단위: 년)
      - 각 (i,j) 측정은 확률 missing_rate 로 결측
    각 시뮬에서 피험자별 OLS 기울기를 구해 집단 평균 기울기의 1-표본 검정
    → 유의비율 = 검정력.

    subject_slope_cv: 피험자 간 진행률 표준편차를 |true_slope| 대비 비율로 가정
      (기본 0.6 — 진행률 개인차가 크다는 보수적 가정).

    반환: dict { power, ... }
    """
    rng = np.random.default_rng(seed)
    t_plan = np.array(visit_months, dtype=float) / 12.0  # 년
    window_years = (window_weeks / 52.0)
    n_visit = len(t_plan)
    z_alpha = _z_from_alpha(alpha)
    tau = abs(subject_slope_cv * true_slope_per_year)
    # MDC 기준선 (반복 1회 가정) — 피험자 단위 변화 분류용
    mdc = compute_mdc(sem, n_repeat=1)["mdc"]

    sig_count = 0
    valid_sim = 0
    # 피험자 단위 지표 누적
    classify_hits = 0          # 관찰 변화 |Δ| > MDC 인 피험자 수
    classify_total = 0         # baseline+종료 둘 다 측정된 피험자 수
    analyzable_total = 0       # 기울기 추정 가능(2회 이상 측정) 피험자 수
    analyzable_planned = 0     # 계획상 전체 피험자 수

    for _ in range(n_sim):
        slopes = []
        for _i in range(n_subjects):
            analyzable_planned += 1
            b0 = rng.normal(0.0, sem)
            slope_i = rng.normal(true_slope_per_year, tau) if tau > 0 else true_slope_per_year
            # visit window 지터
            jitter = rng.uniform(-window_years, window_years, size=n_visit)
            t_obs = t_plan + jitter
            # 결측 마스크
            keep = rng.random(n_visit) >= missing_rate
            true_y = b0 + slope_i * t_obs
            obs_y = true_y + rng.normal(0.0, sem, size=n_visit)
            # --- 피험자 단위 변화 분류 (baseline=visit0, 종료=마지막 visit) ---
            if keep[0] and keep[-1]:
                classify_total += 1
                delta = obs_y[-1] - obs_y[0]
                if abs(delta) > mdc:
                    classify_hits += 1
            # --- 기울기 추정 ---
            if keep.sum() < 2:
                continue  # 기울기 추정 불가
            analyzable_total += 1
            tt = t_obs[keep]
            y = obs_y[keep]
            tbar = tt.mean()
            sxx = np.sum((tt - tbar) ** 2)
            if sxx <= 0:
                continue
            beta = np.sum((tt - tbar) * (y - y.mean())) / sxx
            slopes.append(beta)
        if len(slopes) < 2:
            continue
        slopes = np.array(slopes)
        # 집단 평균 기울기에 대한 1-표본 z 검정 (피험자 간 분산 사용)
        mean_slope = slopes.mean()
        sd_slope = slopes.std(ddof=1)
        se_mean = sd_slope / math.sqrt(len(slopes))
        if se_mean <= 0:
            continue
        zstat = mean_slope / se_mean
        valid_sim += 1
        if abs(zstat) > z_alpha:
            sig_count += 1

    power = (sig_count / valid_sim) if valid_sim > 0 else 0.0
    # 해석적 기준선(결측·window 없음)
    analytic = slope_power_analytic(visit_months, n_subjects, sem,
                                    true_slope_per_year, alpha=alpha)
    # 피험자 단위 변화 검출률 (스케줄 품질에 민감한 핵심 지표)
    detect_rate = (classify_hits / classify_total) if classify_total > 0 else 0.0
    endpoint_completeness = (classify_total / analyzable_planned) if analyzable_planned > 0 else 0.0
    analyzable_frac = (analyzable_total / analyzable_planned) if analyzable_planned > 0 else 0.0
    # 이상적(결측 0, window 0) 변화 검출률 — 동일 모형 무결측 근사
    ideal_detect = _ideal_detect_rate(t_plan, sem, tau, true_slope_per_year, mdc)
    return {
        "power_simulated": power,
        "power_analytic_ideal": analytic,
        "power_loss": max(analytic - power, 0.0),
        "subject_detect_rate": detect_rate,
        "subject_detect_ideal": ideal_detect,
        "detect_rate_loss": max(ideal_detect - detect_rate, 0.0),
        "endpoint_completeness": endpoint_completeness,
        "analyzable_fraction": analyzable_frac,
        "mdc": mdc,
        "n_sim_valid": valid_sim,
        "window_weeks": window_weeks,
        "missing_rate": missing_rate,
        "n_subjects": n_subjects,
    }


def optimize_visit_window(visit_months, n_subjects, sem, true_slope_per_year,
                          window_candidates=(2, 4, 8, 12),
                          missing_rate=0.15, n_sim=300, alpha=0.05, seed=42):
    """
    여러 visit window 폭 후보에 대해 검정력·피험자 변화검출률을 비교, 권장값 산출.

    좁은 window는 측정 시점이 계획에 가까워 baseline~종료 간격(span)이 일정하지만
    실무상 결측 위험이 커지고, 넓은 window는 결측 여유가 크나 span 변동이 커져
    피험자 단위 변화 검출이 흔들린다는 tradeoff를 Monte-Carlo로 평가.
    여기서는 동일 missing_rate 가정 하에서 window 폭 자체의 영향을 비교한다.

    비교 기준은 '피험자 단위 변화 검출률'(subject_detect_rate) — 집단 평균
    기울기 검정력은 표본이 크면 포화되어 스케줄 차이에 둔감하기 때문이다.
    반환: dict { results: [...], recommended_window_weeks: }
    """
    results = []
    for w in window_candidates:
        sim = simulate_missingness(
            visit_months, n_subjects, sem, true_slope_per_year,
            window_weeks=w, missing_rate=missing_rate,
            n_sim=n_sim, alpha=alpha, seed=seed,
        )
        results.append({
            "window_weeks": w,
            "power": sim["power_simulated"],
            "power_loss": sim["power_loss"],
            "detect_rate": sim["subject_detect_rate"],
            "detect_rate_loss": sim["detect_rate_loss"],
            "endpoint_completeness": sim["endpoint_completeness"],
        })
    # 피험자 변화검출률 최대 window 추천 (동률이면 더 좁은 window = span 일관성↑)
    best = max(results, key=lambda r: (round(r["detect_rate"], 4), -r["window_weeks"]))
    return {
        "results": results,
        "recommended_window_weeks": best["window_weeks"],
        "recommended_detect_rate": best["detect_rate"],
        "recommended_power": best["power"],
    }


# ---------------------------------------------------------------------------
# 기능 5: NIT 조합·비용 비교 (비용효율 frontier)
# ---------------------------------------------------------------------------
def cost_efficiency_frontier(library, nit_keys, n_subjects, visit_months,
                             true_slope_map=None, n_repeat=1, alpha=0.05):
    """
    여러 NIT 각각에 대해 1인당 총 검사비와 종단 검정력을 계산, 비용효율 frontier 산출.

    1인당 비용 = 단가 x visit 수 x 반복횟수.
    검정력은 slope_power_analytic 사용. 진행률은 해당 NIT subgroup 기본값 사용
    (true_slope_map 으로 NIT별 override 가능).

    frontier: 같은 검정력 대비 비용이 가장 낮은 NIT(또는 그 역) 집합.
    반환: dict { rows: [...], frontier_keys: [...] }
    """
    n_visit = len(visit_months)
    rows = []
    for key in nit_keys:
        nit = get_nit(library, key)
        sem, unit, baseline = resolve_sem(nit)
        # 진행률
        prog = nit["progression"]
        slope = prog["estimate_per_year"]
        sg = prog.get("subgroups", {})
        if "T2DM_or_MASH" in sg:
            slope = sg["T2DM_or_MASH"]
        if true_slope_map and key in true_slope_map:
            slope = true_slope_map[key]
        # 반복측정 시 유효 SEM 감소
        eff_sem = sem / math.sqrt(n_repeat)
        power = slope_power_analytic(visit_months, n_subjects, eff_sem,
                                     slope, alpha=alpha)
        cost_per_visit = nit.get("cost_krw", 0)
        cost_per_subject = cost_per_visit * n_visit * n_repeat
        rows.append({
            "nit": key,
            "label": nit.get("label", key),
            "category": nit.get("category", ""),
            "unit": unit,
            "cost_per_subject": cost_per_subject,
            "cost_total": cost_per_subject * n_subjects,
            "power": power,
            "slope_per_year": slope,
            "cost_per_power": (cost_per_subject / power) if power > 0 else float("inf"),
        })

    # frontier: 비용 오름차순 정렬 후, 검정력이 단조 증가하는 점만 유지
    by_cost = sorted(rows, key=lambda r: r["cost_per_subject"])
    frontier = []
    best_power = -1.0
    for r in by_cost:
        if r["power"] > best_power:
            frontier.append(r["nit"])
            best_power = r["power"]
    return {
        "rows": rows,
        "frontier_keys": frontier,
        "n_subjects": n_subjects,
        "n_visits": n_visit,
        "n_repeat": n_repeat,
    }


# ---------------------------------------------------------------------------
# 기능 6: 설계 리포트 텍스트 생성
# ---------------------------------------------------------------------------
def generate_design_report(library, nit_key, n_subjects=120, n_followups=2,
                           interval_months=12, subgroup="T2DM_or_MASH",
                           n_repeat=1, window_weeks=4, missing_rate=0.15,
                           budget_krw=None, n_sim=300):
    """
    권장 스케줄·MDC·기간·예산을 관찰연구 프로토콜 측정계획 섹션(국문) 텍스트로 생성.
    반환: 문자열(리포트).
    """
    nit = get_nit(library, nit_key)
    sem, unit, baseline = resolve_sem(nit)
    sched = build_schedule(0, n_followups, interval_months)
    mdc = compute_mdc(sem, n_repeat=n_repeat)["mdc"]

    # 진행률
    prog = nit["progression"]
    slope = prog["estimate_per_year"]
    if subgroup and subgroup in prog.get("subgroups", {}):
        slope = prog["subgroups"][subgroup]

    tradeoff = tradeoff_curve(library, nit_key, subgroup=subgroup,
                              n_repeat=n_repeat, max_months=120, step_months=6)
    sim = simulate_missingness(sched["visit_months"], n_subjects, sem, slope,
                               window_weeks=window_weeks,
                               missing_rate=missing_rate, n_sim=n_sim)

    # 비용
    cost_per_visit = nit.get("cost_krw", 0)
    cost_per_subject = cost_per_visit * sched["n_visits"] * n_repeat
    cost_total = cost_per_subject * n_subjects

    # 의미 변화 기준
    rt = library.get("responder_thresholds", {})
    meaningful_map = {
        "VCTE_LSM": rt.get("VCTE_LSM_meaningful_kPa"),
        "MRE": rt.get("MRE_meaningful_kPa"),
        "FIB-4": rt.get("FIB-4_meaningful_index"),
        "ELF": rt.get("ELF_meaningful_score"),
        "MRI_PDFF": None,
    }
    meaningful = meaningful_map.get(nit_key)
    is_warn, warn_msg = mdc_warning(mdc, meaningful, unit)

    lines = []
    lines.append("=" * 72)
    lines.append("관찰연구 프로토콜 — 측정계획 섹션 (NIT 종단 측정 스케줄)")
    lines.append("MASLD/MASH 비침습검사 설계: MASLDNITSchedule-Kor 자동 생성")
    lines.append("=" * 72)
    lines.append("")
    lines.append(f"[1] 대상 비침습검사(NIT): {nit.get('label', nit_key)} ({nit_key})")
    lines.append(f"    - 범주: {nit.get('category', '-')}")
    lines.append(f"    - 측정 단위: {unit}")
    lines.append(f"    - 측정-재측정 표준오차(SEM, 추정): {sem:.3f} {unit}")
    qc = nit.get("quality_criteria", {})
    if qc.get("rule"):
        lines.append(f"    - 검사 품질기준: {qc['rule']}")
    lines.append("")
    lines.append("[2] 측정 visit 스케줄")
    visit_str = ", ".join(f"{m:.0f}개월" for m in sched["visit_months"])
    lines.append(f"    - 총 visit 수: {sched['n_visits']}회 (baseline 1 + 추적 {n_followups})")
    lines.append(f"    - visit 시점: {visit_str}")
    lines.append(f"    - 총 추적기간: {sched['total_months']:.0f}개월")
    lines.append(f"    - visit window: 계획 시점 +-{window_weeks}주 이내 측정 권장")
    lines.append(f"    - 각 visit당 반복측정: {n_repeat}회 (평균값 분석 사용)")
    lines.append("")
    lines.append("[3] 최소검출변화(MDC) 및 검출 가능성")
    lines.append(f"    - 95% MDC (반복 {n_repeat}회 평균 기준): {mdc:.3f} {unit}")
    lines.append(f"      공식: MDC = 1.96 x sqrt(2) x SEM / sqrt(n_repeat)")
    if meaningful is not None:
        lines.append(f"    - 임상적 의미 변화 기준(추정): {meaningful:.2f} {unit}")
    lines.append(f"    - 판정: {warn_msg}")
    lines.append("")
    lines.append("[4] 연구기간 vs 자연 진행률 tradeoff")
    lines.append(f"    - 적용 하위집단: {subgroup}")
    lines.append(f"    - 가정 연간 진행률(추정): {slope:+.3f} {unit}/년")
    ttm = tradeoff["time_to_mdc_months"]
    if math.isinf(ttm):
        lines.append("    - MDC 도달 예상기간: 진행률 0 가정 — 산출 불가")
    else:
        lines.append(f"    - MDC 만큼 변화 도달 예상기간: 약 {ttm:.1f}개월")
        if ttm > sched["total_months"]:
            lines.append(f"      ※ 계획 추적기간({sched['total_months']:.0f}개월)이 "
                         f"MDC 도달기간보다 짧습니다 — 기간 연장 검토 권장")
        else:
            lines.append(f"      ※ 계획 추적기간({sched['total_months']:.0f}개월) 내 "
                         f"검출 가능 수준 도달 예상")
    lines.append("")
    lines.append("[5] 결측·visit window 반영 검정력 (Monte-Carlo)")
    lines.append(f"    - 표본 수: {n_subjects}명, 가정 결측률: {missing_rate*100:.0f}%, "
                 f"Monte-Carlo {n_sim}회")
    lines.append("    (5-1) 집단 평균 진행 검정력 — 집단 수준 가설")
    lines.append(f"      - 이상적(결측 없음) 해석적 검정력: {sim['power_analytic_ideal']*100:.1f}%")
    lines.append(f"      - 결측·window 반영 시뮬 검정력: {sim['power_simulated']*100:.1f}%")
    lines.append(f"      - 검정력 손실: {sim['power_loss']*100:.1f}%p")
    if sim["power_simulated"] < 0.80:
        lines.append("        ※ 검정력 80% 미만 — 표본 수 증대, 추적기간 연장 필요")
    else:
        lines.append("        ※ 검정력 80% 이상 — 집단 수준 진행 검출은 적정")
    lines.append("    (5-2) 피험자 단위 변화 검출률 — 스케줄 품질 핵심 지표")
    lines.append(f"      - baseline~종료 측정 완료율(엔드포인트 완전성): "
                 f"{sim['endpoint_completeness']*100:.1f}%")
    lines.append(f"      - 이상적(결측 0) 변화 검출률: {sim['subject_detect_ideal']*100:.1f}%")
    lines.append(f"      - 결측·window 반영 변화 검출률: {sim['subject_detect_rate']*100:.1f}%")
    lines.append(f"      - 변화 검출률 손실: {sim['detect_rate_loss']*100:.1f}%p")
    if sim["detect_rate_loss"] > 0.05:
        lines.append("        ※ 결측·window로 인한 검출률 손실 5%p 초과 — "
                     "결측 최소화 및 visit window 단축 검토")
    lines.append("")
    lines.append("[6] 예산 (검사비)")
    lines.append(f"    - 1회 검사 단가(추정): {cost_per_visit:,}원")
    lines.append(f"    - 1인당 총 검사비: {cost_per_subject:,}원 "
                 f"({sched['n_visits']}회 x {n_repeat}반복)")
    lines.append(f"    - 전체 예상 검사비({n_subjects}명): {cost_total:,}원")
    if budget_krw is not None:
        diff = budget_krw - cost_total
        if diff >= 0:
            lines.append(f"    - 가용 예산({budget_krw:,}원) 대비: 여유 {diff:,}원")
        else:
            lines.append(f"    - 가용 예산({budget_krw:,}원) 대비: 부족 {abs(diff):,}원 — 설계 조정 필요")
    lines.append("")
    lines.append("[7] 권장 요약")
    recs = []
    if is_warn:
        recs.append("MDC가 의미 변화보다 큼 → 반복측정 횟수 증가 또는 더 정밀한 NIT 검토")
    if sim["power_simulated"] < 0.80:
        recs.append("집단 검정력 부족 → 표본 수 증대 또는 추적기간 연장")
    if sim["detect_rate_loss"] > 0.05:
        recs.append("결측·window로 인한 피험자 변화 검출률 손실 큼 → "
                    "결측 최소화 전략 및 visit window 단축")
    if sim["endpoint_completeness"] < 0.80:
        recs.append("baseline~종료 측정 완료율 80% 미만 → "
                    "추적 강화 및 결측 보완 visit 추가 검토")
    if not math.isinf(ttm) and ttm > sched["total_months"]:
        recs.append("진행 검출 기간 부족 → 추적기간 연장")
    if not recs:
        recs.append("현 스케줄은 MDC·검정력·검출률·기간 측면에서 적정 수준으로 평가됨")
    for i, r in enumerate(recs, 1):
        lines.append(f"    {i}) {r}")
    lines.append("")
    lines.append("-" * 72)
    lines.append(DISCLAIMER)
    lines.append("출처: AASLD MASLD Practice Guidance 2023; "
                 "EASL-EASD-EASO MASLD CPG 2024 (NIT 권고 참조)")
    lines.append("-" * 72)
    return "\n".join(lines)


if __name__ == "__main__":  # 간단 자가검증
    lib = load_nit_library()
    print("NIT 목록:", list_nits(lib))
    print(generate_design_report(lib, "VCTE_LSM", n_sim=100))