diff --git a/CHANGELOG.md b/CHANGELOG.md index 2aa7d1c..470188f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,10 +5,18 @@ ### Fixed - `hypyp.sync`: asking for the Metal backend on a metric that has no Metal kernel (PLV, CCorr, Coh, ImCoh, EnvCorr, PowCorr) used to run in NumPy silently, while the metric object reported `metal`. The computation still runs in NumPy, so no computed value changes, but it now warns and the metric object reports `numpy`. This holds for `optimization='metal'` and for a `priority` list that reaches Metal on a machine where Metal is available: the backends that follow in the list are not tried, which is left to 0.7.0. Only PLI, wPLI and ACCorr have a Metal kernel; the default selection of `optimization='auto'` is unchanged for the nine built-in metrics (#299) - `hypyp.sync`: when a backend of a `priority` list can neither run on the machine nor be computed by the metric, and no other GPU backend of the list can be used, the warning now says so (for example "'plv' has no Metal implementation") instead of "No GPU backend available" -- `hypyp.sync`: a metric whose `_backend` was set to a backend it does not implement used to compute in NumPy without notice. It still computes in NumPy, now with a warning. A `_backend` that is not the name of a backend raises a `ValueError` naming the metric and the backends it implements; called through `compute_sync`, that error is still reworded as an unsupported metric (#306) +- `hypyp.sync`: a metric whose `_backend` was set to a backend it does not implement used to compute in NumPy without notice. It still computes in NumPy, now with a warning. A `_backend` that is not the name of a backend raises a `ValueError` naming the metric and the backends it implements (#306) - Tests: the six Metal tests that compared NumPy with NumPy are replaced by tests of the fallback itself (warning, backend and result), the dispatch of every metric to every backend it implements is checked without a GPU, and the four tests that run a Metal kernel now assert that the Metal method was called (#300) +- `compute_sync` and `pair_connectivity` reported every `ValueError` as `Metric type "..." not supported`, so an unknown `optimization` was announced as an unsupported metric. Only an unknown metric name gets that message now; any other error keeps its own (#306) +- `pair_connectivity` crashed with an `IndexError` when the frequencies were given as a list and the data had a single epoch, a single channel per participant or a single frequency. Results for data without such a dimension are unchanged +- `utils.normalizing` raised an `UnboundLocalError` for a `type` name other than `'Zscore'` or `'Logratio'`; it now raises a `ValueError` naming the two valid types +- `utils.generate_virtual_epoch` failed inside the solver for an odd number of channels; it now raises a `ValueError` saying that the number must be even, and its docstring example uses four oscillators +- `utils.epochs_from_tasks` failed with `min() iterable argument is empty` when the onset event of a task was absent from the recording; the `ValueError` now names the task and the event +- `utils.create_epochs` warns when the two participants give different numbers of epochs, or, when both are given as lists, different numbers of recordings. It still returns what it returned before +- `stats.statscluster` raised an `UnboundLocalError` for an unknown `test`; it now raises a `ValueError` listing the valid names ### Added +- `stats.metaconn_matrix` has a `plot` argument. It defaults to `True`, the earlier behaviour; `plot=False` computes the matrices without drawing - `BaseMetric.supports(backend)` tells whether a metric implements a backend, for example `PLI.supports('metal')` ### Changed diff --git a/hypyp/analyses.py b/hypyp/analyses.py index 73d94ce..53f77dc 100644 --- a/hypyp/analyses.py +++ b/hypyp/analyses.py @@ -32,7 +32,7 @@ from mne.time_frequency import EpochsSpectrum from .mvarica import MVAR, connectivity_mvarica -from .sync import get_metric +from .sync import METRICS, get_metric def pow( @@ -457,10 +457,9 @@ def pair_connectivity( # compute instantaneous analytic signal from EEG data if type(frequencies) == list: - # average over tapers - values = np.mean( - compute_single_freq(data, sampling_rate, frequencies), 3 - ).squeeze() + # average over tapers; no squeeze, which would also drop a single + # epoch, channel or frequency and break the 5-D layout compute_sync expects + values = np.mean(compute_single_freq(data, sampling_rate, frequencies), 3) elif type(frequencies) == dict: values = compute_freq_bands(data, sampling_rate, frequencies) else: @@ -590,15 +589,14 @@ def compute_sync( } mode_normalized = mode_map.get(mode_lower, mode_lower) - # Get the metric from the sync module - try: - metric = get_metric( - mode_normalized, optimization=optimization, priority=priority - ) - con = metric.compute(complex_signal, n_samp, transpose_axes) - except ValueError: + # Only an unknown metric is reported as unsupported; any other ValueError + # (an unknown backend, for instance) keeps its own message + if mode_normalized not in METRICS: raise ValueError(f'Metric type "{mode}" not supported.') + metric = get_metric(mode_normalized, optimization=optimization, priority=priority) + con = metric.compute(complex_signal, n_samp, transpose_axes) + con = con.swapaxes(0, 1) # n_freq x n_epoch x 2*n_ch x 2*n_ch if epochs_average: con = np.nanmean(con, axis=1) diff --git a/hypyp/stats.py b/hypyp/stats.py index 293665e..d1ca111 100644 --- a/hypyp/stats.py +++ b/hypyp/stats.py @@ -324,6 +324,7 @@ def metaconn_matrix( electrodes: List[Tuple[int, int]], ch_con: scipy.sparse.csr_matrix, freqs_mean: List[float], + plot: bool = True, ) -> namedtuple: """ Compute a priori connectivity between pairs of sensors within one brain. @@ -343,6 +344,10 @@ def metaconn_matrix( freqs_mean : List[float] List of frequencies in the frequency-band-of-interest + plot : bool, optional + Whether to draw the meta-connectivity matrix in the current figure + (default=True, the historical behaviour) + Returns ------- metaconn_matrixTuple : namedtuple @@ -401,9 +406,9 @@ def metaconn_matrix( metaconn_mult = np.tile(metaconn, (l_freq, l_freq)) metaconn_freq = np.multiply(init, metaconn_mult) - # TODO: option with verbose - # vizualising the array - plt.spy(metaconn_freq) + if plot: + # vizualising the array + plt.spy(metaconn_freq) metaconn_matrixTuple = namedtuple("metaconn_matrix", ["metaconn", "metaconn_freq"]) @@ -678,6 +683,11 @@ def stat_fun(*arg): effects="all", pvalue=alpha, ) + else: + raise ValueError( + f"Unknown test '{test}'. Use 'ind ttest', 'rel ttest', 'f oneway' " + "or 'f multipleway'." + ) # computing the cluster permutation t test Stat_obs, clusters, cluster_p_values, H0 = permutation_cluster_test( diff --git a/hypyp/utils.py b/hypyp/utils.py index d742176..de54e33 100644 --- a/hypyp/utils.py +++ b/hypyp/utils.py @@ -15,6 +15,7 @@ import math import random import string +import warnings import numpy as np import pandas as pd @@ -26,8 +27,8 @@ def create_epochs( - raw_S1: mne.io.Raw, raw_S2: mne.io.Raw, duration: float -) -> Tuple[mne.Epochs, mne.Epochs]: + raw_S1: List[mne.io.Raw], raw_S2: List[mne.io.Raw], duration: float +) -> Tuple[List[mne.Epochs], List[mne.Epochs]]: """ Create epochs from continuous raw EEG data for two participants. @@ -37,11 +38,11 @@ def create_epochs( Parameters ---------- - raw_S1 : mne.io.Raw - Raw EEG data for participant 1 + raw_S1 : List[mne.io.Raw] + Raw EEG recordings for participant 1, one per run - raw_S2 : mne.io.Raw - Raw EEG data for participant 2 + raw_S2 : List[mne.io.Raw] + Raw EEG recordings for participant 2, in the same order duration : float Duration of each epoch in seconds @@ -73,6 +74,18 @@ def create_epochs( epoch_S1: List[mne.Epochs] = [] epoch_S2: List[mne.Epochs] = [] + # zip stops at the shorter list, so the extra recordings are not epoched + if ( + hasattr(raw_S1, "__len__") + and hasattr(raw_S2, "__len__") + and len(raw_S1) != len(raw_S2) + ): + warnings.warn( + "The two participants have different numbers of recordings " + f"({len(raw_S1)} and {len(raw_S2)}). Only the first " + f"{min(len(raw_S1), len(raw_S2))} of each are epoched." + ) + for raw1, raw2 in zip(raw_S1, raw_S2): # creating fixed events fixed_events1 = mne.make_fixed_length_events( @@ -130,6 +143,14 @@ def create_epochs( epoch2, reset_bads=True, mode="accurate", origin="auto", verbose=None ) + if len(epoch1) != len(epoch2): + warnings.warn( + "The two participants have different numbers of epochs " + f"({len(epoch1)} and {len(epoch2)}). Hyperscanning analyses need " + "epochs that correspond in time between the two participants; " + "check that the two recordings cover the same period." + ) + epoch_S1.append(epoch1) epoch_S2.append(epoch2) @@ -209,8 +230,6 @@ def merge(epoch_S1: mne.Epochs, epoch_S2: mne.Epochs) -> mne.Epochs: epoch_S1.info["highpass"] != epoch_S2.info["highpass"] or epoch_S1.info["lowpass"] != epoch_S2.info["lowpass"] ): - import warnings - warnings.warn("Filter settings differ between participants. Using S1 settings.") # checking wether data have the same size @@ -463,9 +482,13 @@ def normalizing(baseline: np.ndarray, task: np.ndarray, type: str) -> np.ndarray if type == "Zscore": s = np.subtract(m_task, m_baseline) Normed_task = np.divide(s, std_baseline) - if type == "Logratio": + elif type == "Logratio": d = np.divide(m_task, m_baseline) Normed_task = np.log10(d) + else: + raise ValueError( + f"Unknown normalization type '{type}'. Use 'Zscore' or 'Logratio'." + ) return Normed_task @@ -550,7 +573,8 @@ def generate_virtual_epoch( Template Epochs object to copy structure from W : np.ndarray - Coupling matrix between oscillators, with shape (n_channels, n_channels) + Coupling matrix between oscillators, with shape (n_channels, n_channels). + The number of channels must be even, otherwise a ValueError is raised. frequency_mean : float, optional Mean frequency of oscillators in Hz (default=10) @@ -589,11 +613,12 @@ def generate_virtual_epoch( Examples -------- - >>> # Create a simple coupling matrix (3 oscillators) + >>> # Create a simple coupling matrix (4 oscillators, as many as channels) >>> W = np.array([ - ... [0, 0.2, 0], - ... [0.2, 0, 0.2], - ... [0, 0.2, 0] + ... [0, 0.2, 0, 0], + ... [0.2, 0, 0.2, 0], + ... [0, 0.2, 0, 0.2], + ... [0, 0, 0.2, 0] ... ]) >>> # Generate simulated epochs in the alpha band >>> sim_epochs = generate_virtual_epoch( @@ -608,6 +633,12 @@ def generate_virtual_epoch( n_epo, n_chan, n_samp = epoch.get_data().shape sfreq = epoch.info["sfreq"] + # the initial phases are drawn for two halves of equal size + if n_chan % 2 != 0: + raise ValueError( + f"generate_virtual_epoch needs an even number of channels, got {n_chan}." + ) + Nt = n_samp * n_epo tmax = n_samp / sfreq * n_epo # s tv = np.linspace(0.0, tmax, Nt) @@ -786,6 +817,12 @@ def epochs_from_tasks( t_durations.append(t_duration) task_events.append([t_start, 0, onset_event_id]) + if len(task_events) == 0: + raise ValueError( + f'Cannot find onset of task "{task_key}": no event with ' + f'event_id "{onset_event_id}" in the annotations of the recording' + ) + all_epochs.append( mne.Epochs( raw, diff --git a/tests/test_stats.py b/tests/test_stats.py index 5cb567b..5fe8347 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -351,3 +351,112 @@ def test_pair_connectivity_accorr(epochs): # Diagonal should have high values (self-correlation) diag_values = np.array([con[0, i, i] for i in range(con.shape[1])]) assert np.mean(diag_values) > 0.5 + + +def _synthetic_pair(n_epochs, n_channels, n_times=512, seed=0): + """Random data shaped (2, n_epochs, n_channels, n_times), no download needed.""" + rng = np.random.default_rng(seed) + return rng.standard_normal((2, n_epochs, n_channels, n_times)) + + +def _plv_from_definition(data, sampling_rate, frequencies): + """PLV written from its definition, without compute_sync or hypyp.sync. + + Returns an array shaped (n_freq, n_epochs, 2*n_channels, 2*n_channels), + with the channels of participant 1 first. + """ + # (2, n_epochs, n_channels, n_freq, n_times), tapers averaged + signal = np.mean(analyses.compute_single_freq(data, sampling_rate, frequencies), 3) + phase = signal / np.abs(signal) + both = np.concatenate([phase[0], phase[1]], axis=1) + n_times = both.shape[-1] + plv = np.abs(np.einsum("ecft,edft->fecd", both, both.conj())) / n_times + return plv + + +@pytest.mark.parametrize( + "n_epochs, n_channels, frequencies, n_freq", + [ + (3, 4, [8, 12], 4), # ordinary case, already worked + (1, 4, [8, 12], 4), # single epoch + (3, 1, [8, 12], 4), # single channel per participant + (1, 1, [8, 12], 4), + (3, 4, [10, 11], 1), # single frequency + ], +) +def test_pair_connectivity_frequency_list_singleton_dims( + n_epochs, n_channels, frequencies, n_freq +): + """A frequency list must not crash when a dimension has length one, and + must give the PLV of the definition with the documented layout.""" + data = _synthetic_pair(n_epochs, n_channels) + expected = _plv_from_definition(data, 256, frequencies) + assert expected.shape == (n_freq, n_epochs, 2 * n_channels, 2 * n_channels) + + con = analyses.pair_connectivity( + data, + sampling_rate=256, + frequencies=frequencies, + mode="plv", + epochs_average=False, + ) + assert con.shape == expected.shape + np.testing.assert_allclose(con, expected, rtol=0, atol=1e-12) + + con = analyses.pair_connectivity( + data, sampling_rate=256, frequencies=frequencies, mode="plv" + ) + assert con.shape == (n_freq, 2 * n_channels, 2 * n_channels) + np.testing.assert_allclose(con, expected.mean(axis=1), rtol=0, atol=1e-12) + + +def test_compute_sync_reports_the_real_error(): + """Only an unknown metric is reported as an unsupported metric.""" + complex_signal = analyses.compute_freq_bands( + _synthetic_pair(3, 4), 256, {"alpha": [8, 12]} + ) + + with pytest.raises(ValueError, match='Metric type "nope" not supported.'): + analyses.compute_sync(complex_signal, "nope") + + with pytest.raises(ValueError) as excinfo: + analyses.compute_sync(complex_signal, "plv", optimization="bogus") + assert "not supported" not in str(excinfo.value) + assert "bogus" in str(excinfo.value) + + +def _cluster_data(seed=0): + """Two groups of 12 observations over 6 features; group 2 is larger on + the first three features.""" + rng = np.random.default_rng(seed) + low = rng.standard_normal((12, 6)) + high = rng.standard_normal((12, 6)) + high[:, :3] += 5.0 + adjacency = scipy.sparse.csr_matrix(np.eye(6) + np.eye(6, k=1) + np.eye(6, k=-1)) + return low, high, adjacency + + +def test_statscluster_unknown_test_name(): + low, high, adjacency = _cluster_data() + with pytest.raises(ValueError, match="bogus"): + stats.statscluster([low, high], "bogus", None, adjacency, 0, 50) + + +def test_metaconn_matrix_plot_argument(): + import matplotlib.pyplot as plt + + ch_con = scipy.sparse.csr_matrix(np.eye(4) + np.eye(4, k=1) + np.eye(4, k=-1)) + electrodes = [(0, 1), (1, 2), (2, 3)] + + plt.close("all") + quiet = stats.metaconn_matrix(electrodes, ch_con, [10, 11], plot=False) + assert plt.get_fignums() == [] + + # the default still draws, as before + drawn = stats.metaconn_matrix(electrodes, ch_con, [10, 11]) + assert len(plt.get_fignums()) == 1 + plt.close("all") + + np.testing.assert_array_equal(quiet.metaconn, drawn.metaconn) + np.testing.assert_array_equal(quiet.metaconn_freq, drawn.metaconn_freq) + assert quiet.metaconn_freq.shape == (6, 6) diff --git a/tests/test_utils.py b/tests/test_utils.py index 9b40b7a..45621bd 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -330,3 +330,91 @@ def test_key_dict_as_data_frame(): assert df.columns[2] == "c" assert df["a"][0] == "foo" assert df["b"][0] == "1" + + +def test_normalizing_unknown_type(): + # per channel and frequency: baseline mean 2 and standard deviation 1, + # task mean 5 + baseline = np.stack([np.full((2, 2), 1.0), np.full((2, 2), 3.0)]) + task = np.stack([np.full((2, 2), 4.0), np.full((2, 2), 6.0)]) + + with pytest.raises(ValueError, match="zscore"): + utils.normalizing(baseline, task, "zscore") + + # the two documented types still work + np.testing.assert_allclose( + utils.normalizing(baseline, task, "Zscore"), np.full((2, 2), 3.0) + ) + np.testing.assert_allclose( + utils.normalizing(baseline, task, "Logratio"), np.full((2, 2), np.log10(2.5)) + ) + + +def _fake_epochs(n_channels, n_epochs=2, n_times=50, sfreq=100): + info = mne.create_info( + [f"Foo {i}" for i in range(n_channels)], sfreq=sfreq, ch_types="eeg" + ) + data = np.random.randn(n_epochs, n_channels, n_times) + return mne.EpochsArray(data, info) + + +def test_generate_virtual_epoch_odd_channel_count(): + epochs = _fake_epochs(n_channels=3) + with pytest.raises(ValueError, match="even number of channels"): + utils.generate_virtual_epoch(epochs, W=np.zeros((3, 3))) + + +def test_generate_virtual_epoch_even_channel_count(): + epochs = _fake_epochs(n_channels=4) + simulated = utils.generate_virtual_epoch(epochs, W=np.zeros((4, 4))) + assert simulated.get_data(copy=True).shape == epochs.get_data(copy=True).shape + + +def test_epochs_from_tasks_missing_onset_event(): + raw = get_fake_raw() + raw.set_annotations( + mne.Annotations(onset=[2, 5], duration=[0, 0], description=[1, 2]) + ) + tasks = [utils.Task("ghost", onset_event_id=7, offset_event_id=2)] + + with pytest.raises(ValueError, match='"ghost".*7'): + utils.epochs_from_tasks(raw, tasks) + + +def _fake_raw_of_duration(duration, sfreq=100, n_channels=5): + info = mne.create_info( + [f"Foo {i}" for i in range(n_channels)], sfreq=sfreq, ch_types="eeg" + ) + return mne.io.RawArray(np.random.randn(n_channels, int(sfreq * duration)), info) + + +def test_create_epochs_warns_on_different_epoch_counts(): + raw_long = _fake_raw_of_duration(10) + raw_short = _fake_raw_of_duration(8) + + with pytest.warns(UserWarning, match="different numbers of epochs"): + epo1, epo2 = utils.create_epochs([raw_long], [raw_short], duration=1.0) + # a warning, not an error: both lists are still returned as before + assert len(epo1[0]) != len(epo2[0]) + + +def test_create_epochs_silent_on_equal_epoch_counts(): + import warnings + + raw1 = _fake_raw_of_duration(10) + raw2 = _fake_raw_of_duration(10) + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*different numbers of.*") + epo1, epo2 = utils.create_epochs([raw1], [raw2], duration=1.0) + # iterables without a length were accepted before and still are + utils.create_epochs(iter([raw1]), iter([raw2]), duration=1.0) + assert len(epo1[0]) == len(epo2[0]) + + +def test_create_epochs_warns_on_different_recording_counts(): + raws = [_fake_raw_of_duration(10) for _ in range(3)] + + with pytest.warns(UserWarning, match="different numbers of recordings"): + epo1, epo2 = utils.create_epochs(raws[:2], raws[2:], duration=1.0) + # as before, only the recordings that have a counterpart are epoched + assert len(epo1) == len(epo2) == 1