"""
deconvolution.py — C-peptide 디컨볼루션 기반 인슐린 분비율(ISR) 추정

방법: Van Cauter E, Mestrez F, Sturis J, Polonsky KS.
"Estimation of insulin secretion rates from C-peptide levels: comparison of
individual and standard kinetic parameters for C-peptide clearance."
Diabetes 1992;41:368-377.

C-peptide 동태는 2-구획 모델로 기술되며, 표준 kinetic 파라미터는 연령·성별·BSA·
비만/당뇨 여부로 결정된다. 본 구현은 해당 논문 Table에 보고된 표준 파라미터
회귀식을 사용한다:

  short half-life T1/2(short) = 4.95 분 (모든 대상 공통, 근사)
  long  half-life T1/2(long):
      - 정상/비만:  T1/2(long) = 0.14*age + 29.2   (분)
      - NIDDM:      T1/2(long) = 0.transformed... (본 구현은 정상/비만 식 사용 + 옵션)
  fraction(short) F = 0.76 (근사, 논문 평균)

2-구획 모델 파라미터 변환:
  a = ln(2)/T1/2(short)   (빠른 소실 rate, /min)
  b = ln(2)/T1/2(long)    (느린 소실 rate, /min)
  F = fraction of short component
  분포용적 V = 보고된 표준식 (BSA 기반)

디컨볼루션: 측정 C-peptide 시계열 C(t)와 2-지수 임펄스응답으로부터
ISR(t)를 산출. 본 구현은 Van Cauter 표준 2-구획 질량보존식의
이산(piecewise) 형태를 사용:
  미분형:  dC/dt = ISR/V - (k01+k12)*C + k21*Y   (Y: 말초구획)
간단·견고함을 위해 well-validated한 "2-exponential impulse response +
비음수 최소제곱 디컨볼루션" 방식을 구현하되, scipy 미설치 시
근사 디컨볼루션(매 구간 정상상태 가정)으로 graceful fallback 한다.

ISR 단위: pmol/min (V를 L, C를 pmol/L로 두면 ISR=pmol/min).
입력 C-peptide가 ng/mL이면 units 모듈로 nmol/L 변환 후 pmol/L(*1000) 사용 권장.
"""
from __future__ import annotations
import math
import numpy as np

try:
    from scipy.optimize import nnls  # type: ignore
    _HAVE_SCIPY = True
except Exception:  # pragma: no cover - scipy 미설치 환경
    _HAVE_SCIPY = False


# ---------------------------------------------------------------------------
# 표준 kinetic 파라미터 (Van Cauter 1992)
# ---------------------------------------------------------------------------
def bsa_dubois(weight_kg: float, height_cm: float) -> float:
    """Du Bois & Du Bois 체표면적 (m^2)."""
    return 0.007184 * (weight_kg ** 0.425) * (height_cm ** 0.725)


def standard_kinetics(age_years: float,
                      bsa_m2: float,
                      diabetic: bool = False,
                      obese: bool = False):
    """
    Van Cauter 1992 표준 C-peptide kinetic 파라미터 반환.

    반환 dict:
      F            : short component fraction (무차원)
      t_short_min  : short half-life (분)
      t_long_min   : long half-life (분)
      a, b         : 두 지수 소실 rate (/min)
      V_L          : 분포용적 (L)
    """
    # short half-life: 고정 근사
    t_short = 4.95
    # long half-life: 그룹별 회귀 (Van Cauter 1992, Table 2 근사식)
    if diabetic:
        t_long = 0.14 * age_years + 29.2   # NIDDM 군 근사
    elif obese:
        t_long = 0.11 * age_years + 30.0   # 비만 군 근사
    else:
        t_long = 0.14 * age_years + 29.2   # 정상 군

    F = 0.76  # 평균 fraction(short)

    a = math.log(2) / t_short
    b = math.log(2) / t_long

    # 분포용적: V = 1.92 * BSA + 0.64  (L) 근사 (논문 보고 범위 내)
    V_L = 1.92 * bsa_m2 + 0.64

    return {
        "F": F,
        "t_short_min": t_short,
        "t_long_min": t_long,
        "a": a,
        "b": b,
        "V_L": V_L,
    }


