#!/usr/bin/env python3 """ EEG Deep Analysis: Spectral Connectivity + Graph Metrics + Temporal CNN Features ================================================================================ Advanced 3-way (AD vs FTD vs Control) and binary (AD vs non-AD) classification using graph-theoretic features from coherence networks, Hjorth parameters, envelope statistics, and ensemble learning. Dataset: OpenNeuro ds004504 (88 subjects: 36 AD, 23 FTD, 29 Control) Author: Satyawan Singh — Infonova Solutions """ import os import json import time import pickle import warnings import numpy as np import pandas as pd warnings.filterwarnings('ignore') # NumPy 2.0 compat: np.trapz -> np.trapezoid if not hasattr(np, 'trapz'): np.trapz = np.trapezoid import mne from scipy import signal from scipy.stats import kurtosis, skew from scipy.ndimage import uniform_filter1d from sklearn.model_selection import StratifiedKFold, cross_val_predict from sklearn.preprocessing import StandardScaler, LabelEncoder from sklearn.ensemble import ( GradientBoostingClassifier, RandomForestClassifier, VotingClassifier, ) from sklearn.svm import SVC from sklearn.metrics import ( classification_report, confusion_matrix, roc_auc_score, accuracy_score, ) from sklearn.feature_selection import SelectKBest, f_classif from sklearn.pipeline import Pipeline # ── Paths ── BASE = '/Users/satyawansingh/Documents/alzheimer-research-complete/data/openneuro_ad_eeg' OUTPUT_DIR = '/Users/satyawansingh/Documents/alzheimer-research-complete/models/eeg_deep_analysis' os.makedirs(OUTPUT_DIR, exist_ok=True) CHANNELS = [ 'Fp1', 'Fp2', 'F3', 'F4', 'C3', 'C4', 'P3', 'P4', 'O1', 'O2', 'F7', 'F8', 'T3', 'T4', 'T5', 'T6', 'Fz', 'Cz', 'Pz', ] N_CH = len(CHANNELS) BANDS = { 'theta': (4, 8), 'alpha': (8, 13), 'beta': (13, 30), 'gamma': (30, 45), } # ══════════════════════════════════════════════════════════════ # FEATURE EXTRACTION FUNCTIONS # ══════════════════════════════════════════════════════════════ def compute_coherence_matrix(data, sfreq, band): """Compute pairwise coherence matrix for a given frequency band.""" fmin, fmax = band n_ch = min(data.shape[0], N_CH) coh_matrix = np.zeros((n_ch, n_ch)) for i in range(n_ch): for j in range(i, n_ch): if i == j: coh_matrix[i, j] = 1.0 continue freqs, coh = signal.coherence( data[i], data[j], fs=sfreq, nperseg=min(1024, len(data[i])), ) band_mask = (freqs >= fmin) & (freqs <= fmax) if band_mask.sum() > 0: val = np.mean(coh[band_mask]) else: val = 0.0 coh_matrix[i, j] = val coh_matrix[j, i] = val return coh_matrix def graph_metrics_from_matrix(adj, band_name): """ Extract graph-theoretic metrics from a coherence adjacency matrix. Uses only numpy/scipy — no networkx dependency. """ features = {} n = adj.shape[0] # Threshold to create binary graph (top 30% connections) upper = adj[np.triu_indices(n, k=1)] if len(upper) == 0 or np.std(upper) < 1e-12: # Degenerate matrix — return zeros for key in [ 'mean_coh', 'std_coh', 'global_efficiency', 'clustering_coeff', 'char_path_length', 'small_worldness', 'modularity_approx', 'degree_entropy', 'hub_score_max', 'hub_score_std', 'assortativity_approx', ]: features[f'{band_name}_{key}'] = 0.0 return features threshold = np.percentile(upper, 70) binary = (adj >= threshold).astype(float) np.fill_diagonal(binary, 0) # --- Weighted metrics --- features[f'{band_name}_mean_coh'] = np.mean(upper) features[f'{band_name}_std_coh'] = np.std(upper) # --- Degree & hub scores (weighted) --- strength = adj.sum(axis=1) - 1 # subtract self-connection degree = binary.sum(axis=1) max_degree = degree.max() if degree.max() > 0 else 1 hub_scores = degree / max_degree features[f'{band_name}_hub_score_max'] = hub_scores.max() features[f'{band_name}_hub_score_std'] = hub_scores.std() # Degree entropy degree_norm = degree / (degree.sum() + 1e-10) degree_norm = degree_norm[degree_norm > 0] features[f'{band_name}_degree_entropy'] = -np.sum(degree_norm * np.log2(degree_norm + 1e-15)) # --- Clustering coefficient (binary) --- # C_i = 2 * triangles_i / (k_i * (k_i - 1)) triangles = np.diag(binary @ binary @ binary) / 2 k = degree denom = k * (k - 1) denom[denom == 0] = 1 cc = 2 * triangles / denom features[f'{band_name}_clustering_coeff'] = np.mean(cc) # --- Characteristic path length via Floyd-Warshall on binary --- # Distance matrix: 1/weight for connected, inf for disconnected dist = np.full((n, n), np.inf) np.fill_diagonal(dist, 0) connected = binary > 0 # Use inverse coherence as distance for weighted path dist[connected] = 1.0 / (adj[connected] + 1e-10) # Floyd-Warshall for k_node in range(n): new_dist = dist[:, k_node, None] + dist[None, k_node, :] dist = np.minimum(dist, new_dist) finite_dists = dist[np.triu_indices(n, k=1)] finite_dists = finite_dists[np.isfinite(finite_dists)] if len(finite_dists) > 0: cpl = np.mean(finite_dists) else: cpl = np.inf features[f'{band_name}_char_path_length'] = cpl if np.isfinite(cpl) else 100.0 # --- Global efficiency: mean of 1/d_ij --- inv_dist = 1.0 / (dist + 1e-10) np.fill_diagonal(inv_dist, 0) features[f'{band_name}_global_efficiency'] = inv_dist.sum() / (n * (n - 1)) # --- Small-worldness approximation --- # sigma = (C / C_rand) / (L / L_rand) # For Erdos-Renyi with same density: C_rand ~ p, L_rand ~ ln(n) / ln(k_mean) p = binary.sum() / (n * (n - 1)) c_rand = max(p, 1e-10) k_mean = degree.mean() l_rand = np.log(n) / (np.log(k_mean + 1) + 1e-10) if k_mean > 1 else 10.0 C_real = features[f'{band_name}_clustering_coeff'] L_real = features[f'{band_name}_char_path_length'] sigma = (C_real / (c_rand + 1e-10)) / (L_real / (l_rand + 1e-10) + 1e-10) features[f'{band_name}_small_worldness'] = sigma # --- Modularity approximation (spectral bisection) --- # Q = 1/(2m) * sum( A_ij - k_i*k_j/(2m) ) * delta(c_i, c_j) m = binary.sum() / 2 if m > 0: B = binary - np.outer(degree, degree) / (2 * m + 1e-10) eigvals, eigvecs = np.linalg.eigh(B) # Partition based on sign of leading eigenvector partition = (eigvecs[:, -1] > 0).astype(int) same_comm = np.outer(partition, partition) + np.outer(1 - partition, 1 - partition) Q = np.sum(B * same_comm) / (4 * m + 1e-10) features[f'{band_name}_modularity_approx'] = Q else: features[f'{band_name}_modularity_approx'] = 0.0 # --- Assortativity approximation --- # Correlation of degrees at each end of edges edges_i, edges_j = np.where(np.triu(binary, k=1) > 0) if len(edges_i) > 2: d_i = degree[edges_i] d_j = degree[edges_j] if np.std(d_i) > 0 and np.std(d_j) > 0: features[f'{band_name}_assortativity_approx'] = np.corrcoef(d_i, d_j)[0, 1] else: features[f'{band_name}_assortativity_approx'] = 0.0 else: features[f'{band_name}_assortativity_approx'] = 0.0 return features def compute_temporal_cnn_features(data, sfreq): """ Extract temporal / signal-morphology features per channel: - Hjorth mobility & complexity - Signal envelope statistics - Zero-crossing rate - Line length - Higuchi fractal dimension approximation """ features = {} n_ch = min(data.shape[0], N_CH) for ch_idx in range(n_ch): ch = CHANNELS[ch_idx] x = data[ch_idx] # --- Basic stats --- features[f'{ch}_mean'] = np.mean(x) features[f'{ch}_std'] = np.std(x) features[f'{ch}_kurtosis'] = kurtosis(x) features[f'{ch}_skewness'] = skew(x) features[f'{ch}_rms'] = np.sqrt(np.mean(x ** 2)) # --- Hjorth parameters --- diff1 = np.diff(x) diff2 = np.diff(diff1) activity = np.var(x) mobility = np.sqrt(np.var(diff1) / (activity + 1e-10)) complexity = np.sqrt(np.var(diff2) / (np.var(diff1) + 1e-10)) / (mobility + 1e-10) features[f'{ch}_hjorth_activity'] = activity features[f'{ch}_hjorth_mobility'] = mobility features[f'{ch}_hjorth_complexity'] = complexity # --- Zero-crossing rate --- zcr = np.sum(np.diff(np.sign(x)) != 0) / len(x) features[f'{ch}_zcr'] = zcr # --- Signal envelope (analytic signal) --- analytic = signal.hilbert(x) envelope = np.abs(analytic) features[f'{ch}_env_mean'] = np.mean(envelope) features[f'{ch}_env_std'] = np.std(envelope) features[f'{ch}_env_skew'] = skew(envelope) features[f'{ch}_env_kurtosis'] = kurtosis(envelope) # --- Line length (sum of absolute differences) --- features[f'{ch}_line_length'] = np.mean(np.abs(diff1)) # --- Higuchi fractal dimension (fast approximation, k_max=10) --- k_max = 10 N_pts = len(x) lk = [] for k in range(1, k_max + 1): lengths = [] for m in range(1, k + 1): idx = np.arange(m - 1, N_pts, k) if len(idx) < 2: continue seg = x[idx] L_m = np.sum(np.abs(np.diff(seg))) * (N_pts - 1) / (k * len(seg) * k) lengths.append(L_m) if lengths: lk.append(np.mean(lengths)) if len(lk) > 2: ks = np.arange(1, len(lk) + 1) log_k = np.log(ks) log_lk = np.log(np.array(lk) + 1e-15) # Linear fit slope = fractal dimension slope, _ = np.polyfit(log_k, log_lk, 1) features[f'{ch}_hfd'] = -slope else: features[f'{ch}_hfd'] = 0.0 return features def compute_band_power_features(data, sfreq): """Compute per-channel PSD band powers and key spectral ratios.""" features = {} n_ch = min(data.shape[0], N_CH) for ch_idx in range(n_ch): ch = CHANNELS[ch_idx] x = data[ch_idx] freqs, psd = signal.welch(x, fs=sfreq, nperseg=min(2048, len(x))) total = np.trapz(psd, freqs) + 1e-10 for bname, (fmin, fmax) in BANDS.items(): mask = (freqs >= fmin) & (freqs <= fmax) bp = np.trapz(psd[mask], freqs[mask]) features[f'{ch}_{bname}_rel'] = bp / total # Delta band too delta_mask = (freqs >= 0.5) & (freqs <= 4) delta_power = np.trapz(psd[delta_mask], freqs[delta_mask]) features[f'{ch}_delta_rel'] = delta_power / total # Spectral ratios alpha_mask = (freqs >= 8) & (freqs <= 13) theta_mask = (freqs >= 4) & (freqs <= 8) alpha_p = np.trapz(psd[alpha_mask], freqs[alpha_mask]) theta_p = np.trapz(psd[theta_mask], freqs[theta_mask]) features[f'{ch}_theta_alpha_ratio'] = theta_p / (alpha_p + 1e-10) features[f'{ch}_delta_alpha_ratio'] = delta_power / (alpha_p + 1e-10) # Peak alpha frequency alpha_freqs = freqs[alpha_mask] alpha_psd = psd[alpha_mask] if len(alpha_psd) > 0: features[f'{ch}_peak_alpha_freq'] = alpha_freqs[np.argmax(alpha_psd)] else: features[f'{ch}_peak_alpha_freq'] = 0 # Spectral entropy psd_norm = psd / (psd.sum() + 1e-10) psd_pos = psd_norm[psd_norm > 0] features[f'{ch}_spectral_entropy'] = -np.sum(psd_pos * np.log2(psd_pos)) return features # ══════════════════════════════════════════════════════════════ # STEP 1: Load participants # ══════════════════════════════════════════════════════════════ print("=" * 70) print(" EEG DEEP ANALYSIS: Connectivity Graphs + Temporal + Ensemble") print("=" * 70) participants = pd.read_csv(os.path.join(BASE, 'participants.tsv'), sep='\t') print(f"\nParticipants: {len(participants)}") print(f"Groups: {dict(participants['Group'].value_counts())}") label_map = {'A': 0, 'C': 1, 'F': 2} label_names = {0: 'AD', 1: 'Control', 2: 'FTD'} # ══════════════════════════════════════════════════════════════ # STEP 2: Feature extraction # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" STEP 2: Extracting features (coherence graphs + temporal + PSD)") print("=" * 70) all_features = [] all_labels = [] all_subjects = [] failed = [] t0 = time.time() for idx, row in participants.iterrows(): sub_id = row['participant_id'] group = row['Group'] label = label_map[group] eeg_file = os.path.join(BASE, sub_id, 'eeg', f'{sub_id}_task-eyesclosed_eeg.set') if not os.path.exists(eeg_file): failed.append(sub_id) continue try: raw = mne.io.read_raw_eeglab(eeg_file, preload=True, verbose=False) raw.filter(0.5, 45, verbose=False) data = raw.get_data() sfreq = raw.info['sfreq'] feats = {} # 1) Spectral connectivity graph metrics per band for band_name, band_range in BANDS.items(): coh_mat = compute_coherence_matrix(data, sfreq, band_range) graph_feats = graph_metrics_from_matrix(coh_mat, band_name) feats.update(graph_feats) # Also store upper-triangle coherence stats per band upper = coh_mat[np.triu_indices(N_CH, k=1)] feats[f'{band_name}_coh_median'] = np.median(upper) feats[f'{band_name}_coh_q25'] = np.percentile(upper, 25) feats[f'{band_name}_coh_q75'] = np.percentile(upper, 75) # 2) Temporal / signal-morphology features feats.update(compute_temporal_cnn_features(data, sfreq)) # 3) Band power features feats.update(compute_band_power_features(data, sfreq)) # 4) Demographics feats['age'] = row['Age'] feats['gender'] = 1 if row['Gender'] == 'M' else 0 feats['mmse'] = row['MMSE'] all_features.append(feats) all_labels.append(label) all_subjects.append(sub_id) elapsed = time.time() - t0 print(f" [{idx+1:2d}/{len(participants)}] {sub_id} [{label_names[label]:>7s}] " f"— {len(feats)} features ({elapsed:.0f}s)") except Exception as e: print(f" [{idx+1:2d}/{len(participants)}] {sub_id} FAILED: {e}") failed.append(sub_id) print(f"\nExtracted: {len(all_features)} subjects | Failed: {len(failed)}") print(f"Total time: {time.time() - t0:.0f}s") # Convert to matrix X = pd.DataFrame(all_features).fillna(0) y = np.array(all_labels) # Replace inf with large finite value X = X.replace([np.inf, -np.inf], 0) print(f"Feature matrix: {X.shape}") print(f"Labels: AD={sum(y==0)}, Control={sum(y==1)}, FTD={sum(y==2)}") # Save raw features X.to_csv(os.path.join(OUTPUT_DIR, 'deep_features.csv'), index=False) np.save(os.path.join(OUTPUT_DIR, 'labels.npy'), y) # ══════════════════════════════════════════════════════════════ # STEP 3: Feature selection & scaling # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" STEP 3: Feature selection") print("=" * 70) k_features = min(120, X.shape[1]) selector = SelectKBest(f_classif, k=k_features) X_selected = selector.fit_transform(X, y) selected_mask = selector.get_support() selected_features = X.columns[selected_mask].tolist() print(f"Selected {len(selected_features)} / {X.shape[1]} features") # Top 25 features by F-score scores = selector.scores_[selected_mask] top_idx = np.argsort(scores)[::-1][:25] print("\nTop 25 most discriminative features:") for i, ix in enumerate(top_idx): print(f" {i+1:2d}. {selected_features[ix]:45s} F={scores[ix]:.1f}") scaler = StandardScaler() X_scaled = scaler.fit_transform(X_selected) # ══════════════════════════════════════════════════════════════ # STEP 4: 3-way classification (AD vs FTD vs Control) # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" STEP 4: 3-way classification (AD vs FTD vs Control)") print("=" * 70) cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) gb = GradientBoostingClassifier( n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8, random_state=42, ) rf = RandomForestClassifier( n_estimators=400, max_depth=6, min_samples_leaf=2, random_state=42, ) svm = SVC( kernel='rbf', C=10, gamma='scale', probability=True, random_state=42, ) models_3c = { 'GradientBoosting': gb, 'RandomForest': rf, 'SVM_RBF': svm, } # Also build voting ensemble voting = VotingClassifier( estimators=[('gb', gb), ('rf', rf), ('svm', svm)], voting='soft', ) models_3c['VotingEnsemble'] = voting results_3c = {} for name, model in models_3c.items(): y_pred = cross_val_predict(model, X_scaled, y, cv=cv) y_prob = cross_val_predict(model, X_scaled, y, cv=cv, method='predict_proba') acc = accuracy_score(y, y_pred) # One-vs-rest AUC try: auc = roc_auc_score(y, y_prob, multi_class='ovr', average='weighted') except Exception: auc = 0.0 results_3c[name] = {'accuracy': acc, 'auc': auc, 'y_pred': y_pred, 'y_prob': y_prob} print(f"\n{'─' * 50}") print(f" {name}: Accuracy = {acc:.1%} | AUC(OvR) = {auc:.3f}") print(f"{'─' * 50}") print(classification_report(y, y_pred, target_names=['AD', 'Control', 'FTD'])) print("Confusion matrix:") cm = confusion_matrix(y, y_pred) print(f" {'':>10s} pred_AD pred_Ctrl pred_FTD") for i, lbl in enumerate(['AD', 'Control', 'FTD']): print(f" {lbl:>10s} {cm[i,0]:7d} {cm[i,1]:9d} {cm[i,2]:8d}") best_3c = max(results_3c, key=lambda k: results_3c[k]['accuracy']) print(f"\n>>> Best 3-class model: {best_3c} " f"(Acc={results_3c[best_3c]['accuracy']:.1%}, " f"AUC={results_3c[best_3c]['auc']:.3f})") # ══════════════════════════════════════════════════════════════ # STEP 5: Binary classification (AD vs non-AD) # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" STEP 5: Binary classification (AD vs non-AD)") print("=" * 70) y_binary = (y == 0).astype(int) # AD=1, non-AD=0 gb_b = GradientBoostingClassifier( n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8, random_state=42, ) rf_b = RandomForestClassifier( n_estimators=400, max_depth=6, min_samples_leaf=2, random_state=42, ) svm_b = SVC( kernel='rbf', C=10, gamma='scale', probability=True, random_state=42, ) voting_b = VotingClassifier( estimators=[('gb', gb_b), ('rf', rf_b), ('svm', svm_b)], voting='soft', ) models_bin = { 'GradientBoosting': gb_b, 'RandomForest': rf_b, 'SVM_RBF': svm_b, 'VotingEnsemble': voting_b, } cv_bin = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) results_bin = {} for name, model in models_bin.items(): y_pred = cross_val_predict(model, X_scaled, y_binary, cv=cv_bin) y_prob = cross_val_predict(model, X_scaled, y_binary, cv=cv_bin, method='predict_proba') acc = accuracy_score(y_binary, y_pred) auc = roc_auc_score(y_binary, y_prob[:, 1]) sens = np.sum((y_pred == 1) & (y_binary == 1)) / (np.sum(y_binary == 1) + 1e-10) spec = np.sum((y_pred == 0) & (y_binary == 0)) / (np.sum(y_binary == 0) + 1e-10) results_bin[name] = {'accuracy': acc, 'auc': auc, 'sensitivity': sens, 'specificity': spec} print(f"\n{'─' * 50}") print(f" {name}: Acc={acc:.1%} AUC={auc:.3f} Sens={sens:.1%} Spec={spec:.1%}") print(f"{'─' * 50}") print(classification_report(y_binary, y_pred, target_names=['non-AD', 'AD'])) best_bin = max(results_bin, key=lambda k: results_bin[k]['auc']) print(f"\n>>> Best binary model: {best_bin} " f"(AUC={results_bin[best_bin]['auc']:.3f}, " f"Acc={results_bin[best_bin]['accuracy']:.1%})") # ══════════════════════════════════════════════════════════════ # STEP 6: Train final models on full data & save # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" STEP 6: Training final models & saving artifacts") print("=" * 70) # Final 3-class final_3c = VotingClassifier( estimators=[ ('gb', GradientBoostingClassifier( n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8, random_state=42)), ('rf', RandomForestClassifier( n_estimators=400, max_depth=6, min_samples_leaf=2, random_state=42)), ('svm', SVC(kernel='rbf', C=10, gamma='scale', probability=True, random_state=42)), ], voting='soft', ) final_3c.fit(X_scaled, y) # Final binary final_bin = VotingClassifier( estimators=[ ('gb', GradientBoostingClassifier( n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8, random_state=42)), ('rf', RandomForestClassifier( n_estimators=400, max_depth=6, min_samples_leaf=2, random_state=42)), ('svm', SVC(kernel='rbf', C=10, gamma='scale', probability=True, random_state=42)), ], voting='soft', ) final_bin.fit(X_scaled, y_binary) # Feature importance from the GradientBoosting inside the ensemble gb_inside = final_3c.named_estimators_['gb'] importances = gb_inside.feature_importances_ top_fi = np.argsort(importances)[::-1][:20] print("\nTop 20 feature importances (GradientBoosting, 3-class):") for i, ix in enumerate(top_fi): print(f" {i+1:2d}. {selected_features[ix]:45s} imp={importances[ix]:.4f}") # Save feature importance fi_df = pd.DataFrame({ 'feature': selected_features, 'importance': importances, }).sort_values('importance', ascending=False) fi_df.to_csv(os.path.join(OUTPUT_DIR, 'feature_importance.csv'), index=False) # Save models artifacts = { 'model_3class': final_3c, 'model_binary': final_bin, 'scaler': scaler, 'selector': selector, 'feature_names': list(X.columns), 'selected_features': selected_features, 'label_names_3class': {0: 'AD', 1: 'Control', 2: 'FTD'}, 'label_names_binary': {0: 'non-AD', 1: 'AD'}, 'channels': CHANNELS, 'bands': BANDS, 'results_3class': {k: {kk: vv for kk, vv in v.items() if kk not in ('y_pred', 'y_prob')} for k, v in results_3c.items()}, 'results_binary': results_bin, } model_path = os.path.join(OUTPUT_DIR, 'eeg_deep_analysis.pkl') with open(model_path, 'wb') as f: pickle.dump(artifacts, f) print(f"\nModel saved: {model_path} ({os.path.getsize(model_path)/1e6:.1f} MB)") # Save results summary as JSON summary = { 'dataset': 'OpenNeuro ds004504', 'n_subjects': len(all_features), 'n_failed': len(failed), 'n_features_total': X.shape[1], 'n_features_selected': len(selected_features), 'classification_3way': { name: { 'accuracy': float(f"{v['accuracy']:.4f}"), 'auc_ovr': float(f"{v['auc']:.4f}"), } for name, v in results_3c.items() }, 'classification_binary_AD_vs_nonAD': { name: {k: float(f"{vv:.4f}") for k, vv in v.items()} for name, v in results_bin.items() }, 'best_3class_model': best_3c, 'best_binary_model': best_bin, } with open(os.path.join(OUTPUT_DIR, 'results_summary.json'), 'w') as f: json.dump(summary, f, indent=2) # ══════════════════════════════════════════════════════════════ # FINAL SUMMARY # ══════════════════════════════════════════════════════════════ print(f"\n{'=' * 70}") print(" EEG DEEP ANALYSIS — COMPLETE") print("=" * 70) print(f" Subjects: {len(all_features)} ({sum(y==0)} AD, {sum(y==1)} Ctrl, {sum(y==2)} FTD)") print(f" Features: {X.shape[1]} total -> {len(selected_features)} selected") print(f" 3-class best: {best_3c} Acc={results_3c[best_3c]['accuracy']:.1%} AUC={results_3c[best_3c]['auc']:.3f}") print(f" Binary best: {best_bin} AUC={results_bin[best_bin]['auc']:.3f} Acc={results_bin[best_bin]['accuracy']:.1%}") print(f" Output dir: {OUTPUT_DIR}") print(f" Author: Satyawan Singh — Infonova Solutions") print("=" * 70)