"""
flux_core.py — HepatoLipidFlux 핵심 계산 모듈 (순수 Python)

대사성간질환(MASLD/MASH) 동물모델에서 in vivo 간 지질 플럭스를 정량한다.
Streamlit 의존성 없이 단독으로 import / 테스트 가능하도록 설계했다.

포함 기능:
  1. VLDL-TG 분비율   : tyloxapol/poloxamer-407 주사 후 혈장 TG 시계열의
                       선형(초기) 구간 자동 탐지 + 회귀 → 분비율 계산
  2. DNL 플럭스 (MIDA): 전구체 풀 enrichment + 지방산 mass isotopomer
                       분포 → fractional DNL (신생 분율) 및 합성률
  3. 지방산 산화 플럭스: [13C]palmitate 추적자 호기 13CO2 / 케톤체
                       enrichment → 회수율 보정 산화 플럭스
  4. MASH 모델 참조범위 + QC 경고
  5. 코호트 통계 (ANOVA / 혼합효과) + 기전 분류

연구용·참고용 (for research / reference use only). 임상 진단 목적 사용 금지.
"""

from __future__ import annotations

import math
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Sequence, Tuple

import numpy as np
import pandas as pd

try:
    from scipy import optimize, stats as scipy_stats
    _HAS_SCIPY = True
except Exception:  # pragma: no cover - scipy는 requirements에 포함
    _HAS_SCIPY = False

try:
    import statsmodels.formula.api as smf
    import statsmodels.api as sm
    _HAS_STATSMODELS = True
except Exception:  # pragma: no cover
    _HAS_STATSMODELS = False


# ---------------------------------------------------------------------------
# 0. MASH 동물모델 참조범위 (built-in reference)
# ---------------------------------------------------------------------------
# 값은 마우스 문헌 기반의 대표적 범위(연구용 참고치)이며 절대 기준이 아니다.
# vldl_secretion: mg TG / kg body weight / h
# dnl_fraction  : 간 TG 지방산 중 신생합성 분율 (0-1)
# oxidation     : umol palmitate equivalent / kg / h
MASH_MODEL_REFERENCE: Dict[str, Dict[str, Tuple[float, float]]] = {
    "Chow (정상대조)": {
        "vldl_secretion": (180.0, 320.0),
        "dnl_fraction": (0.05, 0.18),
        "oxidation": (220.0, 420.0),
    },
    "HFD (고지방식이)": {
        "vldl_secretion": (260.0, 480.0),
        "dnl_fraction": (0.10, 0.30),
        "oxidation": (180.0, 380.0),
    },
    "GAN diet": {
        "vldl_secretion": (300.0, 560.0),
        "dnl_fraction": (0.18, 0.42),
        "oxidation": (140.0, 320.0),
    },
    "CDAHFD": {
        "vldl_secretion": (120.0, 300.0),
        "dnl_fraction": (0.12, 0.34),
        "oxidation": (120.0, 300.0),
    },
    "MCD (메티오닌·콜린 결핍)": {
        "vldl_secretion": (60.0, 200.0),
        "dnl_fraction": (0.04, 0.16),
        "oxidation": (110.0, 280.0),
    },
    "AMLN": {
        "vldl_secretion": (300.0, 540.0),
        "dnl_fraction": (0.16, 0.40),
        "oxidation": (150.0, 330.0),
    },
    "db/db": {
        "vldl_secretion": (340.0, 620.0),
        "dnl_fraction": (0.20, 0.46),
        "oxidation": (160.0, 360.0),
    },
    "ob/ob": {
        "vldl_secretion": (320.0, 600.0),
        "dnl_fraction": (0.22, 0.48),
        "oxidation": (150.0, 340.0),
    },
}


# ---------------------------------------------------------------------------
# 결과 컨테이너
# ---------------------------------------------------------------------------
@dataclass
class VLDLResult:
    animal_id: str
    group: str
    secretion_rate_raw: float          # mg TG / dL / h  (혈장 농도 기준 기울기)
    secretion_rate_norm: float         # mg TG / kg BW / h (정규화)
    secretion_rate_lean: Optional[float]  # mg TG / kg lean mass / h
    r_squared: float
    n_points_used: int
    window_minutes: Tuple[float, float]
    intercept: float
    nonlinear_warning: bool
    warnings: List[str] = field(default_factory=list)


