"""breath.py — ¹³C 호기검사(¹³C breath test) 분석.

δ¹³C(‰, vs PDB)의 시점별 측정값으로부터 ¹³C 기질 산화 동태를 산출한다:
  - δ¹³C → 시점별 PDR (Percent Dose Recovery /h, %/h)
  - cumulative PDR (cPDR) 40분/120분
  - 산화율(oxidation rate)
  - 제거 속도상수 kel·반감기 (선형 또는 비선형 fit)

────────────────────────────────────────────────────────────────────────
δ¹³C → PDR 변환식 (표준 호기검사 방법론)
────────────────────────────────────────────────────────────────────────
호기 ¹³CO₂ 농축도를 δ 표기로 측정한다. 시점 t 에서 기저(δ0) 대비 변화량을
DOB(t) = δ(t) - δ0  (Delta Over Baseline, ‰) 라 하자.

기질 투여량 D (mmol ¹³C, 즉 mol substrate × C수 × 동위원소 표지비) 에 대해
회수율(시점별 PDR, %/h)은:

  PDR(t) [%/h]
    = ( VCO2 * (DOB(t)/1000) * R_PDB ) / D * 100
      ──────────  (여기서 0.011237 = R_PDB, PDB 표준 ¹³C/¹²C 비)

  실무 표준식(예: Schadewaldt, Sanaka, Ghoos 등 ¹³C-octanoate/acetate breath
  test 문헌)에서는 다음 형태로 정리한다:

  PDR(t) = DOB(t) * VCO2 * R_PDB / D * 100   [percent of dose per hour]

여기서
  VCO2 : CO₂ 생성률 [mmol/h]
  DOB  : delta over baseline [‰] → 비율로 환산 시 /1000
  R_PDB: 0.0112372 (PDB 국제표준 ¹³C/¹²C 비; 절대비)
  D    : 투여 ¹³C 용량 [mmol]

VCO2 입력 옵션:
  (a) 측정값 직접 입력 [mmol/h], 또는
  (b) 체표면적(BSA) 기반 추정: VCO2 = 300 mmol/m²/h × BSA(m²).
      (널리 쓰이는 표준 가정값 300 mmol CO₂ /m² /h; 가정값은 메타데이터로 강제 기록)

BSA(Du Bois): BSA(m²) = 0.007184 × W(kg)^0.425 × H(cm)^0.725

cumulative PDR (cPDR, %): PDR(t) 를 시간에 대해 적분(사다리꼴) → 누적 회수율.
산화율(oxidation rate) 은 시점별 PDR 그 자체(%/h)로 보고하거나, 절대 기질
산화량 [mmol/h] = PDR(t)/100 × D 로 환산.

kel·반감기: 상승–하강 곡선의 하강부(또는 전체 cPDR 접근)에서 1차 동태 가정.
  본 모듈은 (i) PDR 하강부 단일지수 비선형 fit (scipy.optimize.curve_fit), 또는
  (ii) 로그-선형 회귀(폴백)로 kel 추정, t½ = ln(2)/kel.

면책: 연구용·참고용. 임상 의사결정에 직접 사용 금지.
"""
from __future__ import annotations

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

import numpy as np

import auc as _auc

# PDB(Pee Dee Belemnite) 국제표준 ¹³C/¹²C 절대비
R_PDB = 0.0112372

# BSA 기반 CO2 생성률 표준 가정값 [mmol/m²/h]
DEFAULT_VCO2_PER_BSA = 300.0


def du_bois_bsa(weight_kg: float, height_cm: float) -> float:
    """Du Bois & Du Bois 체표면적(m²)."""
    return 0.007184 * (weight_kg ** 0.425) * (height_cm ** 0.725)


def estimate_vco2(weight_kg: float, height_cm: float,
                  vco2_per_bsa: float = DEFAULT_VCO2_PER_BSA) -> float:
    """BSA 기반 VCO2 추정 [mmol/h] = vco2_per_bsa × BSA."""
    bsa = du_bois_bsa(weight_kg, height_cm)
    return vco2_per_bsa * bsa


