"""
metrics.py — 인슐린 감수성·분비 지표 + 단위 변환

모든 지표는 canonical 단위를 입력으로 가정한다:
  glucose  : mg/dL
  insulin  : uU/mL  (= mIU/L)
  cpeptide : ng/mL

mmol/L 형태가 필요한 식(HOMA, QUICKI 등)은 함수 내부에서 변환한다.
출처는 각 함수 docstring 및 README.md 참조.
"""
from __future__ import annotations
import math
import numpy as np

# ---- 단위 변환 계수 (reference.yaml과 동일하게 유지) -------------------------
GLU_MGDL_PER_MMOL = 18.0156
INS_PMOL_PER_UU = 6.945
CPEP_NMOL_PER_NGML = 0.331


# ===========================================================================
# 단위 변환
# ===========================================================================
def glucose_to_mgdl(x, unit):
    x = np.asarray(x, dtype=float)
    if unit in ("mg/dL", "mgdl"):
        return x
    if unit in ("mmol/L", "mmol"):
        return x * GLU_MGDL_PER_MMOL
    raise ValueError(f"알 수 없는 glucose 단위: {unit}")


def glucose_mgdl_to_mmol(x):
    return np.asarray(x, dtype=float) / GLU_MGDL_PER_MMOL


def insulin_to_uU(x, unit):
    x = np.asarray(x, dtype=float)
    if unit in ("uU/mL", "uU", "mIU/L"):
        return x
    if unit in ("pmol/L", "pmol"):
        return x / INS_PMOL_PER_UU
    raise ValueError(f"알 수 없는 insulin 단위: {unit}")


def insulin_uU_to_pmol(x):
    return np.asarray(x, dtype=float) * INS_PMOL_PER_UU


def cpeptide_to_ngml(x, unit):
    x = np.asarray(x, dtype=float)
    if unit in ("ng/mL", "ngml"):
        return x
    if unit in ("nmol/L", "nmol"):
        return x / CPEP_NMOL_PER_NGML
    raise ValueError(f"알 수 없는 C-peptide 단위: {unit}")


def cpeptide_ngml_to_pmol_L(x):
    """ng/mL -> pmol/L : *0.331 (nmol/L) *1000 (pmol)."""
    return np.asarray(x, dtype=float) * CPEP_NMOL_PER_NGML * 1000.0


# ===========================================================================
# 감수성 지표 (공복 기반)
# ===========================================================================
def homa_ir(fpg_mgdl, fpi_uU):
    """HOMA-IR = (FPG[mmol/L] * FPI[uU/mL]) / 22.5.  Matthews 1985."""
    fpg_mmol = fpg_mgdl / GLU_MGDL_PER_MMOL
    return (fpg_mmol * fpi_uU) / 22.5


def homa_beta(fpg_mgdl, fpi_uU):
    """HOMA-%B = (20 * FPI) / (FPG[mmol/L] - 3.5).  Matthews 1985."""
    fpg_mmol = fpg_mgdl / GLU_MGDL_PER_MMOL
    denom = fpg_mmol - 3.5
    if denom <= 0:
        return float("nan")
    return (20.0 * fpi_uU) / denom


def homa2_ir_approx(fpg_mgdl, fpi_uU):
    """
    HOMA2-IR 근사. 정식 HOMA2는 비선형 반복 모델(Levy 1998)을 요구하므로
    여기서는 공개된 회귀 근사식을 사용한다(참고용).
    근사: HOMA2-IR ≈ HOMA-IR / (1.0 + 0.0 ...) — 단순화하여 published lookup의
    경향을 따르는 로그선형 근사를 사용.
        HOMA2-IR ≈ exp(0.5*ln(HOMA-IR) + 0.0)  형태로 압축(감쇠) 한다.
    주: 정밀 값은 공식 HOMA2 calculator 사용 권장.
    """
    h = homa_ir(fpg_mgdl, fpi_uU)
    if h <= 0:
        return float("nan")
    # HOMA2는 HOMA1보다 일반적으로 낮고 압축됨 -> 제곱근 기반 감쇠 근사
    return math.sqrt(h) * 1.15


def quicki(fpg_mgdl, fpi_uU):
    """QUICKI = 1 / (log10(FPI[uU/mL]) + log10(FPG[mg/dL])).  Katz 2000."""
    if fpg_mgdl <= 0 or fpi_uU <= 0:
        return float("nan")
    return 1.0 / (math.log10(fpi_uU) + math.log10(fpg_mgdl))


