"""
drug_comsol_params.py — 약물별 COMSOL 입력 파라미터 (v2, 물리공식 기반)
============================================================
약물 선택 시 해당 약물의 10개 COMSOL 파라미터를 반환한다.

★ v2 변경: BA 역산(순환 논리) 폐기 → 물리 공식·문헌 기반 재유도
  기존 v1은 k_decay를 "문헌 BA가 나오도록 역산"해 넣어, 대리모델→BA
  계산이 순환(BA를 넣어 BA를 얻음)이었다. v2는 이를 제거하고 각
  파라미터를 물리적 근거로 독립 결정한다. 그 결과 BA는 순방향 예측값이 된다.

파라미터별 결정 방법:
  Lp_ve, K          조직 성질(약물 무관) → 공통 상수 (Lp_ve=8e-12, K=1e-15)
  sigma_ve, sigma_le  Curry 반사계수 σ=1-(1-(a/R)²)², a=0.0515·MW^0.392 nm
                      (Erickson 2009). 벽 공극반경 R_ve=6.22, R_le=24.5 nm는
                      파일 3종(IgG/INS/ALB) 실측 best-fit.
  D_gel             Stokes-Einstein 자유확산 D=kBT/(6πηr) × ECM 방해인자
                      f=0.376·exp(0.1439·r_h[nm]) (파일 3종 실측 피팅, 37°C)
  p_ve, p_le        벽 확산투과도. 파일 3종 실측을 log(MW) 선형 보간.
  kf_m, kr_m        ECM 결합 상수 → 파일(250415_Drug_list) 실측 그대로.
  k_decay           ★ 목표 BA 재현값(calibration). 알려진 약물은 문헌 SC BA를
                      재현하도록 역산된 값이다(순방향 '예측'이 아니라 '재현').
                      LODO 검증상 물성으로 예측 불가(27.9%p)이므로 미지 약물엔
                      k_decay_physics(MW)(MW 로그스케일 물리추정)를 쓸 것.
                      물리값 원본은 K_DECAY_PHYSICS dict에 보존.
                      TRA는 추천모델5(흡수 왜곡)라 model4 기준으로 재현함.

파일 실측 3종(IgG/INS/ALB, 251103_revised)은 sigma·p·D_gel을 실측 그대로 사용.

────────────────────────────────────────────────────────────
★ 파라미터 근거 등급 (생체물리 문헌 감사, 상세: comsol_params_biophysics_refs_ko.md)
  주의: 이 10개는 임상 PK 상수(CL·V·BA)가 아니라 조직-분자 수송 생체물리량이라
  FDA/PopPK DB에서 약물별로 찾을 수 없다. 미세혈관 생리·ECM 확산·공극이론으로 유도.

  🟩 물리공식/문헌실측 (5개):
    sigma_ve, sigma_le : Curry-Michel 공극이론 + Rippe & Haraldsson two-pore
                         (two-pore eqs; 소공극~4nm/대공극~24nm). 파일 3종 재현.
    D_gel              : Stokes-Einstein + ECM 방해(Pluen 2001 PNAS, Ramanujan 2002 BiophysJ).
    p_ve, p_le         : 미세혈관 단백질 투과도 실측과 자릿수 일치
                         (알부민 1.8e-9, IgG 1.9~4.6e-9 m/s; Yuan 1993 PMID:8321142, PMC2363402).
  ⬜ 조직상수 고정 (약물 무관, 2개):
    Lp_ve              : 사람 미세혈관 수리전도도 범위 고정. Levick & Michel 2010.
    K                  : 사람 간질 Darcy 투과도 범위 고정. Levick 1987, Q J Exp Physiol 72:409.
  🟧 기전은 문헌확인, 절대 속도상수 미검증 (3개 — 파일 실측 의존):
    kf_m, kr_m         : ECM-하이알루로난 정전기 결합은 실재 기전(J Pharm Sci 2020)이나
                         약물별 kon/koff 표준값 없음 → 파일(250415) 실측 사용.
    k_decay            : SC presystemic 분해. 기전·방향(펩타이드↑ mAb↓)·BA범위(52~80%)는
                         문헌확인(PMC5613179; FcRn PBPK antibodies 2024 doi:10.3390/antib13030070;
                         FcRn결손 BA 82.5%→28.3% PMC6150636)이나 약물별 속도상수 실측 없음
                         → MW 대리스케일. BA 역산 아님(순환 회피).
────────────────────────────────────────────────────────────

파라미터 순서:
  Lp_ve, K, p_le, sigma_ve, sigma_le, p_ve, D_gel, k_decay, kf_m, kr_m
"""

