{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# =========================================================\n# PhysioNet ECG Image Digitization — v4m+ (safe template-first blend)\n# ResNet34-H + FiLM + Temporal Head  (3-channel input)\n# ROI: heuristic(train) + UNet(train/val), test = grid-only (stable)\n# Tilt-aware TTA, multi-scale freq loss, cosine LR, motion/perspective augs\n# Safe blend: lag-aligned + sign-check; template default; model weight gated by NCC\n# Output: /kaggle/working/submission.parquet  &  submission.csv\n# Includes: submission validator + offline LB estimate (official metric)\n# =========================================================\nimport os, gc, cv2, glob, math, random, warnings, pickle\nwarnings.filterwarnings(\"ignore\")\nimport numpy as np, pandas as pd\nfrom tqdm import tqdm\nimport scipy.signal, scipy.optimize\n\nimport torch, torch.nn as nn, torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.model_selection import train_test_split\n\n# ---------------------------\n# Config (safer blending + longer train)\n# ---------------------------\nclass CFG:\n    DATA_DIR = \"/kaggle/input/physionet-ecg-image-digitization\"\n    SEED = 1337\n    IMG_H = 352\n    IMG_W = 1280\n    T_TRAIN = 1536\n    BATCH_SIZE = 2\n    EPOCHS_MAIN = 10                 # ↑ train a bit longer\n    LR = 3e-4\n    WD = 1e-4\n    NUM_WORKERS = 2\n    AMP = True\n    CHANNELS_LAST = True\n    VALID_FRAC = 0.05\n    SAVE_BEST = True\n\n    # ROI modes\n    ROI_MODE_TRAIN = \"heuristic\"\n    ROI_MODE_VAL   = \"unet\"\n    ROI_MODE_TESTS = [\"grid\"]        # ← safest & most stable for LB\n\n    GRID_USE_RHYTHM_FOR_II = True\n    TRAIN_ROI_UNET = True\n    ROI_EPOCHS = 1\n    ROI_IN_TRAIN = True\n    LEAD_BAND_FRAC = 0.06\n    NMS_MIN_DIST_FRAC = 0.04\n    ROI_BAND_MARGIN = 0.015\n\n    # Augs\n    AUG_BRIGHTNESS = 0.12\n    AUG_CONTRAST = 0.12\n    AUG_ROT_DEG = 0.7\n    AUG_PERSPECTIVE = 0.015\n    AUG_MOTION_BLUR_P = 0.25\n\n    # Inference\n    TTA_BASE_DEGS = [-0.4, 0.0, 0.4]\n    TTA_EXTRA = True\n\n    # Filters (we’ll fix cutoff to 15 Hz)\n    BANDPASS_LOW = 0.5\n    BANDPASS_HIGH = 40.0\n    BANDPASS_ORDER = 3\n\n    # Templates\n    TEMPLATE_LEN = 500\n    TEMPLATE_MAX_SAMPLES = 4000\n    TEMPLATE_CACHE = \"/kaggle/working/lead_templates.pkl\"\n\n    # Blend (safe: template default)\n    ENSEMBLE_BETA_BASE = 0.0         # ← no unconditional model weight\n    ENSEMBLE_BETA_GAIN = 0.50        # scales with correlation\n    ENSEMBLE_BETA_CLIP = (0.00, 0.65)\n\n    OUT_MIN = 0.00\n    OUT_MAX = 0.09                   # ← match stronger baseline\n\n    # Train knobs\n    GRAD_ACCUM_TARGET = 8            # larger effective batch\n\n    # Local sanity\n    DO_LOCAL_SNR = True\n\ndef set_seed(seed=1337):\n    random.seed(seed); np.random.seed(seed)\n    torch.manual_seed(seed); torch.cuda.manual_seed_all(seed)\nset_seed(CFG.SEED)\n\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\ntorch.backends.cuda.matmul.allow_tf32 = True\ntorch.backends.cudnn.benchmark = True\n\ntrain_csv = pd.read_csv(f\"{CFG.DATA_DIR}/train.csv\")\ntest_csv  = pd.read_csv(f\"{CFG.DATA_DIR}/test.csv\")\nLEADS = ['I','II','III','aVR','aVL','aVF','V1','V2','V3','V4','V5','V6']\nLEAD_TO_IDX = {l:i for i,l in enumerate(LEADS)}\nVARIANTS = [\"0001\",\"0003\",\"0004\",\"0005\",\"0006\",\"0009\",\"0010\",\"0011\",\"0012\"]\nprint(\"Device:\", device, \"GPUs:\", torch.cuda.device_count())\nprint(\"Train records:\", len(train_csv), \" Test records:\", len(test_csv))\n\n# ---------------------------\n# IO / preprocess\n# ---------------------------\ndef pick_image_path(base_id, prefer_first=True):\n    folder = f\"{CFG.DATA_DIR}/train/{base_id}\"\n    if not os.path.isdir(folder): return None\n    paths = [f\"{folder}/{base_id}-{s}.png\" for s in VARIANTS if os.path.exists(f\"{folder}/{base_id}-{s}.png\")]\n    if not paths: paths = glob.glob(f\"{folder}/*.png\")\n    if not paths: return None\n    return paths[0] if prefer_first else random.choice(paths)\n\ndef apply_motion_blur(img01):\n    if np.random.rand() > CFG.AUG_MOTION_BLUR_P: return img01\n    k = np.random.choice([3,5,7])\n    kernel = np.zeros((k, k), dtype=np.float32)\n    idx = k//2\n    kernel[idx, :] = 1.0 / k\n    return cv2.filter2D(img01, -1, kernel)\n\ndef perspective_jitter(img01):\n    h,w = img01.shape\n    d = CFG.AUG_PERSPECTIVE\n    if d <= 1e-6: return img01\n    dx, dy = d*w, d*h\n    src = np.float32([[0,0],[w,0],[w,h],[0,h]])\n    dst = np.float32([\n        [np.random.uniform(-dx, dx), np.random.uniform(-dy, dy)],\n        [w+np.random.uniform(-dx, dx), np.random.uniform(-dy, dy)],\n        [w+np.random.uniform(-dx, dx), h+np.random.uniform(-dy, dy)],\n        [np.random.uniform(-dx, dx), h+np.random.uniform(-dy, dy)]\n    ])\n    M = cv2.getPerspectiveTransform(src, dst)\n    out = cv2.warpPerspective((img01*255).astype(np.uint8), M, (w,h),\n                              flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REPLICATE)\n    return out.astype(np.float32)/255.0\n\ndef load_img_gray(path):\n    img = cv2.imread(path, cv2.IMREAD_GRAYSCALE)\n    img = cv2.resize(img, (CFG.IMG_W, CFG.IMG_H), interpolation=cv2.INTER_AREA)\n    img = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8)).apply(img)\n    return img.astype(np.float32) / 255.0\n\ndef rotate_image(img01, deg):\n    if abs(deg) < 1e-6: return img01\n    M = cv2.getRotationMatrix2D((CFG.IMG_W/2, CFG.IMG_H/2), deg, 1.0)\n    rot = cv2.warpAffine((img01*255).astype(np.uint8), M, (CFG.IMG_W, CFG.IMG_H),\n                         flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REPLICATE)\n    return rot.astype(np.float32)/255.0\n\ndef to_three_channel(img01):\n    sx = cv2.Sobel(img01, cv2.CV_32F, 1, 0, ksize=3)\n    sy = cv2.Sobel(img01, cv2.CV_32F, 0, 1, ksize=3)\n    sx = (sx - sx.min()) / (sx.max() - sx.min() + 1e-6)\n    sy = (sy - sy.min()) / (sy.max() - sy.min() + 1e-6)\n    return np.stack([img01, sx, sy], axis=0)  # (3,H,W)\n\ndef random_augment(img01):\n    if CFG.AUG_BRIGHTNESS>0:\n        b=np.random.uniform(-CFG.AUG_BRIGHTNESS, CFG.AUG_BRIGHTNESS); img01=np.clip(img01+b,0,1)\n    if CFG.AUG_CONTRAST>0:\n        c=np.random.uniform(1-CFG.AUG_CONTRAST,1+CFG.AUG_CONTRAST); m=img01.mean(); img01=np.clip((img01-m)*c+m,0,1)\n    if CFG.AUG_ROT_DEG>0:\n        deg=np.random.uniform(-CFG.AUG_ROT_DEG, CFG.AUG_ROT_DEG); img01=rotate_image(img01,deg)\n    img01 = perspective_jitter(img01)\n    img01 = apply_motion_blur(img01)\n    return img01\n\ndef resample_1d(x, new_len):\n    new_len = int(new_len)\n    if new_len <= 0:\n        return np.zeros(0, dtype=np.float32)\n    x = np.asarray(x, dtype=np.float32).reshape(-1)\n    n = x.shape[0]\n    if n == 0:\n        return np.zeros(new_len, dtype=np.float32)\n    if n == 1:\n        return np.full(new_len, float(x[0]), dtype=np.float32)\n    xp = np.linspace(0.0, n - 1.0, num=n, dtype=np.float64)\n    fp = np.nan_to_num(x, nan=0.0, posinf=0.0, neginf=0.0).astype(np.float64)\n    xnew = np.linspace(0.0, n - 1.0, num=new_len, dtype=np.float64)\n    y = np.interp(xnew, xp, fp).astype(np.float32)\n    return y\n\n# Tiny tilt estimator (for grid-aligned pages)\ndef estimate_tilt_deg(img01):\n    im = (img01*255).astype(np.uint8)\n    edges = cv2.Canny(im, 50, 150, apertureSize=3)\n    lines = cv2.HoughLines(edges, 1, np.pi/180, threshold=200)\n    if lines is None: return 0.0\n    angs=[]\n    for rho,theta in lines[:,0,:]:\n        a = (theta - np.pi/2) * 180/np.pi  # horizontal grid → near 0°\n        if -3.0 <= a <= 3.0: angs.append(a)\n    if not angs: return 0.0\n    a = float(np.median(angs))\n    return float(np.clip(a, -1.2, 1.2))\n\n# ---------------------------\n# ROI helpers (heuristic / UNet / grid)\n# ---------------------------\ndef heuristic_peak_bands(img01):\n    h,w=img01.shape\n    edges=cv2.Canny((img01*255).astype(np.uint8),40,120)\n    proj=edges.sum(axis=1).astype(np.float32)\n    p=(proj-proj.min())/(proj.max()-proj.min()+1e-6)\n    min_dist=max(1,int(CFG.NMS_MIN_DIST_FRAC*h))\n    peaks=[]\n    for y in range(h):\n        l=max(0,y-min_dist); r=min(h,y+min_dist+1)\n        if p[y]==p[l:r].max() and p[y]>0.1: peaks.append((p[y],y))\n    peaks=sorted(peaks,key=lambda t:t[0],reverse=True)\n    bands=[]; used=np.zeros(h,bool); band_half=int(CFG.LEAD_BAND_FRAC*h/2)\n    for _,y in peaks:\n        if used[max(0,y-min_dist):min(h,y+min_dist+1)].any(): continue\n        y0=max(0,y-band_half-int(CFG.ROI_BAND_MARGIN*h)); y1=min(h,y+band_half+int(CFG.ROI_BAND_MARGIN*h))\n        bands.append((y0,y1)); used[y0:y1]=True\n        if len(bands)==12: break\n    bands=sorted(bands,key=lambda b:(b[0]+b[1])/2)\n    while len(bands)<12: bands.append((int(len(bands)*h/12), int((len(bands)+1)*h/12)))\n    return bands[:12]\n\ndef crops_from_heuristic(img01):\n    return [img01[y0:y1,:] for (y0,y1) in heuristic_peak_bands(img01)]\n\ndef stack_crop_to_x3(crop):\n    Hc=max(48,int(CFG.LEAD_BAND_FRAC*CFG.IMG_H))\n    c=cv2.resize(crop,(CFG.IMG_W,Hc),interpolation=cv2.INTER_AREA)\n    return to_three_channel(c)  # (3,Hc,W)\n\nclass TinyUNet(nn.Module):\n    def __init__(self,in_ch=1,out_ch=12,base=16):\n        super().__init__()\n        self.e1=nn.Sequential(nn.Conv2d(in_ch,base,3,1,1),nn.ReLU(True),nn.Conv2d(base,base,3,1,1),nn.ReLU(True))\n        self.e2=nn.Sequential(nn.Conv2d(base,base*2,3,2,1),nn.ReLU(True))\n        self.e3=nn.Sequential(nn.Conv2d(base*2,base*4,3,2,1),nn.ReLU(True))\n        self.d2=nn.Sequential(nn.ConvTranspose2d(base*4,base*2,2,2),nn.ReLU(True))\n        self.d1=nn.Sequential(nn.ConvTranspose2d(base*2,base,2,2),nn.ReLU(True))\n        self.outc=nn.Conv2d(base,out_ch,1)\n    def forward(self,x):\n        x1=self.e1(x); x2=self.e2(x1); x3=self.e3(x2)\n        y=self.d2(x3); y=self.d1(y); return self.outc(y)\n\ndef crops_from_mask(img01, mask12):\n    H,W=img01.shape; crops=[]\n    for i in range(12):\n        m=mask12[i]; col=m.mean(axis=1)\n        y=int(np.clip(np.average(np.arange(H),weights=(col+1e-6)),0,H-1))\n        band=max(4,int(CFG.LEAD_BAND_FRAC*H)); y0=max(0,y-band//2); y1=min(H,y0+band)\n        if y1-y0<band: y0=max(0,y1-band)\n        crops.append(img01[y0:y1,:])\n    return crops\n\nGRID_ROWS = [\n    ['I','II','III'],\n    ['aVR','aVL','aVF'],\n    ['V1','V2','V3'],\n    ['V4','V5','V6'],\n]\ndef crops_from_grid(img01, use_rhythm_for_ii=True):\n    H, W = img01.shape\n    rhythm_h = int(0.22 * H); rhythm_h = max(int(0.15*H), min(int(0.30*H), rhythm_h))\n    grid_h = H - rhythm_h\n    row_edges = np.linspace(0, grid_h, 5).astype(int)\n    col_edges = np.linspace(0, W, 4).astype(int)\n    lead_to_crop = {}\n    for r in range(4):\n        y0, y1 = row_edges[r], row_edges[r+1]\n        for c in range(3):\n            x0, x1 = col_edges[c], col_edges[c+1]\n            lead = GRID_ROWS[r][c]\n            lead_to_crop[lead] = img01[y0:y1, x0:x1]\n    rhythm = img01[grid_h:H, :]\n    if use_rhythm_for_ii:\n        lead_to_crop['II'] = rhythm\n    return [lead_to_crop[l] for l in LEADS]\n\ndef get_crops(img01, roi_mode=\"heuristic\", roi_net=None):\n    if roi_mode == \"unet\" and (roi_net is not None):\n        with torch.no_grad():\n            # CPU tensor; DO NOT .to(device) in a worker\n            x = torch.from_numpy(img01).unsqueeze(0).unsqueeze(0).float()\n            if CFG.CHANNELS_LAST:\n                x = x.to(memory_format=torch.channels_last)\n            m = roi_net(x).sigmoid()[0].detach().numpy()\n        return crops_from_mask(img01, m)\n    elif roi_mode == \"grid\":\n        return crops_from_grid(img01, use_rhythm_for_ii=CFG.GRID_USE_RHYTHM_FOR_II)\n    else:\n        return crops_from_heuristic(img01)\n        \n# ---------------------------\n# Model: ResNet34-H + FiLM + Temporal\n# ---------------------------\nclass BasicBlockH(nn.Module):\n    expansion=1\n    def __init__(self,inplanes,planes,stride=(1,1),downsample=None):\n        super().__init__()\n        self.conv1=nn.Conv2d(inplanes,planes,3,stride=stride,padding=1,bias=False)\n        self.bn1=nn.BatchNorm2d(planes); self.relu=nn.ReLU(True)\n        self.conv2=nn.Conv2d(planes,planes,3,1,1,bias=False); self.bn2=nn.BatchNorm2d(planes)\n        self.downsample=downsample\n    def forward(self,x):\n        idt=x\n        out=self.relu(self.bn1(self.conv1(x)))\n        out=self.bn2(self.conv2(out))\n        if self.downsample is not None: idt=self.downsample(x)\n        return self.relu(out+idt)\n\ndef _make_layer_h(inplanes,planes,blocks,stride_h):\n    down=None\n    if stride_h!=1 or inplanes!=planes:\n        down=nn.Sequential(nn.Conv2d(inplanes,planes,1,(stride_h,1),bias=False), nn.BatchNorm2d(planes))\n    layers=[BasicBlockH(inplanes,planes,(stride_h,1),down)]\n    for _ in range(1,blocks): layers.append(BasicBlockH(planes,planes,(1,1)))\n    return nn.Sequential(*layers)\n\nclass ResNet34H(nn.Module):\n    def __init__(self,in_ch=3):   # ← 3 channels\n        super().__init__()\n        self.conv1=nn.Conv2d(in_ch,64,7,(2,1),3,bias=False)\n        self.bn1=nn.BatchNorm2d(64); self.relu=nn.ReLU(True)\n        self.maxpool=nn.MaxPool2d((3,3),(2,1),1)\n        self.layer1=_make_layer_h(64,64,3,1)\n        self.layer2=_make_layer_h(64,128,4,2)\n        self.layer3=_make_layer_h(128,256,6,2)\n        self.layer4=_make_layer_h(256,512,3,2)\n        self.proj=nn.Conv2d(512,256,1)\n    def forward(self,x):\n        x=self.relu(self.bn1(self.conv1(x))); x=self.maxpool(x)\n        x=self.layer1(x); x=self.layer2(x); x=self.layer3(x); x=self.layer4(x)\n        x=self.proj(x).mean(dim=2)  # (B,256,W)\n        return x\n\nclass DilatedTemporalBlock(nn.Module):\n    def __init__(self,c,d):\n        super().__init__(); self.conv=nn.Conv1d(c,c,3,padding=d,dilation=d,bias=False)\n        self.bn=nn.BatchNorm1d(c); self.act=nn.ReLU(True)\n    def forward(self,x): return x + self.act(self.bn(self.conv(x)))\n\nclass TemporalHead(nn.Module):\n    def __init__(self,c=256,h=128):\n        super().__init__()\n        self.pre=nn.Sequential(nn.Conv1d(c,c,7,padding=3,bias=False), nn.BatchNorm1d(c), nn.ReLU(True))\n        self.d1=DilatedTemporalBlock(c,1); self.d2=DilatedTemporalBlock(c,2)\n        self.d4=DilatedTemporalBlock(c,4); self.d8=DilatedTemporalBlock(c,8)\n        self.rnn=nn.GRU(c,h,batch_first=True,bidirectional=True)\n        self.fuse=nn.Conv1d(2*h,c,1)\n    def forward(self,x):\n        x=self.pre(x); x=self.d1(x); x=self.d2(x); x=self.d4(x); x=self.d8(x)\n        seq=x.transpose(1,2); self.rnn.flatten_parameters()\n        y,_=self.rnn(seq); return self.fuse(y.transpose(1,2))\n\nclass FiLMLead(nn.Module):\n    def __init__(self,n=12,d=64,c=256):\n        super().__init__(); self.emb=nn.Embedding(n,d); self.g=nn.Linear(d,c); self.b=nn.Linear(d,c)\n    def forward(self,x,lead_idx):\n        e=self.emb(lead_idx); gamma=self.g(e).unsqueeze(-1); beta=self.b(e).unsqueeze(-1)\n        return x*(1+gamma)+beta\n\nclass Head1D(nn.Module):\n    def __init__(self,c=256): super().__init__(); self.net=nn.Sequential(nn.Conv1d(c,64,3,1,1), nn.ReLU(True), nn.Conv1d(64,1,3,1,1))\n    def forward(self,x): return self.net(x).squeeze(1)\n\nclass ECGNetV4M(nn.Module):\n    def __init__(self,in_ch=3):\n        super().__init__(); self.enc=ResNet34H(in_ch); self.tem=TemporalHead(); self.film=FiLMLead(); self.head=Head1D()\n    def forward(self,x3,lead_idx):\n        f=self.enc(x3); t=self.tem(f); t=self.film(t,lead_idx); return self.head(t)\n\n# --- robust FFT + multiscale\ndef fft_mag(x: torch.Tensor) -> torch.Tensor:\n    with torch.cuda.amp.autocast(enabled=False):\n        x32 = x.float()\n        X = torch.fft.rfft(x32, dim=-1)\n        mag = torch.abs(X)\n        return mag\n\ndef ms_downsample(x, k):\n    # avgpool downsample by k along last dim\n    return F.avg_pool1d(x.unsqueeze(1), kernel_size=k, stride=k, ceil_mode=False).squeeze(1)\n\n# --- safer CSV loader ---\ndef load_wave_csv_matrix(path):\n    df = pd.read_csv(path)\n    vals = []\n    for c in LEADS:\n        v = pd.to_numeric(df[c], errors='coerce').values.astype(np.float32)\n        v = np.nan_to_num(v, nan=0.0, posinf=0.0, neginf=0.0)\n        vals.append(v)\n    arr = np.stack(vals, axis=0)\n    return arr\n\ndef first_deriv(x): return x[...,1:]-x[...,:-1]\n\ndef main_loss(pred,target):\n    # base L1 + derivative + multiscale spectral L1 (three scales)\n    T=min(pred.shape[-1],target.shape[-1]); p=pred[...,:T]; t=target[...,:T]\n    l = F.l1_loss(p,t) + 0.05*F.l1_loss(first_deriv(p),first_deriv(t))\n    l += 0.10*F.l1_loss(fft_mag(p), fft_mag(t))\n    p2,t2 = ms_downsample(p,2), ms_downsample(t,2)\n    p4,t4 = ms_downsample(p,4), ms_downsample(t,4)\n    l += 0.06*F.l1_loss(fft_mag(p2), fft_mag(t2))\n    l += 0.04*F.l1_loss(fft_mag(p4), fft_mag(t4))\n    return l\n\n# ---------------------------\n# Dataset\n# ---------------------------\nclass LeadDataset(Dataset):\n    def __init__(self, ids, roi_net=None, roi_mode=\"heuristic\", augment=False, mix_unet=False):\n        self.ids = list(ids); self.roi_net = roi_net; self.roi_mode = roi_mode; self.augment = augment\n        self.mix_unet = mix_unet\n        self.samples = [(id_, li) for id_ in self.ids for li in range(12)]\n    def __len__(self): return len(self.samples)\n    def __getitem__(self, idx):\n        base_id, li = self.samples[idx]\n        ipath = pick_image_path(base_id, prefer_first=True)\n        if ipath is None: raise FileNotFoundError(f\"no image for {base_id}\")\n        img01 = load_img_gray(ipath)\n        if self.augment: img01 = random_augment(img01)\n        # occasionally use UNet ROI even during \"heuristic\" train\n        roi_mode_eff = self.roi_mode\n        if self.mix_unet and (self.roi_net is not None) and (np.random.rand() < 0.5):\n            roi_mode_eff = \"unet\"\n        crops = get_crops(img01, roi_mode=roi_mode_eff, roi_net=self.roi_net if roi_mode_eff==\"unet\" else None)\n        crop = crops[li]\n        x3 = stack_crop_to_x3(crop)                      # (3,Hc,W)\n        csv_path=f\"{CFG.DATA_DIR}/train/{base_id}/{base_id}.csv\"\n        M=load_wave_csv_matrix(csv_path)\n        y=resample_1d(M[li], CFG.T_TRAIN)\n        y = np.nan_to_num(y, nan=0.0, posinf=0.0, neginf=0.0)\n        return (torch.tensor(x3, dtype=torch.float32),\n                torch.tensor(y, dtype=torch.float32),\n                torch.tensor(li, dtype=torch.long))\n\n# ---------------------------\n# Metric (exact) — used for quick local sanity on Lead II\n# ---------------------------\nMAX_TIME_SHIFT = 0.2; PERFECT_SCORE = 384\nclass ParticipantVisibleError(Exception): pass\ndef compute_power_metric(label, prediction):\n    if label.ndim!=1 or prediction.ndim!=1: raise ParticipantVisibleError('Inputs must be 1D.')\n    finite_mask=np.isfinite(prediction)\n    if not np.any(finite_mask): raise ParticipantVisibleError(\"prediction all NaN/inf\")\n    prediction=prediction.copy(); prediction[~np.isfinite(prediction)]=0\n    noise=label-prediction; return np.sum(label**2), np.sum(noise**2)\ndef compute_snr_metric(signal, noise):\n    if noise==0: return PERFECT_SCORE\n    if signal==0: return 0\n    return min((signal/noise), PERFECT_SCORE)\ndef align_signals_metric(label,pred,max_shift=float('inf')):\n    if np.any(~np.isfinite(label)): raise ParticipantVisibleError('label finite')\n    if np.sum(np.isfinite(pred))==0: raise ParticipantVisibleError('prediction finite')\n    la=np.asarray(label,dtype=np.float64); pr=np.asarray(pred,dtype=np.float64)\n    lc=la-np.mean(la); pc=pr-np.mean(pr)\n    corr=scipy.signal.correlate(lc,pc,mode='full'); nl=la.size; npd=pr.size\n    lags=scipy.signal.correlation_lags(nl,npd,mode='full'); valid=(lags>=-max_shift)&(lags<=max_shift)\n    maxc=np.nanmax(corr[valid]); allmax=np.flatnonzero(corr==maxc); best=min(allmax,key=lambda i:abs(lags[i])); shift=lags[best]\n    sp=max(shift,0); s=max(-shift,0); e=min(nl-shift,npd); ep=max(nl-npd-shift,0)\n    aligned=np.concatenate((np.full(sp,np.nan), pr[s:e], np.full(ep,np.nan)))\n    def obj(v): return np.nansum((la-(aligned-v))**2)\n    if np.any(np.isfinite(la)&np.isfinite(aligned)):\n        res=scipy.optimize.minimize_scalar(obj,method='Brent'); aligned-=res.x\n    return aligned\ndef _calc_image_score(group):\n    ufs=group['fs'].unique()\n    if len(ufs)!=1: raise ParticipantVisibleError('fs consistent')\n    fs=ufs[0]\n    if fs!=int(len(group[group['lead']=='II'])/10): raise ParticipantVisibleError('len check')\n    sum_sig=sum_noise=0.0\n    for lead in LEADS:\n        sub=group[group['lead']==lead]\n        label=sub['value_true'].values; pred=sub['value_pred'].values\n        aligned=align_signals_metric(label,pred,int(fs*MAX_TIME_SHIFT))\n        ps,pn=compute_power_metric(label,aligned); sum_sig+=ps; sum_noise+=pn\n    return compute_snr_metric(sum_sig,sum_noise)\ndef competition_score(solution, submission, row_id_column_name=\"id\"):\n    for df in [solution,submission]:\n        if row_id_column_name not in df.columns: raise ParticipantVisibleError(f\"'{row_id_column_name}' missing\")\n        if df['value'].isna().any(): raise ParticipantVisibleError('NaN in dfs')\n        if not np.isfinite(df['value']).all(): raise ParticipantVisibleError('Inf in dfs')\n    sub=submission[['id','value']]\n    m=pd.merge(solution, sub, on=row_id_column_name, suffixes=('_true','_pred'))\n    m['image_id']=m[row_id_column_name].str.split('_').str[0]\n    m['row_id']=m[row_id_column_name].str.split('_').str[1].astype('int64')\n    m['lead']=m[row_id_column_name].str.split('_').str[2]\n    m.sort_values(by=['image_id','row_id','lead'], inplace=True)\n    img_scores=m.groupby('image_id').apply(_calc_image_score, include_groups=False)\n    return max(float(10*np.log10(img_scores.mean())), -PERFECT_SCORE)\n\n# ---------------------------\n# ROI quick train (small subset)\n# ---------------------------\nroi_net = None\nif CFG.TRAIN_ROI_UNET:\n    roi_net = TinyUNet().to(device)\n    if CFG.CHANNELS_LAST:\n        roi_net = roi_net.to(memory_format=torch.channels_last)\n    opt_roi = torch.optim.AdamW(roi_net.parameters(), lr=1e-3, weight_decay=1e-4)\n    scaler_roi = torch.cuda.amp.GradScaler(enabled=CFG.AMP)\n\n    def build_pseudo_mask(img01: np.ndarray) -> np.ndarray:\n        H, W = img01.shape\n        mask = np.zeros((12, H, W), dtype=np.float32)\n        for k, (y0, y1) in enumerate(heuristic_peak_bands(img01)):\n            if k >= 12: break\n            mask[k, y0:y1, :] = 1.0\n        if mask.shape[0] < 12:\n            pad = np.zeros((12 - mask.shape[0], H, W), dtype=np.float32)\n            mask = np.concatenate([mask, pad], axis=0)\n        return mask\n\n    ids_all = train_csv[\"id\"].values.copy()\n    np.random.shuffle(ids_all)\n    ids_small = ids_all[:max(200, int(0.15 * len(ids_all)))]\n    bs = 8\n\n    for ep in range(1, CFG.ROI_EPOCHS + 1):\n        roi_net.train()\n        losses = []\n        pbar = tqdm(range(0, len(ids_small), bs), desc=f\"ROI ep{ep}\")\n        for i in pbar:\n            batch_ids = ids_small[i:i+bs]\n            imgs, masks = [], []\n            for bid in batch_ids:\n                p = pick_image_path(bid, prefer_first=True)\n                if p is None: continue\n                img01 = load_img_gray(p)\n                imgs.append(img01)\n                masks.append(build_pseudo_mask(img01))\n            if not imgs: continue\n\n            x = torch.tensor(np.stack(imgs)[:, None, ...],\n                             dtype=torch.float16 if CFG.AMP else torch.float32,\n                             device=device)\n            y = torch.tensor(np.stack(masks),\n                             dtype=torch.float16 if CFG.AMP else torch.float32,\n                             device=device)\n            if CFG.CHANNELS_LAST:\n                x = x.to(memory_format=torch.channels_last)\n\n            opt_roi.zero_grad(set_to_none=True)\n            with torch.cuda.amp.autocast(enabled=CFG.AMP):\n                logits = roi_net(x)\n                loss = F.binary_cross_entropy_with_logits(logits, y)\n            scaler_roi.scale(loss).backward()\n            scaler_roi.step(opt_roi)\n            scaler_roi.update()\n\n            losses.append(loss.item())\n            pbar.set_postfix(loss=np.mean(losses))\n        print(f\"ROI epoch {ep} loss={np.mean(losses):.4f}\")\n\n    torch.cuda.empty_cache()\n\nif roi_net is not None:\n    roi_net.eval()\n    roi_net = roi_net.to('cpu')   # <<< important: avoid CUDA in worker\n\n# ---------------------------\n# Train main model\n# ---------------------------\nmodel=ECGNetV4M().to(device)\nif CFG.CHANNELS_LAST: model=model.to(memory_format=torch.channels_last)\nopt=torch.optim.AdamW(model.parameters(), lr=CFG.LR, weight_decay=CFG.WD)\n# cosine LR\nscheduler = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=CFG.EPOCHS_MAIN, eta_min=CFG.LR*0.05)\nscaler=torch.cuda.amp.GradScaler(enabled=CFG.AMP)\n\ntrain_ids, val_ids = train_test_split(train_csv[\"id\"].values, test_size=CFG.VALID_FRAC, random_state=CFG.SEED)\ntrain_ds=LeadDataset(train_ids, roi_net=roi_net if CFG.ROI_IN_TRAIN else None,\n                     roi_mode=CFG.ROI_MODE_TRAIN, augment=True, mix_unet=CFG.ROI_IN_TRAIN)\nval_roi_mode = \"unet\" if roi_net is not None else \"grid\"\nval_ds  =LeadDataset(val_ids,   roi_net=roi_net, roi_mode=val_roi_mode, augment=False)\n\ntrain_loader=DataLoader(train_ds, batch_size=CFG.BATCH_SIZE, shuffle=True,\n                        num_workers=CFG.NUM_WORKERS, pin_memory=True, persistent_workers=False)\nval_loader = DataLoader(val_ds, batch_size=CFG.BATCH_SIZE, shuffle=False,\n                        num_workers=0, pin_memory=False, persistent_workers=False)\n\naccum_steps=max(1, CFG.GRAD_ACCUM_TARGET//CFG.BATCH_SIZE)\n\ndef train_one_epoch():\n    model.train()\n    losses = []\n    opt.zero_grad(set_to_none=True)\n    pbar = tqdm(enumerate(train_loader), total=len(train_loader), desc=\"Train\", leave=False)\n    for step, (x3, y, li) in pbar:\n        x3 = x3.to(device, non_blocking=True)\n        y  = y.to(device,  non_blocking=True)\n        li = li.to(device, non_blocking=True)\n        if CFG.CHANNELS_LAST:\n            x3 = x3.to(memory_format=torch.channels_last)\n        with torch.cuda.amp.autocast(enabled=CFG.AMP):\n            pred = model(x3, li)\n            if pred.shape[-1] != CFG.T_TRAIN:\n                pred = F.interpolate(pred.unsqueeze(1), size=CFG.T_TRAIN, mode=\"linear\", align_corners=False).squeeze(1)\n        pred = torch.nan_to_num(pred, nan=0.0, posinf=0.0, neginf=0.0)\n        y    = torch.nan_to_num(y,    nan=0.0, posinf=0.0, neginf=0.0)\n        with torch.cuda.amp.autocast(enabled=False):\n            loss = main_loss(pred.float(), y.float()) / accum_steps\n        scaler.scale(loss).backward()\n        if (step + 1) % accum_steps == 0:\n            torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n            scaler.step(opt); scaler.update(); opt.zero_grad(set_to_none=True)\n        losses.append(loss.item() * accum_steps)\n        pbar.set_postfix(loss=float(np.mean(losses)))\n    scheduler.step()\n    return float(np.mean(losses))\n\n@torch.no_grad()\ndef valid_epoch():\n    model.eval()\n    losses = []\n    for x3, y, li in val_loader:\n        x3 = x3.to(device, non_blocking=True)\n        y  = y.to(device,  non_blocking=True)\n        li = li.to(device, non_blocking=True)\n        if CFG.CHANNELS_LAST:\n            x3 = x3.to(memory_format=torch.channels_last)\n        with torch.cuda.amp.autocast(enabled=CFG.AMP):\n            pred = model(x3, li)\n            if pred.shape[-1] != CFG.T_TRAIN:\n                pred = F.interpolate(pred.unsqueeze(1), size=CFG.T_TRAIN, mode=\"linear\", align_corners=False).squeeze(1)\n        pred = torch.nan_to_num(pred, nan=0.0, posinf=0.0, neginf=0.0)\n        y    = torch.nan_to_num(y,    nan=0.0, posinf=0.0, neginf=0.0)\n        with torch.cuda.amp.autocast(enabled=False):\n            loss = main_loss(pred.float(), y.float())\n        losses.append(loss.item())\n    return float(np.mean(losses))\n\nbest_val=1e9; best_path=\"ecgnet_v4m_plus_best.pt\"\nfor ep in range(1, CFG.EPOCHS_MAIN+1):\n    tr=train_one_epoch(); va=valid_epoch()\n    print(f\"[Main {ep}/{CFG.EPOCHS_MAIN}] train={tr:.5f}  valid={va:.5f}\")\n    if CFG.SAVE_BEST and va<best_val:\n        best_val=va\n        torch.save(model.state_dict(), best_path); print(\"  -> saved\", best_path)\n    torch.cuda.empty_cache(); gc.collect()\n\nif os.path.exists(best_path):\n    model.load_state_dict(torch.load(best_path, map_location=device)); print(\"Loaded best weights.\")\ntorch.cuda.empty_cache(); gc.collect()\n\n# ---------------------------\n# Filters / helpers (fixed 15 Hz; stable)\n# ---------------------------\ndef lowpass(sig, fs, order=2, cutoff_hz=15.0):\n    fs = int(fs)\n    sig = np.asarray(sig, dtype=np.float32)\n    nyq = 0.5 * fs\n    w = min(cutoff_hz / nyq, 0.999)\n    b,a = scipy.signal.butter(order, w, btype='low')\n    if len(sig) > 3 * max(len(a), len(b)):\n        try:\n            y = scipy.signal.filtfilt(b, a, sig).astype(np.float32)\n            return np.nan_to_num(y, nan=0.0, posinf=0.0, neginf=0.0)\n        except:\n            pass\n    return np.nan_to_num(sig, nan=0.0, posinf=0.0, neginf=0.0)\n\ndef scale_to_range(x, a=0.0, b=0.09):\n    x = np.asarray(x, dtype=np.float32)\n    xmin, xmax = float(np.min(x)), float(np.max(x))\n    if not np.isfinite(xmin) or not np.isfinite(xmax):\n        return np.full_like(x, (a+b)/2, dtype=np.float32)\n    if xmax - xmin < 1e-8:\n        return np.full_like(x, (a+b)/2, dtype=np.float32)\n    y = (x - xmin) / (xmax - xmin)\n    return (a + y * (b - a)).astype(np.float32)\n\n# ---------------------------\n# NCC + Safe dynamic blend (lag-align + sign-check)\n# ---------------------------\ndef ncc(a, b):\n    a = np.asarray(a, dtype=np.float32); b=np.asarray(b, dtype=np.float32)\n    if len(a)!=len(b):\n        b = resample_1d(b, len(a))\n    ac = a - a.mean(); bc = b - b.mean()\n    da = float(np.sqrt((ac*ac).sum()) + 1e-8)\n    db = float(np.sqrt((bc*bc).sum()) + 1e-8)\n    return float((ac*bc).sum() / (da*db))\n\ndef safe_blend(w_net, w_tpl, fs):\n    # inputs same length & unscaled\n    w_net = np.asarray(w_net, dtype=np.float32)\n    w_tpl = np.asarray(w_tpl, dtype=np.float32)\n    w_net = np.nan_to_num(w_net, nan=0.0, posinf=0.0, neginf=0.0)\n    w_tpl = np.nan_to_num(w_tpl, nan=0.0, posinf=0.0, neginf=0.0)\n\n    # Lag search within ±0.2s\n    max_shift = int(fs * 0.2)\n    a = (w_net - w_net.mean()) / (w_net.std() + 1e-8)\n    b = (w_tpl - w_tpl.mean()) / (w_tpl.std() + 1e-8)\n    corr = scipy.signal.correlate(a, b, mode='full')\n    lags = scipy.signal.correlation_lags(len(a), len(b), mode='full')\n    mask = (lags >= -max_shift) & (lags <= max_shift)\n    lag = int(lags[mask][np.argmax(corr[mask])])\n\n    if lag > 0:\n        w_net_aligned = np.concatenate([np.full(lag, w_net[0], dtype=np.float32), w_net[:-lag]])\n    elif lag < 0:\n        lag2 = -lag\n        w_net_aligned = np.concatenate([w_net[lag2:], np.full(lag2, w_net[-1], dtype=np.float32)])\n    else:\n        w_net_aligned = w_net\n\n    # Sign check\n    r = ncc(w_net_aligned, w_tpl)\n    r_flip = ncc(-w_net_aligned, w_tpl)\n    if r_flip > r:\n        w_net_aligned = -w_net_aligned\n        r = r_flip\n\n    # Convert correlation to weight (template default)\n    if r < 0.15:\n        beta = 0.0\n    else:\n        beta = CFG.ENSEMBLE_BETA_GAIN * min(1.0, max(0.0, r))\n        beta = float(np.clip(beta, CFG.ENSEMBLE_BETA_CLIP[0], CFG.ENSEMBLE_BETA_CLIP[1]))\n\n    return beta, w_net_aligned\n\n# ---------------------------\n# Predict 12 leads for one image (ROI ensemble + tilt-aware TTA)\n# ---------------------------\n@torch.no_grad()\ndef predict_lead_waves_for_image(img01, roi_modes=(\"grid\",\"unet\"), roi_net_use=True):\n    # detect tilt\n    base_tilt = estimate_tilt_deg(img01)\n    ttas = [t+base_tilt for t in CFG.TTA_BASE_DEGS]\n    if CFG.TTA_EXTRA:\n        ttas += [base_tilt-0.8, base_tilt+0.8]\n    preds=[]\n    for deg in ttas:\n        im=rotate_image(img01,deg)\n        for mode in roi_modes:\n            use_net = (mode==\"unet\") and roi_net_use and (roi_net is not None)\n            crops = get_crops(im, roi_mode=(\"unet\" if use_net else mode), roi_net=(roi_net if use_net else None))\n            xs=[stack_crop_to_x3(c) for c in crops]\n            X=torch.tensor(np.stack(xs),dtype=torch.float32).to(device)\n            li=torch.arange(12, device=device)\n            if CFG.CHANNELS_LAST: X=X.to(memory_format=torch.channels_last)\n            with torch.cuda.amp.autocast(enabled=CFG.AMP):\n                out=model(X, li)  # (12,W)\n            preds.append(out.detach().cpu().numpy())\n    return np.mean(np.stack(preds,axis=0),axis=0).astype(np.float32)  # (12,W)\n\n# ---------------------------\n# Per-lead sign/gain calibration (small subset) — robust median\n# ---------------------------\n@torch.no_grad()\ndef compute_lead_calibration(sample_ids, max_images=32):\n    gains = {lead: [] for lead in LEADS}\n    for base_id in tqdm(list(sample_ids)[:max_images], desc=\"Calibrating (lead gains)\"):\n        p = pick_image_path(base_id, True)\n        if p is None: \n            continue\n        img01 = load_img_gray(p)\n        waves = predict_lead_waves_for_image(img01, roi_modes=CFG.ROI_MODE_TESTS, roi_net_use=True)  # (12, Wpred)\n        fs = int(train_csv.loc[train_csv.id==base_id, 'fs'].values[0])\n        M  = load_wave_csv_matrix(f\"{CFG.DATA_DIR}/train/{base_id}/{base_id}.csv\")\n        for lead in LEADS:\n            li = LEAD_TO_IDX[lead]\n            L  = int(fs*10) if lead==\"II\" else int(fs*2.5)\n            gt = resample_1d(M[li], L).astype(np.float32)\n            pr = resample_1d(waves[li], L).astype(np.float32)\n            gtc = gt - gt.mean()\n            prc = pr - pr.mean()\n            denom = float((prc**2).sum()) + 1e-6\n            a = float((prc*gtc).sum() / denom)\n            gains[lead].append(a)\n    lead_gain = {}\n    for lead, arr in gains.items():\n        if len(arr)==0:\n            lead_gain[lead] = 1.0\n        else:\n            med = float(np.median(np.clip(arr, -3.0, 3.0)))\n            lead_gain[lead] = float(np.clip(med, -2.5, 2.5))\n    return lead_gain\n\ncalib_ids = list(train_ids)[:32]\nLEAD_GAIN = compute_lead_calibration(calib_ids, max_images=32)\nprint(\"Per-lead gains:\", {k: round(v,3) for k,v in LEAD_GAIN.items()})\n\n# ---------------------------\n# Templates\n# ---------------------------\ndef build_lead_templates(template_len=500, max_samples=4000, cache_path=CFG.TEMPLATE_CACHE):\n    if os.path.exists(cache_path):\n        with open(cache_path, \"rb\") as f:\n            obj = pickle.load(f)\n        print(\"Loaded templates from cache.\")\n        return obj[\"templates\"], obj[\"counts\"]\n\n    counts = {l:0 for l in LEADS}\n    sums   = {l:np.zeros(template_len, dtype=np.float64) for l in LEADS}\n\n    ids = train_csv[\"id\"].values\n    np.random.shuffle(ids)\n    used = 0\n    for base_id in tqdm(ids, total=min(len(ids), max_samples), desc=\"Building lead templates\"):\n        if used >= max_samples: break\n        csv_path = f\"{CFG.DATA_DIR}/train/{base_id}/{base_id}.csv\"\n        if not os.path.exists(csv_path): continue\n        try:\n            df = pd.read_csv(csv_path)\n        except:\n            continue\n        used += 1\n        for lead in LEADS:\n            if lead not in df.columns: continue\n            s = pd.to_numeric(df[lead], errors=\"coerce\").dropna().values.astype(np.float32)\n            if len(s) < 50: continue\n            s = (s - s.mean()) / (s.std() + 1e-8)\n            s_resamp = np.interp(np.linspace(0,1,template_len), np.linspace(0,1,len(s)), s)\n            sums[lead] += s_resamp.astype(np.float64)\n            counts[lead] += 1\n\n    templates = {}\n    for lead in LEADS:\n        if counts[lead] > 0:\n            templates[lead] = (sums[lead] / counts[lead]).astype(np.float32)\n        else:\n            t = np.linspace(0,1,template_len, dtype=np.float32)\n            templates[lead] = np.sin(2*np.pi*t).astype(np.float32)\n\n    with open(cache_path, \"wb\") as f:\n        pickle.dump({\"templates\": templates, \"counts\": counts}, f)\n    print(\"Saved templates to cache.\")\n    return templates, counts\n\nLEAD_TEMPLATES, TEMPLATE_COUNTS = build_lead_templates(\n    template_len=CFG.TEMPLATE_LEN,\n    max_samples=CFG.TEMPLATE_MAX_SAMPLES,\n    cache_path=CFG.TEMPLATE_CACHE\n)\nprint(\"Template counts:\", TEMPLATE_COUNTS)\n\n# ---------------------------\n# Local SNR sanity (Lead II)\n# ---------------------------\nif CFG.DO_LOCAL_SNR:\n    try:\n        model.eval()\n        sample_ids = list(train_ids)[:4]\n        rows_sol, rows_sub = [], []\n        for base_id in tqdm(sample_ids, desc=\"Local SNR (Lead II only)\"):\n            p = pick_image_path(base_id, True)\n            if p is None: continue\n            img01 = load_img_gray(p)\n            waves = predict_lead_waves_for_image(img01, roi_modes=CFG.ROI_MODE_TESTS, roi_net_use=True)  # (12, Wpred)\n\n            fs = int(train_csv[train_csv.id == base_id]['fs'].values[0])\n            M  = load_wave_csv_matrix(f\"{CFG.DATA_DIR}/train/{base_id}/{base_id}.csv\")\n\n            lead = \"II\"; li = LEAD_TO_IDX[lead]; L = int(fs * 10)\n            gt = resample_1d(M[li], L)\n            pr = waves[li].astype(np.float32)\n            pr = LEAD_GAIN.get(lead, 1.0) * pr\n            pr = resample_1d(pr, L)\n\n            # lowpass & scale for quick local check\n            pr = lowpass(pr, fs, order=2, cutoff_hz=15.0)\n            pr = scale_to_range(pr, CFG.OUT_MIN, CFG.OUT_MAX)\n\n            for j in range(L):\n                rid = f\"{base_id}_{j}_{lead}\"\n                rows_sol.append((rid, fs, gt[j]))\n                rows_sub.append((rid, pr[j]))\n        sol = pd.DataFrame(rows_sol, columns=[\"id\",\"fs\",\"value\"])\n        sub = pd.DataFrame(rows_sub, columns=[\"id\",\"value\"])\n        local_snr = competition_score(sol, sub, \"id\")\n        print(\"Local SNR (Lead II only, dB):\", round(local_snr, 4))\n    except Exception as e:\n        print(\"Local SNR eval skipped:\", repr(e))\n\n# ---------------------------\n# Test Inference + Safe Template Blend + Submission\n# ---------------------------\nmodel.eval()\nrows=[]\nfor base_id in tqdm(test_csv[\"id\"].unique(), desc=\"Test inference\"):\n    p=f\"{CFG.DATA_DIR}/test/{base_id}.png\"\n    if not os.path.exists(p): continue\n    img01=load_img_gray(p)\n\n    # network prediction (12 x Wpred)\n    waves=predict_lead_waves_for_image(img01, roi_modes=CFG.ROI_MODE_TESTS, roi_net_use=True)\n\n    subset=test_csv[test_csv[\"id\"]==base_id]\n    for _,r in subset.iterrows():\n        lead=r[\"lead\"]; fs=int(r[\"fs\"]); nrows=int(r[\"number_of_rows\"])\n        li=LEAD_TO_IDX[lead]\n\n        # --- network path (unscaled first) ---\n        w = waves[li].astype(np.float32)\n        w = LEAD_GAIN.get(lead, 1.0) * w\n        w = np.nan_to_num(w, nan=0.0, posinf=0.0, neginf=0.0).astype(np.float32)\n        w = resample_1d(w, nrows)  # resample to per-lead length\n\n        # --- template path (unscaled) ---\n        tpl = LEAD_TEMPLATES.get(lead, LEAD_TEMPLATES['II']).copy().astype(np.float32)\n        tpl = np.interp(np.linspace(0,1,nrows), np.linspace(0,1,len(tpl)), tpl).astype(np.float32)\n\n        # --- safe blend (align + sign safety; template default) ---\n        beta, w_aligned = safe_blend(w, tpl, fs)\n\n        # --- Low-pass & scale after alignment ---\n        w_final   = scale_to_range(lowpass(w_aligned, fs, order=2, cutoff_hz=15.0), CFG.OUT_MIN, CFG.OUT_MAX)\n        tpl_final = scale_to_range(lowpass(tpl,       fs, order=2, cutoff_hz=15.0), CFG.OUT_MIN, CFG.OUT_MAX)\n\n        final = beta * w_final + (1.0 - beta) * tpl_final\n        final = np.clip(final, CFG.OUT_MIN, CFG.OUT_MAX).astype(np.float32)\n\n        for j in range(nrows):\n            rows.append((f\"{base_id}_{j}_{lead}\", float(final[j])))\n\n# --- Write submissions to Kaggle working output and verify ---\nout_dir = \"/kaggle/working\"\nos.makedirs(out_dir, exist_ok=True)\n\nsubmission = pd.DataFrame(rows, columns=[\"id\",\"value\"]).astype({\"id\":\"string\",\"value\":\"float32\"})\nsubmission_parquet = os.path.join(out_dir, \"submission.parquet\")\nsubmission_csv     = os.path.join(out_dir, \"submission.csv\")\nsubmission.to_parquet(submission_parquet, index=False)\nsubmission.to_csv(submission_csv, index=False)\n\nprint(\"✅ Saved:\", submission_parquet, submission_csv)\nprint(\"CWD:\", os.getcwd())\nprint(\"Working dir listing:\", os.listdir(out_dir)[:20])\nprint(submission.head())\nprint(\"✅ Saved submission.parquet; rows:\", len(submission))\n\n# ===========================\n# Submission validation + OFFLINE LB estimate (official metric)\n# ===========================\nprint(\"\\n=== Submission validation ===\")\nsub = submission.copy()\n# Basic checks\nn_ids_expected = len(test_csv)\nn_rows_expected = int(test_csv['number_of_rows'].sum())\nprint(\"Unique test ECGs:\", n_ids_expected)\nprint(\"Expected rows:\", n_rows_expected, \" | Submission rows:\", len(sub))\n\nassert len(sub) == n_rows_expected, \"Row count mismatch.\"\nassert sub['id'].is_unique, \"IDs are not unique.\"\nassert sub['value'].isna().sum() == 0, \"NaNs in submission.\"\nassert np.isfinite(sub['value']).all(), \"Non-finite values in submission.\"\nassert sub['value'].between(CFG.OUT_MIN - 1e-6, CFG.OUT_MAX + 1e-6).all(), \"Values out of range.\"\nprint(\"✅ Submission looks valid.\")\n\n# ---------------------------\n# OFFLINE LB estimate on held-out train (uses the official scoring functions)\n# ---------------------------\nfrom typing import Tuple\nimport scipy.optimize, scipy.signal\n\nLEADS = ['I','II','III','aVR','aVL','aVF','V1','V2','V3','V4','V5','V6']\nMAX_TIME_SHIFT = 0.2\nPERFECT_SCORE = 384\n\nclass ParticipantVisibleError(Exception): pass\n\ndef compute_power(label: np.ndarray, prediction: np.ndarray) -> Tuple[float, float]:\n    if label.ndim != 1 or prediction.ndim != 1:\n        raise ParticipantVisibleError('Inputs must be 1-dimensional arrays.')\n    finite_mask = np.isfinite(prediction)\n    if not np.any(finite_mask):\n        raise ParticipantVisibleError(\"The 'prediction' array contains no finite values (all NaN or inf).\")\n    prediction = prediction.copy()\n    prediction[~np.isfinite(prediction)] = 0\n    noise = label - prediction\n    p_signal = np.sum(label**2)\n    p_noise = np.sum(noise**2)\n    return p_signal, p_noise\n\ndef compute_snr(signal: float, noise: float) -> float:\n    if noise == 0:\n        return PERFECT_SCORE\n    if signal == 0:\n        return 0\n    return min((signal / noise), PERFECT_SCORE)\n\ndef align_signals(label: np.ndarray, pred: np.ndarray, max_shift: float = float('inf')) -> np.ndarray:\n    if np.any(~np.isfinite(label)):\n        raise ParticipantVisibleError('values in label should all be finite')\n    if np.sum(np.isfinite(pred)) == 0:\n        raise ParticipantVisibleError('prediction can not all be infinite')\n    label_arr = np.asarray(label, dtype=np.float64)\n    pred_arr = np.asarray(pred, dtype=np.float64)\n    label_mean = np.mean(label_arr)\n    pred_mean = np.mean(pred_arr)\n    label_arr_centered = label_arr - label_mean\n    pred_arr_centered = pred_arr - pred_mean\n    correlation = scipy.signal.correlate(label_arr_centered, pred_arr_centered, mode='full')\n    n_label = np.size(label_arr)\n    n_pred = np.size(pred_arr)\n    lags = scipy.signal.correlation_lags(n_label, n_pred, mode='full')\n    valid_lags_mask = (lags >= -max_shift) & (lags <= max_shift)\n    max_correlation = np.nanmax(correlation[valid_lags_mask])\n    all_max_indices = np.flatnonzero(correlation == max_correlation)\n    best_idx = min(all_max_indices, key=lambda i: abs(lags[i]))\n    time_shift = lags[best_idx]\n    start_padding_len = max(time_shift, 0)\n    pred_slice_start = max(-time_shift, 0)\n    pred_slice_end = min(n_label - time_shift, n_pred)\n    end_padding_len = max(n_label - n_pred - time_shift, 0)\n    aligned_pred = np.concatenate((np.full(start_padding_len, np.nan), pred_arr[pred_slice_start:pred_slice_end], np.full(end_padding_len, np.nan)))\n    def objective_func(v_shift):\n        return np.nansum((label_arr - (aligned_pred - v_shift)) ** 2)\n    if np.any(np.isfinite(label_arr) & np.isfinite(aligned_pred)):\n        results = scipy.optimize.minimize_scalar(objective_func, method='Brent')\n        vertical_shift = results.x\n        aligned_pred -= vertical_shift\n    return aligned_pred\n\ndef _calculate_image_score(group: pd.DataFrame) -> float:\n    unique_fs_values = group['fs'].unique()\n    if len(unique_fs_values) != 1:\n        raise ParticipantVisibleError('Sampling frequency should be consistent across each ecg')\n    sampling_frequency = unique_fs_values[0]\n    if sampling_frequency != int(len(group[group['lead']=='II']) / 10):\n        raise ParticipantVisibleError('The sequence_length should be sampling frequency * 10s')\n    sum_signal = 0.0\n    sum_noise = 0.0\n    for lead in LEADS:\n        sub = group[group['lead'] == lead]\n        label = sub['value_true'].values\n        pred  = sub['value_pred'].values\n        aligned_pred = align_signals(label, pred, int(sampling_frequency * MAX_TIME_SHIFT))\n        p_signal, p_noise = compute_power(label, aligned_pred)\n        sum_signal += p_signal\n        sum_noise  += p_noise\n    return compute_snr(sum_signal, sum_noise)\n\ndef score(solution: pd.DataFrame, submission: pd.DataFrame, row_id_column_name: str) -> float:\n    for df in [solution, submission]:\n        if row_id_column_name not in df.columns:\n            raise ParticipantVisibleError(f\"'{row_id_column_name}' column not found in DataFrame.\")\n        if df['value'].isna().any():\n            raise ParticipantVisibleError('NaN exists in solution/submission')\n        if not np.isfinite(df['value']).all():\n            raise ParticipantVisibleError('Infinity exists in solution/submission')\n    submission = submission[['id', 'value']]\n    merged_df = pd.merge(solution, submission, on=row_id_column_name, suffixes=('_true', '_pred'))\n    merged_df['image_id'] = merged_df[row_id_column_name].str.split('_').str[0]\n    merged_df['row_id']   = merged_df[row_id_column_name].str.split('_').str[1].astype('int64')\n    merged_df['lead']     = merged_df[row_id_column_name].str.split('_').str[2]\n    merged_df.sort_values(by=['image_id','row_id','lead'], inplace=True)\n    image_scores = merged_df.groupby('image_id').apply(_calculate_image_score, include_groups=False)\n    return max(float(10 * np.log10(image_scores.mean())), -PERFECT_SCORE)\n\n# ---- Build an offline validation set from held-out train (small sample) ----\nnp.random.seed(CFG.SEED)\nval_eval_ids = list(val_ids)[:12]   # evaluate ~12 images for speed\n\nrows_sol, rows_sub = [], []\nfor base_id in tqdm(val_eval_ids, desc=\"Offline LB (val subset)\"):\n    ip = pick_image_path(base_id, True)\n    if ip is None: \n        continue\n    img01 = load_img_gray(ip)\n    waves = predict_lead_waves_for_image(img01, roi_modes=CFG.ROI_MODE_TESTS, roi_net_use=False)  # model inference only; blend below\n\n    fs = int(train_csv.loc[train_csv.id==base_id, 'fs'].values[0])\n    M  = load_wave_csv_matrix(f\"{CFG.DATA_DIR}/train/{base_id}/{base_id}.csv\")\n\n    for lead in LEADS:\n        li = LEAD_TO_IDX[lead]\n        nrows = int(fs * (10 if lead == \"II\" else 2.5))\n\n        # network (unscaled)\n        w = waves[li].astype(np.float32)\n        w = LEAD_GAIN.get(lead, 1.0) * w\n        w = resample_1d(np.nan_to_num(w, nan=0.0, posinf=0.0, neginf=0.0), nrows)\n\n        # template (unscaled)\n        tpl = LEAD_TEMPLATES.get(lead, LEAD_TEMPLATES['II']).copy().astype(np.float32)\n        tpl = np.interp(np.linspace(0,1,nrows), np.linspace(0,1,len(tpl)), tpl).astype(np.float32)\n\n        # safe blend\n        beta, w_aligned = safe_blend(w, tpl, fs)\n        w_final   = scale_to_range(lowpass(w_aligned, fs, order=2, cutoff_hz=15.0), CFG.OUT_MIN, CFG.OUT_MAX)\n        tpl_final = scale_to_range(lowpass(tpl,       fs, order=2, cutoff_hz=15.0), CFG.OUT_MIN, CFG.OUT_MAX)\n        pred = np.clip(beta * w_final + (1.0 - beta) * tpl_final, CFG.OUT_MIN, CFG.OUT_MAX).astype(np.float32)\n\n        # GT\n        gt = resample_1d(M[li], nrows).astype(np.float32)\n\n        for j in range(nrows):\n            rows_sol.append((f\"{base_id}_{j}_{lead}\", fs, float(gt[j])))\n            rows_sub.append((f\"{base_id}_{j}_{lead}\", float(pred[j])))\n\nsol_df = pd.DataFrame(rows_sol, columns=[\"id\",\"fs\",\"value\"])\nsub_df = pd.DataFrame(rows_sub, columns=[\"id\",\"value\"])\ntry:\n    offline_snr_db = score(sol_df, sub_df, \"id\")\n    print(f\"📈 Offline LB estimate (val subset, dB): {offline_snr_db:.4f}\")\nexcept Exception as e:\n    print(\"Offline scoring skipped:\", repr(e))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null}]}