@dataclass
class DNLResult:
    animal_id: str
    group: str
    fractional_dnl: float              # 신생합성 분율 (0-1)
    synthesis_rate: Optional[float]    # mg newly-made TG-FA / (조직 g) (선택)
    precursor_enrichment: float        # 전구체 풀 enrichment 가정값 (0-1)
    n_label_assumed: int               # MIDA에 가정한 라벨 가능 위치 수
    sensitivity: Dict[str, float] = field(default_factory=dict)
    warnings: List[str] = field(default_factory=list)


@dataclass
class OxidationResult:
    animal_id: str
    group: str
    oxidation_flux: float              # umol / kg / h (회수율 보정)
    recovery_fraction: float
    raw_flux: float
    warnings: List[str] = field(default_factory=list)


# ---------------------------------------------------------------------------
# 1. VLDL-TG 분비율
# ---------------------------------------------------------------------------
def _linregress(x: np.ndarray, y: np.ndarray) -> Tuple[float, float, float]:
    """단순 선형회귀. (slope, intercept, r_squared) 반환. scipy 없이 동작."""
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    n = len(x)
    if n < 2:
        return float("nan"), float("nan"), float("nan")
    xm, ym = x.mean(), y.mean()
    sxx = np.sum((x - xm) ** 2)
    sxy = np.sum((x - xm) * (y - ym))
    if sxx == 0:
        return float("nan"), float("nan"), float("nan")
    slope = sxy / sxx
    intercept = ym - slope * xm
    ss_tot = np.sum((y - ym) ** 2)
    ss_res = np.sum((y - (slope * x + intercept)) ** 2)
    r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
    return float(slope), float(intercept), float(r2)


def detect_linear_window(
    times: Sequence[float],
    concentrations: Sequence[float],
    min_points: int = 3,
    r2_threshold: float = 0.95,
) -> Tuple[int, int, float, float, float]:
    """
    tyloxapol/poloxamer 주사 후 혈장 TG 시계열에서 초기 선형 구간을 자동 탐지한다.

    전략: 시작점(0번 인덱스)부터 가능한 한 많은 연속 포인트를 포함하되,
    r^2 >= r2_threshold 를 유지하는 가장 긴 윈도우를 선택한다.
    이는 후기 포화(saturation) 구간을 자동으로 배제한다.

    반환: (start_idx, end_idx_inclusive, slope, intercept, r_squared)
    """
    t = np.asarray(times, dtype=float)
    c = np.asarray(concentrations, dtype=float)
    order = np.argsort(t)
    t, c = t[order], c[order]
    n = len(t)
    if n < min_points:
        slope, intercept, r2 = _linregress(t, c)
        return 0, n - 1, slope, intercept, r2

    # 0번 시점부터 시작하는 윈도우를 늘려가며, r2 임계를 만족하는
    # 가장 긴 구간을 선택한다. 잡음에 의한 일시적 r2 하락에 견고하도록
    # 첫 실패에서 즉시 중단하지 않고, 연속 2회 실패 시에만 중단한다
    # (포화 구간은 r2를 지속적으로 떨어뜨리므로 구분된다).
    s0, i0, r0 = _linregress(t[:min_points], c[:min_points])
    best = (0, min_points - 1, s0, i0, r0)
    best_end_ok = min_points - 1 if (not math.isnan(r0) and r0 >= r2_threshold) else -1
    consecutive_fail = 0
    for end in range(min_points, n):
        slope, intercept, r2 = _linregress(t[: end + 1], c[: end + 1])
        if not math.isnan(r2) and r2 >= r2_threshold:
            best = (0, end, slope, intercept, r2)
            best_end_ok = end
            consecutive_fail = 0
        else:
            consecutive_fail += 1
            if consecutive_fail >= 2:
                break
    # 임계를 한 번도 만족 못했으면 가능한 최장 구간(전체)을 반환
    if best_end_ok < 0:
        slope, intercept, r2 = _linregress(t, c)
        return 0, n - 1, slope, intercept, r2
    return best


