Merge pull request #776 from neuropsychology/dev

0.2.4

AUTHORS.rst

@@ -30,6 +30,7 @@ Core contributors
 Contributors
 -------------
 
+* `Gansheng Tan <https://github.com/GanshengT>`_ *(Washington University, USA)*
 * `Hung Pham <https://github.com/hungpham2511>`_ *(Eureka Robotics, Singapore)*
 * `Christopher Schölzel <https://github.com/CSchoel>`_ *(THM University of Applied Sciences, Germany)*
 * `Duy Le <https://github.com/duylp>`_ *(Hubble, Singapore)*
@@ -53,7 +54,7 @@ Contributors
 * `Marek Sokol <https://github.com/sokolmarek>`_ *(Faculty of Biomedical Engineering of the CTU in Prague, Czech Republic)*
 
 
-Thanks also to `Gansheng Tan <https://github.com/GanshengT>`_, `Chuan-Peng Hu <https://github.com/hcp4715>`_, `@ucohen <https://github.com/ucohen>`_, `Anthony Gatti <https://github.com/gattia>`_, `Julien Lamour <https://github.com/lamourj>`_, `@renatosc <https://github.com/renatosc>`_, `Nicolas Beaudoin-Gagnon <https://github.com/Fegalf>`_ and `@rubinovitz <https://github.com/rubinovitz>`_ for their contribution in `NeuroKit 1 <https://github.com/neuropsychology/NeuroKit.py>`_.
+Thanks also to `Chuan-Peng Hu <https://github.com/hcp4715>`_, `@ucohen <https://github.com/ucohen>`_, `Anthony Gatti <https://github.com/gattia>`_, `Julien Lamour <https://github.com/lamourj>`_, `@renatosc <https://github.com/renatosc>`_, `Nicolas Beaudoin-Gagnon <https://github.com/Fegalf>`_ and `@rubinovitz <https://github.com/rubinovitz>`_ for their contribution in `NeuroKit 1 <https://github.com/neuropsychology/NeuroKit.py>`_.
 
 
 

NEWS.rst

@@ -1,7 +1,16 @@
 News
 =====
 
+0.2.4
+-------------------
+Fixes
++++++++++++++
 
+* `eda_sympathetic()` has been reviewed: low-pass filter and resampling have been added to be in
+  line with the original paper
+* `eda_findpeaks()` using methods proposed in nabian2018 is reviewed and improved. Differentiation
+  has been added before smoothing. Skin conductance response criteria have been revised based on
+  the original paper.
 
 
 

docs/functions/eda.rst

@@ -31,14 +31,6 @@ Preprocessing
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_phasic
 
-*eda_autocor()*
-"""""""""""""""
-.. autofunction:: neurokit2.eda.eda_autocor
-
-*eda_changepoints()*
-"""""""""""""""""""""
-.. autofunction:: neurokit2.eda.eda_changepoints
-
 *eda_peaks()*
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_peaks
@@ -47,9 +39,6 @@ Preprocessing
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_fixpeaks
 
-*eda_sympathetic()*
-"""""""""""""""""""""
-.. autofunction:: neurokit2.eda.eda_sympathetic
 
 
 Analysis
@@ -62,6 +51,18 @@ Analysis
 """""""""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_intervalrelated
 