def delta_to_pdr(times_min: Sequence[float],
                 delta13c: Sequence[float],
                 dose_mmol_13c: float,
                 vco2_mmol_h: float,
                 baseline_delta: Optional[float] = None) -> Tuple[np.ndarray, np.ndarray]:
    """δ¹³C(‰) 시점값 → 시점별 PDR(%/h).

    Returns (times_h, pdr_pct_per_h).

    PDR(t) = DOB(t)/1000 * VCO2 * R_PDB / D * 100
    DOB(t) = δ(t) - δ0
    """
    t_min = np.asarray(times_min, dtype=float)
    d = np.asarray(delta13c, dtype=float)
    if baseline_delta is None:
        baseline_delta = float(d[0])
    dob = d - baseline_delta  # ‰
    pdr = (dob / 1000.0) * vco2_mmol_h * R_PDB / dose_mmol_13c * 100.0
    t_h = t_min / 60.0
    return t_h, pdr


def cumulative_pdr(times_min: Sequence[float],
                   pdr_pct_per_h: Sequence[float],
                   t_end_min: float) -> float:
    """cPDR(%) = 0..t_end 의 PDR(%/h) 사다리꼴 적분.

    PDR 단위가 %/h 이고 시간 적분은 h 단위여야 하므로 분→시 변환 후 적분.
    """
    t_h = np.asarray(times_min, dtype=float) / 60.0
    return _auc.cumulative_auc_to(t_h, pdr_pct_per_h, t_end_min / 60.0)


def _mono_exp(t, A, k, C):
    """단일지수 감쇠 모델: y = A*exp(-k*t) + C."""
    return A * np.exp(-k * t) + C


def fit_kel_halflife(times_h: Sequence[float],
                     pdr_pct_per_h: Sequence[float],
                     use_nonlinear: bool = True) -> Dict[str, object]:
    """PDR 하강부에서 제거 속도상수 kel(/h)·반감기(h) 추정.

    use_nonlinear=True 이고 scipy 가용하면 단일지수 비선형 fit,
    실패하거나 미설치면 로그-선형 회귀로 폴백.

    반환 dict: {kel_per_h, half_life_h, method, r_or_note}
    """
    t = np.asarray(times_h, dtype=float)
    y = np.asarray(pdr_pct_per_h, dtype=float)

    # peak 이후 하강부만 사용
    if len(t) < 3:
        return {"kel_per_h": float("nan"), "half_life_h": float("nan"),
                "method": "none", "note": "시점 부족(<3)"}
    peak_idx = int(np.argmax(y))
    t_d = t[peak_idx:]
    y_d = y[peak_idx:]
    if len(t_d) < 3:
        # 하강부 부족 → 전체 사용 시도
        t_d, y_d = t, y

    method = "none"
    note = ""
    kel = float("nan")

    if use_nonlinear:
        try:
            from scipy.optimize import curve_fit  # noqa
            A0 = max(float(np.max(y_d) - np.min(y_d)), 1e-6)
            p0 = [A0, 1.0, max(float(np.min(y_d)), 0.0)]
            popt, _ = curve_fit(_mono_exp, t_d - t_d[0], y_d, p0=p0, maxfev=20000)
            kel = float(abs(popt[1]))
            method = "nonlinear_monoexp"
            note = f"A={popt[0]:.4g}, C={popt[2]:.4g}"
        except Exception as e:  # 폴백
            note = f"비선형 fit 실패→선형폴백 ({type(e).__name__})"
            method = None  # 아래 선형으로

    if method in (None, "none"):
        # 로그-선형 회귀: ln(y - C0) ~ -kel * t ; C0 baseline 추정값
        c0 = max(float(np.min(y_d)) - 1e-6, 0.0)
        yy = y_d - c0
        valid = yy > 0
        if valid.sum() >= 2:
            slope, _intercept = np.polyfit(t_d[valid] - t_d[0], np.log(yy[valid]), 1)
            kel = float(abs(slope))
            method = "loglinear" if method is None else "loglinear"
            if not note:
                note = "loglinear regression"
            else:
                note = note + "; loglinear"
        else:
            method = "failed"
            note = (note + "; 유효 양수점 부족") if note else "유효 양수점 부족"

    half_life = float(np.log(2) / kel) if (kel and kel > 0 and np.isfinite(kel)) else float("nan")
    return {"kel_per_h": kel, "half_life_h": half_life, "method": method, "note": note}


