diff --git a/doc/changes/devel/12326.other.rst b/doc/changes/devel/12326.other.rst new file mode 100644 index 00000000000..b8f2966bbf9 --- /dev/null +++ b/doc/changes/devel/12326.other.rst @@ -0,0 +1 @@ +Updated the text in the preprocessing tutorial to use :class:`mne.io.Raw.pick()` instead of the legacy :class:`mne.io.Raw.pick_types()`, by :newcontrib:`btkcodedev`. diff --git a/doc/changes/devel/12371.newfeature.rst b/doc/changes/devel/12371.newfeature.rst new file mode 100644 index 00000000000..4d28ff1f5ce --- /dev/null +++ b/doc/changes/devel/12371.newfeature.rst @@ -0,0 +1 @@ +Speed up :func:`mne.io.read_raw_neuralynx` on large datasets with many gaps, by `Kristijan Armeni`_. \ No newline at end of file diff --git a/doc/changes/devel/12383.newfeature.rst b/doc/changes/devel/12383.newfeature.rst new file mode 100644 index 00000000000..f896572eb93 --- /dev/null +++ b/doc/changes/devel/12383.newfeature.rst @@ -0,0 +1 @@ +Add ability to detect minima peaks found in :class:`mne.Evoked` if data is all positive and maxima if data is all negative. \ No newline at end of file diff --git a/doc/changes/devel/12389.bugfix.rst b/doc/changes/devel/12389.bugfix.rst new file mode 100644 index 00000000000..85892df97a8 --- /dev/null +++ b/doc/changes/devel/12389.bugfix.rst @@ -0,0 +1 @@ +Fix bug where :func:`mne.preprocessing.regress_artifact` projection check was not specific to the channels being processed, by `Eric Larson`_. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 811029ddaa7..0389f75e83e 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -68,6 +68,8 @@ .. _Bruno Nicenboim: https://bnicenboim.github.io +.. _btkcodedev: https://github.com/btkcodedev + .. _buildqa: https://github.com/buildqa .. _Carlos de la Torre-Ortiz: https://ctorre.me diff --git a/mne/_fiff/pick.py b/mne/_fiff/pick.py index 2722e91f517..9e2e369ab71 100644 --- a/mne/_fiff/pick.py +++ b/mne/_fiff/pick.py @@ -649,7 +649,8 @@ def pick_info(info, sel=(), copy=True, verbose=None): return info elif len(sel) == 0: raise ValueError("No channels match the selection.") - n_unique = len(np.unique(np.arange(len(info["ch_names"]))[sel])) + ch_set = set(info["ch_names"][k] for k in sel) + n_unique = len(ch_set) if n_unique != len(sel): raise ValueError( "Found %d / %d unique names, sel is not unique" % (n_unique, len(sel)) @@ -687,6 +688,15 @@ def pick_info(info, sel=(), copy=True, verbose=None): if info.get("custom_ref_applied", False) and not _electrode_types(info): with info._unlock(): info["custom_ref_applied"] = FIFF.FIFFV_MNE_CUSTOM_REF_OFF + # remove unused projectors + if info.get("projs", False): + projs = list() + for p in info["projs"]: + if any(ch_name in ch_set for ch_name in p["data"]["col_names"]): + projs.append(p) + if len(projs) != len(info["projs"]): + with info._unlock(): + info["projs"] = projs info._check_consistency() return info diff --git a/mne/_fiff/tests/test_pick.py b/mne/_fiff/tests/test_pick.py index 5494093cd23..fa503a04ab3 100644 --- a/mne/_fiff/tests/test_pick.py +++ b/mne/_fiff/tests/test_pick.py @@ -558,11 +558,17 @@ def test_clean_info_bads(): # simulate the bad channels raw.info["bads"] = eeg_bad_ch + meg_bad_ch + assert len(raw.info["projs"]) == 3 + raw.set_eeg_reference(projection=True) + assert len(raw.info["projs"]) == 4 + # simulate the call to pick_info excluding the bad eeg channels info_eeg = pick_info(raw.info, picks_eeg) + assert len(info_eeg["projs"]) == 1 # simulate the call to pick_info excluding the bad meg channels info_meg = pick_info(raw.info, picks_meg) + assert len(info_meg["projs"]) == 3 assert info_eeg["bads"] == eeg_bad_ch assert info_meg["bads"] == meg_bad_ch diff --git a/mne/datasets/config.py b/mne/datasets/config.py index b7780778f24..238b61998d6 100644 --- a/mne/datasets/config.py +++ b/mne/datasets/config.py @@ -87,8 +87,12 @@ # To update the `testing` or `misc` datasets, push or merge commits to their # respective repos, and make a new release of the dataset on GitHub. Then # update the checksum in the MNE_DATASETS dict below, and change version -# here: ↓↓↓↓↓ ↓↓↓ -RELEASES = dict(testing="0.151", misc="0.27") +# here: ↓↓↓↓↓↓↓↓ +RELEASES = dict( + testing="0.151", + misc="0.27", + phantom_kit="0.2", +) TESTING_VERSIONED = f'mne-testing-data-{RELEASES["testing"]}' MISC_VERSIONED = f'mne-misc-data-{RELEASES["misc"]}' @@ -176,9 +180,9 @@ ) MNE_DATASETS["phantom_kit"] = dict( - archive_name="MNE-phantom-KIT-24bit.zip", - hash="md5:CAF82EE978DD473C7DE6C1034D9CCD45", - url="https://osf.io/download/svnt3/", + archive_name="MNE-phantom-KIT-data.tar.gz", + hash="md5:7bfdf40bbeaf17a66c99c695640e0740", + url="https://osf.io/fb6ya/download?version=1", folder_name="MNE-phantom-KIT-data", config_key="MNE_DATASETS_PHANTOM_KIT_PATH", ) diff --git a/mne/datasets/phantom_kit/phantom_kit.py b/mne/datasets/phantom_kit/phantom_kit.py index a4eac6c4a50..d57ca875f2c 100644 --- a/mne/datasets/phantom_kit/phantom_kit.py +++ b/mne/datasets/phantom_kit/phantom_kit.py @@ -10,7 +10,7 @@ def data_path( ): # noqa: D103 return _download_mne_dataset( name="phantom_kit", - processor="unzip", + processor="untar", path=path, force_update=force_update, update_path=update_path, diff --git a/mne/evoked.py b/mne/evoked.py index c988368c314..1f694f7c11b 100644 --- a/mne/evoked.py +++ b/mne/evoked.py @@ -914,6 +914,8 @@ def get_peak( time_as_index=False, merge_grads=False, return_amplitude=False, + *, + strict=True, ): """Get location and latency of peak amplitude. @@ -941,6 +943,12 @@ def get_peak( If True, return also the amplitude at the maximum response. .. versionadded:: 0.16 + strict : bool + If True, raise an error if values are all positive when detecting + a minimum (mode='neg'), or all negative when detecting a maximum + (mode='pos'). Defaults to True. + + .. versionadded:: 1.7 Returns ------- @@ -1032,7 +1040,14 @@ def get_peak( data, _ = _merge_ch_data(data, ch_type, []) ch_names = [ch_name[:-1] + "X" for ch_name in ch_names[::2]] - ch_idx, time_idx, max_amp = _get_peak(data, self.times, tmin, tmax, mode) + ch_idx, time_idx, max_amp = _get_peak( + data, + self.times, + tmin, + tmax, + mode, + strict=strict, + ) out = (ch_names[ch_idx], time_idx if time_as_index else self.times[time_idx]) @@ -1949,7 +1964,7 @@ def _write_evokeds(fname, evoked, check=True, *, on_mismatch="raise", overwrite= end_block(fid, FIFF.FIFFB_MEAS) -def _get_peak(data, times, tmin=None, tmax=None, mode="abs"): +def _get_peak(data, times, tmin=None, tmax=None, mode="abs", *, strict=True): """Get feature-index and time of maximum signal from 2D array. Note. This is a 'getter', not a 'finder'. For non-evoked type @@ -1970,6 +1985,10 @@ def _get_peak(data, times, tmin=None, tmax=None, mode="abs"): values will be considered. If 'neg' only negative values will be considered. If 'abs' absolute values will be considered. Defaults to 'abs'. + strict : bool + If True, raise an error if values are all positive when detecting + a minimum (mode='neg'), or all negative when detecting a maximum + (mode='pos'). Defaults to True. Returns ------- @@ -2008,12 +2027,12 @@ def _get_peak(data, times, tmin=None, tmax=None, mode="abs"): maxfun = np.argmax if mode == "pos": - if not np.any(data[~mask] > 0): + if strict and not np.any(data[~mask] > 0): raise ValueError( "No positive values encountered. Cannot " "operate in pos mode." ) elif mode == "neg": - if not np.any(data[~mask] < 0): + if strict and not np.any(data[~mask] < 0): raise ValueError( "No negative values encountered. Cannot " "operate in neg mode." ) diff --git a/mne/io/neuralynx/neuralynx.py b/mne/io/neuralynx/neuralynx.py index 55de7579d67..ab768d57b13 100644 --- a/mne/io/neuralynx/neuralynx.py +++ b/mne/io/neuralynx/neuralynx.py @@ -223,11 +223,8 @@ def __init__( [np.full(shape=(n,), fill_value=i) for i, n in enumerate(sizes_sorted)] ) - # construct Annotations() - gap_seg_ids = np.unique(sample2segment)[gap_indicator] - gap_start_ids = np.array( - [np.where(sample2segment == seg_id)[0][0] for seg_id in gap_seg_ids] - ) + # get the start sample index for each gap segment () + gap_start_ids = np.cumsum(np.hstack([[0], sizes_sorted[:-1]]))[gap_indicator] # recreate time axis for gap annotations mne_times = np.arange(0, len(sample2segment)) / info["sfreq"] diff --git a/mne/preprocessing/_regress.py b/mne/preprocessing/_regress.py index 31a842f7d4f..260796a221d 100644 --- a/mne/preprocessing/_regress.py +++ b/mne/preprocessing/_regress.py @@ -6,7 +6,7 @@ import numpy as np -from .._fiff.pick import _picks_to_idx +from .._fiff.pick import _picks_to_idx, pick_info from ..defaults import _BORDER_DEFAULT, _EXTRAPOLATE_DEFAULT, _INTERPOLATION_DEFAULT from ..epochs import BaseEpochs from ..evoked import Evoked @@ -178,9 +178,7 @@ def fit(self, inst): reference (see :func:`mne.set_eeg_reference`) before performing EOG regression. """ - self._check_inst(inst) - picks = _picks_to_idx(inst.info, self.picks, none="data", exclude=self.exclude) - picks_artifact = _picks_to_idx(inst.info, self.picks_artifact) + picks, picks_artifact = self._check_inst(inst) # Calculate regression coefficients. Add a row of ones to also fit the # intercept. @@ -232,9 +230,7 @@ def apply(self, inst, copy=True): """ if copy: inst = inst.copy() - self._check_inst(inst) - picks = _picks_to_idx(inst.info, self.picks, none="data", exclude=self.exclude) - picks_artifact = _picks_to_idx(inst.info, self.picks_artifact) + picks, picks_artifact = self._check_inst(inst) # Check that the channels are compatible with the regression weights. ref_picks = _picks_to_idx( @@ -324,19 +320,25 @@ def _check_inst(self, inst): _validate_type( inst, (BaseRaw, BaseEpochs, Evoked), "inst", "Raw, Epochs, Evoked" ) - if _needs_eeg_average_ref_proj(inst.info): + picks = _picks_to_idx(inst.info, self.picks, none="data", exclude=self.exclude) + picks_artifact = _picks_to_idx(inst.info, self.picks_artifact) + all_picks = np.unique(np.concatenate([picks, picks_artifact])) + use_info = pick_info(inst.info, all_picks) + del all_picks + if _needs_eeg_average_ref_proj(use_info): raise RuntimeError( - "No reference for the EEG channels has been " - "set. Use inst.set_eeg_reference() to do so." + "No average reference for the EEG channels has been " + "set. Use inst.set_eeg_reference(projection=True) to do so." ) if self.proj and not inst.proj: inst.apply_proj() - if not inst.proj and len(inst.info.get("projs", [])) > 0: + if not inst.proj and len(use_info.get("projs", [])) > 0: raise RuntimeError( "Projections need to be applied before " "regression can be performed. Use the " ".apply_proj() method to do so." ) + return picks, picks_artifact def __repr__(self): """Produce a string representation of this object.""" diff --git a/mne/preprocessing/tests/test_regress.py b/mne/preprocessing/tests/test_regress.py index 8050b6bebf7..48d960e0464 100644 --- a/mne/preprocessing/tests/test_regress.py +++ b/mne/preprocessing/tests/test_regress.py @@ -42,6 +42,19 @@ def test_regress_artifact(): epochs, betas = regress_artifact(epochs, picks="eog", picks_artifact="eog") assert np.ptp(epochs.get_data("eog")) < 1e-15 # constant value assert_allclose(betas, 1) + # proj should only be required of channels being processed + raw = read_raw_fif(raw_fname).crop(0, 1).load_data() + raw.del_proj() + raw.set_eeg_reference(projection=True) + model = EOGRegression(proj=False, picks="meg", picks_artifact="eog") + model.fit(raw) + model.apply(raw) + model = EOGRegression(proj=False, picks="eeg", picks_artifact="eog") + with pytest.raises(RuntimeError, match="Projections need to be applied"): + model.fit(raw) + raw.del_proj() + with pytest.raises(RuntimeError, match="No average reference for the EEG"): + model.fit(raw) @testing.requires_testing_data diff --git a/mne/tests/test_evoked.py b/mne/tests/test_evoked.py index 2c5f064606d..31110596be6 100644 --- a/mne/tests/test_evoked.py +++ b/mne/tests/test_evoked.py @@ -589,6 +589,24 @@ def test_get_peak(): with pytest.raises(ValueError, match="No positive values"): evoked_all_neg.get_peak(mode="pos") + # Test finding minimum and maximum values + evoked_all_neg_outlier = evoked_all_neg.copy() + evoked_all_pos_outlier = evoked_all_pos.copy() + + # Add an outlier to the data + evoked_all_neg_outlier.data[0, 15] = -1e-20 + evoked_all_pos_outlier.data[0, 15] = 1e-20 + + ch_name, time_idx, max_amp = evoked_all_neg_outlier.get_peak( + mode="pos", return_amplitude=True, strict=False + ) + assert max_amp == -1e-20 + + ch_name, time_idx, min_amp = evoked_all_pos_outlier.get_peak( + mode="neg", return_amplitude=True, strict=False + ) + assert min_amp == 1e-20 + # Test interaction between `mode` and `tmin` / `tmax` # For the test, create an Evoked where half of the values are negative # and the rest is positive diff --git a/mne/viz/misc.py b/mne/viz/misc.py index 3d8a9469620..c9be85f8f9e 100644 --- a/mne/viz/misc.py +++ b/mne/viz/misc.py @@ -130,7 +130,7 @@ def plot_cov( fig_cov : instance of matplotlib.figure.Figure The covariance plot. fig_svd : instance of matplotlib.figure.Figure | None - The SVD spectra plot of the covariance. + The SVD plot of the covariance (i.e., the eigenvalues or "matrix spectrum"). See Also -------- diff --git a/tutorials/inverse/95_phantom_KIT.py b/tutorials/inverse/95_phantom_KIT.py index 6a07658e13a..75e0025a9c2 100644 --- a/tutorials/inverse/95_phantom_KIT.py +++ b/tutorials/inverse/95_phantom_KIT.py @@ -15,9 +15,8 @@ # Copyright the MNE-Python contributors. # %% -import matplotlib.pyplot as plt +import mne_bids import numpy as np -from scipy.signal import find_peaks import mne @@ -25,14 +24,33 @@ actual_pos, actual_ori = mne.dipole.get_phantom_dipoles("oyama") actual_pos, actual_ori = actual_pos[:49], actual_ori[:49] # only 49 of 50 dipoles -raw = mne.io.read_raw_kit(data_path / "002_phantom_11Hz_100uA.con") -# cut from ~800 to ~300s for speed, and also at convenient dip stim boundaries -# chosen by examining MISC 017 by eye. -raw.crop(11.5, 302.9).load_data() -raw.filter(None, 40) # 11 Hz stimulation, no need to keep higher freqs +bids_path = mne_bids.BIDSPath( + root=data_path, + subject="01", + task="phantom", + run="01", + datatype="meg", +) +# ignore warning about misc units +raw = mne_bids.read_raw_bids(bids_path).load_data() + +# Let's apply a little bit of preprocessing (temporal filtering and reference +# regression) +picks_artifact = ["MISC 001", "MISC 002", "MISC 003"] +picks = np.r_[ + mne.pick_types(raw.info, meg=True), + mne.pick_channels(raw.info["ch_names"], picks_artifact), +] +raw.filter(None, 40, picks=picks) +mne.preprocessing.regress_artifact( + raw, picks="meg", picks_artifact=picks_artifact, copy=False, proj=False +) plot_scalings = dict(mag=5e-12) # large-amplitude sinusoids raw_plot_kwargs = dict(duration=15, n_channels=50, scalings=plot_scalings) -raw.plot(**raw_plot_kwargs) +events, event_id = mne.events_from_annotations(raw) +raw.plot(events=events, **raw_plot_kwargs) +n_dip = len(event_id) +assert n_dip == 49 # sanity check # %% # We can also look at the power spectral density to see the phantom oscillations at @@ -46,81 +64,10 @@ fig.axes[0].axvline(dip_freq, color="r", ls="--", lw=2, zorder=4) # %% -# To find the events, we can look at the MISC channel that recorded the activations. -# Here we use a very simple thresholding approach to find the events. -# The MISC 017 channel holds the dipole activations, which are 2-cycle 11 Hz sinusoidal -# bursts with the initial sinusoidal deflection downward, so we do a little bit of -# signal manipulation to help :func:`~scipy.signal.find_peaks`. - -# Figure out events -dip_act, dip_t = raw["MISC 017"] -dip_act = dip_act[0] # 2D to 1D array -dip_act -= dip_act.mean() # remove DC offset -dip_act *= -1 # invert so first deflection is positive -thresh = np.percentile(dip_act, 90) -min_dist = raw.info["sfreq"] / dip_freq * 0.9 # 90% of period, to be safe -peaks = find_peaks(dip_act, height=thresh, distance=min_dist)[0] -assert len(peaks) % 2 == 0 # 2-cycle modulations -peaks = peaks[::2] # take only first peaks of each 2-cycle burst - -fig, ax = plt.subplots(layout="constrained", figsize=(12, 4)) -stop = int(15 * raw.info["sfreq"]) # 15 sec -ax.plot(dip_t[:stop], dip_act[:stop], color="k", lw=1) -ax.axhline(thresh, color="r", ls="--", lw=1) -peak_idx = peaks[peaks < stop] -ax.plot(dip_t[peak_idx], dip_act[peak_idx], "ro", zorder=5, ms=5) -ax.set(xlabel="Time (s)", ylabel="Dipole activation (AU)\n(MISC 017 adjusted)") -ax.set(xlim=dip_t[[0, stop - 1]]) - -# We know that there are 32 dipoles, so mark the first ones as well -n_dip = 49 -assert len(peaks) % n_dip == 0 # we found them all (hopefully) -ax.plot(dip_t[peak_idx[::n_dip]], dip_act[peak_idx[::n_dip]], "bo", zorder=4, ms=10) - -# Knowing we've caught the top of the first cycle of a 11 Hz sinusoid, plot onsets -# with red X's. -onsets = peaks - np.round(raw.info["sfreq"] / dip_freq / 4.0).astype( - int -) # shift to start -onset_idx = onsets[onsets < stop] -ax.plot(dip_t[onset_idx], dip_act[onset_idx], "rx", zorder=5, ms=5) - -# %% -# Given the onsets are now stored in ``peaks``, we can create our events array and plot -# on our raw data. +# Now we can figure out our epoching parameters and epoch the data and plot it. -n_rep = len(peaks) // n_dip -events = np.zeros((len(peaks), 3), int) -events[:, 0] = onsets + raw.first_samp -events[:, 2] = np.tile(np.arange(1, n_dip + 1), n_rep) -raw.plot(events=events, **raw_plot_kwargs) - -# %% -# Now we can figure out our epoching parameters and epoch the data, sanity checking -# some values along the way knowing how the stimulation was done. - -# Sanity check and determine epoching params -deltas = np.diff(events[:, 0], axis=0) -group_deltas = deltas[n_dip - 1 :: n_dip] / raw.info["sfreq"] # gap between 49 and 1 -assert (group_deltas > 0.8).all() -assert (group_deltas < 0.9).all() -others = np.delete(deltas, np.arange(n_dip - 1, len(deltas), n_dip)) # remove 49->1 -others = others / raw.info["sfreq"] -assert (others > 0.25).all() -assert (others < 0.3).all() -tmax = 1 / dip_freq * 2.0 # 2 cycles -tmin = tmax - others.min() -assert tmin < 0 -epochs = mne.Epochs( - raw, - events, - tmin=tmin, - tmax=tmax, - baseline=(None, 0), - decim=10, - picks="data", - preload=True, -) +tmin, tmax = -0.08, 0.18 +epochs = mne.Epochs(raw, tmin=tmin, tmax=tmax, decim=10, picks="data", preload=True) del raw epochs.plot(scalings=plot_scalings) @@ -131,7 +78,7 @@ t_peak = 1.0 / dip_freq / 4.0 data = np.zeros((len(epochs.ch_names), n_dip)) for di in range(n_dip): - data[:, [di]] = epochs[str(di + 1)].average().crop(t_peak, t_peak).data + data[:, [di]] = epochs[f"dip{di + 1:02d}"].average().crop(t_peak, t_peak).data evoked = mne.EvokedArray(data, epochs.info, tmin=0, comment="KIT phantom activations") evoked.plot_joint() @@ -141,22 +88,12 @@ trans = mne.transforms.Transform("head", "mri", np.eye(4)) sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08) cov = mne.compute_covariance(epochs, tmax=0, method="empirical") -# We need to correct the ``dev_head_t`` because it's incorrect for these data! -# relative to the helmet: hleft, forward, up -translation = mne.transforms.translation(x=0.01, y=-0.015, z=-0.088) -# pitch down (rot about x/R), roll left (rot about y/A), yaw left (rot about z/S) -rotation = mne.transforms.rotation( - x=np.deg2rad(5), - y=np.deg2rad(-1), - z=np.deg2rad(-3), -) -evoked.info["dev_head_t"]["trans"][:] = translation @ rotation dip, residual = mne.fit_dipole(evoked, cov, sphere, n_jobs=None) # %% # Finally let's look at the results. -# sphinx_gallery_thumbnail_number = 7 +# sphinx_gallery_thumbnail_number = 5 print(f"Average amplitude: {np.mean(dip.amplitude) * 1e9:0.1f} nAm") print(f"Average GOF: {np.mean(dip.gof):0.1f}%") diff --git a/tutorials/preprocessing/15_handling_bad_channels.py b/tutorials/preprocessing/15_handling_bad_channels.py index daac97976a5..06e9ffd6e53 100644 --- a/tutorials/preprocessing/15_handling_bad_channels.py +++ b/tutorials/preprocessing/15_handling_bad_channels.py @@ -238,8 +238,9 @@ fig.suptitle(title, size="xx-large", weight="bold") # %% -# Note that we used the ``exclude=[]`` trick in the call to -# :meth:`~mne.io.Raw.pick_types` to make sure the bad channels were not +# Note that the method :meth:`~mne.io.Raw.pick()` default +# arguments includes ``exclude=()`` which ensures that bad +# channels are not # automatically dropped from the selection. Here is the corresponding example # with the interpolated gradiometer channel; since there are more channels # we'll use a more transparent gray color this time: