diff --git a/doc/changes/devel/0000.enhancement.rst b/doc/changes/devel/0000.enhancement.rst new file mode 100644 index 00000000000..ff668050076 --- /dev/null +++ b/doc/changes/devel/0000.enhancement.rst @@ -0,0 +1 @@ +Add example :ref:`ex-ica-comp` comparing ICA algorithms on clean vs noisy MEG data, by :newcontrib:`Ganasekhar Kalla`. diff --git a/doc/changes/devel/0000.newfeature.rst b/doc/changes/devel/0000.newfeature.rst new file mode 100644 index 00000000000..7b9675e8dd1 --- /dev/null +++ b/doc/changes/devel/0000.newfeature.rst @@ -0,0 +1 @@ +Add new example :ref:`ex-ica-comp` that compares performance of multiple ICA algorithms (FastICA, Picard, Infomax, Extended Infomax) on clean vs noisy MEG data showing differences in runtime, component detection and EOG-related ICs, by `Ganasekhar Kalla`_. diff --git a/doc/changes/devel/13215.enhancement.rst b/doc/changes/devel/13215.enhancement.rst new file mode 100644 index 00000000000..ff668050076 --- /dev/null +++ b/doc/changes/devel/13215.enhancement.rst @@ -0,0 +1 @@ +Add example :ref:`ex-ica-comp` comparing ICA algorithms on clean vs noisy MEG data, by :newcontrib:`Ganasekhar Kalla`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 2e4027b82b7..5f2de74282c 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -96,6 +96,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 diff --git a/examples/preprocessing/ica_comparison.py b/examples/preprocessing/ica_comparison.py index d4246b80362..3cf083573ef 100644 --- a/examples/preprocessing/ica_comparison.py +++ b/examples/preprocessing/ica_comparison.py @@ -1,23 +1,34 @@ """ .. _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 + +# authors : Ganasekhar Kalla # # 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 @@ -25,26 +36,37 @@ 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 + )