@dataclass
class BreathResult:
    """¹³C 호기검사 단일 곡선 분석 결과."""
    vco2_mmol_h: float
    dose_mmol_13c: float
    baseline_delta: float
    peak_pdr: float
    tmax_h: float
    cpdr_40min: float
    cpdr_120min: float
    cpdr_total: float
    kel_per_h: float
    half_life_h: float
    fit_method: str
    times_h: List[float] = field(default_factory=list)
    pdr: List[float] = field(default_factory=list)
    flags: List[str] = field(default_factory=list)
    note: str = ""

    def as_dict(self) -> Dict[str, object]:
        return {
            "vco2_mmol_h": self.vco2_mmol_h,
            "dose_mmol_13c": self.dose_mmol_13c,
            "baseline_delta_permil": self.baseline_delta,
            "peak_pdr_pct_per_h": self.peak_pdr,
            "tmax_h": self.tmax_h,
            "cpdr_40min_pct": self.cpdr_40min,
            "cpdr_120min_pct": self.cpdr_120min,
            "cpdr_total_pct": self.cpdr_total,
            "kel_per_h": self.kel_per_h,
            "half_life_h": self.half_life_h,
            "fit_method": self.fit_method,
            "flags": "; ".join(self.flags) if self.flags else "",
            "note": self.note,
        }


def analyze_breath(times_min: Sequence[float],
                   delta13c: Sequence[float],
                   dose_mmol_13c: float,
                   vco2_mmol_h: Optional[float] = None,
                   weight_kg: Optional[float] = None,
                   height_cm: Optional[float] = None,
                   vco2_per_bsa: float = DEFAULT_VCO2_PER_BSA,
                   baseline_delta: Optional[float] = None,
                   use_nonlinear: bool = True,
                   expected_times_min: Optional[Sequence[float]] = None) -> BreathResult:
    """¹³C 호기검사 단일 곡선 전체 분석.

    VCO2: vco2_mmol_h 직접 입력 우선. 없으면 weight/height + vco2_per_bsa 추정.
    """
    note_parts: List[str] = []
    if vco2_mmol_h is None:
        if weight_kg is None or height_cm is None:
            raise ValueError("VCO2 미입력 시 weight_kg, height_cm 필요")
        vco2_mmol_h = estimate_vco2(weight_kg, height_cm, vco2_per_bsa)
        note_parts.append(
            f"VCO2 추정: {vco2_per_bsa} mmol/m²/h × BSA "
            f"{du_bois_bsa(weight_kg, height_cm):.3f} m² = {vco2_mmol_h:.1f} mmol/h"
        )
    else:
        note_parts.append(f"VCO2 측정값 입력: {vco2_mmol_h:.1f} mmol/h")

    t_h, pdr = delta_to_pdr(times_min, delta13c, dose_mmol_13c,
                            vco2_mmol_h, baseline_delta=baseline_delta)
    base_delta = float(delta13c[0]) if baseline_delta is None else float(baseline_delta)

    peak_pdr, tmax = _auc.peak_and_tmax(t_h, pdr)
    cpdr40 = cumulative_pdr(times_min, pdr, 40.0)
    cpdr120 = cumulative_pdr(times_min, pdr, 120.0)
    cpdr_total = _auc.total_auc(t_h, pdr)

    fit = fit_kel_halflife(t_h, pdr, use_nonlinear=use_nonlinear)

    flags: List[str] = []
    # QC: 음수 PDR
    if np.any(np.asarray(pdr) < -1e-9):
        flags.append("음수 PDR 존재(기저 이상/측정오류 가능)")
    # QC: δ¹³C 기저 이상 (PDB 대비 통상 -10 ~ -30‰ 범위 밖이면 플래그)
    if not (-40.0 <= base_delta <= 0.0):
        flags.append(f"δ¹³C 기저 이상치 {base_delta:.1f}‰ (통상 -30~-10‰)")
    # QC: 시점 누락
    if expected_times_min is not None:
        present = set(round(float(x), 3) for x in times_min)
        for et in expected_times_min:
            if round(float(et), 3) not in present:
                flags.append(f"누락시점 {et}분")
    if fit.get("method") in ("failed", "none"):
        flags.append("kel/반감기 fit 실패")

    return BreathResult(
        vco2_mmol_h=float(vco2_mmol_h),
        dose_mmol_13c=float(dose_mmol_13c),
        baseline_delta=base_delta,
        peak_pdr=peak_pdr,
        tmax_h=tmax,
        cpdr_40min=cpdr40,
        cpdr_120min=cpdr120,
        cpdr_total=cpdr_total,
        kel_per_h=float(fit["kel_per_h"]),
        half_life_h=float(fit["half_life_h"]),
        fit_method=str(fit["method"]),
        times_h=list(map(float, t_h)),
        pdr=list(map(float, pdr)),
        flags=flags,
        note="; ".join(note_parts + ([str(fit["note"])] if fit.get("note") else [])),
    )