Release/0.4.6

docs/source/overview/ovw_whatsnew.rst

@@ -3,23 +3,41 @@
 What's new
 ==========
 
-v0.4.5
+v0.4.6
 ------
 
 Bug fixes
 +++++++++
-* Compatibility with NumPy 2 (:commit:`0e85f55d`, :commit:`bc45cab7`)
-* Support mne-python versions newer than 1.5.0, up to the current release (:commit:`9a81a6dc`)
-* Restore a compatibility shim so mne keeps working with NumPy >= 2.4, which removed ``np.in1d`` while mne still relies on it internally (:commit:`ad7cfe65`)
-* Fix type comparison in :func:`frites.dataset.ds_utils.multi_to_uni_conditions` and in :func:`frites.plot.plot_conn_circle` (:commit:`cee7ed4a`)
-* Fix :func:`frites.conn.conn_spec` in multitaper mode : the cross- and auto-spectra are now averaged over tapers (the complex coefficients used to be averaged before, cancelling most of the signal). Coherence and PLV values change for every multitaper user (:pull:`69`)
+* Fix :func:`frites.conn.conn_spec` in multitaper mode : the cross- and auto-spectra are now averaged over tapers (the complex coefficients used to be averaged before, cancelling most of the signal). **Coherence and PLV values change for every multitaper user** (:pull:`69`)
 * Fix the hanning smoothing kernel of :func:`frites.conn.conn_spec` (``sm_times`` of 2 samples gave a NaN output, 3 samples no smoothing) (:pull:`69`)
 * The cross-spectrum (``metric='sxy'``) of :func:`frites.conn.conn_spec` is now complex (the imaginary part used to be silently dropped) (:pull:`69`)
-* :func:`frites.conn.conn_spec` sets ``zero_mean=False`` explicitly in the time-frequency decomposition : mne changed its default to ``True``, which silently changed the results between mne versions
+* :func:`frites.conn.conn_spec` sets ``zero_mean=False`` explicitly in the time-frequency decomposition : mne changed its default to ``True``, which silently changed the results between mne versions (:pull:`70`)
+* :func:`frites.set_mpl_style` no longer depends on ``pkg_resources`` (absent from Python >= 3.12 environments) (:pull:`70`)
+* Remove the ``np.in1d`` compatibility shim added in v0.4.5 : no supported mne release uses ``np.in1d`` and the shim modified the NumPy namespace globally
+* :func:`frites.stats.trial_swap_surrogates` ignored its ``random_state`` (the trials were shuffled with the global random generator) ; surrogates are now reproducible
+* New ``random_state`` parameter for :func:`frites.simulations.sim_local_ccd_ms` (the conditions were drawn from the global random generator)
+* The copula normalization (:func:`frites.core.copnorm_nd` and the GCMI estimators built on it) now raises a ``ValueError`` on NaN values instead of silently ranking them and returning a finite, meaningless mutual information
+* :class:`frites.workflow.WfMi` and :class:`frites.workflow.WfConnComod` raise an explicit ``ValueError`` when ``n_perm=0`` is combined with a permutation-based correction (used to fail with a bare ``AssertionError``) ; use ``mcp='noperm'`` to only compute the mutual information
+* :func:`frites.conn.conn_reshape_undirected` and :func:`frites.conn.conn_reshape_directed` no longer rely on deprecated xarray behaviours (``xr.concat`` default ``coords`` and implicit ``pandas.MultiIndex`` promotion) that were about to change their results
+* :func:`frites.conn.conn_spec` no longer parallelizes both the time-frequency decomposition and the loop over pairs (nested parallelism oversubscribed the CPU with the default ``n_jobs=-1``)
+* :func:`frites.plot.plot_conn_circle` uses ``Colormap.with_extremes`` (``set_bad`` is pending deprecation in matplotlib)
 
 Dependencies
 ++++++++++++
-* Tested against mne 1.5 up to 1.13 (mne >= 1.13 requires Python >= 3.11) and Python 3.10 to 3.12 ; the continuous integration matrix moved from Python 3.8/3.9 to 3.10/3.11/3.12
+* Python >= 3.10 is now required (``python_requires``) ; the continuous integration matrix moved from Python 3.8/3.9 to 3.10/3.11/3.12
+* Tested against mne 1.5 up to 1.13 (mne >= 1.13 requires Python >= 3.11) and NumPy 2.x up to 2.5
+* Documentation build : ``sphinxcontrib-bibtex >= 2.6`` and ``sphinx >= 7.4`` (:pull:`70`)
+* ``xarray >= 2023.8`` (``xarray.Coordinates.from_pandas_multiindex``)
+
+v0.4.5
+------
+
+Bug fixes
++++++++++
+* Compatibility with NumPy 2 (:commit:`0e85f55d`, :commit:`bc45cab7`)
+* Support mne-python versions newer than 1.5.0, up to the current release (:commit:`9a81a6dc`)
+* Restore a compatibility shim so mne keeps working with NumPy >= 2.4, which removed ``np.in1d`` while mne still relies on it internally (:commit:`ad7cfe65`)
+* Fix type comparison in :func:`frites.dataset.ds_utils.multi_to_uni_conditions` and in :func:`frites.plot.plot_conn_circle` (:commit:`cee7ed4a`)
 
 v0.4.4
 ------
@@ -30,7 +48,7 @@ New Features
 * New function :func:`frites.core.ent_nd_g` to compute entropy on tensors (:commit:`17587a15`)
 * New function :func:`frites.conn.conn_ii` to estimate the interaction information (:commit:`10938b46`)
 * New function :func:`frites.conn.conn_pid` to estimate the partial information decomposition (:commit:`ac9798dd`)
-- New function :func:`frites.conn.conn_fit` to estimate the feature specific information transfer (:pull:`59`) - :author:`aopy`
+* New function :func:`frites.conn.conn_fit` to estimate the feature specific information transfer (:pull:`59`) - :author:`aopy`
 
 Bug fixes
 +++++++++

frites/__init__.py

@@ -6,26 +6,11 @@
 """
 import logging
 
-import numpy as np
-
-# -----------------------------------------------------------------------------
-# numpy / mne transitional compatibility shim
-#
-# numpy >= 2.4 removed `np.in1d` (deprecated since numpy 2.0 in favor of
-# `np.isin`), but mne (up to and including 1.12.1, its latest release as of
-# writing) still calls `np.in1d` internally (e.g. in mne/epochs.py) when
-# building Epochs. mne's development branch has already dropped that call,
-# but no release with the fix is out yet. Restore the alias so mne keeps
-# working until such a release is available.
-# TODO: drop this once the minimum supported mne version no longer needs it.
-if not hasattr(np, 'in1d'):
-    np.in1d = np.isin
-
 from frites import (
     io, core, conn, plot, stats, utils, workflow, simulations,  estimator  # noqa
 )
 
-__version__ = "0.4.5"
+__version__ = "0.4.6"
 
 # -----------------------------------------------------------------------------
 # Set 'info' as the default logging level

frites/conn/conn_spec.py

@@ -325,7 +325,7 @@ def conn_spec(
         w = _tf_decomp(
             x[tr, ...], sfreq, freqs, n_cycles=n_cycles, decim=decim,
             mode=mode, mt_bandwidth=mt_bandwidth, kw_cwt=kw_cwt, kw_mt=kw_mt,
-            n_jobs=n_jobs)
+            n_jobs=1)  # parallelism is at the pair level (below)
 
         # ----------------------------- CONN TRIALS ---------------------------
         # give indication about computed trials

frites/conn/conn_utils.py

@@ -474,10 +474,13 @@ def conn_reshape_directed(
     # transpose, reindex and reorder (if needed)
     if 'direction' in list(da.dims):
         da_xy, da_yx = da.sel(direction='x->y'), da.sel(direction='y->x')
+        # (explicit `coords` / `compat` : xarray is changing its defaults)
         if net:
-            da = xr.concat((da_xy - da_yx, da_xy - da_yx), axis)
+            da = xr.concat((da_xy - da_yx, da_xy - da_yx), axis,
+                           coords='different', compat='equals')
         else:
-            da = xr.concat((da_xy, da_yx), axis)
+            da = xr.concat((da_xy, da_yx), axis, coords='different',
+                           compat='equals')
         s_, t_ = sources + targets, targets + sources
     else:
         s_, t_ = sources, targets
@@ -516,17 +519,19 @@ def _dataarray_unstack(
     dim_names = dim_names[:cut_at] + ['sources', 'targets'] + dim_names[
         cut_at + 1:]
 
-    # build the multi-index
-    da[axis] = pd.MultiIndex.from_arrays(
+    # build the multi-index (xarray no longer promotes a raw pandas
+    # MultiIndex assigned to a coordinate)
+    midx = pd.MultiIndex.from_arrays(
         [sources, targets], names=['sources', 'targets'])
+    da = da.assign_coords(xr.Coordinates.from_pandas_multiindex(midx, axis))
 
     # test for duplicated entries
     st_names = pd.Series([f"{s}-{t}" for s, t in zip(sources, targets)])
     duplicates = np.array(list(st_names.duplicated(keep='first')))
     if duplicates.any():
         logger.warning(f"Duplicated entries found and removed : "
                        f"{da[axis].data[duplicates]}")
-        da = da.sel(roi=~duplicates)
+        da = da.isel({axis: ~duplicates})
 
     # unstack to be 2D/3D
     da = da.unstack(fill_value=fill_value)

frites/conn/tests/test_conn_utils.py

@@ -1,4 +1,6 @@
 """Test window functions."""
+import warnings
+
 import numpy as np
 import xarray as xr
 
@@ -51,6 +53,11 @@ def test_conn_reshape_directed(self):
         # reshape it with the time dimension
         gc_times = conn_reshape_directed(gc.copy())
         assert gc_times.shape == (n_roi, n_roi, len(gc['times']))
+        # the reshape must not rely on deprecated xarray behaviours
+        with warnings.catch_warnings():
+            warnings.simplefilter('error', FutureWarning)
+            conn_reshape_directed(gc.copy())
+            conn_reshape_directed(gc.copy(), net=True)
         # try the reorder
         gc_order = conn_reshape_directed(gc.copy(), order=order)
         assert gc_order.shape == (2, 2, len(gc['times']))

frites/core/copnorm.py

@@ -19,6 +19,12 @@ def ctransform(x):
         Empirical CDF value along the last axis of x. Data is ranked and scaled
         within [0 1] (open interval)
     """