def compute_vldl_secretion(
    df: pd.DataFrame,
    body_weight_g: Optional[float] = None,
    lean_mass_g: Optional[float] = None,
    animal_id: str = "NA",
    group: str = "NA",
    r2_threshold: float = 0.95,
    time_col: str = "time_min",
    conc_col: str = "plasma_tg_mg_dl",
    plasma_volume_ml_per_g: float = 0.04,
) -> VLDLResult:
    """
    한 마리의 혈장 TG 시계열로부터 VLDL-TG 분비율을 계산한다.

    분비율(raw) 단위: mg TG / dL plasma / h  (농도 증가 기울기)
    정규화 분비율: 혈장량을 체중에 비례한다고 가정하여 체중(kg)당 mg/h 로 변환.
        plasma_volume_ml = body_weight_g * plasma_volume_ml_per_g
        secretion(mg/h) = slope(mg/dL/min) * 60 * plasma_volume_dL
        norm(mg/kg/h)  = secretion(mg/h) / (body_weight_g/1000)
    """
    warnings: List[str] = []
    times = df[time_col].to_numpy(dtype=float)
    conc = df[conc_col].to_numpy(dtype=float)

    if len(times) < 2:
        raise ValueError(f"{animal_id}: 시계열 포인트가 부족하다 (>=2 필요).")

    start, end, slope, intercept, r2 = detect_linear_window(
        times, conc, r2_threshold=r2_threshold
    )
    n_used = end - start + 1
    n_total = len(times)

    nonlinear = False
    if n_used < n_total:
        nonlinear = True
        warnings.append(
            f"비선형 포화 구간 감지: 전체 {n_total}점 중 초기 {n_used}점만 회귀에 사용 "
            f"(후기 {n_total - n_used}점 제외). 포화 구간은 분비율을 과소평가한다."
        )
    if not math.isnan(r2) and r2 < r2_threshold:
        warnings.append(
            f"선형구간 r^2={r2:.3f} < {r2_threshold}: 데이터 잡음/포화로 회귀 신뢰도 낮음."
        )
    if slope <= 0:
        warnings.append(
            "기울기가 0 이하: TG 농도가 증가하지 않음 — 추적자 주사 실패 또는 채혈 오류 의심."
        )

    sorted_t = np.sort(times)
    window = (float(sorted_t[start]), float(sorted_t[end]))

    # raw: mg/dL/h
    secretion_raw = slope * 60.0

    secretion_norm = float("nan")
    secretion_lean = None
    if body_weight_g is not None and body_weight_g > 0:
        plasma_volume_ml = body_weight_g * plasma_volume_ml_per_g
        plasma_volume_dl = plasma_volume_ml / 100.0
        secretion_mg_per_h = slope * 60.0 * plasma_volume_dl
        secretion_norm = secretion_mg_per_h / (body_weight_g / 1000.0)
        if lean_mass_g is not None and lean_mass_g > 0:
            secretion_lean = secretion_mg_per_h / (lean_mass_g / 1000.0)
    else:
        warnings.append("체중 미입력: 체중 정규화 분비율을 계산할 수 없음 (raw 값만 유효).")

    return VLDLResult(
        animal_id=animal_id,
        group=group,
        secretion_rate_raw=float(secretion_raw),
        secretion_rate_norm=float(secretion_norm),
        secretion_rate_lean=secretion_lean,
        r_squared=float(r2),
        n_points_used=int(n_used),
        window_minutes=window,
        intercept=float(intercept),
        nonlinear_warning=nonlinear,
        warnings=warnings,
    )