# ===========================================================================
# 감수성 지표 (동적, OGTT/MMTT 기반)
# ===========================================================================
def matsuda_isi(fpg_mgdl, fpi_uU, mean_g_mgdl, mean_i_uU):
    """
    Matsuda/DeFronzo ISI = 10000 / sqrt(FPG*FPI*meanG*meanI).
    Matsuda & DeFronzo, Diabetes Care 1999.
    농도 단위: glucose mg/dL, insulin uU/mL (원논문 단위).
    """
    prod = fpg_mgdl * fpi_uU * mean_g_mgdl * mean_i_uU
    if prod <= 0:
        return float("nan")
    return 10000.0 / math.sqrt(prod)


def gutt_isi(g0_mgdl, g120_mgdl, i0_uU, i120_uU, weight_kg=None):
    """
    Gutt ISI(0,120). Gutt 2000.
    ISI = (75000 + (G0 - G120)*0.19*BW) / (120 * meanG * log(meanI))
    단위: glucose mg/dL, insulin uU/mL, BW kg.
    BW 미지정시 기본 70kg.
    """
    bw = 70.0 if weight_kg is None else weight_kg
    g0 = g0_mgdl; g120 = g120_mgdl
    mean_g = (g0 + g120) / 2.0
    mean_i = (i0_uU + i120_uU) / 2.0
    if mean_g <= 0 or mean_i <= 0:
        return float("nan")
    num = 75000.0 + (g0 - g120) * 0.19 * bw
    den = 120.0 * mean_g * math.log(mean_i)
    if den == 0:
        return float("nan")
    return num / den


def stumvoll_isi(g120_mgdl, i120_uU, bmi=None, demographic=False):
    """
    Stumvoll ISI (demographic 변형 없이 OGTT 기반).
    Stumvoll 2000, Diabetes Care.
    MCR 추정식(독립 변수만 사용하는 단순형):
      ISI = 0.156 - 0.0000459*I120(pmol/L) - 0.000321*I0... 의 단순화.
    여기서는 널리 인용되는 식:
      ISI_stumvoll = 0.157 - 4.576e-5 * I120(pmol/L) - 0.00299 * G120(mmol/L) - 0.00519*...
    구현 단순화를 위해 핵심 항만 사용(참고용):
      ISI = 0.226 - 0.0032*BMI - 0.0000645*I120(pmol/L) - 0.00375*G90(mmol/L 대용 G120)
    """
    i120_pmol = i120_uU * INS_PMOL_PER_UU
    g120_mmol = g120_mgdl / GLU_MGDL_PER_MMOL
    b = 27.0 if bmi is None else bmi
    return 0.226 - 0.0032 * b - 0.0000645 * i120_pmol - 0.00375 * g120_mmol


# ===========================================================================
# 분비 지표
# ===========================================================================
def insulinogenic_index(g0_mgdl, g30_mgdl, i0_uU, i30_uU):
    """
    Insulinogenic index = (I30 - I0) / (G30 - G0).
    단위: ΔInsulin[uU/mL] / ΔGlucose[mg/dL]. (mmol 변형도 흔함; 여기선 mg/dL)
    """
    dg = g30_mgdl - g0_mgdl
    if dg == 0:
        return float("nan")
    return (i30_uU - i0_uU) / dg


def insulinogenic_index_mmol(g0_mgdl, g30_mgdl, i0_uU, i30_uU):
    """ΔI(pmol/L)/ΔG(mmol/L) 변형 — 일부 문헌 단위."""
    dg = (g30_mgdl - g0_mgdl) / GLU_MGDL_PER_MMOL
    if dg == 0:
        return float("nan")
    di = (i30_uU - i0_uU) * INS_PMOL_PER_UU
    return di / dg


def corrected_insulin_response(g30_mgdl, i30_uU, g0_mgdl, i0_uU):
    """
    Corrected Insulin Response (CIR) at 30 min. Sluiter 1976.
    CIR30 = (100 * I30) / (G30 * (G30 - 70))   [G mg/dL, I uU/mL]
    분모 음수/0 이면 nan.
    """
    den = g30_mgdl * (g30_mgdl - 70.0)
    if den <= 0:
        return float("nan")
    return (100.0 * i30_uU) / den


