{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# ============================================================\n# PhysioNet ECG Digitization — Full Pipeline (time-calibrated + phase-align)\n# + Paper speed(25/50) detection & override\n# + TTA with shared time-scale (spp_ref) and clamp\n# + Light 1D-CNN denoiser (optional training & inference)\n# ============================================================\n\nimport os, glob, cv2, math, warnings\nwarnings.filterwarnings(\"ignore\")\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.auto import tqdm\nfrom scipy.signal import butter, filtfilt, find_peaks\n\n# -----------------------\n# Paths / Config\n# -----------------------\nTRAIN_DIR    = '/kaggle/input/physionet-ecg-image-digitization/train/'\nTRAIN_CSV    = '/kaggle/input/physionet-ecg-image-digitization/train.csv'\nTEST_DIR     = '/kaggle/input/physionet-ecg-image-digitization/test/'\nTEST_CSV     = '/kaggle/input/physionet-ecg-image-digitization/test.csv'\nWORK_DIR     = '/kaggle/working'\n\nTEMPLATE_NPZ = os.path.join(WORK_DIR, 'lead_templates_beats.npz')  # 저장만, 로딩은 안 함\nVIS_DIR      = os.path.join(WORK_DIR, 'train_vis')\nos.makedirs(VIS_DIR, exist_ok=True)\n\nSUBMISSION_CSV = os.path.join(WORK_DIR, 'submission.csv')\n\n# 리드 순서 (3x4 패널 가정)\nLEAD_GRID = [\n    [\"I\",\"II\",\"III\",\"aVR\"],\n    [\"aVL\",\"aVF\",\"V1\",\"V2\"],\n    [\"V3\",\"V4\",\"V5\",\"V6\"],\n]\nLEADS = sum(LEAD_GRID, [])\n\n# 제출 스케일\nMIN_VAL, MAX_VAL = 0.0, 0.07\n\n# ROI/DP/블렌딩 파라미터 (약간 튜닝)\nINK_GRAY_THR     = 48\nLOCAL_INK_THR    = 0.06\nMARGIN_COLS      = 8\n\nDP_LAMBDA        = 1.25\nDP_WIN_FRAC      = 0.10\nDP_EDGE_GAIN     = 0.45\n\nCONF_BAND        = 3\nUSE_IMG_MIN_CONF = 0.20\n\n# 템플릿 (beat-wise)\nTEMPLATE_BEAT_LEN = 360\nR_PRE_S           = 0.22\nR_POST_S          = 0.42\n\n# ---- Paper speed override (환경변수로 강제 가능: '25' 또는 '50') ----\nPAPER_SPEED_OVERRIDE = os.getenv(\"ECG_PAPER_SPEED\", \"\").strip()  # \"25\" or \"50\" or \"\"\n\n# -----------------------\n# TTA (증강) 설정\n# -----------------------\nTTA_ENABLED = bool(int(os.getenv(\"ECG_TTA\", \"1\")))   # 1=켜기, 0=끄기\nMAX_AUG     = int(os.getenv(\"ECG_TTA_N\", \"6\"))\nTTA_AGG     = os.getenv(\"ECG_TTA_AGG\", \"weighted_mean\")  # 'weighted_mean' or 'median'\n\n# 미세 지오메트릭 + 포토메트릭 (시간스케일 안정 위해 각/시어는 매우 작게)\nAUG_PRESETS = [\n    {\"angle\": -0.2}, {\"angle\": +0.2},\n    {\"shear\": -0.3}, {\"shear\": +0.3},\n    {\"tx\": -6.0}, {\"tx\": +6.0},\n    {\"alpha\": 1.10}, {\"alpha\": 0.90},\n    {\"gamma\": 0.90}, {\"gamma\": 1.10},\n]\n\n# -----------------------\n# Denoiser(학습/추론) 설정\n# -----------------------\nDENOISER_ENABLE      = bool(int(os.getenv(\"ECG_DENOISER\", \"1\")))   # 1=사용, 0=비사용\nDENOISER_TRAIN       = bool(int(os.getenv(\"ECG_DENOISER_TRAIN\", \"0\"))) # 1=학습 실행\nDENOISER_EPOCHS      = int(os.getenv(\"ECG_DENOISER_EPOCHS\", \"1\"))\nDENOISER_LR          = float(os.getenv(\"ECG_DENOISER_LR\", \"1e-3\"))\nDENOISER_PATH        = os.getenv(\"ECG_DENOISER_PATH\", os.path.join(WORK_DIR, \"denoiser1d.pt\"))\nDENOISER_FREQ_LOSS_W = float(os.getenv(\"ECG_DENOISER_FREQW\", \"0.2\"))  # 주파수 손실 가중치\nDENOISER_USE_TPL_CH  = bool(int(os.getenv(\"ECG_DENOISER_TPLCH\", \"1\"))) # 템플릿 채널 사용\n\n# -----------------------\n# Utils\n# -----------------------\ndef lowpass(x, fs, cutoff_hz=15.0, order=2):\n    x = np.asarray(x, dtype=np.float32)\n    if x.size <= 10: return x\n    nyq = 0.5*float(fs); wn = min(cutoff_hz/max(nyq,1e-6), 0.99)\n    b,a = butter(order, wn, btype='low')\n    return filtfilt(b,a,x).astype(np.float32)\n\ndef zscore(x):\n    x = np.asarray(x, dtype=np.float32)\n    return (x - np.mean(x)) / (np.std(x)+1e-8)\n\ndef rescale_range(x, lo=MIN_VAL, hi=MAX_VAL):\n    x = np.asarray(x, dtype=np.float32)\n    mn, mx = float(np.min(x)), float(np.max(x))\n    if not np.isfinite(mn) or not np.isfinite(mx) or mx <= mn:\n        return np.full_like(x, (lo+hi)/2, dtype=np.float32)\n    y = (x - mn) / (mx - mn)\n    return (lo + y * (hi-lo)).astype(np.float32)\n\ndef tukey_window(n, alpha=0.25):\n    if n <= 1: return np.ones(n, np.float32)\n    w = np.ones(n, np.float32)\n    e = int(alpha*(n-1)/2.0)\n    if e > 0:\n        ramp = (1 - np.cos(np.linspace(0, np.pi, e*2, dtype=np.float32)))/2.0\n        w[:e] = ramp[:e]; w[-e:] = ramp[-e:]\n    return w\n\ndef sigmoid_blend(alpha, k=8.0, bias=-0.10, lo=0.12, hi=0.92):\n    s = 1.0/(1.0 + np.exp(-k*(alpha + bias)))\n    return float(np.clip(lo + (hi-lo)*s, lo, hi))\n\n# -----------------------\n# Augmentation helpers\n# -----------------------\ndef _photometric_bgr(img_bgr, alpha=1.0, beta=0.0, gamma=1.0):\n    x = img_bgr.astype(np.float32)\n    if alpha is None: alpha = 1.0\n    if beta  is None: beta  = 0.0\n    if gamma is None: gamma = 1.0\n    x = x * float(alpha) + float(beta)\n    x = np.clip(x, 0, 255)\n    if abs(float(gamma) - 1.0) > 1e-6:\n        x = (x / 255.0) ** (1.0 / float(gamma))\n        x = np.clip(x * 255.0, 0, 255)\n    return x.astype(np.uint8)\n\ndef _affine_shear_rotate_translate(img_bgr, angle=0.0, shear=0.0, tx=0.0, ty=0.0, scale=1.0):\n    H, W = img_bgr.shape[:2]\n    cx, cy = (W-1)*0.5, (H-1)*0.5\n    def _to33(M23):\n        M33 = np.eye(3, dtype=np.float32)\n        M33[:2,:3] = M23\n        return M33\n    C    = np.array([[1,0,-cx],[0,1,-cy],[0,0,1]], np.float32)\n    Cinv = np.array([[1,0, cx],[0,1, cy],[0,0,1]], np.float32)\n    R23  = cv2.getRotationMatrix2D((0,0), float(angle), float(scale))\n    R    = _to33(R23)\n    sh = math.tan(math.radians(float(shear)))\n    S  = np.array([[1, sh, 0],[0,1,0],[0,0,1]], np.float32)\n    T  = np.array([[1,0,float(tx)],[0,1,float(ty)],[0,0,1]], np.float32)\n    M  = T @ Cinv @ R @ S @ C\n    M23 = M[:2,:]\n    return cv2.warpAffine(img_bgr, M23, (W, H), flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REPLICATE)\n\ndef augment_panel(panel_bgr, angle=0.0, shear=0.0, tx=0.0, ty=0.0, scale=1.0,\n                  alpha=1.0, beta=0.0, gamma=1.0):\n    out = _affine_shear_rotate_translate(panel_bgr, angle=angle, shear=shear, tx=tx, ty=ty, scale=scale)\n    out = _photometric_bgr(out, alpha=alpha, beta=beta, gamma=gamma)\n    return out\n\n# -----------------------\n# 0) 패널 분할 (3x4 균등)\n# -----------------------\ndef split_3x4_panels(img_bgr, trim=6):\n    H, W = img_bgr.shape[:2]\n    w = W // 4; h = H // 3\n    panels = {}\n    for r in range(3):\n        for c in range(4):\n            y0, y1 = r*h, (r+1)*h\n            x0, x1 = c*w, (c+1)*w\n            roi = img_bgr[y0:y1, x0:x1]\n            if trim>0: roi = roi[trim:-trim, trim:-trim]\n            lead = LEAD_GRID[r][c]\n            panels[lead] = roi\n    return panels\n\n# -----------------------\n# 1) 격자 억제 + 활성구간\n# -----------------------\ndef degrid_gray(bgr):\n    hsv = cv2.cvtColor(bgr, cv2.COLOR_BGR2HSV)\n    h,s,v = cv2.split(hsv)\n    red = ((h < 12) | (h > 168)) & (s > 60)\n    gray = v.copy()\n    gray[red] = np.median(v)\n    gray = cv2.GaussianBlur(gray, (3,3), 0)\n    gray = cv2.normalize(gray, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)\n    return gray\n\ndef find_active_columns_longest(gray, ink_thr=INK_GRAY_THR, min_ratio=LOCAL_INK_THR, margin=MARGIN_COLS):\n    H, W = gray.shape\n    inv = 255 - gray\n    ink = (inv > ink_thr).astype(np.uint8)\n    col_sum = ink.sum(axis=0)\n    thr = max(int(min_ratio*H), 5)\n    mask = col_sum >= thr\n\n    best_len, best = 0, (0, W-1); s = None\n    for i, v in enumerate(mask):\n        if v and s is None: s = i\n        if ((not v) or i==W-1) and s is not None:\n            e = i if not v else i\n            if e - s > best_len:\n                best_len = e - s\n                best = (max(0, s-margin), min(W-1, e+margin))\n            s = None\n    lo, hi = best\n    return lo, hi\n\n# -----------------------\n# 1.5) 격자 간격 기반 보정 (NEW)\n# -----------------------\ndef _red_mask(bgr):\n    hsv = cv2.cvtColor(bgr, cv2.COLOR_BGR2HSV)\n    h,s,v = cv2.split(hsv)\n    m1 = (h <= 10) & (s >= 50) & (v >= 40)\n    m2 = (h >= 170) & (s >= 50) & (v >= 40)\n    mask = (m1 | m2).astype(np.uint8)*255\n    mask = cv2.medianBlur(mask, 3)\n    mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, np.ones((3,3), np.uint8))\n    return mask\n\ndef _autocorr_period(sig, min_lag=6, max_lag=None):\n    sig = np.asarray(sig, np.float32)\n    sig = (sig - sig.mean())/(sig.std()+1e-6)\n    ac = np.correlate(sig, sig, mode='full')[len(sig)-1:]\n    if max_lag is None: max_lag = len(ac)//2\n    min_lag = max(min_lag, 1); max_lag = max(max_lag, min_lag+1)\n    if max_lag <= min_lag: return None\n    k = np.argmax(ac[min_lag:max_lag]) + min_lag\n    return int(k)\n\ndef measure_grid_spacing_bgr(panel_bgr):\n    H, W = panel_bgr.shape[:2]\n    mask = _red_mask(panel_bgr)\n    col = mask.sum(axis=0)\n    row = mask.sum(axis=1)\n\n    dx_small = _autocorr_period(col, min_lag=max(6, W//200), max_lag=W//2)\n    dy_small = _autocorr_period(row, min_lag=max(6, H//200), max_lag=H//2)\n\n    dx_big = int(dx_small*5) if dx_small else None\n    dy_big = int(dy_small*5) if dy_small else None\n    if dx_big is not None and (dx_big <= 0 or dx_big > W): dx_big = None\n    if dy_big is not None and (dy_big <= 0 or dy_big > H): dy_big = None\n    return dx_big, dy_big\n\n# ---- 속도표기(25/50) 감지 ----\ndef detect_paper_speed_label(img_bgr):\n    \"\"\"\n    returns: 25, 50, or None\n    - 하단 30%에서 '25' / '50' 템플릿매칭\n    - 빨간 글자 우선, 약하면 흑백 adaptive threshold\n    - 환경변수 ECG_PAPER_SPEED로 강제 가능\n    \"\"\"\n    if PAPER_SPEED_OVERRIDE in (\"25\",\"50\"):\n        return int(PAPER_SPEED_OVERRIDE)\n\n    H, W = img_bgr.shape[:2]\n    roi = img_bgr[int(H*0.70):, :]\n    hsv = cv2.cvtColor(roi, cv2.COLOR_BGR2HSV)\n    h,s,v = cv2.split(hsv)\n    red = (((h < 15) | (h > 165)) & (s > 60) & (v > 40)).astype(np.uint8)*255\n\n    gray = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)\n    mask = red.copy()\n    if mask.sum() < 0.002 * mask.size:\n        mask = cv2.adaptiveThreshold(gray, 255, cv2.ADAPTIVE_THRESH_MEAN_C,\n                                     cv2.THRESH_BINARY_INV, 31, 5)\n    vis = cv2.bitwise_and(gray, mask)\n\n    def best_score_for(label):\n        base = np.zeros((40, 90), np.uint8)\n        cv2.putText(base, label, (2, 34), cv2.FONT_HERSHEY_SIMPLEX, 1.2, 255, 2, cv2.LINE_AA)\n        best = -1.0\n        for sc in (0.8, 0.9, 1.0, 1.1, 1.2):\n            t = cv2.resize(base, (int(base.shape[1]*sc), int(base.shape[0]*sc)), interpolation=cv2.INTER_AREA)\n            if vis.shape[0] < t.shape[0] or vis.shape[1] < t.shape[1]:\n                continue\n            res = cv2.matchTemplate(vis, t, cv2.TM_CCOEFF_NORMED)\n            if res.size:\n                best = max(best, float(res.max()))\n        return best\n\n    s25 = best_score_for(\"25\")\n    s50 = best_score_for(\"50\")\n    thr = 0.35\n    if s25 < thr and s50 < thr:\n        return None\n    return 25 if s25 >= s50 else 50\n\n# -----------------------\n# (A) Hough fallback (기존)\n# -----------------------\ndef sec_per_pixel_from_grid_calibrated(gray, x_lo, x_hi, T_known):\n    g = cv2.equalizeHist(gray)\n    edges = cv2.Canny(g, 40, 100)\n    lines = cv2.HoughLinesP(edges, 1, np.pi/180, threshold=120, minLineLength=40, maxLineGap=6)\n\n    d = None\n    if lines is not None:\n        xs = [int(round((x1+x2)/2)) for x1,y1,x2,y2 in lines[:,0]\n              if abs(x1-x2) < 3 and abs(y1-y2) > 25]\n        xs = np.array(sorted(set(xs)))\n        if len(xs) >= 3:\n            diffs = np.diff(xs)\n            diffs = diffs[(diffs>5) & (diffs<2000)]\n            if len(diffs): d = float(np.median(diffs))\n\n    cands = []\n    if d and d > 0:\n        cands += [0.2/d, 0.1/d]\n\n    Wroi = max(1, (x_hi - x_lo + 1))\n    if not cands:\n        return T_known / Wroi\n\n    best = None; best_err = 1e9\n    for spp in cands:\n        err = abs(Wroi*spp - T_known)\n        if err < best_err:\n            best_err, best = err, spp\n    return best\n\n# ---- s/px 산출 (속도표기 반영) ----\ndef sec_per_pixel_from_grid(panel_bgr, x_lo, x_hi, T_known, gray_fallback=None, paper_speed=None):\n    dx_big, _ = measure_grid_spacing_bgr(panel_bgr)\n    Wroi = max(1, (x_hi - x_lo + 1))\n\n    cands = []\n    if dx_big:\n        if paper_speed == 25:\n            cands = [0.2/float(dx_big)]\n        elif paper_speed == 50:\n            cands = [0.1/float(dx_big)]\n        else:\n            cands = [0.2/float(dx_big), 0.1/float(dx_big)]\n\n    if not cands and gray_fallback is not None:\n        spp_hough = sec_per_pixel_from_grid_calibrated(gray_fallback, x_lo, x_hi, T_known)\n        if spp_hough is not None:\n            cands.append(spp_hough)\n\n    if not cands:\n        return T_known / float(Wroi)\n\n    best_spp, best_err = None, 1e18\n    for spp in cands:\n        err = abs(Wroi*spp - T_known)\n        if err < best_err:\n            best_err, best_spp = err, spp\n\n    if len(cands) == 2:\n        e0 = abs(Wroi * cands[0] - T_known)\n        e1 = abs(Wroi * cands[1] - T_known)\n        if abs(e0 - e1) < 0.06 * T_known:\n            best_spp = cands[0]  # 25 mm/s 선호\n\n    return best_spp\n\n# -----------------------\n# 2) DP centerline\n# -----------------------\ndef dp_trace_center(gray, x_lo, x_hi):\n    H, W = gray.shape\n    x_lo = max(0, int(x_lo)); x_hi = min(W-1, int(x_hi))\n    if x_hi <= x_lo: x_lo, x_hi = 0, W-1\n\n    g = gray.astype(np.float32)/255.0\n    inv = 1.0 - g\n    gx = cv2.Sobel(g, cv2.CV_32F, 1, 0, ksize=3)\n    gy = cv2.Sobel(g, cv2.CV_32F, 0, 1, ksize=3)\n    edge = cv2.magnitude(gx, gy); edge = edge / (edge.max()+1e-6)\n\n    cost_img = (1.0 - (inv[:, x_lo:x_hi+1] + DP_EDGE_GAIN*edge[:, x_lo:x_hi+1]))\n    cost_img = np.clip(cost_img, 0.0, 2.0).astype(np.float32)\n\n    Hc, Wc = cost_img.shape\n    win = max(1, int(DP_WIN_FRAC*Hc)); lam = DP_LAMBDA\n    C = np.full((Hc, Wc), np.inf, np.float32)\n    P = np.full((Hc, Wc), -1, np.int32)\n    C[:,0] = cost_img[:,0]\n\n    for x in range(1, Wc):\n        prev = C[:, x-1]\n        for dy in range(-win, win+1):\n            trans = 0.0 if dy==0 else lam*abs(dy)\n            if dy < 0:\n                src = prev[-dy:Hc]; dst_rows = slice(0, Hc+dy)\n            elif dy > 0:\n                src = prev[0:Hc-dy]; dst_rows = slice(dy, Hc)\n            else:\n                src = prev; dst_rows = slice(0, Hc)\n            cand = src + cost_img[dst_rows, x] + trans\n            better = cand < C[dst_rows, x]\n            C[dst_rows, x][better] = cand[better]\n            P[dst_rows, x][better] = np.arange(len(src), dtype=np.int32)[better]\n\n    y = int(np.argmin(C[:, -1]))\n    ys = [y]\n    for x in range(Wc-1, 0, -1):\n        py = P[ys[-1], x]\n        if py < 0: py = ys[-1]\n        ys.append(py)\n    ys = ys[::-1]\n    xs = np.arange(x_lo, x_hi+1, dtype=np.int32)\n    if len(xs) != len(ys): xs = np.linspace(x_lo, x_hi, len(ys), dtype=np.int32)\n    return xs, np.array(ys, np.int32)\n\ndef path_confidence(gray, xs, ys, band=CONF_BAND):\n    H, W = gray.shape\n    xs = np.clip(np.asarray(xs, np.int32), 0, W-1)\n    ys = np.clip(np.asarray(ys, np.int32), 0, H-1)\n\n    inv = 255 - gray\n    gx = cv2.Sobel(gray, cv2.CV_32F, 1,0,ksize=3)\n    gy = cv2.Sobel(gray, cv2.CV_32F, 0,1,ksize=3)\n    mag = cv2.magnitude(gx, gy)\n\n    ink_cnt, ink_hit, edge_acc = 0, 0, 0.0\n    for x,y in zip(xs,ys):\n        y0, y1 = max(0,y-band), min(H-1,y+band)\n        line = inv[y0:y1+1, x]; edge = mag[y0:y1+1, x]\n        ink_cnt += line.size; ink_hit += int((line > INK_GRAY_THR).sum())\n        edge_acc += float(np.mean(edge))\n    ink_ratio = 0.0 if ink_cnt==0 else ink_hit/float(ink_cnt)\n    edge_norm = (edge_acc / max(len(xs),1)) / 50.0\n    edge_norm = float(np.clip(edge_norm, 0.0, 1.0))\n    conf = 0.35 * np.clip((ink_ratio - 0.10)/0.35, 0.0, 1.0) + 0.65 * edge_norm\n    return float(np.clip(conf, 0.0, 1.0))\n\n# ROI → 값 시퀀스 (시간보정 포함)\ndef roi_series_from_path_t(gray, xs, ys, spp, n_out, fs):\n    inv = 255 - gray\n    band = 2\n    vals = []\n    for x,y in zip(xs,ys):\n        y0, y1 = max(0,y-band), min(gray.shape[0]-1,y+band)\n        vals.append(np.mean(inv[y0:y1+1, x]))\n    vals = np.array(vals, np.float32)\n    vals = zscore(vals); vals = lowpass(vals, fs, 20.0, 2)\n\n    T = (n_out-1)/float(fs)\n    t_pix = (np.arange(len(vals), dtype=np.float32)) * float(spp)\n    t_pix -= t_pix[0]\n    t_pix = np.clip(t_pix, 0.0, T)\n\n    t_out = np.linspace(0.0, T, n_out, dtype=np.float32)\n    out = np.interp(t_out, t_pix, vals).astype(np.float32)\n    return out\n\n# -----------------------\n# 서로상관 기반 시프트 정렬\n# -----------------------\ndef phase_align_to_template(y_img, y_tpl):\n    a = zscore(y_img); b = zscore(y_tpl)\n    L = min(len(a), len(b))\n    if L < 8: return y_img\n    c = np.correlate(a[:L], b[:L], mode=\"full\")\n    lag = int(np.argmax(c) - (L-1))\n    if lag > 0:\n        y_img = np.r_[np.zeros(lag, dtype=y_img.dtype), y_img[:-lag]]\n    elif lag < 0:\n        lag = -lag\n        y_img = np.r_[y_img[lag:], np.zeros(lag, dtype=y_img.dtype)]\n    return y_img\n\n# -----------------------\n# 3) Beat-wise Templates (train)\n# -----------------------\ndef bandpass_ecg(x, fs, lo=5.0, hi=25.0, order=2):\n    nyq = 0.5*fs\n    lo = max(lo/nyq, 1e-3); hi = min(hi/nyq, 0.99)\n    b,a = butter(order, [lo,hi], btype='band')\n    return filtfilt(b,a,x).astype(np.float32)\n\ndef build_lead_templates_beatwise(train_csv_path, train_dir,\n                                  leads=LEADS,\n                                  beat_len=TEMPLATE_BEAT_LEN,\n                                  pre_R=R_PRE_S, post_R=R_POST_S):\n    meta = pd.read_csv(train_csv_path)\n    beats = {ld: [] for ld in leads}\n    used  = {ld: 0 for ld in leads}\n    for row in tqdm(meta.itertuples(index=False), total=len(meta), desc=\"Build beatwise templates\"):\n        rid = str(row.id); fs = int(row.fs)\n        csvp = os.path.join(train_dir, rid, f\"{rid}.csv\")\n        if not os.path.exists(csvp): continue\n        try: df = pd.read_csv(csvp)\n        except: continue\n        for ld in leads:\n            if ld not in df.columns: continue\n            y = df[ld].dropna().to_numpy(np.float32)\n            if y.size < 50: continue\n            y0 = zscore(y)\n            yf = bandpass_ecg(y0, fs, 5, 25, 2)\n            peaks, _ = find_peaks(yf, distance=int(0.35*fs), prominence=max(0.5*np.std(yf), 0.2))\n            if len(peaks) < 2: continue\n            n_pre  = int(round(pre_R*fs)); n_post = int(round(post_R*fs))\n            for pk in peaks:\n                a = pk - n_pre; b = pk + n_post\n                if a < 0 or b >= len(y0): continue\n                seg = y0[a:b+1]\n                seg_rs = np.interp(np.linspace(0,1,beat_len, dtype=np.float32),\n                                   np.linspace(0,1,len(seg), dtype=np.float32),\n                                   seg).astype(np.float32)\n                beats[ld].append(seg_rs); used[ld]+=1\n    templates = {}\n    for ld in leads:\n        if beats[ld]:\n            arr = np.vstack(beats[ld])\n            tpl = np.median(arr, axis=0).astype(np.float32)\n        else:\n            t = np.linspace(0,1,beat_len, dtype=np.float32)\n            tpl = np.sin(2*np.pi*t).astype(np.float32)\n        templates[ld] = zscore(tpl)\n    return templates, used\n\ndef stretch_template_for_bpm(template, fs, bpm):\n    beat_samp = max(6, int(round(fs * 60.0 / max(bpm, 1e-3))))\n    x = np.linspace(0, 1, beat_samp, dtype=np.float32)\n    tx = np.linspace(0, 1, len(template), dtype=np.float32)\n    return np.interp(x, tx, template).astype(np.float32)\n\ndef template_series(template, fs, n_out, bpm, amp=1.0):\n    one = stretch_template_for_bpm(template, fs, bpm)\n    reps = int(np.ceil(n_out / len(one)))\n    y = np.tile(one, reps)[:n_out]\n    y = zscore(y) * float(amp)\n    y = lowpass(y, fs, 20.0, 2)\n    return y\n\ndef estimate_bpm_from_series(y, fs):\n    y = zscore(y); y = lowpass(y, fs, 15.0, 2)\n    ac = np.correlate(y, y, mode='full')[len(y)-1:]\n    ac[:int(0.30*fs)] = 0\n    pk = int(np.argmax(ac[:max(int(2.0*fs),1)]))\n    if pk <= 0: return 75.0\n    rr = pk/float(fs)\n    return float(np.clip(60.0/max(rr, 1e-3), 40.0, 160.0))\n\n# -----------------------\n# 4) 이미지 패널 → 시리즈 (시간보정 + 시프트정렬 + 블렌딩)\n#     - 원본에서 s/px 한 번만 추정(spp_ref) → 모든 TTA에 공유\n#     - s/px를 T_known/Wroi의 ±5%로 클램프\n#     - paper_speed(25/50) 지정 시 해당 후보만 사용\n# -----------------------\ndef panel_to_series(panel_bgr, fs, n_out, template, debug_path=None,\n                    tta=TTA_ENABLED, aug_presets=AUG_PRESETS, aggregate=TTA_AGG,\n                    paper_speed=None):\n\n    def _decode_one(img_bgr, spp_override=None, dbg=False):\n        gray = degrid_gray(img_bgr)\n        x_lo, x_hi = find_active_columns_longest(gray)\n        xs, ys = dp_trace_center(gray, x_lo, x_hi)\n        conf = path_confidence(gray, xs, ys, band=CONF_BAND)\n\n        T = (n_out-1)/float(fs)\n        if spp_override is None:\n            spp = sec_per_pixel_from_grid(img_bgr, x_lo, x_hi, T_known=T,\n                                          gray_fallback=gray, paper_speed=paper_speed)\n        else:\n            spp = float(spp_override)\n\n        # 가드레일: 전체 시간 길이 유지\n        Wroi  = max(1, (x_hi - x_lo + 1))\n        ideal = T / float(Wroi)\n        spp   = float(np.clip(float(spp), 0.95*ideal, 1.05*ideal))\n\n        # ROI → 시리즈 (시간보정)\n        y_img = roi_series_from_path_t(gray, xs, ys, spp, n_out, fs)\n\n        # BPM 추정 + 템플릿\n        bpm_est = estimate_bpm_from_series(y_img, fs) if conf > 0.15 else 75.0\n        y_tpl   = template_series(template, fs, n_out, bpm=bpm_est, amp=1.0)\n\n        # 서로상관 기반 시프트 정렬\n        y_img = phase_align_to_template(y_img, y_tpl)\n\n        # 신뢰도 기반 블렌딩\n        a = sigmoid_blend(conf)\n        if conf < USE_IMG_MIN_CONF: a *= 0.5\n        y_mix = rescale_range(a*y_img + (1.0-a)*y_tpl, MIN_VAL, MAX_VAL)\n\n        if dbg and (debug_path is not None):\n            H,W = gray.shape\n            color_dbg = panel_bgr.copy()\n            cv2.rectangle(color_dbg, (x_lo,0), (x_hi,H-1), (0,255,255), 2)\n            for x,y in zip(xs,ys):\n                cv2.circle(color_dbg, (int(x),int(y)), 1, (0,0,255), -1)\n            cv2.imwrite(debug_path, color_dbg)\n\n        return y_mix, float(conf), float(bpm_est), float(spp)\n\n    # 1) 원본으로 먼저 복원 → 참조 s/px 확보\n    y0, c0, b0, spp_ref = _decode_one(panel_bgr, spp_override=None, dbg=True)\n\n    if not tta:\n        return y0, c0, b0\n\n    # 2) 증강들은 모두 spp_ref를 강제 사용(시간 스케일 고정)\n    ys, confs, bpms = [y0], [c0], [b0]\n    use_presets = (aug_presets or [])[:max(0, int(MAX_AUG))]\n\n    for cfg in use_presets:\n        try:\n            aug = augment_panel(\n                panel_bgr,\n                angle=cfg.get(\"angle\", 0.0),\n                shear=cfg.get(\"shear\", 0.0),\n                tx=cfg.get(\"tx\", 0.0),\n                ty=cfg.get(\"ty\", 0.0),\n                scale=cfg.get(\"scale\", 1.0),\n                alpha=cfg.get(\"alpha\", 1.0),\n                beta=cfg.get(\"beta\", 0.0),\n                gamma=cfg.get(\"gamma\", 1.0),\n            )\n            y, c, b, _ = _decode_one(aug, spp_override=spp_ref, dbg=False)\n            ys.append(y); confs.append(c); bpms.append(b)\n        except Exception:\n            continue\n\n    Y = np.vstack(ys)\n    C = np.asarray(confs, np.float32)\n    W = np.maximum(C**2, 1e-3)\n\n    if aggregate == \"median\":\n        y_agg = np.median(Y, axis=0)\n    else:\n        y_agg = np.average(Y, axis=0, weights=W)\n\n    y_agg   = rescale_range(y_agg, MIN_VAL, MAX_VAL)\n    conf_ag = float(np.average(C, weights=W))\n    bpm_ag  = float(np.average(np.asarray(bpms, np.float32), weights=W))\n    return y_agg, conf_ag, bpm_ag\n\n# ============================================================\n# 5) Light 1D-CNN Denoiser (옵션): 학습/로드/적용\n# ============================================================\ntry:\n    import torch, torch.nn as nn\n    TORCH_OK = True\nexcept Exception:\n    TORCH_OK = False\n\nclass Denoiser1D(nn.Module):\n    def __init__(self, in_ch=2):\n        super().__init__()\n        self.net = nn.Sequential(\n            nn.Conv1d(in_ch, 32, 9, padding=4), nn.ReLU(),\n            nn.Conv1d(32, 32, 9, padding=4), nn.ReLU(),\n            nn.Conv1d(32, 16, 5, padding=2), nn.ReLU(),\n            nn.Conv1d(16, 1,  1)\n        )\n    def forward(self, x):\n        # residual: 입력의 첫 채널(y_mix_norm)에 잔차를 더해 보정값 산출\n        res = self.net(x)\n        return x[:, :1, :] + res\n\ndef spectral_mse(a, b):\n    # a,b: (B,1,N)\n    A = torch.fft.rfft(a, dim=-1)\n    B = torch.fft.rfft(b, dim=-1)\n    return torch.mean((A.real - B.real)**2 + (A.imag - B.imag)**2)\n\ndef train_denoiser(templates, train_meta_path, train_dir, leads=LEADS,\n                   epochs=1, lr=1e-3, use_tpl=True, device=\"cpu\",\n                   freq_w=0.2, paper_speed_fn=detect_paper_speed_label):\n    model = Denoiser1D(in_ch=2 if use_tpl else 1).to(device)\n    opt = torch.optim.Adam(model.parameters(), lr=lr)\n    huber = nn.SmoothL1Loss()\n\n    meta = pd.read_csv(train_meta_path)\n    ids = meta[\"id\"].astype(int).tolist()\n\n    model.train()\n    for ep in range(epochs):\n        np.random.shuffle(ids)\n        running = 0.0; steps = 0\n        for rid in tqdm(ids, desc=f\"Denoiser epoch {ep+1}/{epochs}\"):\n            csvp = os.path.join(train_dir, str(rid), f\"{rid}.csv\")\n            imgs = sorted(glob.glob(os.path.join(train_dir, str(rid), f\"{rid}-*.png\")))\n            if not (os.path.exists(csvp) and imgs): continue\n            df  = pd.read_csv(csvp)\n            fs  = int(meta.loc[meta['id']==rid,'fs'].iloc[0])\n            img = cv2.imread(imgs[0])\n            paper_speed = paper_speed_fn(img)\n            panels = split_3x4_panels(img)\n\n            for ld in leads:\n                if ld not in df.columns or ld not in panels: continue\n                y_gt = df[ld].dropna().to_numpy(np.float32)\n                n    = len(y_gt)\n                # 파이프라인 복원\n                y_mix, conf, bpm = panel_to_series(panels[ld], fs, n,\n                                                   templates.get(ld, templates['II']),\n                                                   debug_path=None, paper_speed=paper_speed)\n                # 입력/타깃 표준화\n                yy = zscore(y_gt).astype(np.float32)\n                xx = zscore(y_mix).astype(np.float32)\n                if use_tpl:\n                    ytpl = zscore(template_series(templates.get(ld, templates['II']), fs, n, bpm, 1.0)).astype(np.float32)\n                    x_np = np.stack([xx, ytpl])[None,...]   # (1,2,n)\n                else:\n                    x_np = xx[None,None,:]                  # (1,1,n)\n\n                t_np = yy[None,None,:]                      # (1,1,n)\n\n                x = torch.tensor(x_np, dtype=torch.float32, device=device)\n                t = torch.tensor(t_np, dtype=torch.float32, device=device)\n\n                opt.zero_grad()\n                y_hat = model(x)               # (1,1,n)\n                loss  = huber(y_hat, t) + freq_w * spectral_mse(y_hat, t)\n                loss.backward()\n                opt.step()\n\n                running += float(loss.detach().cpu().item()); steps += 1\n\n        print(f\"[Denoiser] epoch {ep+1}: avg loss={running/max(1,steps):.5f}\")\n\n    # 저장\n    os.makedirs(os.path.dirname(DENOISER_PATH), exist_ok=True)\n    torch.save({\"state_dict\": model.state_dict(),\n                \"in_ch\": 2 if use_tpl else 1}, DENOISER_PATH)\n    print(f\"[Denoiser] saved -> {DENOISER_PATH}\")\n    return model\n\ndef load_denoiser(path=DENOISER_PATH, device=\"cpu\"):\n    if not TORCH_OK or not os.path.exists(path):\n        return None\n    ckpt = torch.load(path, map_location=device)\n    in_ch = ckpt.get(\"in_ch\", 2)\n    model = Denoiser1D(in_ch=in_ch).to(device)\n    model.load_state_dict(ckpt[\"state_dict\"], strict=True)\n    model.eval()\n    return model\n\ndef apply_denoiser(y_mix, y_tpl, model, device=\"cpu\"):\n    if (model is None) or (not DENOISER_ENABLE) or (not TORCH_OK):\n        return y_mix\n    # 표준화 입력 → 보정 → 제출 스케일로 재스케일\n    xx = zscore(y_mix).astype(np.float32)\n    if y_tpl is not None:\n        xt = zscore(y_tpl).astype(np.float32)\n        x_np = np.stack([xx, xt])[None,...]   # (1,2,n)\n    else:\n        x_np = xx[None,None,:]\n    x = torch.tensor(x_np, dtype=torch.float32, device=device)\n    with torch.no_grad():\n        y_hat = model(x).cpu().numpy()[0,0]\n    return rescale_range(y_hat, MIN_VAL, MAX_VAL)\n\n# -----------------------\n# 6) Train 시각화 (GT vs Pred) with optional denoiser\n# -----------------------\ndef plot_train_gt_vs_pred(rec_id, templates, leads=('II','V2','V5'), denoiser=None, device=\"cpu\"):\n    rid = str(int(rec_id))\n    csvp = os.path.join(TRAIN_DIR, rid, f\"{rid}.csv\")\n    imgs = sorted(glob.glob(os.path.join(TRAIN_DIR, rid, f\"{rid}-*.png\")))\n    if not os.path.exists(csvp) or not imgs:\n        print(f\"[Skip] {rid} (csv/image missing)\")\n        return\n    df = pd.read_csv(csvp)\n    fs = int(pd.read_csv(TRAIN_CSV).loc[lambda d: d['id']==int(rid),'fs'].iloc[0])\n\n    img = cv2.imread(imgs[0])\n    paper_speed = detect_paper_speed_label(img)\n    panels = split_3x4_panels(img)\n\n    rows = len(leads)\n    fig, axes = plt.subplots(rows, 1, figsize=(12, 2.8*rows), squeeze=False)\n    axes = axes.ravel()\n\n    for i, ld in enumerate(leads):\n        if ld not in df.columns or ld not in panels:\n            axes[i].set_axis_off(); continue\n        y_gt = df[ld].dropna().to_numpy(np.float32)\n        n = len(y_gt)\n        y_pred, conf, bpm = panel_to_series(\n            panels[ld], fs, n, templates.get(ld, templates['II']),\n            debug_path=None, paper_speed=paper_speed\n        )\n        # (옵션) 디노이저 적용\n        y_tpl = template_series(templates.get(ld, templates['II']), fs, n, bpm, 1.0)\n        y_pred = apply_denoiser(y_pred, y_tpl if DENOISER_USE_TPL_CH else None, denoiser, device=device)\n\n        yy = zscore(y_gt); pp = zscore(y_pred)\n        pear = np.corrcoef(yy, pp)[0,1] if (np.std(yy)>0 and np.std(pp)>0) else 0.0\n        rmse = float(np.sqrt(np.mean((yy-pp)**2)))\n        ss_tot = float(np.sum((yy - np.mean(yy))**2)); r2 = float(1.0 - np.sum((yy-pp)**2)/ss_tot) if ss_tot>0 else 0.0\n\n        t = np.arange(n)/float(fs)\n        ax = axes[i]\n        ax.plot(t, yy, color='k', lw=1.0, label=f'GT {ld}')\n        ax.plot(t, pp, color='crimson', lw=1.0, label=f'Pred {ld}')\n        ax.set_title(f\"Train ID {rid} — Lead {ld}   r={pear:.3f}  RMSE={rmse:.3f}  R²={r2:.3f}\")\n        ax.set_xlabel(\"Time (s)\"); ax.set_ylabel(\"Norm amp\")\n        ax.grid(True); ax.legend(loc='upper right', fontsize=9)\n\n    plt.tight_layout()\n    outp = os.path.join(VIS_DIR, f\"gt_pred_{rid}.png\")\n    plt.savefig(outp, dpi=200, bbox_inches='tight')\n    print(f\"[Saved] {outp}\")\n    plt.show()\n\n# -----------------------\n# 7) TEST Inference & Submission (with optional denoiser)\n# -----------------------\ndef run_test_submission(templates, denoiser=None, device=\"cpu\"):\n    test = pd.read_csv(TEST_CSV)\n    id2img = {}\n    for tid in test['id'].unique():\n        p = sorted(glob.glob(os.path.join(TEST_DIR, str(int(tid)), f\"{int(tid)}-*.png\")))\n        id2img[int(tid)] = p[0] if p else None\n\n    sub_rows = []\n    for r in tqdm(test.itertuples(index=False), total=len(test), desc=\"Predict(test)\"):\n        base_id, lead, fs, n = int(r.id), str(r.lead), int(r.fs), int(r.number_of_rows)\n        imgp = id2img.get(base_id, None)\n        if (imgp is None) or (not os.path.exists(imgp)) or (lead not in LEADS):\n            y = template_series(templates.get(lead, templates['II']), fs, n, bpm=75.0, amp=1.0)\n            y = rescale_range(y, MIN_VAL, MAX_VAL)\n        else:\n            img = cv2.imread(imgp)\n            paper_speed = detect_paper_speed_label(img)\n            panels = split_3x4_panels(img)\n            if lead not in panels:\n                y = template_series(templates.get(lead, templates['II']), fs, n, bpm=75.0, amp=1.0)\n                y = rescale_range(y, MIN_VAL, MAX_VAL)\n            else:\n                y, conf, bpm = panel_to_series(\n                    panels[lead], fs, n, templates.get(lead, templates['II']),\n                    debug_path=None, paper_speed=paper_speed\n                )\n                # (옵션) 디노이저\n                if DENOISER_ENABLE and TORCH_OK:\n                    y_tpl = template_series(templates.get(lead, templates['II']), fs, n, bpm, 1.0)\n                    y = apply_denoiser(y, y_tpl if DENOISER_USE_TPL_CH else None, denoiser, device=device)\n\n        ids = [f\"{base_id}_{i}_{lead}\" for i in range(n)]\n        sub_rows.extend(zip(ids, y.tolist()))\n    sub = pd.DataFrame(sub_rows, columns=['id','value'])\n    sub.to_csv(SUBMISSION_CSV, index=False)\n    print(sub.head(10))\n    print(f\"[OK] Wrote submission: {SUBMISSION_CSV} (rows={len(sub)})\")\n    return sub\n\n# =======================\n# RUN\n# =======================\n# 1) 템플릿을 항상 새로 학습\ntemplates, used = build_lead_templates_beatwise(TRAIN_CSV, TRAIN_DIR, leads=LEADS)\n# (옵션) 저장만 수행 (로딩 안 함)\nnp.savez_compressed(TEMPLATE_NPZ, **templates)\nprint(\"[OK] (re)built templates and saved ->\", TEMPLATE_NPZ)\nfor ld in LEADS:\n    print(f\"  {ld:>3}: beats={used.get(ld,0)} tpl_len={len(templates[ld])}\")\n\n# 2) (옵션) 디노이저 학습 또는 로드\ndevice = \"cuda\" if TORCH_OK and (os.getenv(\"CUDA_VISIBLE_DEVICES\",\"\") != \"\") and (getattr(__import__('torch'), 'cuda', None) and __import__('torch').cuda.is_available()) else \"cpu\"\ndenoiser = None\nif DENOISER_ENABLE and TORCH_OK:\n    if DENOISER_TRAIN:\n        denoiser = train_denoiser(templates, TRAIN_CSV, TRAIN_DIR, leads=LEADS,\n                                  epochs=DENOISER_EPOCHS, lr=DENOISER_LR,\n                                  use_tpl=DENOISER_USE_TPL_CH, device=device,\n                                  freq_w=DENOISER_FREQ_LOSS_W)\n    else:\n        denoiser = load_denoiser(DENOISER_PATH, device=device)\n        if denoiser is None:\n            print(\"[Denoiser] no weights found; run with ECG_DENOISER_TRAIN=1 to train.\")\n\n# 3) Train 시각화 (예시 2개)\ntrain_meta = pd.read_csv(TRAIN_CSV)\nexample_ids = [int(train_meta.iloc[i]['id']) for i in range(min(2, len(train_meta)))]\nfor rid in example_ids:\n    plot_train_gt_vs_pred(rid, templates, leads=('II','V2','V5'), denoiser=denoiser, device=device)\n\n# 4) Test → submission.csv\n_ = run_test_submission(templates, denoiser=denoiser, device=device)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-24T07:57:44.779380Z","iopub.execute_input":"2025-10-24T07:57:44.780774Z","iopub.status.idle":"2025-10-24T07:58:28.039094Z","shell.execute_reply.started":"2025-10-24T07:57:44.780736Z","shell.execute_reply":"2025-10-24T07:58:28.038265Z"}},"outputs":[],"execution_count":null}]}