# ---------------------------------------------------------------------------
# 2. DNL 플럭스 — MIDA (Mass Isotopomer Distribution Analysis)
# ---------------------------------------------------------------------------
def binomial_isotopomer_distribution(n: int, p: float) -> np.ndarray:
    """
    n개의 라벨 가능 위치, 각 위치 라벨 확률 p 일 때
    질량 isotopomer 분포 M0, M1, ..., Mn 를 이항분포로 반환한다.

    MIDA 기본 모델: 신생합성된 분자는 전구체 풀 enrichment(p)로부터
    n개 위치에 독립적으로 라벨을 받는다고 가정한다.
    """
    if not (0.0 <= p <= 1.0):
        raise ValueError("precursor enrichment p 는 0..1 범위여야 한다.")
    k = np.arange(n + 1)
    # 이항계수 * p^k * (1-p)^(n-k)
    coeff = np.array([math.comb(n, int(ki)) for ki in k], dtype=float)
    dist = coeff * (p ** k) * ((1.0 - p) ** (n - k))
    return dist / dist.sum()


def compute_fractional_dnl(
    measured_distribution: Sequence[float],
    precursor_enrichment: float,
    n_label_sites: int = 22,
    natural_abundance: Optional[Sequence[float]] = None,
) -> Tuple[float, np.ndarray]:
    """
    측정된 지방산 mass isotopomer 분포로부터 fractional DNL 을 계산한다.

    모델: 측정분포 = (1 - f) * 기존(natural abundance) 분포
                    + f * 신생합성(MIDA 이항) 분포
    f (fractional DNL) 를 최소제곱으로 추정한다.

    palmitate(C16)의 경우 2H2O 라벨 모델에서 라벨 가능 H 위치는
    통상 22개 내외로 가정한다 (n_label_sites).

    반환: (fractional_dnl, 적합된_분포)
    """
    measured = np.asarray(measured_distribution, dtype=float)
    measured = measured / measured.sum()  # 정규화
    m_len = len(measured)

    new_dist_full = binomial_isotopomer_distribution(n_label_sites, precursor_enrichment)
    new_dist = new_dist_full[:m_len]
    if new_dist.sum() > 0:
        new_dist = new_dist / new_dist.sum()

    if natural_abundance is None:
        # 자연존재비: M0 우세, 미세한 M1/M2 (13C 자연존재 1.1%)
        nat = np.zeros(m_len)
        nat[0] = 0.955
        if m_len > 1:
            nat[1] = 0.040
        if m_len > 2:
            nat[2] = 0.005
        nat = nat / nat.sum()
    else:
        nat = np.asarray(natural_abundance, dtype=float)[:m_len]
        nat = nat / nat.sum()

    # measured = (1-f)*nat + f*new_dist  →  measured - nat = f*(new_dist - nat)
    design = (new_dist - nat)
    target = (measured - nat)
    denom = float(np.dot(design, design))
    if denom == 0:
        f = 0.0
    else:
        f = float(np.dot(design, target) / denom)
    f = max(0.0, min(1.0, f))  # 0..1 클리핑
    fitted = (1.0 - f) * nat + f * new_dist
    return f, fitted


def compute_dnl(
    measured_distribution: Sequence[float],
    precursor_enrichment: float,
    animal_id: str = "NA",
    group: str = "NA",
    n_label_sites: int = 22,
    hepatic_tg_fa_mass_mg: Optional[float] = None,
    label_duration_h: Optional[float] = None,
) -> DNLResult:
    """
    한 마리의 지방산 isotopomer 분포로부터 fractional DNL 과
    (선택적으로) 합성률을 계산하고, 전구체 풀 enrichment 가정에 대한
    민감도 분석을 함께 제공한다.
    """
    warnings: List[str] = []
    if precursor_enrichment <= 0:
        warnings.append(
            "전구체 풀 enrichment 가정값이 0 이하: MIDA 계산 불가능 — 2H2O/13C "
            "전구체 enrichment 측정값을 반드시 입력해야 함."
        )
        precursor_enrichment = 1e-6

    f, _ = compute_fractional_dnl(
        measured_distribution, precursor_enrichment, n_label_sites
    )

    # 민감도: 전구체 enrichment ±20% 변동 시 fractional DNL 변화
    sensitivity: Dict[str, float] = {}
    for label, factor in (("precursor_-20%", 0.8), ("precursor_+20%", 1.2)):
        p_alt = min(1.0, max(1e-6, precursor_enrichment * factor))
        f_alt, _ = compute_fractional_dnl(
            measured_distribution, p_alt, n_label_sites
        )
        sensitivity[label] = float(f_alt)
    spread = abs(sensitivity["precursor_+20%"] - sensitivity["precursor_-20%"])
    if f > 0 and spread / max(f, 1e-9) > 0.3:
        warnings.append(
            f"전구체 enrichment 가정에 민감 (±20% 가정 시 DNL이 {spread:.3f} 변동). "
            "전구체 풀 enrichment 측정 정확도를 확인할 것."
        )

    synthesis_rate = None
    if (
        hepatic_tg_fa_mass_mg is not None
        and hepatic_tg_fa_mass_mg > 0
        and label_duration_h is not None
        and label_duration_h > 0
    ):
        # 신생합성 TG-FA 질량 / 라벨 노출시간
        synthesis_rate = (f * hepatic_tg_fa_mass_mg) / label_duration_h

    warnings.append(
        f"전구체 풀 enrichment 가정 = {precursor_enrichment:.4f} "
        f"(2H2O body water 또는 13C 전구체 plateau enrichment). 결과는 이 가정에 의존함."
    )

    return DNLResult(
        animal_id=animal_id,
        group=group,
        fractional_dnl=float(f),
        synthesis_rate=synthesis_rate,
        precursor_enrichment=float(precursor_enrichment),
        n_label_assumed=int(n_label_sites),
        sensitivity=sensitivity,
        warnings=warnings,
    )


