ENH: Add eSSS

doc/changes/latest.inc

@@ -107,6 +107,8 @@ Changelog
 
 - Add automatic T3 magnetometer detection and application of :meth:`mne.io.Raw.fix_mag_coil_types` to :func:`mne.preprocessing.maxwell_filter` by `Eric Larson`_
 
+- Add extended SSS (eSSS) support to :func:`mne.preprocessing.maxwell_filter` by `Eric Larson`_
+
 - Add ``'auto'`` option to :meth:`mne.preprocessing.ICA.find_bads_ecg` to automatically determine the threshold for CTPS method by `Yu-Han Luo`_
 
 - Add a ``notebook`` 3d backend for visualization in jupyter notebook with :func:`mne.viz.set_3d_backend` by `Guillaume Favelier`_

mne/chpi.py

@@ -517,7 +517,7 @@ def _setup_ext_proj(info, ext_order):
     ext = _sss_basis(
         dict(origin=(0., 0., 0.), int_order=0, ext_order=ext_order),
         mf_coils).T
-    out_removes = _regularize_out(0, 1, mag_or_fine)
+    out_removes = _regularize_out(0, 1, mag_or_fine, [])
     ext = ext[~np.in1d(np.arange(len(ext)), out_removes)]
     ext = linalg.orth(ext.T).T
     assert ext.shape[1] == len(meg_picks)

mne/fixes.py

@@ -20,6 +20,7 @@
 import warnings
 
 import numpy as np
+import scipy
 from scipy import linalg
 from scipy.linalg import LinAlgError
 
@@ -314,6 +315,22 @@ def _get_dpss():
     from numpy.fft import fft, ifft, fftfreq, rfft, irfft, rfftfreq, ifftshift
 
 
+###############################################################################
+# Orth with rcond argument (SciPy 1.1)
+
+if LooseVersion(scipy.__version__) >= '1.1':
+    from scipy.linalg import orth
+else:
+    def orth(A, rcond=None):  # noqa
+        u, s, vh = linalg.svd(A, full_matrices=False)
+        M, N = u.shape[0], vh.shape[1]
+        if rcond is None:
+            rcond = numpy.finfo(s.dtype).eps * max(M, N)
+        tol = np.amax(s) * rcond
+        num = np.sum(s > tol, dtype=int)
+        Q = u[:, :num]
+        return Q
+
 ###############################################################################
 # NumPy Generator (NumPy 1.17)
 

mne/proj.py

@@ -18,7 +18,6 @@
 from .forward import (is_fixed_orient, _subject_from_forward,
                       convert_forward_solution)
 from .source_estimate import _make_stc
-from .rank import _get_rank_sss
 
 
 def read_proj(fname):
@@ -82,8 +81,6 @@ def _compute_proj(data, info, n_grad, n_mag, n_eeg, desc_prefix,
 
     _check_option('meg', meg, ['separate', 'combined'])
     if meg == 'combined':
-        _get_rank_sss(info, msg='meg="combined" can only be used with '
-                      'Maxfiltered data', verbose=False)
         if n_grad != n_mag:
             raise ValueError('n_grad (%d) must be equal to n_mag (%d) when '
                              'using meg="combined"')

mne/tests/test_proj.py

@@ -376,8 +376,6 @@ def test_sss_proj():
     raw = read_raw_fif(raw_fname)
     raw.crop(0, 1.0).load_data().pick_types(meg=True, exclude=())
     raw.pick_channels(raw.ch_names[:51]).del_proj()
-    with pytest.raises(ValueError, match='can only be used with Maxfiltered'):
-        compute_proj_raw(raw, meg='combined')
     raw_sss = maxwell_filter(raw, int_order=5, ext_order=2)
     sss_rank = 21  # really low due to channel picking
     assert len(raw_sss.info['projs']) == 0

mne/utils/docs.py

@@ -522,6 +522,13 @@
     or :meth:`mne.io.Raw.append`, or separated during acquisition.
     To disable, provide an empty list.
 """
+docdict['maxwell_extended'] = """
+extended_proj : list
+    The empty-room projection vectors used to extend the external
+    SSS basis (i.e., use eSSS).
+
+    .. versionadded:: 0.21
+"""
 
 # Rank
 docdict['rank'] = """