# ---------------------------------------------------------------------------
# 임펄스 응답 (2-exponential)
# ---------------------------------------------------------------------------
def impulse_response(t, params):
    """
    단위 임펄스(1 pmol 순간 분비)에 대한 plasma C-peptide 농도 응답 (pmol/L).
    h(t) = (1/V) * [F*exp(-a t) + (1-F)*exp(-b t)]
    """
    F = params["F"]; a = params["a"]; b = params["b"]; V = params["V_L"]
    t = np.asarray(t, dtype=float)
    return (1.0 / V) * (F * np.exp(-a * t) + (1.0 - F) * np.exp(-b * t))


# ---------------------------------------------------------------------------
# ISR 디컨볼루션
# ---------------------------------------------------------------------------
def deconvolve_isr(times_min, cpep_pmol_L, params, dt=1.0):
    """
    측정 C-peptide(pmol/L) 시계열로부터 ISR(pmol/min) 산출.

    방법:
      1) 측정점을 dt(기본 1분) 그리드로 선형보간
      2) 임펄스응답 h로 컨볼루션 행렬 H 구성 (이산, dt 가중)
      3) C ≈ H @ ISR 를 비음수 최소제곱(nnls)으로 풀이
         (scipy 없으면 정상상태 근사로 fallback)

    반환 dict:
      grid_min   : 분 그리드
      isr_pmol_min: ISR 시계열 (pmol/min)
      isr_basal  : 기저 ISR (첫 구간 평균)
      isr_total_auc: ISR AUC (pmol) = 총 분비량 추정
      method     : 'nnls' | 'steadystate_fallback'
    """
    times_min = np.asarray(times_min, dtype=float)
    cpep = np.asarray(cpep_pmol_L, dtype=float)
    ok = ~(np.isnan(times_min) | np.isnan(cpep))
    times_min, cpep = times_min[ok], cpep[ok]
    if times_min.size < 2:
        return {
            "grid_min": np.array([]), "isr_pmol_min": np.array([]),
            "isr_basal": float("nan"), "isr_total_auc": float("nan"),
            "method": "insufficient_data",
        }

    order = np.argsort(times_min)
    times_min, cpep = times_min[order], cpep[order]

    t0, t1 = float(times_min[0]), float(times_min[-1])
    grid = np.arange(t0, t1 + dt, dt)
    Cg = np.interp(grid, times_min, cpep)
    n = len(grid)

    a = params["a"]; b = params["b"]; F = params["F"]; V = params["V_L"]

    use_nnls = _HAVE_SCIPY and n >= 3
    if use_nnls:
        try:
            # 컨볼루션 행렬: H[i,j] = h((i-j)*dt)*dt  for i>=j else 0
            tau = np.arange(n) * dt
            h = impulse_response(tau, params)  # pmol/L per pmol-impulse
            H = np.zeros((n, n))
            for j in range(n):
                length = n - j
                H[j:, j] = h[:length] * dt
            # 열 스케일링으로 조건수 개선: 각 열을 그 L2-norm으로 정규화 후 복원
            col_norm = np.linalg.norm(H, axis=0)
            col_norm[col_norm == 0] = 1.0
            Hs = H / col_norm
            # Tikhonov 1차 차분 평활 정규화 (스케일된 공간에서)
            lam = 5e-2
            D = np.zeros((n - 1, n))
            for i in range(n - 1):
                D[i, i] = -1.0; D[i, i + 1] = 1.0
            A = np.vstack([Hs, lam * D])
            rhs = np.concatenate([Cg, np.zeros(n - 1)])
            with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
                isr_s, _ = nnls(A, rhs, maxiter=10 * n)
            isr = isr_s / col_norm  # 원래 스케일 복원
            if not np.all(np.isfinite(isr)):
                raise FloatingPointError("nnls non-finite")
            method = "nnls"
        except Exception:
            use_nnls = False
    if not use_nnls:
        # Fallback: 각 시점에서 정상상태 가정 C ≈ ISR * (F/a + (1-F)/b)/V
        gain = (F / a + (1.0 - F) / b) / V  # pmol/L per (pmol/min)
        isr = Cg / gain if gain > 0 else np.full(n, np.nan)
        method = "steadystate_fallback"

    # 기저 ISR: 처음 5분(또는 첫 점) 평균
    basal_mask = grid <= (t0 + 5)
    isr_basal = float(np.nanmean(isr[basal_mask])) if basal_mask.any() else float(isr[0])
    isr_total_auc = float(np.trapezoid(isr, grid))

    return {
        "grid_min": grid,
        "isr_pmol_min": isr,
        "isr_basal": isr_basal,
        "isr_total_auc": isr_total_auc,
        "method": method,
    }