Skip to content

Add MEG support#76

Draft
vferat wants to merge 27 commits into
mainfrom
meg
Draft

Add MEG support#76
vferat wants to merge 27 commits into
mainfrom
meg

Conversation

@vferat
Copy link
Copy Markdown
Collaborator

@vferat vferat commented Apr 15, 2026

No description provided.

neurotuning-personal and others added 27 commits February 28, 2026 12:51
Implement frequency-specific epoch sizing (MATLAB port) and refactor band fitting. Adds compute_epoch_sizes_per_band to compute per-band epoch durations from a cycles-per-band parameter, ensures even sample counts and marks bands too long for the data. Introduces Gedai.epoch_size_in_cycles, a new _fit_single_band helper that encapsulates the per-band SENSAI fitting logic, and refactors fit_epochs to use it. Adds _fit_raw_frequency_specific: when epoch_size_in_cycles is set, fit_raw now computes a broadband reference covariance, determines wavelet bands, computes per-band epoch durations, epochs/decomposes the raw data per band, skips or zeroes bands that are too short or below the low cutoff, and fits each band individually. Preserves original fixed-duration path when epoch_size_in_cycles is None and includes logging and small argument-check cleanups.
Introduce a fast frequency-specific (spectral) processing path that computes a single MODWT/MRA pass and then fits/applies per-band GEDAI filters. Key changes:

- Frequency-specific path no longer forces duration adjustment for broadband reference epochs (broadband_duration uses the user-supplied value).
- Retain 2^level duration adjustment only for the fixed-duration (original) path.
- Implement single MODWT (modwt) + MRA (modwtmra) over the full signal, with padding to the nearest multiple of 2^level when needed.
- Compute frequency bands directly and compute per-band epoch sizes; slice the band's MRA into frequency-appropriate epochs to call _fit_single_band, avoiding repeated MODWT passes.
- Add fast transform: _transform_raw_frequency_specific applies fitted per-band spatial filters directly in MRA space and reconstructs cleaned signal by summing bands.
- Add fit_transform_raw and _fit_transform_raw_frequency_specific for single-pass fit+transform that incrementally fits each band and accumulates the cleaned output to reduce memory footprint.
- Preserve behavior for bands that are too short or below the low-frequency cutoff by creating pass-through / zeroed fits.
- Minor API/logging adjustments and use of compute_epoch_sizes_per_band, _compute_refcov, and _process_single_epoch as appropriate.
- Rename example: opm_tutorial_gedai_full.py -> opm_tutorial_gedai_broadband.py and add a new opm_tutorial_gedai_spectral.py demonstrating the spectral GEDAI workflow.

These changes improve performance and memory usage for spectral GEDAI by avoiding O(n_bands) MODWT repetitions and enabling single-pass fit+transform for wavelet-based, frequency-specific epoching.
Introduce a fast analytical SENSAI scoring path and use it in the spectral fit path to avoid repeated GEVD/time-series reconstruction and speed up threshold search.

Changes:
- gedai/gedai/gedai.py: compute and cache per-epoch GEVD eigenvalues and eigenvectors, precompute reference subspace, remove EpochsArray creation for the fast path, and call _sensai_gridsearch_fast/_sensai_optimize_fast instead of the legacy epoch-based routines. Docstrings updated to note retained API args for compatibility.
- gedai/sensai/__init__.py: export new fast helpers (_sensai_gridsearch_fast, _sensai_optimize_fast, _sensai_score_fast) alongside existing API.
- gedai/sensai/sensai.py: add _sensai_score_fast which analytically computes signal/noise covariances from GEVD results (no time-series reconstruction), plus _sensai_gridsearch_fast and _sensai_optimize_fast that use the fast scorer (with parallel support). Legacy epoch-based _sensai_score/_sensai_gridsearch/_sensai_optimize are preserved with deprecation notes. Minor cleanup in subspace_angles and eigen<->sensai mapping helper spacing.

Motivation: the fast path reuses eigenvectors/eigenvalues and performs analytic covariance operations, significantly reducing cost of evaluating many thresholds (especially for low-frequency bands). Backwards compatibility maintained by keeping legacy epoch-based functions for the non-fast code paths.
Introduce signal_type support and automatic detection to differentiate EEG vs MEG workflows. Adds _detect_signal_type(info) and a signal_type parameter ("eeg"/"meg"/"auto") to Gedai; validates and propagates the resolved type through fit/transform paths. EEG now applies average referencing and uses fixed 3 PCs / 98th-percentile mapping; MEG uses adaptive reference-PC selection (≥85% variance), 4 PCs for SSI, and 99th-percentile mapping, with a warning if reference PCs are too few. sensai conversion and optimization functions (_sensai_to_eigen, _eigen_to_sensai, _sensai_optimize_fast) gain a percentile argument (default 95 for backward compatibility) so eigenvalue↔SENSAI mapping can match MATLAB behavior. Also update reference-eigenvector selection to use refcov_n_pc and thread signal_type through relevant internal methods. Finally, adjust the tutorial example noise_multiplier from 3.0 to 1.5.
Introduce a highpass_cutoff option to Gedai (default 0.1 Hz) with validation and a new internal _modwt_highpass() helper that uses MODWT/modwtmra to zero the approximation band and reconstruct a high-passed Raw copy. Apply this pre-processing step before GEDAI fitting/transforming. Improve MEG reference-PC handling for very small sensor arrays (adjust refcov_n_pc/ssi_n_pc and log info instead of warning). Update opm_tutorial_gedai_spectral.py: adjust filtering, apply notch before GEDAI, import the compare plotting helper, increase noise_multiplier, and add an interactive before/after comparison. Add a new example script source_power_spectrum_opm.py that computes PSDs and demonstrates GEDAI denoising, source localization, and band-wise visualization for VectorView and OPM data.
Introduce a preliminary_broadband_noise_multiplier option to Gedai (default 6.0) that runs a single-band broadband GEDAI pass on the full signal before the per-band wavelet analysis, mirroring the MATLAB GEDAI.m pipeline to strip large-amplitude artifacts. The ctor, docs and internal code paths were updated to run a temporary Gedai (wavelet_level=0) preliminary pass (disabled recursively) in both raw- and epoch-based flows when the multiplier is set. Also update the spectral tutorial to use noise_multiplier=3.0 (was 2.0).
Validate and support an 'auto' wavelet_level for Gedai. Adds input validation to require a non-negative int or 'auto', and implements _resolve_wavelet_level(sfreq) which computes the decomposition level when 'auto' is set using L = ceil(log2(sfreq / cutoff)) - 1 (fallback cutoff=0.5). Replaces direct uses of self.wavelet_level with the resolved level in band computation and preprocessing, and prints a brief info line showing the chosen level and its frequency band. Update example scripts to use wavelet_level='auto' (adjusted wavelet_low_cutoff and noise_multiplier in the tutorial) so the level is chosen automatically based on sampling rate.
Add parallelization and efficiency improvements to MODWT/MRA and surface per-band ENOVA reporting. _modwt.py: introduce _swt/_mra helpers, add joblib Parallel (loky) support and n_jobs parameter; implement an O(L) cascade for MRA reconstruction (fewer iSWT calls) and avoid GIL contention. gedai.py: call modwt/modwtmra with n_jobs=-1 at top-level and use n_jobs=1 in inner/already-parallel contexts; compute per-band ENOVA and preliminary broadband ENOVA, collect a band table and print a concise summary. transform.py: pass n_jobs=1 for epoch-level calls to prevent oversubscription. Also small f-string formatting and memory-release tweaks.
Correct wavelet MRA ordering and frequency-band mapping, plus small behavioral fixes and a demo.

- gedai/wavelet/_modwt.py: Rework MRA computation so bands follow [D1..DL, S_L]; compute details as A_j - A_{j-1} and place the smooth approximation in the last index.
- gedai/wavelet/transform.py: Make epochs_to_wavelet produce freq_bands in MRA order (D1..DL, S_L).
- gedai/gedai/gedai.py: Adjust usages to match new MRA ordering (zero approximation band at mra[-1]) and fix noisy debug prints / preliminary-broadband calls to use the actual noise_multiplier variable.
- source_power_spectrum_opm.py: Increase GEDAI noise_multiplier from 1.5 to 2.0 for several MEG calls.
- Added demo_caueeg_gedai.py and bundled CAUEEG.set test data to provide a quick spectral GEDAI denoising demo.

