fix: clear errors instead of crashes in analyses, stats and utils

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

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)

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(

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,

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)