""" verify.py — 검수 스크립트 (QA.md 근거 생성) 핵심 지표를 합성 데이터에 대해 실제 호출하고, 손계산과 비교한다. 종료코드 0 = 통과, 1 = 실패. """ import math import sys import numpy as np import pandas as pd import auc import metrics import deconvolution as dc import demo_data FAILS = [] WARNS = [] def check(name, cond, detail=""): status = "PASS" if cond else "FAIL" if not cond: FAILS.append(f"{name}: {detail}") print(f"[{status}] {name} {detail}") def approx(a, b, tol=1e-6): if a is None or b is None: return False if isinstance(a, float) and math.isnan(a): return False return abs(a - b) <= tol * max(1.0, abs(b)) print("=" * 60) print("GlucoDynamics-Kor 검수") print("=" * 60) # --------------------------------------------------------------------------- # 1) HOMA-IR 손계산 검증 # FPG=90 mg/dL -> 4.9957 mmol/L, FPI=10 uU -> HOMA = 4.9957*10/22.5 = 2.2203 # --------------------------------------------------------------------------- h = metrics.homa_ir(90.0, 10.0) expected = (90.0 / 18.0156 * 10.0) / 22.5 check("HOMA-IR 손계산", approx(h, expected, 1e-9), f"got={h:.6f} expected={expected:.6f}") # HOMA-IR 단순 케이스: FPG=126(=6.9939mmol), FPI=22.5/6.9939 -> HOMA=1 fpi_for_unit = 22.5 / (126.0 / 18.0156) check("HOMA-IR=1 단위검증", approx(metrics.homa_ir(126.0, fpi_for_unit), 1.0, 1e-9), f"got={metrics.homa_ir(126.0, fpi_for_unit):.6f}") # --------------------------------------------------------------------------- # 2) QUICKI 손계산: FPG=90, FPI=10 -> 1/(log10(10)+log10(90)) = 1/(1+1.9542) # --------------------------------------------------------------------------- q = metrics.quicki(90.0, 10.0) exp_q = 1.0 / (math.log10(10.0) + math.log10(90.0)) check("QUICKI 손계산", approx(q, exp_q, 1e-9), f"got={q:.6f} expected={exp_q:.6f}") # --------------------------------------------------------------------------- # 3) Matsuda 손계산: FPG=90,FPI=10,meanG=120,meanI=40 # 10000/sqrt(90*10*120*40)=10000/sqrt(4320000)=10000/2078.46=4.8113 # --------------------------------------------------------------------------- m = metrics.matsuda_isi(90.0, 10.0, 120.0, 40.0) exp_m = 10000.0 / math.sqrt(90.0 * 10.0 * 120.0 * 40.0) check("Matsuda 손계산", approx(m, exp_m, 1e-9), f"got={m:.6f} expected={exp_m:.6f}") # --------------------------------------------------------------------------- # 4) Insulinogenic index: (I30-I0)/(G30-G0) = (60-6)/(150-88) # --------------------------------------------------------------------------- ig = metrics.insulinogenic_index(88.0, 150.0, 6.0, 60.0) exp_ig = (60.0 - 6.0) / (150.0 - 88.0) check("Insulinogenic index 손계산", approx(ig, exp_ig, 1e-9), f"got={ig:.6f} expected={exp_ig:.6f}") # --------------------------------------------------------------------------- # 5) AUC trapezoidal: glucose 직사각형/삼각형 검증 # times=[0,60,120], values=[100,100,100] total = 100*120 = 12000 # --------------------------------------------------------------------------- a_flat = auc.trapz_auc([0, 60, 120], [100, 100, 100], mode="total") check("AUC total(평탄)", approx(a_flat, 12000.0, 1e-9), f"got={a_flat}") # incremental of flat (baseline=100) = 0 a_inc = auc.trapz_auc([0, 60, 120], [100, 100, 100], mode="incremental") check("AUC incremental(평탄)=0", approx(a_inc, 0.0, 1e-9), f"got={a_inc}") # triangle: [0,60,120], [0,100,0] total trapz = (0+100)/2*60 + (100+0)/2*60 = 6000 a_tri = auc.trapz_auc([0, 60, 120], [0, 100, 0], mode="total") check("AUC total(삼각)", approx(a_tri, 6000.0, 1e-9), f"got={a_tri}") # 불균등 간격: [0,30,120],[10,20,10] total = (10+20)/2*30 + (20+10)/2*90 = 450+1350=1800 a_uneven = auc.trapz_auc([0, 30, 120], [10, 20, 10], mode="total") check("AUC 불균등간격", approx(a_uneven, 1800.0, 1e-9), f"got={a_uneven}") # positive incremental: baseline=10, values [10,20,5] over [0,60,120] # inc = [0,10,-5]; zero crossing between t=60(10) and t=120(-5): # frac=10/15=0.6667 -> tc=60+0.6667*60=100; pos area: 0..60 (0->10) tri=300; # 60..100 (10->0) tri = 10/2*40=200; 100..120 (0->-5 clipped 0)=0 => 500 a_pos = auc.trapz_auc([0, 60, 120], [10, 20, 5], mode="positive_incremental") check("AUC positive-incremental(영점교차)", approx(a_pos, 500.0, 1e-6), f"got={a_pos}") # window 분리: 0-120 vs 0-180 on extended curve times_ext = [0, 30, 60, 90, 120, 150, 180] vals_ext = [100, 100, 100, 100, 100, 100, 100] a120 = auc.trapz_auc(times_ext, vals_ext, 0, 120, "total") a180 = auc.trapz_auc(times_ext, vals_ext, 0, 180, "total") check("AUC window 0-120", approx(a120, 12000.0, 1e-9), f"got={a120}") check("AUC window 0-180", approx(a180, 18000.0, 1e-9), f"got={a180}") # mean_over_curve: flat 100 -> 100 mean_flat = auc.mean_over_curve([0, 60, 120], [100, 100, 100]) check("mean_over_curve(평탄)=100", approx(mean_flat, 100.0, 1e-9), f"got={mean_flat}") # --------------------------------------------------------------------------- # 6) ISR deconvolution: 합성 C-peptide 곡선에 대해 NaN/에러 없이 산출 # --------------------------------------------------------------------------- df = demo_data.make_ogtt_normal() times = df["time_min"].to_numpy(dtype=float) cpep_ngml = df["cpeptide"].to_numpy(dtype=float) cpep_pmol = metrics.cpeptide_ngml_to_pmol_L(cpep_ngml) bsa = dc.bsa_dubois(75.0, 170.0) params = dc.standard_kinetics(age_years=45, bsa_m2=bsa, diabetic=False, obese=False) res = dc.deconvolve_isr(times, cpep_pmol, params, dt=1.0) isr = res["isr_pmol_min"] check("ISR 산출 비어있지않음", isr.size > 0, f"n={isr.size} method={res['method']}") check("ISR 유한값", np.all(np.isfinite(isr)) and isr.size > 0, f"any_nan={np.any(~np.isfinite(isr)) if isr.size else 'NA'}") check("ISR 비음수", isr.size > 0 and np.all(isr >= -1e-6), "음수 ISR 발견" ) check("ISR total AUC 유한", math.isfinite(res["isr_total_auc"]), f"auc={res['isr_total_auc']:.2f} basal={res['isr_basal']:.3f}") print(f" [info] scipy 사용여부: {dc._HAVE_SCIPY}, ISR 방법={res['method']}") # kinetics sanity: a > b (short rate > long rate), V > 0 check("kinetics a>b>0", params["a"] > params["b"] > 0, f"a={params['a']:.4f} b={params['b']:.5f} V={params['V_L']:.3f}") # --------------------------------------------------------------------------- # 7) 전체 패널: 3개 OGTT 데모에 대해 NaN-안전 산출 # --------------------------------------------------------------------------- for label, maker in [("NGT", demo_data.make_ogtt_normal), ("IGT", demo_data.make_ogtt_igt), ("T2DM", demo_data.make_ogtt_t2dm)]: d = maker() panel = metrics.compute_panel( d["time_min"].to_numpy(float), d["glucose"].to_numpy(float), d["insulin"].to_numpy(float), d["cpeptide"].to_numpy(float), weight_kg=75, bmi=26, auc_module=auc) core = ["HOMA-IR", "Matsuda/ISI", "Insulinogenic index", "Oral Disposition Index"] finite_core = all(math.isfinite(panel[k]) for k in core) check(f"패널 핵심지표 유한 [{label}]", finite_core, " ".join(f"{k}={panel[k]:.3f}" for k in core)) # 생리학적 방향성: T2DM HOMA-IR > NGT HOMA-IR, NGT Matsuda > T2DM Matsuda pn = metrics.compute_panel(*[demo_data.make_ogtt_normal()[c].to_numpy(float) for c in ["time_min", "glucose", "insulin", "cpeptide"]], weight_kg=75, bmi=24, auc_module=auc) pt = metrics.compute_panel(*[demo_data.make_ogtt_t2dm()[c].to_numpy(float) for c in ["time_min", "glucose", "insulin", "cpeptide"]], weight_kg=85, bmi=30, auc_module=auc) check("방향성 Matsuda NGT>T2DM", pn["Matsuda/ISI"] > pt["Matsuda/ISI"], f"NGT={pn['Matsuda/ISI']:.2f} T2DM={pt['Matsuda/ISI']:.2f}") # --------------------------------------------------------------------------- # 8) Clamp M-value # --------------------------------------------------------------------------- dcl = demo_data.make_clamp() mval = metrics.clamp_m_value(dcl["gir"].to_numpy(float), dcl["time_min"].to_numpy(float), (90, 120)) moi = metrics.clamp_m_over_i(mval, float(dcl["insulin"].iloc[-1])) check("Clamp M-value 유한&양수", math.isfinite(mval) and mval > 0, f"M={mval:.3f} mg/kg/min, M/I={moi:.4f}") # --------------------------------------------------------------------------- # 9) 단위 변환 왕복 # --------------------------------------------------------------------------- g = 100.0 back = metrics.glucose_to_mgdl(metrics.glucose_mgdl_to_mmol(g), "mmol/L") check("glucose 단위 왕복", approx(back, g, 1e-9), f"got={back}") ins = 12.0 back_i = metrics.insulin_to_uU(metrics.insulin_uU_to_pmol(ins), "pmol/L") check("insulin 단위 왕복", approx(back_i, ins, 1e-9), f"got={back_i}") # --------------------------------------------------------------------------- # 10) CSV 로드 테스트 # --------------------------------------------------------------------------- import os for f in ["demo_ogtt_normal.csv", "demo_ogtt_igt.csv", "demo_ogtt_t2dm.csv", "demo_clamp.csv"]: p = os.path.join("data", f) try: dd = pd.read_csv(p) check(f"CSV 로드 [{f}]", "time_min" in dd.columns and len(dd) >= 5, f"rows={len(dd)} cols={list(dd.columns)}") except Exception as e: check(f"CSV 로드 [{f}]", False, f"err={e}") print("=" * 60) if FAILS: print(f"실패 {len(FAILS)}건:") for f in FAILS: print(" -", f) sys.exit(1) else: print("모든 검수 통과") sys.exit(0)