import math

# 약물 무관 조직 공통값
_COMMON = dict(Lp_ve=8e-12, K=1e-15)

# 약물별 COMSOL 파라미터
#   IgG/INS/ALB: 파일 실측(sigma·p·D_gel) + 문헌 k_decay
#   그 외 7종  : 물리공식(sigma·D_gel·p) + 파일 kf/kr + 문헌 k_decay
# ── k_decay 물리값(MW 스케일, BA 미보정) — 미지 약물 폴백·참조용 보존 ──
K_DECAY_PHYSICS = {
    "INS": 3.000e-05,
    "EXE": 3.925e-05,
    "ETA": 2.000e-06,
    "TRA": 2.051e-06,
    "RIT": 2.070e-06,
    "IgG": 2.000e-06,
    "BOC": 2.056e-06,
}
# 위 DRUG_COMSOL의 k_decay는 아래 목표 BA를 재현하도록 보정된 값(calibration)이다.
# ⚠ 이는 순방향 '예측'이 아니라 알려진 문헌 BA의 '재현'값이다. LODO 검증 결과
#   물성(MW·sigma_ve)으로 예측 불가(평균오차 27.9%p)하므로, 미지 약물에는
#   k_decay_physics()(MW 로그스케일)를 써야 한다. 상세: pipeline_validation_report_ko.md
CALIBRATION_BA = {"INS": 66, "EXE": 100, "ETA": 58, "TRA": 77, "RIT": 64, "IgG": 70, "BOC": 46.5}

DRUG_COMSOL = {
    "INS": dict(MW=5.80, sigma_ve=0.1000, sigma_le=0.0000, p_ve=1.000e-08, p_le=7.200e-08, D_gel=9.400e-11, k_decay=1.494e-05, kf_m=1.450, kr_m=6.82, **_COMMON),
    "ALB": dict(MW=66.50, sigma_ve=0.7000, sigma_le=0.0500, p_ve=1.000e-09, p_le=1.000e-08, D_gel=6.200e-11, k_decay=3.937e-06, kf_m=0.170, kr_m=4.77, **_COMMON),
    "IgG": dict(MW=150.00, sigma_ve=0.9000, sigma_le=0.1000, p_ve=1.000e-09, p_le=1.000e-08, D_gel=4.500e-11, k_decay=4.102e-06, kf_m=0.480, kr_m=4.20, **_COMMON),
    "EXE": dict(MW=4.20, sigma_ve=0.0927, sigma_le=0.0061, p_ve=1.171e-08, p_le=8.241e-08, D_gel=1.093e-10, k_decay=1.361e-06, kf_m=0.624, kr_m=11.14, **_COMMON),
    "LYS": dict(MW=14.30, sigma_ve=0.2328, sigma_le=0.0159, p_ve=4.600e-09, p_le=3.700e-08, D_gel=7.628e-11, k_decay=1.415e-05, kf_m=1.451, kr_m=99.70, **_COMMON),
    "OVA": dict(MW=45.00, sigma_ve=0.5168, sigma_le=0.0389, p_ve=1.920e-09, p_le=1.749e-08, D_gel=5.821e-11, k_decay=5.449e-06, kf_m=0.371, kr_m=6.36, **_COMMON),
    "RIT": dict(MW=143.90, sigma_ve=0.9417, sigma_le=0.0954, p_ve=7.913e-10, p_le=8.182e-09, D_gel=4.909e-11, k_decay=6.238e-06, kf_m=0.296, kr_m=6.35, **_COMMON),
    "BOC": dict(MW=145.10, sigma_ve=0.9440, sigma_le=0.0960, p_ve=7.863e-10, p_le=8.137e-09, D_gel=4.905e-11, k_decay=1.140e-05, kf_m=0.311, kr_m=5.99, **_COMMON),
    "TRA": dict(MW=145.50, sigma_ve=0.9448, sigma_le=0.0962, p_ve=7.846e-10, p_le=8.123e-09, D_gel=4.904e-11, k_decay=2.666e-06, kf_m=0.517, kr_m=8.77, **_COMMON),
    "ETA": dict(MW=149.99, sigma_ve=0.9531, sigma_le=0.0985, p_ve=7.667e-10, p_le=7.963e-09, D_gel=4.892e-11, k_decay=6.981e-06, kf_m=0.756, kr_m=13.40, **_COMMON),
}

# 파일에서 직접 측정된 약물 (실측 sigma·p·D_gel, 신뢰도 높음)
FILE_MEASURED = {"IgG", "INS", "ALB"}

PARAM_ORDER = ["Lp_ve", "K", "p_le", "sigma_ve", "sigma_le",
               "p_ve", "D_gel", "k_decay", "kf_m", "kr_m"]


def get_comsol_params(drug="IgG"):
    if drug not in DRUG_COMSOL:
        raise KeyError(f"약물 '{drug}' 없음. 지원: {list(DRUG_COMSOL)}")
    return DRUG_COMSOL[drug]


def get_param_vector(drug="IgG"):
    p = get_comsol_params(drug)
    return [p[name] for name in PARAM_ORDER]


def is_file_measured(drug):
    """파일에서 직접 측정된 약물인지 (vs 물리공식 계산)."""
    return drug in FILE_MEASURED


def k_decay_physics(mw_kda):
    """미지 약물용 k_decay (MW 로그스케일 물리 추정, 1/s).
    BA calibration이 없는 약물에 사용. 펩타이드(~5kDa)~3e-5, mAb(~150kDa)~2e-6.
    ⚠ 물리 순방향 추정값 — BA 재현 보장 안 됨."""
    import math
    return 10**(-0.8326*math.log10(mw_kda*1000.0) + -1.3896)


def get_k_decay(drug, mw_kda=None):
    """약물 k_decay 반환. 알려진 약물은 calibration(목표 BA 재현)값,
    미지 약물은 k_decay_physics(MW) 물리 추정값."""
    if drug in DRUG_COMSOL:
        return DRUG_COMSOL[drug]["k_decay"]
    if mw_kda is None:
        raise ValueError("미지 약물은 mw_kda 필요")
    return k_decay_physics(mw_kda)


# ── 미지 약물: 분자량(MW)만으로 10개 파라미터 추정 ───────────────────────
# 검증(2026-07): sigma·D_gel 공식은 물리유도 7종을 완벽 재현, p_ve/p_le는
# 자릿수 일치. kf_m·kr_m은 MW와 무관하므로 MW 최근접 등록약물에서 차용한다.
_MW_R_VE_NM = 6.22     # 미세혈관 소공극 반경 (Rippe & Haraldsson 1994)
_MW_R_LE_NM = 24.5     # 림프관 대공극 반경 (파일 3종 실측 best-fit)
_MW_KB = 1.380649e-23  # Boltzmann (J/K)
_MW_T  = 310.15        # 37°C (K)
_MW_ETA = 0.0007       # ECM 점도 (Pa·s)

def _solute_radius_nm(mw_da):
    """용질 Stokes 반경 (nm), Erickson 2009: a = 0.0515·MW_Da^0.392."""
    return 0.0515 * mw_da ** 0.392

def _sigma_curry(mw_da, R_nm):
    """Curry-Michel 반사계수 σ = 1−(1−(a/R)²)²."""
    r = _solute_radius_nm(mw_da) / R_nm
    return 1.0 - (1.0 - r * r) ** 2 if r < 1.0 else 1.0

def _d_gel_stokes(mw_da):
    """Stokes-Einstein 자유확산 × ECM 방해인자 (m²/s)."""
    r_nm = _solute_radius_nm(mw_da)
    D = _MW_KB * _MW_T / (6.0 * math.pi * _MW_ETA * r_nm * 1e-9)
    f = 0.376 * math.exp(0.144 * r_nm)      # ECM 방해(감소)인자
    return D * f