# ---------------------------------------------------------------------------
# 3. 지방산 산화 플럭스
# ---------------------------------------------------------------------------
def compute_oxidation_flux(
    tracer_dose_umol: float,
    co2_label_recovery_fraction: float,
    ketone_label_fraction: float = 0.0,
    body_weight_g: float = 25.0,
    bicarbonate_recovery: float = 0.80,
    measurement_duration_h: float = 1.0,
    animal_id: str = "NA",
    group: str = "NA",
) -> OxidationResult:
    """
    [13C]palmitate 추적자 호기 13CO2 (및 케톤체) enrichment 로부터
    회수율 보정 지방산 산화 플럭스를 계산한다.

    raw_flux  = (투여 추적자 umol * (CO2 회수분율 + 케톤체 분율)) /
                (체중kg * 측정시간h)
    보정 플럭스 = raw_flux / bicarbonate_recovery
        (호기 13CO2는 체내 bicarbonate pool 에 일부 격리되므로
         bicarbonate_recovery(<1) 로 보정하여 과소평가를 교정한다.)
    """
    warnings: List[str] = []
    if not (0.0 <= co2_label_recovery_fraction <= 1.0):
        warnings.append("CO2 회수분율이 0..1 범위를 벗어남 — 입력값 확인 필요.")
    if not (0.0 < bicarbonate_recovery <= 1.0):
        warnings.append("bicarbonate 회수율은 0..1 범위여야 함; 기본 0.80 사용 권장.")
        bicarbonate_recovery = 0.80
    if tracer_dose_umol <= 0:
        raise ValueError("추적자 투여량(umol)은 양수여야 한다.")
    if body_weight_g <= 0:
        raise ValueError("체중(g)은 양수여야 한다.")

    total_recovered = co2_label_recovery_fraction + ketone_label_fraction
    bw_kg = body_weight_g / 1000.0
    raw_flux = (tracer_dose_umol * total_recovered) / (
        bw_kg * measurement_duration_h
    )
    corrected_flux = raw_flux / bicarbonate_recovery

    if total_recovered < 0.02:
        warnings.append(
            "표지 회수율이 매우 낮음 (<2%): 산화 경로 활성 저하 또는 측정 누락 가능."
        )
    if total_recovered > 0.9:
        warnings.append(
            "표지 회수율이 비정상적으로 높음 (>90%): 회수율 보정 가정 재검토 필요."
        )

    return OxidationResult(
        animal_id=animal_id,
        group=group,
        oxidation_flux=float(corrected_flux),
        recovery_fraction=float(total_recovered),
        raw_flux=float(raw_flux),
        warnings=warnings,
    )