These changes fix incorrect band indexing that could produce wrong high-/low-frequency assignments and ensure consistent usage of the noise_multiplier parameter.
Switch demo to a spectral GEDAI configuration by enabling automatic wavelet level selection and adding wavelet_low_cutoff and epoch_size_in_cycles parameters. Simplify epoch-size logic to treat any band with fmin==0 as the approximation band (fixed 1s epoch). Add explicit handling to zero-out and log the approximation band in the GEDAI processing pipeline (both when building wavelet fits and when finalizing bands), ensuring the approximation band is excluded from spectral output. Minor message/format tweaks and preservation of existing low-frequency cutoff behavior.
Assert minimum 0.5 overlap across Gedai API (update docs and error messages) to reflect mathematical constraints of the squared-cosine cross-fade windowing. Restructure transform_raw to handle spectral vs broadband paths consistently, add timing and a Total metrics summary (Total SENSAI, Total ENOVA, elapsed time) printed after processing, and compute SENSAI using sig_ss - 1.0*noise_ss (matching MATLAB SENSAI_basic). Minor defensive fixes: guard attribute access on wavelets_fits and handle threshold=None when selecting the chosen run. In the demo, add preliminary_broadband_noise_multiplier=3.0 with a comment and remove the now-duplicated manual metrics/printing block since the library prints the summary.
Replace manual SENSAI averaging logic with a call to Gedai.score_sensai_basic to compute the total SENSAI score (using noise_multiplier=1.0) and fall back to 0.0 on errors. Add handling for reference_cov == "leadfield" to resolve and load the packaged leadfield .mat file via _compute_refcov. Also update demo_caueeg_gedai.py to set preliminary_broadband_noise_multiplier from 3.0 to 6.0.
Add a pure-Python bounded scalar minimizer and use it across SENSAI optimization, plus related API and behavior tweaks. Changes include:

- Added _minimize_scalar_bounded and SENSAI_fminbnd.py wrapper to provide a Brent-like bounded optimizer (removes reliance on scipy.optimize.minimize_scalar in the fast path).
- Updated gedai/sensai conversions (_sensai_to_eigen and _eigen_to_sensai) to use natural logs, an auto scale-factor for very small eigenvalues, and adjusted offsets to recover EEGLAB microvolt scaling.
- Replaced calls to scipy minimize in _sensai_optimize_fast with the new _minimize_scalar_bounded implementation.
- Adjusted Gedai behavior: coarsened sensai threshold step (0.25), changed epoch/window stepping to use full window (no overlap) in several places, and propagate/print mean ENOVA and noise_multiplier in status output.
- Updated MATLAB runner to use 'precomputed' refCOV argument and removed inline refCOV loading.
- Added test_opt.py to exercise gridsearch vs. optimizer and help validate the new optimizer and scaling.

These changes aim to remove the SciPy dependency for the bounded minimizer, improve numeric robustness for tiny eigenvalues, and simplify epoch stepping for faster deterministic runs. Please review for numerical equivalence and performance implications on existing pipelines.
Compute and choose PCA components dynamically in Gedai: for EEG default to 3 PCs; otherwise compute eigenvalues of the reference covariance and pick the number of PCs that capture 85% variance (bounded between 1 and n_channels-1). The selected ssi_n_pc and refcov_n_pc are printed. This logic was inserted in two mirrored locations where reference covariance is regularized.

Add interactive coregistration visualizations in source_power_spectrum_opm.py: show VectorView and OPM alignment figures, wait for user input to close each, and then continue. Also increase Gedai.fit_transform_raw noise_multiplier from 2.0 to 3.0 for grad, mag, and opm processing to apply stronger denoising.
Signed-off-by: Victor Férat <victor.ferat@live.Fr>
@vferat vferat marked this pull request as draft May 6, 2026 10:12
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants