Delete archive

archive/classification/aggregate_subject_metrics.py

@@ -1,96 +0,0 @@
-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

@@ -1,66 +0,0 @@
-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

@@ -1,93 +0,0 @@
-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

@@ -1,71 +0,0 @@
-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

@@ -1,494 +0,0 @@
-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

@@ -1,64 +0,0 @@
-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

@@ -1,34 +0,0 @@
-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

@@ -1,297 +0,0 @@
-"""
-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

@@ -1,128 +0,0 @@
-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

@@ -1,427 +0,0 @@
-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

@@ -1,21 +0,0 @@
-# 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

@@ -1,302 +0,0 @@
-# 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)