# ---------------------------------------------------------------------------
# 4. QC — 참조범위 대비 경고
# ---------------------------------------------------------------------------
def qc_against_reference(
    metric: str,
    value: float,
    model: str,
) -> List[str]:
    """
    측정값을 MASH 모델 참조범위와 비교하여 경고 메시지를 반환한다.
    metric: 'vldl_secretion' | 'dnl_fraction' | 'oxidation'
    """
    msgs: List[str] = []
    if model not in MASH_MODEL_REFERENCE:
        msgs.append(f"'{model}' 모델은 내장 참조범위에 없음 — QC 비교 생략.")
        return msgs
    ref = MASH_MODEL_REFERENCE[model].get(metric)
    if ref is None:
        return msgs
    lo, hi = ref
    if math.isnan(value):
        msgs.append(f"{metric}: 값이 NaN — 계산 실패 또는 입력 누락.")
        return msgs
    if value < lo:
        msgs.append(
            f"{metric}={value:.3g} 가 {model} 참조범위({lo:g}-{hi:g}) 하한 미만."
        )
    elif value > hi:
        msgs.append(
            f"{metric}={value:.3g} 가 {model} 참조범위({lo:g}-{hi:g}) 상한 초과."
        )
    return msgs


# ---------------------------------------------------------------------------
# 5. 코호트 통계 + 기전 분류
# ---------------------------------------------------------------------------
def group_summary(df: pd.DataFrame, value_col: str, group_col: str = "group") -> pd.DataFrame:
    """그룹별 기술통계 (n, mean, sd, sem)."""
    rows = []
    for g, sub in df.groupby(group_col):
        vals = sub[value_col].dropna().to_numpy(dtype=float)
        n = len(vals)
        mean = float(np.mean(vals)) if n else float("nan")
        sd = float(np.std(vals, ddof=1)) if n > 1 else float("nan")
        sem = sd / math.sqrt(n) if n > 1 else float("nan")
        rows.append({"group": g, "n": n, "mean": mean, "sd": sd, "sem": sem})
    return pd.DataFrame(rows)


def one_way_anova(df: pd.DataFrame, value_col: str, group_col: str = "group") -> Dict[str, float]:
    """
    일원배치 ANOVA. scipy 가 있으면 scipy.stats.f_oneway 사용,
    없으면 순수 Python 으로 F, p 근사 계산.
    반환: {'F': ..., 'p_value': ..., 'df_between': ..., 'df_within': ...}
    """
    groups = [
        sub[value_col].dropna().to_numpy(dtype=float)
        for _, sub in df.groupby(group_col)
    ]
    groups = [g for g in groups if len(g) > 0]
    if len(groups) < 2:
        return {"F": float("nan"), "p_value": float("nan"),
                "df_between": float("nan"), "df_within": float("nan")}

    if _HAS_SCIPY:
        F, p = scipy_stats.f_oneway(*groups)
        k = len(groups)
        n_total = sum(len(g) for g in groups)
        return {
            "F": float(F),
            "p_value": float(p),
            "df_between": float(k - 1),
            "df_within": float(n_total - k),
        }

    # 순수 Python fallback
    grand = np.concatenate(groups)
    grand_mean = grand.mean()
    k = len(groups)
    n_total = len(grand)
    ss_between = sum(len(g) * (g.mean() - grand_mean) ** 2 for g in groups)
    ss_within = sum(np.sum((g - g.mean()) ** 2) for g in groups)
    df_b = k - 1
    df_w = n_total - k
    if df_w <= 0 or ss_within == 0:
        return {"F": float("nan"), "p_value": float("nan"),
                "df_between": float(df_b), "df_within": float(df_w)}
    ms_b = ss_between / df_b
    ms_w = ss_within / df_w
    F = ms_b / ms_w
    return {
        "F": float(F),
        "p_value": float("nan"),  # scipy 없으면 p 계산 생략
        "df_between": float(df_b),
        "df_within": float(df_w),
    }


