Upload folder using huggingface_hub

.gitattributes

@@ -33,3 +33,65 @@ saved_model/**/* filter=lfs diff=lfs merge=lfs -text
 *.zip filter=lfs diff=lfs merge=lfs -text
 *.zst filter=lfs diff=lfs merge=lfs -text
 *tfevents* filter=lfs diff=lfs merge=lfs -text
+assets/logo.png filter=lfs diff=lfs merge=lfs -text
+example_inputs/test.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/confusion_matrix-checkpoint.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/confusion_matrix_normalized-checkpoint.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/labels_correlogram-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/P_curve-checkpoint.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/results-checkpoint.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/val_batch0_labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/val_batch0_pred-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/.ipynb_checkpoints/val_batch1_labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/confusion_matrix.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/confusion_matrix_normalized.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/F1_curve.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/labels_correlogram.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/P_curve.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/PR_curve.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/R_curve.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/results.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch0.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch1.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch2.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch237330.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch237331.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/train_batch237332.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch0_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch0_pred.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch1_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch1_pred.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch2_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_lead/val_batch2_pred.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/.ipynb_checkpoints/labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/.ipynb_checkpoints/train_batch0-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/.ipynb_checkpoints/train_batch1-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/.ipynb_checkpoints/train_batch2-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/labels_correlogram.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/train_batch0.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/train_batch1.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_patch/train_batch2.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/results-checkpoint.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/train_batch50402-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/val_batch0_pred-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/val_batch1_labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/.ipynb_checkpoints/val_batch2_labels-checkpoint.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/labels_correlogram.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/results.png filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch0.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch1.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch2.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch50400.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch50401.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/train_batch50402.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch0_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch0_pred.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch1_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch1_pred.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch2_labels.jpg filter=lfs diff=lfs merge=lfs -text
+models/digitization_models/yolo11_pulse/val_batch2_pred.jpg filter=lfs diff=lfs merge=lfs -text

.gitignore

@@ -0,0 +1,15 @@
+# Ignore model checkpoints (heavy weights)
+models/
+
+# Streamlit run output directories
+output/
+
+# Python cache files
+__pycache__/
+*.pyc
+
+# Local Streamlit config files
+.streamlit/
+
+# Legacy/archive code directories (retained locally but excluded from version control)
+archive/

.streamlit/config.toml

@@ -0,0 +1,15 @@
+[theme]
+primaryColor = "#E63946"
+backgroundColor = "#FFFFFF"
+secondaryBackgroundColor = "#F1F5F9"
+textColor = "#1E293B"
+font = "sans serif"
+
+[server]
+headless = true
+enableCORS = false
+enableXsrfProtection = false
+maxUploadSize = 200
+
+[runner]
+fastReruns = true

LICENSE

