{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":159854820,"sourceType":"kernelVersion"}],"dockerImageVersionId":31011,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom xgboost import XGBClassifier\nfrom sklearn.model_selection import train_test_split, learning_curve\nfrom sklearn.metrics import log_loss, confusion_matrix, ConfusionMatrixDisplay, accuracy_score, precision_recall_curve, auc\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.decomposition import PCA\nfrom scipy.stats import skew, kurtosis\nfrom tqdm import tqdm\nimport multiprocessing\nimport pickle\nimport warnings\nimport random\nimport matplotlib.pyplot as plt\n\nwarnings.filterwarnings('ignore')\n\n# Configuration\nBASE_PATH = \"/kaggle/input/hms-harmful-brain-activity-classification\"\nPREPROCESSED_PATH = \"/kaggle/input/how-to-make-spectrogram-from-eeg/EEG_Spectrograms/\"  # Path to preprocessed spectrogram .npy files\nTRAIN_LABELS_PATH = os.path.join(BASE_PATH, \"train.csv\")\nMODEL_OUTPUT_PATH = os.path.join(\"/kaggle/working/\", \"models\")\nFEATURE_CACHE_PATH = os.path.join(\"/kaggle/working/\", \"feature_cache\")\nEXPECTED_CHANNELS = 128  # Adjusted to match Chris Deotte’s 4 montages (LL, LP, RR, RP)\nFEATURES_PER_CHANNEL = 9  # Mean, var, skew, kurtosis, 4 frequency bands, spectral centroid\n\n# Create output directories\nos.makedirs(MODEL_OUTPUT_PATH, exist_ok=True)\nos.makedirs(FEATURE_CACHE_PATH, exist_ok=True)\n\n# Define classes\nCLASSES = ['Seizure', 'LPD', 'GPD', 'LRDA', 'GRDA', 'Other']\nN_CLASSES = len(CLASSES)\nTARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']\n\n# Metrics functions\ndef kl_divergence(y_true, y_pred):\n    epsilon = 1e-10\n    y_pred = np.clip(y_pred, epsilon, 1 - epsilon)\n    y_true = np.clip(y_true, epsilon, 1.0)\n    kl_div = np.sum(y_true * np.log(y_true / y_pred), axis=1)\n    return np.mean(kl_div)\n\n# Check GPU availability\ndef check_gpu_availability():\n    try:\n        test_model = XGBClassifier(tree_method='gpu_hist', n_estimators=1)\n        test_model.fit(np.zeros((2, 2)), [0, 1])\n        return 'gpu_hist'\n    except:\n        print(\"GPU not available, falling back to CPU histogram\")\n        return 'hist'\n\n# Check if files exist\nif not os.path.exists(TRAIN_LABELS_PATH):\n    print(f\"Error: Training labels file not found at {TRAIN_LABELS_PATH}\")\n    print(\"Available files in BASE_PATH:\", os.listdir(BASE_PATH) if os.path.exists(BASE_PATH) else \"BASE_PATH does not exist\")\n    exit()\n\n# Load and preprocess training labels\ntry:\n    df = pd.read_csv(TRAIN_LABELS_PATH)\n    print(f\"Loaded {len(df)} annotations with {len(df['eeg_id'].unique())} unique EEG IDs\")\nexcept Exception as e:\n    print(f\"Error loading training labels: {e}\")\n    exit()\n\n# Non-overlapping EEG ID processing\ntrain = df.groupby('eeg_id')[['spectrogram_id', 'spectrogram_label_offset_seconds']].agg(\n    {'spectrogram_id': 'first', 'spectrogram_label_offset_seconds': 'min'})\ntrain.columns = ['spec_id', 'min']\n\ntmp = df.groupby('eeg_id')[['spectrogram_label_offset_seconds']].agg('max')\ntrain['max'] = tmp\n\ntmp = df.groupby('eeg_id')[['patient_id']].agg('first')\ntrain['patient_id'] = tmp\n\ntmp = df.groupby('eeg_id')[TARGETS].agg('sum')\nfor t in TARGETS:\n    train[t] = tmp[t].values\n\ny_data = train[TARGETS].values\ny_data = y_data / y_data.sum(axis=1, keepdims=True)\ntrain[TARGETS] = y_data\n\ntmp = df.groupby('eeg_id')[['expert_consensus']].agg('first')\ntrain['target'] = tmp\n\ntrain = train.reset_index()\nprint('Train non-overlap eeg_id shape:', train.shape)\nprint(train.head())\n\n# Encode hard labels\nlabel_encoder = LabelEncoder()\nlabel_encoder.fit(CLASSES)\ntrain['label'] = label_encoder.transform(train['target'])\nprint(f\"Computed labels for {len(train)} non-overlapping EEG IDs\")\n\n# Display original class distribution\nprint(\"Original class distribution:\")\nprint(train['target'].value_counts(normalize=True))\n\n# Data augmentation function for spectrograms\ndef augment_spectrogram_data(spec_data, prob=0.5):\n    augmented = spec_data.copy()\n    for channel in range(spec_data.shape[0]):  # Loop over channels\n        if random.random() < prob:\n            # Add Gaussian noise\n            noise = np.random.normal(0, 0.1 * np.std(augmented[channel, :, :]), augmented[channel, :, :].shape)\n            augmented[channel, :, :] += noise\n        if random.random() < prob:\n            # Time shift (random shift within 10% of time axis)\n            shift = int(random.uniform(-0.1, 0.1) * augmented.shape[2])\n            augmented[channel, :, :] = np.roll(augmented[channel, :, :], shift, axis=1)\n        if random.random() < prob:\n            # Scale amplitude (random scaling between 0.9 and 1.1)\n            scale = random.uniform(0.9, 1.1)\n            augmented[channel, :, :] *= scale\n    return augmented\n\n# Feature loading and caching for spectrograms\ndef load_and_cache_features(eeg_id, augment=False):\n    cache_file = os.path.join(FEATURE_CACHE_PATH, f\"{eeg_id}_aug.npy\" if augment else f\"{eeg_id}.npy\")\n    if os.path.exists(cache_file):\n        try:\n            features = np.load(cache_file)\n            if features.shape[0] == EXPECTED_CHANNELS * FEATURES_PER_CHANNEL and not np.any(np.isnan(features)) and not np.any(np.isinf(features)):\n                return features\n            else:\n                print(f\"Warning: Cached features for {eeg_id} have incorrect shape {features.shape} or contain NaN/Inf\")\n        except Exception as e:\n            print(f\"Error loading cached features for {eeg_id}: {e}\")\n    \n    try:\n        spec_path = os.path.join(PREPROCESSED_PATH, f\"{eeg_id}.npy\")\n        if not os.path.exists(spec_path):\n            print(f\"Spectrogram file not found for {eeg_id} at {spec_path}\")\n            return None\n        spec_data = np.load(spec_path).astype(np.float32)  # Expected shape: (channels, freq, time)\n        \n        if augment:\n            spec_data = augment_spectrogram_data(spec_data, prob=0.5)\n        \n        # Check shape and print for debugging\n        #print(f\"Spectrogram {eeg_id} shape: {spec_data.shape}\")\n        if spec_data.shape[0] != EXPECTED_CHANNELS:\n            print(f\"Error: Spectrogram data for {eeg_id} has {spec_data.shape[0]} channels, expected {EXPECTED_CHANNELS}\")\n            return None\n        \n        if np.any(np.isnan(spec_data)) or np.any(np.isinf(spec_data)):\n            print(f\"Warning: NaN or Inf values in spectrogram data for {eeg_id}\")\n            return None\n        \n        features = []\n        # Assume frequency bins cover 0-100 Hz, adjust if different\n        freqs = np.linspace(0, 50, spec_data.shape[1])  # Frequency axis (axis 1)\n        for channel in range(min(spec_data.shape[0], EXPECTED_CHANNELS)):  # Limit to expected channels\n            signal = spec_data[channel, :, :].flatten()  # Flatten each channel for stats\n            # Compute statistical features\n            mean_val = np.mean(signal)\n            var_val = np.var(signal)\n            skew_val = skew(signal)\n            kurt_val = kurtosis(signal)\n            if np.any(np.isnan([mean_val, var_val, skew_val, kurt_val])):\n                print(f\"Warning: NaN in statistical features for spectrogram {eeg_id}, channel {channel}\")\n                return None\n            features.extend([mean_val, var_val, skew_val, kurt_val])\n            # Compute frequency band powers (frequency axis is axis 1)\n            freq_bands = [(0, 4), (4, 8), (8, 13), (13, 30)]  # Delta, Theta, Alpha, Beta\n            for low, high in freq_bands:\n                band_indices = (freqs >= low) & (freqs < high)\n                band_power = np.mean(spec_data[channel, band_indices, :])\n                if np.isnan(band_power) or np.isinf(band_power):\n                    print(f\"Warning: NaN/Inf in band power for spectrogram {eeg_id}, channel {channel}, band ({low},{high})\")\n                    return None\n                features.append(band_power)\n            # Compute spectral centroid (weighted mean of frequencies)\n            spec_channel = np.mean(spec_data[channel, :, :], axis=1)  # Average over time\n            spec_sum = np.sum(spec_channel)\n            if spec_sum == 0 or np.isnan(spec_sum) or np.isinf(spec_sum):\n                print(f\"Warning: Zero or invalid sum for spectrogram {eeg_id}, channel {channel}, using default centroid 0\")\n                spectral_centroid = 0.0  # Default value\n            else:\n                spectral_centroid = np.sum(freqs * spec_channel) / (spec_sum + 1e-10)\n                if np.isnan(spectral_centroid) or np.isinf(spectral_centroid):\n                    print(f\"Warning: NaN/Inf in spectral centroid for spectrogram {eeg_id}, channel {channel}, using default 0\")\n                    spectral_centroid = 0.0\n            features.append(spectral_centroid)\n        \n        features = np.array(features, dtype=np.float32)\n        expected_feature_size = EXPECTED_CHANNELS * FEATURES_PER_CHANNEL\n        if features.shape[0] != expected_feature_size:\n            print(f\"Error: Feature vector for {eeg_id} has shape {features.shape}, expected ({expected_feature_size},)\")\n            return None\n        \n        if np.any(np.isnan(features)) or np.any(np.isinf(features)):\n            print(f\"Warning: NaN or Inf in final features for spectrogram {eeg_id}\")\n            return None\n        \n        try:\n            np.save(cache_file, features)\n        except Exception as e:\n            print(f\"Error saving cached features for {eeg_id}: {e}\")\n        return features\n    except Exception as e:\n        print(f\"Error loading features for {eeg_id}: {e}\")\n        return None\n\n# Check preprocessed data\nif not os.path.exists(PREPROCESSED_PATH):\n    print(f\"Error: Preprocessed spectrogram path not found at {PREPROCESSED_PATH}\")\n    print(\"Available files in /kaggle/working/:\", os.listdir(\"/kaggle/working/\") if os.path.exists(\"/kaggle/working/\") else \"Path does not exist\")\n    exit()\n\n# Load features\nprint(\"Loading features...\")\nX = []\ny_soft = []\ny_hard = []\neeg_ids = []\nskipped_samples = 0\nfor idx, row in train.iterrows():\n    eeg_id = str(row['eeg_id'])\n    # Original features\n    features = load_and_cache_features(eeg_id, augment=False)\n    if features is not None:\n        X.append(features)\n        y_soft.append(row[TARGETS].values)\n        y_hard.append(row['label'])\n        eeg_ids.append(eeg_id)\n    else:\n        skipped_samples += 1\n    # Augmented features\n    features_aug = load_and_cache_features(eeg_id, augment=True)\n    if features_aug is not None:\n        X.append(features_aug)\n        y_soft.append(row[TARGETS].values)\n        y_hard.append(row['label'])\n        eeg_ids.append(eeg_id + \"_aug\")\n    else:\n        skipped_samples += 1\n\nif skipped_samples > 0:\n    print(f\"Warning: Skipped {skipped_samples} samples due to missing or invalid spectrogram data\")\n\nif len(X) == 0:\n    print(\"Error: No features could be loaded. Check preprocessed data in\", PREPROCESSED_PATH)\n    exit()\n\ntry:\n    X = np.array(X)\n    y_soft = np.array(y_soft, dtype=np.float64)\n    y_hard = np.array(y_hard, dtype=int)\nexcept ValueError as e:\n    print(f\"Error converting X to array: {e}\")\n    print(\"Debugging feature shapes:\")\n    for i, (features, eeg_id) in enumerate(zip(X, eeg_ids)):\n        print(f\"Spectrogram {eeg_id}: Feature shape {np.array(features).shape}\")\n    exit()\n\n# Apply PCA\npca = PCA(n_components=0.95)  # Retain 95% of variance\nX_pca = pca.fit_transform(X)\nprint(f\"PCA reduced feature dimension from {X.shape[1]} to {X_pca.shape[1]}\")\nX = X_pca  # Replace X with PCA-transformed features\n\nprint(f\"Loaded features for {len(X)} samples\")\nprint(f\"Feature vector shape: {X.shape}\")\n\n# Train-validation split\nX_train, X_val, y_train, y_val, y_soft_train, y_soft_val, train_eeg_ids, val_eeg_ids = train_test_split(\n    X, y_hard, y_soft, eeg_ids, test_size=0.2, random_state=42, stratify=y_hard\n)\nprint(f\"Training samples: {len(X_train)}, Validation samples: {len(X_val)}\")\n\n# Verify class distribution\ntrain_class_counts = pd.Series(y_train).value_counts(normalize=True)\ntrain_class_counts.index = [CLASSES[i] for i in train_class_counts.index]\nval_class_counts = pd.Series(y_val).value_counts(normalize=True)\nval_class_counts.index = [CLASSES[i] for i in val_class_counts.index]\nprint(\"\\nClass distribution in training set:\")\nprint(train_class_counts)\nprint(\"\\nClass distribution in validation set:\")\nprint(val_class_counts)\n\n# Plot class distribution\nplt.figure(figsize=(8, 5))\nplt.bar(train_class_counts.index, train_class_counts.values, alpha=0.5, label='Training')\nplt.bar(val_class_counts.index, val_class_counts.values, alpha=0.5, label='Validation')\nplt.title('Class Distribution')\nplt.xlabel('Class')\nplt.ylabel('Proportion')\nplt.legend()\nplt.tight_layout()\nplt.savefig(os.path.join(MODEL_OUTPUT_PATH, 'class_distribution.png'), dpi=100, bbox_inches='tight')\nplt.close()\n\n# Set up multiprocessing\ntry:\n    max_cores = min(2, multiprocessing.cpu_count() - 1)\nexcept:\n    max_cores = 1\nprint(f\"Using {max_cores} cores\")\n\n# Fixed hyperparameters for XGBoost\nfixed_params = {\n    'n_estimators': 50,\n    'max_depth': 4,\n    'learning_rate': 0.05,\n    'subsample': 0.8,\n    'lambda': 1.0,\n    'alpha': 0.1,\n    'random_state': 42,\n    'n_jobs': max_cores,\n    'objective': 'multi:softprob',\n    'num_class': N_CLASSES,\n    'eval_metric': 'mlogloss',\n    'tree_method': check_gpu_availability(),\n    'max_bin': 64\n}\n\n# Compute class weights for training\nclass_counts = pd.Series(y_train).value_counts()\nclass_weights = {i: len(y_train) / (N_CLASSES * count) for i, count in class_counts.items()}\nsample_weights = np.array([class_weights[label] for label in y_train])\n\n# Train XGBoost model\nprint(\"Training XGBoost model...\")\nxgb_model = XGBClassifier(**fixed_params, early_stopping_rounds=10)\ntry:\n    xgb_model.fit(\n        X_train, y_train,\n        eval_set=[(X_val, y_val)],\n        sample_weight=sample_weights,\n        verbose=False\n    )\n    print(f\"Number of boosting rounds used: {xgb_model.get_booster().num_boosted_rounds()}\")\n    \n    # Evaluate XGBoost on training set\n    y_train_pred = xgb_model.predict(X_train)\n    y_train_pred_proba = xgb_model.predict_proba(X_train)\n    train_acc = accuracy_score(y_train, y_train_pred)\n    train_kl = kl_divergence(y_soft_train, y_train_pred_proba)\n    train_ll = log_loss(y_train, y_train_pred_proba)\n    print(f\"XGBoost Training Accuracy: {train_acc:.6f}\")\n    print(f\"XGBoost Training KL Divergence: {train_kl:.6f}\")\n    print(f\"XGBoost Training Log Loss: {train_ll:.6f}\")\n    \n    # Evaluate XGBoost on validation set\n    y_val_pred = xgb_model.predict(X_val)\n    y_val_pred_proba = xgb_model.predict_proba(X_val)\n    val_acc = accuracy_score(y_val, y_val_pred)\n    val_kl = kl_divergence(y_soft_val, y_val_pred_proba)\n    val_ll = log_loss(y_val, y_val_pred_proba)\n    print(f\"XGBoost Validation Accuracy: {val_acc:.6f}\")\n    print(f\"XGBoost Validation KL Divergence: {val_kl:.6f}\")\n    print(f\"XGBoost Validation Log Loss: {val_ll:.6f}\")\n    \n    # Save XGBoost model\n    xgb_model_path = os.path.join(MODEL_OUTPUT_PATH, \"xgb_model_final.pkl\")\n    with open(xgb_model_path, \"wb\") as f:\n        pickle.dump(xgb_model, f)\n    print(\"XGBoost model saved to:\", xgb_model_path)\nexcept Exception as e:\n    print(f\"Error training XGBoost model: {e}\")\n    exit()\n\n# Feature importance visualization for XGBoost\ntry:\n    feature_names = []\n    for channel in range(EXPECTED_CHANNELS):\n        feature_names.extend([\n            f\"Channel_{channel}_Mean\", f\"Channel_{channel}_Var\",\n            f\"Channel_{channel}_Skew\", f\"Channel_{channel}_Kurtosis\",\n            f\"Channel_{channel}_Delta\", f\"Channel_{channel}_Theta\",\n            f\"Channel_{channel}_Alpha\", f\"Channel_{channel}_Beta\",\n            f\"Channel_{channel}_Centroid\"\n        ])\n    feature_importance_df = pd.DataFrame({\n        'Feature': feature_names[:len(xgb_model.feature_importances_)],\n        'Importance': xgb_model.feature_importances_\n    })\n    top_features = feature_importance_df.sort_values('Importance', ascending=False)\n    print(\"Top feature importances for XGBoost:\")\n    print(top_features)\n    # Plot feature importance\n    plt.figure(figsize=(10, 6))\n    plt.barh(top_features['Feature'][:10], top_features['Importance'][:10])\n    plt.title('Top 10 Feature Importances for XGBoost')\n    plt.xlabel('Importance')\n    plt.tight_layout()\n    plt.savefig(os.path.join(MODEL_OUTPUT_PATH, 'xgb_feature_importance.png'), dpi=100, bbox_inches='tight')\n    plt.close()\nexcept Exception as e:\n    print(f\"Error creating XGBoost feature importance: {e}\")\n\n# Confusion matrix for XGBoost validation set\ntry:\n    cm = confusion_matrix(y_val, y_val_pred, normalize='true')\n    disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=CLASSES)\n    disp.plot(cmap='Blues')\n    plt.title(f'Normalized Confusion Matrix - XGBoost Validation (KL: {val_kl:.6f})')\n    plt.tight_layout()\n    plt.savefig(os.path.join(MODEL_OUTPUT_PATH, 'xgb_confusion_matrix.png'), dpi=100, bbox_inches='tight')\n    plt.close()\nexcept Exception as e:\n    print(f\"Error creating XGBoost confusion matrix: {e}\")\n\n# Learning Curves for XGBoost\ntry:\n    train_sizes, train_scores, val_scores = learning_curve(\n        XGBClassifier(**fixed_params),\n        X, y_hard,\n        cv=5,\n        scoring='neg_log_loss',\n        train_sizes=np.linspace(0.1, 1.0, 10),\n        n_jobs=max_cores\n    )\n    plt.figure(figsize=(8, 5))\n    plt.plot(train_sizes, -train_scores.mean(axis=1), label='Train Log Loss')\n    plt.plot(train_sizes, -val_scores.mean(axis=1), label='Validation Log Loss')\n    plt.title('XGBoost Learning Curves')\n    plt.xlabel('Training Samples')\n    plt.ylabel('Log Loss')\n    plt.legend()\n    plt.tight_layout()\n    plt.savefig(os.path.join(MODEL_OUTPUT_PATH, 'xgb_learning_curves.png'), dpi=100, bbox_inches='tight')\n    plt.close()\nexcept Exception as e:\n    print(f\"Error creating XGBoost learning curves: {e}\")\n\n# Precision-Recall Curves for XGBoost\ntry:\n    plt.figure(figsize=(10, 6))\n    for i, cls in enumerate(CLASSES):\n        y_val_binary = (y_val == i).astype(int)\n        y_val_pred_proba_cls = y_val_pred_proba[:, i]\n        precision, recall, _ = precision_recall_curve(y_val_binary, y_val_pred_proba_cls)\n        auc_pr = auc(recall, precision)\n        plt.plot(recall, precision, label=f'{cls} (AUC = {auc_pr:.2f})')\n    plt.title('XGBoost Precision-Recall Curves by Class')\n    plt.xlabel('Recall')\n    plt.ylabel('Precision')\n    plt.legend()\n    plt.tight_layout()\n    plt.savefig(os.path.join(MODEL_OUTPUT_PATH, 'xgb_pr_curves.png'), dpi=100, bbox_inches='tight')\n    plt.close()\nexcept Exception as e:\n    print(f\"Error creating XGBoost PR curves: {e}\")\n\nprint(\"\\nTraining completed!\")\nprint(f\"XGBoost Training KL Divergence: {train_kl:.6f}\")\nprint(f\"XGBoost Validation KL Divergence: {val_kl:.6f}\")\nprint(f\"XGBoost Training Accuracy: {train_acc:.6f}\")\nprint(f\"XGBoost Validation Accuracy: {val_acc:.6f}\")\nprint(f\"XGBoost Training Log Loss: {train_ll:.6f}\")\nprint(f\"XGBoost Validation Log Loss: {val_ll:.6f}\")\nprint(f\"Files saved in: {MODEL_OUTPUT_PATH}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T13:06:57.449486Z","iopub.execute_input":"2025-05-30T13:06:57.450123Z"}},"outputs":[],"execution_count":null}]}