"""calorimetry.py — 간접열량측정 핵심 계산 로직.

본 도구는 연구용·참고용이며 임상 의사결정에 직접 사용 금지.

단위 규약 (반드시 일관 유지):
- VO2, VCO2 는 함수 인자로 L/min 단위를 받는 것을 표준으로 한다.
- mL/min 데이터는 호출 전에 /1000 으로 변환하여 L/min 으로 넘긴다(importers.py 가 담당).
- Weir 식 표준형:
      REE(kcal/day) = [3.941 * VO2(L/min) + 1.106 * VCO2(L/min)] * 1440  - 2.17 * UN(g/day)
  여기서 mL/min 으로 계산하면 동일 결과를 1000 으로 나눈 형태가 된다:
      REE = [3.941 * VO2(mL/min) + 1.106 * VCO2(mL/min)] * 1440 / 1000 - 2.17 * UN
  (1440 = 분/일)
- 단백 보정(UN)은 24h 요소질소(g/day). 미입력/off 면 0.
"""

from __future__ import annotations

import numpy as np

# Frayn(1983) 기질산화 상수 (단백 보정형, mol→g 환산 포함된 고전적 계수)
# c(carbohydrate oxidation, g/min) = 4.55*VCO2 - 3.21*VO2 - 2.87*n
# f(fat oxidation, g/min)          = 1.67*(VO2 - VCO2) - 1.92*n
# 단, VO2/VCO2 는 L/min, n = 단백 산화율(g/min) = UN(g/day) * 6.25 / 1440
FRAYN_CHO_VCO2 = 4.55
FRAYN_CHO_VO2 = 3.21
FRAYN_CHO_N = 2.87
FRAYN_FAT_DIFF = 1.67
FRAYN_FAT_N = 1.92

# 에너지 밀도 (kcal/g) — 산화 기질 에너지 환산
KCAL_PER_G_CHO = 4.0
KCAL_PER_G_FAT = 9.0
KCAL_PER_G_PRO = 4.0

# 생리적 RQ 범위 (QC)
RQ_PHYS_MIN = 0.67
RQ_PHYS_MAX = 1.00


def urinary_nitrogen_to_protein_ox_per_min(un_g_day: float) -> float:
    """24h 요소질소(g/day) -> 단백 산화율 n (g protein / min).

    단백질 1 g 당 질소 0.16 g (즉 1/6.25). 따라서 단백질(g) = N(g) * 6.25.
    """
    if un_g_day is None or np.isnan(un_g_day):
        return 0.0
    protein_g_day = un_g_day * 6.25
    return protein_g_day / 1440.0


def rq(vo2_lmin, vco2_lmin):
    """호흡지수 RQ = VCO2 / VO2. 0 나눗셈 방지."""
    vo2 = np.asarray(vo2_lmin, dtype=float)
    vco2 = np.asarray(vco2_lmin, dtype=float)
    with np.errstate(divide="ignore", invalid="ignore"):
        out = np.where(vo2 > 0, vco2 / vo2, np.nan)
    return out if out.shape else float(out)


def weir_ree(vo2_lmin, vco2_lmin, un_g_day=0.0, protein_correction=True):
    """Weir(1949) REE (kcal/day). VO2/VCO2 는 L/min.

    REE = [3.941*VO2 + 1.106*VCO2] * 1440 - 2.17*UN(g/day)
    protein_correction=False 이면 UN 항 무시(약식 Weir).
    """
    vo2 = np.asarray(vo2_lmin, dtype=float)
    vco2 = np.asarray(vco2_lmin, dtype=float)
    base = (3.941 * vo2 + 1.106 * vco2) * 1440.0
    if protein_correction and un_g_day:
        base = base - 2.17 * float(un_g_day)
    return base if base.shape else float(base)


def frayn_substrate_oxidation(vo2_lmin, vco2_lmin, un_g_day=0.0, protein_correction=True):
    """Frayn(1983) 탄수화물·지방 산화율.

    반환 dict (스칼라 입력 시 스칼라):
      cho_g_min, fat_g_min, pro_g_min,
      cho_kcal_day, fat_kcal_day, pro_kcal_day,
      np_rq (비단백 RQ)
    VO2/VCO2 L/min.
    """
    vo2 = np.asarray(vo2_lmin, dtype=float)
    vco2 = np.asarray(vco2_lmin, dtype=float)

    if protein_correction:
        n = urinary_nitrogen_to_protein_ox_per_min(un_g_day)
    else:
        n = 0.0

    cho = FRAYN_CHO_VCO2 * vco2 - FRAYN_CHO_VO2 * vo2 - FRAYN_CHO_N * n
    fat = FRAYN_FAT_DIFF * (vo2 - vco2) - FRAYN_FAT_N * n
    pro = np.full_like(vo2, n)  # g/min

    # 비단백 RQ: 단백 산화에 쓰인 가스교환 제거 후 RQ
    # 단백 1 g 당 VO2 ~0.966 L, VCO2 ~0.782 L (고전 값)
    vo2_pro = 0.966 * n
    vco2_pro = 0.782 * n
    npr_vo2 = vo2 - vo2_pro
    npr_vco2 = vco2 - vco2_pro
    with np.errstate(divide="ignore", invalid="ignore"):
        np_rq = np.where(npr_vo2 > 0, npr_vco2 / npr_vo2, np.nan)

    result = {
        "cho_g_min": cho,
        "fat_g_min": fat,
        "pro_g_min": pro,
        "cho_kcal_day": cho * 1440.0 * KCAL_PER_G_CHO,
        "fat_kcal_day": fat * 1440.0 * KCAL_PER_G_FAT,
        "pro_kcal_day": pro * 1440.0 * KCAL_PER_G_PRO,
        "np_rq": np_rq,
    }
    # 스칼라 입력이면 스칼라로 반환
    if vo2.shape == ():
        result = {k: float(v) for k, v in result.items()}
    return result