@@ -0,0 +1,18 @@
+NON-COMMERCIAL ACADEMIC AND RESEARCH LICENSE AGREEMENT
+
+Copyright (c) 2026 SCAI Lab (Spinal Cord Injury and Artificial Intelligence Lab). All rights reserved.
+
+Redistribution and use of this software (in source or binary forms, with or without modification) and its associated trained model weights/checkpoints are permitted solely for academic, educational, and non-profit research purposes, subject to the following terms and conditions:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. The software, source code, binary forms, and all associated model checkpoints/weights (located under the `models/` directory or downloaded separately) may NOT be used for any commercial purpose, including but not limited to:
+   - Integration into a commercial or profit-generating product, platform, or service.
+   - Use in paid consulting, corporate research, clinical trial services, or commercial patient diagnostics.
+   - Sublicensing, renting, leasing, selling, or distributing the software or derivatives for financial gain.
+
+4. Any publications, patents, research papers, presentations, or software derivatives that utilize this software or the associated digitization/classification pipeline must credit the original authors and cite this repository.
+
+THIS SOFTWARE AND MODEL WEIGHTS ARE PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES (INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT) ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -0,0 +1,375 @@
+# ⚡ ECGLight: ECG Digitization & Classification Dashboard
+
+<!--
+<table align="center" border="0" style="border-collapse: collapse; border: none;">
+  <tr style="border: none;">
+    <td align="center" valign="middle" style="border: none; padding-right: 40px;">
+      <img src="assets/logo.png" alt="ECG Digitization & Classification Logo" width="180px" style="border-radius: 14px; box-shadow: 0 4px 16px rgba(15, 23, 42, 0.08);" />
+    </td>
+    <td align="center" valign="middle" style="border: none;">
+      <img src="assets/scai_lab_logo.svg" alt="SCAI Lab Logo" height="100px" />
+    </td>
+  </tr>
+</table>
+-->
+
+An advanced, interactive Streamlit web workstation designed to convert printed/photographed 12-lead paper ECG reports into high-resolution digitized signals and carry out different types of classification over them. The suite is engineered to run with low computational resource requirements, operating seamlessly on a standard consumer laptop GPU (via CUDA) or running completely on CPU.
+
+## <a id="table-of-contents"></a>📌 Table of Contents
+
+- [🖥️ Web Dashboard Workstation Overview](#web-dashboard-workstation-overview)
+- [🚀 Installation & Setup](#installation-setup)
+- [🧠 Pre-Trained Classifiers & Tasks](#pre-trained-classifiers-tasks)
+- [🚀 Command Line Usage](#command-line-usage)
+- [📁 Repository Structure & Directory Organization](#repository-structure-directory-organization)
+- [📷 How It Works: Signal Digitization](#how-it-works-signal-digitization)
+- [📈 How It Works: Signal Analysis & Visualization](#how-it-works-signal-analysis-visualization)
+- [⚡ How It Works: Heartbeat Segmentation](#how-it-works-heartbeat-segmentation)
+- [🧠 How It Works: Cardiac Classification](#how-it-works-cardiac-classification)
+- [🤝 Collaborating Institutions](#collaborating-institutions)
+- [📄 Citation](#citation)
+- [👥 Authors & Contact](#authors-contact)
+- [📄 License](#license)
+
+## <a id="web-dashboard-workstation-overview"></a>🖥️ Web Dashboard Workstation Overview
+
+The dashboard provides a premium, responsive user interface designed for research, education, and clinical workflow exploration. It coordinates the digitization and classification pipelines into a unified, lightweight web application.
+
+### Workstation Modules
+
+1. **📷 ECG Image Digitizer**:
+   - Upload any scanned or photographed ECG image (`.png`, `.jpg`, `.jpeg`).
+   - Run the sequential YOLOv11 pipeline step-by-step with real-time progress indicators.
+   - Outputs a summary of detected leads and total samples.
+   - Automatically saves the digitized CSV to disk under `output/digitization/latest_digitized.csv` for downstream consumption.
+   
+2. **📈 ECG Signal Viewer**:
+   - Visualizes multi-channel ECG signals interactively using client-side native Streamlit line charts (supporting zoom, pan, and hover tooltips).
+   - Supports stacked subplots (with distinct clinical colors for each lead: clinical red, teal, deep blue, yellow, purple, etc.) or overlaid graphs.
+   - Displays statistical summaries (mean, standard deviation, min/max, range) and allows row-by-row signal previewing.
+
+3. **❤️ ECG Classification**:
+   - Predicts cardiac conditions using pre-trained ensemble and deep learning classifiers.
+   - Automatically segments raw signals into heartbeats around R-peaks using the Pan-Tompkins algorithm before running inference.
+   - **Inference Mode (No Ground-Truth)**: If the uploaded CSV lacks diagnostic labels, the dashboard displays a downloadable **Predictions Table** detailing predicted class diagnoses and model confidence probabilities.
+   - **Evaluation Mode (With Ground-Truth)**: If labels are present, the page calculates and plots performance metrics (Accuracy, F1-Score, Sensitivity, Specificity, Confusion Matrix).
+
+### 🔄 Workflow
+
+The workstation coordinates the pipeline through four distinct steps: **Digitization**, **Analysis**, **Segmentation**, and **Classification**. The technical workflows for each process are documented below in their respective **How It Works** sections.
+
+---
+
+## <a id="installation-setup"></a>🚀 Installation & Setup
+
+### Prerequisites
+- Python 3.9
+- CUDA-capable GPU recommended (automatically falls back to CPU if unavailable).
+
+### Conda Environment & Model Setup
+
+1. Clone the repository and navigate to the project directory:
+   ```bash
+   git clone https://github.com/scai-lab/ECG-Digitization-Classification.git
+   cd ECG-Digitization-Classification
+   ```
+
+2. Create the conda environment using the provided `environment.yml` configuration:
+   ```bash
+   conda env create -f environment.yml
+   conda activate infer
+   ```
+
+   > [!IMPORTANT]
+   > **Windows Compatibility & TensorFlow Setup**:
+   > If you are on Windows and encounter native runtime loading failures (`ImportError: DLL load failed while importing _pywrap_tensorflow_internal: A dynamic link library (DLL) initialization routine failed`), you need to install a stable version pairing of TensorFlow and Protobuf:
+   > ```bash
+   > pip install tensorflow==2.15.0 protobuf==4.25.3
+   > ```
+   > *(Make sure no background Streamlit or Python tasks are running when executing this command, to prevent file locking issues on `.pyd` libraries).*
+
+3. **Download Pre-Trained Model Weights**:
+   Due to their file sizes, the YOLO detection checkpoints and pre-trained classifiers are hosted externally. Download the `models/` directory from the link below and place it directly in the root of the project:
+   
+   👉 **[Download Pre-Trained Models Directory (ETH Zürich Polybox)](https://polybox.ethz.ch/index.php/s/GDACstPtsoTrrWH)**
+   
+   Once extracted, verify that the weights are located inside the directory tree structure:
+   ```text
+   models/
+   ├── digitization_models/
+   │   ├── yolo11_full/weights/best.pt
+   │   ├── yolo11_lead/weights/best.pt
+   │   ├── yolo11_pulse/weights/best.pt
+   │   └── yolo11_patch/weights/best.pt
+   └── classifier_models/
+       ├── mi_vs_normal_segmented/
+       ├── omi_vs_nonomi/
+       └── ecg_surgery/
+   ```
+
+Key packages installed by the environment: `torch 2.7`, `ultralytics 8.3`, `opencv-python 4.11`, `scikit-image 0.24`, `wfdb 4.3`, `patched-yolo-infer 1.3.8`, `sktime`, `streamlit`.
+
+---
+
+## <a id="pre-trained-classifiers-tasks"></a>🧠 Pre-Trained Classifiers & Tasks
+
+The classification engine supports three diagnostic tasks using the pre-trained weights in `classifier_models/`:
+
+| Classification Task | Model Type | Expected Input Shape | Test Accuracy | Positive Class |
+| :--- | :--- | :--- | :---: | :--- |
+| **Normal vs Myocardial Infarction (MI) - Segmented** | Arsenal | 12 leads × 140 timesteps | **92.3%** | `MYOCARDIAL_INFARCTION` |
+| **Occlusive MI (OMI) vs non-OMI** | Rocket | 12 leads × 141 timesteps | **88.9%** | `OMI` |
+| **Pre-Procedural vs Post-Procedural MI** | InceptionTime | 12 leads × 140 timesteps | **91.4%** | `pre-procedural MI` |
+
+- **Arsenal**: An ensemble of ROCKET classifiers utilizing random convolutional kernels to extract feature representations combined with ridge regression.
+- **Rocket**: Random Omni-directional Kernel Extraction (ROCKET) classifier, computing kernel convolutions quickly for high-dimensional time-series data.
+- **InceptionTime**: A deep convolutional network ensemble modeled on the Inception architecture, extracting multi-scale temporal features.
+
+---
+
+## <a id="command-line-usage"></a>🚀 Command Line Usage
+
+### Run Batch Digitization (`run_org.py`)
+
+The batch processing script processes nested hospital directories, exporting structured folders of digitized CSVs:
+
+1. Configure path variables at the top of `run_org.py`:
+   ```python
+   ORGANIZED_DIR = "../ecg_files/ECG_organized_all"   # Input dataset root
+   OUTPUT_DIR    = "../ecg_files/ECG_digitized"       # Mirrored CSV directory
+   CATEGORIES    = ["pre", "index", "post"]           # Categories to process
+   ```
+
+2. Run the script:
+   ```bash
+   python run_org.py
+   ```
+
+### Run Model Inference (`run_inference.py`)
+
+Execute predictions directly on digitized data from the command line using `run_inference.py` located in `archive/classification/`:
+
+```bash
+# MI vs Normal Segmented heartbeat classification
+python archive/classification/run_inference.py --model mi_vs_normal_segmented --input data/ptb_xl/segmented_heartbeats.csv
+
+# OMI vs non-OMI classification
+python archive/classification/run_inference.py --model omi_vs_nonomi --input data/ecg_matrix_omi_segmented_50_150_90.csv
+
+# Custom output file path
+python archive/classification/run_inference.py --model ecg_surgery --input data/ecg_surgery_segmented_50_150_70.csv --output results/surgery_preds.csv
+```
+
+---
+
+## <a id="repository-structure-directory-organization"></a>📁 Repository Structure & Directory Organization
+
+The repository is structured to maintain a clean root directory, moving utility runners, UI views, model checkpoints, and legacy/training scripts into distinct modules:
+
+```
+.
+├── app.py                          # Streamlit application main router
+├── config.py                       # Centralized configuration and model registry
+├── digitization.py                 # Core ECGImage extraction pipeline class
+├── environment.yml                 # Conda environment dependency file
+├── README.md                       # Comprehensive repository documentation
+│
+├── backend/                        # Dashboard background execution adapters
+│   ├── __init__.py                 # Backend package declaration
+│   ├── digitization_runner.py      # YOLO loader and single-image processor
+│   └── classification_runner.py    # Pre-trained model loader and preprocessor
+│
+├── utils/                          # Streamlit front-end page components
+│   ├── __init__.py                 # Utils package declaration
+│   ├── branding.py                 # Sidebar titles, headers, and footer logos
+│   ├── css.py                      # Custom clinical theme and grid background CSS
+│   ├── hardware.py                 # Displays CPU/GPU hardware properties (cached)
+│   ├── page_digitizer.py           # Front-end for the ECG Digitizer page
+│   ├── page_csv_viewer.py          # Front-end for the interactive Signal Viewer
+│   └── page_classifier.py          # Front-end for the Classification workstation
+│
+├── models/                         # Relocated YOLO checkpoints and classifiers
+│   ├── digitization_models/        # YOLO v11 checkpoints for digitization
+│   │   ├── yolo11_full/            # YOLO Bounding boxes
+│   │   ├── yolo11_lead/            # YOLO Lead names
+│   │   ├── yolo11_pulse/           # YOLO Reference pulses
+│   │   └── yolo11_patch/           # YOLO Waveform segmentations
+│   │
+│   └── classifier_models/          # Bundled pre-trained diagnostic classifiers
+│       ├── mi_vs_normal_segmented/ # Pre-trained Arsenal model (segmented beats)
+│       ├── omi_vs_nonomi/          # Pre-trained Rocket model (segmented beats)
+│       └── ecg_surgery/            # Pre-trained InceptionTime model (segmented beats)
+│
+└── archive/                        # Archived developer, training, and legacy scripts
+    └── classification/             
+        ├── train_and_save_models.py# Script used to compile pre-trained models
+        ├── run_inference.py        # Independent CLI inference execution script
+        ├── run_classification.py   # Baseline MLP classifier pipeline
+        ├── run_benchmarking.py     # Comparison benchmarking suite
+        ├── run_lead_importance_test.py# Individual lead performance evaluator
+        ├── feature_analysis.py     # Original all-in-one analysis script
+        ├── aggregate_subject_metrics.py# Multi-subject performance aggregator
+        ├── dataset_curate.py       # Local dataset curation utility
+        └── re_plotter.py           # Advanced Gaussian signal generator & visualizer
+```
+
+---
+
+## <a id="how-it-works-signal-digitization"></a>📷 How It Works: Signal Digitization
+
+The core class [digitization.py](file:///d:/Projects/ECGLight/digitization.py) operates a multi-stage sequential computer vision pipeline to translate raster images into digitized signals:
+
+```mermaid
+graph TD
+    A[ECG Image Upload] --> B[Preprocessing: Otsu & Blurring]
+    B --> C[YOLOv11 Detection & Segmentation]
+    subgraph YOLOv11 Models
+        C1[yolo11_full: Lead Boundaries]
+        C2[yolo11_lead: Text Name Labels]
+        C3[yolo11_pulse: Calibration Pulses]
+        C4[yolo11_patch: Waveform Segments]
+    end
+    C --> C1 & C2 & C3 & C4
+    C1 & C2 & C3 & C4 --> D[Hough Lines Calibration]
+    D --> E[K-Means Row & Column Grid Construction]
+    E --> F[Anti-Leakage Connected Components Filter]
+    F --> G[Centroid Trace & Resampling to 500Hz]
+    G --> H[Export latest_digitized.csv]
+```
+
+1. **Preprocessing**: Cleans the scanned image using shadow-removal masks, Otsu binarization, and Gaussian blurring to isolate ink lines from paper textures.
+2. **YOLO Segmentation**: Applies a patched YOLO segmentation model at three crop scales (`4×`, `4.5×`, and `5×` height) to isolate individual lead waveform contours.
+3. **Sequential Detections**: Runs three YOLO models in parallel:
+   - `yolo11_full`: Bounding boxes for the 12 lead channels.
+   - `yolo11_lead`: Text labels representing lead names (I, II, aVR...).
+   - `yolo11_pulse`: Bounding boxes for the calibration reference pulses (typically 1mV high, representing vertical scale).
+4. **Scale Calibration**: Fits Hough lines to the calibration pulse boundaries. The pixel height determines the voltage scale (`volt/pixel`), while the width determines the time scale (`time/pixel`).
+5. **Grid Construction**: Employs K-Means clustering on lead coordinates to map rows and columns, automatically parsing standard Cabrera orders and grid formats (3×4, 4×3, 6×2, 12×1).
+6. **Signal Extraction & Post-Processing**: Traces contours to extract raw pixel centroids, performs baseline correction, applies linear interpolation to bridge gaps, and resamples to a standard **500 Hz** frequency calibrated in **millivolts (mV)**. In addition, an **anti-leakage component filter** is executed per cell crop using connected components analysis to automatically identify the primary waveform trace and strip out smaller, boundary-adjacent components (leaked signals from neighboring leads) that sit far from the row baseline.
+
+---
+
+## <a id="how-it-works-signal-analysis-visualization"></a>📈 How It Works: Signal Analysis & Visualization
+
+Once continuous signals are extracted, the dashboard runs analytical tasks and displays interactive previews:
+
+```mermaid
+graph TD
+    A[Upload Digitized CSV] --> B[Parse Lead Voltages & Timestamps]
+    B --> C[Vega-Lite Interactive Visualizer]
+    C --> C1[Render Stacked Leads]
+    C --> C2[Render Overlaid Signals]
+    B --> D[Compute Signal Statistics: Mean, SD, Min/Max]
+    D --> E[Display Summary Dataframes & Row Previews]
+```
+
+1. **Parse Signals**: Reads digitized CSV format, validating lead names and timestamps.
+2. **Vega-Lite Visualization**: Renders interactive charts supporting native client-side zoom, pan, and hover tooltips for all channels.
+3. **Signal Statistics**: Automatically computes statistical characteristics (mean, standard deviation, min/max values) for each lead.
+
+---
+
+## <a id="how-it-works-heartbeat-segmentation"></a>⚡ How It Works: Heartbeat Segmentation
+
+To prepare continuous digitized signals for the classification models, the pipeline runs the Pan-Tompkins R-peak detection algorithm:
+
+```mermaid
+graph TD
+    A[Digitized 500Hz Signal] --> B[Bandpass Filter 5-15Hz]
+    B --> C[Derivative Filter]
+    C --> D[Squaring Operation]
+    D --> E[Moving Window Integration]
+    E --> F[Adaptive Thresholding & R-Peak Search]
+    F --> G[Extract 140-sample Beats: 50ms pre-R, 150ms post-R]
+    G --> H[Max-Absolute Voltage Normalization]
+```
+
+1. **Filtering**: The Lead II signal is filtered via a bandpass filter (5-15 Hz) to suppress muscle noise, baseline wander, and T-wave interference.
+2. **Differentiation**: Computes the slope of the signal to highlight the rapid change in the QRS complex.
+3. **Squaring**: Performs point-by-point squaring to amplify QRS slopes while attenuating smaller waves.
+4. **Integration**: A moving window integrator (typically 150ms wide) compiles the slope information into a peak window.
+5. **Adaptive Thresholding & Peak Search**: Dynamically computes threshold constants based on average noise and signal levels, locating R-peaks.
+6. **Beat Windowing**: Extracts a localized heartbeat around each R-peak (typically extending 50ms before and 150ms after the peak), normalizes the voltage per heartbeat using max-absolute scaling, and truncates/pads the resulting segments to the target model input width (e.g. 140 or 141 timesteps).
+
+---
+
+## <a id="how-it-works-cardiac-classification"></a>🧠 How It Works: Cardiac Classification
+
+The heartbeat segment tensors are evaluated using pre-trained time-series classification models:
+
+```mermaid
+graph TD
+    A[Segmented Heartbeats] --> B[Numpy3D Reshaping: N_instances × 12_leads × N_timesteps]
+    B --> C[Select Classification Task]
+    subgraph Model Registry
+        C1[Normal vs MI: Arsenal]
+        C2[OMI vs non-OMI: Rocket]
+        C3[Pre vs Post-Procedural MI: InceptionTime]
+    end
+    C --> C1 & C2 & C3
+    C1 & C2 & C3 --> D[Load Pre-Trained Pickled Estimator]
+    D --> E[Predict Class Labels & Probabilities]
+    E --> F[Generate Downloadable Predictions CSV]
+```
+
+1. **Numpy3D Formatting**: Formats the heartbeat segments into a standard `sktime` `Numpy3D` tensor with shape `(N_instances, 12_leads, N_timesteps)`.
+2. **Dynamic Task Selection**: Loads the pre-trained pickled model corresponding to the selected classification task.
+3. **Model Inference**: Evaluates the model to compute class predictions and probability confidences.
+4. **Result Generation**: Automatically builds downloadable prediction tables and calculates performance metrics if ground-truth labels are present in the dataset.
+
+---
+
+## <a id="collaborating-institutions"></a>🤝 Collaborating Institutions
+
+This project was developed in collaboration with:
+
+- [ETH Zürich](https://ethz.ch)
+- [Istituto Cardiocentro Ticino (EOC)](https://www.cardiocentro.org)
+- [Università della Svizzera italiana (USI)](https://www.usi.ch)
+- [Università della Campania Luigi Vanvitelli](https://www.unicampania.it)
+
+<p align="center">
+  <a href="https://ethz.ch" target="_blank">
+    <img src="assets/ETH_Zürich_Logo_black.svg.png" alt="ETH Zürich" height="30px" style="vertical-align: middle; margin: 0 15px;" />
+  </a>
+  &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
+  <a href="https://www.cardiocentro.org" target="_blank">
+    <img src="assets/eoc_logo.png" alt="Istituto Cardiocentro Ticino (EOC)" height="35px" style="vertical-align: middle; margin: 0 15px;" />
+  </a>
+  &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
+  <a href="https://www.usi.ch" target="_blank">
+    <img src="assets/usi_logo.png" alt="USI" height="35px" style="vertical-align: middle; margin: 0 15px;" />
+  </a>
+  &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
+  <a href="https://www.unicampania.it" target="_blank">
+    <img src="assets/Logo_Vanvitelli_university.svg.png" alt="Università della Campania Luigi Vanvitelli" height="35px" style="vertical-align: middle; margin: 0 15px;" />
+  </a>
+</p>
+
+---
+
+## <a id="citation"></a>📄 Citation
+
+```bibtex
+@article{natraj2026ecglight,
+  title={ECGLight: Compute-Light Framework For Paper ECG Digitization and Myocardial Infarction Screening},
+  author={Natraj, Shreyasvi and Achtari, Cyrus and Gragnano, Felice and Milzi, Andrea and Valgimigli, Marco and Paez-Granados, Diego},
+  journal={arXiv preprint arXiv:2607.07683},
+  year={2026},
+  url={https://arxiv.org/abs/2607.07683},
+  doi={10.48550/arXiv.2607.07683}
+}
+```
+
+---
+
+## <a id="authors-contact"></a>👥 Authors & Contact
+
+- **Shreyasvi Natraj** — [snatraj@ethz.ch](mailto:snatraj@ethz.ch)
+- **Cyrus Achtari**
+
+---
+
+## <a id="license"></a>📄 License
+
+This project is released under the **Non-Commercial Academic and Research License Agreement**. Please refer to the [LICENSE](file:///d:/Projects/ECGLight/LICENSE) file in the repository root for the full licensing terms. The codebase and trained model weights are provided free of charge for personal, academic, and non-profit research use only. Commercial use is strictly prohibited.

app.py

@@ -0,0 +1,114 @@
+# app.py
+"""
+Main entry point for the ECGLight Dashboard.
+Slim router that delegates to page modules in utils/.
+"""
+
+import os
+import sys
+import streamlit as st
+
+# Add project root to path
+sys.path.append(os.path.abspath(os.path.dirname(__file__)))
+
+# Import configuration
+import config
+
+# Import utility modules (CSS, branding, hardware, pages)
+from utils.css import inject_custom_css
+from utils.branding import render_sidebar_logo, render_logo_footer
+from utils.hardware import render_sidebar_hardware
+from utils import page_digitizer, page_csv_viewer, page_classifier
+
+# Backend runners are imported lazily inside page routing blocks below
+# to avoid loading torch/sklearn/scipy when the user only needs the CSV viewer.
+
+# ==============================================================================
+# PAGE CONFIG (must be the first Streamlit command)
+# ==============================================================================
+st.set_page_config(
+    page_title="ECGLight",
+    page_icon="⚡",
+    layout="wide",
+    initial_sidebar_state="expanded"
+)
+
+# ==============================================================================
+# INJECT CSS THEME
+# ==============================================================================
+inject_custom_css()
+
+# ==============================================================================
+# SIDEBAR
+# ==============================================================================
+# SCAI Lab logo at the top of the sidebar
+render_sidebar_logo()
+
+st.sidebar.markdown(
+    '<h1 style="font-size: 1.8rem; margin-bottom: 5px; font-weight: 800; letter-spacing: -0.5px; text-align: center; padding-right: 1.5rem;">'
+    '⚡ <span class="glow-text">ECGLight</span></h1>',
+    unsafe_allow_html=True
+)
+st.sidebar.markdown(
+    '<p class="section-subtitle" style="text-align: center;">Lightweight Digitization & Classification Suite</p>',
+    unsafe_allow_html=True
+)
+
+st.sidebar.markdown("---")
+
+# Navigation (CSV Signal Viewer above Classification)
+page = st.sidebar.radio(
+    "🧭 Navigation Pages",
+    ["📷 ECG Image Digitizer", "📈 ECG Signal Viewer", "❤️ ECG Classification"],
+    index=0
+)
+
+# Keep track of active page to reset session state on page navigation
+if "current_page" not in st.session_state:
+    st.session_state["current_page"] = page
+
+if st.session_state["current_page"] != page:
+    # Clear session state keys to reset state on page change to prevent slowing things down,
+    # but we DO NOT clear "latest_digitized_path" to share the default digitized CSV.
+    keys_to_clear = [
+        "digitized_df", "digitized_name", "last_uploaded_name",
+        "clf_results", "clf_logs", "selected_task",
+        "csv_viewer_upload", "classifier_csv_upload"
+    ]
+    for key in keys_to_clear:
+        st.session_state.pop(key, None)
+    
+    st.session_state["current_page"] = page
+    st.rerun()
+
+st.sidebar.markdown("---")
+
+# Hardware info (cached — only runs once per session)
+render_sidebar_hardware()
+
+# ==============================================================================
+# YOLO MODEL CHECK (used by Page 1)
+# ==============================================================================
+yolo_configured = all([
+    config.YOLO_BOX_MODEL_PATH,
+    config.YOLO_LEAD_NAME_MODEL_PATH,
+    config.YOLO_PULSE_MODEL_PATH,
+    config.YOLO_SEGMENTATION_MODEL_PATH
+])
+
+# ==============================================================================
+# PAGE ROUTING (lazy backend imports — only load heavy modules when needed)
+# ==============================================================================
+if page == "📷 ECG Image Digitizer":
+    from backend import digitization_runner
+    page_digitizer.render(config, digitization_runner, yolo_configured)
+    render_logo_footer()
+
+elif page == "📈 ECG Signal Viewer":
+    page_csv_viewer.render()
+    render_logo_footer()
+
+elif page == "❤️ ECG Classification":
+    from backend import classification_runner
+    page_classifier.render(config, classification_runner)
+    render_logo_footer()

archive/classification/aggregate_subject_metrics.py

@@ -0,0 +1,96 @@
+import pandas as pd
+import numpy as np
+import os
+import re
+from sklearn.metrics import accuracy_score, precision_score, f1_score
+
+def get_base_subject_id(subject_id_with_beat):
+    """
+    Extracts the base subject ID by removing suffixes like _beatN,  beatN, -index_beatN, etc.
+    """
+    # Regex to catch variations like _beat1,  beat1, -index_beat1, INDEX_beat1, etc.
+    # It looks for 'beat' preceded by common delimiters and captures everything before it.
+    match = re.search(r'^(.*?)(?:[_\s-](?:index|INDEX|post|POST|checkpoint|POST PCI|POST PCI \d|POST \d)*)?(?:[_\s-]beat\d+)$', str(subject_id_with_beat))
+    if match:
+        return match.group(1).strip()
+    
+    # Fallback: just split by 'beat' if regex is too specific
+    if 'beat' in str(subject_id_with_beat):
+        return str(subject_id_with_beat).split('beat')[0].rstrip('_ -').strip()
+    
+    return str(subject_id_with_beat).strip()
+
+def calculate_metrics(df_subject):
+    """
+    Calculates metrics for subject-level predictions.
+    """
+    # Ground truth: since all beats for a subject have the same label, we just take the first
+    y_true = df_subject['true_label_numeric'].values
+    
+    # Prediction based on averaged probability for Pre-Procedural MI
+    # Assuming probability_Pre-Procedural MI >= 0.5 is Pre-Procedural MI (1)
+    y_pred = (df_subject['avg_prob_pre_procedural'] >= 0.5).astype(int)
+    
+    acc = accuracy_score(y_true, y_pred)
+    prec = precision_score(y_true, y_pred, zero_division=0)
+    f1 = f1_score(y_true, y_pred, zero_division=0)
+    
+    return acc, prec, f1
+
+def main():
+    prediction_dir = r'output/benchmark_results/ecg_matrix_dataset_segmented/model_predictions'
+    
+    if not os.path.exists(prediction_dir):
+        print(f"Directory not found: {prediction_dir}")
+        return
+
+    # Filter out the combined file if it exists
+    files = [f for f in os.listdir(prediction_dir) if f.endswith('.csv') and f != 'all_models_predictions_combined.csv']
+    
+    results = []
+    
+    print(f"{'Model':<20} | {'Accuracy':<10} | {'Precision':<10} | {'F1 Score':<10}")
+    print("-" * 60)
+
+    for file in files:
+        model_name = file.replace('_predictions.csv', '')
+        file_path = os.path.join(prediction_dir, file)
+        
+        try:
+            df = pd.read_csv(file_path)
+            
+            # 1. Extract base subject ID
+            df['base_subject_id'] = df['subject_id'].apply(get_base_subject_id)
+            
+            # 2. Group by base subject ID
+            # We average the probability and take the first true label (since it should be invariant for a subject)
+            subject_group = df.groupby('base_subject_id').agg({
+                'probability_Pre-Procedural MI': 'mean',
+                'true_label_numeric': 'first'
+            }).reset_index()
+            
+            subject_group.rename(columns={'probability_Pre-Procedural MI': 'avg_prob_pre_procedural'}, inplace=True)
+            
+            # 3. Calculate Metrics
+            acc, prec, f1 = calculate_metrics(subject_group)
+            
+            results.append({
+                'Model': model_name,
+                'Accuracy': acc,
+                'Precision': prec,
+                'F1 Score': f1
+            })
+            
+            print(f"{model_name:<20} | {acc:<10.4f} | {prec:<10.4f} | {f1:<10.4f}")
+            
+        except Exception as e:
+            print(f"Error processing {file}: {e}")
+
+    # Save summary to CSV
+    summary_df = pd.DataFrame(results)
+    summary_path = os.path.join(prediction_dir, 'subject_level_metrics_summary.csv')
+    summary_df.to_csv(summary_path, index=False)
+    print(f"\nSummary saved to: {summary_path}")
+
+if __name__ == "__main__":
+    main()

archive/classification/dataset_curate.py

@@ -0,0 +1,66 @@
+import os
+import pandas as pd
+from pathlib import Path
+
+def process_ecg_data_wide_format(root_folder):
+    """
+    Process ECG data in wide format (one row per timestamp with all leads as columns)
+    """
+    lead_names = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+    
+    all_data = []
+    base_path = Path(root_folder)
+    
+    # Process both classes
+    for class_name, class_folder in [("Pre-Procedural MI", "Pre"), 
+                                    ("Post-Procedural MI", "Post")]:
+        class_path = base_path / class_folder
+        
+        if class_path.exists():
+            for csv_file in class_path.glob("*.csv"):
+                try:
+                    df = pd.read_csv(csv_file)
+                    subject_id = csv_file.stem.replace('_hr', '')
+                    
+                    # Add metadata columns
+                    df['subject_id'] = subject_id
+                    df['class'] = class_name
+                    df['filename'] = csv_file.name
+                    df['timestamp'] = range(len(df))
+                    
+                    all_data.append(df)
+                    print(f"Processed {csv_file.name}")
+                    print(class_name)
+                    
+                except Exception as e:
+                    print(f"Error processing {csv_file.name}: {str(e)}")
+    
+    if all_data:
+        # Combine all DataFrames
+        combined_df = pd.concat(all_data, ignore_index=True)
+        
+        # Reorder columns to have metadata first
+        metadata_cols = ['subject_id', 'timestamp', 'class', 'filename']
+        signal_cols = [col for col in combined_df.columns if col not in metadata_cols]
+        combined_df = combined_df[metadata_cols + signal_cols]
+        
+        return combined_df
+    return None
+
+# Usage for wide format:
+if __name__ == "__main__":
+    folder_path = "."
+    
+    if os.path.exists(folder_path):
+        print("Processing ECG data in wide format...")
+        wide_df = process_ecg_data_wide_format(folder_path)
+        
+        if wide_df is not None:
+            wide_df.to_csv("combined_ecg_data_wide.csv", index=False)
+            print("Wide format data saved to combined_ecg_data_wide.csv")
+            print(f"Shape: {wide_df.shape}")
+            print(f"Columns: {list(wide_df.columns)}")
+            print("\nSample data:")
+            print(wide_df.head())
+    else:
+        print(f"Folder '{folder_path}' not found.")

archive/classification/feature_analysis.py

@@ -0,0 +1,93 @@
+import os
+import pandas as pd
+import numpy as np
+from utils.mlp_feature_downsample import (
+    perform_comprehensive_heartbeat_analysis,
+    perform_full_signal_analysis,
+    DEFAULT_HEARTBEAT_OUTPUT_DIR,
+    DEFAULT_TARGET_FS,
+)
+
+if __name__ == "__main__":
+    
+    # Create output directory
+    os.makedirs('output', exist_ok=True)
+    
+    # Use the processed CSV from the previous step
+    combined_input_csv = 'data/ecg_dataset_processed.csv'
+    
+    # Load and verify the data
+    df_combined = pd.read_csv(combined_input_csv)
+    print(f"Loaded combined dataset shape: {df_combined.shape}")
+    print(f"Unique subjects: {df_combined['subject_id'].nunique()}")
+
+    print("\n" + "="*50)
+    print("RUNNING SEGMENTED HEARTBEAT ANALYSIS (500 Hz)")
+    print("="*50)
+    
+    heartbeat_results = perform_comprehensive_heartbeat_analysis(
+        input_csv=combined_input_csv,
+        output_dir=DEFAULT_HEARTBEAT_OUTPUT_DIR,
+        pca_components=50,
+        shap_samples=100,
+        shap_background_size=50,
+        hidden_layer_sizes=(100, 50),
+        max_iter=1000,
+        target_fs=500,
+    )
+    
+    print("\n" + "="*50)
+    print("RUNNING FULL SIGNAL ANALYSIS (500 Hz)")
+    print("="*50)
+    
+    full_signal_results = perform_full_signal_analysis(
+        input_csv=combined_input_csv,
+        output_dir='output/full_signal_analysis',
+        pca_components=50,
+        shap_samples=100,
+        shap_background_size=50,
+        hidden_layer_sizes=(100, 50),
+        max_iter=1000,
+        target_fs=500,
+    )
+
+    # Print Summary Results
+    print("\n" + "#" * 60)
+    print("FINAL SUMMARY RESULTS")
+    print("#" * 60)
+    
+    if heartbeat_results is not None and 'metrics' in heartbeat_results:
+        m = heartbeat_results['metrics']
+        print("\n--- Segmented Heartbeat Analysis ---")
+        print(f"  Accuracy:    {m['accuracy']:.4f}")
+        print(f"  Precision:   {m['precision']:.4f}")
+        print(f"  Recall:      {m['recall']:.4f}")
+        print(f"  F1-score:    {m['f1']:.4f}")
+        print(f"  Sensitivity: {m['sensitivity']:.4f}")
+        print(f"  Specificity: {m['specificity']:.4f}")
+        print("=" * 50)
+    '''
+    ############################################## Supervised Classifier Benchmarking ##################################
+    #print("\nPerforming classifier benchmarking...")
+    benchmark_results, models = perform_benchmark(combined_input_csv)
+    
+    print("\nBenchmarking Results:")
+    print(benchmark_results)
+
+    # Save all results to a summary file
+    results_summary = {
+        #'kmeans_results': kmeans_results,
+        'mlp_results': mlp_results,
+        'benchmark_results': benchmark_results,
+    }
+    
+    # Save summary to file
+    summary_df = pd.DataFrame([
+        #{'model': 'K-Means', 'details': str(kmeans_results)},
+        #{'model': 'MLP', 'details': str(mlp_results)},
+        {'model': 'Benchmark', 'details': str(benchmark_results)},
+    ])
+    
+    summary_df.to_csv('output/analysis_summary.csv', index=False)
+    print("\nAnalysis summary saved to 'output/analysis_summary.csv'")
+    '''

archive/classification/prepare_multiclass_data.py

@@ -0,0 +1,71 @@
+import pandas as pd
+import os
+
+if __name__ == "__main__":
+    # ── CONFIGURATION ──
+    # Input CSV file path
+    input_csv = 'data/combined_ecg_matrix_data_segmented.csv'
+    
+    # Output CSV file path
+    output_csv = 'data/combined_ecg_matrix_data_filtered.csv'
+
+    # Option 1: Only include specific classes (set to None to use all)
+    # e.g. ['Pre', 'Post', 'Index']
+    include = None
+
+    # Option 2: Exclude specific classes (set to None to skip)
+    # e.g. ['Dimi']
+    exclude = None
+
+    # Option 3: Combine/merge classes under a new label (set to None to skip)
+    # e.g. {'Surgery': ['Pre', 'Post'], 'Other': ['Dimi', 'Discharge']}
+    combine = {'Post-Procedural MI': ['Discharge', 'Dimi', 'Post'], 'Pre-Procedural MI': ['Index', 'Pre']}
+
+    print(f"Loading data from '{input_csv}'...")
+    if not os.path.exists(input_csv):
+        print(f"Error: Input file '{input_csv}' not found.")
+        exit(1)
+
+    df = pd.read_csv(input_csv)
+    initial_rows = len(df)
+    print(f"Loaded {initial_rows} rows.")
+    
+    if 'class' not in df.columns:
+        print("Error: 'class' column not found in input CSV.")
+        exit(1)
+
+    print(f"\nOriginal classes: {sorted(df['class'].dropna().unique())}")
+    for cls in sorted(df['class'].dropna().unique()):
+        print(f"  {cls}: {len(df[df['class'] == cls])} rows")
+
+    # ── 1. Combine/Merge Classes ──
+    if combine:
+        print("\nApplying class combinations...")
+        for new_name, old_names in combine.items():
+            df.loc[df['class'].isin(old_names), 'class'] = new_name
+            print(f"  Combined {old_names} -> '{new_name}'")
+
+    # ── 2. Include Classes ──
+    if include:
+        print(f"\nIncluding only: {include}")
+        df = df[df['class'].isin(include)]
+
+    # ── 3. Exclude Classes ──
+    if exclude:
+        print(f"\nExcluding: {exclude}")
+        df = df[~df['class'].isin(exclude)]
+
+    if len(df) == 0:
+        print("\nERROR: No data remaining after class filtering!")
+        exit(1)
+
+    print(f"\nData processing complete. Remaining rows: {len(df)} ({(len(df)/initial_rows)*100:.1f}%)")
+    
+    final_classes = sorted(df['class'].dropna().unique())
+    print(f"\nFinal classes hierarchy: {final_classes}")
+    for cls in final_classes:
+        print(f"  {cls}: {len(df[df['class'] == cls])} rows")
+
+    os.makedirs(os.path.dirname(output_csv), exist_ok=True)
+    df.to_csv(output_csv, index=False)
+    print(f"\nSaved filtered data to '{output_csv}'")

archive/classification/re_plotter.py

@@ -0,0 +1,494 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import os
+from matplotlib.gridspec import GridSpec
+
+def replot_idealized_shap_all_leads(shap_values_csv, top_regions_csv, output_dir, n_timesteps=40, fs=500):
+    """
+    Replot idealized heartbeat with SHAP values for all leads below it in order of importance
+    """
+    
+    # Create output directory
+    os.makedirs(output_dir, exist_ok=True)
+    
+    # Load SHAP values and top regions
+    print("Loading SHAP values and top regions...")
+    shap_df = pd.read_csv(shap_values_csv)
+    top_regions_df = pd.read_csv(top_regions_csv)
+    
+    print(f"SHAP DataFrame shape: {shap_df.shape}")
+    print(f"Number of timesteps inferred from columns: {len([col for col in shap_df.columns if '_t' in col]) // 12}")  # 12 leads
+    
+    # Calculate mean absolute SHAP values for each feature
+    mean_shap_abs = np.abs(shap_df).mean(axis=0)
+    
+    # Extract unique leads and their overall importance
+    leads = []
+    lead_importance = {}
+    
+    for feature_name in mean_shap_abs.index:
+        if '_t' in feature_name:
+            lead = feature_name.split('_t')[0]
+            if lead not in lead_importance:
+                lead_importance[lead] = 0
+            lead_importance[lead] += mean_shap_abs[feature_name]
+    
+    # Sort leads by overall importance
+    sorted_leads = sorted(lead_importance.items(), key=lambda x: x[1], reverse=True)
+    lead_names = [lead for lead, importance in sorted_leads]
+    lead_importances = [importance for lead, importance in sorted_leads]
+    
+    print("Lead importance ranking:")
+    for i, (lead, importance) in enumerate(sorted_leads):
+        print(f"{i+1}. {lead}: {importance:.4f}")
+    
+    # Dynamically determine n_timesteps from the data
+    actual_timesteps = infer_timesteps_from_data(shap_df, lead_names[0])
+    print(f"Using {actual_timesteps} timesteps based on data")
+    
+    # Create time points in percentage (0-100%)
+    time_points_pct = np.linspace(0, 100, actual_timesteps)
+    
+    # Create main plot with idealized ECG at top and all leads below
+    create_main_comprehensive_plot(lead_names, shap_df, time_points_pct, output_dir, lead_importances, actual_timesteps)
+    
+    return lead_names, lead_importances
+
+def infer_timesteps_from_data(shap_df, sample_lead):
+    """Infer number of timesteps from the data"""
+    lead_columns = [col for col in shap_df.columns if col.startswith(sample_lead + '_t')]
+    return len(lead_columns)
+
+def create_main_comprehensive_plot(lead_names, shap_df, time_points_pct, output_dir, lead_importances, n_timesteps):
+    """
+    Create main comprehensive plot with idealized ECG at top and all leads' SHAP values below
+    """
+    n_leads = len(lead_names)
+    
+    # Create larger figure with better spacing
+    fig = plt.figure(figsize=(18, 4 + n_leads * 1.5))  # Increased size
+    
+    # Create grid with more space for idealized ECG and better spacing for leads
+    gs = GridSpec(n_leads + 1, 1, figure=fig, height_ratios=[3] + [1.2] * n_leads, hspace=0.4)
+    
+    # Plot 1: Idealized ECG at the top
+    ax_ideal = fig.add_subplot(gs[0])
+    plot_idealized_heartbeat(ax_ideal, time_points_pct)
+    
+    # Plot 2: SHAP values for each lead below (in order of importance)
+    for i, lead in enumerate(lead_names):
+        ax_shap = fig.add_subplot(gs[i + 1])
+        
+        # Get SHAP values for this lead
+        lead_shap_values = extract_lead_shap_values(shap_df, lead, n_timesteps)
+        
+        # Plot SHAP bars for this lead with cardiac cycle regions
+        plot_lead_shap_bars(ax_shap, time_points_pct, lead_shap_values, lead, lead_importances[i], i + 1, n_leads)
+    
+    plt.tight_layout()
+    plt.savefig(f'{output_dir}/idealized_with_all_leads_shap.png', 
+                dpi=300, bbox_inches='tight')
+    plt.savefig(f'{output_dir}/idealized_with_all_leads_shap.pdf', 
+                bbox_inches='tight')
+    plt.close()
+    
+    print(f"Main comprehensive plot saved: {output_dir}/idealized_with_all_leads_shap.png")
+
+def extract_lead_shap_values(shap_df, lead, n_timesteps):
+    """
+    Extract SHAP values for a specific lead across all timesteps
+    """
+    lead_shap_values = np.zeros(n_timesteps)
+    found_features = 0
+    
+    for t in range(n_timesteps):
+        feature_name = f"{lead}_t{t}"
+        if feature_name in shap_df.columns:
+            # Take mean absolute SHAP value for this time point
+            lead_shap_values[t] = np.abs(shap_df[feature_name]).mean()
+            found_features += 1
+        else:
+            print(f"Warning: Feature {feature_name} not found in SHAP DataFrame")
+            lead_shap_values[t] = 0
+    
+    print(f"Lead {lead} - Found {found_features}/{n_timesteps} features, SHAP range: [{lead_shap_values.min():.6f}, {lead_shap_values.max():.6f}]")
+    return lead_shap_values
+
+def plot_idealized_heartbeat(ax, time_points_pct):
+    """
+    Plot accurate idealized ECG heartbeat on given axis with proper waveform characteristics
+    """
+    # Create time array
+    t = time_points_pct
+    
+    # Initialize ECG signal
+    ideal_ecg = np.zeros_like(t, dtype=float)
+    
+    # Define standard ECG regions with PR segment removed and P-wave extended
+    ecg_regions = [
+        {'name': 'P Wave', 'description': 'Atrial Depolarization', 'start': 0, 'end': 25, 'center': 12.5, 'color': 'lightblue', 'alpha': 0.5},
+        {'name': 'QRS Complex', 'description': 'Ventricular Depolarization', 'start': 25, 'end': 42, 'center': 33.5, 'color': 'lightgreen', 'alpha': 0.5},
+        {'name': 'ST Segment', 'description': 'Ventricular Plateau', 'start': 42, 'end': 58, 'center': 50, 'color': 'orange', 'alpha': 0.5},
+        {'name': 'T Wave', 'description': 'Ventricular Repolarization', 'start': 58, 'end': 83, 'center': 70.5, 'color': 'violet', 'alpha': 0.5},
+        {'name': 'TP Segment', 'description': 'Ventricular Diastole', 'start': 83, 'end': 100, 'center': 91.5, 'color': 'saddlebrown', 'alpha': 0.5}
+    ]
+    
+    # Define precise landmark positions and values FIRST to ensure consistency
+    # We align the signal generation to these landmarks
+    landmark_data = [
+        # P wave peak - Center of P region
+        {'time': 12.5, 'label': 'P', 'value': 0.25, 'description': 'Atrial Depolarization'},
+        # Q wave nadir - Early QRS
+        {'time': 27.0, 'label': 'Q', 'value': -0.15, 'description': 'Start of Ventricular Depolarization'},
+        # R wave peak - Center of QRS
+        {'time': 33.5, 'label': 'R', 'value': 1.0, 'description': 'Ventricular Depolarization Peak'},
+        # S wave nadir - Late QRS
+        {'time': 39.0, 'label': 'S', 'value': -0.25, 'description': 'End of Ventricular Depolarization'},
+        # T wave peak - Center of T region
+        {'time': 70.5, 'label': 'T', 'value': 0.35, 'description': 'Ventricular Repolarization'}
+    ]
+    
+    # Helper to create Gaussian waves
+    def gaussian(x, mu, sig, amp):
+        return amp * np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
+
+    # Create accurate ECG waveform components using Gaussians for natural shape
+    for i, time_val in enumerate(t):
+        val = 0.0
+        
+        # P Wave: centered at 12.5, wide
+        val += gaussian(time_val, 12.5, 4.0, 0.25)
+        
+        # Q Wave: centered at 27.0, narrow
+        val += gaussian(time_val, 27.0, 1.2, -0.15)
+        
+        # R Wave: centered at 33.5, narrow and tall
+        val += gaussian(time_val, 33.5, 1.5, 1.0)
+        
+        # S Wave: centered at 39.0, narrow
+        val += gaussian(time_val, 39.0, 1.2, -0.25)
+        
+        # T Wave: centered at 70.5, wide
+        val += gaussian(time_val, 70.5, 6.0, 0.35)
+        
+        # Assign to array
+        ideal_ecg[i] = val
+    
+    # Plot cardiac cycle regions with increased alpha
+    for region in ecg_regions:
+        ax.axvspan(region['start'], region['end'], alpha=region['alpha'], color=region['color'], zorder=0)
+        
+        # Add region labels - name above and description below
+        ax.text(region['center'], -0.45, region['name'], 
+                ha='center', va='center', fontsize=9, fontweight='bold', 
+                bbox=dict(boxstyle="round,pad=0.2", facecolor='white', alpha=0.9, edgecolor=region['color']),
+                zorder=10)
+        ax.text(region['center'], -0.65, region['description'], 
+                ha='center', va='center', fontsize=8, 
+                bbox=dict(boxstyle="round,pad=0.2", facecolor='white', alpha=0.8, edgecolor=region['color']),
+                zorder=10)
+    
+    # Plot idealized ECG
+    ax.plot(t, ideal_ecg, 'k-', linewidth=2, label='Idealized ECG', zorder=5)
+    ax.set_title('Idealized ECG with Standard Regions', fontsize=16, fontweight='bold', pad=20)
+    ax.set_ylabel('Normalized Amplitude', fontsize=12)
+    ax.grid(True, alpha=0.3, zorder=1)
+    
+    # Plot ECG landmarks with annotations
+    for landmark in landmark_data:
+        # Calculate actual value on the curve at this time
+        # This ensures the point is exactly ON the line, even if Gaussian summation shifts it slightly
+        curve_val = 0.0
+        t_lm = landmark['time']
+        curve_val += gaussian(t_lm, 12.5, 4.0, 0.25)
+        curve_val += gaussian(t_lm, 27.0, 1.2, -0.15)
+        curve_val += gaussian(t_lm, 33.5, 1.5, 1.0)
+        curve_val += gaussian(t_lm, 39.0, 1.2, -0.25)
+        curve_val += gaussian(t_lm, 70.5, 6.0, 0.35)
+
+        ax.plot(t_lm, curve_val, 'ro', markersize=8, zorder=10)
+        
+        # Position text based on whether it's a peak or trough
+        if landmark['value'] > 0:
+            text_y_offset = 0.15
+            va = 'bottom'
+        else:
+            text_y_offset = -0.15
+            va = 'top'
+        
+        ax.annotate(landmark['label'], 
+                   (t_lm, curve_val), 
+                   xytext=(0, text_y_offset * 30),
+                   textcoords='offset points',
+                   ha='center', va=va,
+                   fontweight='bold', 
+                   fontsize=12,
+                   bbox=dict(boxstyle="round,pad=0.3", facecolor='white', alpha=0.9, edgecolor='red'),
+                   zorder=15)
+        
+    ax.set_ylim(-0.8, 1.5)  # Increased range to accommodate annotations
+    ax.set_xlim(0, 100)
+    
+    # Remove x-axis labels and ticks for the top plot
+    ax.set_xticklabels([])
+    ax.set_xticks([])
+    ax.set_xlabel('')
+    
+    ax.tick_params(axis='y', which='major', labelsize=10)
+    
+    # Add a legend for the waveform
+    ax.legend(loc='upper right', framealpha=0.9)
+
+def plot_lead_shap_bars(ax, time_points_pct, shap_importance, lead, lead_importance, rank, total_leads):
+    """
+    Plot SHAP importance bars for a single lead with cardiac cycle regions
+    """
+    if len(shap_importance) == len(time_points_pct):
+        # Define cardiac cycle regions with PR segment removed and P-wave extended
+        ecg_regions = [
+            {'name': 'P Wave', 'start': 0, 'end': 25, 'color': 'lightblue', 'alpha': 0.5},
+            {'name': 'QRS Complex', 'start': 25, 'end': 42, 'color': 'lightgreen', 'alpha': 0.5},
+            {'name': 'ST Segment', 'start': 42, 'end': 58, 'color': 'orange', 'alpha': 0.5},
+            {'name': 'T Wave', 'start': 58, 'end': 83, 'color': 'violet', 'alpha': 0.5},
+            {'name': 'TP Segment', 'start': 83, 'end': 100, 'color': 'saddlebrown', 'alpha': 0.5}
+        ]
+        
+        # Plot cardiac cycle regions as background
+        for region in ecg_regions:
+            ax.axvspan(region['start'], region['end'], alpha=region['alpha'], color=region['color'], zorder=0)
+        
+        bar_width = 100 / len(time_points_pct) * 0.8
+        
+        # Debug info
+        max_shap = np.max(shap_importance) if len(shap_importance) > 0 else 0
+        print(f"  {rank}. {lead} - Max SHAP: {max_shap:.6f}")
+        
+        # Create bars with dark grey color
+        if max_shap > 0:
+            bars = ax.bar(time_points_pct, shap_importance, alpha=0.8, color='darkslategray', 
+                         width=bar_width, edgecolor='black', linewidth=0.5, zorder=5)
+            
+            # Highlight top region
+            top_idx = np.argmax(shap_importance)
+            top_value = shap_importance[top_idx]
+            top_time = time_points_pct[top_idx]
+            
+            # Highlight the top bar with red
+            bars[top_idx].set_color('red')
+            bars[top_idx].set_alpha(1.0)
+            bars[top_idx].set_edgecolor('darkred')
+            bars[top_idx].set_linewidth(1)
+            
+            # Add annotation for top region
+            if top_value > max_shap * 0.1:  # Only annotate if significant
+                ax.text(top_time, top_value * 0.9, 
+                       f'{top_time:.1f}%', 
+                       ha='center', va='top', fontsize=8, fontweight='bold', zorder=10,
+                       bbox=dict(boxstyle="round,pad=0.2", facecolor='white', alpha=0.9))
+        else:
+            # If no SHAP values, plot zeros with different color
+            bars = ax.bar(time_points_pct, shap_importance, alpha=0.3, color='gray', 
+                         width=bar_width, edgecolor='darkgray', linewidth=0.5, zorder=5)
+        
+        # Set labels with better formatting
+        ax.set_ylabel(f'{lead}\n({lead_importance:.3f})\nSHAP Value', fontsize=10, rotation=0, 
+                     ha='right', va='center', labelpad=10)
+        ax.grid(True, alpha=0.3, zorder=1)
+        ax.set_xlim(0, 100)
+        
+        # Remove x-axis labels and ticks for all but bottom plot
+        if rank < total_leads:
+            ax.set_xticklabels([])
+            ax.set_xticks([])
+        else:
+            ax.set_xlabel('Cardiac Cycle (%)', fontsize=11)
+            ax.tick_params(axis='x', labelsize=9)
+        
+        # Set y-axis limit based on max SHAP value with some margin
+        max_shap_val = max_shap if max_shap > 0 else 0.001
+        ax.set_ylim(0, max_shap_val * 1.3)
+        ax.tick_params(axis='y', labelsize=8)
+        
+        # Remove y-axis ticks and labels (the numbers 1-12)
+        ax.set_yticklabels([])
+        ax.set_yticks([])
+
+def create_lead_importance_table(lead_names, lead_importances, output_dir):
+    """
+    Create a simple table showing lead importance ranking
+    """
+    # Create importance DataFrame
+    importance_df = pd.DataFrame({
+        'Rank': range(1, len(lead_names) + 1),
+        'Lead': lead_names,
+        'Total_SHAP_Importance': lead_importances
+    })
+    
+    # Save to CSV
+    importance_df.to_csv(f'{output_dir}/lead_importance_ranking.csv', index=False)
+    
+    # Create simple text summary
+    with open(f'{output_dir}/lead_importance_summary.txt', 'w') as f:
+        f.write("Lead Importance Ranking (by Total SHAP Value)\n")
+        f.write("=" * 50 + "\n")
+        for i, (lead, importance) in enumerate(zip(lead_names, lead_importances)):
+            f.write(f"{i+1:2d}. {lead:4s}: {importance:.4f}\n")
+    
+    print(f"Lead importance table saved: {output_dir}/lead_importance_ranking.csv")
+    return importance_df
+
+# Function to fix other plots (avg_with_shap plots)
+def fix_other_plots(analysis_output_dir, n_timesteps=40):
+    """
+    Fix the other plots that are not showing SHAP values
+    """
+    print("\nFixing other plots...")
+    
+    # Load the necessary data
+    shap_values_csv = f'{analysis_output_dir}/shap/shap_values_detailed.csv'
+    top_regions_csv = f'{analysis_output_dir}/shap/top_shap_regions.csv'
+    
+    if not os.path.exists(shap_values_csv) or not os.path.exists(top_regions_csv):
+        print("Required files not found. Skipping other plots.")
+        return
+    
+    shap_df = pd.read_csv(shap_values_csv)
+    top_regions_df = pd.read_csv(top_regions_csv)
+    
+    # Get lead names from top regions
+    lead_names = top_regions_df['feature'].unique()
+    
+    # Create time points
+    time_points_pct = np.linspace(0, 100, n_timesteps)
+    
+    # Recreate avg_with_shap plots for each lead
+    for lead in lead_names:
+        create_fixed_avg_shap_plot(lead, shap_df, time_points_pct, analysis_output_dir, n_timesteps)
+
+def create_fixed_avg_shap_plot(lead, shap_df, time_points_pct, output_dir, n_timesteps):
+    """
+    Create fixed avg_with_shap plot for a single lead
+    """
+    # Extract SHAP values
+    lead_shap_values = extract_lead_shap_values(shap_df, lead, n_timesteps)
+    
+    # Create figure
+    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 8))
+    
+    # Plot 1: Placeholder for average heartbeat (you would need the original data for this)
+    ax1.text(0.5, 0.5, f'Average Heartbeats - {lead}\n(Original data required)', 
+             ha='center', va='center', transform=ax1.transAxes, fontsize=14)
+    ax1.set_title(f'Average Heartbeats - {lead}', fontsize=16)
+    ax1.set_ylabel('Normalized Amplitude')
+    ax1.grid(True, alpha=0.3)
+    
+    # Plot 2: SHAP importance with updated regions
+    # Define cardiac cycle regions with PR segment removed and P-wave extended
+    ecg_regions = [
+        {'name': 'P Wave', 'start': 0, 'end': 25, 'color': 'lightblue', 'alpha': 0.3},
+        {'name': 'QRS Complex', 'start': 25, 'end': 42, 'color': 'lightgreen', 'alpha': 0.3},
+        {'name': 'ST Segment', 'start': 42, 'end': 58, 'color': 'orange', 'alpha': 0.3},
+        {'name': 'T Wave', 'start': 58, 'end': 83, 'color': 'violet', 'alpha': 0.3},
+        {'name': 'TP Segment', 'start': 83, 'end': 100, 'color': 'saddlebrown', 'alpha': 0.3}
+    ]
+    
+    # Plot cardiac cycle regions as background
+    for region in ecg_regions:
+        ax2.axvspan(region['start'], region['end'], alpha=region['alpha'], color=region['color'], zorder=0)
+    
+    bar_width = 100 / len(time_points_pct) * 0.8
+    ax2.bar(time_points_pct, lead_shap_values, alpha=0.7, color='darkslategray', 
+            width=bar_width, label='SHAP Importance', zorder=5)
+    ax2.set_title(f'SHAP Importance - {lead}', fontsize=16)
+    ax2.set_xlabel('Cardiac Cycle (%)')
+    ax2.set_ylabel('SHAP Value')
+    ax2.grid(True, alpha=0.3)
+    ax2.legend()
+    
+    plt.tight_layout()
+    plt.savefig(f'{output_dir}/heartbeat_plots/avg_with_shap_{lead}_FIXED.png', 
+                dpi=300, bbox_inches='tight')
+    plt.close()
+    
+    print(f"Fixed plot saved: {output_dir}/heartbeat_plots/avg_with_shap_{lead}_FIXED.png")
+
+# Additional function to create standalone idealized ECG plot
+def create_standalone_idealized_ecg_plot(output_dir, n_timesteps=40):
+    """
+    Create a standalone idealized ECG plot with the updated regions
+    """
+    time_points_pct = np.linspace(0, 100, n_timesteps)
+    
+    fig, ax = plt.subplots(figsize=(15, 6))
+    plot_idealized_heartbeat(ax, time_points_pct)
+    
+    plt.tight_layout()
+    plt.savefig(f'{output_dir}/standalone_idealized_ecg.png', 
+                dpi=300, bbox_inches='tight')
+    plt.savefig(f'{output_dir}/standalone_idealized_ecg.pdf', 
+                bbox_inches='tight')
+    plt.close()
+    
+    print(f"Standalone idealized ECG plot saved: {output_dir}/standalone_idealized_ecg.png")
+
+# Main execution function
+def main():
+    """
+    Main function to run the SHAP replotting
+    """
+    # Define paths
+    analysis_output_dir = 'results/heartbeat_analysis_500Hz_good'
+    shap_values_csv = f'{analysis_output_dir}/shap/shap_values_detailed.csv'
+    top_regions_csv = f'{analysis_output_dir}/shap/top_shap_regions.csv'
+    output_dir = f'{analysis_output_dir}/shap_replots'
+    
+    # Check if files exist
+    if not os.path.exists(shap_values_csv):
+        print(f"Error: SHAP values file not found at {shap_values_csv}")
+        return
+    if not os.path.exists(top_regions_csv):
+        print(f"Error: Top regions file not found at {top_regions_csv}")
+        return
+    
+    print("Starting SHAP replotting...")
+    print(f"SHAP values file: {shap_values_csv}")
+    print(f"Top regions file: {top_regions_csv}")
+    print(f"Output directory: {output_dir}")
+    
+    # Replot all leads
+    lead_names, lead_importances = replot_idealized_shap_all_leads(
+        shap_values_csv, top_regions_csv, output_dir
+    )
+    
+    # Create lead importance table
+    importance_df = create_lead_importance_table(lead_names, lead_importances, output_dir)
+    
+    # Create standalone idealized ECG plot
+    create_standalone_idealized_ecg_plot(output_dir)
+    
+    # Fix other plots
+    fix_other_plots(analysis_output_dir)
+    
+    print("\n" + "="*60)
+    print("SHAP REPLOTTING COMPLETE!")
+    print("="*60)
+    print(f"Results saved to: {output_dir}")
+    print(f"Lead ranking (by importance):")
+    for i, (lead, importance) in enumerate(zip(lead_names, lead_importances)):
+        print(f"  {i+1:2d}. {lead}: {importance:.4f}")
+
+# Alternative: Function to use in your existing code
+def replot_from_existing_analysis(analysis_output_dir):
+    """
+    Replot from existing analysis output directory
+    """
+    shap_values_csv = f'{analysis_output_dir}/shap/shap_values_detailed.csv'
+    top_regions_csv = f'{analysis_output_dir}/shap/top_shap_regions.csv'
+    output_dir = f'{analysis_output_dir}/shap_replots'
+    
+    return replot_idealized_shap_all_leads(shap_values_csv, top_regions_csv, output_dir)
+
+# Example usage:
+if __name__ == "__main__":
+    main()

archive/classification/run_benchmarking.py

@@ -0,0 +1,64 @@
+import argparse
+import os
+from utils.classifier_benchmarking import perform_benchmark
+
+def run_benchmarking(input_file, output_folder=None):
+    """
+    Run benchmarking on the specified input file.
+    
+    Args:
+        input_file: Path to input CSV
+        output_folder: Directory to save results. If None, derives from input filename.
+    """
+    if not os.path.exists(input_file):
+        print(f"Error: Input file '{input_file}' not found.")
+        return
+
+    # Determine output folder if not specified
+    if output_folder is None:
+        base_name = os.path.splitext(os.path.basename(input_file))[0]
+        output_folder = f'output/benchmark_results/{base_name}'
+
+    print(f"\n{'='*60}")
+    print(f"STARTING BENCHMARKING")
+    print(f"Input:  {input_file}")
+    print(f"Output: {output_folder}")
+    print(f"{'='*60}")
+    
+    try:
+        perform_benchmark(input_file, output_dir=output_folder)
+        print(f"\nBenchmarking complete for {input_file}")
+        print(f"Results saved to: {output_folder}")
+    except Exception as e:
+        print(f"\nError running benchmark: {e}")
+        import traceback
+        traceback.print_exc()
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Run ECG Classifier Benchmarking")
+    parser.add_argument('--input', type=str, help='Path to input CSV file (optional)')
+    parser.add_argument('--all', action='store_true', help='Run on both standard datasets (segmented and full)')
+    
+    args = parser.parse_args()
+    
+    # Default files
+    segmented_file = 'data/segmented_heartbeats.csv'
+    processed_file = 'data/ecg_dataset_processed.csv'
+    
+    # Define explicit output paths for standard runs
+    output_segmented = 'output/benchmark_results/segmented_heartbeats'
+    output_processed = 'output/benchmark_results/ecg_dataset_processed'
+    
+    if args.input:
+        # Run on user specified file (folder derived from filename)
+        run_benchmarking(args.input)
+    elif args.all:
+        # Run on both with specific folders
+        print("Running comprehensive benchmarking on ALL datasets...")
+        run_benchmarking(segmented_file, output_folder=output_segmented)
+        run_benchmarking(processed_file, output_folder=output_processed)
+    else:
+        # Default run
+        print("No input specified. Defaulting to processed ECG data.")
+        print(f"Usage: python run_benchmarking.py --input <path> OR --all")
+        run_benchmarking(processed_file, output_folder=output_processed)

archive/classification/run_classification.py

@@ -0,0 +1,34 @@
+from utils.ecg_classification import perform_classification_analysis
+import os
+# Import the universally defined target frequency
+from run_segmentation import TARGET_FS
+
+if __name__ == "__main__":
+    # Define paths
+    input_csv = 'data/ecg_matrix_dataset_segmented.csv'
+    output_dir = 'output/heartbeat_analysis_matrix'
+    
+    # Check if input exists
+    if not os.path.exists(input_csv):
+        print(f"Error: Input file {input_csv} not found.")
+        print("Please run run_segmentation.py first.")
+        exit(1)
+        
+    print(f"Starting classification analysis...")
+    print(f"Input: {input_csv}")
+    print(f"Output Directory: {output_dir}")
+    print(f"Target Frequency: {TARGET_FS} Hz (imported from run_segmentation.py)")
+    
+    # Run analysis
+    # Normalization happens inside this function
+    results = perform_classification_analysis(
+        segmented_csv=input_csv,
+        output_dir=output_dir,
+        target_fs=TARGET_FS
+    )
+    
+    if results:
+        print("\nAnalysis successful!")
+        print(f"Results available in {output_dir}")
+    else:
+        print("\nAnalysis failed.")

archive/classification/run_inference.py

@@ -0,0 +1,297 @@
+"""
+Run inference on new ECG data using a previously saved model.
+
+Loads the serialised model and its metadata, applies identical preprocessing
+(per-subject normalisation, reshape to numpy3D), and outputs per-subject
+predictions with confidence scores.
+
+Usage:
+    # MI vs Normal
+    python run_inference.py --model mi_vs_normal --input data/new_ecg.csv
+
+    # OMI vs non-OMI
+    python run_inference.py --model omi_vs_nonomi --input data/new_ecg.csv
+
+    # Custom output path
+    python run_inference.py --model mi_vs_normal --input data/new_ecg.csv --output results/preds.csv
+"""
+
+import os
+import argparse
+import pickle
+import json
+import numpy as np
+import pandas as pd
+from tqdm.auto import tqdm
+
+# ─── Configuration ───────────────────────────────────────────────────────────
+
+MODELS_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '../../models/classifier_models'))
+
+# ─── Model Loading ───────────────────────────────────────────────────────────
+
+
+def load_model(use_case):
+    """
+    Load a saved model and its metadata from the models directory.
+
+    Args:
+        use_case: One of 'mi_vs_normal' or 'omi_vs_nonomi'.
+
+    Returns:
+        (model, metadata) tuple.
+    """
+    model_dir = os.path.join(MODELS_DIR, use_case)
+
+    # Metadata
+    metadata_path = os.path.join(model_dir, 'model_metadata.json')
+    if not os.path.exists(metadata_path):
+        raise FileNotFoundError(
+            f"Model metadata not found: {metadata_path}\n"
+            f"Have you run train_and_save_models.py first?")
+
+    with open(metadata_path, 'r') as f:
+        metadata = json.load(f)
+
+    # Model pickle
+    model_path = os.path.join(model_dir, metadata['model_file'])
+    if not os.path.exists(model_path):
+        raise FileNotFoundError(
+            f"Model file not found: {model_path}\n"
+            f"Have you run train_and_save_models.py first?")
+
+    print(f"Loading {metadata['model_name']} model for {metadata['use_case']} ...")
+    with open(model_path, 'rb') as f:
+        model = pickle.load(f)
+
+    print(f"  Model loaded successfully.")
+    print(f"  Trained on:      {metadata['training_data']}")
+    print(f"  Training date:   {metadata['training_date']}")
+    print(f"  Expected input:  {metadata['n_features']} leads x "
+          f"{metadata['n_timesteps']} timesteps")
+    print(f"  Performance metrics on held-out test set:")
+    metrics = metadata.get('test_metrics', {})
+    print(f"    - Accuracy:    {metrics.get('accuracy', 0.0):.4f}")
+    print(f"    - Precision:   {metrics.get('precision', 0.0):.4f}")
+    print(f"    - Recall:      {metrics.get('recall', 0.0):.4f}")
+    print(f"    - F1 Score:    {metrics.get('f1', 0.0):.4f}")
+    print(f"    - Sensitivity: {metrics.get('sensitivity', 0.0):.4f}")
+    print(f"    - Specificity: {metrics.get('specificity', 0.0):.4f}")
+
+    return model, metadata
+
+
+# ─── Preprocessing ───────────────────────────────────────────────────────────
+
+
+def preprocess_data(input_csv, metadata):
+    """
+    Preprocess new ECG data for inference using the same pipeline as training.
+
+    Steps (mirrors prepare_data() in classifier_benchmarking.py):
+      1. Per-subject normalisation (divide each lead by its max absolute value)
+      2. Reshape to numpy3D  (n_subjects, n_features, n_timesteps)
+      3. Truncate or zero-pad each subject to the training n_timesteps
+
+    Args:
+        input_csv:  Path to the input CSV.  Must contain columns
+                    'subject_id', 'timestamp', and the 12 ECG lead columns.
+        metadata:   The model_metadata dict loaded from JSON.
+
+    Returns:
+        (X, valid_subject_ids)  where X has shape
+        (n_subjects, n_features, n_timesteps).
+    """
+    features    = metadata['features']
+    n_timesteps = metadata['n_timesteps']
+    n_features  = metadata['n_features']
+
+    # ── Load ──────────────────────────────────────────────────────────────
+    print(f"\nLoading input data from {input_csv} ...")
+    df = pd.read_csv(input_csv, index_col=['subject_id', 'timestamp'])
+
+    missing = [f for f in features if f not in df.columns]
+    if missing:
+        raise ValueError(f"Missing required lead columns in input: {missing}")
+
+    # ── Normalise ─────────────────────────────────────────────────────────
+    def safe_normalize(x):
+        max_val = x.abs().max()
+        return x / max_val if max_val > 0 else x
+
+    print("Normalising features per subject ...")
+    df_X = df[features].copy()
+    for col in tqdm(df_X.columns, desc="Normalising"):
+        df_X[col] = df_X.groupby(level='subject_id')[col].transform(safe_normalize)
+
+    # ── Reshape ───────────────────────────────────────────────────────────
+    subject_ids = df_X.index.get_level_values('subject_id').unique()
+    timesteps_per_subject = df_X.groupby(level='subject_id').size()
+
+    print(f"  Subjects found:          {len(subject_ids)}")
+    print(f"  Timesteps per subject — "
+          f"min: {timesteps_per_subject.min()}, "
+          f"max: {timesteps_per_subject.max()}")
+    print(f"  Target timesteps (train): {n_timesteps}")
+
+    X_list = []
+    valid_subject_ids = []
+    skipped = 0
+
+    for subject_id in tqdm(subject_ids, desc="Reshaping subjects"):
+        try:
+            subject_data = df_X.xs(subject_id, level='subject_id')
+            actual_len = len(subject_data)
+
+            if actual_len >= n_timesteps:
+                # Truncate to training length
+                arr = subject_data.iloc[:n_timesteps].values.T
+            else:
+                # Zero-pad short subjects
+                arr = np.zeros((n_features, n_timesteps))
+                arr[:, :actual_len] = subject_data.values.T
+                print(f"  [warn] Subject {subject_id}: {actual_len} timesteps "
+                      f"→ padded to {n_timesteps}")
+
+            if arr.shape == (n_features, n_timesteps):
+                X_list.append(arr)
+                valid_subject_ids.append(subject_id)
+            else:
+                skipped += 1
+        except Exception as e:
+            print(f"  [warn] Skipped subject {subject_id}: {e}")
+            skipped += 1
+
+    if not X_list:
+        raise ValueError("No valid subjects found in the input data!")
+
+    X = np.array(X_list)
+    print(f"\n  Preprocessed shape: {X.shape}  "
+          f"({len(valid_subject_ids)} subjects, {skipped} skipped)")
+
+    return X, valid_subject_ids
+
+
+# ─── Inference ───────────────────────────────────────────────────────────────
+
+
+def run_predictions(model, X, metadata, subject_ids):
+    """
+    Run model predictions and package them into a DataFrame.
+
+    Returns:
+        pd.DataFrame with columns:
+            subject_id, predicted_class, confidence,
+            probability_<class0>, probability_<class1>
+    """
+    print(f"\nRunning inference on {len(subject_ids)} subjects ...")
+
+    # Predict labels and probabilities in threading backend to bypass Numba read-only array issue
+    from joblib import parallel_backend
+    with parallel_backend('threading', n_jobs=-1):
+        y_pred = model.predict(X)
+
+        # Predict probabilities (Arsenal & Rocket both support this)
+        has_proba = False
+        y_proba = None
+        try:
+            y_proba = model.predict_proba(X)
+            has_proba = True
+        except (AttributeError, NotImplementedError):
+            pass
+
+    # ── Build output DataFrame ────────────────────────────────────────────
+    results = {
+        'subject_id':      subject_ids,
+        'predicted_class':  y_pred,
+    }
+
+    if has_proba:
+        class_labels = list(model.classes_)
+        for i, cls in enumerate(class_labels):
+            results[f'probability_{cls}'] = y_proba[:, i]
+
+        # Confidence = probability of the predicted class
+        confidence = []
+        for j in range(len(y_pred)):
+            pred_idx = class_labels.index(y_pred[j])
+            confidence.append(float(y_proba[j, pred_idx]))
+        results['confidence'] = confidence
+
+    return pd.DataFrame(results)
+
+
+# ─── Entry Point ──────────────────────────────────────────────────────────────
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description='Run ECG classification inference with a saved model',
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog="""
+Examples:
+  python run_inference.py --model mi_vs_normal_processed --input data/ptb_xl/ecg_dataset_processed.csv
+  python run_inference.py --model mi_vs_normal_segmented --input data/ptb_xl/segmented_heartbeats.csv
+  python run_inference.py --model omi_vs_nonomi          --input data/ecg_matrix_omi_segmented_50_150_90.csv
+  python run_inference.py --model ecg_surgery            --input data/ecg_surgery_segmented_50_150_70.csv
+""")
+
+    parser.add_argument('--model', type=str, required=True,
+                        choices=['mi_vs_normal_processed', 'mi_vs_normal_segmented', 'omi_vs_nonomi', 'ecg_surgery'],
+                        help='Which saved model to use')
+    parser.add_argument('--input', type=str, required=True,
+                        help='Path to input CSV (subject_id, timestamp, 12 leads)')
+    parser.add_argument('--output', type=str, default=None,
+                        help='Path to save predictions CSV '
+                             '(default: results/<model>_predictions.csv)')
+
+    args = parser.parse_args()
+
+    # Default output
+    if args.output is None:
+        os.makedirs('results', exist_ok=True)
+        args.output = f'results/{args.model}_predictions.csv'
+
+    print('=' * 60)
+    print('  ECG CLASSIFICATION - INFERENCE')
+    print('=' * 60)
+
+    # ── Load model ────────────────────────────────────────────────────────
+    model, metadata = load_model(args.model)
+
+    # ── Preprocess ────────────────────────────────────────────────────────
+    X, subject_ids = preprocess_data(args.input, metadata)
+
+    # ── Predict ───────────────────────────────────────────────────────────
+    results_df = run_predictions(model, X, metadata, subject_ids)
+
+    # ── Save ──────────────────────────────────────────────────────────────
+    out_dir = os.path.dirname(args.output)
+    if out_dir:
+        os.makedirs(out_dir, exist_ok=True)
+    results_df.to_csv(args.output, index=False)
+
+    # ── Summary ───────────────────────────────────────────────────────────
+    print(f"\n{'=' * 60}")
+    print(f"  INFERENCE RESULTS")
+    print(f"{'=' * 60}")
+    print(f"  Model:          {metadata['model_name']} ({metadata['use_case']})")
+    print(f"  Total subjects: {len(results_df)}")
+    print(f"\n  Predicted class distribution:")
+    for cls, count in results_df['predicted_class'].value_counts().items():
+        print(f"    {cls}: {count}")
+
+    if 'confidence' in results_df.columns:
+        conf = results_df['confidence']
+        print(f"\n  Confidence statistics:")
+        print(f"    Mean:   {conf.mean():.4f}")
+        print(f"    Median: {conf.median():.4f}")
+        print(f"    Min:    {conf.min():.4f}")
+        print(f"    Max:    {conf.max():.4f}")
+
+    print(f"\n  Results saved to: {args.output}")
+    print(f"{'=' * 60}")
+
+
+if __name__ == '__main__':
+    main()

archive/classification/run_lead_importance_test.py

@@ -0,0 +1,128 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import os
+from sklearn.model_selection import train_test_split
+from sklearn.neural_network import MLPClassifier
+from sklearn.preprocessing import StandardScaler
+from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
+from tqdm.auto import tqdm
+
+# Import common configuration
+try:
+    from run_segmentation import TARGET_FS
+except ImportError:
+    TARGET_FS = 500  # Default fallback
+
+def run_lead_importance_test(input_csv, output_dir):
+    """
+    Test the accuracy of identifying MI using each ECG lead individually.
+    """
+    os.makedirs(output_dir, exist_ok=True)
+    
+    # Load segmented data
+    print(f"Loading data from {input_csv}...")
+    df_segmented = pd.read_csv(input_csv, index_col=['subject_id', 'timestamp'])
+    
+    if df_segmented.empty:
+        print("Error: Input data is empty.")
+        return
+    
+    # Define features/leads
+    leads = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+    label_map = {'Post-Procedural MI': 0, 'Pre-Procedural MI': 1}
+    
+    results = []
+    
+    print(f"\nEvaluating {len(leads)} leads individually at {TARGET_FS} Hz...")
+    
+    for lead in leads:
+        print(f"\n>>> Testing Lead {lead}...")
+        
+        # 1. Prepare data for this lead
+        df_lead = df_segmented[[lead, 'class']].copy()
+        
+        # Normalize per heartbeat
+        df_lead[lead] = df_lead.groupby(level='subject_id')[lead].transform(
+            lambda x: x / x.abs().max() if x.abs().max() > 0 else x
+        )
+        
+        # Reshape to (heartbeats, timesteps)
+        beat_ids = df_lead.index.get_level_values('subject_id').unique()
+        timesteps_per_beat = df_lead.groupby(level='subject_id').size().min()
+        
+        X_list = []
+        y_list = []
+        
+        for beat_id in beat_ids:
+            beat_data = df_lead.xs(beat_id, level='subject_id')
+            if len(beat_data) >= timesteps_per_beat:
+                X_list.append(beat_data[lead].values[:timesteps_per_beat])
+                y_list.append(label_map[beat_data['class'].iloc[0]])
+        
+        X = np.array(X_list)
+        y = np.array(y_list)
+        
+        # 2. Train and Evaluate
+        X_train, X_test, y_train, y_test = train_test_split(
+            X, y, test_size=0.25, random_state=42, stratify=y
+        )
+        
+        scaler = StandardScaler()
+        X_train_scaled = scaler.fit_transform(X_train)
+        X_test_scaled = scaler.transform(X_test)
+        
+        clf = MLPClassifier(
+            hidden_layer_sizes=(100,),
+            max_iter=1000,
+            random_state=42
+        )
+        clf.fit(X_train_scaled, y_train)
+        
+        y_pred = clf.predict(X_test_scaled)
+        
+        # Collect metrics
+        metrics = {
+            'Lead': lead,
+            'Accuracy': accuracy_score(y_test, y_pred),
+            'Precision': precision_score(y_test, y_pred, zero_division=0),
+            'Recall': recall_score(y_test, y_pred, zero_division=0),
+            'F1_Score': f1_score(y_test, y_pred, zero_division=0)
+        }
+        results.append(metrics)
+        print(f"  Accuracy: {metrics['Accuracy']:.4f}")
+    
+    # 3. Save and Visualize Results
+    results_df = pd.DataFrame(results).sort_values(by='Accuracy', ascending=False)
+    results_df.to_csv(f'{output_dir}/lead_accuracies.csv', index=False)
+    
+    print("\n" + "="*40)
+    print("LEAD IMPORTANCE TEST RESULTS")
+    print("="*40)
+    print(results_df.to_string(index=False))
+    
+    # Plot results
+    plt.figure(figsize=(12, 6))
+    plt.bar(results_df['Lead'], results_df['Accuracy'], color='skyblue', edgecolor='navy')
+    plt.axhline(y=0.5, color='red', linestyle='--', label='Baseline (Random)')
+    plt.title(f'MLP Accuracy per Single ECG Lead ({TARGET_FS} Hz)', fontsize=15, fontweight='bold')
+    plt.xlabel('ECG Lead', fontsize=12)
+    plt.ylabel('Accuracy', fontsize=12)
+    plt.ylim(0, 1.0)
+    plt.grid(axis='y', alpha=0.3)
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(f'{output_dir}/lead_accuracies_plot.png', dpi=300)
+    plt.close()
+    
+    print(f"\nResults saved to {output_dir}")
+
+if __name__ == "__main__":
+    #input_csv = 'data/segmented_heartbeats.csv'
+    input_csv = 'data/ecg_matrix_dataset_segmented.csv'
+    output_dir = 'output/lead_importance_test'
+    
+    if not os.path.exists(input_csv):
+        print(f"Error: {input_csv} not found. Run segmentation first.")
+    else:
+        run_lead_importance_test(input_csv, output_dir)

archive/classification/train_and_save_models.py

@@ -0,0 +1,427 @@
+import os
+import argparse
+import pickle
+import json
+import time
+import numpy as np
+import pandas as pd
+from sklearn.model_selection import train_test_split
+from imblearn.combine import SMOTETomek
+from sklearn.metrics import (accuracy_score, precision_score,
+                             recall_score, f1_score, confusion_matrix)
+from sktime.classification.kernel_based import Arsenal, RocketClassifier
+from sktime.classification.deep_learning.inceptiontime import InceptionTimeClassifier
+from tqdm.auto import tqdm
+import warnings
+warnings.filterwarnings('ignore')
+
+
+# ─── Configuration ───────────────────────────────────────────────────────────
+
+MODELS_DIR = 'models'
+
+FEATURES = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5',
+            'III', 'aVF', 'V3', 'V6']
+
+USE_CASES = {
+    'mi_vs_normal_processed': {
+        'name': 'MI vs Normal (Processed)',
+        'input_csv': 'data/ptb_xl/ecg_dataset_processed.csv',
+        'model_class': Arsenal,
+        'model_name': 'Arsenal',
+        'model_file': 'arsenal_processed_model.pkl',
+        'positive_class': 'MYOCARDIAL_INFARCTION',
+        'negative_class': 'NORMAL',
+        'hyperparameters': {
+            'num_kernels': 10000,
+            'random_state': 42,
+            'n_jobs': 1,
+        },
+    },
+    'mi_vs_normal_segmented': {
+        'name': 'MI vs Normal (Segmented)',
+        'input_csv': 'data/ptb_xl/segmented_heartbeats.csv',
+        'model_class': Arsenal,
+        'model_name': 'Arsenal',
+        'model_file': 'arsenal_segmented_model.pkl',
+        'positive_class': 'MYOCARDIAL_INFARCTION',
+        'negative_class': 'NORMAL',
+        'hyperparameters': {
+            'num_kernels': 20000,
+            'random_state': 42,
+            'n_jobs': 1,
+        },
+    },
+    'omi_vs_nonomi': {
+        'name': 'OMI vs non-OMI',
+        'input_csv': 'data/ecg_matrix_omi_segmented_50_150_90.csv',
+        'model_class': RocketClassifier,
+        'model_name': 'Rocket',
+        'model_file': 'rocket_model.pkl',
+        'positive_class': 'OMI',
+        'negative_class': 'non-OMI',
+        'hyperparameters': {
+            'num_kernels': 50000,
+            'random_state': 42,
+            'n_jobs': -1,
+        },
+    },
+    'ecg_surgery': {
+        'name': 'Pre vs Post Procedural MI (Surgery)',
+        'input_csv': 'data/ecg_surgery_segmented_50_150_70.csv',
+        'model_class': InceptionTimeClassifier,
+        'model_name': 'InceptionTime',
+        'model_file': 'inceptiontime_model.pkl',
+        'positive_class': 'pre-procedural MI',
+        'negative_class': 'post-procedural MI',
+        'hyperparameters': {
+            'n_epochs': 200,
+            'batch_size': 64,
+            'random_state': 42,
+            'verbose': True,
+        },
+    },
+}
+
+# ─── Data Preprocessing ──────────────────────────────────────────────────────
+
+def prepare_data(df):
+    """Prepare data for benchmarking in sktime compatible format (numpy3D)"""
+    features = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+    
+    def safe_normalize(x):
+        max_val = x.abs().max()
+        # Prevent division by zero which yields NaNs
+        return x / max_val if max_val > 0 else x
+
+    print(" Normalizing features...")
+    # Normalize features within each subject
+    df_X = df[features].copy()
+    for col in tqdm(df_X.columns, desc="Normalizing features"):
+        df_X[col] = df_X.groupby(level='subject_id')[col].transform(safe_normalize)
+    
+    # Get unique subjects
+    subject_ids = df_X.index.get_level_values('subject_id').unique()
+    
+    # Check timesteps per subject
+    timesteps_per_subject = df_X.groupby(level='subject_id').size()
+    print(f" Timesteps per subject stats:")
+    print(f"   Min: {timesteps_per_subject.min()}, Max: {timesteps_per_subject.max()}, Mean: {timesteps_per_subject.mean():.1f}")
+    
+    # Use minimum timesteps to ensure consistent shape
+    n_timesteps = timesteps_per_subject.min()
+    n_features = len(features)
+    
+    print(f" Using {n_timesteps} timesteps per subject (minimum)")
+    
+    # Manually reshape to numpy3D format: (n_instances, n_features, n_timesteps)
+    print(" Reshaping data to numpy3D format...")
+    X_list = []
+    y_list = []
+    valid_subject_ids = []
+    
+    for subject_id in tqdm(subject_ids, desc="Processing subjects"):
+        try:
+            subject_data = df_X.xs(subject_id, level='subject_id')
+            
+            # Ensure we have enough timesteps
+            if len(subject_data) >= n_timesteps:
+                # Take first n_timesteps and transpose to (n_features, n_timesteps)
+                subject_array = subject_data.iloc[:n_timesteps].values.T
+                
+                # Check shape
+                if subject_array.shape == (n_features, n_timesteps):
+                    X_list.append(subject_array)
+                    
+                    # Get label for this subject
+                    subject_label = df.loc[subject_id].iloc[0]['class']
+                    y_list.append(subject_label)
+                    valid_subject_ids.append(subject_id)
+            else:
+                continue
+                
+        except Exception as e:
+            continue
+    
+    # Convert to numpy arrays
+    X = np.array(X_list)  # Shape: (n_instances, n_features, n_timesteps)
+    y = pd.Series(y_list, index=valid_subject_ids)
+    
+    print(f" Final data shape: {X.shape}")
+    print(f" Class distribution:")
+    print(y.value_counts())
+    
+    # Split into train/test
+    print(" Splitting data into train/test sets...")
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.25, random_state=42, stratify=y
+    )
+    
+    # --- SMOTE-TOMEK INJECTION ---
+    print(" Executing SMOTE-Tomek to isolate and synthetically balance training geometries...")
+    try:
+        # Flatten 3D arrays to 2D for SMOTETomek: (N, features, timesteps) -> (N, features * timesteps)
+        n_train, n_feat, n_time = X_train.shape
+        X_train_flat = X_train.reshape(n_train, n_feat * n_time)
+        
+        # Balance classes
+        smote_tomek = SMOTETomek(random_state=42)
+        X_resampled_flat, y_resampled = smote_tomek.fit_resample(X_train_flat, y_train)
+        
+        # Deflate back out to precise 3D format for Sktime models
+        n_new = X_resampled_flat.shape[0]
+        X_train = X_resampled_flat.reshape(n_new, n_feat, n_time)
+        y_train = y_resampled
+        
+        print(f" -> SMOTE Success! Training geometries augmented from {n_train} to {n_new} samples.")
+    except Exception as e:
+        print(f" [!] SMOTETomek failed (are you missing imbalanced-learn?): {e}")
+    # -----------------------------
+    
+    print(f" Final Training set: {X_train.shape}")
+    print(f" Final Test set: {X_test.shape}")
+    print(f" y_train balanced distribution: {y_train.value_counts().to_dict()}")
+    print(f" y_test strictly isolated distribution: {y_test.value_counts().to_dict()}")
+    
+    return X_train, X_test, y_train, y_test
+
+
+# ─── Training Logic ──────────────────────────────────────────────────────────
+
+
+def train_and_save(use_case_key, config):
+    """
+    Train a single model for the given use case, evaluate it, and persist
+    both the serialised model and a metadata JSON to disk.
+
+    Returns True on success, False on failure.
+    """
+
+    print(f"\n{'=' * 60}")
+    print(f"  TRAINING: {config['name']}")
+    print(f"  Model:    {config['model_name']}")
+    print(f"{'=' * 60}")
+
+    input_csv = config['input_csv']
+    if not os.path.exists(input_csv):
+        print(f"[ERROR] Input file not found: {input_csv}")
+        return False
+
+    # ── Load data ──────────────────────────────────────────────────────────
+    print(f"\nLoading data from {input_csv} ...")
+    df = pd.read_csv(input_csv, index_col=['subject_id', 'timestamp'])
+    
+    # Map labels to Pre/Post Procedural MI if it's the surgery dataset
+    if use_case_key == 'ecg_surgery':
+        print("  Mapping class labels to 'pre-procedural MI' and 'post-procedural MI' ...")
+        df['class'] = df['class'].replace({
+            'Pre-Surgery MI': 'pre-procedural MI',
+            'Post-Surgery MI': 'post-procedural MI'
+        })
+        
+    print(f"  DataFrame shape: {df.shape}")
+
+    subject_class = df.groupby(level='subject_id').first()['class']
+    print(f"  Class distribution (subjects):\n{subject_class.value_counts().to_string()}")
+
+    # ── Prepare data (normalize → reshape → split → SMOTE-Tomek) ──────────
+    print("\nPreparing data (normalize, reshape, split, SMOTE-Tomek) ...")
+    X_train, X_test, y_train, y_test = prepare_data(df)
+
+    n_timesteps = X_train.shape[2]
+    n_features = X_train.shape[1]
+    print(f"  n_timesteps = {n_timesteps}")
+    print(f"  n_features  = {n_features}")
+    print(f"  X_train shape: {X_train.shape}  |  X_test shape: {X_test.shape}")
+
+    # ── Initialise & train model ──────────────────────────────────────────
+    model = config['model_class'](**config['hyperparameters'])
+
+    params_str = ""
+    if 'num_kernels' in config['hyperparameters']:
+        params_str = f" (num_kernels={config['hyperparameters']['num_kernels']})"
+    elif 'n_epochs' in config['hyperparameters']:
+        params_str = f" (n_epochs={config['hyperparameters']['n_epochs']})"
+        if 'batch_size' in config['hyperparameters']:
+            params_str += f" (batch_size={config['hyperparameters']['batch_size']})"
+
+    print(f"\nTraining {config['model_name']}{params_str} ...")
+    train_start = time.time()
+    
+    from joblib import parallel_backend
+    with parallel_backend('threading'):
+        model.fit(X_train, y_train)
+        train_time = time.time() - train_start
+        print(f"  Training completed in {train_time:.1f}s")
+
+        # ── Evaluate on test set ──────────────────────────────────────────────
+        print("\nEvaluating on held-out test set ...")
+        y_pred = model.predict(X_test)
+
+    positive = config['positive_class']
+    y_test_num = (y_test == positive).astype(int)
+    y_pred_num = (y_pred == positive).astype(int)
+
+    accuracy    = accuracy_score(y_test_num, y_pred_num)
+    precision   = precision_score(y_test_num, y_pred_num, zero_division=0)
+    recall      = recall_score(y_test_num, y_pred_num, zero_division=0)
+    f1          = f1_score(y_test_num, y_pred_num, zero_division=0)
+
+    cm = confusion_matrix(y_test_num, y_pred_num)
+    if cm.shape == (2, 2):
+        tn, fp, fn, tp = cm.flatten()
+        sensitivity = tp / (tp + fn) if (tp + fn) > 0 else 0
+        specificity = tn / (tn + fp) if (tn + fp) > 0 else 0
+    else:
+        sensitivity = specificity = 0.0
+
+    test_metrics = {
+        'accuracy':         float(accuracy),
+        'precision':        float(precision),
+        'recall':           float(recall),
+        'f1':               float(f1),
+        'sensitivity':      float(sensitivity),
+        'specificity':      float(specificity),
+        'confusion_matrix': cm.tolist(),
+    }
+
+    print(f"\n{'-' * 40}")
+    print(f"  TEST-SET RESULTS")
+    print(f"{'-' * 40}")
+    print(f"  Accuracy:    {accuracy:.4f}")
+    print(f"  Precision:   {precision:.4f}")
+    print(f"  Recall:      {recall:.4f}")
+    print(f"  F1 Score:    {f1:.4f}")
+    print(f"  Sensitivity: {sensitivity:.4f}")
+    print(f"  Specificity: {specificity:.4f}")
+    print(f"  Confusion Matrix:\n    {cm}")
+
+    # ── Save model (pickle) ──────────────────────────────────────────────
+    output_dir = os.path.join(MODELS_DIR, use_case_key)
+    os.makedirs(output_dir, exist_ok=True)
+
+    model_path = os.path.join(output_dir, config['model_file'])
+    print(f"\nSaving model to {model_path} ...")
+    with open(model_path, 'wb') as f:
+        pickle.dump(model, f)
+
+    model_size_mb = os.path.getsize(model_path) / (1024 * 1024)
+    print(f"  Model saved  ({model_size_mb:.1f} MB)")
+
+    # ── Save metadata (JSON) ─────────────────────────────────────────────
+    metadata = {
+        'model_name':       config['model_name'],
+        'model_file':       config['model_file'],
+        'use_case':         config['name'],
+        'positive_class':   config['positive_class'],
+        'negative_class':   config['negative_class'],
+        'class_labels':     [config['negative_class'], config['positive_class']],
+        'features':         FEATURES,
+        'n_features':       int(n_features),
+        'n_timesteps':      int(n_timesteps),
+        'hyperparameters':  config['hyperparameters'],
+        'training_data':    input_csv,
+        'training_date':    time.strftime('%Y-%m-%d %H:%M:%S'),
+        'training_time_seconds': round(train_time, 2),
+        'preprocessing': {
+            'normalization': 'per_subject_max_abs',
+            'smote_tomek':   True,
+            'test_size':     0.25,
+            'random_state':  42,
+        },
+        'training_samples': int(X_train.shape[0]),
+        'test_samples':     int(X_test.shape[0]),
+        'test_metrics':     test_metrics,
+    }
+
+    metadata_path = os.path.join(output_dir, 'model_metadata.json')
+    with open(metadata_path, 'w') as f:
+        json.dump(metadata, f, indent=2)
+    print(f"  Metadata saved to {metadata_path}")
+
+    print(f"\n[OK] {config['name']} - done.")
+    return True
+
+
+# ─── Entry Point ──────────────────────────────────────────────────────────────
+
+def show_saved_metrics():
+    """
+    Search the MODELS_DIR directory for any saved model metadata JSON files
+    and print a beautiful summary of their performance metrics.
+    """
+    print('=' * 60)
+    print('  TRAINED ECG MODEL METRICS SUMMARY')
+    print('=' * 60)
+    
+    found = False
+    if os.path.exists(MODELS_DIR):
+        for item in sorted(os.listdir(MODELS_DIR)):
+            meta_path = os.path.join(MODELS_DIR, item, 'model_metadata.json')
+            if os.path.exists(meta_path):
+                try:
+                    with open(meta_path, 'r') as f:
+                        meta = json.load(f)
+                    
+                    found = True
+                    print(f"\nUse Case:      {meta.get('use_case', 'Unknown')}")
+                    print(f"Model Class:   {meta.get('model_name', 'Unknown')}")
+                    print(f"Model File:    {os.path.join(MODELS_DIR, item, meta.get('model_file', ''))}")
+                    print(f"Trained On:    {meta.get('training_data', 'Unknown')}")
+                    print(f"Training Date: {meta.get('training_date', 'Unknown')}")
+                    if 'training_time_seconds' in meta:
+                        print(f"Train Time:    {meta['training_time_seconds']:.1f}s")
+                    
+                    metrics = meta.get('test_metrics', {})
+                    print("Test Metrics:")
+                    print(f"  - Accuracy:    {metrics.get('accuracy', 0.0):.4f}")
+                    print(f"  - Precision:   {metrics.get('precision', 0.0):.4f}")
+                    print(f"  - Recall:      {metrics.get('recall', 0.0):.4f}")
+                    print(f"  - F1 Score:    {metrics.get('f1', 0.0):.4f}")
+                    print(f"  - Sensitivity: {metrics.get('sensitivity', 0.0):.4f}")
+                    print(f"  - Specificity: {metrics.get('specificity', 0.0):.4f}")
+                    
+                    cm = metrics.get('confusion_matrix', None)
+                    if cm:
+                        print(f"  - Confusion Matrix:\n      {np.array(cm)}")
+                    print("-" * 60)
+                except Exception as e:
+                    print(f"[Error loading metadata from {meta_path}]: {e}")
+                    
+    if not found:
+        print(f"\nNo saved model metadata files found under '{MODELS_DIR}/'.")
+        print("Please train models first to generate performance metrics.")
+    print('=' * 60)
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        description='Train and save the best ECG classification models')
+    parser.add_argument(
+        '--use-case', type=str,
+        choices=['mi_vs_normal_processed', 'mi_vs_normal_segmented', 'omi_vs_nonomi', 'ecg_surgery', 'all'],
+        default='all',
+        help='Which use case to train (default: all)')
+    parser.add_argument(
+        '--show-metrics', action='store_true',
+        help='Display performance metrics of already trained/saved models and exit')
+    args = parser.parse_args()
+
+    if args.show_metrics:
+        show_saved_metrics()
+    else:
+        cases = list(USE_CASES.keys()) if args.use_case == 'all' else [args.use_case]
+
+        print('=' * 60)
+        print('  ECG MODEL TRAINING & SERIALISATION')
+        print('=' * 60)
+
+        for key in cases:
+            success = train_and_save(key, USE_CASES[key])
+            if not success:
+                print(f"\n[FAIL] Could not train {USE_CASES[key]['name']}")
+
+        print(f"\n{'=' * 60}")
+        print(f"  ALL TRAINING COMPLETE")
+        print(f"  Models saved to: {MODELS_DIR}/")
+        print(f"{'=' * 60}")

archive/digitization/README.md

@@ -0,0 +1,21 @@
+# archive/digitization
+
+Legacy ECG image digitization scripts that have been **superseded** by the
+Streamlit dashboard in the workspace root.
+
+## Contents
+
+| File | Status | Reason for archival |
+| --- | --- | --- |
+| `run_ecg.py` | Legacy batch digitizer | Processes a directory tree of ECG paper scans (one CSV per image). Uses Linux-only `SIGALRM` timeouts and hard-coded model paths under `Submission code/...`. Replaced by the interactive pipeline in `backend/digitization_runner.py` invoked from `app.py`. |
+
+## Current digitization entry point
+
+Use the Streamlit dashboard:
+
+```bash
+streamlit run app.py
+```
+
+Then open the **"📷 Image Digitizer"** (now **"ECG Image Digitizer"**) page and
+either upload an ECG scan or check the sample image option.

archive/digitization/run_ecg.py

@@ -0,0 +1,302 @@
+# Force all BLAS / OpenMP thread pools to a single thread before any other
+# import. This must happen before numpy, torch, or cv2 are imported, otherwise
+# the libraries have already spawned their thread pools and the env vars have
+# no effect. Without this, BLAS and PyTorch worker threads compete for CPU
+# cores and can deadlock on multi-core servers.
+import os
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+import glob
+import signal
+import torch
+import pandas as pd
+from ultralytics import YOLO
+from digitization import ECGImage
+
+# -------------------------
+# Configuration
+# -------------------------
+
+# Root of the organised input dataset. Each immediate subdirectory is treated
+# as one hospital and must itself contain subfolders named after CATEGORIES.
+ORGANIZED_DIR = "../ecg_files/ECG_organized_all"
+
+# Root where output CSVs are written, mirroring the hospital/category structure
+# of ORGANIZED_DIR. Created automatically if it does not exist.
+OUTPUT_DIR    = "../ecg_files/ECG_digitized"
+
+# Subfolders to process within each hospital directory.
+# Any subfolder not listed here is silently skipped.
+# Add or remove entries to match your dataset's naming convention.
+CATEGORIES    = ["pre", "index", "post", "discharge", "dimi"]
+
+# Per-image wall-clock time limit in seconds, enforced via SIGALRM.
+# Images that exceed this budget are logged as timeouts and skipped.
+# Increase for very large images or slow CPU-only machines.
+# Note: SIGALRM is only available on Linux.
+TIMEOUT       = 600
+
+# -------------------------
+# Torch / GPU settings
+# -------------------------
+
+if torch.cuda.is_available():
+    # benchmark=True lets cuDNN auto-tune convolution algorithms at the first
+    # forward pass. Speeds up subsequent passes when input sizes are fixed,
+    # which is the case here (images are resized to a fixed height).
+    torch.backends.cudnn.benchmark     = True
+
+    # deterministic=False allows cuDNN to use non-deterministic (but faster)
+    # algorithms. Set to True if exact reproducibility is required.
+    torch.backends.cudnn.deterministic = False
+
+    print(f"GPU: {torch.cuda.get_device_name(0)}", flush=True)
+    print(f"VRAM: {torch.cuda.get_device_properties(0).total_memory / 1e9:.1f} GB\n",
+          flush=True)
+else:
+    print("No GPU found — running on CPU.\n", flush=True)
+
+# -------------------------
+# Timeout (Linux only)
+# -------------------------
+
+class TimeoutException(Exception):
+    """Raised by timeout_handler when SIGALRM fires."""
+    pass
+
+def timeout_handler(signum, frame):
+    """Signal handler that converts SIGALRM into a TimeoutException.
+
+    Registered below with signal.signal() so that any blocking call
+    inside process_folder() can be interrupted after TIMEOUT seconds.
+    """
+    raise TimeoutException("Processing timed out")
+
+# Register the handler for SIGALRM. signal.alarm(n) starts a countdown;
+# if the process has not called signal.alarm(0) within n seconds, this
+# handler fires and raises TimeoutException.
+signal.signal(signal.SIGALRM, timeout_handler)
+
+# -------------------------
+# Load models once (on GPU if available)
+# -------------------------
+
+# All four models are loaded once before the processing loop so that the
+# weights are only read from disk once, regardless of how many images are
+# processed. Loading inside the loop would add several seconds of overhead
+# per image.
+DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
+
+print("Loading models...", flush=True)
+
+# Detection model for full lead bounding boxes (one box per lead region).
+box_model       = YOLO("Submission code/3. Model training/runs/yolo11_full/weights/best.pt")
+
+# Classification model for lead name labels (I, II, III, aVR, aVL, aVF, V1–V6).
+lead_name_model = YOLO("Submission code/3. Model training/runs/yolo11_lead/weights/best.pt")
+
+# Detection model for calibration (reference) pulse boxes.
+pulse_model     = YOLO("Submission code/3. Model training/runs/yolo11_pulse/weights/best.pt")
+
+box_model.to(DEVICE)
+lead_name_model.to(DEVICE)
+pulse_model.to(DEVICE)
+
+# The segmentation model is kept as a file path rather than a loaded model
+# object because it is consumed by patched_yolo_infer.MakeCropsDetectThem,
+# which loads it internally on each call. This avoids a second copy of the
+# weights living in VRAM alongside the three detection models above.
+segmentation_model = "Submission code/3. Model training/runs/yolo11_patch/weights/best.pt"
+
+print(f"Models loaded on {DEVICE.upper()}.\n", flush=True)
+
+# -------------------------
+# Process a single folder
+# -------------------------
+
+def process_folder(input_dir, output_dir):
+    """Digitize all ECG images in input_dir and save CSVs to output_dir.
+
+    Iterates over all .png / .jpg / .jpeg files in input_dir (sorted
+    alphabetically). Each image is processed by ECGImage.run_full_pipeline()
+    and the resulting signals are exported as a CSV file with the same base
+    name as the image.
+
+    Parameters
+    ----------
+    input_dir : str
+        Directory containing the input ECG images.
+    output_dir : str
+        Directory where output CSV files are written (created if absent).
+
+    Returns
+    -------
+    success : int
+        Number of images successfully digitized.
+    total : int
+        Total number of images found in input_dir.
+    """
+    # Collect all supported image formats; sort for reproducible ordering.
+    files = []
+    for ext in ("*.png", "*.jpg", "*.jpeg"):
+        files.extend(glob.glob(os.path.join(input_dir, ext)))
+    files.sort()
+
+    # Return early without creating the output directory if the folder is empty.
+    if not files:
+        return 0, 0
+
+    os.makedirs(output_dir, exist_ok=True)
+    success = errors = timeouts = 0
+
+    for i, file_path in enumerate(files):
+        base_name = os.path.splitext(os.path.basename(file_path))[0]
+        print(f"    [{i+1}/{len(files)}] {base_name}", flush=True)
+        try:
+            # Start the per-image countdown. If the block below takes longer
+            # than TIMEOUT seconds, SIGALRM fires and TimeoutException is raised.
+            signal.alarm(TIMEOUT)
+
+            # Instantiate the digitizer for this image. Models are passed in
+            # so they are not reloaded on every iteration.
+            ecg = ECGImage(
+                box_model=box_model,
+                segmentation_model=segmentation_model,
+                lead_name_model=lead_name_model,
+                pulse_model=pulse_model,
+                image_path=file_path,
+            )
+
+            # Run all pipeline stages: preprocessing → segmentation → YOLO
+            # detection (×3) → scale calibration → grid construction →
+            # signal extraction. Internally retries at different image sizes
+            # if any detection returns empty.
+            ecg.run_full_pipeline()
+
+            # Write one CSV per image: columns = lead names, rows = time samples.
+            ecg.save_signals_as_csv(base_name, directory=output_dir)
+            success += 1
+
+        except TimeoutException:
+            # The image took longer than TIMEOUT seconds. Counted separately
+            # from errors so they can be investigated independently.
+            timeouts += 1
+            print(f"    ⏰ TIMEOUT: {base_name}", flush=True)
+
+        except Exception as e:
+            # Any other failure (bad image, detection failure, calibration
+            # error, etc.) is caught here so the loop continues with the
+            # remaining images.
+            errors += 1
+            print(f"    ❌ ERROR: {base_name} → {e}", flush=True)
+
+        finally:
+            # Always cancel the alarm, even if an exception was raised, to
+            # prevent a stale countdown from interrupting the next image.
+            signal.alarm(0)
+
+            # Release VRAM after each image to prevent fragmentation across
+            # a long batch run. synchronize() ensures all GPU ops are complete
+            # before empty_cache() frees the memory.
+            if torch.cuda.is_available():
+                try:
+                    torch.cuda.synchronize()
+                    torch.cuda.empty_cache()
+                except Exception:
+                    pass
+
+    total = len(files)
+    print(f"    → {success}/{total} ok | {errors} errors | {timeouts} timeouts\n",
+          flush=True)
+    return success, total
+
+# -------------------------
+# Main loop
+# -------------------------
+
+# Discover hospital directories. Each immediate child directory of
+# ORGANIZED_DIR is treated as one hospital; files at the root level are ignored.
+hospitals = sorted([
+    name for name in os.listdir(ORGANIZED_DIR)
+    if os.path.isdir(os.path.join(ORGANIZED_DIR, name))
+])
+
+# Accumulate per-hospital/category counts for the summary report.
+log_rows = []
+
+print(f"Found {len(hospitals)} hospital(s). Starting digitization...\n", flush=True)
+
+for hospital in hospitals:
+    print(f"{'='*60}", flush=True)
+    print(f"Hospital: {hospital}", flush=True)
+    print(f"{'='*60}", flush=True)
+
+    for cat in CATEGORIES:
+        input_dir  = os.path.join(ORGANIZED_DIR, hospital, cat)
+        output_dir = os.path.join(OUTPUT_DIR,    hospital, cat)
+
+        # Skip silently if this category does not exist for this hospital.
+        if not os.path.isdir(input_dir):
+            continue
+
+        print(f"  [{cat.upper()}]", flush=True)
+        success, total = process_folder(input_dir, output_dir)
+
+        # Only log categories that had at least one image, to keep the
+        # summary table free of zero-row entries.
+        if total > 0:
+            log_rows.append({
+                "hospital": hospital,
+                "category": cat,
+                "total"   : total,
+                "success" : success,
+                "failed"  : total - success,
+                "rate"    : f"{100 * success / total:.1f}%",
+            })
+
+# -------------------------
+# Summary report
+# -------------------------
+
+print("\n" + "="*60, flush=True)
+print("DIGITIZATION SUMMARY", flush=True)
+print("="*60, flush=True)
+
+df = pd.DataFrame(log_rows)
+
+if df.empty:
+    print("No images were processed.", flush=True)
+else:
+    # Flat per-hospital/category breakdown.
+    print("\nPer hospital / category:\n")
+    print(df.to_string(index=False))
+
+    # Pivot to a hospital × category matrix of successful digitizations,
+    # with a TOTAL row and column for quick inspection.
+    print("\n\nSuccessful digitizations — distribution:\n")
+    success_pivot = (
+        df.pivot_table(index="hospital", columns="category",
+                       values="success", aggfunc="sum", fill_value=0)
+          .reindex(columns=CATEGORIES, fill_value=0)  # preserve CATEGORIES column order
+    )
+    success_pivot["TOTAL"] = success_pivot.sum(axis=1)   # row totals
+    success_pivot.loc["TOTAL"] = success_pivot.sum(axis=0)  # grand total row
+    print(success_pivot.to_string())
+
+    # Overall counts across the entire dataset.
+    total_imgs   = df["total"].sum()
+    total_ok     = df["success"].sum()
+    total_failed = df["failed"].sum()
+    print(f"\n✅ Overall: {total_ok}/{total_imgs} successful "
+          f"({100*total_ok/total_imgs:.1f}%) | "
+          f"{total_failed} failed", flush=True)
+
+    # Persist the pivot table as a CSV for later analysis.
+    summary_path = os.path.join(OUTPUT_DIR, "digitization_summary.csv")
+    os.makedirs(OUTPUT_DIR, exist_ok=True)
+    success_pivot.to_csv(summary_path)
+    print(f"📄 Summary saved to {summary_path}", flush=True)
+
+print("\nAll processing complete.", flush=True)

backend/__init__.py

@@ -0,0 +1 @@
+# backend package

backend/classification_runner.py

@@ -0,0 +1,571 @@
+# backend/classification_runner.py
+"""
+Backend adapter for ECG Classification.
+Handles heartbeat segmentation, data preparation for sktime models, and model evaluation.
+"""
+
+import os
+import sys
+import numpy as np
+import pandas as pd
+from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
+
+# Add parent directory to path to ensure digitization and other local modules are importable
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
+
+# Custom Exceptions for dependencies
+class DependencyMissingError(Exception):
+    """Raised when a required third-party library is not installed."""
+    pass
+
+class EmptyDatasetError(Exception):
+    """Raised when the uploaded dataset is empty or invalid."""
+    pass
+
+# Helper to check dependencies
+def check_classification_dependencies():
+    """
+    Checks if sktime is installed in the current environment.
+    Raises DependencyMissingError if it is not.
+    """
+    missing = []
+    try:
+        import sktime
+    except ImportError:
+        missing.append("sktime")
+        
+    if missing:
+        raise DependencyMissingError(
+            f"The following libraries are required for classification but are missing in this environment:\n"
+            f"{', '.join(missing)}\n\n"
+            f"Please run: pip install sktime"
+        )
+
+# ==============================================================================
+# HEARTBEAT SEGMENTATION
+# ==============================================================================
+# ==============================================================================
+# LOCAL HEARTBEAT SEGMENTATION IMPLEMENTATION (copied from ecg_segmentation.py)
+# ==============================================================================
+
+def bandpass_filter(signal_data, lowcut=5, highcut=15, fs=500, order=2):
+    """Bandpass filter for ECG signal"""
+    from scipy.signal import butter, filtfilt
+    try:
+        nyquist = 0.5 * fs
+        low = lowcut / nyquist
+        high = highcut / nyquist
+        b, a = butter(order, [low, high], btype='band')
+        return filtfilt(b, a, signal_data)
+    except Exception:
+        return signal_data
+
+
+def find_r_peaks(ecg_signal, min_heart_rate=40, max_heart_rate=180, fs=500):
+    """Find R-peaks in ECG signal"""
+    from scipy.signal import find_peaks
+    # Calculate expected distance between peaks based on heart rate
+    min_distance = int((60 / max_heart_rate) * fs)  # samples
+    max_distance = int((60 / min_heart_rate) * fs)  # samples
+    
+    # Find peaks
+    peaks, properties = find_peaks(ecg_signal, 
+                                  height=np.percentile(ecg_signal, 80),
+                                  distance=min_distance,
+                                  prominence=np.std(ecg_signal) * 0.3)
+    
+    # Filter peaks based on reasonable heart rate
+    if len(peaks) == 0:
+        return peaks
+        
+    valid_peaks = [peaks[0]]
+    for i in range(1, len(peaks)):
+        distance = peaks[i] - valid_peaks[-1]
+        if min_distance <= distance <= max_distance:
+            valid_peaks.append(peaks[i])
+    
+    return np.array(valid_peaks)
+
+
+def downsample_beat(beat_data, features, downsampling_factor):
+    """
+    Downsample a single heartbeat by taking every nth sample.
+    Uses vectorized iloc slicing instead of row-by-row iteration.
+    """
+    if downsampling_factor <= 1:
+        return beat_data
+    try:
+        downsampled = beat_data.iloc[::downsampling_factor]
+        return downsampled if len(downsampled) > 0 else None
+    except Exception as e:
+        print(f"Error in downsampling: {e}")
+        return None
+
+
+def segment_beats_around_r_peaks(subject_data, r_peaks, features, fs=500, target_fs=500, pre_r_ms=150, post_r_ms=250):
+    """
+    Segment heartbeats around R-peaks with downsampling
+    """
+    heartbeats = []
+    total_samples = len(subject_data)
+    
+    # Calculate downsampling factor
+    downsampling_factor = fs // target_fs
+    if downsampling_factor < 1:
+        downsampling_factor = 1
+    
+    # Convert milliseconds to samples at original and target frequencies
+    pre_r_original = int((pre_r_ms / 1000) * fs)
+    post_r_original = int((post_r_ms / 1000) * fs)
+    
+    min_beat_len = int((pre_r_original + post_r_original) * 0.7)
+    
+    for r_peak in r_peaks:
+        start_idx = max(0, r_peak - pre_r_original)
+        end_idx = min(total_samples, r_peak + post_r_original)
+        
+        # Ensure we have enough samples for a complete beat
+        if (end_idx - start_idx) >= min_beat_len:
+            beat_data_original = subject_data.iloc[start_idx:end_idx]
+            
+            # Downsample the beat
+            beat_data_downsampled = downsample_beat(beat_data_original, features, downsampling_factor)
+            
+            if beat_data_downsampled is not None:
+                heartbeats.append((beat_data_downsampled, r_peak))
+    
+    return heartbeats
+
+
+def local_segment_heartbeats(input_csv, output_csv=None, min_heart_rate=40, max_heart_rate=180, 
+                             fs=500, target_fs=500):
+    """
+    Segment ECG data into individual heartbeats for each subject with downsampling.
+    Self-contained implementation.
+    Optimized: uses vectorized pandas concat instead of row-by-row dict assembly.
+    """
+    print("Loading ECG data for heartbeat segmentation...")
+    df = pd.read_csv(input_csv, index_col=['subject_id', 'timestamp'])
+    features = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+    
+    # Fill NaN values in features using linear interpolation and fallback to 0.0
+    existing_features = [f for f in features if f in df.columns]
+    if existing_features:
+        df[existing_features] = df[existing_features].interpolate(method='linear', limit_direction='both').fillna(0.0)
+            
+    # Calculate downsampling factor
+    downsampling_factor = fs // target_fs
+    print(f"Downsampling from {fs} Hz to {target_fs} Hz (factor: {downsampling_factor})")
+    
+    # Collect beat DataFrames in a list for a single pd.concat at the end
+    beat_frames = []
+    
+    # Get unique subjects
+    subject_ids = df.index.get_level_values('subject_id').unique()
+    
+    for subject_id in subject_ids:
+        # Get subject data
+        subject_data = df.xs(subject_id, level='subject_id')
+        subject_class = subject_data['class'].iloc[0]
+        
+        # Use Lead II for R-peak detection (commonly used)
+        ecg_signal = subject_data['II'].values
+        
+        # Skip if signal is too short
+        if len(ecg_signal) < fs * 2:  # Less than 2 seconds
+            continue
+            
+        # Preprocess signal (apply bandpass filter for robust R-peak detection)
+        try:
+            ecg_filtered = bandpass_filter(ecg_signal, lowcut=5, highcut=15, fs=fs, order=2)
+        except Exception:
+            ecg_filtered = ecg_signal
+        
+        # Find R-peaks
+        r_peaks = find_r_peaks(ecg_filtered, min_heart_rate=min_heart_rate, max_heart_rate=max_heart_rate, fs=fs)
+        
+        if len(r_peaks) < 2:
+            # Fallback: if filtered failed or returned no peaks, try on unfiltered signal
+            try:
+                r_peaks = find_r_peaks(ecg_signal, min_heart_rate=min_heart_rate, max_heart_rate=max_heart_rate, fs=fs)
+            except Exception:
+                pass
+            
+        if len(r_peaks) < 2:
+            print(f"  Warning: Only {len(r_peaks)} R-peaks found for subject {subject_id}")
+            continue
+        
+        # Segment heartbeats around R-peaks
+        heartbeats = segment_beats_around_r_peaks(subject_data, r_peaks, features, fs=fs, target_fs=target_fs)
+        
+        # Vectorized beat assembly: build each beat as a small DataFrame, then concat once
+        for i, (beat_data, r_peak_pos) in enumerate(heartbeats):
+            beat_id = f"{subject_id}_beat{i+1}"
+            n_rows = len(beat_data)
+            
+            # Build the beat frame from the existing beat_data slice
+            beat_frame = beat_data[existing_features].copy()
+            beat_frame = beat_frame.reset_index(drop=True)
+            beat_frame['subject_id'] = beat_id
+            beat_frame['original_subject'] = subject_id
+            beat_frame['beat_number'] = i + 1
+            beat_frame['r_peak_position'] = int(r_peak_pos)
+            beat_frame['class'] = subject_class
+            beat_frame['timestamp'] = np.arange(n_rows)
+            
+            beat_frames.append(beat_frame)
+        
+    # Single concat at the end instead of thousands of dict appends
+    if not beat_frames:
+        return pd.DataFrame()
+        
+    segmented_df = pd.concat(beat_frames, ignore_index=True)
+    
+    # Set multi-index
+    if not segmented_df.empty:
+        segmented_df.set_index(['subject_id', 'timestamp'], inplace=True)
+    
+    # Save if output path provided
+    if output_csv and not segmented_df.empty:
+        os.makedirs(os.path.dirname(output_csv), exist_ok=True)
+        segmented_df.to_csv(output_csv)
+        print(f"Segmented data saved to: {output_csv}")
+    
+    print(f"\nSegmentation complete:")
+    print(f"Total heartbeats: {segmented_df.index.get_level_values('subject_id').nunique() if not segmented_df.empty else 0}")
+    
+    return segmented_df
+
+
+def segment_uploaded_csv(input_csv_path, output_csv_path, target_fs=500):
+    """
+    Segments a raw digitized dataset CSV into individual heartbeats.
+    Uses the local self-contained heartbeat segmentation implementation.
+    """
+    # Read the first few lines to inspect columns
+    df_raw = pd.read_csv(input_csv_path)
+    if df_raw.empty:
+        raise EmptyDatasetError("The uploaded CSV file is empty.")
+        
+    # Standardize input formatting for the segmentation module.
+    # The segmentation module expects 'subject_id', 'timestamp', 'class' columns.
+    required_cols = ['subject_id', 'timestamp', 'class']
+    needs_formatting = not all(col in df_raw.columns for col in required_cols)
+    
+    formatted_csv_path = input_csv_path
+    if needs_formatting:
+        print("Formatting raw CSV columns to match expected segmentation schema...")
+        # Check if we have lead columns
+        lead_cols = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+        existing_leads = [col for col in lead_cols if col in df_raw.columns]
+        if not existing_leads:
+            raise ValueError(
+                f"The uploaded CSV does not contain standard ECG lead columns.\n"
+                f"Expected some of: {lead_cols}"
+            )
+            
+        # Format df_raw
+        df_formatted = df_raw.copy()
+        if 'subject_id' not in df_formatted.columns:
+            df_formatted['subject_id'] = 'subject_1'
+        if 'class' not in df_formatted.columns:
+            # Default placeholder label
+            df_formatted['class'] = 'Pre-Procedural MI'
+        if 'timestamp' not in df_formatted.columns:
+            df_formatted['timestamp'] = range(len(df_formatted))
+            
+        # Ensure subject_id and class are strings
+        df_formatted['subject_id'] = df_formatted['subject_id'].astype(str)
+        df_formatted['class'] = df_formatted['class'].astype(str)
+        
+        # Save to temporary path for segmentation input
+        formatted_csv_path = input_csv_path + "_formatted.csv"
+        df_formatted.to_csv(formatted_csv_path, index=False)
+        
+    # Run the local segmentation module
+    try:
+        segmented_df = local_segment_heartbeats(
+            formatted_csv_path,
+            output_csv=output_csv_path,
+            fs=target_fs,
+            target_fs=target_fs
+        )
+    finally:
+        # Clean up temporary formatted file if created
+        if needs_formatting and os.path.exists(formatted_csv_path):
+            try:
+                os.remove(formatted_csv_path)
+            except Exception:
+                pass
+                
+    if segmented_df is None or segmented_df.empty:
+        raise ValueError(
+            "Heartbeat segmentation returned no beats. This can happen if the signal is too short or R-peaks were not detected."
+        )
+        
+    return segmented_df
+
+
+# 
+
+
+# ==============================================================================
+# PRE-TRAINED MODEL LOADING & INFERENCE SUPPORT
+# ==============================================================================
+
+def _load_model_and_metadata_impl(model_key):
+    """
+    Internal implementation: loads a pre-trained model and its metadata from
+    models/classifier_models/{model_key}.
+    """
+    import json
+    import pickle
+
+    model_dir = os.path.join("models", "classifier_models", model_key)
+    metadata_path = os.path.join(model_dir, "model_metadata.json")
+
+    if not os.path.exists(metadata_path):
+        raise FileNotFoundError(
+            f"Pre-trained model metadata not found: {metadata_path}\n"
+            f"Please verify models/classifier_models/ directory contents."
+        )
+
+    with open(metadata_path, 'r') as f:
+        metadata = json.load(f)
+
+    model_path = os.path.join(model_dir, metadata['model_file'])
+    if not os.path.exists(model_path):
+        raise FileNotFoundError(
+            f"Pre-trained model pickle not found: {model_path}\n"
+            f"Please verify models/classifier_models/ directory contents."
+        )
+
+    print(f"Loading {metadata['model_name']} pre-trained model for {metadata['use_case']}...")
+    with open(model_path, 'rb') as f:
+        model = pickle.load(f)
+
+    # Restore missing classes_ on unpickled scikit-learn models due to version inconsistencies
+    def _fix_unpickled_classes(obj, visited=None):
+        if obj is None:
+            return
+        if visited is None:
+            visited = set()
+        obj_id = id(obj)
+        if obj_id in visited:
+            return
+        visited.add(obj_id)
+        
+        if hasattr(obj, '_label_binarizer') and hasattr(obj._label_binarizer, 'classes_'):
+            if not hasattr(obj, 'classes_') or obj.classes_ is None:
+                try:
+                    obj.classes_ = obj._label_binarizer.classes_
+                except Exception:
+                    pass
+                    
+        if isinstance(obj, (list, tuple)):
+            for item in obj:
+                _fix_unpickled_classes(item, visited)
+        elif isinstance(obj, dict):
+            for val in obj.values():
+                _fix_unpickled_classes(val, visited)
+        else:
+            for attr in ['estimator', 'estimator_', 'classifier', 'classifier_', 
+                         'estimators', 'estimators_', 'transformers', 'transformers_', 
+                         'steps', 'steps_']:
+                if hasattr(obj, attr):
+                    try:
+                        val = getattr(obj, attr)
+                        _fix_unpickled_classes(val, visited)
+                    except Exception:
+                        pass
+
+    _fix_unpickled_classes(model)
+
+    return model, metadata
+
+
+# Try to use Streamlit cache; fall back to plain function if running outside Streamlit
+try:
+    import streamlit as st
+
+    @st.cache_resource(show_spinner="Loading pre-trained model...")
+    def load_pretrained_model_and_metadata(model_key):
+        """Cached wrapper — model pickle is only deserialized once per model_key."""
+        return _load_model_and_metadata_impl(model_key)
+except Exception:
+    # Fallback when Streamlit is not available (e.g. tests, scripts)
+    load_pretrained_model_and_metadata = _load_model_and_metadata_impl
+
+
+def preprocess_dataframe_for_inference(df, metadata):
+    """
+    Preprocesses raw or segmented DataFrame into numpy3D array for the pre-trained model.
+    Reshapes each subject's signal to match the exact target length (n_timesteps),
+    applies max absolute voltage normalization per subject, and returns the 3D numpy array.
+    Optimized: single grouped normalization instead of per-column loop.
+    """
+    features = metadata['features']
+    n_timesteps = metadata['n_timesteps']
+    n_features = metadata['n_features']
+
+    # Standardize input formatting
+    # If index is multi-index with subject_id, or if subject_id/timestamp are columns
+    df_reset = df.copy()
+    if 'subject_id' in df_reset.index.names:
+        df_reset = df_reset.reset_index()
+    elif 'subject_id' not in df_reset.columns:
+        df_reset['subject_id'] = 'subject_1'
+
+    if 'timestamp' in df_reset.index.names:
+        df_reset = df_reset.reset_index()
+    elif 'timestamp' not in df_reset.columns:
+        df_reset['timestamp'] = df_reset.groupby('subject_id').cumcount()
+
+    df_reset = df_reset.set_index(['subject_id', 'timestamp'])
+
+    # Verify required features exist
+    missing = [f for f in features if f not in df_reset.columns]
+    if missing:
+        # Fallback to standard channels if names differ slightly (e.g. casing)
+        df_reset.columns = [c.strip() for c in df_reset.columns]
+        missing = [f for f in features if f not in df_reset.columns]
+        if missing:
+            raise ValueError(f"Missing required lead columns in input: {missing}")
+
+    # Vectorized normalization: single grouped abs-max across all feature columns at once
+    df_X = df_reset[features].copy()
+    grouped_max = df_X.abs().groupby(level='subject_id').transform('max')
+    # Avoid division by zero
+    grouped_max = grouped_max.replace(0, 1.0)
+    df_X = df_X / grouped_max
+
+    # Reshape to 3D array: (n_subjects, n_features, n_timesteps)
+    subject_ids = df_X.index.get_level_values('subject_id').unique()
+    X_list = []
+    valid_subject_ids = []
+    y_true_str = []
+
+    for subject_id in subject_ids:
+        try:
+            subject_data = df_X.xs(subject_id, level='subject_id')
+            actual_len = len(subject_data)
+
+            if actual_len >= n_timesteps:
+                # Truncate
+                arr = subject_data.iloc[:n_timesteps].values.T
+            else:
+                # Zero-pad
+                arr = np.zeros((n_features, n_timesteps))
+                arr[:, :actual_len] = subject_data.values.T
+
+            if arr.shape == (n_features, n_timesteps):
+                X_list.append(arr)
+                valid_subject_ids.append(subject_id)
+                
+                # Capture ground truth class label if it exists in the source DataFrame
+                if 'class' in df_reset.columns:
+                    lbl = df_reset.loc[subject_id].iloc[0]['class']
+                    y_true_str.append(str(lbl))
+        except Exception:
+            continue
+
+    if not X_list:
+        raise ValueError(
+            f"No valid subjects/heartbeats found after preprocessing. "
+            f"Ensure data contains at least one subject with some timesteps."
+        )
+
+    return np.array(X_list), valid_subject_ids, y_true_str
+
+
+def run_pretrained_inference(model, X, metadata, subject_ids):
+    """
+    Runs pre-trained model predictions and returns a results dictionary.
+    """
+    from joblib import parallel_backend
+
+    with parallel_backend('threading', n_jobs=-1):
+        y_pred = model.predict(X)
+
+        has_proba = False
+        y_proba = None
+        try:
+            y_proba = model.predict_proba(X)
+            has_proba = True
+        except (AttributeError, NotImplementedError):
+            pass
+
+    results = {
+        'subject_id': subject_ids,
+        'predicted_class': y_pred.tolist(),
+    }
+
+    if has_proba and y_proba is not None:
+        class_labels = list(model.classes_)
+        results['class_labels'] = class_labels
+        results['probabilities'] = y_proba.tolist()
+
+        # Calculate confidence as the probability of the predicted class
+        confidence = []
+        for j in range(len(y_pred)):
+            try:
+                pred_idx = class_labels.index(y_pred[j])
+                confidence.append(float(y_proba[j, pred_idx]))
+            except ValueError:
+                confidence.append(1.0)
+        results['confidence'] = confidence
+
+    return results
+
+
+def calculate_evaluation_metrics(y_true_str, y_pred_str, class_labels, positive_class):
+    """
+    Calculates accuracy, precision, recall, f1, sensitivity, specificity, and confusion matrix.
+    Maps ground-truth and prediction string labels to binary (0/1).
+    """
+    # Map strings to binary values (0/1) based on positive class
+    pos_str = str(positive_class).strip().lower()
+    
+    # We do a lenient comparison (check substring or exact match)
+    y_true_num = [1 if pos_str in str(y).strip().lower() else 0 for y in y_true_str]
+    y_pred_num = [1 if pos_str in str(y).strip().lower() else 0 for y in y_pred_str]
+
+    # Handle cases where all targets are of one class
+    if len(np.unique(y_true_num)) <= 1 and len(np.unique(y_pred_num)) <= 1:
+        # Avoid standard metrics crashing
+        accuracy = float(accuracy_score(y_true_num, y_pred_num))
+        return {
+            "Accuracy": accuracy,
+            "Precision": 1.0 if accuracy == 1.0 else 0.0,
+            "Recall": 1.0 if accuracy == 1.0 else 0.0,
+            "F1": 1.0 if accuracy == 1.0 else 0.0,
+            "Sensitivity": 1.0 if accuracy == 1.0 else 0.0,
+            "Specificity": 1.0 if accuracy == 1.0 else 0.0,
+            "Confusion Matrix": confusion_matrix(y_true_num, y_pred_num, labels=[0, 1]).tolist(),
+            "TN": int(accuracy * len(y_true_num)) if y_true_num[0] == 0 else 0,
+            "FP": 0,
+            "FN": 0,
+            "TP": int(accuracy * len(y_true_num)) if y_true_num[0] == 1 else 0
+        }
+
+    accuracy = float(accuracy_score(y_true_num, y_pred_num))
+    precision = float(precision_score(y_true_num, y_pred_num, zero_division=0))
+    recall = float(recall_score(y_true_num, y_pred_num, zero_division=0))
+    f1 = float(f1_score(y_true_num, y_pred_num, zero_division=0))
+
+    cm = confusion_matrix(y_true_num, y_pred_num, labels=[0, 1])
+    tn, fp, fn, tp = cm.flatten()
+    
+    sensitivity = float(tp / (tp + fn) if (tp + fn) > 0 else 0)
+    specificity = float(tn / (tn + fp) if (tn + fp) > 0 else 0)
+
+    return {
+        "Accuracy": accuracy,
+        "Precision": precision,
+        "Recall": recall,
+        "F1": f1,
+        "Sensitivity": sensitivity,
+        "Specificity": specificity,
+        "Confusion Matrix": cm.tolist(),
+        "TN": int(tn), "FP": int(fp), "FN": int(fn), "TP": int(tp)
+    }

backend/digitization_runner.py

@@ -0,0 +1,174 @@
+# backend/digitization_runner.py
+"""
+Backend adapter for ECG Digitization.
+Wraps the ECGImage pipeline from digitization.py and manages model caching.
+"""
+
+import os
+import sys
+import pandas as pd
+import numpy as np
+import streamlit as st
+import torch
+
+# Add parent directory to path to ensure digitization and other local modules are importable
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
+
+try:
+    from digitization import ECGImage
+    from ultralytics import YOLO
+except ImportError as e:
+    # We will log the error but keep running so the dashboard can load and display elegant warnings
+    pass
+
+class ModelPathNotSetError(Exception):
+    """Raised when one or more YOLO model paths are not set in the configuration."""
+    pass
+
+class ModelFileNotFoundError(Exception):
+    """Raised when a specified YOLO model file does not exist on disk."""
+    pass
+
+def check_models_configured(config):
+    """
+    Verify that all YOLO model paths are set in the configuration.
+    """
+    missing_configs = []
+    if not config.YOLO_BOX_MODEL_PATH:
+        missing_configs.append("YOLO Box Model (YOLO_BOX_MODEL_PATH)")
+    if not config.YOLO_LEAD_NAME_MODEL_PATH:
+        missing_configs.append("YOLO Lead Name Model (YOLO_LEAD_NAME_MODEL_PATH)")
+    if not config.YOLO_PULSE_MODEL_PATH:
+        missing_configs.append("YOLO Pulse Model (YOLO_PULSE_MODEL_PATH)")
+    if not config.YOLO_SEGMENTATION_MODEL_PATH:
+        missing_configs.append("YOLO Segmentation Model (YOLO_SEGMENTATION_MODEL_PATH)")
+        
+    if missing_configs:
+        raise ModelPathNotSetError(
+            f"The following YOLO model paths are not set in config.py:\n" + 
+            "\n".join([f"- {m}" for m in missing_configs]) + 
+            "\n\nPlease configure these paths in config.py before running digitization."
+        )
+
+def check_model_files_exist(config):
+    """
+    Verify that the configured model files actually exist on disk.
+    """
+    missing_files = []
+    for name, path in [
+        ("Box Model", config.YOLO_BOX_MODEL_PATH),
+        ("Lead Name Model", config.YOLO_LEAD_NAME_MODEL_PATH),
+        ("Pulse Model", config.YOLO_PULSE_MODEL_PATH),
+        ("Segmentation Model", config.YOLO_SEGMENTATION_MODEL_PATH)
+    ]:
+        if path and not os.path.exists(path):
+            missing_files.append(f"{name} at path: '{path}'")
+            
+    if missing_files:
+        raise ModelFileNotFoundError(
+            f"The following model files were not found on disk:\n" + 
+            "\n".join([f"- {f}" for f in missing_files]) + 
+            "\n\nPlease verify that the files are downloaded and located at these paths."
+        )
+
+@st.cache_resource
+def load_yolo_models(box_path, lead_name_path, pulse_path):
+    """
+    Load YOLO models once and cache them in Streamlit memory.
+    Returns a dict containing the loaded model objects.
+    """
+    device = "cuda" if torch.cuda.is_available() else "cpu"
+    
+    # Load each model
+    box_model = YOLO(box_path)
+    lead_name_model = YOLO(lead_name_path)
+    pulse_model = YOLO(pulse_path)
+    
+    # Move to GPU if available
+    box_model.to(device)
+    lead_name_model.to(device)
+    pulse_model.to(device)
+    
+    return {
+        "box_model": box_model,
+        "lead_name_model": lead_name_model,
+        "pulse_model": pulse_model,
+        "device": device
+    }
+
+def run_digitization_pipeline(image_path, config):
+    """
+    Runs the full digitization pipeline on a single ECG image.
+    
+    Parameters
+    ----------
+    image_path : str
+        Path to the uploaded image file.
+    config : module
+        The config module containing paths and settings.
+        
+    Returns
+    -------
+    df : pandas.DataFrame
+        Digitized ECG signals with lead names as columns and timesteps as rows.
+    ecg_instance : ECGImage
+        The raw ECGImage instance for accessing signal_grid, masks, or metadata.
+    """
+    # 1. Validation
+    check_models_configured(config)
+    check_model_files_exist(config)
+    
+    # 2. Load models (cached via st.cache_resource)
+    models = load_yolo_models(
+        config.YOLO_BOX_MODEL_PATH,
+        config.YOLO_LEAD_NAME_MODEL_PATH,
+        config.YOLO_PULSE_MODEL_PATH
+    )
+    
+    # 3. Instantiate ECGImage
+    ecg = ECGImage(
+        box_model=models["box_model"],
+        segmentation_model=config.YOLO_SEGMENTATION_MODEL_PATH,
+        lead_name_model=models["lead_name_model"],
+        pulse_model=models["pulse_model"],
+        image_path=image_path
+    )
+    
+    # 4. Run the pipeline (without Linux-specific SIGALRM timeouts)
+    ecg.run_full_pipeline()
+    
+    # 5. Extract DataFrame from signal_grid
+    df = extract_dataframe_from_ecg(ecg)
+    
+    # 6. Clear GPU cache if needed
+    if torch.cuda.is_available():
+        try:
+            torch.cuda.synchronize()
+            torch.cuda.empty_cache()
+        except Exception:
+            pass
+            
+    return df, ecg
+
+def extract_dataframe_from_ecg(ecg_instance):
+    """
+    Extracts signal data from an ECGImage instance and converts it to a pandas DataFrame.
+    """
+    signals, lead_names, max_length = [], [], 0
+    for row in ecg_instance.signal_grid:
+        for cell in row:
+            if 'signal' in cell and 'lead' in cell:
+                sig = np.asarray(cell['signal'], dtype=float)
+                signals.append(sig)
+                lead_names.append(cell['lead'])
+                max_length = max(max_length, len(sig))
+                
+    if not signals:
+        raise ValueError("No valid lead signals were extracted from the ECG image.")
+        
+    # Pad shorter signals with NaN to make them uniform length
+    padded = [np.pad(s, (0, max_length - len(s)), mode='constant', constant_values=np.nan) for s in signals]
+    
+    # Build DataFrame
+    df = pd.DataFrame(np.vstack(padded).T, columns=lead_names)
+    return df

config.py

@@ -0,0 +1,64 @@
+# config.py
+"""
+Centralized configuration settings for the ECG Digitization and Classification Dashboard.
+This holds paths, hyperparameter defaults, and metadata for classification tasks.
+"""
+
+import os
+
+# ==============================================================================
+# YOLO Models Settings (For Digitization)
+# ==============================================================================
+# These paths are initially left empty as requested. The user will set them in their environment 
+# or directly in this file. The Streamlit app will warn if they are empty.
+YOLO_BOX_MODEL_PATH = "models/digitization_models/yolo11_full/weights/best.pt"          # Bounding boxes for full lead regions
+YOLO_LEAD_NAME_MODEL_PATH = "models/digitization_models/yolo11_lead/weights/best.pt"    # Classification for lead labels (I, II, etc.)
+YOLO_PULSE_MODEL_PATH = "models/digitization_models/yolo11_pulse/weights/best.pt"        # Detection for calibration reference pulses
+YOLO_SEGMENTATION_MODEL_PATH = "models/digitization_models/yolo11_patch/weights/best.pt" # Path to the YOLO patch segmentation model checkpoint
+
+# ==============================================================================
+# ECG Signal Metadata
+# ==============================================================================
+# Target sampling frequency for digitization and classification (in Hz)
+TARGET_FS = 500
+
+# Ordered standard list of 12 ECG Leads
+LEAD_NAMES = ['I', 'aVR', 'V1', 'V4', 'II', 'aVL', 'V2', 'V5', 'III', 'aVF', 'V3', 'V6']
+
+# ==============================================================================
+# Classification Task Mapping
+# ==============================================================================
+# Maps the task selected in the Streamlit UI to its associated model and target class labels.
+CLASSIFICATION_TASKS = {
+    "Normal vs Myocardial Infarction (MI) - Segmented": {
+        "model_dir": "mi_vs_normal_segmented",
+        "model": "Arsenal",
+        "labels": {0: "NORMAL", 1: "MYOCARDIAL_INFARCTION"},
+        "description": "Differentiates normal healthy heartbeats from Myocardial Infarction using the pre-trained Arsenal ensemble classifier (Segmented Heartbeats)."
+    },
+    "Occlusive Myocardial Infarction (OMI) vs non-OMI": {
+        "model_dir": "omi_vs_nonomi",
+        "model": "Rocket",
+        "labels": {0: "non-OMI", 1: "OMI"},
+        "description": "Identifies Occlusive Myocardial Infarction (OMI) vs Non-Occlusive Myocardial Infarction/controls using the pre-trained Rocket classifier."
+    },
+    "Pre-Procedural vs Post-Procedural MI": {
+        "model_dir": "ecg_surgery",
+        "model": "InceptionTime",
+        "labels": {0: "post-procedural MI", 1: "pre-procedural MI"},
+        "description": "Classifies pre-procedural MI vs post-procedural MI using the pre-trained InceptionTime deep learning time-series classifier."
+    }
+}
+
+
+# ==============================================================================
+# File Management
+# ==============================================================================
+# Base output directories for storing run results
+OUTPUT_BASE_DIR = "output"
+DIGITIZATION_OUTPUT_DIR = os.path.join(OUTPUT_BASE_DIR, "digitization")
+CLASSIFICATION_OUTPUT_DIR = os.path.join(OUTPUT_BASE_DIR, "classification")
+
+# Ensure base directories exist
+os.makedirs(DIGITIZATION_OUTPUT_DIR, exist_ok=True)
+os.makedirs(CLASSIFICATION_OUTPUT_DIR, exist_ok=True)

digitization.py

@@ -0,0 +1,1784 @@
+"""
+digitization.py
+===============
+Core ECG image-to-signal conversion module.
+
+The public interface is the :class:`ECGImage` class.  Instantiate it with
+four pre-loaded YOLO models and the path to an ECG image, then call
+:meth:`ECGImage.run_full_pipeline` followed by
+:meth:`ECGImage.save_signals_as_csv` (or :meth:`ECGImage.save_signals_as_wfdb`).
+
+Module-level helpers
+--------------------
+plot_image        — Quick matplotlib preview of a grayscale array.
+shadow_removal    — Morphological background subtraction.
+line_length       — Euclidean distance between two 2-D points.
+parse_layout_from_folder — Parse rows/cols/layout flags from a folder name.
+"""
+
+import os
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["OMP_NUM_THREADS"]      = "1"
+os.environ["MKL_NUM_THREADS"]      = "1"
+
+import numpy as np
+import cv2
+cv2.setNumThreads(0)   # disable OpenCV's internal thread pool
+
+import matplotlib.pyplot as plt
+from ultralytics import YOLO
+from PIL import Image
+from sklearn.cluster import KMeans
+from sklearn.metrics import silhouette_score, mean_squared_error
+from scipy.interpolate import interp1d
+import wfdb
+from scipy import signal
+from skimage import morphology, segmentation
+from scipy.signal import savgol_filter, find_peaks
+from scipy.stats import pearsonr
+from skimage.filters import threshold_multiotsu
+from concurrent.futures import ThreadPoolExecutor
+import re
+import pandas as pd
+import torch
+import time
+
+
+from patched_yolo_infer import (
+    MakeCropsDetectThem,
+    CombineDetections,
+    visualize_results,
+)
+
+
+def plot_image(img, title="Image Plot", size=(12, 12), show_axis=False):
+    """Display a grayscale image with matplotlib.
+
+    Parameters
+    ----------
+    img : np.ndarray
+        Grayscale image array (H × W).
+    title : str, optional
+        Figure title (default ``"Image Plot"``).
+    size : tuple[int, int], optional
+        Figure size in inches ``(width, height)`` (default ``(12, 12)``).
+    show_axis : bool, optional
+        Whether to draw axis ticks (default ``False``).
+    """
+    plt.figure(figsize=size)
+    plt.imshow(img, cmap='gray')
+    plt.title(title)
+    if not show_axis:
+        plt.axis('off')
+    plt.show()
+
+
+def shadow_removal(img):
+    """Remove uneven illumination / shadow from a grayscale image.
+
+    Uses morphological dilation followed by median blur to estimate the
+    background, then subtracts it from the original to yield a
+    normalised, shadow-free image.
+
+    Parameters
+    ----------
+    img : np.ndarray
+        Grayscale uint8 image.
+
+    Returns
+    -------
+    np.ndarray
+        Shadow-corrected uint8 image, intensity range [0, 255].
+    """
+    dilated_img = cv2.dilate(img, np.ones((7, 7), np.uint8))
+    bg_img = cv2.medianBlur(dilated_img, 15)
+    diff_img = 255 - cv2.absdiff(img, bg_img)
+    norm_img = cv2.normalize(diff_img, None, alpha=0, beta=255,
+                             norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8UC1)
+    return norm_img
+
+
+def line_length(x1, y1, x2, y2):
+    """Return the Euclidean length of a line segment.
+
+    Parameters
+    ----------
+    x1, y1 : float
+        Coordinates of the first endpoint.
+    x2, y2 : float
+        Coordinates of the second endpoint.
+
+    Returns
+    -------
+    float
+        Length of the segment in the same units as the inputs.
+    """
+    return np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
+
+
+def parse_layout_from_folder(folder_path):
+    """Parse ECG layout metadata encoded in a folder name.
+
+    Expected naming convention::
+
+        ecg_signals_<rows>x<cols>_<Rythm|None>_<Cabrera|Normal>
+
+    Parameters
+    ----------
+    folder_path : str
+        Absolute or relative path whose final component encodes the layout.
+
+    Returns
+    -------
+    layout_key : tuple[int, int, bool] or None
+        ``(rows, cols, is_cabrera)``.  ``None`` if the name does not match.
+    calibration : bool or None
+        ``True`` when a rhythm (calibration) strip is present, ``None`` on
+        parse failure.
+    """
+    base_name = folder_path.split('/')[-1]
+    match = re.search(r'ecg_signals_(\d+)x(\d+)_(None|Rythm)_(Cabrera|Normal)', base_name)
+    if match:
+        rows = int(match.group(1))
+        cols = int(match.group(2))
+        calibration = (match.group(3) == 'Rythm')
+        cabrera_flag = (match.group(4) == 'Cabrera')
+        layout_key = (rows, cols, cabrera_flag)
+        return layout_key, calibration
+    else:
+        return None, None
+
+
+class ECGImage:
+    """End-to-end digitizer for a single 12-lead ECG image.
+
+    The class encapsulates every processing stage — from raw pixel
+    loading to per-lead voltage time-series export.  All intermediate
+    results (masks, bounding boxes, calibration constants, signal grid)
+    are stored as instance attributes so that individual stages can be
+    inspected or visualised after the pipeline runs.
+
+    Typical usage
+    -------------
+    >>> ecg = ECGImage(box_model, seg_model, name_model, pulse_model, "ecg.png")
+    >>> ecg.run_full_pipeline()
+    >>> ecg.save_signals_as_csv("record_001", directory="./output")
+
+    Parameters
+    ----------
+    box_model : ultralytics.YOLO
+        Pre-loaded YOLO model that detects lead bounding boxes.
+    segmentation_model : str
+        Path to the YOLO segmentation model checkpoint.  Loaded
+        internally by ``patched_yolo_infer.MakeCropsDetectThem`` so
+        that it can operate on image patches.
+    lead_name_model : ultralytics.YOLO
+        Pre-loaded YOLO model that classifies lead name labels
+        (I, II, III, aVR, aVL, aVF, V1–V6).
+    pulse_model : ultralytics.YOLO
+        Pre-loaded YOLO model that detects calibration pulse boxes.
+    image_path : str
+        Path to the ECG image (.png, .jpg, or .jpeg).
+    wfdb_path : str, optional
+        Path to a WFDB record (without extension) used for ground-truth
+        comparison in :meth:`calculate_metrics_ptb`.  Not required for
+        normal digitization.
+
+    Attributes set by the pipeline
+    --------------------------------
+    image : np.ndarray
+        Loaded and padded grayscale image.
+    processed_image : np.ndarray
+        Shadow-removed, blurred version used for all model inference.
+    lead_segmentation : list
+        ``CombineDetections`` objects from the three segmentation passes.
+    mask_image : np.ndarray
+        Binary mask (H × W, uint8) combining all segmented lead polygons.
+    row_centers : np.ndarray
+        Y-pixel positions of detected lead row centres.
+    roi : tuple[float, float]
+        ``(min_y, max_y)`` bounding the active ECG area.
+    lead_bboxes : list[list[float]]
+        ``[x1, y1, x2, y2]`` boxes from ``box_model``.
+    lead_name_bboxes : list[dict]
+        ``{'bbox': [...], 'class_name': str}`` from ``lead_name_model``.
+    reference_pulses : list[dict]
+        ``{'bbox': [...], 'image': np.ndarray}`` from ``pulse_model``.
+    volt_per_pixel : float
+        Calibration constant: millivolts per pixel (vertical axis).
+    time_per_pixel : float
+        Calibration constant: seconds per pixel (horizontal axis).
+    is_cabrera : bool
+        Whether the ECG uses Cabrera lead ordering.
+    has_calibration_pulse : bool
+        Whether a rhythm/calibration strip row is present.
+    layout : list[list[str]]
+        2-D grid of lead names matching the detected row/column layout.
+    grid : list[list[dict]]
+        Mask grid — each cell contains ``'lead'`` (str) and ``'signal'``
+        (np.ndarray of the binary mask slice).
+    signal_grid : list[list[dict]]
+        After :meth:`extract_signals`: cells additionally contain
+        ``'signal'`` as a list of float mV values.
+    """
+
+    def __init__(self, box_model, segmentation_model, lead_name_model, pulse_model,
+                 image_path, wfdb_path=""):
+        self.image_path = image_path          # stored for retry in run_full_pipeline
+        self.load_image(image_path)
+        self.wfdb_path = wfdb_path
+
+        self.box_model = box_model
+        self.segmentation_model = segmentation_model
+        self.lead_name_model = lead_name_model
+        self.pulse_model = pulse_model
+
+        self.is_cabrera = None
+        self.has_calibration_pulse = None
+
+        self.standard_layouts = {
+            # Standard layouts
+            (12, 1, False): [['I'], ['II'], ['III'], ['aVR'], ['aVL'], ['aVF'],
+                             ['V1'], ['V2'], ['V3'], ['V4'], ['V5'], ['V6']],
+            (6, 2, False):  [['I', 'V1'], ['II', 'V2'], ['III', 'V3'],
+                             ['aVR', 'V4'], ['aVL', 'V5'], ['aVF', 'V6']],
+            (4, 3, False):  [['I', 'II', 'III'], ['aVR', 'aVL', 'aVF'],
+                             ['V1', 'V2', 'V3'], ['V4', 'V5', 'V6']],
+            (3, 4, False):  [['I', 'aVR', 'V1', 'V4'], ['II', 'aVL', 'V2', 'V5'],
+                             ['III', 'aVF', 'V3', 'V6']],
+            # Cabrera layouts
+            (12, 1, True):  [['aVL'], ['I'], ['aVR'], ['II'], ['aVF'], ['III'],
+                             ['V1'], ['V2'], ['V3'], ['V4'], ['V5'], ['V6']],
+            (6, 2, True):   [['aVL', 'V1'], ['I', 'V2'], ['aVR', 'V3'],
+                             ['II', 'V4'], ['aVF', 'V5'], ['III', 'V6']],
+            (4, 3, True):   [['aVL', 'I', 'aVR'], ['II', 'aVF', 'III'],
+                             ['V1', 'V2', 'V3'], ['V4', 'V5', 'V6']],
+            (3, 4, True):   [['aVL', 'II', 'V1', 'V4'], ['I', 'aVF', 'V2', 'V5'],
+                             ['aVR', 'III', 'V3', 'V6']],
+        }
+
+    # ------------------------------------------------------------------
+    # Image loading & preprocessing
+    # ------------------------------------------------------------------
+
+    def load_image(self, path, target_size=2100):
+        """Load, resize, and pad an ECG image.
+
+        The image is scaled so that its height equals *target_size* (cubic
+        interpolation) and then padded with a 20-pixel white border on all
+        sides.  The result is stored in ``self.image``.
+
+        Parameters
+        ----------
+        path : str
+            Path to the image file.
+        target_size : int, optional
+            Target image height in pixels (default ``2100``).
+
+        Raises
+        ------
+        FileNotFoundError
+            If OpenCV cannot open the file at *path*.
+        """
+        img = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
+        if img is None:
+            raise FileNotFoundError(f"Cannot read image: {path}")
+        resample_factor = target_size / img.shape[0]
+        img = cv2.resize(img,
+                         (int(img.shape[1] * resample_factor),
+                          int(img.shape[0] * resample_factor)),
+                         interpolation=cv2.INTER_CUBIC)
+        self.image = cv2.copyMakeBorder(img, 20, 20, 20, 20,
+                                        cv2.BORDER_CONSTANT, value=[255, 255, 255])
+
+    def preprocess_image(self):
+        """Apply shadow removal and Gaussian smoothing to ``self.image``.
+
+        The result is stored in ``self.processed_image`` and is used as
+        input for all subsequent YOLO inference calls.
+        """
+        rem = shadow_removal(self.image)
+        self.processed_image = cv2.GaussianBlur(rem, (3, 3), 0)
+
+    # ------------------------------------------------------------------
+    # Segmentation
+    # ------------------------------------------------------------------
+
+    def segment_leads(self):
+        """Run patched YOLO instance segmentation at three crop scales.
+
+        The segmentation model is applied to overlapping image patches at
+        crop factors of 4, 4.5, and 5 (relative to image height).
+        Detections from each pass are combined with NMS
+        (``nms_threshold=0.5``) and stored in ``self.lead_segmentation``
+        as a list of three ``CombineDetections`` objects.
+        """
+        segmentations = []
+        for shape in [4, 4.5, 5]:
+            element_crops = MakeCropsDetectThem(
+                image=cv2.cvtColor(self.processed_image, cv2.COLOR_GRAY2BGR),
+                model_path=self.segmentation_model,
+                segment=True,
+                show_crops=False,
+                shape_x=int(self.processed_image.shape[0] // shape),
+                shape_y=int(self.processed_image.shape[0] // shape),
+                overlap_x=50,
+                overlap_y=50,
+                conf=0.8,
+                iou=0.7,
+                classes_list=[0],
+            )
+            segmentations.append(CombineDetections(element_crops, nms_threshold=0.5))
+        self.lead_segmentation = segmentations
+
+    def make_segmentation_mask(self):
+        """Rasterise all segmented polygons into a single binary mask.
+
+        Polygons from all three segmentation passes in
+        ``self.lead_segmentation`` are filled and merged into a combined
+        mask.  A morphological opening (5 × 5 kernel) removes small
+        spurious regions.  The result is stored in ``self.mask_image``
+        (uint8, values 0 or 255).
+        """
+        height, width = self.image.shape[:2]
+        combined_mask = np.zeros((height, width), dtype=np.uint8)
+
+        for segmentation in self.lead_segmentation:
+            for poly in segmentation.filtered_polygons:
+                pts = np.array(poly, dtype=np.int32)
+                if pts.ndim == 2:
+                    pts = [pts]
+                cv2.fillPoly(combined_mask, pts, color=255)
+
+        kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5))
+        self.mask_image = cv2.morphologyEx(combined_mask, cv2.MORPH_OPEN, kernel)
+
+    # ------------------------------------------------------------------
+    # Row detection
+    # ------------------------------------------------------------------
+
+    def find_row_centers(self):
+        """Detect the vertical centre of each ECG lead row.
+
+        Projects the binary mask onto the vertical axis (row sums) and
+        applies ``scipy.signal.find_peaks`` twice: once to identify all
+        peaks, and a second time with an adaptive minimum distance
+        (⅔ × mean inter-peak spacing) to consolidate merged rows.
+
+        Sets
+        ----
+        self.row_centers : np.ndarray
+            Y-pixel indices of each row centre.
+        self.first_peak_start : int
+            Y-pixel where the topmost row's signal begins.
+        self.last_peak_end : int
+            Y-pixel where the bottommost row's signal ends.
+        """
+        image_height, image_width = self.mask_image.shape[:2]
+        proj = np.sum(self.mask_image, 1)
+
+        height   = (image_width // 10) * 255
+        distance = image_height // 30
+
+        peaks, _ = find_peaks(proj, height=height, distance=distance)
+        row_centers, _ = find_peaks(
+            proj, height=height,
+            distance=int(np.mean(np.diff(peaks)) * (2 / 3))
+        )
+        self.row_centers = row_centers
+
+        if len(peaks) > 0:
+            first_peak  = peaks[0]
+            start_index = first_peak
+            zero_gap    = 0
+            for i in range(first_peak - 1, -1, -1):
+                if proj[i] == 0:
+                    zero_gap += 1
+                    if zero_gap > 2:
+                        break
+                else:
+                    zero_gap    = 0
+                    start_index = i
+            self.first_peak_start = start_index
+
+            last_peak = peaks[-1]
+            end_index = last_peak
+            zero_gap  = 0
+            for i in range(last_peak + 1, image_height):
+                if proj[i] == 0:
+                    zero_gap += 1
+                    if zero_gap > 2:
+                        break
+                else:
+                    zero_gap  = 0
+                    end_index = i
+            self.last_peak_end = end_index
+
+    def get_roi(self):
+        """Compute the vertical region-of-interest (ROI) spanning all lead rows.
+
+        Extends ⅔ of the mean row spacing above the first row centre and
+        below the last.  Stores the result as ``self.roi = (min_y, max_y)``.
+        """
+        spacing = 2 / 3 * np.mean(np.diff(self.row_centers))
+        min_y = max(0, self.row_centers[0] - spacing)
+        max_y = min(self.image.shape[0], self.row_centers[-1] + spacing)
+        self.roi = (min_y, max_y)
+
+    # ------------------------------------------------------------------
+    # YOLO inference helpers — GPU-safe wrappers
+    # ------------------------------------------------------------------
+
+    def _predict_safe(self, model, image_bgr, **kwargs):
+        """Run a YOLO model prediction in a thread-safe manner.
+
+        On CUDA: wraps the call in a dedicated ``torch.cuda.Stream`` and
+        synchronises before returning, preventing contention on the
+        default stream when multiple threads call this simultaneously.
+        On CPU: wraps with ``torch.no_grad()`` only.
+
+        Parameters
+        ----------
+        model : ultralytics.YOLO
+            The detection/classification model to run.
+        image_bgr : np.ndarray
+            BGR image array as expected by Ultralytics.
+        **kwargs
+            Forwarded to ``model.predict()`` (e.g. ``conf``, ``iou``).
+
+        Returns
+        -------
+        list
+            Ultralytics ``Results`` objects.
+        """
+        if torch.cuda.is_available():
+            stream = torch.cuda.Stream()
+            with torch.no_grad():
+                with torch.cuda.stream(stream):
+                    results = model.predict(image_bgr, verbose=False, **kwargs)
+                stream.synchronize()
+        else:
+            with torch.no_grad():
+                results = model.predict(image_bgr, verbose=False, **kwargs)
+        return results
+
+    # ------------------------------------------------------------------
+    # Timing helper
+    # ------------------------------------------------------------------
+
+    @staticmethod
+    def _timed(label, fn, *args, **kwargs):
+        """Call *fn* with *args*/*kwargs*, print its wall-clock duration, and return its result.
+
+        Parameters
+        ----------
+        label : str
+            Human-readable stage name printed alongside the elapsed time.
+        fn : callable
+            The stage function to time.
+        *args, **kwargs
+            Forwarded to *fn*.
+
+        Returns
+        -------
+        Any
+            Whatever *fn* returns.
+        """
+        t = time.time()
+        result = fn(*args, **kwargs)
+        print(f"    ⏱ {label}: {time.time() - t:.1f}s", flush=True)
+        return result
+
+    # ------------------------------------------------------------------
+    # Box & name extraction
+    # ------------------------------------------------------------------
+
+    def extract_lead_boxes(self):
+        """Detect lead bounding boxes with ``box_model`` and filter to the ROI.
+
+        Runs ``_predict_safe`` on the preprocessed BGR image at confidence
+        0.8, then discards boxes whose ``y1``/``y2`` coordinates fall
+        outside ``self.roi``.  Sets ``self.lead_bboxes`` as a list of
+        ``[x1, y1, x2, y2]`` float lists.
+        """
+        image_bgr = cv2.cvtColor(self.processed_image, cv2.COLOR_GRAY2BGR)
+        results   = self._predict_safe(self.box_model, image_bgr, conf=0.8)
+
+        min_y, max_y = self.roi
+        lead_boxes   = []
+        for r in results:
+            for box in r.boxes:
+                x1, y1, x2, y2 = box.xyxy.cpu().numpy()[0].tolist()
+                if y1 >= min_y and y2 <= max_y:
+                    lead_boxes.append([x1, y1, x2, y2])
+        self.lead_bboxes = lead_boxes
+
+    def extract_lead_name_boxes(self):
+        """Detect and classify lead name labels with ``lead_name_model``.
+
+        Runs at confidence 0.8, filters to the ROI, and stores results in
+        ``self.lead_name_bboxes`` as a list of
+        ``{'bbox': [x1, y1, x2, y2], 'class_name': str}`` dicts.
+        """
+        image_bgr = cv2.cvtColor(self.processed_image, cv2.COLOR_GRAY2BGR)
+        results   = self._predict_safe(self.lead_name_model, image_bgr, conf=0.8)
+
+        min_y, max_y = self.roi
+        name_boxes   = []
+        for r in results:
+            for box in r.boxes:
+                x1, y1, x2, y2 = box.xyxy.cpu().numpy()[0].tolist()
+                if y1 >= min_y and y2 <= max_y:
+                    cls_id   = int(box.cls.cpu().numpy()[0])
+                    cls_name = self.lead_name_model.names[cls_id]
+                    name_boxes.append({'bbox': [x1, y1, x2, y2], 'class_name': cls_name})
+        self.lead_name_bboxes = name_boxes
+
+    def extract_reference_pulses(self):
+        """Detect calibration pulse boxes with ``pulse_model``.
+
+        Runs at confidence 0.7 (lower than lead detection to improve
+        recall on small pulses).  For each detection, stores the
+        corresponding image crop alongside the bounding-box coordinates.
+        Sets ``self.reference_pulses`` as a list of
+        ``{'bbox': [x1, y1, x2, y2], 'image': np.ndarray}`` dicts.
+        """
+        image_bgr  = cv2.cvtColor(self.processed_image, cv2.COLOR_GRAY2BGR)
+        results    = self._predict_safe(self.pulse_model, image_bgr, conf=0.7)
+
+        pulse_boxes = []
+        for r in results:
+            for box in r.boxes:
+                coord = box.xyxy.cpu().numpy()[0].tolist()
+                pulse_boxes.append({
+                    'image': self.image[int(coord[1]) - 5:int(coord[3]) + 5,
+                                        int(coord[0]) - 5:int(coord[2]) + 5],
+                    'bbox': coord,
+                })
+        self.reference_pulses = pulse_boxes
+
+    # ------------------------------------------------------------------
+    # Visualisation helpers
+    # ------------------------------------------------------------------
+
+    def visualize_boxes(self, task='Lead name', show_axis=False):
+        """Overlay detection bounding boxes on the original image.
+
+        Parameters
+        ----------
+        task : {'Lead name', 'Lead box', 'Reference pulse'}
+            Which set of boxes to draw.
+        show_axis : bool, optional
+            Whether to display axis ticks (default ``False``).
+        """
+        img_copy = cv2.cvtColor(self.image.copy(), cv2.COLOR_GRAY2BGR)
+
+        if task == 'Lead name':
+            if not self.lead_name_bboxes:
+                print("Lead name boxes not extracted.")
+                return
+            for box in self.lead_name_bboxes:
+                x1, y1, x2, y2 = map(int, box['bbox'])
+                cv2.rectangle(img_copy, (x1, y1), (x2, y2), (255, 0, 0), 2)
+                cv2.putText(img_copy, box['class_name'], (x1, y1 - 10),
+                            cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 0, 0), 3)
+
+        elif task == 'Lead box':
+            if not self.lead_bboxes:
+                print("Lead boxes not extracted.")
+                return
+            for bbox in self.lead_bboxes:
+                x1, y1, x2, y2 = map(int, bbox)
+                cv2.rectangle(img_copy, (x1, y1), (x2, y2), (255, 0, 0), 2)
+
+        elif task == 'Reference pulse':
+            if not self.reference_pulses:
+                print("Reference pulses not extracted.")
+                return
+            for bbox in self.reference_pulses:
+                x1, y1, x2, y2 = map(int, bbox['bbox'])
+                cv2.rectangle(img_copy, (x1, y1), (x2, y2), (255, 0, 0), 2)
+
+        else:
+            print(f"Unknown task: {task}")
+            return
+
+        plt.figure(figsize=(12, 10))
+        plt.imshow(cv2.cvtColor(img_copy, cv2.COLOR_BGR2RGB))
+        if not show_axis:
+            plt.axis('off')
+        plt.show()
+
+    def visualize_segmentation(self, show_boxes=False, show_axis=False,
+                               fill_mask=True, thickness=1):
+        """Render all segmented lead polygons on the original image.
+
+        Aggregates results from all three segmentation passes in
+        ``self.lead_segmentation`` and passes them to
+        ``patched_yolo_infer.visualize_results``.
+
+        Parameters
+        ----------
+        show_boxes : bool, optional
+            Also draw bounding boxes around each polygon (default ``False``).
+        show_axis : bool, optional
+            Display axis ticks (default ``False``).
+        fill_mask : bool, optional
+            Fill polygon interiors with a translucent colour (default ``True``).
+        thickness : int, optional
+            Polygon outline thickness in pixels (default ``1``).
+        """
+        all_confidences, all_boxes, all_polygons = [], [], []
+        all_classes_ids, all_classes_names       = [], []
+
+        for seg in self.lead_segmentation:
+            all_confidences.extend(seg.filtered_confidences)
+            all_boxes.extend(seg.filtered_boxes)
+            all_polygons.extend(seg.filtered_polygons)
+            all_classes_ids.extend(seg.filtered_classes_id)
+            all_classes_names.extend(seg.filtered_classes_names)
+
+        visualize_results(
+            img=cv2.cvtColor(self.image, cv2.COLOR_GRAY2RGB),
+            confidences=all_confidences,
+            boxes=all_boxes,
+            polygons=all_polygons,
+            classes_ids=all_classes_ids,
+            classes_names=all_classes_names,
+            segment=True,
+            thickness=thickness,
+            fill_mask=fill_mask,
+            show_boxes=show_boxes,
+            show_class=False,
+            axis_off=(not show_axis),
+        )
+
+    # ------------------------------------------------------------------
+    # Calibration (volt/pixel & time/pixel)
+    # ------------------------------------------------------------------
+
+    def get_reference_scale(self):
+        """Derive pixel-to-physical calibration constants from the calibration pulses.
+
+        For each pulse in ``self.reference_pulses``, the method:
+
+        1. Enhances contrast with CLAHE and binarises with multi-Otsu thresholding.
+        2. Detects horizontal and vertical lines via morphological filtering and
+           Hough transform to measure the pulse amplitude (voltage) and width (time).
+        3. Refines the time calibration using the Line Segment Detector (LSD).
+
+        Sets
+        ----
+        self.volt_per_pixel : float
+            Millivolts per pixel (vertical axis).  Calibrated to the standard
+            1 mV / 10 mm calibration pulse.
+        self.time_per_pixel : float
+            Seconds per pixel (horizontal axis).  Calibrated to 0.2 s between
+            the two inner vertical edges of the calibration pulse.
+
+        Raises
+        ------
+        RuntimeError
+            If no valid pulse yields voltage or time measurements.
+        """
+        voltages, times, dist = [], [], []
+
+        def _line_length(x1, y1, x2, y2):
+            return np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
+
+        for pulse_idx, pulse in enumerate(self.reference_pulses):
+            try:
+                img = pulse['image']
+                if img is None or img.size == 0:
+                    continue
+                h, w = img.shape
+                if h < 10 or w < 10:
+                    continue
+
+                clahe    = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
+                enhanced = clahe.apply(img)
+                blurred  = cv2.GaussianBlur(enhanced, (3, 3), 0)
+
+                thresholds = threshold_multiotsu(blurred, classes=2)
+                regions    = np.digitize(blurred, bins=thresholds)
+                binary     = (regions == 0).astype(np.uint8) * 255
+
+                h_kernel    = cv2.getStructuringElement(cv2.MORPH_RECT, (25, 1))
+                v_kernel    = cv2.getStructuringElement(cv2.MORPH_RECT, (1, 25))
+                h_lines_img = cv2.morphologyEx(binary, cv2.MORPH_OPEN, h_kernel)
+                v_lines_img = cv2.morphologyEx(binary, cv2.MORPH_OPEN, v_kernel)
+
+                h_lines_raw = cv2.HoughLinesP(
+                    h_lines_img, rho=1, theta=np.pi / 180,
+                    threshold=w // 4, minLineLength=w // 2, maxLineGap=1
+                )
+                v_lines_raw = cv2.HoughLinesP(
+                    v_lines_img, rho=1, theta=np.pi / 180,
+                    threshold=h // 4, minLineLength=h // 2, maxLineGap=1
+                )
+
+                combined_lines = []
+                if h_lines_raw is not None:
+                    combined_lines.extend(h_lines_raw)
+                if v_lines_raw is not None:
+                    combined_lines.extend(v_lines_raw)
+
+                horizontal_lines, vertical_lines = [], []
+                for line in combined_lines:
+                    x1, y1, x2, y2 = line[0]
+                    angle = np.arctan2(abs(y2 - y1), abs(x2 - x1)) * 180 / np.pi
+                    if angle < 10 or angle > 170:
+                        horizontal_lines.append((x1, y1, x2, y2))
+                    elif 80 < angle < 100:
+                        vertical_lines.append((x1, y1, x2, y2))
+
+                if len(vertical_lines) >= 2:
+                    longest_v = sorted(vertical_lines,
+                                       key=lambda l: _line_length(*l), reverse=True)[:2]
+                    xc1 = (longest_v[0][0] + longest_v[0][2]) / 2
+                    xc2 = (longest_v[1][0] + longest_v[1][2]) / 2
+                    vertical_spacing = abs(xc1 - xc2)
+
+                    volt_vals = sorted([_line_length(*l) for l in vertical_lines],
+                                       reverse=True)[:2]
+                    if volt_vals:
+                        voltages.append(volt_vals)
+                    if vertical_spacing > 5:
+                        times.append(vertical_spacing)
+
+                # LSD for more precise time calibration
+                lsd = cv2.createLineSegmentDetector(refine=2)
+                lsd_lines, _, _, _ = lsd.detect(v_lines_img)
+
+                min_length    = h / 3
+                angle_tol_rad = np.deg2rad(5)
+                filtered      = []
+
+                if lsd_lines is not None:
+                    for line in lsd_lines:
+                        x1, y1, x2, y2 = line[0]
+                        length = np.hypot(x2 - x1, y2 - y1)
+                        if length >= min_length:
+                            angle = np.arctan2(y2 - y1, x2 - x1)
+                            if np.abs(np.abs(angle) - np.pi / 2) <= angle_tol_rad:
+                                filtered.append([[x1, y1, x2, y2]])
+
+                if len(filtered) >= 4:
+                    filtered.sort(
+                        key=lambda l: np.hypot(l[0][2] - l[0][0], l[0][3] - l[0][1]),
+                        reverse=True
+                    )
+                    top4 = filtered[:4]
+                    top4.sort(key=lambda l: (l[0][0] + l[0][2]) / 2)
+
+                    midpoints = []
+                    for l1, l2 in [(top4[0], top4[1]), (top4[2], top4[3])]:
+                        xc1 = (l1[0][0] + l1[0][2]) / 2
+                        xc2 = (l2[0][0] + l2[0][2]) / 2
+                        midpoints.append((xc1 + xc2) / 2)
+
+                    dist.append(midpoints[1] - midpoints[0])
+
+            except Exception as e:
+                print(f"    ⚠ Skipping pulse {pulse_idx}: {e}", flush=True)
+                continue
+
+        if not voltages:
+            raise RuntimeError("No valid reference pulses found for voltage calibration.")
+        if not dist:
+            raise RuntimeError("No valid reference pulses found for time calibration.")
+
+        self.volt_per_pixel = 1 / np.mean(voltages)
+        self.time_per_pixel = 0.2 / np.mean(dist)
+
+    # ------------------------------------------------------------------
+    # Layout helpers
+    # ------------------------------------------------------------------
+
+    def make_bounding_box_features(self, box, axis):
+        """Extract three positional features from a bounding box along one axis.
+
+        Parameters
+        ----------
+        box : list[float]
+            Bounding box as ``[x1, y1, x2, y2]``.
+        axis : {'x', 'y'}
+            Axis along which to compute features.
+
+        Returns
+        -------
+        list[float]
+            ``[axis_min, axis_center, axis_max]``.
+        """
+        if axis == 'y':
+            axis_min, axis_max = box[1], box[3]
+        else:
+            axis_min, axis_max = box[0], box[2]
+        axis_center = (axis_min + axis_max) / 2
+        return [axis_min, axis_center, axis_max]
+
+    def bounding_boxes_kmeans(self, bounding_boxes, axis='y',
+                              k_min=1, k_max=13, return_model=True):
+        """Cluster bounding boxes along one axis using K-Means.
+
+        Tries all values of *k* in ``[k_min, k_max]`` and selects the one
+        that maximises the silhouette score.  Cluster labels are sorted by
+        mean position along *axis* (ascending).
+
+        Parameters
+        ----------
+        bounding_boxes : list[list[float]]
+            List of ``[x1, y1, x2, y2]`` boxes.
+        axis : {'x', 'y'}, optional
+            Axis along which to cluster (default ``'y'``).
+        k_min : int, optional
+            Minimum number of clusters (default ``1``).
+        k_max : int, optional
+            Maximum number of clusters (default ``13``).
+        return_model : bool, optional
+            If ``True``, also return the fitted ``KMeans`` model and the
+            label remapping dict (default ``True``).
+
+        Returns
+        -------
+        sorted_labels : np.ndarray of int
+            Cluster index (sorted) for each input box.
+        best_k : int
+            Optimal number of clusters.
+        sorted_centers : np.ndarray
+            Mean positional value of each cluster, sorted ascending.
+        label_map : dict
+            Mapping from original K-Means label → sorted label.
+        best_model : sklearn.cluster.KMeans
+            Fitted model (only when ``return_model=True``).
+
+        Raises
+        ------
+        ValueError
+            If *axis* is not ``'x'`` or ``'y'``, or if there are fewer
+            boxes than *k_min*.
+        """
+        if axis not in ('x', 'y'):
+            raise ValueError("Axis must be 'x' or 'y'")
+        if len(bounding_boxes) < k_min:
+            raise ValueError("Not enough bounding boxes to cluster")
+
+        features = np.array([self.make_bounding_box_features(b, axis)
+                             for b in bounding_boxes])
+
+        best_score, best_k = -1, k_min
+        best_labels = best_centers = best_model = None
+
+        for k in range(k_min, min(k_max + 1, len(bounding_boxes))):
+            kmeans = KMeans(n_clusters=k, random_state=42, n_init="auto")
+            labels = kmeans.fit_predict(features)
+            score  = silhouette_score(features, labels)
+            if score > best_score:
+                best_score   = score
+                best_k       = k
+                best_labels  = labels
+                best_centers = kmeans.cluster_centers_
+                best_model   = kmeans
+
+        cluster_avgs   = best_centers.mean(axis=1)
+        sorted_indices = np.argsort(cluster_avgs)
+        label_map      = {old: new for new, old in enumerate(sorted_indices)}
+        sorted_labels  = np.array([label_map[l] for l in best_labels])
+        sorted_centers = cluster_avgs[sorted_indices]
+
+        if return_model:
+            return sorted_labels, best_k, sorted_centers, label_map, best_model
+        return sorted_labels, best_k, sorted_centers, label_map
+
+    def check_cabrera(self, num_rows, num_cols):
+        """Determine whether the ECG uses Cabrera lead ordering.
+
+        The heuristic varies by layout:
+
+        * **12- or 6-row layouts:** compares the vertical spacing of
+          augmented-limb leads (aVR/aVL/aVF) against precordial leads.
+          Cabrera re-orders the limb leads so their spacing differs.
+        * **4-row, 3-col / 5-row layouts:** checks the vertical spread
+          of augmented leads (high std → interleaved → Cabrera).
+        * **3-row / 4-col layouts:** checks the horizontal spread of
+          augmented leads.
+
+        Sets ``self.is_cabrera`` as a side-effect.
+
+        Parameters
+        ----------
+        num_rows : int
+            Number of detected lead rows.
+        num_cols : int
+            Number of detected lead columns.
+
+        Returns
+        -------
+        bool
+            ``True`` if Cabrera ordering is detected.
+        """
+        if num_rows in [13, 12, 7, 6]:
+            av_leads = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'aVR', 'aVL', 'aVF'}]
+            v_leads  = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'V1', 'V2', 'V3', 'V4', 'V5', 'V6'}]
+
+            if not av_leads or not v_leads:
+                return False
+
+            y_v  = sorted([(b['bbox'][1] + b['bbox'][3]) / 2 for b in v_leads])
+            y_av = sorted([(b['bbox'][1] + b['bbox'][3]) / 2 for b in av_leads])
+
+            threshold = 30
+            diff_v    = np.diff(y_v)
+            diff_av   = np.diff(y_av)
+            diff_v    = diff_v[np.abs(diff_v) > threshold]
+            diff_av   = diff_av[np.abs(diff_av) > threshold]
+
+            if len(diff_v) == 0 or len(diff_av) == 0:
+                return False
+
+            if abs(np.min(diff_v) - np.min(diff_av)) > (0.25 * np.min(diff_v)):
+                self.is_cabrera = True
+                return True
+            self.is_cabrera = False
+            return False
+
+        elif (num_rows == 4 and num_cols == 3) or num_rows == 5:
+            av_leads = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'aVR', 'aVL', 'aVF'}]
+            if not av_leads:
+                return False
+            y_coords = [(b['bbox'][1] + b['bbox'][3]) / 2 for b in av_leads]
+            self.is_cabrera = np.std(y_coords) > 25
+            return self.is_cabrera
+
+        elif (num_rows == 4 and num_cols == 4) or num_rows == 3:
+            av_leads = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'aVR', 'aVL', 'aVF'}]
+            if not av_leads:
+                return False
+            x_coords = [(b['bbox'][0] + b['bbox'][2]) / 2 for b in av_leads]
+            self.is_cabrera = np.std(x_coords) > 25
+            return self.is_cabrera
+
+    def get_layout(self, num_rows):
+        """Map the detected row count to a standard ECG layout.
+
+        Calls :meth:`check_cabrera` to determine lead ordering, sets
+        ``self.layout`` (2-D list of lead name strings) and
+        ``self.has_calibration_pulse``, and returns the number of columns.
+
+        Supported row counts
+        --------------------
+        * 13 → 12×1 with calibration pulse
+        * 12 → 12×1 without calibration pulse
+        * 7  → 6×2 with calibration pulse
+        * 6  → 6×2 without calibration pulse
+        * 5  → 4×3 with calibration pulse
+        * 4  → 3×4 or 4×3 (disambiguated by V-lead spatial distribution)
+        * 3  → 3×4 without calibration pulse
+
+        Parameters
+        ----------
+        num_rows : int
+            Number of detected lead rows.
+
+        Returns
+        -------
+        int
+            Number of lead columns in the detected layout.
+        """
+        if num_rows == 13:
+            num_cols = 1
+            self.has_calibration_pulse = True
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(12, 1, cabrera)]
+
+        elif num_rows == 12:
+            num_cols = 1
+            self.has_calibration_pulse = False
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(12, 1, cabrera)]
+
+        elif num_rows == 7:
+            num_cols = 2
+            self.has_calibration_pulse = True
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(6, 2, cabrera)]
+
+        elif num_rows == 6:
+            num_cols = 2
+            self.has_calibration_pulse = False
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(6, 2, cabrera)]
+
+        elif num_rows == 5:
+            num_cols = 3
+            self.has_calibration_pulse = True
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(4, 3, cabrera)]
+
+        elif num_rows == 4:
+            v_leads1 = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'V1', 'V2', 'V3'}]
+            v_leads2 = [b for b in self.lead_name_bboxes
+                        if b['class_name'] in {'V4', 'V5', 'V6'}]
+
+            def centers(boxes, axis):
+                if axis == 'x':
+                    return [(b['bbox'][0] + b['bbox'][2]) / 2 for b in boxes]
+                return [(b['bbox'][1] + b['bbox'][3]) / 2 for b in boxes]
+
+            x_std1 = np.std(centers(v_leads1, 'x'))
+            y_std1 = np.std(centers(v_leads1, 'y'))
+            x_std2 = np.std(centers(v_leads2, 'x'))
+            y_std2 = np.std(centers(v_leads2, 'y'))
+
+            if x_std1 < y_std1 and x_std2 < y_std2:
+                num_cols = 4
+                self.has_calibration_pulse = True
+            elif y_std1 < x_std1 and y_std2 < x_std2:
+                num_cols = 3
+                self.has_calibration_pulse = False
+            else:
+                num_cols = 4
+                self.has_calibration_pulse = True
+
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            key = (3, 4, cabrera) if num_cols == 4 else (4, 3, cabrera)
+            self.layout = self.standard_layouts[key]
+
+        elif num_rows == 3:
+            num_cols = 4
+            self.has_calibration_pulse = False
+            cabrera = self.check_cabrera(num_rows, num_cols)
+            self.layout = self.standard_layouts[(3, 4, cabrera)]
+
+        return num_cols
+
+    def clean_lead_mask(self, img, baseline_y):
+        """Remove leaked segments of neighboring leads from a single lead cell binary mask.
+
+        Parameters
+        ----------
+        img : np.ndarray
+            Binary mask for a single lead cell (H × W).
+        baseline_y : float
+            Y-coordinate of the row baseline relative to the cropped row mask.
+
+        Returns
+        -------
+        np.ndarray
+            Cleaned binary mask.
+        """
+        height, width = img.shape
+        if height < 10 or width < 10:
+            return img
+
+        num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(img)
+        if num_labels <= 2:
+            return img
+
+        # Find the primary component (largest horizontal span)
+        primary_c = -1
+        max_w = -1
+        for c in range(1, num_labels):
+            w = stats[c, cv2.CC_STAT_WIDTH]
+            if w > max_w:
+                max_w = w
+                primary_c = c
+
+        if primary_c == -1:
+            return img
+
+        primary_area = stats[primary_c, cv2.CC_STAT_AREA]
+        cleaned_img = img.copy()
+
+        for c in range(1, num_labels):
+            if c == primary_c:
+                continue
+
+            c_top = stats[c, cv2.CC_STAT_TOP]
+            c_height = stats[c, cv2.CC_STAT_HEIGHT]
+            c_bottom = c_top + c_height
+            c_width = stats[c, cv2.CC_STAT_WIDTH]
+            c_area = stats[c, cv2.CC_STAT_AREA]
+            c_centroid_y = centroids[c][1]
+
+            # Heuristics to identify a leak:
+            # 1. It is relatively small (either narrow or small area)
+            is_small = (c_width < 0.35 * width) or (c_area < 0.25 * primary_area)
+
+            # 2. It is located near the vertical boundaries of the cell
+            near_boundary = (c_top < 0.15 * height) or (c_bottom > 0.85 * height)
+
+            # 3. It doesn't cross the baseline (with safety margin)
+            baseline_margin = max(5, int(0.08 * height))
+            crosses_baseline = (c_top - baseline_margin <= baseline_y <= c_bottom + baseline_margin)
+
+            # 4. It is far from the baseline
+            far_from_baseline = abs(c_centroid_y - baseline_y) > 0.2 * height
+
+            if is_small and near_boundary and not crosses_baseline and far_from_baseline:
+                cleaned_img[labels == c] = 0
+
+        return cleaned_img
+
+    # ------------------------------------------------------------------
+    # Grid construction
+    # ------------------------------------------------------------------
+
+    def make_grid(self, padding=0):
+        """Build the lead mask grid from segmentation polygons and lead boxes.
+
+        Steps:
+
+        1. Assigns each segmented polygon to the nearest row centre by
+           polygon centroid Y-coordinate.
+        2. Crops each row mask to its vertical extent.
+        3. Clusters lead bounding boxes into *num_cols* columns with K-Means.
+        4. Slices each (row, col) cell from the corresponding row mask.
+        5. Applies a morphological opening to each cell slice.
+
+        Sets
+        ----
+        self.grid : list[list[dict]]
+            ``grid[row][col]`` contains ``{'lead': str, 'signal': np.ndarray}``,
+            where ``'signal'`` is currently a binary mask slice.
+        self.baseline : list[float]
+            Relative Y-position of the isoelectric baseline within each row
+            (pixels from the top of the cropped row mask).
+        self.row : dict[int, np.ndarray]
+            Cropped binary mask for each row index.
+
+        Parameters
+        ----------
+        padding : int, optional
+            Extra pixel padding (currently unused, reserved for future use).
+        """
+        image_height, image_width = self.image.shape[:2]
+        num_rows = len(self.row_centers)
+        num_cols = self.get_layout(num_rows)
+
+        row_masks    = {i: np.zeros((image_height, image_width), dtype=np.uint8)
+                        for i in range(num_rows)}
+        row_polygons = {i: [] for i in range(num_rows)}
+        row_limits   = {i: [] for i in range(num_rows)}
+
+        for seg in self.lead_segmentation:
+            for box, polygon in zip(seg.filtered_boxes, seg.filtered_polygons):
+                poly_np   = np.array(polygon, dtype=np.int32).reshape((-1, 1, 2))
+                temp_mask = np.zeros((image_height, image_width), dtype=np.uint8)
+                cv2.fillPoly(temp_mask, [poly_np], color=1)
+
+                temp_proj      = np.sum(temp_mask, axis=1)
+                total_weight_y = np.sum(temp_proj)
+                if total_weight_y == 0:
+                    continue
+
+                centroid_y = int(np.sum(np.arange(temp_proj.shape[0]) * temp_proj)
+                                 / total_weight_y)
+
+                y_vals      = [pt[1] for pt in polygon]
+                min_y, max_y = min(y_vals), max(y_vals)
+
+                if max_y < self.first_peak_start or min_y > self.last_peak_end:
+                    continue
+
+                closest_idx = int(np.argmin(np.abs(self.row_centers - centroid_y)))
+
+                if self.has_calibration_pulse and closest_idx == num_rows - 1:
+                    continue
+
+                row_polygons[closest_idx].append(polygon)
+                cv2.fillPoly(row_masks[closest_idx], [poly_np], color=1)
+
+                if not row_limits[closest_idx]:
+                    row_limits[closest_idx] = [min_y, max_y]
+                else:
+                    row_limits[closest_idx][0] = min(row_limits[closest_idx][0], min_y)
+                    row_limits[closest_idx][1] = max(row_limits[closest_idx][1], max_y)
+
+        cropped_row_masks = {}
+        for i in range(num_rows):
+            if self.has_calibration_pulse and i == num_rows - 1:
+                continue
+            if not row_limits[i]:
+                continue
+            min_y, max_y = row_limits[i]
+            cropped_row_masks[i] = row_masks[i][min_y:max_y + 1, :]
+
+        lead_boxes = [
+            box for box in self.lead_bboxes
+            if not (
+                self.has_calibration_pulse and
+                int(np.argmin(np.abs(self.row_centers - (box[1] + box[3]) / 2)))
+                == num_rows - 1
+            )
+        ] if self.has_calibration_pulse else self.lead_bboxes
+
+        self.row = cropped_row_masks
+
+        if num_cols != 1:
+            labels_cols, _, _, _ = self.bounding_boxes_kmeans(
+                lead_boxes, axis='x', k_min=num_cols, k_max=num_cols, return_model=False
+            )
+        else:
+            labels_cols = np.zeros(len(lead_boxes), dtype=int)
+
+        boxes_arr = np.array(lead_boxes)
+        min_x_per_col, max_x_per_col = [], []
+        for col_label in range(num_cols):
+            col_boxes = boxes_arr[labels_cols == col_label]
+            min_x_per_col.append(col_boxes[:, 0].min())
+            max_x_per_col.append(col_boxes[:, 2].max())
+
+        mask_grid, relative_baselines = [], []
+        for row_idx, row_slice in cropped_row_masks.items():
+            row_cells = []
+            for col_idx in range(num_cols):
+                x_min      = max(0, int(min_x_per_col[col_idx]))
+                x_max      = min(image_width, int(max_x_per_col[col_idx]))
+                cell_slice = row_slice[:, x_min:x_max]
+
+                kernel     = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))
+                cell_slice = cv2.morphologyEx(cell_slice, cv2.MORPH_OPEN, kernel)
+
+                # Post-processing: Clean up leaked segments from neighboring leads
+                baseline_y = self.row_centers[row_idx] - row_limits[row_idx][0]
+                cell_slice = self.clean_lead_mask(cell_slice, baseline_y)
+
+                row_cells.append({
+                    'lead':   self.layout[row_idx][col_idx],
+                    'signal': cell_slice,
+                })
+
+            relative_baselines.append(
+                self.row_centers[row_idx] - row_limits[row_idx][0]
+            )
+            mask_grid.append(row_cells)
+
+        self.grid     = mask_grid
+        self.baseline = relative_baselines
+
+    def visualize_grid(self, figsize=(15, 10)):
+        """Display the lead mask grid — one subplot per (row, col) cell.
+
+        Parameters
+        ----------
+        figsize : tuple[int, int], optional
+            Overall figure size in inches (default ``(15, 10)``).
+
+        Raises
+        ------
+        ValueError
+            If :meth:`make_grid` has not been called yet.
+        """
+        if not hasattr(self, 'grid') or not self.grid:
+            raise ValueError("Grid not generated. Call make_grid() first.")
+
+        num_rows = len(self.grid)
+        num_cols = len(self.grid[0]) if isinstance(self.grid[0], list) else 1
+
+        fig, axes = plt.subplots(num_rows, num_cols, figsize=figsize)
+        if num_cols == 1:
+            axes = np.atleast_2d(axes).T
+        elif num_rows == 1:
+            axes = np.atleast_2d(axes)
+
+        for row_idx in range(num_rows):
+            row = self.grid[row_idx] if isinstance(self.grid[row_idx], list) \
+                  else [self.grid[row_idx]]
+            for col_idx, cell in enumerate(row):
+                ax = axes[row_idx][col_idx]
+                ax.imshow(cell['signal'], cmap='gray', aspect='auto')
+                ax.set_title(cell['lead'], fontsize=10)
+                ax.axis('off')
+
+        plt.tight_layout()
+        plt.show()
+
+    # ------------------------------------------------------------------
+    # Signal extraction
+    # ------------------------------------------------------------------
+
+    def binarize_signal(self, img, window_length=11, polyorder=2):
+        """Extract a 1-D signal trace from a binary lead mask column by column.
+
+        For each image column, computes the intensity-weighted centroid of
+        non-zero pixels.  Uses the second derivative of the initial centroid
+        trace (Savitzky-Golay filtered) to disambiguate ambiguous columns:
+        upward concavity → take the topmost pixels; downward concavity →
+        take the bottommost pixels; otherwise use all non-zero pixels.
+
+        Parameters
+        ----------
+        img : np.ndarray
+            Binary uint8 mask for a single lead cell (H × W).
+        window_length : int, optional
+            Savitzky-Golay filter window (must be odd, default ``11``).
+        polyorder : int, optional
+            Savitzky-Golay polynomial order (default ``2``).
+
+        Returns
+        -------
+        x_coords : np.ndarray of int
+            Column indices where a signal was detected.
+        final_signal : np.ndarray of float
+            Weighted-centroid Y-coordinate for each column in *x_coords*.
+        """
+        height, width = img.shape
+        x_coords = []
+        initial_signal = []
+
+        for col in range(width):
+            column = img[:, col]
+            if np.sum(column) == 0:
+                continue
+
+            y_indices = np.arange(height)
+            weights = column.astype(float)
+            centroid = np.average(y_indices, weights=weights)
+            initial_signal.append(centroid)
+            x_coords.append(col)
+
+        x_coords = np.array(x_coords)
+        initial_signal = np.array(initial_signal)
+
+        second_deriv = np.gradient(np.gradient(savgol_filter(initial_signal, window_length, polyorder)))
+        second_deriv = savgol_filter(second_deriv, 11, 3)
+
+        x_out, final_signal = [], []
+        for idx, col in enumerate(x_coords):
+            column = img[:, col]
+            nz_idx = np.where(column > 0)[0]
+            if nz_idx.size == 0:
+                continue
+
+            if second_deriv[idx] > 0.5:
+                sel_idx = nz_idx[:5]
+            elif second_deriv[idx] < -0.5:
+                sel_idx = nz_idx[-5:]
+            else:
+                sel_idx = nz_idx
+
+            weights  = column[sel_idx].astype(float)
+            centroid = np.average(sel_idx, weights=weights)
+            x_out.append(col)
+            final_signal.append(centroid)
+
+        return np.array(x_out), np.array(final_signal)
+
+    def fill_gaps(self, x_coords, y_coords, method='linear'):
+        """Interpolate missing columns to produce a dense, evenly-spaced signal.
+
+        Parameters
+        ----------
+        x_coords : np.ndarray of int
+            Sparse column indices from :meth:`binarize_signal`.
+        y_coords : np.ndarray of float
+            Corresponding centroid Y-values.
+        method : str, optional
+            Interpolation kind passed to ``scipy.interpolate.interp1d``
+            (default ``'linear'``).
+
+        Returns
+        -------
+        x_full : np.ndarray of int
+            Dense array ``[x_coords[0], …, x_coords[-1]]``.
+        y_interp : np.ndarray of float
+            Interpolated Y-values at every integer column.
+        """
+        x_full = np.arange(x_coords[0], x_coords[-1] + 1)
+        interpolator = interp1d(x_coords, y_coords, kind=method, fill_value="extrapolate")
+        y_interp = interpolator(x_full)
+        return x_full, y_interp
+
+    def smooth_signal(self, sig, window_length=7, polyorder=4):
+        """Apply a Savitzky-Golay low-pass filter to a 1-D signal.
+
+        Parameters
+        ----------
+        sig : np.ndarray
+            Input signal array.
+        window_length : int, optional
+            Filter window size in samples (must be odd, default ``7``).
+        polyorder : int, optional
+            Polynomial order (default ``4``).
+
+        Returns
+        -------
+        np.ndarray
+            Smoothed signal of the same length as *sig*.
+        """
+        return savgol_filter(sig, window_length, polyorder)
+
+    def extract_signals(self):
+        """Convert mask grid cells into calibrated voltage time-series.
+
+        For each cell in ``self.grid``:
+
+        1. Calls :meth:`binarize_signal` to get column-centroid Y-values.
+        2. Trims 5 edge samples, fills gaps, and smooths with Savitzky-Golay.
+        3. Converts pixel Y-coordinates to millivolts using
+           ``self.baseline`` and ``self.volt_per_pixel``.
+        4. Converts column indices to seconds using ``self.time_per_pixel``.
+        5. Resamples uniformly at 500 Hz via linear interpolation.
+
+        Sets ``self.signal_grid`` with the same structure as ``self.grid``
+        but each cell additionally contains:
+
+        * ``'time'``   – np.ndarray of time values in seconds.
+        * ``'signal'`` – np.ndarray of voltage values in millivolts.
+        """
+        signal_grid = []
+        sample_rate = 500
+
+        for row_idx, row in enumerate(self.grid):
+            baseline_y = self.baseline[row_idx]
+            row_signals = []
+
+            for cell in row:
+                x_coords, sig = self.binarize_signal(cell['signal'])
+
+                x_coords = x_coords[5:-5]
+                sig = sig[5:-5]
+
+                x_coords, sig = self.fill_gaps(x_coords, sig, method='linear')
+                sig = self.smooth_signal(sig)
+
+                signal_volts = (baseline_y - sig) * self.volt_per_pixel
+                x_seconds = x_coords * self.time_per_pixel
+
+                duration = x_seconds[-1] - x_seconds[0]
+                num_samples = round(duration * sample_rate) + 1
+                resampled_time = np.linspace(x_seconds[0], x_seconds[-1], num_samples)
+
+                interpolator = interp1d(x_seconds, signal_volts, kind='linear',
+                                        fill_value="extrapolate")
+                resampled_signal = interpolator(resampled_time)
+
+                row_signals.append({
+                    'lead': cell['lead'],
+                    'time': resampled_time,
+                    'signal': resampled_signal
+                })
+
+            signal_grid.append(row_signals)
+
+        self.signal_grid = signal_grid
+
+    # ------------------------------------------------------------------
+    # Metrics & comparison
+    # ------------------------------------------------------------------
+
+    def sliding_metrics(self, signal_a, signal_b, return_aligned_signals=False):
+        """Compute alignment-aware similarity metrics between two signals.
+
+        Slides the shorter signal across the longer one and finds the offset
+        that maximises the Pearson correlation.  Reports RMSE and SNR at
+        the best-aligned position.
+
+        Parameters
+        ----------
+        signal_a, signal_b : array-like
+            The two signals to compare (need not be the same length).
+        return_aligned_signals : bool, optional
+            If ``True``, also return the aligned signal pair (default ``False``).
+
+        Returns
+        -------
+        max_corr : float
+            Maximum Pearson correlation found across all offsets.
+        rmse : float
+            Root mean squared error at the best offset.
+        snr : float
+            Signal-to-noise ratio in dB at the best offset
+            (``inf`` if noise power is zero).
+        short_signal : np.ndarray
+            The shorter signal (only when *return_aligned_signals* is ``True``).
+        best_window : np.ndarray
+            The aligned window from the longer signal
+            (only when *return_aligned_signals* is ``True``).
+        """
+        len_a, len_b = len(signal_a), len(signal_b)
+        if len_a > len_b:
+            long_signal, short_signal = signal_a, signal_b
+        else:
+            long_signal, short_signal = signal_b, signal_a
+
+        len_long, len_short          = len(long_signal), len(short_signal)
+        max_corr, best_corr_offset   = -1, 0
+
+        for i in range(len_long - len_short + 1):
+            window = long_signal[i:i + len_short]
+            corr, _ = pearsonr(window, short_signal)
+            if corr > max_corr:
+                max_corr         = corr
+                best_corr_offset = i
+
+        best_window  = long_signal[best_corr_offset:best_corr_offset + len_short]
+        rmse         = mean_squared_error(best_window, short_signal)
+        signal_power = np.mean(np.square(short_signal))
+        noise_power  = np.mean(np.square(best_window - short_signal))
+        snr          = (10 * np.log10(signal_power / noise_power)
+                        if noise_power > 0 else np.inf)
+
+        if return_aligned_signals:
+            return max_corr, rmse, snr, short_signal, best_window
+        return max_corr, rmse, snr
+
+    def calculate_metrics_ptb(self, plot_signals=True, per_lead_scores=None):
+        """Evaluate digitization quality against a WFDB ground-truth record.
+
+        For each lead in ``self.signal_grid``, aligns the extracted signal
+        to the corresponding WFDB signal using :meth:`sliding_metrics` and
+        records Pearson r, RMSE, SNR, and p-value.  Requires ``self.wfdb_path``
+        to point to a valid WFDB record.
+
+        Sets the following instance attributes after completion:
+
+        * ``self.average_pearson`` – mean Pearson r across all leads.
+        * ``self.average_rmse``    – mean RMSE (mV).
+        * ``self.average_snr``     – mean SNR (dB).
+        * ``self.average_pval``    – mean p-value.
+
+        Parameters
+        ----------
+        plot_signals : bool, optional
+            Unused placeholder (default ``True``).
+        per_lead_scores : dict or None, optional
+            If provided, per-lead metrics are appended here for leads where
+            Pearson r > 0.60.  Expected structure:
+            ``{lead_name: {'pearson': [], 'rmse': [], 'snr': [], 'pval': []}}``.
+        """
+        record = wfdb.rdrecord(self.wfdb_path)
+        avg_pearson, avg_rmse, avg_snr, avg_pval = [], [], [], []
+
+        for row in self.signal_grid:
+            for cell in row:
+                if 'lead' not in cell or 'signal' not in cell:
+                    continue
+                try:
+                    lead_index = record.sig_name.index(cell['lead'])
+                except ValueError:
+                    continue
+
+                wfdb_signal    = record.p_signal[:, lead_index]
+                wfdb_signal    = wfdb_signal[~np.isnan(wfdb_signal)]
+                voltage_signal = np.array(cell['signal'])
+
+                pearson_val, rmse, snr, sig1, sig2 = self.sliding_metrics(
+                    voltage_signal, wfdb_signal, return_aligned_signals=True
+                )
+
+                try:
+                    pearson, pval = pearsonr(sig1, sig2)
+                except Exception:
+                    pearson, pval = np.nan, np.nan
+
+                cell.update({'pearson': pearson, 'rmse': rmse,
+                             'snr': snr,         'pval': pval})
+
+                avg_pearson.append(pearson)
+                avg_rmse.append(rmse)
+                avg_snr.append(snr)
+                avg_pval.append(pval)
+
+                if per_lead_scores is not None and pearson > 0.60:
+                    lead = cell['lead']
+                    if lead not in per_lead_scores:
+                        per_lead_scores[lead] = {
+                            'pearson': [], 'rmse': [], 'snr': [], 'pval': []
+                        }
+                    per_lead_scores[lead]['pearson'].append(pearson)
+                    per_lead_scores[lead]['rmse'].append(rmse)
+                    per_lead_scores[lead]['snr'].append(snr)
+                    per_lead_scores[lead]['pval'].append(pval)
+
+        self.average_pearson = np.mean(avg_pearson)
+        self.average_rmse    = np.mean(avg_rmse)
+        self.average_snr     = np.mean(avg_snr)
+        self.average_pval    = np.mean(avg_pval)
+
+    def plot_signals(self, title='', plot_wfdb=False):
+        """Plot extracted signals, optionally overlaid with the WFDB ground truth.
+
+        Generates one figure per lead.  When *plot_wfdb* is ``True``, the
+        extracted signal and the best-aligned WFDB window are plotted
+        together for visual comparison.
+
+        Parameters
+        ----------
+        title : str, optional
+            Figure title applied to every subplot (default ``''``).
+        plot_wfdb : bool, optional
+            If ``True``, overlay the ground-truth WFDB signal.  Requires
+            ``self.wfdb_path`` and a loaded ``record`` in scope
+            (default ``False``).
+        """
+        for row in self.signal_grid:
+            for cell in row:
+                voltage_signal = cell['signal']
+                plt.figure(figsize=(10, 4))
+                if not plot_wfdb:
+                    plt.plot(voltage_signal, linewidth=1.5)
+                else:
+                    lead_index  = record.sig_name.index(cell['lead'])
+                    wfdb_signal = record.p_signal[:, lead_index]
+                    wfdb_signal = [x for x in wfdb_signal if not np.isnan(x)]
+                    _, _, _, sig1, sig2 = self.sliding_metrics(
+                        voltage_signal, wfdb_signal, return_aligned_signals=True
+                    )
+                    plt.plot(sig1, label='Extracted Signal', linewidth=1.5)
+                    plt.plot(sig2, label='Ground Truth',     linewidth=1.5)
+                plt.title(title)
+                plt.legend()
+                plt.xlabel("Time (ms)")
+                plt.ylabel("Voltage (mV)")
+                plt.tight_layout()
+                plt.show()
+
+    # ------------------------------------------------------------------
+    # Export
+    # ------------------------------------------------------------------
+
+    def save_signals_as_wfdb(self, record_name, directory='.'):
+        """Export extracted signals as a WFDB record (.hea + .dat files).
+
+        Pads all lead signals to the same length with NaN (converted to 0),
+        then writes a 500 Hz, 16-bit WFDB record in millivolts.
+
+        Parameters
+        ----------
+        record_name : str
+            Base name for the output files (no extension).
+        directory : str, optional
+            Destination directory (default ``'.'``).
+        """
+        signals, lead_names, max_length = [], [], 0
+        for row in self.signal_grid:
+            for cell in row:
+                if 'signal' in cell and 'lead' in cell:
+                    sig = np.array(cell['signal'])
+                    signals.append(sig)
+                    lead_names.append(cell['lead'])
+                    max_length = max(max_length, len(sig))
+
+        padded = [np.pad(s, (0, max_length - len(s)),
+                         mode='constant', constant_values=np.nan)
+                  for s in signals]
+        signal_array = np.nan_to_num(np.array(padded).T, nan=0.0)
+
+        wfdb.wrsamp(
+            record_name=record_name,
+            fs=500,
+            units=['mV'] * len(lead_names),
+            sig_name=lead_names,
+            p_signal=signal_array,
+            fmt=['16'] * len(lead_names),
+            write_dir=directory,
+        )
+
+    def save_signals_as_csv(self, record_name, directory='.'):
+        """
+        Save ECG signals to CSV. Each column = one lead, each row = one time sample.
+        """
+        signals, lead_names, max_length = [], [], 0
+        for row in self.signal_grid:
+            for cell in row:
+                if 'signal' in cell and 'lead' in cell:
+                    sig = np.asarray(cell['signal'], dtype=float)
+                    signals.append(sig)
+                    lead_names.append(cell['lead'])
+                    max_length = max(max_length, len(sig))
+
+        if not signals:
+            raise ValueError("No valid signals found in signal_grid.")
+
+        padded = [np.pad(s, (0, max_length - len(s)),
+                         mode='constant', constant_values=np.nan)
+                  for s in signals]
+        df = pd.DataFrame(np.vstack(padded).T, columns=lead_names)
+
+        os.makedirs(directory, exist_ok=True)
+        csv_path = os.path.join(directory, f"{record_name}.csv")
+        df.to_csv(csv_path, index=False, float_format="%.6f")
+
+    # ------------------------------------------------------------------
+    # Full pipeline
+    # ------------------------------------------------------------------
+
+    def run_full_pipeline(self):
+        """Run all digitization stages end-to-end for a single ECG image.
+
+        Executes the following stages in order, with per-stage wall-clock
+        timing printed to stdout via :meth:`_timed`:
+
+        1. **Image loading & preprocessing** — ``load_image`` + ``preprocess_image``.
+        2. **Segmentation** — ``segment_leads`` → ``make_segmentation_mask``
+           → ``find_row_centers`` → ``get_roi``.
+        3. **Sequential YOLO detections** — ``extract_lead_boxes``,
+           ``extract_lead_name_boxes``, ``extract_reference_pulses``.
+           (Sequential, not parallel, to avoid a PyTorch thread-pool deadlock.)
+        4. **Retry loop** — if any of the three detections is empty, the
+           image is reloaded at the next fallback target size from
+           ``[2000, 2100, 1900, 2200, 1800, 1700]`` px and stages 1–3 are
+           re-run.  Raises ``RuntimeError`` if all sizes fail.
+        5. **Calibration** — ``get_reference_scale``.
+        6. **Grid & signal extraction** — ``make_grid`` → ``extract_signals``.
+
+        Raises
+        ------
+        RuntimeError
+            If all fallback target sizes yield at least one empty detection.
+        """
+        FALLBACK_SIZES = [2000, 2100, 1900, 2200, 1800, 1700]
+
+        for target_size in FALLBACK_SIZES:
+
+            # -- Image loading & preprocessing --
+            self._timed("load_image",       self.load_image, self.image_path,
+                        target_size=target_size)
+            self._timed("preprocess",       self.preprocess_image)
+
+            # -- Segmentation (sequential: each step depends on the previous) --
+            self._timed("segment_leads",    self.segment_leads)
+            self._timed("make_mask",        self.make_segmentation_mask)
+            self._timed("find_row_centers", self.find_row_centers)
+            self._timed("get_roi",          self.get_roi)
+
+            # -- Three YOLO calls sequentially --
+            self._timed("extract_lead_boxes",       self.extract_lead_boxes)
+            self._timed("extract_lead_name_boxes",  self.extract_lead_name_boxes)
+            self._timed("extract_reference_pulses", self.extract_reference_pulses)
+
+            # -- Check all three detections; retry at next size if any is empty --
+            missing = []
+            if len(self.reference_pulses) == 0:
+                missing.append("pulses")
+            if len(self.lead_bboxes) == 0:
+                missing.append("lead boxes")
+            if len(self.lead_name_bboxes) == 0:
+                missing.append("lead names")
+
+            if not missing:
+                if target_size != FALLBACK_SIZES[0]:
+                    print(f"    ↳ All detections found at target_size={target_size}",
+                          flush=True)
+                break
+            else:
+                print(f"    ↳ Missing {', '.join(missing)} at target_size={target_size},"
+                      f" retrying...", flush=True)
+
+        else:
+            raise RuntimeError(
+                f"Detections incomplete at all target sizes {FALLBACK_SIZES}. "
+                f"Last missing: {', '.join(missing)}."
+            )
+
+        # -- Remainder of pipeline (sequential) --
+        self._timed("get_reference_scale", self.get_reference_scale)
+        self._timed("make_grid",           self.make_grid)
+        self._timed("extract_signals",     self.extract_signals)