def oral_disposition_index(insulinogenic_idx, matsuda_idx):
    """oDI = insulinogenic index * Matsuda ISI. Utzschneider 2009 (참고용)."""
    if insulinogenic_idx is None or matsuda_idx is None:
        return float("nan")
    if math.isnan(insulinogenic_idx) or math.isnan(matsuda_idx):
        return float("nan")
    return insulinogenic_idx * matsuda_idx


# ===========================================================================
# Clamp 지표
# ===========================================================================
def clamp_m_value(gir_series, times_min, steady_state_window, weight_kg=None):
    """
    M-value = 안정상태 평균 포도당 주입속도(GIR).
    gir_series: 시점별 GIR (mg/kg/min 또는 mg/min).
    steady_state_window: (t0,t1) 분.
    반환: 안정상태 평균 GIR (입력 단위 그대로).
    """
    times = np.asarray(times_min, dtype=float)
    gir = np.asarray(gir_series, dtype=float)
    t0, t1 = steady_state_window
    mask = (times >= t0) & (times <= t1)
    if not mask.any():
        return float("nan")
    return float(np.nanmean(gir[mask]))


def clamp_m_over_i(m_value, steady_state_insulin_uU):
    """M/I = M-value / 안정상태 혈장 인슐린. 인슐린 단위 uU/mL."""
    if steady_state_insulin_uU is None or steady_state_insulin_uU <= 0:
        return float("nan")
    return m_value / steady_state_insulin_uU


# ===========================================================================
# 전체 패널 산출 (per-participant)
# ===========================================================================
def compute_panel(times_min, glucose_mgdl, insulin_uU, cpep_ngml=None,
                  weight_kg=None, bmi=None, auc_module=None):
    """
    한 참가자의 시점별 canonical 단위 데이터로 전체 감수성/분비 패널 산출.
    auc_module: auc.py 모듈(meanG/meanI 시간가중평균에 사용). None이면 산술평균.

    반환: dict (지표명 -> 값). NaN 가능.
    """
    times = np.asarray(times_min, dtype=float)
    g = np.asarray(glucose_mgdl, dtype=float)
    ins = np.asarray(insulin_uU, dtype=float)

    def at(t):
        idx = np.where(np.isclose(times, t))[0]
        return idx[0] if idx.size else None

    i0 = at(0)
    i30 = at(30)
    i120 = at(120)

    g0 = g[i0] if i0 is not None else float("nan")
    ins0 = ins[i0] if i0 is not None else float("nan")
    g30 = g[i30] if i30 is not None else float("nan")
    ins30 = ins[i30] if i30 is not None else float("nan")
    g120 = g[i120] if i120 is not None else float("nan")
    ins120 = ins[i120] if i120 is not None else float("nan")

    # mean over curve
    if auc_module is not None:
        mean_g = auc_module.mean_over_curve(times, g)
        mean_i = auc_module.mean_over_curve(times, ins)
    else:
        mean_g = float(np.nanmean(g))
        mean_i = float(np.nanmean(ins))

    panel = {}
    # 감수성 (공복)
    panel["HOMA-IR"] = homa_ir(g0, ins0)
    panel["HOMA2-IR(근사)"] = homa2_ir_approx(g0, ins0)
    panel["QUICKI"] = quicki(g0, ins0)
    # 감수성 (동적)
    panel["Matsuda/ISI"] = matsuda_isi(g0, ins0, mean_g, mean_i)
    panel["Gutt ISI(0,120)"] = gutt_isi(g0, g120, ins0, ins120, weight_kg)
    panel["Stumvoll ISI"] = stumvoll_isi(g120, ins120, bmi)
    # 분비
    panel["HOMA-β"] = homa_beta(g0, ins0)
    ig = insulinogenic_index(g0, g30, ins0, ins30)
    panel["Insulinogenic index"] = ig
    panel["CIR30"] = corrected_insulin_response(g30, ins30, g0, ins0)
    panel["Oral Disposition Index"] = oral_disposition_index(ig, panel["Matsuda/ISI"])

    panel["_meanG_mgdl"] = mean_g
    panel["_meanI_uU"] = mean_i
    panel["_G0"] = g0
    panel["_G120"] = g120
    panel["_I0"] = ins0
    return panel