def qc_flags(vo2_lmin, vco2_lmin, rq_value=None):
    """생리범위 밖 RQ, 음수/결측 등 QC 플래그 리스트 반환(스칼라)."""
    flags = []
    vo2 = float(vo2_lmin) if vo2_lmin is not None else np.nan
    vco2 = float(vco2_lmin) if vco2_lmin is not None else np.nan
    if np.isnan(vo2) or np.isnan(vco2):
        flags.append("MISSING_GAS")
    if not np.isnan(vo2) and vo2 <= 0:
        flags.append("NONPOS_VO2")
    if not np.isnan(vco2) and vco2 <= 0:
        flags.append("NONPOS_VCO2")
    r = rq_value
    if r is None and not np.isnan(vo2) and vo2 > 0:
        r = vco2 / vo2
    if r is not None and not (np.isnan(r)):
        if r < RQ_PHYS_MIN:
            flags.append("RQ_BELOW_PHYS")
        elif r > RQ_PHYS_MAX:
            flags.append("RQ_ABOVE_PHYS")
    return flags


def ee_timeseries_kcal_day(df, vo2_col="VO2_Lmin", vco2_col="VCO2_Lmin",
                           un_g_day=0.0, protein_correction=True):
    """각 시점의 순간 EE(kcal/day) 시계열(ndarray) 반환 (Weir)."""
    vo2 = df[vo2_col].values.astype(float)
    vco2 = df[vco2_col].values.astype(float)
    return weir_ree(vo2, vco2, un_g_day, protein_correction)


def dit_analysis(df, time_col="time_min", vo2_col="VO2_Lmin", vco2_col="VCO2_Lmin",
                 baseline_end_min=30.0, meal_energy_kcal=None,
                 un_g_day=0.0, protein_correction=True):
    """식이유발열생성(DIT) 분석.

    - baseline: time < baseline_end_min 구간 평균 EE(kcal/day) 를 기저로.
    - incremental EE(kcal/day) = 순간 EE - baseline.
    - incremental AUC: 식후(baseline 이후) 구간을 사다리꼴 적분.
        시간은 분 단위, EE 는 kcal/day 이므로 kcal 환산: EE(kcal/day) * (dt_min/1440).
    - meal_energy_kcal 주어지면 % thermogenesis = DIT_kcal / meal * 100.

    반환 dict: baseline_ee_kcal_day, ee_series(kcal/day ndarray),
               incr_ee(kcal/day ndarray), dit_auc_kcal, pct_thermogenesis(or None),
               peak_incr_kcal_day, peak_time_min
    """
    t = df[time_col].values.astype(float)
    ee = ee_timeseries_kcal_day(df, vo2_col, vco2_col, un_g_day, protein_correction)

    base_mask = t < baseline_end_min
    baseline = float(np.nanmean(ee[base_mask])) if base_mask.any() else float(np.nanmean(ee))
    incr = ee - baseline

    post_mask = t >= baseline_end_min
    tp = t[post_mask]
    incp = incr[post_mask]
    if tp.size >= 2:
        # 사다리꼴: kcal = sum( mean_incr_EE(kcal/day) * dt_min / 1440 )
        dt = np.diff(tp)
        seg = (incp[:-1] + incp[1:]) / 2.0
        dit_auc_kcal = float(np.nansum(seg * dt / 1440.0))
    else:
        dit_auc_kcal = float("nan")

    pct = None
    if meal_energy_kcal:
        pct = float(dit_auc_kcal / meal_energy_kcal * 100.0)

    if incp.size:
        pk_idx = int(np.nanargmax(incp))
        peak_incr = float(incp[pk_idx])
        peak_time = float(tp[pk_idx])
    else:
        peak_incr = float("nan")
        peak_time = float("nan")

    return {
        "baseline_ee_kcal_day": baseline,
        "ee_series": ee,
        "incr_ee": incr,
        "dit_auc_kcal": dit_auc_kcal,
        "pct_thermogenesis": pct,
        "peak_incr_kcal_day": peak_incr,
        "peak_time_min": peak_time,
    }


def summarize_steady(df, vo2_col="VO2_Lmin", vco2_col="VCO2_Lmin",
                     un_g_day=0.0, protein_correction=True):
    """안정구간 DataFrame 평균에 대해 핵심 지표를 산출한 dict 반환."""
    vo2_mean = float(np.nanmean(df[vo2_col].values))
    vco2_mean = float(np.nanmean(df[vco2_col].values))
    r = float(vco2_mean / vo2_mean) if vo2_mean > 0 else float("nan")
    ree = weir_ree(vo2_mean, vco2_mean, un_g_day, protein_correction)
    frayn = frayn_substrate_oxidation(vo2_mean, vco2_mean, un_g_day, protein_correction)
    flags = qc_flags(vo2_mean, vco2_mean, r)
    out = {
        "vo2_lmin": vo2_mean,
        "vco2_lmin": vco2_mean,
        "rq": r,
        "ree_kcal_day": float(ree),
        "n_samples": int(len(df)),
        "qc_flags": flags,
    }
    out.update({k: (float(v) if np.isscalar(v) or getattr(v, "shape", None) == () else v)
                for k, v in frayn.items()})
    return out