Plot ica comparison

doc/changes/dev/13215.other.rst

@@ -0,0 +1 @@
+Extend :ref:`ex-ica-comp` example on comparing ICA algorithms with clean vs noisy MEG data, by :newcontrib:`Ganasekhar Kalla`.

doc/changes/names.inc

@@ -98,6 +98,7 @@
 .. _Florin Pop: https://github.com/florin-pop
 .. _Frederik Weber: https://github.com/Frederik-D-Weber
 .. _Fu-Te Wong: https://github.com/zuxfoucault
+.. _Ganasekhar Kalla: https://github.com/Ganasekhar-gif
 .. _Gennadiy Belonosov: https://github.com/Genuster
 .. _Geoff Brookshire: https://github.com/gbrookshire
 .. _George O'Neill: https://georgeoneill.github.io

examples/preprocessing/ica_comparison.py

@@ -1,50 +1,72 @@
 """
 .. _ex-ica-comp:
 
-===========================================
-Compare the different ICA algorithms in MNE
-===========================================
+===========================================================
+Compare the performance of different ICA algorithms in MNE
+===========================================================
 
-Different ICA algorithms are fit to raw MEG data, and the corresponding maps
-are displayed.
+This example compares various ICA algorithms (FastICA, Picard, Infomax,
+Extended Infomax) on the same raw MEG data. For each algorithm:
 
+- The ICA fit time (speed) is shown
+- All components (up to 20) are visualized
+- The EOG-related component from each method is detected and compared
+- Comparison on clean vs noisy data is done
+
+Note: In typical preprocessing, only one ICA algorithm is used.
+This example is for educational purposes.
 """
 # Authors: Pierre Ablin <[email protected]>
+#          Ganasekhar Kalla <[email protected]>
 #
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
 # %%
 
+from pathlib import Path
 from time import time
 
+import numpy as np
+
 import mne
 from mne.datasets import sample
 from mne.preprocessing import ICA
 
 print(__doc__)
 
 # %%
+
 # Read and preprocess the data. Preprocessing consists of:
 #
 # - MEG channel selection
 # - 1-30 Hz band-pass filter
 
-data_path = sample.data_path()
-meg_path = data_path / "MEG" / "sample"
-raw_fname = meg_path / "sample_audvis_filt-0-40_raw.fif"
+# Load sample dataset
+data_path = Path(sample.data_path())
+raw_file = data_path / "MEG" / "sample" / "sample_audvis_raw.fif"
+raw = mne.io.read_raw_fif(raw_file, preload=True)
+raw.pick_types(meg=True, eeg=False, eog=True)
+raw.crop(0, 60)
 
-raw = mne.io.read_raw_fif(raw_fname).crop(0, 60).pick("meg").load_data()
+# %%
 
-reject = dict(mag=5e-12, grad=4000e-13)
-raw.filter(1, 30, fir_design="firwin")
+# Copy for clean and noisy
+raw_clean = raw.copy()
+raw_noisy = raw_clean.copy()
+noise = 1e-12 * np.random.randn(*raw_noisy._data.shape)
+raw_noisy._data += noise
 
+# Rejection thresholds
+reject_clean = dict(mag=5e-12, grad=4000e-13)
+reject_noisy = dict(mag=1e-11, grad=8000e-13)
 
 # %%
-# Define a function that runs ICA on the raw MEG data and plots the components
 
 
-def run_ica(method, fit_params=None):
+# Run ICA
+def run_ica(raw_input, method, fit_params=None, reject=None):
+    print(f"\nRunning ICA with: {method}")
     ica = ICA(
         n_components=20,
         method=method,
@@ -53,24 +75,82 @@ def run_ica(method, fit_params=None):
         random_state=0,
     )
     t0 = time()
-    ica.fit(raw, reject=reject)
+    ica.fit(raw_input, reject=reject)
     fit_time = time() - t0
-    title = f"ICA decomposition using {method} (took {fit_time:.1f}s)"
+    print(f"Fitting ICA took {fit_time:.1f}s.")
+
+    # Updated code with broken long line
+    title = (
+        f"ICA decomposition using {method} on "
+        f"{'noisy' if raw_input is raw_noisy else 'clean'} data\n"
+        f"(took {fit_time:.1f}s)"
+    )
     ica.plot_components(title=title)
 
+    return ica, fit_time
+
 
 # %%
-# FastICA
-run_ica("fastica")
+
+
+# Run all ICA methods
+def run_all_ica(raw_input, label, reject):
+    icas = {}
+    fit_times = {}
+    eog_components = {}
+    for method, params in [
+        ("fastica", None),
+        ("picard", None),
+        ("infomax", None),
+        ("infomax", {"extended": True}),
+    ]:
+        name = f"{method}" if not params else f"{method}_extended"
+        full_label = f"{label}_{name}"
+        ica, t = run_ica(raw_input, method, params, reject)
+        icas[full_label] = ica
+        fit_times[full_label] = t
+
+        eog_inds, _ = ica.find_bads_eog(raw_input, threshold=3.0, verbose="ERROR")
+        if eog_inds:
+            eog_components[full_label] = eog_inds[0]
+            print(f"{full_label}:Detected EOG comp at index {eog_inds[0]}")
+        else:
+            eog_components[full_label] = None
+            print(f"{full_label}: No EOG component detected")
+
+    return icas, fit_times, eog_components
+
 
 # %%
-# Picard
-run_ica("picard")
+
+
+# Run on both raw versions
+icas_clean, times_clean, eog_clean = run_all_ica(raw_clean, "clean", reject_clean)
+icas_noisy, times_noisy, eog_noisy = run_all_ica(raw_noisy, "noisy", reject_noisy)
+
+# Combine results
+icas = {**icas_clean, **icas_noisy}
+times = {**times_clean, **times_noisy}
+eog_comps = {**eog_clean, **eog_noisy}
 
 # %%
-# Infomax
-run_ica("infomax")
+
+# Clean EOG components for each algorithm (Column 1)
+for method in ["fastica", "picard", "infomax", "infomax_extended"]:
+    key = f"clean_{method}"
+    comp = eog_comps.get(key)
+    if comp is not None:
+        icas[key].plot_components(
+            picks=[comp], title=f"{key} - EOG Component (Clean Data)", show=True
+        )
 
 # %%
-# Extended Infomax
-run_ica("infomax", fit_params=dict(extended=True))
+
+# Noisy EOG components for each algorithm (Column 2)
+for method in ["fastica", "picard", "infomax", "infomax_extended"]:
+    key = f"noisy_{method}"
+    comp = eog_comps.get(key)
+    if comp is not None:
+        icas[key].plot_components(
+            picks=[comp], title=f"{key} - EOG Component (Noisy Data)", show=True
+        )