+    # NaN would silently be ranked like any other value
+    if np.isnan(x).any():
+        raise ValueError(
+            "The data contain NaN values. The copula normalization ranks NaN "
+            "like any other value, which silently corrupts the mutual "
+            "information estimate. Remove or impute NaN values first.")
     xr = np.argsort(np.argsort(x)).astype(float)
     xr += 1.
     xr /= float(xr.shape[-1] + 1)

frites/core/tests/test_copnorm.py

@@ -12,6 +12,16 @@ def test_copnorm_1d(self):
         arr = np.random.randint(0, 10, (20,))
         copnorm_1d(arr)
 
+    def test_copnorm_nan(self):
+        """NaN values must raise instead of being silently ranked."""
+        import pytest
+        x = np.random.rand(50)
+        x[3] = np.nan
+        with pytest.raises(ValueError, match="NaN"):
+            copnorm_1d(x)
+        with pytest.raises(ValueError, match="NaN"):
+            copnorm_nd(np.tile(x, (2, 1)), axis=-1)
+
     def test_copnorm_cat_1d(self):
         """Test function copnorm_cat_1d."""
         arr = np.random.randint(0, 10, (20,))

frites/plot/plot_conn.py

@@ -87,15 +87,15 @@ def _prepare_plot_conn(
     if (nodes_cmap is None) or isinstance(nodes_cmap, str):
         nodes_cmap = plt.get_cmap(nodes_cmap).copy()
     if nodes_bad:
-        nodes_cmap.set_bad(color=nodes_bad)
+        nodes_cmap = nodes_cmap.with_extremes(bad=nodes_bad)
     cfg['nodes_color'] = [nodes_cmap(k) for k in cfg['nodes_data']]
 
     # _________________________________ COLOR _________________________________
     # colormap
     if (cmap is None) or isinstance(cmap, str):
         cmap = plt.get_cmap(cmap).copy()
     if bad:
-        cmap.set_bad(color=bad)
+        cmap = cmap.with_extremes(bad=bad)
     cfg['cmap'] = cmap
 
     # vmin, vmax trick

frites/simulations/sim_local_mi.py

@@ -304,7 +304,7 @@ def sim_local_cd_ss(n_conditions=2, n_epochs=10, n_times=100, n_roi=1,
 """
 
 
-def sim_local_ccd_ms(n_subjects, **kwargs):
+def sim_local_ccd_ms(n_subjects, random_state=None, **kwargs):
     """Multi-subjects simulations for computing local MI (CCD).
 
     This function can be used for simulating local representations of mutual
@@ -315,6 +315,9 @@ def sim_local_ccd_ms(n_subjects, **kwargs):
     ----------
     n_subjects : int
         Number of subjects
+    random_state : int | None
+        Fix the random state of the machine (use it for reproducibility). If
+        None, a random state is randomly assigned.
     kwargs : dict | {}
         Additional arguments are send to the function :func:`sim_local_ccd_ss`
 
@@ -337,9 +340,13 @@ def sim_local_ccd_ms(n_subjects, **kwargs):
     n_c = kwargs.get('n_conditions', 2)
     if 'n_conditions' in list(kwargs.keys()):
         kwargs.pop('n_conditions')
-    x, y, roi, times = sim_local_cc_ms(n_subjects, **kwargs)
+    if not isinstance(random_state, int):
+        random_state = np.random.randint(1000)
+    rnd = np.random.RandomState(random_state)
+    x, y, roi, times = sim_local_cc_ms(
+        n_subjects, random_state=random_state, **kwargs)
     n_e = len(y[0])
-    z = [np.random.randint(0, n_c, (n_e,)) for k in range(n_subjects)]
+    z = [rnd.randint(0, n_c, (n_e,)) for k in range(n_subjects)]
 
     return x, y, z, roi, times
 

frites/simulations/tests/test_sim_local_mi.py

@@ -73,6 +73,19 @@ def test_sim_local_ccd_ms(self):
         assert times.shape == (n_times,)
         assert roi.shape == (n_roi,)
 
+    def test_sim_local_ccd_ms_random_state(self):
+        """Test the reproducibility of sim_local_ccd_ms."""
+        kw = dict(n_epochs=n_epochs, n_roi=n_roi, n_times=n_times,
+                  n_conditions=n_conditions, cl_index=cl_index)
+        x_0, y_0, z_0, _, _ = sim_local_ccd_ms(3, random_state=0, **kw)
+        x_0b, y_0b, z_0b, _, _ = sim_local_ccd_ms(3, random_state=0, **kw)
+        _, _, z_1, _, _ = sim_local_ccd_ms(3, random_state=1, **kw)
+        for k in range(3):
+            np.testing.assert_array_equal(x_0[k], x_0b[k])
+            np.testing.assert_array_equal(y_0[k], y_0b[k])
+            np.testing.assert_array_equal(z_0[k], z_0b[k])
+        assert not all([np.array_equal(z_0[k], z_1[k]) for k in range(3)])
+
     def test_sim_local_ccd_ss(self):
         """Test function sim_local_ccd_ss."""
         x, y, z, roi, times = sim_local_ccd_ms(

frites/stats/stats_nonparam.py

@@ -346,7 +346,7 @@ def trial_swap_surrogates(x, random_state=0, verbose=False):
 
     for c in range(n_nodes):
         # destroy roi-roi trial relation
-        np.random.shuffle(trials)
+        rnd.shuffle(trials)
         x_surr[:, c, :] = x[trials, c, :]
         pbar.update_with_increment_value(1)
 

frites/stats/tests/test_stats_nonparam.py

@@ -4,7 +4,7 @@
 
 from frites.stats import (permute_mi_vector, permute_mi_trials,
                           bootstrap_partitions, dist_to_ci,
-                          confidence_interval)
+                          confidence_interval, trial_swap_surrogates)
 
 rnd = np.random.RandomState(0)
 
@@ -64,6 +64,20 @@ def test_bootstrap_partitions(self):
             assert (2 <= k[2:4].min()) and (k[2:4].max() <= 3)
             assert (4 <= k[4::].min()) and (k[4::].max() <= 5)
 
+    def test_trial_swap_surrogates(self):
+        """Test function trial_swap_surrogates (reproducibility)."""
+        x = np.random.RandomState(0).rand(10, 3, 20)
+        s_0 = trial_swap_surrogates(x, random_state=0, verbose=False)
+        s_0b = trial_swap_surrogates(x, random_state=0, verbose=False)
+        s_1 = trial_swap_surrogates(x, random_state=1, verbose=False)
+        # same seed -> same surrogates ; different seed -> different ones
+        np.testing.assert_array_equal(s_0, s_0b)
+        assert not np.array_equal(s_0, s_1)
+        # every channel is a permutation of the trials of that channel
+        for c in range(x.shape[1]):
+            np.testing.assert_array_equal(
+                np.sort(s_0[:, c, :], axis=0), np.sort(x[:, c, :], axis=0))
+
     def test_dist_to_ci(self):
         """Test function dist_to_ci."""
         # standard definition

frites/workflow/tests/test_wf_mi.py

@@ -77,6 +77,18 @@ def test_mi_ccd(self):
             WfMi(mi_type='ccd', inference='rfx', estimator=estimator,
                  verbose=False).fit(dt, **kw_mi)
 
+    def test_n_perm_zero(self):
+        """n_perm=0 is only valid without a permutation-based correction."""
+        import pytest
+        y, gt = sim_mi_cc(x.copy(), snr=1.)
+        dt = DatasetEphy(x.copy(), y=y, roi=roi, times=time)
+        with pytest.raises(ValueError, match="n_perm=0"):
+            WfMi('cc', 'ffx', verbose=False).fit(
+                dt, mcp='cluster', n_perm=0, n_jobs=1)
+        mi, pv = WfMi('cc', 'ffx', verbose=False).fit(
+            dt, mcp='noperm', n_perm=0, n_jobs=1)
+        assert pv.min() == pv.max() == 1.
+
     def test_no_stat(self):
         """Test on no stats / no permutations / don't repeat computations."""
         y, gt = sim_mi_cc(x.copy(), snr=1.)

frites/workflow/wf_conn_comod.py

@@ -200,6 +200,12 @@ def fit(self, dataset, mcp='cluster', n_perm=1000, cluster_th=None,
         # don't compute permutations if mcp is either nostat / None
         if mcp in ['noperm', None]:
             n_perm = 0
+        elif (n_perm < 1) and (mcp != 'nostat'):
+            raise ValueError(
+                f"n_perm={n_perm} : permutations are required for the "
+                f"multiple comparison correction mcp='{mcp}'. Use "
+                "mcp='noperm' (or None) to only compute the mutual "
+                "information")
         # get important dataset's variables
         self._times = dataset.times
 

frites/workflow/wf_mi.py

@@ -226,6 +226,12 @@ def fit(self, dataset=None, mcp='cluster', n_perm=1000, cluster_th=None,
             # don't compute permutations if mcp is either nostat / None
             if mcp in ['noperm', None]:
                 n_perm = 0
+            elif (n_perm < 1) and (mcp != 'nostat'):
+                raise ValueError(
+                    f"n_perm={n_perm} : permutations are required for the "
+                    f"multiple comparison correction mcp='{mcp}'. Use "
+                    "mcp='noperm' (or None) to only compute the mutual "
+                    "information")
 
             # get needed dataset's informations
             self._times, self._roi = dataset.times, dataset.roi_names

frites/workflow/wf_stats.py

@@ -101,7 +101,8 @@ def fit(self, effect, perms, inference='rfx', mcp='cluster', tail=1,
         assert all([isinstance(k, np.ndarray) and k.ndim >= 2 for k in effect])
         n_roi, n_times, tvalues = len(effect), effect[0].shape[-1], None
         # don't compute statistics if `mcp` is None
-        if (mcp in [None, 'noperm']) or not len(perms):
+        no_perm = (not len(perms)) or all([k.size == 0 for k in perms])
+        if (mcp in [None, 'noperm']) or no_perm:
             return np.ones((n_times, n_roi), dtype=float), tvalues
         assert all([isinstance(k, np.ndarray) and k.ndim >= 3 for k in perms])
         assert len(effect) == len(perms)

requirements.txt

@@ -2,7 +2,7 @@ numpy
 scipy
 mne>=1.5.0
 joblib
-xarray
+xarray>=2023.8
 netCDF4
 h5netcdf
 scikit-learn

setup.py

@@ -2,9 +2,18 @@
 # -*- coding: utf-8 -*-
 # License: 3-clause BSD
 import os
+import re
 from setuptools import setup, find_packages
 
-__version__ = "0.4.5"
+
+def get_version():
+    """Read the version from frites/__init__.py (single source of truth)."""
+    init = os.path.join(os.path.dirname(__file__), 'frites', '__init__.py')
+    with open(init) as f:
+        return re.search(r'^__version__ = "([^"]+)"', f.read(), re.M).group(1)
+
+
+__version__ = get_version()
 NAME = 'frites'
 AUTHOR = "BraiNets"
 MAINTAINER = "Andrea Brovelli"
@@ -39,6 +48,7 @@ def read(fname):
 setup(
     name=NAME,
     version=__version__,
+    python_requires='>=3.10',
     packages=find_packages(),
     package_dir={'frites': 'frites'},
     package_data=PACKAGE_DATA,