+*eda_sympathetic()*
+"""""""""""""""""""""
+.. autofunction:: neurokit2.eda.eda_sympathetic
+
+*eda_autocor()*
+"""""""""""""""
+.. autofunction:: neurokit2.eda.eda_autocor
+
+*eda_changepoints()*
+"""""""""""""""""""""
+.. autofunction:: neurokit2.eda.eda_changepoints
+
 
 
 Miscellaneous

docs/functions/hrv.rst

@@ -46,5 +46,5 @@ Intervals
 
 .. automodule:: neurokit2.hrv
    :members:
-   :exclude-members: hrv, hrv_time, hrv_frequency, hrv_nonlinear, hrv_rqa, hrv_rsa
+   :exclude-members: hrv, hrv_time, hrv_frequency, hrv_nonlinear, hrv_rqa, hrv_rsa, intervals_process, intervals_to_peaks
 

docs/installation.rst

@@ -23,7 +23,7 @@ Then, at the top of each of your Python script, you should be able to import the
 
     .. code-block:: console
 
-        pip install https://github.com/neuropsychology/neurokit/zipball/dev
+        pip install https://github.com/neuropsychology/neurokit/zipball/dev --upgrade
 
 
 

docs/readme/README_examples.py

@@ -10,6 +10,7 @@
 # Setup matplotlib with Agg to run on server
 matplotlib.use("Agg")
 plt.rcParams["figure.figsize"] = (10, 6.5)
+plt.rcParams["savefig.facecolor"] = "white"
 
 # =============================================================================
 # Quick Example

neurokit2/__init__.py

@@ -32,7 +32,7 @@
 from .video import *
 
 # Info
-__version__ = "0.2.3"
+__version__ = "0.2.4"
 
 
 # Maintainer info

neurokit2/complexity/__init__.py

@@ -8,6 +8,7 @@
 from .complexity_lyapunov import complexity_lyapunov
 from .complexity_relativeroughness import complexity_relativeroughness
 from .complexity_rqa import complexity_rqa
+from .entropy_angular import entropy_angular
 from .entropy_approximate import entropy_approximate
 from .entropy_attention import entropy_attention
 from .entropy_bubble import entropy_bubble
@@ -155,6 +156,7 @@
     "complexity_dfa",
     "complexity_relativeroughness",
     "complexity_rqa",
+    "entropy_angular",
     "entropy_maximum",
     "entropy_shannon",
     "entropy_shannon_joint",

neurokit2/complexity/entropy_angular.py

@@ -0,0 +1,147 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from .utils_complexity_embedding import complexity_embedding
+
+
+def entropy_angular(signal, delay=1, dimension=2, show=False, **kwargs):
+    """**Angular entropy (AngEn)**
+
+    The Angular Entropy (AngEn) is the name that we use in NeuroKit to refer to the complexity
+    method described in Nardelli et al. (2022), referred as comEDA due to its application to EDA
+    signal. The method comprises the following steps: 1) Phase space reconstruction, 2) Calculation
+    of the angular distances between all the pairs of points in the phase space; 3) Computation of
+    the probability density function (PDF) of the distances; 4) Quadratic Rényi entropy of the PDF.
+
+    Parameters
+    ----------
+    signal : Union[list, np.array, pd.Series]
+        The signal (i.e., a time series) in the form of a vector of values.
+    delay : int
+        Time delay (often denoted *Tau* :math:`\\tau`, sometimes referred to as *lag*) in samples.
+        See :func:`complexity_delay` to estimate the optimal value for this parameter.
+    dimension : int
+        Embedding Dimension (*m*, sometimes referred to as *d* or *order*). See
+        :func:`complexity_dimension` to estimate the optimal value for this parameter.
+    **kwargs : optional
+        Other arguments.
+
+    Returns
+    --------
+    angen : float
+        The Angular Entropy (AngEn) of the signal.
+    info : dict
+        A dictionary containing additional information regarding the parameters used
+        to compute the index.
+
+    See Also
+    --------
+    entropy_renyi
+
+    Examples
+    ----------
+    .. ipython:: python
+
+      import neurokit2 as nk
+
+      # Simulate a Signal with Laplace Noise
+      signal = nk.signal_simulate(duration=2, frequency=[5, 3], noise=0.1)
+
+      # Compute Angular Entropy
+      @savefig p_entropy_angular1.png scale=100%
+      angen, info = nk.entropy_angular(signal, delay=1, dimension=3, show=True)
+      @suppress
+      plt.close()
+
+
+    References
+    -----------
+    * Nardelli, M., Greco, A., Sebastiani, L., & Scilingo, E. P. (2022). ComEDA: A new tool for
+      stress assessment based on electrodermal activity. Computers in Biology and Medicine, 150,
+      106144.
+
+    """
+    # Sanity checks
+    if isinstance(signal, (np.ndarray, pd.DataFrame)) and signal.ndim > 1:
+        raise ValueError(
+            "Multidimensional inputs (e.g., matrices or multichannel data) are not supported yet."
+        )
+
+    # 1. Phase space reconstruction (time-delay embeddings)
+    embedded = complexity_embedding(signal, delay=delay, dimension=dimension)
+
+    # 2. Angular distances between all the pairs of points in the phase space
+    angles = _angular_distance(embedded)
+
+    # 3. Compute the probability density function (PDF) of the upper triangular matrix
+    bins, pdf = _kde_sturges(angles)
+
+    # 4. Apply the quadratic Rényi entropy to the PDF
+    angen = -np.log2(np.sum(pdf**2))
+
+    # Normalize to the range [0, 1] by the log of the number of bins
+
+    # Note that in the paper (eq. 4 page 4) there is a minus sign, but adding it would give
+    # negative values, plus the linked code does not seem to do that
+    # https://github.com/NardelliM/ComEDA/blob/main/comEDA.m#L103
+    angen = angen / np.log2(len(bins))
+
+    if show is True:
+        # Plot the PDF as a bar chart
+        plt.bar(bins[:-1], pdf, width=bins[1] - bins[0], align="edge", alpha=0.5)
+        # Set the x-axis limits to the range of the data
+        plt.xlim([np.min(angles), np.max(angles)])
+        # Print titles
+        plt.suptitle(f"Angular Entropy (AngEn) = {angen:.3f}")
+        plt.title("Distribution of Angular Distances:")
+
+    return angen, {"bins": bins, "pdf": pdf}
+
+
+def _angular_distance(m):
+    """
+    Compute angular distances between all the pairs of points.
+    """
+    # Get index of upper triangular to avoid double counting
+    idx = np.triu_indices(m.shape[0], k=1)
+
+    # compute the magnitude of each vector
+    magnitudes = np.linalg.norm(m, axis=1)
+
+    # compute the dot product between all pairs of vectors using np.matmul function, which is
+    # more efficient than np.dot for large matrices; and divide the dot product matrix by the
+    # product of the magnitudes to get the cosine of the angle
+    cos_angles = np.matmul(m, m.T)[idx] / np.outer(magnitudes, magnitudes)[idx]
+
+    # clip the cosine values to the range [-1, 1] to avoid any numerical errors and compute angles
+    return np.arccos(np.clip(cos_angles, -1, 1))
+
+
+def _kde_sturges(x):
+    """
+    Computes the PDF of a vector x using a kernel density estimator based on linear diffusion
+    processes with a Gaussian kernel. The number of bins of the PDF is chosen applying the Sturges
+    method.
+    """
+    # Estimate the bandwidth
+    iqr = np.percentile(x, 75) - np.percentile(x, 25)
+    bandwidth = 0.9 * iqr / (len(x) ** 0.2)
+
+    # Compute the number of bins using the Sturges method
+    nbins = int(np.ceil(np.log2(len(x)) + 1))
+
+    # Compute the bin edges
+    bins = np.linspace(np.min(x), np.max(x), nbins + 1)
+
+    # Compute the kernel density estimate
+    xi = (bins[:-1] + bins[1:]) / 2
+    pdf = np.sum(
+        scipy.stats.norm.pdf((xi.reshape(-1, 1) - x.reshape(1, -1)) / bandwidth), axis=1
+    ) / (len(x) * bandwidth)
+
+    # Normalize the PDF
+    pdf = pdf / np.sum(pdf)
+
+    return bins, pdf

neurokit2/complexity/entropy_differential.py

@@ -17,18 +17,17 @@ def entropy_differential(signal, base=2, **kwargs):
     ----------
     signal : Union[list, np.array, pd.Series]
         The signal (i.e., a time series) in the form of a vector of values.
-
+    base: float
+        The logarithmic base to use, defaults to ``2``, giving a unit in *bits*. Note that ``scipy.
+        stats.entropy()`` uses Euler's number (``np.e``) as default (the natural logarithm), giving
+        a measure of information expressed in *nats*.
     **kwargs : optional
         Other arguments passed to ``scipy.stats.differential_entropy()``.
 
     Returns
     --------
     diffen : float
         The Differential entropy of the signal.
-    base: float
-        The logarithmic base to use, defaults to ``2``, giving a unit in *bits*. Note that ``scipy.
-        stats.entropy()`` uses Euler's number (``np.e``) as default (the natural logarithm), giving
-        a measure of information expressed in *nats*.
     info : dict
         A dictionary containing additional information regarding the parameters used
         to compute Differential entropy.

neurokit2/complexity/entropy_sample.py

@@ -58,7 +58,7 @@ def entropy_sample(signal, delay=1, dimension=2, tolerance="sd", **kwargs):
 
       signal = nk.signal_simulate(duration=2, frequency=5)
 
-      sampen, parameters = nk.entropy_sample(signal)
+      sampen, parameters = nk.entropy_sample(signal, delay=1, dimension=2)
       sampen
 
     """

