{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","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":"gpu","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# # ECG Image Digitization — **ML‑Hybrid Baseline** (Fast, Deterministic, Competition‑Ready)\n# \n# **Upgrades vs pure signal baseline**\n# - Per‑lead morphology via **PCA templates** (learned from beats)\n# - **BPM regressor** (RandomForest) using rhythm/spectral features\n# - Same deterministic synthesis (tiling + gentle filtering)\n# - Optional **Einthoven blending** for limb leads\n# - Optional tiny **CNN denoiser** (disabled by default to keep runtime short)\n# \n# **I/O** follows the Kaggle PhysioNet ECG Image Digitization dataset structure.\n# \n# If this helps, an upvote ⭐ on Kaggle helps others find it! 🙌","metadata":{"_uuid":"46b7ff09-74c6-4bb0-92e3-967a98fbe808","_cell_guid":"23ba28c6-69ab-41a4-8961-d443255b9f8f","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false}}},{"cell_type":"code","source":"# --- Imports, Reproducibility, and Global Config ---\nimport os, random, math, warnings\nimport numpy as np\nimport pandas as pd\nfrom tqdm.auto import tqdm\nfrom scipy.signal import butter, sosfiltfilt, find_peaks\nfrom sklearn.decomposition import PCA\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import r2_score, mean_absolute_error\nwarnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n\n# Optional: lightweight CNN denoiser (disabled by default)\ntry:\n    import torch\n    import torch.nn as nn\n    import torch.optim as optim\n    TORCH_OK = True\nexcept Exception:\n    TORCH_OK = False\n\n# Reproducibility\nos.environ[\"PYTHONHASHSEED\"] = \"0\"\nrandom.seed(0)\nnp.random.seed(0)\n\n# --- Paths (Kaggle dataset structure) ---\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'\n\n# --- Config ---\nLEADS = ['I','II','III','aVR','aVL','aVF','V1','V2','V3','V4','V5','V6']\nR_PRE_S, R_POST_S = 0.20, 0.40\nBEAT_LEN = 360\nBP_LO_HZ, BP_HI_HZ, BP_ORDER = 5.0, 25.0, 2\nBPM_CANDIDATES = [45, 55, 65, 75, 85, 95, 105]\nMIN_VAL, MAX_VAL = 0.0, 0.09\nEINTHOVEN_BLEND_W = 0.6\nENSEMBLE_W = np.array([0.5, 0.3, 0.2], dtype=np.float32)\nSAFE_MODE = True\n\n# Feature/Model knobs\nPCA_COMPONENTS = 6             # per‑lead PCA size\nRF_TREES = 200                 # BPM regressor trees\nENABLE_CNN_DENOISER = True    # keep False for speed; set True to train tiny CNN\nCNN_EPOCHS = 5                 # very small, to keep runtime reasonable","metadata":{"_uuid":"f4c287a8-38f0-4e91-b495-1551fbd2a7f9","_cell_guid":"2967d776-5bf9-46cc-8fd6-f51227e0516b","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:44:33.612891Z","iopub.execute_input":"2025-10-27T16:44:33.613726Z","iopub.status.idle":"2025-10-27T16:44:33.623502Z","shell.execute_reply.started":"2025-10-27T16:44:33.613695Z","shell.execute_reply":"2025-10-27T16:44:33.622480Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Utilities ---\n\ndef zscore(x):\n    x = np.asarray(x, np.float32)\n    mu, s = np.mean(x), np.std(x)\n    if not np.isfinite(s) or s < 1e-8: s = 1e-8\n    return (x - mu) / s\n\n\ndef bandpass(x, fs, lo=BP_LO_HZ, hi=BP_HI_HZ, order=BP_ORDER):\n    if len(x) < 32:\n        return x\n    nyq = 0.5 * fs\n    lo_n, hi_n = max(lo/nyq, 1e-3), min(hi/nyq, 0.99)\n    if hi_n <= lo_n + 1e-3:\n        return x\n    sos = butter(order, [lo_n, hi_n], btype='band', output='sos')\n    try:\n        return sosfiltfilt(sos, x).astype(np.float32)\n    except Exception:\n        return x.astype(np.float32)\n\n\ndef apply_lowpass(x, fs, cutoff=15.0, order=2):\n    if len(x) < 32:\n        return x\n    nyq = 0.5 * fs\n    wn = min(cutoff/nyq, 0.99)\n    if wn <= 1e-3:\n        return x.astype(np.float32)\n    sos = butter(order, wn, btype='low', output='sos')\n    try:\n        return sosfiltfilt(sos, x).astype(np.float32)\n    except Exception:\n        return x.astype(np.float32)\n\n\ndef _nextpow2(n):\n    return 1 << (int(math.ceil(math.log2(max(1, n)))))\n\n\ndef autocorr_peak_score(y, fs, min_rr_s=0.35, max_rr_s=1.8):\n    y = zscore(y)\n    n = len(y)\n    if n < 16: return 0.0\n    m = _nextpow2(2*n - 1)\n    Y = np.fft.rfft(y, n=m)\n    ac = np.fft.irfft(Y * np.conj(Y), n=m)[:n]\n    lo, hi = int(min_rr_s * fs), int(max_rr_s * fs)\n    if hi <= lo: return 0.0\n    seg = ac[lo:hi]\n    peak = np.max(seg) if seg.size else 0\n    return float(np.clip(peak / (ac[0] + 1e-8), 0.0, 1.0))\n\n\ndef soft_minmax_scale(x, lo=MIN_VAL, hi=MAX_VAL):\n    x = np.asarray(x, np.float32)\n    if x.size == 0:\n        return np.full(1, (lo + hi) / 2, np.float32)\n    mn, mx = np.nanmin(x), np.nanmax(x)\n    if not np.isfinite(mn) or not np.isfinite(mx) or mx <= mn:\n        return np.full_like(x, (lo + hi) / 2, np.float32)\n    y = (x - mn) / (mx - mn)\n    return np.clip(lo + y * (hi - lo), lo, hi).astype(np.float32)\n\n\ndef scale_to_lead_range(x, lead_stat=None, lo=MIN_VAL, hi=MAX_VAL):\n    if lead_stat is not None:\n        mn, mx = lead_stat.get('min', -0.5), lead_stat.get('max', 0.5)\n        x = np.clip(x, mn, mx)\n    return soft_minmax_scale(x, lo, hi)\n\n\ndef resample_to_length(x, n):\n    return np.interp(\n        np.linspace(0, 1, n, dtype=np.float32),\n        np.linspace(0, 1, len(x), dtype=np.float32),\n        x\n    ).astype(np.float32)\n\n\ndef derive_limb_leads_from_I_II(yI, yII):\n    III  = yII - yI\n    aVR  = -(yI + yII) / 2.0\n    aVL  = yI - 0.5 * yII\n    aVF  = yII - 0.5 * yI\n    return {'III': III, 'aVR': aVR, 'aVL': aVL, 'aVF': aVF}\n\n\ndef soft_blend(a, b, w):\n    return (1 - w) * a + w","metadata":{"_uuid":"6a5663b0-47ca-4b5f-bd39-a29b15def127","_cell_guid":"8c6b25b1-1fdc-4a0d-9c6c-02048d025a05","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:44:33.624863Z","iopub.execute_input":"2025-10-27T16:44:33.625193Z","iopub.status.idle":"2025-10-27T16:44:33.655952Z","shell.execute_reply.started":"2025-10-27T16:44:33.625171Z","shell.execute_reply":"2025-10-27T16:44:33.654786Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step A: Scan training set to collect stats, beats, and BPM supervision ---\n\ndef build_per_lead_assets(train_csv, train_dir, leads=LEADS):\n    meta = pd.read_csv(train_csv)\n    lead_vals = {ld: [] for ld in leads}\n    lead_beats = {ld: [] for ld in leads}\n    lead_bpm_samples = {ld: [] for ld in leads}\n\n    # For BPM regressor supervision (features/targets per record/lead)\n    bpm_feats, bpm_targets = [], []\n\n    for row in tqdm(meta.itertuples(index=False), total=len(meta), desc=\"Scan train\"):\n        rid = str(row.id)\n        fs  = int(row.fs)\n        csvp = os.path.join(train_dir, rid, f\"{rid}.csv\")\n        if not os.path.exists(csvp):\n            continue\n        try:\n            df = pd.read_csv(csvp)\n        except Exception:\n            continue\n\n        for ld in leads:\n            if ld not in df.columns:\n                continue\n            y = df[ld].dropna().to_numpy(np.float32)\n            if len(y) < 200:\n                continue\n\n            # stats pool (raw)\n            lead_vals[ld].append(y)\n\n            # R-peak detection on band-passed signal\n            y_bp = bandpass(zscore(y), fs)\n            iqr = np.subtract(*np.percentile(y_bp, [75, 25]))\n            scale = iqr if np.isfinite(iqr) and iqr > 0 else np.std(y_bp)\n            prominence = max(0.25 * scale, 0.08)\n            distance = int(max(0.28 * fs, 1))\n            pks, _ = find_peaks(y_bp, distance=distance, prominence=prominence)\n            if len(pks) < 2:\n                continue\n\n            rr = np.diff(pks) / float(fs)\n            rr = rr[(rr > 0.3) & (rr < 2.0)]\n            if rr.size >= 1:\n                bpm = float(np.clip(60.0 / np.median(rr), 40.0, 160.0))\n                lead_bpm_samples[ld].append(bpm)\n\n                # --- BPM features per record/lead ---\n                # Autocorr and spectral/simple stats\n                y_norm = zscore(y)\n                ac_sc = autocorr_peak_score(y_norm, fs)\n                dur_s = len(y) / fs\n                var = float(np.var(y_norm))\n                mad = float(np.mean(np.abs(np.diff(y_norm))))\n                # spectral proxy: power in 0.5–15 Hz band vs 15–40 Hz\n                def band_power(x, fs, lo, hi):\n                    X = np.fft.rfft(x)\n                    freqs = np.fft.rfftfreq(len(x), d=1.0/fs)\n                    m = (freqs>=lo) & (freqs<=hi)\n                    return float(np.sum(np.abs(X[m])**2) + 1e-8)\n                p_lo = band_power(y_norm, fs, 0.5, 15.0)\n                p_hi = band_power(y_norm, fs, 15.0, 40.0)\n                ratio = p_lo / (p_hi + 1e-8)\n                bpm_feats.append([fs, dur_s, ac_sc, var, mad, ratio])\n                bpm_targets.append(bpm)\n\n            # beat windows around R\n            n_pre  = int(round(R_PRE_S * fs))\n            n_post = int(round(R_POST_S * fs))\n            for pk in pks:\n                a, b = pk - n_pre, pk + n_post\n                if a < 0 or b >= len(y):\n                    continue\n                seg = y[a:b+1].astype(np.float32)\n                seg_rs = resample_to_length(seg, BEAT_LEN)\n                lead_beats[ld].append(seg_rs)\n\n    # Lead stats & trimmed median beat\n    lead_stats, lead_templates_raw = {}, {}\n    for ld in leads:\n        if len(lead_vals[ld]) == 0:\n            lead_stats[ld] = {'mean':0.0,'std':0.1,'min':-0.5,'max':0.5}\n        else:\n            vals = np.concatenate(lead_vals[ld]).astype(np.float32)\n            lead_stats[ld] = {\n                'mean': float(np.mean(vals)),\n                'std':  float(np.std(vals) if vals.size>1 else 0.1),\n                'min':  float(np.min(vals)),\n                'max':  float(np.max(vals))\n            }\n        if len(lead_beats[ld])>0:\n            arr = np.vstack(lead_beats[ld]).astype(np.float32)\n            if arr.shape[0] > 5:\n                provisional = np.median(arr, axis=0)\n                d = np.linalg.norm(arr - provisional, axis=1)\n                keep = d < np.percentile(d, 90)\n                arr = arr[keep]\n            lead_templates_raw[ld] = np.median(arr, axis=0).astype(np.float32)\n        else:\n            t = np.linspace(0,1,BEAT_LEN,dtype=np.float32)\n            lead_templates_raw[ld] = np.sin(2*np.pi*t).astype(np.float32)\n\n    bpm_feats = np.array(bpm_feats, dtype=np.float32)\n    bpm_targets = np.array(bpm_targets, dtype=np.float32)\n    return lead_stats, lead_templates_raw, lead_beats, lead_bpm_samples, bpm_feats, bpm_targets","metadata":{"_uuid":"0020dc17-da37-412b-ae68-ba0945890a8c","_cell_guid":"b3fe6f06-6eb6-48cf-a232-b5e20fcf3834","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:44:33.657418Z","iopub.execute_input":"2025-10-27T16:44:33.658168Z","iopub.status.idle":"2025-10-27T16:44:33.686057Z","shell.execute_reply.started":"2025-10-27T16:44:33.658143Z","shell.execute_reply":"2025-10-27T16:44:33.685057Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"[1/5] Building assets from train…\")\nlead_stats, lead_templates_raw, lead_beats, lead_bpms, bpm_feats, bpm_targets = build_per_lead_assets(TRAIN_CSV, TRAIN_DIR, LEADS)\nprint(\"  BPM samples:\", {k: len(v) for k,v in lead_bpms.items()})\nprint(\"  BPM reg feats:\", bpm_feats.shape)","metadata":{"_uuid":"9f0b99ed-1c92-434a-a36b-85f3a20c080d","_cell_guid":"f0b86199-78f2-4e02-9ae2-4cda9fd99094","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:44:33.882167Z","iopub.execute_input":"2025-10-27T16:44:33.883496Z","iopub.status.idle":"2025-10-27T16:45:19.349260Z","shell.execute_reply.started":"2025-10-27T16:44:33.883448Z","shell.execute_reply":"2025-10-27T16:45:19.348109Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step B: Per‑lead PCA templates ---\nlead_pca = {}\nlead_template = {}\nfor ld in LEADS:\n    beats = lead_beats[ld]\n    if len(beats) >= max(20, PCA_COMPONENTS+5):\n        X = np.vstack(beats).astype(np.float32)\n        pca = PCA(n_components=PCA_COMPONENTS, random_state=0).fit(X)\n        lead_pca[ld] = pca\n        mean_code = np.zeros(PCA_COMPONENTS, dtype=np.float32)\n        tpl = pca.inverse_transform(mean_code)\n        lead_template[ld] = zscore(tpl)\n    else:\n        lead_pca[ld] = None\n        lead_template[ld] = zscore(lead_templates_raw[ld])\nprint(\"[2/5] PCA templates built.\")","metadata":{"_uuid":"690609e8-8d55-403f-9d43-2d226a00d1cf","_cell_guid":"f19ee874-3950-4212-a2af-cda6b0798876","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:45:19.360351Z","iopub.execute_input":"2025-10-27T16:45:19.360718Z","iopub.status.idle":"2025-10-27T16:45:20.688990Z","shell.execute_reply.started":"2025-10-27T16:45:19.360691Z","shell.execute_reply":"2025-10-27T16:45:20.685915Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step C: BPM Regressor (RandomForest) ---\nif len(bpm_targets) >= 50:\n    rf = RandomForestRegressor(n_estimators=RF_TREES, random_state=0, n_jobs=-1)\n    rf.fit(bpm_feats, bpm_targets)\n    pred = rf.predict(bpm_feats)\n    print(\"[3/5] BPM RF: R2=%.3f, MAE=%.2f bpm\" % (r2_score(bpm_targets, pred), mean_absolute_error(bpm_targets, pred)))\nelse:\n    rf = None\n    print(\"[3/5] Not enough BPM supervision, skipping RF.\")\n\n# Simple amplitude scaler: per‑lead linear model from (std,iqr) → (clip range)\n# It refines the per‑lead min/max from training to be a bit more adaptive at test.\nlead_scalers = {}\nfor ld in LEADS:\n    # gather training stats per record for this lead (approximate using beats)\n    if len(lead_beats[ld]) >= 20:\n        Xs, ys_min, ys_max = [], [], []\n        # sample up to 200 beats for statistics\n        arr = np.vstack(lead_beats[ld])\n        step = max(1, arr.shape[0] // 200)\n        for seg in arr[::step]:\n            s = np.std(seg); i = np.subtract(*np.percentile(seg, [75,25]))\n            Xs.append([s, i])\n            ys_min.append(np.min(seg))\n            ys_max.append(np.max(seg))\n        Xs = np.array(Xs, np.float32)\n        ys_min = np.array(ys_min, np.float32)\n        ys_max = np.array(ys_max, np.float32)\n        if len(Xs) > 10:\n            lr_min = LinearRegression().fit(Xs, ys_min)\n            lr_max = LinearRegression().fit(Xs, ys_max)\n            lead_scalers[ld] = (lr_min, lr_max)\n        else:\n            lead_scalers[ld] = None\n    else:\n        lead_scalers[ld] = None\nprint(\"[3.5/5] Per‑lead amplitude scalers ready.\")","metadata":{"_uuid":"085624e8-4581-48ac-95a0-d372c949e3eb","_cell_guid":"9c568c2f-51fd-423a-94dd-9d07efcb9f71","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:45:20.690210Z","iopub.execute_input":"2025-10-27T16:45:20.690522Z","iopub.status.idle":"2025-10-27T16:45:25.741030Z","shell.execute_reply.started":"2025-10-27T16:45:20.690502Z","shell.execute_reply":"2025-10-27T16:45:25.739791Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step D: Optional Tiny CNN denoiser (disabled by default) ---\nclass SmallDenoiser(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.net = nn.Sequential(\n            nn.Conv1d(1, 8, 9, padding=4), nn.ReLU(),\n            nn.Conv1d(8, 8, 9, padding=4), nn.ReLU(),\n            nn.Conv1d(8, 1, 9, padding=4)\n        )\n    def forward(self, x):\n        return self.net(x)\n\ndenoiser = None\nif ENABLE_CNN_DENOISER and TORCH_OK:\n    print(\"[4/5] Training tiny CNN denoiser (quick)…\")\n    # Build a tiny train set: synthesized (noisy) → true beat median\n    X_train, Y_train = [], []\n    for ld in LEADS:\n        beats = lead_beats[ld]\n        if len(beats) < 40: continue\n        tpl = lead_template[ld]\n        # corrupt with mild noise to simulate synthesis artifacts\n        for b in beats[:200]:\n            x = tpl + 0.15*zscore(np.random.randn(len(tpl)).astype(np.float32))\n            X_train.append(x.astype(np.float32))\n            Y_train.append(zscore(b).astype(np.float32))\n    if len(X_train) > 1000:\n        X = torch.tensor(np.stack(X_train)[:,None,:])\n        Y = torch.tensor(np.stack(Y_train)[:,None,:])\n        denoiser = SmallDenoiser()\n        opt = optim.Adam(denoiser.parameters(), lr=1e-3)\n        loss_fn = nn.MSELoss()\n        denoiser.train()\n        for ep in range(CNN_EPOCHS):\n            perm = torch.randperm(X.size(0))\n            for i in range(0, X.size(0), 64):\n                idx = perm[i:i+64]\n                xb, yb = X[idx], Y[idx]\n                out = denoiser(xb)\n                loss = loss_fn(out, yb)\n                opt.zero_grad(); loss.backward(); opt.step()\n            print(f\"  epoch {ep+1}/{CNN_EPOCHS} — loss {loss.item():.5f}\")\n        denoiser.eval()\n    else:\n        print(\"  Not enough data for CNN; skipping.\")\nelse:\n    print(\"[4/5] CNN denoiser disabled (set ENABLE_CNN_DENOISER=True to try).\")","metadata":{"_uuid":"c21aa85e-75a5-45ef-bbd1-f1ed6b3a5fed","_cell_guid":"bd30cf92-6dc9-49f9-afec-842678dc9778","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:45:25.742488Z","iopub.execute_input":"2025-10-27T16:45:25.743384Z","iopub.status.idle":"2025-10-27T16:45:25.755668Z","shell.execute_reply.started":"2025-10-27T16:45:25.743357Z","shell.execute_reply":"2025-10-27T16:45:25.754539Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step E: Synthesis helpers (with ML upgrades) ---\n\ndef tile_template(template_beat, fs, n_out, bpm, amp=1.0):\n    beat_samples = max(4, int(round((60.0 / max(bpm, 1e-6)) * fs)))\n    one = resample_to_length(template_beat, beat_samples)\n    reps = int(np.ceil(n_out / len(one)))\n    y = np.tile(one, reps)[:n_out]\n    return zscore(y) * float(amp)\n\n\ndef choose_bpm(template_beat, fs, n_out):\n    # ML path if RF exists\n    if rf is not None:\n        # Build the same features used in training but on the raw template signal of length n_out\n        # Use autocorr best tiling as a candidate to compute features\n        # First, a quick autocorr sweep to generate a provisional y\n        best_bpm, best_sc = None, -1\n        for bpm in BPM_CANDIDATES:\n            y = tile_template(template_beat, fs, n_out, bpm, amp=1.0)\n            sc = autocorr_peak_score(y, fs)\n            if sc > best_sc:\n                best_sc, best_bpm = sc, bpm\n        # Feature vector approximated from provisional realization\n        y0 = tile_template(template_beat, fs, n_out, best_bpm, amp=1.0)\n        y0 = zscore(y0)\n        ac_sc = autocorr_peak_score(y0, fs)\n        dur_s = n_out / fs\n        var = float(np.var(y0))\n        mad = float(np.mean(np.abs(np.diff(y0))))\n        def band_power(x, fs, lo, hi):\n            X = np.fft.rfft(x)\n            freqs = np.fft.rfftfreq(len(x), d=1.0/fs)\n            m = (freqs>=lo) & (freqs<=hi)\n            return float(np.sum(np.abs(X[m])**2) + 1e-8)\n        p_lo = band_power(y0, fs, 0.5, 15.0)\n        p_hi = band_power(y0, fs, 15.0, 40.0)\n        ratio = p_lo / (p_hi + 1e-8)\n        feat = np.array([[fs, dur_s, ac_sc, var, mad, ratio]], dtype=np.float32)\n        bpm_pred = float(np.clip(rf.predict(feat)[0], 40.0, 160.0))\n        # also try nearby candidates and pick best by autocorr\n        candidates = sorted(set(BPM_CANDIDATES + [int(round(bpm_pred))]))\n        best_bpm, best_y, best_sc = None, None, -1\n        for bpm in candidates:\n            y = tile_template(template_beat, fs, n_out, bpm, amp=1.0)\n            sc = autocorr_peak_score(y, fs)\n            if sc > best_sc:\n                best_sc, best_bpm, best_y = sc, bpm, y\n        return best_bpm, best_y\n    else:\n        # fallback: pure autocorr sweep\n        best_bpm, best_score, best_y = None, -1.0, None\n        for bpm in BPM_CANDIDATES:\n            y = tile_template(template_beat, fs, n_out, bpm, amp=1.0)\n            sc = autocorr_peak_score(y, fs)\n            if sc > best_score:\n                best_bpm, best_score, best_y = bpm, sc, y\n        return best_bpm, best_y","metadata":{"_uuid":"0fac8856-799c-472b-b066-297c9a4e90b5","_cell_guid":"815f8555-9fe3-4f14-af38-bb398b1e46c0","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:45:25.757685Z","iopub.execute_input":"2025-10-27T16:45:25.758072Z","iopub.status.idle":"2025-10-27T16:45:25.778481Z","shell.execute_reply.started":"2025-10-27T16:45:25.758049Z","shell.execute_reply":"2025-10-27T16:45:25.777380Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Predict Test with Micro‑Ensemble + Einthoven blending + optional denoiser ---\nprint(\"[5/5] Predicting test…\")\n\ntest = pd.read_csv(TEST_CSV)\nrecords = {}\nfor r in test.itertuples(index=False):\n    records.setdefault(int(r.id), []).append(r)\n\npredictions = {}\n\nfor rid, items in tqdm(records.items(), desc=\"Records\"):\n    tmp_store = {}\n    scales_store = {}\n\n    for r in items:\n        lead = str(r.lead)\n        fs   = int(r.fs)\n        n    = int(r.number_of_rows)\n\n        # PCA template if available, else median template\n        tpl_beat = lead_template.get(lead, zscore(lead_templates_raw['II']))\n\n        # BPM selection (RF‑augmented)\n        best_bpm, y_best = choose_bpm(tpl_beat, fs, n)\n\n        # Prior BPM (median from training for this lead) branch\n        bpm_fixed = float(np.median(lead_bpms.get(lead, [75.0]))) if len(lead_bpms.get(lead, []))>0 else 75.0\n        y_fixed   = tile_template(tpl_beat, fs, n, bpm_fixed, amp=1.0)\n\n        # Plain mean template branch = PCA template stretched to n\n        y_mean = resample_to_length(tpl_beat, n)\n\n        # Light low‑pass\n        y_best  = apply_lowpass(y_best, fs, cutoff=15.0, order=2)\n        y_fixed = apply_lowpass(y_fixed, fs, cutoff=15.0, order=2)\n        y_mean  = apply_lowpass(y_mean, fs,  cutoff=15.0, order=2)\n\n        # Normalize to shape space\n        B, F, M = zscore(y_best), zscore(y_fixed), zscore(y_mean)\n\n        # Micro‑ensemble\n        w = ENSEMBLE_W / (np.sum(ENSEMBLE_W) + 1e-8)\n        y_syn = (w[0]*B + w[1]*F + w[2]*M).astype(np.float32)\n\n        # Optional CNN denoiser to subtly correct shape\n        if 'torch' in globals() and ENABLE_CNN_DENOISER and ('denoiser' in globals()) and (denoiser is not None):\n            with torch.no_grad():\n                inp = torch.tensor(y_syn[None,None,:])\n                y_syn = denoiser(inp).numpy().reshape(-1).astype(np.float32)\n\n        tmp_store[lead] = y_syn\n        scales_store[lead] = (fs, n)\n\n    # Einthoven blending\n    if EINTHOVEN_BLEND_W > 0.0 and ('I' in tmp_store) and ('II' in tmp_store):\n        for dlead in ['III', 'aVR', 'aVL', 'aVF']:\n            if dlead not in scales_store:\n                continue\n            _, n_d = scales_store[dlead]\n            yI_rs  = resample_to_length(tmp_store['I'],  n_d)\n            yII_rs = resample_to_length(tmp_store['II'], n_d)\n            derived_all = derive_limb_leads_from_I_II(yI_rs, yII_rs)\n            ydrv = zscore(derived_all[dlead])\n            tmp_store[dlead] = soft_blend(tmp_store.get(dlead, ydrv), ydrv, EINTHOVEN_BLEND_W)\n\n    # Final scaling per lead with adaptive min/max (if learned), then map to [MIN,MAX]\n    for lead, y in tmp_store.items():\n        scaler = lead_scalers.get(lead)\n        if scaler is not None:\n            lr_min, lr_max = scaler\n            s = np.std(y); i = np.subtract(*np.percentile(y, [75,25]))\n            est_min = float(lr_min.predict([[s,i]])[0])\n            est_max = float(lr_max.predict([[s,i]])[0])\n            ls = {'min': min(est_min, est_max), 'max': max(est_min, est_max)}\n            y_scaled = scale_to_lead_range(y, ls, MIN_VAL, MAX_VAL)\n        else:\n            y_scaled = scale_to_lead_range(y, lead_stats.get(lead), MIN_VAL, MAX_VAL)\n        predictions[(rid, lead)] = y_scaled.astype(np.float32)\n\n# Write submission\nrows = []\nfor r in test.itertuples(index=False):\n    rid   = int(r.id)\n    lead  = str(r.lead)\n    n     = int(r.number_of_rows)\n    y     = predictions[(rid, lead)]\n    if len(y) != n:\n        y = resample_to_length(y, n)\n    for i in range(n):\n        rows.append((f\"{rid}_{i}_{lead}\", float(y[i])))\n\nsub = pd.DataFrame(rows, columns=['id','value'])\nassert set(sub.columns) == {'id','value'}\nassert sub['value'].between(MIN_VAL-1e-8, MAX_VAL+1e-8).all()\nassert sub['id'].nunique() == len(sub)\nsub.to_csv('submission.csv', index=False)\nprint(\"✅ submission.csv written successfully!\")\n\n# Peek\nprint(sub.head(10).to_string(index=False))","metadata":{"_uuid":"6b662475-5e95-4446-9f56-1244148cb67c","_cell_guid":"295017a2-ed96-4225-9d2f-7c720885557d","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-10-27T16:45:25.779745Z","iopub.execute_input":"2025-10-27T16:45:25.780114Z","iopub.status.idle":"2025-10-27T16:45:27.495953Z","shell.execute_reply.started":"2025-10-27T16:45:25.780080Z","shell.execute_reply":"2025-10-27T16:45:27.494971Z"}},"outputs":[],"execution_count":null}]}