Remove tensorpac dependency (#239)

Re-implement ndPAC

Author: raphaelvallat

docs/conf.py

@@ -585,7 +585,6 @@
     "scipy": ("https://docs.scipy.org/doc/scipy", None),
     "seaborn": ("https://seaborn.pydata.org", None),
     "sklearn": ("https://scikit-learn.org/stable", None),
-    "tensorpac": ("https://etiennecmb.github.io/tensorpac", None),
 }
 
 # The maximum number of days to cache remote inventories.

docs/index.rst

@@ -70,7 +70,6 @@ The core dependencies of YASA are:
 * `lspopt <https://github.com/hbldh/lspopt>`_ >= 1.4
 
 Some features require optional dependencies (`SleepECG <https://sleepecg.readthedocs.io/>`_,
-`TensorPAC <https://etiennecmb.github.io/tensorpac/>`_,
 `PyRiemann <https://pyriemann.readthedocs.io/>`_,
 `ipywidgets <https://ipywidgets.readthedocs.io/>`_).
 

pyproject.toml

@@ -45,10 +45,9 @@ dependencies = [
 
 [project.optional-dependencies]
 heart = ["sleepecg>=0.5.0"]
-pac   = ["tensorpac>=0.6.5"]
 art   = ["pyriemann>=0.2.7"]
 plot  = ["ipywidgets"]
-full  = ["yasa[heart,pac,art,plot]"]
+full  = ["yasa[heart,art,plot]"]
 
 [dependency-groups]
 dev = [

src/yasa/detection.py

@@ -21,11 +21,12 @@
 from scipy.stats import circmean
 from sklearn.ensemble import IsolationForest
 
-from .io import is_pyriemann_installed, is_tensorpac_installed, set_log_level
+from .io import is_pyriemann_installed, set_log_level
 from .numba import _detrend, _rms
 from .others import (
     _index_to_events,
     _merge_close,
+    _norm_direct_pac,
     _zerocrossings,
     get_centered_indices,
     moving_transform,
@@ -1576,9 +1577,8 @@ def sw_detect(
           be calculated for each slow-waves using a 2-seconds epoch centered around the negative
           peak of the slow-waves (i.e. 1 second on each side).
 
-        * ``p`` is a parameter passed to the :py:func:`tensorpac.methods.norm_direct_pac``
-          function. It represents the p-value to use for thresholding of unreliable coupling
-          values. Sub-threshold PAC values will be set to 0. To disable this behavior (no masking),
+        * ``p`` is the p-value used for thresholding of unreliable coupling values (ndPAC).
+          Sub-threshold PAC values will be set to 0. To disable this behavior (no masking),
           use ``p=1`` or ``p=None``.
 
         .. versionadded:: 0.6.0
@@ -1710,9 +1710,6 @@ def sw_detect(
 
     # Extract the spindles-related sigma signal for coupling
     if coupling:
-        is_tensorpac_installed()
-        import tensorpac.methods as tpm
-
         # The width of the transition band is set to 1.5 Hz on each side,
         # meaning that for freq_sp = (12, 15 Hz), the -6 dB points are located
         # at 11.25 and 15.75 Hz. The frequency band for the amplitude signal
@@ -1925,9 +1922,7 @@ def sw_detect(
             # 3) Normalized Direct PAC, with thresholding
             # Unreliable values are set to 0
             ndp = np.squeeze(
-                tpm.norm_direct_pac(
-                    sw_pha_ev[None, ...], sp_amp_ev[None, ...], p=coupling_params["p"]
-                )
+                _norm_direct_pac(sw_pha_ev[None, ...], sp_amp_ev[None, ...], p=coupling_params["p"])
             )
             sw_params["ndPAC"] = np.ones(n_peaks) * np.nan
             sw_params["ndPAC"][idx_valid] = ndp

src/yasa/io.py

@@ -35,14 +35,6 @@ def set_log_level(verbose=None):
             raise ValueError("verbose must be in %s" % ", ".join(LOGGING_TYPES))
 
 
-def is_tensorpac_installed():
-    """Test if tensorpac is installed."""
-    try:
-        import tensorpac  # noqa
-    except ImportError:  # pragma: no cover
-        raise ImportError("tensorpac needs to be installed. Please use `pip install yasa[pac]`.")
-
-
 def is_pyriemann_installed():
     """Test if pyRiemann is installed."""
     try:

src/yasa/others.py

@@ -7,6 +7,7 @@
 import numpy as np
 from numpy.lib.stride_tricks import as_strided
 from scipy.interpolate import interp1d
+from scipy.special import erfinv
 
 from .numba import _corr, _covar, _rms, _slope_lstsq
 
@@ -482,3 +483,36 @@ def rng(x):
     idx_ep_nomask = np.unique(idx_ep.nonzero()[0])
     idx_ep = np.ma.compress_rows(idx_ep)
     return idx_ep, idx_ep_nomask
+
+
+def _norm_direct_pac(pha, amp, p=0.05):
+    """Normalized direct PAC (ndPAC).
+
+    Re-implementation of tensorpac's ``norm_direct_pac`` (Ozkurt et al. 2012).
+
+    Parameters
+    ----------
+    pha : array_like
+        Phase array of shape (n_pha, ..., n_times).
+    amp : array_like
+        Amplitude array of shape (n_amp, ..., n_times).
+    p : float | .05
+        P-value threshold. Sub-threshold PAC values are set to 0.
+        Use ``p=1`` or ``p=None`` to disable thresholding.
+
+    Returns
+    -------
+    pac : array_like
+        Phase-amplitude coupling array of shape (n_amp, n_pha, ...).
+    """
+    n_times = amp.shape[-1]
+    amp = np.subtract(amp, np.mean(amp, axis=-1, keepdims=True))
+    amp = np.divide(amp, np.std(amp, ddof=1, axis=-1, keepdims=True))
+    pac = np.abs(np.einsum("i...j, k...j->ik...", amp, np.exp(1j * pha)))
+    if p == 1.0 or p is None:
+        return pac / n_times
+    s = pac**2
+    pac /= n_times
+    xlim = n_times * erfinv(1 - p) ** 2
+    pac[s <= 2 * xlim] = 0.0
+    return pac

tests/test_io.py

@@ -8,7 +8,6 @@
 from yasa.io import (
     is_pyriemann_installed,
     is_sleepecg_installed,
-    is_tensorpac_installed,
     set_log_level,
 )
 
@@ -38,6 +37,5 @@ def test_logger(self):
 
     def test_dependence(self):
         """Test dependancies."""
-        is_tensorpac_installed()
         is_pyriemann_installed()
         is_sleepecg_installed()