neurokit2/ecg/ecg_delineate.py

@@ -164,6 +164,11 @@ def ecg_delineate(
             "NeuroKit error: ecg_delineate(): 'method' should be one of 'peak'," "'cwt' or 'dwt'."
         )
 
+    # Ensure that all indices are not larger than ECG signal indices
+    for _, value in waves.items():
+        if value[-1] >= len(ecg_cleaned):
+            value[-1] = np.nan
+
     # Remove NaN in Peaks, Onsets, and Offsets
     waves_noNA = waves.copy()
     for feature in waves_noNA.keys():
@@ -947,11 +952,6 @@ def _ecg_delineator_peak(ecg, rpeaks=None, sampling_rate=1000):
         "ECG_T_Offsets": T_offsets,
     }
 
-    # Ensure that all indices are not larger than ECG signal indices
-    for _, value in info.items():
-        if value[-1] >= len(ecg):
-            value[-1] = np.nan
-
     # Return info dictionary
     return info
 

neurokit2/ecg/ecg_segment.py

@@ -59,6 +59,11 @@ def ecg_segment(ecg_cleaned, rpeaks=None, sampling_rate=1000, show=False):
         epochs_end=epochs_end,
     )
 
+    # pad last heartbeat with nan so that segments are equal length
+    last_heartbeat_key = str(np.max(np.array(list(heartbeats.keys()), dtype=int)))
+    after_last_index = heartbeats[last_heartbeat_key]["Index"] < len(ecg_cleaned)
+    heartbeats[last_heartbeat_key].loc[after_last_index, "Signal"] = np.nan
+
     if show:
         heartbeats_plot = epochs_to_df(heartbeats)
         heartbeats_pivoted = heartbeats_plot.pivot(index="Time", columns="Label", values="Signal")