def _p_interp_from_anchors(mw_da, key):
    """IgG/INS/ALB 실측 3점을 log(MW)–log(p) 선형보간 (외삽은 양끝 고정)."""
    anc = sorted((DRUG_COMSOL[d]["MW"] * 1000.0, DRUG_COMSOL[d][key])
                 for d in ("INS", "ALB", "IgG"))
    xs = [math.log10(m) for m, _ in anc]
    ys = [math.log10(p) for _, p in anc]
    x = math.log10(mw_da)
    if x <= xs[0]:   return 10 ** ys[0]
    if x >= xs[-1]:  return 10 ** ys[-1]
    for i in range(len(xs) - 1):
        if xs[i] <= x <= xs[i + 1]:
            t = (x - xs[i]) / (xs[i + 1] - xs[i])
            return 10 ** (ys[i] + t * (ys[i + 1] - ys[i]))
    return 10 ** ys[-1]

def _nearest_mw_drug(mw_kda):
    """MW(kDa)가 가장 가까운 등록 약물 이름."""
    return min(DRUG_COMSOL, key=lambda d: abs(DRUG_COMSOL[d]["MW"] - mw_kda))

def params_from_mw(mw_kda, kf_m=None, kr_m=None, return_detail=False):
    """분자량(kDa)만으로 미지 약물의 COMSOL 10개 파라미터를 추정한다.

    추정 방식:
      Lp_ve, K          : 조직 상수 (_COMMON, 약물 무관)
      sigma_ve/sigma_le : Curry-Michel 공극이론 (MW → Stokes 반경 → σ). 검증 완료.
      D_gel             : Stokes-Einstein × ECM 방해인자. 검증 완료.
      p_ve, p_le        : IgG/INS/ALB 실측 log(MW) 보간 (자릿수 신뢰).
      k_decay           : k_decay_physics(MW) 로그스케일 물리추정 (근사, BA 미보정).
      kf_m, kr_m        : MW로 예측 불가 → 인자로 주면 그 값, 없으면 MW 최근접
                          등록약물의 값을 차용.

    인자:
      mw_kda      : 분자량 (kDa)
      kf_m, kr_m  : 알면 직접 지정 (None이면 최근접 약물에서 차용)
      return_detail : True면 (vector, detail dict) 반환
    반환:
      PARAM_ORDER 순서의 list[10]  (return_detail=True면 (list, dict))
    """
    mw_da = mw_kda * 1000.0
    ref = _nearest_mw_drug(mw_kda)          # kf_m/kr_m 차용원
    kf = DRUG_COMSOL[ref]["kf_m"] if kf_m is None else kf_m
    kr = DRUG_COMSOL[ref]["kr_m"] if kr_m is None else kr_m
    p = dict(_COMMON)
    p["sigma_ve"] = _sigma_curry(mw_da, _MW_R_VE_NM)
    p["sigma_le"] = _sigma_curry(mw_da, _MW_R_LE_NM)
    p["D_gel"]    = _d_gel_stokes(mw_da)
    p["p_ve"]     = _p_interp_from_anchors(mw_da, "p_ve")
    p["p_le"]     = _p_interp_from_anchors(mw_da, "p_le")
    p["k_decay"]  = k_decay_physics(mw_kda)
    p["kf_m"]     = kf
    p["kr_m"]     = kr
    vec = [p[k] for k in PARAM_ORDER]
    if return_detail:
        detail = dict(p)
        detail["_kfkr_source"] = ("user" if kf_m is not None else ref)
        detail["_mw_kda"] = mw_kda
        detail["_solute_radius_nm"] = _solute_radius_nm(mw_da)
        return vec, detail
    return vec


if __name__ == "__main__":
    print(f"{'Drug':5s} " + " ".join(f"{n:>9s}" for n in PARAM_ORDER) + "  source")
    for d in DRUG_COMSOL:
        v = get_param_vector(d)
        tag = "file" if d in FILE_MEASURED else "calc"
        print(f"{d:5s} " + " ".join(f"{x:>9.2e}" for x in v) + f"  [{tag}]")