environment.yml

@@ -0,0 +1,194 @@
+name: infer
+channels:
+  - main
+  - r
+  - conda-forge
+  - bioconda
+  - pytorch
+dependencies:
+  - _libgcc_mutex=0.1=main
+  - _openmp_mutex=5.1=1_gnu
+  - aiohappyeyeballs=2.4.4=py39h06a4308_0
+  - aiohttp=3.11.10=py39h5eee18b_0
+  - aiosignal=1.2.0=pyhd3eb1b0_0
+  - aom=3.6.0=h6a678d5_0
+  - async-timeout=5.0.1=py39h06a4308_0
+  - attrs=24.3.0=py39h06a4308_0
+  - blas=1.0=openblas
+  - blosc=1.21.3=h6a678d5_0
+  - bottleneck=1.4.2=py39ha9d4c09_0
+  - brotli=1.0.9=h5eee18b_9
+  - brotli-bin=1.0.9=h5eee18b_9
+  - brotli-python=1.0.9=py39h6a678d5_9
+  - brunsli=0.1=h2531618_0
+  - bzip2=1.0.8=h5eee18b_6
+  - c-ares=1.19.1=h5eee18b_0
+  - ca-certificates=2025.2.25=h06a4308_0
+  - certifi=2025.4.26=py39h06a4308_0
+  - cffi=1.17.1=py39h1fdaa30_1
+  - cfitsio=3.470=h5893167_7
+  - charls=2.2.0=h2531618_0
+  - dav1d=1.2.1=h5eee18b_0
+  - freetype=2.13.3=h4a9f257_0
+  - frozenlist=1.5.0=py39h5eee18b_0
+  - fsspec=2025.3.2=py39he3bba80_0
+  - gettext=0.21.0=hedfda30_2
+  - giflib=5.2.2=h5eee18b_0
+  - icu=73.1=h6a678d5_0
+  - imagecodecs=2024.9.22=py39h652d938_1
+  - imageio=2.37.0=py39h06a4308_0
+  - importlib_resources=6.4.0=py39h06a4308_0
+  - joblib=1.4.2=py39h06a4308_0
+  - jpeg=9e=h5eee18b_3
+  - jxrlib=1.1=h7b6447c_2
+  - krb5=1.20.1=h143b758_1
+  - lazy_loader=0.4=py39h06a4308_0
+  - lcms2=2.16=h92b89f2_1
+  - ld_impl_linux-64=2.40=h12ee557_0
+  - lerc=4.0.0=h6a678d5_0
+  - libaec=1.1.3=h6a678d5_0
+  - libavif=1.1.1=h5eee18b_0
+  - libbrotlicommon=1.0.9=h5eee18b_9
+  - libbrotlidec=1.0.9=h5eee18b_9
+  - libbrotlienc=1.0.9=h5eee18b_9
+  - libcurl=8.12.1=hc9e6f67_0
+  - libdeflate=1.22=h5eee18b_0
+  - libedit=3.1.20230828=h5eee18b_0
+  - libev=4.33=h7f8727e_1
+  - libffi=3.4.4=h6a678d5_1
+  - libflac=1.3.4=h27087fc_0
+  - libgcc-ng=11.2.0=h1234567_1
+  - libgfortran-ng=11.2.0=h00389a5_1
+  - libgfortran5=11.2.0=h1234567_1
+  - libgomp=11.2.0=h1234567_1
+  - libnghttp2=1.57.0=h2d74bed_0
+  - libogg=1.3.5=h27cfd23_1
+  - libopenblas=0.3.29=ha39b09d_0
+  - libopus=1.3.1=h5eee18b_1
+  - libpng=1.6.39=h5eee18b_0
+  - libsndfile=1.0.31=h9c3ff4c_1
+  - libssh2=1.11.1=h251f7ec_0
+  - libstdcxx-ng=11.2.0=h1234567_1
+  - libtiff=4.7.0=hde9077f_0
+  - libvorbis=1.3.7=h7b6447c_0
+  - libwebp-base=1.3.2=h5eee18b_1
+  - libxml2=2.13.8=hfdd30dd_0
+  - libzopfli=1.0.3=he6710b0_0
+  - lz4-c=1.9.4=h6a678d5_1
+  - matplotlib-base=3.9.2=py39hbfdbfaf_1
+  - multidict=6.1.0=py39h5eee18b_0
+  - ncurses=6.4=h6a678d5_0
+  - networkx=3.2.1=py39h06a4308_0
+  - numexpr=2.10.1=py39hd28fd6d_0
+  - numpy=1.26.4=py39heeff2f4_0
+  - numpy-base=1.26.4=py39h8a23956_0
+  - openjpeg=2.5.2=h0d4d230_1
+  - openssl=3.0.16=h5eee18b_0
+  - pandas=2.2.3=py39h6a678d5_0
+  - pip=25.1=pyhc872135_2
+  - propcache=0.3.1=py39h5eee18b_0
+  - pybind11-abi=4=hd3eb1b0_1
+  - pycparser=2.21=pyhd3eb1b0_0
+  - pysocks=1.7.1=py39h06a4308_0
+  - pysoundfile=0.11.0=pyhd8ed1ab_0
+  - python=3.9.21=he870216_1
+  - python-dateutil=2.9.0post0=py39h06a4308_2
+  - python-tzdata=2025.2=pyhd3eb1b0_0
+  - readline=8.2=h5eee18b_0
+  - requests=2.32.3=py39h06a4308_1
+  - scikit-image=0.24.0=py39h1128e8f_0
+  - scikit-learn=1.6.1=py39h6a678d5_0
+  - scipy=1.13.1=py39heeff2f4_1
+  - setuptools=78.1.1=py39h06a4308_0
+  - six=1.17.0=py39h06a4308_0
+  - snappy=1.2.1=h6a678d5_0
+  - sqlite=3.45.3=h5eee18b_0
+  - threadpoolctl=3.5.0=py39h2f386ee_0
+  - tifffile=2023.4.12=py39h06a4308_0
+  - tk=8.6.14=h39e8969_0
+  - typing_extensions=4.12.2=py39h06a4308_0
+  - tzdata=2025b=h04d1e81_0
+  - unicodedata2=15.1.0=py39h5eee18b_1
+  - wfdb=4.3.0=pyhd8ed1ab_0
+  - wheel=0.45.1=py39h06a4308_0
+  - xz=5.6.4=h5eee18b_1
+  - yarl=1.18.0=py39h5eee18b_0
+  - zfp=1.0.0=h6a678d5_0
+  - zipp=3.21.0=py39h06a4308_0
+  - zlib=1.2.13=h5eee18b_1
+  - zlib-ng=2.0.7=h5eee18b_0
+  - zstd=1.5.6=hc292b87_0
+  - pip:
+      - asttokens==3.0.0
+      - charset-normalizer==3.4.2
+      - comm==0.2.2
+      - contourpy==1.3.0
+      - cycler==0.12.1
+      - debugpy==1.8.14
+      - decorator==5.2.1
+      - exceptiongroup==1.3.0
+      - executing==2.2.0
+      - filelock==3.18.0
+      - fonttools==4.58.0
+      - idna==3.10
+      - importlib-metadata==8.7.0
+      - importlib-resources==6.5.2
+      - ipykernel==6.29.5
+      - ipython==8.18.1
+      - jedi==0.19.2
+      - jinja2==3.1.6
+      - jupyter-client==8.6.3
+      - jupyter-core==5.7.2
+      - kiwisolver==1.4.7
+      - markupsafe==3.0.2
+      - matplotlib==3.9.4
+      - matplotlib-inline==0.1.7
+      - mpmath==1.3.0
+      - nest-asyncio==1.6.0
+      - nvidia-cublas-cu12==12.6.4.1
+      - nvidia-cuda-cupti-cu12==12.6.80
+      - nvidia-cuda-nvrtc-cu12==12.6.77
+      - nvidia-cuda-runtime-cu12==12.6.77
+      - nvidia-cudnn-cu12==9.5.1.17
+      - nvidia-cufft-cu12==11.3.0.4
+      - nvidia-cufile-cu12==1.11.1.6
+      - nvidia-curand-cu12==10.3.7.77
+      - nvidia-cusolver-cu12==11.7.1.2
+      - nvidia-cusparse-cu12==12.5.4.2
+      - nvidia-cusparselt-cu12==0.6.3
+      - nvidia-nccl-cu12==2.26.2
+      - nvidia-nvjitlink-cu12==12.6.85
+      - nvidia-nvtx-cu12==12.6.77
+      - opencv-python==4.11.0.86
+      - packaging==25.0
+      - parso==0.8.4
+      - patched-yolo-infer==1.3.8
+      - pexpect==4.9.0
+      - pillow==11.2.1
+      - platformdirs==4.3.8
+      - prompt-toolkit==3.0.51
+      - psutil==7.0.0
+      - ptyprocess==0.7.0
+      - pure-eval==0.2.3
+      - py-cpuinfo==9.0.0
+      - pygments==2.19.1
+      - pyparsing==3.2.3
+      - pytz==2025.2
+      - pyyaml==6.0.2
+      - pyzmq==26.4.0
+      - seaborn==0.13.2
+      - simpleitk==2.5.0
+      - stack-data==0.6.3
+      - sympy==1.14.0
+      - torch==2.7.0
+      - torchvision==0.22.0
+      - tornado==6.4.2
+      - tqdm==4.67.1
+      - traitlets==5.14.3
+      - triton==3.3.0
+      - typing-extensions==4.13.2
+      - ultralytics==8.3.133
+      - ultralytics-thop==2.0.14
+      - urllib3==2.4.0
+      - wcwidth==0.2.13
+prefix: /cluster/apps/scaimed/cachtari/conda/envs/infer