def mixed_effects_model(
    df: pd.DataFrame,
    value_col: str,
    group_col: str = "group",
    random_col: str = "batch",
) -> Optional[Dict[str, object]]:
    """
    혼합효과 모델: 고정효과=group, 임의효과=batch(또는 cage).
    statsmodels 가 있을 때만 동작. 없으면 None 반환.
    """
    if not _HAS_STATSMODELS:
        return None
    if random_col not in df.columns:
        return None
    work = df[[value_col, group_col, random_col]].dropna()
    if work[group_col].nunique() < 2 or work[random_col].nunique() < 2:
        return None
    try:
        md = smf.mixedlm(
            f"{value_col} ~ C({group_col})", work, groups=work[random_col]
        )
        res = md.fit(reml=True, method="lbfgs")
        return {
            "summary": str(res.summary()),
            "params": res.params.to_dict(),
            "pvalues": res.pvalues.to_dict(),
        }
    except Exception as exc:  # pragma: no cover
        return {"error": str(exc)}


def classify_mechanism(
    vldl_change_pct: float,
    dnl_change_pct: float,
    oxidation_change_pct: float,
    threshold_pct: float = 15.0,
) -> Dict[str, object]:
    """
    대조군 대비 변화율(%)을 받아 개입(약물/식이)의 작용 기전을 분류한다.

    분류 카테고리:
      - 'synthesis-suppressing'  : DNL 분율 감소가 우세
      - 'secretion-altering'     : VLDL 분비율 변화가 우세
      - 'oxidation-promoting'    : 산화 플럭스 증가가 우세
      - 'mixed'                  : 여러 축이 동시에 유의
      - 'no clear effect'        : 모든 축이 임계 미만

    반환: {'classification': ..., 'rationale': ..., 'flags': {...}}
    """
    flags = {
        "synthesis-suppressing": dnl_change_pct <= -threshold_pct,
        "secretion-altering": abs(vldl_change_pct) >= threshold_pct,
        "oxidation-promoting": oxidation_change_pct >= threshold_pct,
    }
    active = [k for k, v in flags.items() if v]

    if not active:
        classification = "no clear effect"
        rationale = (
            f"세 축 모두 변화율이 임계({threshold_pct:g}%) 미만 — 뚜렷한 플럭스 기전 없음."
        )
    elif len(active) == 1:
        classification = active[0]
        magnitudes = {
            "synthesis-suppressing": dnl_change_pct,
            "secretion-altering": vldl_change_pct,
            "oxidation-promoting": oxidation_change_pct,
        }
        rationale = (
            f"{classification} 단일 축 우세 (해당 변화율 {magnitudes[active[0]]:+.1f}%)."
        )
    else:
        # 가장 큰 절대 변화율을 주 기전으로
        magnitudes = {
            "synthesis-suppressing": abs(dnl_change_pct),
            "secretion-altering": abs(vldl_change_pct),
            "oxidation-promoting": abs(oxidation_change_pct),
        }
        primary = max(active, key=lambda k: magnitudes[k])
        classification = "mixed"
        rationale = (
            f"복수 축 동시 변화 ({', '.join(active)}). 주 기전 추정: {primary}."
        )

    return {
        "classification": classification,
        "rationale": rationale,
        "flags": flags,
        "changes_pct": {
            "vldl_secretion": vldl_change_pct,
            "dnl_fraction": dnl_change_pct,
            "oxidation": oxidation_change_pct,
        },
    }


# ---------------------------------------------------------------------------
# 배치 처리 헬퍼 — 데이터프레임 → 결과 테이블
# ---------------------------------------------------------------------------
def batch_vldl(
    ts_df: pd.DataFrame,
    meta_df: pd.DataFrame,
    r2_threshold: float = 0.95,
) -> pd.DataFrame:
    """
    여러 동물의 혈장 TG 시계열(long format)을 일괄 처리한다.

    ts_df 컬럼   : animal_id, group, time_min, plasma_tg_mg_dl
    meta_df 컬럼 : animal_id, group, body_weight_g, lean_mass_g (lean 선택)
    """
    meta = meta_df.set_index("animal_id")
    rows = []
    for animal_id, sub in ts_df.groupby("animal_id"):
        group = str(sub["group"].iloc[0])
        bw = None
        lean = None
        if animal_id in meta.index:
            bw = float(meta.loc[animal_id, "body_weight_g"])
            if "lean_mass_g" in meta.columns and not pd.isna(
                meta.loc[animal_id, "lean_mass_g"]
            ):
                lean = float(meta.loc[animal_id, "lean_mass_g"])
        res = compute_vldl_secretion(
            sub, body_weight_g=bw, lean_mass_g=lean,
            animal_id=str(animal_id), group=group, r2_threshold=r2_threshold,
        )
        rows.append({
            "animal_id": res.animal_id,
            "group": res.group,
            "secretion_rate_norm": res.secretion_rate_norm,
            "secretion_rate_raw": res.secretion_rate_raw,
            "secretion_rate_lean": res.secretion_rate_lean,
            "r_squared": res.r_squared,
            "n_points_used": res.n_points_used,
            "window_start_min": res.window_minutes[0],
            "window_end_min": res.window_minutes[1],
            "nonlinear_warning": res.nonlinear_warning,
            "n_warnings": len(res.warnings),
        })
    return pd.DataFrame(rows)


