diff --git a/docs/source/overview/ovw_whatsnew.rst b/docs/source/overview/ovw_whatsnew.rst index 07e976990..a30cc9e5a 100644 --- a/docs/source/overview/ovw_whatsnew.rst +++ b/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 +++++++++ diff --git a/frites/__init__.py b/frites/__init__.py index 640eed497..ec0855eef 100644 --- a/frites/__init__.py +++ b/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 diff --git a/frites/conn/conn_spec.py b/frites/conn/conn_spec.py index 22dc6f3d5..4a67cd190 100644 --- a/frites/conn/conn_spec.py +++ b/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 diff --git a/frites/conn/conn_utils.py b/frites/conn/conn_utils.py index 89375434d..7150d91b4 100644 --- a/frites/conn/conn_utils.py +++ b/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,9 +519,11 @@ 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)]) @@ -526,7 +531,7 @@ def _dataarray_unstack( 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) diff --git a/frites/conn/tests/test_conn_utils.py b/frites/conn/tests/test_conn_utils.py index 399f1fe79..91d41fb94 100644 --- a/frites/conn/tests/test_conn_utils.py +++ b/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'])) diff --git a/frites/core/copnorm.py b/frites/core/copnorm.py index b6a701118..52ec64323 100644 --- a/frites/core/copnorm.py +++ b/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) diff --git a/frites/core/tests/test_copnorm.py b/frites/core/tests/test_copnorm.py index b2e364985..567ef911c 100644 --- a/frites/core/tests/test_copnorm.py +++ b/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,)) diff --git a/frites/plot/plot_conn.py b/frites/plot/plot_conn.py index c5632a9e9..6f8835e82 100644 --- a/frites/plot/plot_conn.py +++ b/frites/plot/plot_conn.py @@ -87,7 +87,7 @@ 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 _________________________________ @@ -95,7 +95,7 @@ def _prepare_plot_conn( 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 diff --git a/frites/simulations/sim_local_mi.py b/frites/simulations/sim_local_mi.py index 71c218b61..b885afecf 100644 --- a/frites/simulations/sim_local_mi.py +++ b/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 diff --git a/frites/simulations/tests/test_sim_local_mi.py b/frites/simulations/tests/test_sim_local_mi.py index a3ddc5ad4..9668223b3 100644 --- a/frites/simulations/tests/test_sim_local_mi.py +++ b/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( diff --git a/frites/stats/stats_nonparam.py b/frites/stats/stats_nonparam.py index f6f732a2c..44fe96395 100644 --- a/frites/stats/stats_nonparam.py +++ b/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) diff --git a/frites/stats/tests/test_stats_nonparam.py b/frites/stats/tests/test_stats_nonparam.py index aec99fad8..678c20451 100644 --- a/frites/stats/tests/test_stats_nonparam.py +++ b/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 diff --git a/frites/workflow/tests/test_wf_mi.py b/frites/workflow/tests/test_wf_mi.py index 1d5ef038e..636e03789 100644 --- a/frites/workflow/tests/test_wf_mi.py +++ b/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.) diff --git a/frites/workflow/wf_conn_comod.py b/frites/workflow/wf_conn_comod.py index 7dca19296..31ec414e4 100644 --- a/frites/workflow/wf_conn_comod.py +++ b/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 diff --git a/frites/workflow/wf_mi.py b/frites/workflow/wf_mi.py index 350146e57..3e6a73fab 100644 --- a/frites/workflow/wf_mi.py +++ b/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 diff --git a/frites/workflow/wf_stats.py b/frites/workflow/wf_stats.py index f5cd1831e..c85ed60f0 100644 --- a/frites/workflow/wf_stats.py +++ b/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) diff --git a/requirements.txt b/requirements.txt index 85e1bfc20..d3a0db434 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ numpy scipy mne>=1.5.0 joblib -xarray +xarray>=2023.8 netCDF4 h5netcdf scikit-learn diff --git a/setup.py b/setup.py index a28572513..a0a0a7430 100644 --- a/setup.py +++ b/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,