models/classifier_models/ecg_surgery/inceptiontime_model.pkl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f6d52babc1f5c750f89c9a442ccd5fdb7835e6aac2cba31c39c184d59305c282
+size 2153

models/classifier_models/ecg_surgery/model_metadata.json

@@ -0,0 +1,62 @@
+{
+  "model_name": "InceptionTime",
+  "model_file": "inceptiontime_model.pkl",
+  "use_case": "Pre vs Post Procedural MI (Surgery)",
+  "positive_class": "pre-procedural MI",
+  "negative_class": "post-procedural MI",
+  "class_labels": [
+    "post-procedural MI",
+    "pre-procedural MI"
+  ],
+  "features": [
+    "I",
+    "aVR",
+    "V1",
+    "V4",
+    "II",
+    "aVL",
+    "V2",
+    "V5",
+    "III",
+    "aVF",
+    "V3",
+    "V6"
+  ],
+  "n_features": 12,
+  "n_timesteps": 140,
+  "hyperparameters": {
+    "n_epochs": 200,
+    "batch_size": 64,
+    "random_state": 42,
+    "verbose": true
+  },
+  "training_data": "data/ecg_surgery_segmented_50_150_70.csv",
+  "training_date": "2026-06-02 20:26:15",
+  "training_time_seconds": 251.83,
+  "preprocessing": {
+    "normalization": "per_subject_max_abs",
+    "smote_tomek": true,
+    "test_size": 0.25,
+    "random_state": 42
+  },
+  "training_samples": 1104,
+  "test_samples": 338,
+  "test_metrics": {
+    "accuracy": 0.9142011834319527,
+    "precision": 0.8974358974358975,
+    "recall": 0.9150326797385621,
+    "f1": 0.9061488673139159,
+    "sensitivity": 0.9150326797385621,
+    "specificity": 0.9135135135135135,
+    "confusion_matrix": [
+      [
+        169,
+        16
+      ],
+      [
+        13,
+        140
+      ]
+    ]
+  }
+}