def batch_dnl(dnl_df: pd.DataFrame) -> pd.DataFrame:
    """
    여러 동물의 isotopomer 분포(wide format)를 일괄 처리한다.

    dnl_df 컬럼: animal_id, group, precursor_enrichment,
                 n_label_sites(선택), M0, M1, M2, M3, M4 ...
    """
    m_cols = [c for c in dnl_df.columns if c.upper().startswith("M") and c[1:].isdigit()]
    m_cols = sorted(m_cols, key=lambda c: int(c[1:]))
    rows = []
    for _, r in dnl_df.iterrows():
        dist = [float(r[c]) for c in m_cols]
        n_sites = int(r["n_label_sites"]) if "n_label_sites" in dnl_df.columns else 22
        res = compute_dnl(
            dist,
            precursor_enrichment=float(r["precursor_enrichment"]),
            animal_id=str(r["animal_id"]),
            group=str(r["group"]),
            n_label_sites=n_sites,
        )
        rows.append({
            "animal_id": res.animal_id,
            "group": res.group,
            "fractional_dnl": res.fractional_dnl,
            "precursor_enrichment": res.precursor_enrichment,
            "n_label_sites": res.n_label_assumed,
            "sens_minus20": res.sensitivity.get("precursor_-20%", float("nan")),
            "sens_plus20": res.sensitivity.get("precursor_+20%", float("nan")),
        })
    return pd.DataFrame(rows)


def batch_oxidation(ox_df: pd.DataFrame) -> pd.DataFrame:
    """
    여러 동물의 산화 추적자 데이터(wide format)를 일괄 처리한다.

    ox_df 컬럼: animal_id, group, tracer_dose_umol, co2_recovery_fraction,
                ketone_label_fraction(선택), body_weight_g,
                bicarbonate_recovery(선택), duration_h(선택)
    """
    rows = []
    for _, r in ox_df.iterrows():
        res = compute_oxidation_flux(
            tracer_dose_umol=float(r["tracer_dose_umol"]),
            co2_label_recovery_fraction=float(r["co2_recovery_fraction"]),
            ketone_label_fraction=float(r.get("ketone_label_fraction", 0.0) or 0.0),
            body_weight_g=float(r["body_weight_g"]),
            bicarbonate_recovery=float(r.get("bicarbonate_recovery", 0.80) or 0.80),
            measurement_duration_h=float(r.get("duration_h", 1.0) or 1.0),
            animal_id=str(r["animal_id"]),
            group=str(r["group"]),
        )
        rows.append({
            "animal_id": res.animal_id,
            "group": res.group,
            "oxidation_flux": res.oxidation_flux,
            "raw_flux": res.raw_flux,
            "recovery_fraction": res.recovery_fraction,
        })
    return pd.DataFrame(rows)


if __name__ == "__main__":  # 간단한 자체 점검
    print("flux_core 자체 점검")
    dist = binomial_isotopomer_distribution(22, 0.04)
    print("  isotopomer dist sum =", round(dist.sum(), 6))
    f, _ = compute_fractional_dnl([0.70, 0.20, 0.08, 0.015, 0.005], 0.04, 22)
    print("  fractional DNL 예시 =", round(f, 4))
    print("내장 MASH 모델 수 =", len(MASH_MODEL_REFERENCE))