Merge remote-tracking branch 'upstream/master' into refactor

* upstream/master:
  MRG, ENH: Add method to project onto max power ori (mne-tools#7883)
  WIP: Warn if untested NIRX device (mne-tools#7905)
  MRG, BUG: Fix bug with volume morph and subject_to!="fsaverage" (mne-tools#7896)
  MRG, MAINT: Clean up use of bool, float, int (mne-tools#7917)
  ENH: Better error message for incompatible Evoked objects (mne-tools#7910)
  try to fix nullcontext (mne-tools#7908)

doc/changes/latest.inc

@@ -27,12 +27,18 @@ Changelog
 
 - Add :meth:`mne.MixedSourceEstimate.surface` and :meth:`mne.MixedSourceEstimate.volume` methods to allow surface and volume extraction by `Eric Larson`_
 
+- Add :meth:`mne.VectorSourceEstimate.project` to project vector source estimates onto the direction of maximum source power by `Eric Larson`_
+
 - Add support to :func:`mne.extract_label_time_course` for vector-valued and volumetric source estimates by `Eric Larson`_
 
 - Add method :meth:`mne.VolSourceEstimate.in_label` by `Eric Larson`_
 
 - Add support for mixed source spaces to :func:`mne.compute_source_morph` by `Eric Larson`_
 
+- Add support for omitting the SDR step in volumetric morphing by passing ``n_iter_sdr=()`` to `mne.compute_source_morph` by `Eric Larson`_
+
+- Add support for passing a string argument to ``bg_img`` in `mne.viz.plot_volume_source_estimates` by `Eric Larson`_
+
 - Add support for providing the destination surface source space in the ``src_to`` argument of :func:`mne.compute_source_morph` by `Eric Larson`_
 
 - Add ``tol_kind`` option to :func:`mne.compute_rank` by `Eric Larson`_
@@ -89,6 +95,8 @@ Changelog
 
 - Add support for whitening and multiple channel types in :func:`mne.beamformer.make_dics` by `Marijn van Vliet`_
 
+- Add better error message when trying to save incompatible `~mne.Evoked` objects to the same file by `Eric Larson`_
+
 - Add support for loading complex numbers from mat files by `Thomas Hartmann`_
 
 - Add generic reader function :func:`mne.io.read_raw` that loads files based on their extensions (it wraps the underlying specific ``read_raw_xxx`` functions) by `Clemens Brunner`_
@@ -116,6 +124,10 @@ Bug
 
 - Fix bug with :func:`mne.compute_source_morph` when more than one volume source space was present (e.g., when using labels) where only the first label would be interpolated when ``mri_resolution=True`` by `Eric Larson`_
 
+- Fix bug with :func:`mne.compute_source_morph` when morphing to a volume source space when ``src_to`` is used and the destination subject is not ``fsaverage`` by `Eric Larson`_
+
+- Fix bug with :func:`mne.compute_source_morph` where outermost voxels in the destination source space could be errantly omitted by `Eric Larson`_
+
 - Fix bug with :func:`mne.minimum_norm.compute_source_psd_epochs` and :func:`mne.minimum_norm.source_band_induced_power` raised errors when ``method='eLORETA'`` by `Eric Larson`_
 
 - Fix bug with :func:`mne.minimum_norm.apply_inverse` where the explained variance did not work for complex data by `Eric Larson`_
@@ -234,6 +246,8 @@ API
 
 - The method ``stc_mixed.plot_surface`` for a :class:`mne.MixedSourceEstimate` has been deprecated in favor of :meth:`stc.surface().plot(...) <mne.MixedSourceEstimate.surface>` by `Eric Larson`_
 
+- The method ``stc.normal`` for :class:`mne.VectorSourceEstimate` has been deprecated in favor of :meth:`stc.project('nn', src) <mne.VectorSourceEstimate.project>` by `Eric Larson`_
+
 - Add ``use_dev_head_trans`` parameter to :func:`mne.preprocessing.annotate_movement` to allow choosing the device to head transform is used to define the fixed cHPI coordinates by `Luke Bloy`_
 
 - The function ``mne.channels.read_dig_captrack`` will be deprecated in version 0.22 in favor of :func:`mne.channels.read_dig_captrak` to correct the spelling error: "captraCK" -> "captraK", by `Stefan Appelhoff`_

examples/inverse/plot_morph_volume_stc.py

@@ -9,7 +9,8 @@
 :class:`mne.VolSourceEstimate` to a common reference space. We achieve this
 using :class:`mne.SourceMorph`. Pre-computed data will be morphed based on
 an affine transformation and a nonlinear registration method
-known as Symmetric Diffeomorphic Registration (SDR) by Avants et al. [1]_.
+known as Symmetric Diffeomorphic Registration (SDR) by
+:footcite:`AvantsEtAl2008`.
 
 Transformation is estimated from the subject's anatomical T1 weighted MRI
 (brain) to `FreeSurfer's 'fsaverage' T1 weighted MRI (brain)
@@ -18,13 +19,6 @@
 Afterwards the transformation will be applied to the volumetric source
 estimate. The result will be plotted, showing the fsaverage T1 weighted
 anatomical MRI, overlaid with the morphed volumetric source estimate.
-
-References
-----------
-.. [1] Avants, B. B., Epstein, C. L., Grossman, M., & Gee, J. C. (2009).
-       Symmetric Diffeomorphic Image Registration with Cross- Correlation:
-       Evaluating Automated Labeling of Elderly and Neurodegenerative
-       Brain, 12(1), 26-41.
 """
 # Author: Tommy Clausner <tommy.clausner@gmail.com>
 #
@@ -153,3 +147,6 @@
 #
 #     >>> morph.apply(stc)
 #
+# References
+# ----------
+# .. footbibliography::

examples/inverse/plot_vector_mne_solution.py

@@ -21,6 +21,7 @@
 #
 # License: BSD (3-clause)
 
+import numpy as np
 import mne
 from mne.datasets import sample
 from mne.minimum_norm import read_inverse_operator, apply_inverse
@@ -52,6 +53,21 @@
 brain = stc.plot(
     initial_time=peak_time, hemi='lh', subjects_dir=subjects_dir)
 
+###############################################################################
+# Plot the activation in the direction of maximal power for this data:
+
+stc_max, directions = stc.project('pca', src=inv['src'])
+# These directions must by design be close to the normals because this
+# inverse was computed with loose=0.2:
+print('Absolute cosine similarity between source normals and directions: '
+      f'{np.abs(np.sum(directions * inv["source_nn"][2::3], axis=-1)).mean()}')
+brain_max = stc_max.plot(
+    initial_time=peak_time, hemi='lh', subjects_dir=subjects_dir,
+    time_label='Max power')
+brain_normal = stc.project('normal', inv['src'])[0].plot(
+    initial_time=peak_time, hemi='lh', subjects_dir=subjects_dir,
+    time_label='Normal')
+
 ###############################################################################
 # You can also do this with a fixed-orientation inverse. It looks a lot like
 # the result above because the ``loose=0.2`` orientation constraint keeps

examples/simulation/plot_simulate_evoked_data.py

@@ -46,7 +46,7 @@
 ###############################################################################
 # Generate source time courses from 2 dipoles and the correspond evoked data
 
-times = np.arange(300, dtype=np.float) / raw.info['sfreq'] - 0.1
+times = np.arange(300, dtype=np.float64) / raw.info['sfreq'] - 0.1
 rng = np.random.RandomState(42)
 
 

examples/simulation/plot_simulated_raw_data_using_subject_anatomy.py

@@ -143,7 +143,7 @@ def data_fun(times, latency, duration):
 # event, the second is not used. The third one is the event id, which is
 # different for each of the 4 areas.
 
-times = np.arange(150, dtype=np.float) / info['sfreq']
+times = np.arange(150, dtype=np.float64) / info['sfreq']
 duration = 0.03
 rng = np.random.RandomState(7)
 source_simulator = mne.simulation.SourceSimulator(src, tstep=tstep)

examples/time_frequency/plot_time_frequency_simulated.py

@@ -43,7 +43,7 @@
 noise = rng.randn(n_epochs, len(ch_names), n_times)
 
 # Add a 50 Hz sinusoidal burst to the noise and ramp it.
-t = np.arange(n_times, dtype=np.float) / sfreq
+t = np.arange(n_times, dtype=np.float64) / sfreq
 signal = np.sin(np.pi * 2. * 50. * t)  # 50 Hz sinusoid signal
 signal[np.logical_or(t < 0.45, t > 0.55)] = 0.  # Hard windowing
 on_time = np.logical_and(t >= 0.45, t <= 0.55)

mne/beamformer/tests/test_dics.py

@@ -284,7 +284,7 @@ def test_make_dics(tmpdir, _load_forward, idx, whiten):
     assert isinstance(filters, Beamformer)
     assert isinstance(filters_read, Beamformer)
     for key in ['tmin', 'tmax']:  # deal with strictness of object_diff
-        setattr(filters['csd'], key, np.float(getattr(filters['csd'], key)))
+        setattr(filters['csd'], key, np.float64(getattr(filters['csd'], key)))
     assert object_diff(filters, filters_read) == ''
 
 

mne/bem.py

@@ -1341,7 +1341,7 @@ def _read_bem_surface(fid, this, def_coord_frame, s_id=None):
     if tag is None:
         raise ValueError('Vertex data not found')
 
-    res['rr'] = tag.data.astype(np.float)  # XXX : double because of mayavi bug
+    res['rr'] = tag.data.astype(np.float64)  # XXX : double because of mayavi
     if res['rr'].shape[0] != res['np']:
         raise ValueError('Vertex information is incorrect')
 

mne/channels/channels.py

@@ -373,7 +373,7 @@ def _set_channel_positions(self, pos, names):
         if len(pos) != len(names):
             raise ValueError('Number of channel positions not equal to '
                              'the number of names given.')
-        pos = np.asarray(pos, dtype=np.float)
+        pos = np.asarray(pos, dtype=np.float64)
         if pos.shape[-1] != 3 or pos.ndim != 2:
             msg = ('Channel positions must have the shape (n_points, 3) '
                    'not %s.' % (pos.shape,))

mne/channels/interpolation.py

@@ -136,8 +136,8 @@ def _interpolate_bads_eeg(inst, origin, verbose=None):
     inst : mne.io.Raw, mne.Epochs or mne.Evoked
         The data to interpolate. Must be preloaded.
     """
-    bads_idx = np.zeros(len(inst.ch_names), dtype=np.bool)
-    goods_idx = np.zeros(len(inst.ch_names), dtype=np.bool)
+    bads_idx = np.zeros(len(inst.ch_names), dtype=bool)
+    goods_idx = np.zeros(len(inst.ch_names), dtype=bool)
 
     picks = pick_types(inst.info, meg=False, eeg=True, exclude=[])
     inst.info._check_consistency()

mne/channels/layout.py

@@ -126,7 +126,7 @@ def _read_lout(fname):
                 name = chkind + ' ' + nb
             else:
                 cid, x, y, dx, dy, name = splits
-            pos.append(np.array([x, y, dx, dy], dtype=np.float))
+            pos.append(np.array([x, y, dx, dy], dtype=np.float64))
             names.append(name)
             ids.append(int(cid))
 
@@ -147,7 +147,7 @@ def _read_lay(fname):
                 name = chkind + ' ' + nb
             else:
                 cid, x, y, dx, dy, name = splits
-            pos.append(np.array([x, y, dx, dy], dtype=np.float))
+            pos.append(np.array([x, y, dx, dy], dtype=np.float64))
             names.append(name)
             ids.append(int(cid))
 

mne/chpi.py

@@ -557,7 +557,7 @@ def _fit_chpi_amplitudes(raw, time_sl, hpi):
             # loads hpi_stim channel
             chpi_data = raw[hpi['hpi_pick'], time_sl][0]
 
-        ons = (np.round(chpi_data).astype(np.int) &
+        ons = (np.round(chpi_data).astype(np.int64) &
                hpi['on'][:, np.newaxis]).astype(bool)
         n_on = ons.all(axis=-1).sum(axis=0)
         if not (n_on >= 3).all():

mne/conftest.py

@@ -88,6 +88,11 @@ def pytest_configure(config):
     ignore:.*trait.*handler.*deprecated.*:DeprecationWarning
     ignore:.*rich_compare.*metadata.*deprecated.*:DeprecationWarning
     ignore:.*In future, it will be an error for 'np.bool_'.*:DeprecationWarning
+    ignore:.*`np.bool` is a deprecated alias.*:DeprecationWarning
+    ignore:.*`np.int` is a deprecated alias.*:DeprecationWarning
+    ignore:.*`np.float` is a deprecated alias.*:DeprecationWarning
+    ignore:.*`np.object` is a deprecated alias.*:DeprecationWarning
+    ignore:.*`np.long` is a deprecated alias:DeprecationWarning
     ignore:.*Converting `np\.character` to a dtype is deprecated.*:DeprecationWarning
     ignore:.*sphinx\.util\.smartypants is deprecated.*:
     ignore:.*pandas\.util\.testing is deprecated.*:

mne/connectivity/effective.py

@@ -133,7 +133,7 @@ def phase_slope_index(data, indices=None, sfreq=2 * np.pi,
     # allocate space for output
     out_shape = list(cohy.shape)
     out_shape[freq_dim] = n_bands
-    psi = np.zeros(out_shape, dtype=np.float)
+    psi = np.zeros(out_shape, dtype=np.float64)
 
     # allocate accumulator
     acc_shape = copy.copy(out_shape)

mne/connectivity/spectral.py

@@ -980,7 +980,7 @@ def _prepare_connectivity(epoch_block, tmin, tmax, fmin, fmax, sfreq, indices,
             raise ValueError('define frequencies of interest using '
                              'cwt_freqs')
         else:
-            cwt_freqs = cwt_freqs.astype(np.float)
+            cwt_freqs = cwt_freqs.astype(np.float64)
         if any(cwt_freqs > (sfreq / 2.)):
             raise ValueError('entries in cwt_freqs cannot be '
                              'larger than Nyquist (sfreq / 2)')
@@ -1003,7 +1003,7 @@ def _prepare_connectivity(epoch_block, tmin, tmax, fmin, fmax, sfreq, indices,
                                   5. / np.min(fmin), five_cycle_freq))
 
     # create a frequency mask for all bands
-    freq_mask = np.zeros(len(freqs_all), dtype=np.bool)
+    freq_mask = np.zeros(len(freqs_all), dtype=bool)
     for f_lower, f_upper in zip(fmin, fmax):
         freq_mask |= ((freqs_all >= f_lower) & (freqs_all <= f_upper))
 

mne/cov.py

@@ -846,8 +846,8 @@ def _unpack_epochs(epochs):
         # prepare mean covs
         n_epoch_types = len(epochs)
         data_mean = [0] * n_epoch_types
-        n_samples = np.zeros(n_epoch_types, dtype=np.int)
-        n_epochs = np.zeros(n_epoch_types, dtype=np.int)
+        n_samples = np.zeros(n_epoch_types, dtype=np.int64)
+        n_epochs = np.zeros(n_epoch_types, dtype=np.int64)
 
         for ii, epochs_t in enumerate(epochs):
 
@@ -1759,7 +1759,10 @@ def compute_whitener(noise_cov, info=None, picks=None, rank=None,
     nzero = (eig > 0)
     eig[~nzero] = 0.  # get rid of numerical noise (negative) ones
 
-    dtype = np.complex if noise_cov['eigvec'].dtype == np.complex else np.float
+    if noise_cov['eigvec'].dtype.kind == 'c':
+        dtype = np.complex128
+    else:
+        dtype = np.float64
     W = np.zeros((n_chan, 1), dtype)
     W[nzero, 0] = 1.0 / np.sqrt(eig[nzero])
     #   Rows of eigvec are the eigenvectors
@@ -1966,7 +1969,7 @@ def _write_cov(fid, cov):
     else:
         # Store only lower part of covariance matrix
         dim = cov['dim']
-        mask = np.tril(np.ones((dim, dim), dtype=np.bool)) > 0
+        mask = np.tril(np.ones((dim, dim), dtype=bool)) > 0
         vals = cov['data'][mask].ravel()
         write_double(fid, FIFF.FIFF_MNE_COV, vals)
 

mne/datasets/utils.py

@@ -239,7 +239,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.92', misc='0.6')
+    releases = dict(testing='0.93', misc='0.6')
     # And also update the "md5_hashes['testing']" variable below.
 
     # To update any other dataset, update the data archive itself (upload
@@ -326,7 +326,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='12b75d1cb7df9dfb4ad73ed82f61094f',
         somato='ea825966c0a1e9b2f84e3826c5500161',
         spm='9f43f67150e3b694b523a21eb929ea75',
-        testing='42daafd1b882da2ef041de860ca6e771',
+        testing='2eb998a0893a28faedd583973c70b517',
         multimodal='26ec847ae9ab80f58f204d09e2c08367',
         fnirs_motor='c4935d19ddab35422a69f3326a01fef8',
         opm='370ad1dcfd5c47e029e692c85358a374',

mne/decoding/base.py

@@ -224,7 +224,7 @@ def _set_cv(cv, estimator=None, X=None, y=None):
     # Setup CV
     from sklearn import model_selection as models
     from sklearn.model_selection import (check_cv, StratifiedKFold, KFold)
-    if isinstance(cv, (int, np.int)):
+    if isinstance(cv, (int, np.int64)):
         XFold = StratifiedKFold if est_is_classifier else KFold
         cv = XFold(n_splits=cv)
     elif isinstance(cv, str):

mne/decoding/tests/test_receptive_field.py

@@ -108,7 +108,7 @@ def test_time_delay():
             _delay_time_series(X, tmin, tmax, sfreq=[1])
         # Delays must be int/float
         with pytest.raises(TypeError, match='.*complex.*'):
-            _delay_time_series(X, np.complex(tmin), tmax, 1)
+            _delay_time_series(X, np.complex128(tmin), tmax, 1)
         # Make sure swapaxes works
         start, stop = int(round(tmin * isfreq)), int(round(tmax * isfreq)) + 1
         n_delays = stop - start

mne/dipole.py

@@ -1306,7 +1306,7 @@ def fit_dipole(evoked, cov, bem, trans=None, min_dist=5., n_jobs=1,
     # cov = prepare_noise_cov(cov, info_nb, info_nb['ch_names'], verbose=False)
     # nzero = (cov['eig'] > 0)
     # n_chan = len(info_nb['ch_names'])
-    # whitener = np.zeros((n_chan, n_chan), dtype=np.float)
+    # whitener = np.zeros((n_chan, n_chan), dtype=np.float64)
     # whitener[nzero, nzero] = 1.0 / np.sqrt(cov['eig'][nzero])
     # whitener = np.dot(whitener, cov['eigvec'])
 

mne/epochs.py

@@ -904,7 +904,7 @@ def subtract_evoked(self, evoked=None):
         else:
             if self._offset is None:
                 self._offset = np.zeros((len(self.ch_names), len(self.times)),
-                                        dtype=np.float)
+                                        dtype=np.float64)
             self._offset[ep_picks] -= evoked.data[picks]
         logger.info('[done]')
 
@@ -2900,35 +2900,37 @@ def add_channels_epochs(epochs_list, verbose=None):
     return epochs
 
 
-def _compare_epochs_infos(info1, info2, ind):
+def _compare_epochs_infos(info1, info2, name):
     """Compare infos."""
+    if not isinstance(name, str):  # passed epochs index
+        name = f'epochs[{name:d}]'
     info1._check_consistency()
     info2._check_consistency()
     if info1['nchan'] != info2['nchan']:
-        raise ValueError('epochs[%d][\'info\'][\'nchan\'] must match' % ind)
+        raise ValueError(f'{name}.info[\'nchan\'] must match')
     if info1['bads'] != info2['bads']:
-        raise ValueError('epochs[%d][\'info\'][\'bads\'] must match' % ind)
+        raise ValueError(f'{name}.info[\'bads\'] must match')
     if info1['sfreq'] != info2['sfreq']:
-        raise ValueError('epochs[%d][\'info\'][\'sfreq\'] must match' % ind)
+        raise ValueError(f'{name}.info[\'sfreq\'] must match')
     if set(info1['ch_names']) != set(info2['ch_names']):
-        raise ValueError('epochs[%d][\'info\'][\'ch_names\'] must match' % ind)
+        raise ValueError(f'{name}.info[\'ch_names\'] must match')
     if len(info2['projs']) != len(info1['projs']):
-        raise ValueError('SSP projectors in epochs files must be the same')
+        raise ValueError(f'SSP projectors in {name} must be the same')
     if any(not _proj_equal(p1, p2) for p1, p2 in
            zip(info2['projs'], info1['projs'])):
-        raise ValueError('SSP projectors in epochs files must be the same')
+        raise ValueError(f'SSP projectors in {name} must be the same')
     if (info1['dev_head_t'] is None) != (info2['dev_head_t'] is None) or \
             (info1['dev_head_t'] is not None and not
              np.allclose(info1['dev_head_t']['trans'],
                          info2['dev_head_t']['trans'], rtol=1e-6)):
-        raise ValueError('epochs[%d][\'info\'][\'dev_head_t\'] must match. '
-                         'The epochs probably come from different runs, and '
+        raise ValueError(f'{name}.info[\'dev_head_t\'] must match. The '
+                         'instances probably come from different runs, and '
                          'are therefore associated with different head '
                          'positions. Manually change info[\'dev_head_t\'] to '
                          'avoid this message but beware that this means the '
                          'MEG sensors will not be properly spatially aligned. '
                          'See mne.preprocessing.maxwell_filter to realign the '
-                         'runs to a common head position.' % ind)
+                         'runs to a common head position.')
 
 
 def _concatenate_epochs(epochs_list, with_data=True, add_offset=True):