models/classifier_models/mi_vs_normal_segmented/arsenal_segmented_model.pkl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fb22793560c0939a574e1cdab6e59b3e732cd8592be10cad7dd3fc46e46134c8
+size 112024703

models/classifier_models/mi_vs_normal_segmented/model_metadata.json

@@ -0,0 +1,61 @@
+{
+  "model_name": "Arsenal",
+  "model_file": "arsenal_segmented_model.pkl",
+  "use_case": "MI vs Normal (Segmented)",
+  "positive_class": "MYOCARDIAL_INFARCTION",
+  "negative_class": "NORMAL",
+  "class_labels": [
+    "NORMAL",
+    "MYOCARDIAL_INFARCTION"
+  ],
+  "features": [
+    "I",
+    "aVR",
+    "V1",
+    "V4",
+    "II",
+    "aVL",
+    "V2",
+    "V5",
+    "III",
+    "aVF",
+    "V3",
+    "V6"
+  ],
+  "n_features": 12,
+  "n_timesteps": 140,
+  "hyperparameters": {
+    "num_kernels": 20000,
+    "random_state": 42,
+    "n_jobs": 1
+  },
+  "training_data": "data/ptb_xl/segmented_heartbeats.csv",
+  "training_date": "2026-06-03 03:58:21",
+  "training_time_seconds": 10302.47,
+  "preprocessing": {
+    "normalization": "per_subject_max_abs",
+    "smote_tomek": true,
+    "test_size": 0.25,
+    "random_state": 42
+  },
+  "training_samples": 6954,
+  "test_samples": 2349,
+  "test_metrics": {
+    "accuracy": 0.922520221370796,
+    "precision": 0.940864960282436,
+    "recall": 0.9026248941574937,
+    "f1": 0.9213483146067416,
+    "sensitivity": 0.9026248941574937,
+    "specificity": 0.9426369863013698,
+    "confusion_matrix": [
+      [
+        1101,
+        67
+      ],
+      [
+        115,
+        1066
+      ]
+    ]
+  }
+}

models/classifier_models/omi_vs_nonomi/model_metadata.json

@@ -0,0 +1,61 @@
+{
+  "model_name": "Rocket",
+  "model_file": "rocket_model.pkl",
+  "use_case": "OMI vs non-OMI",
+  "positive_class": "OMI",
+  "negative_class": "non-OMI",
+  "class_labels": [
+    "non-OMI",
+    "OMI"
+  ],
+  "features": [
+    "I",
+    "aVR",
+    "V1",
+    "V4",
+    "II",
+    "aVL",
+    "V2",
+    "V5",
+    "III",
+    "aVF",
+    "V3",
+    "V6"
+  ],
+  "n_features": 12,
+  "n_timesteps": 141,
+  "hyperparameters": {
+    "num_kernels": 50000,
+    "random_state": 42,
+    "n_jobs": -1
+  },
+  "training_data": "data/ecg_matrix_omi_segmented_50_150_90.csv",
+  "training_date": "2026-06-02 16:41:15",
+  "training_time_seconds": 132.78,
+  "preprocessing": {
+    "normalization": "per_subject_max_abs",
+    "smote_tomek": true,
+    "test_size": 0.25,
+    "random_state": 42
+  },
+  "training_samples": 2084,
+  "test_samples": 666,
+  "test_metrics": {
+    "accuracy": 0.8888888888888888,
+    "precision": 0.8571428571428571,
+    "recall": 0.9171974522292994,
+    "f1": 0.8861538461538462,
+    "sensitivity": 0.9171974522292994,
+    "specificity": 0.8636363636363636,
+    "confusion_matrix": [
+      [
+        304,
+        48
+      ],
+      [
+        26,
+        288
+      ]
+    ]
+  }
+}

models/classifier_models/omi_vs_nonomi/rocket_model.pkl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d35281f43a57aae562b0f929c7912e5bb8c17efad4dc3a50907fe93cf4f40318
+size 11204282

models/digitization_models/yolo11_full/weights/best.pt

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2a9adc74478e051e36652555cc81e35d787435a4799fa7a8abde6c6fa117a68a
+size 124776986

models/digitization_models/yolo11_lead/.ipynb_checkpoints/args-checkpoint.yaml

@@ -0,0 +1,105 @@
+task: detect
+mode: train
+model: yolo11x.pt
+data: /cluster/home/cachtari/dataset_leadname.yaml
+epochs: 100
+time: null
+patience: 100
+batch: 6
+imgsz: 1500
+save: true
+save_period: -1
+cache: false
+device: 0,1
+workers: 8
+project: /cluster/home/cachtari/runs
+name: yolo11_1500_leadname_class
+exist_ok: false
+pretrained: true
+optimizer: auto
+verbose: true
+seed: 0
+deterministic: true
+single_cls: false
+rect: false
+cos_lr: false
+close_mosaic: 10
+resume: false
+amp: true
+fraction: 1.0
+profile: false
+freeze: null
+multi_scale: false
+overlap_mask: true
+mask_ratio: 4
+dropout: 0.0
+val: true
+split: val
+save_json: false
+conf: null
+iou: 0.7
+max_det: 300
+half: false
+dnn: false
+plots: true
+source: null
+vid_stride: 1
+stream_buffer: false
+visualize: false
+augment: false
+agnostic_nms: false
+classes: null
+retina_masks: false
+embed: null
+show: false
+save_frames: false
+save_txt: false
+save_conf: false
+save_crop: false
+show_labels: true
+show_conf: true
+show_boxes: true
+line_width: null
+format: torchscript
+keras: false
+optimize: false
+int8: false
+dynamic: false
+simplify: true
+opset: null
+workspace: null
+nms: false
+lr0: 0.01
+lrf: 0.01
+momentum: 0.937
+weight_decay: 0.0005
+warmup_epochs: 3.0
+warmup_momentum: 0.8
+warmup_bias_lr: 0.1
+box: 7.5
+cls: 0.5
+dfl: 1.5
+pose: 12.0
+kobj: 1.0
+nbs: 64
+hsv_h: 0.015
+hsv_s: 0.7
+hsv_v: 0.4
+degrees: 0.0
+translate: 0.1
+scale: 0.5
+shear: 0.0
+perspective: 0.0
+flipud: 0.0
+fliplr: 0.5
+bgr: 0.0
+mosaic: 1.0
+mixup: 0.0
+cutmix: 0.0
+copy_paste: 0.0
+copy_paste_mode: flip
+auto_augment: randaugment
+erasing: 0.4
+cfg: null
+tracker: botsort.yaml
+save_dir: /cluster/home/cachtari/runs/yolo11_1500_leadname_class

models/digitization_models/yolo11_lead/.ipynb_checkpoints/results-checkpoint.csv

@@ -0,0 +1,101 @@
+epoch,time,train/box_loss,train/cls_loss,train/dfl_loss,metrics/precision(B),metrics/recall(B),metrics/mAP50(B),metrics/mAP50-95(B),val/box_loss,val/cls_loss,val/dfl_loss,lr/pg0,lr/pg1,lr/pg2
+1,886.991,0.93774,0.93982,0.89298,0.98053,0.95723,0.98816,0.78539,0.80694,0.38397,0.85301,0.00333207,0.00333207,0.00333207
+2,1741.49,0.79997,0.45797,0.84868,0.9877,0.96985,0.99114,0.81433,0.7171,0.35068,0.84086,0.00659942,0.00659942,0.00659942
+3,2587.75,0.79563,0.44381,0.84741,0.9897,0.97185,0.99174,0.82059,0.70896,0.34325,0.83936,0.00980076,0.00980076,0.00980076
+4,3432.96,0.75251,0.41456,0.84102,0.98979,0.97597,0.99211,0.84977,0.62835,0.30703,0.82137,0.009703,0.009703,0.009703
+5,4277.5,0.69831,0.37499,0.8299,0.99061,0.97892,0.99275,0.85163,0.62284,0.29359,0.81764,0.009604,0.009604,0.009604
+6,5120.98,0.6707,0.35552,0.82548,0.99217,0.98008,0.99302,0.86206,0.60218,0.28213,0.813,0.009505,0.009505,0.009505
+7,5963.46,0.65473,0.34407,0.82271,0.9925,0.98065,0.99334,0.86668,0.59208,0.27651,0.81129,0.009406,0.009406,0.009406
+8,6809.02,0.64014,0.33376,0.82023,0.99242,0.98181,0.99349,0.87175,0.57154,0.26426,0.80764,0.009307,0.009307,0.009307
+9,7654.95,0.6281,0.32268,0.81734,0.99407,0.98343,0.99377,0.87428,0.57133,0.25716,0.80706,0.009208,0.009208,0.009208
+10,8498.47,0.61846,0.31632,0.81624,0.99353,0.98488,0.99378,0.87652,0.56228,0.25264,0.80488,0.009109,0.009109,0.009109
+11,9341.04,0.6143,0.3126,0.81357,0.99451,0.98566,0.99391,0.87779,0.55769,0.24658,0.80416,0.00901,0.00901,0.00901
+12,10184.5,0.60103,0.30526,0.81238,0.99462,0.98598,0.99398,0.88037,0.55125,0.24366,0.8027,0.008911,0.008911,0.008911
+13,11031.8,0.59679,0.30283,0.81192,0.99405,0.98645,0.99409,0.87994,0.55084,0.24197,0.8026,0.008812,0.008812,0.008812
+14,11876.7,0.59333,0.29933,0.80987,0.99453,0.98626,0.99402,0.88191,0.54564,0.23862,0.80157,0.008713,0.008713,0.008713
+15,12722.1,0.58581,0.29494,0.80972,0.99523,0.98615,0.99412,0.88326,0.53836,0.23728,0.80025,0.008614,0.008614,0.008614
+16,13565.6,0.58006,0.29061,0.80786,0.99542,0.98692,0.99403,0.88504,0.53513,0.23455,0.79935,0.008515,0.008515,0.008515
+17,14411.3,0.57672,0.28795,0.80764,0.99561,0.98806,0.99427,0.88538,0.53908,0.23224,0.79983,0.008416,0.008416,0.008416
+18,15258.7,0.57058,0.28694,0.80589,0.99513,0.98829,0.99422,0.88413,0.53386,0.23022,0.79868,0.008317,0.008317,0.008317
+19,16103.3,0.56536,0.28123,0.80563,0.99493,0.98831,0.99425,0.8901,0.52304,0.22738,0.79724,0.008218,0.008218,0.008218
+20,16946.7,0.56263,0.28067,0.80593,0.99553,0.98797,0.99426,0.88536,0.53237,0.22661,0.79834,0.008119,0.008119,0.008119
+21,17789.4,0.55984,0.27929,0.80478,0.99557,0.98861,0.99427,0.88764,0.52399,0.226,0.7972,0.00802,0.00802,0.00802
+22,18638.8,0.55793,0.27731,0.80402,0.99598,0.98838,0.99433,0.89168,0.51603,0.22054,0.79583,0.007921,0.007921,0.007921
+23,19485,0.55206,0.27382,0.80311,0.99603,0.98875,0.99435,0.8932,0.51208,0.2185,0.79524,0.007822,0.007822,0.007822
+24,20326.6,0.55068,0.27161,0.8036,0.99586,0.98887,0.99437,0.89249,0.51388,0.21833,0.79516,0.007723,0.007723,0.007723
+25,21169.1,0.54745,0.27037,0.80235,0.99581,0.98928,0.99437,0.89536,0.50525,0.21598,0.79397,0.007624,0.007624,0.007624
+26,22012.1,0.54547,0.26931,0.80247,0.99554,0.99001,0.99433,0.89836,0.49678,0.21205,0.79302,0.007525,0.007525,0.007525
+27,22860.5,0.54133,0.26642,0.80207,0.99579,0.99015,0.99439,0.89703,0.4992,0.21184,0.79319,0.007426,0.007426,0.007426
+28,23703.3,0.53816,0.26553,0.80118,0.99645,0.98994,0.99437,0.89788,0.49825,0.20991,0.79271,0.007327,0.007327,0.007327
+29,24544.8,0.53313,0.2634,0.80064,0.99592,0.99047,0.99434,0.8979,0.49459,0.20961,0.79266,0.007228,0.007228,0.007228
+30,25388.8,0.52973,0.26017,0.79985,0.99619,0.99031,0.99446,0.90008,0.48961,0.20792,0.79182,0.007129,0.007129,0.007129
+31,26233.2,0.5269,0.25997,0.79837,0.99702,0.98979,0.9944,0.89964,0.49105,0.20738,0.79174,0.00703,0.00703,0.00703
+32,27079.8,0.52688,0.25916,0.79927,0.99672,0.99078,0.9944,0.90069,0.48737,0.20495,0.79115,0.006931,0.006931,0.006931
+33,27921.5,0.52129,0.25693,0.79811,0.99683,0.99049,0.99454,0.90165,0.48401,0.20341,0.79081,0.006832,0.006832,0.006832
+34,28764.2,0.52114,0.25471,0.79824,0.99719,0.99018,0.99447,0.90317,0.47991,0.20193,0.79012,0.006733,0.006733,0.006733
+35,29606.6,0.51987,0.25446,0.79883,0.9966,0.99087,0.99438,0.90382,0.47828,0.20192,0.79008,0.006634,0.006634,0.006634
+36,30456.1,0.51496,0.25208,0.79743,0.99719,0.99037,0.99437,0.90508,0.47383,0.20009,0.78931,0.006535,0.006535,0.006535
+37,31300.5,0.511,0.24972,0.79762,0.99658,0.99091,0.99437,0.90613,0.47076,0.19873,0.7887,0.006436,0.006436,0.006436
+38,32143.9,0.51136,0.25042,0.79623,0.99692,0.99084,0.99445,0.90558,0.47097,0.19822,0.78883,0.006337,0.006337,0.006337
+39,32985.8,0.50541,0.24745,0.79565,0.99724,0.99077,0.99448,0.90597,0.46946,0.19736,0.78834,0.006238,0.006238,0.006238
+40,33830.3,0.50479,0.24695,0.79518,0.99724,0.99083,0.99452,0.90742,0.46542,0.19641,0.78804,0.006139,0.006139,0.006139
+41,34677.4,0.5029,0.24614,0.79591,0.99696,0.99136,0.99454,0.90771,0.46354,0.19568,0.78762,0.00604,0.00604,0.00604
+42,35520.5,0.49842,0.24446,0.79503,0.99725,0.99107,0.99452,0.90858,0.46129,0.19448,0.78728,0.005941,0.005941,0.005941
+43,36364.6,0.49839,0.24419,0.7955,0.99732,0.99082,0.99444,0.90947,0.45884,0.19314,0.78682,0.005842,0.005842,0.005842
+44,37206.5,0.49216,0.24015,0.79435,0.99722,0.99107,0.99449,0.91033,0.45611,0.19261,0.78638,0.005743,0.005743,0.005743
+45,38049.1,0.49418,0.24203,0.79385,0.99703,0.99144,0.99451,0.91087,0.4543,0.19214,0.78613,0.005644,0.005644,0.005644
+46,38897.8,0.48824,0.23884,0.79406,0.99717,0.99129,0.99453,0.91115,0.45407,0.19199,0.78622,0.005545,0.005545,0.005545
+47,39741.3,0.48797,0.23914,0.79344,0.99707,0.99147,0.99454,0.91148,0.45203,0.19103,0.78606,0.005446,0.005446,0.005446
+48,40584.9,0.48281,0.23633,0.79246,0.99712,0.9916,0.99454,0.91221,0.44899,0.19046,0.7859,0.005347,0.005347,0.005347
+49,41428.6,0.4808,0.23643,0.79376,0.99757,0.99099,0.99447,0.91283,0.4469,0.19003,0.78564,0.005248,0.005248,0.005248
+50,42275,0.4799,0.23503,0.7923,0.99698,0.99156,0.99457,0.91367,0.44459,0.18939,0.78536,0.005149,0.005149,0.005149
+51,43117.5,0.47796,0.23351,0.79202,0.99712,0.99146,0.99457,0.91408,0.44334,0.18854,0.7851,0.00505,0.00505,0.00505
+52,43961,0.47425,0.23329,0.79199,0.99745,0.99117,0.9946,0.91474,0.44112,0.18791,0.78484,0.004951,0.004951,0.004951
+53,44804,0.46837,0.22882,0.79024,0.99725,0.99152,0.99463,0.91509,0.43921,0.1873,0.78449,0.004852,0.004852,0.004852
+54,45648,0.46714,0.22827,0.79055,0.99736,0.99148,0.99464,0.91561,0.4377,0.18655,0.78424,0.004753,0.004753,0.004753
+55,46495.3,0.46614,0.22808,0.79074,0.99736,0.99148,0.99463,0.91603,0.43667,0.18608,0.78411,0.004654,0.004654,0.004654
+56,47339.5,0.4637,0.2255,0.78996,0.99686,0.99221,0.99465,0.91641,0.43518,0.18521,0.78393,0.004555,0.004555,0.004555
+57,48184.1,0.46273,0.22496,0.78987,0.99748,0.99156,0.99458,0.91676,0.43384,0.18469,0.78368,0.004456,0.004456,0.004456
+58,49029.3,0.45682,0.2228,0.79042,0.99749,0.99158,0.99462,0.91735,0.43217,0.18397,0.78336,0.004357,0.004357,0.004357
+59,49877.6,0.45594,0.22277,0.79048,0.99748,0.99154,0.99461,0.9177,0.43028,0.18332,0.78314,0.004258,0.004258,0.004258
+60,50721.8,0.45355,0.2194,0.78919,0.99744,0.99169,0.99461,0.91808,0.42898,0.18286,0.78299,0.004159,0.004159,0.004159
+61,51566.9,0.44861,0.21824,0.7893,0.99742,0.99164,0.99466,0.91869,0.42741,0.18213,0.78286,0.00406,0.00406,0.00406
+62,52412.3,0.44858,0.21778,0.78838,0.99736,0.99171,0.99467,0.91925,0.42557,0.18148,0.78267,0.003961,0.003961,0.003961
+63,53257.7,0.44843,0.21688,0.78885,0.99726,0.99179,0.99467,0.91961,0.42442,0.18094,0.78249,0.003862,0.003862,0.003862
+64,54107,0.44477,0.21446,0.78753,0.99715,0.99188,0.99467,0.91993,0.42315,0.18063,0.78229,0.003763,0.003763,0.003763
+65,54953,0.44175,0.21357,0.7874,0.99734,0.99172,0.99467,0.92038,0.42182,0.18018,0.78208,0.003664,0.003664,0.003664
+66,55800,0.43928,0.21206,0.78695,0.99751,0.99154,0.99469,0.92063,0.42075,0.17969,0.78192,0.003565,0.003565,0.003565
+67,56645.5,0.43963,0.21207,0.78765,0.9974,0.99168,0.99469,0.92097,0.41954,0.17931,0.78178,0.003466,0.003466,0.003466
+68,57496.2,0.43514,0.21014,0.78565,0.99751,0.99162,0.99468,0.92127,0.41848,0.17882,0.78167,0.003367,0.003367,0.003367
+69,58350.2,0.42996,0.20736,0.78585,0.99755,0.9915,0.99462,0.92163,0.41739,0.17858,0.78156,0.003268,0.003268,0.003268
+70,59199.3,0.42742,0.20551,0.78586,0.99759,0.9915,0.99461,0.9218,0.41656,0.1783,0.78142,0.003169,0.003169,0.003169
+71,60048.6,0.42606,0.20497,0.78549,0.99755,0.99153,0.99469,0.92203,0.41563,0.17809,0.78128,0.00307,0.00307,0.00307
+72,60897.2,0.42453,0.20353,0.78485,0.99738,0.99163,0.99471,0.92237,0.41471,0.17778,0.78115,0.002971,0.002971,0.002971
+73,61751.8,0.42075,0.20147,0.78419,0.99736,0.99168,0.99471,0.92261,0.41389,0.17746,0.781,0.002872,0.002872,0.002872
+74,62622.8,0.41807,0.20066,0.78485,0.99734,0.99175,0.99471,0.92299,0.41298,0.17724,0.7809,0.002773,0.002773,0.002773
+75,63476.1,0.41766,0.20019,0.78364,0.99736,0.99163,0.99472,0.92325,0.41214,0.17706,0.78079,0.002674,0.002674,0.002674
+76,64327.2,0.416,0.19909,0.78397,0.99737,0.99158,0.99473,0.92332,0.41149,0.17701,0.78071,0.002575,0.002575,0.002575
+77,65176.8,0.41124,0.19658,0.78418,0.99716,0.99183,0.99473,0.92365,0.41066,0.17665,0.78061,0.002476,0.002476,0.002476
+78,66035,0.40958,0.19567,0.78388,0.99712,0.99182,0.99473,0.92393,0.40989,0.17629,0.78052,0.002377,0.002377,0.002377
+79,66886.2,0.40919,0.19443,0.78283,0.99715,0.9918,0.99473,0.92408,0.40921,0.1762,0.78042,0.002278,0.002278,0.002278
+80,67738.1,0.40715,0.19381,0.78292,0.99708,0.99192,0.99472,0.92434,0.40846,0.17597,0.78031,0.002179,0.002179,0.002179
+81,68589.2,0.40328,0.19213,0.78175,0.99709,0.99194,0.99472,0.9245,0.40775,0.17579,0.78023,0.00208,0.00208,0.00208
+82,69443.4,0.40319,0.19144,0.78211,0.99719,0.99182,0.99472,0.92463,0.40726,0.17563,0.78012,0.001981,0.001981,0.001981
+83,70289.6,0.39893,0.18903,0.78224,0.99705,0.99196,0.99471,0.92486,0.40658,0.17557,0.78004,0.001882,0.001882,0.001882
+84,71137.2,0.39812,0.18907,0.78186,0.99678,0.99223,0.9947,0.92504,0.40588,0.17558,0.77995,0.001783,0.001783,0.001783
+85,71985.2,0.39451,0.18653,0.78083,0.99717,0.9918,0.9947,0.92519,0.40524,0.17548,0.77985,0.001684,0.001684,0.001684
+86,72835.6,0.39281,0.18653,0.78087,0.99708,0.99194,0.9947,0.92533,0.40454,0.17543,0.77977,0.001585,0.001585,0.001585
+87,73688.6,0.38824,0.18342,0.78001,0.99712,0.99192,0.99469,0.92568,0.40373,0.17536,0.77965,0.001486,0.001486,0.001486
+88,74539.5,0.38613,0.18234,0.78023,0.99715,0.99184,0.9947,0.92593,0.40293,0.17522,0.77955,0.001387,0.001387,0.001387
+89,75389.8,0.38569,0.18154,0.78027,0.99716,0.99181,0.9947,0.92611,0.40243,0.17525,0.77946,0.001288,0.001288,0.001288
+90,76239.2,0.38303,0.17992,0.7801,0.99713,0.99185,0.9947,0.92634,0.40173,0.17533,0.77937,0.001189,0.001189,0.001189
+91,77086.7,0.37095,0.17454,0.77936,0.99716,0.99186,0.9947,0.92656,0.40097,0.17513,0.77928,0.00109,0.00109,0.00109
+92,77938.8,0.36715,0.17267,0.77984,0.9971,0.99194,0.99469,0.9268,0.40027,0.17496,0.77918,0.000991,0.000991,0.000991
+93,78809.8,0.36022,0.16862,0.77824,0.99715,0.99183,0.99469,0.92707,0.39938,0.17471,0.77907,0.000892,0.000892,0.000892
+94,79725.6,0.36091,0.16907,0.77796,0.99718,0.99181,0.99469,0.92743,0.39847,0.17434,0.77896,0.000793,0.000793,0.000793
+95,80641,0.35823,0.16708,0.77799,0.99714,0.99188,0.99469,0.92771,0.39764,0.17414,0.77886,0.000694,0.000694,0.000694
+96,81570.1,0.35413,0.16608,0.77726,0.99715,0.9918,0.99469,0.92805,0.39675,0.17394,0.77875,0.000595,0.000595,0.000595
+97,82477.9,0.35209,0.16355,0.77688,0.99703,0.9919,0.99469,0.92822,0.39607,0.17364,0.77865,0.000496,0.000496,0.000496
+98,83378.5,0.34975,0.1626,0.77636,0.99695,0.99201,0.99469,0.92863,0.39533,0.17326,0.77857,0.000397,0.000397,0.000397
+99,84261.4,0.34865,0.16152,0.77566,0.99699,0.99195,0.99469,0.92895,0.39459,0.1728,0.77846,0.000298,0.000298,0.000298
+100,85165.5,0.34551,0.15998,0.77677,0.99728,0.99157,0.99469,0.92913,0.39388,0.17234,0.77837,0.000199,0.000199,0.000199