Source code for neurokit2.complexity.complexity_rqa

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from ..misc import find_groups
from .entropy_shannon import entropy_shannon
from .optim_complexity_tolerance import complexity_tolerance
from .utils_recurrence_matrix import recurrence_matrix

[docs] def complexity_rqa( signal, dimension=3, delay=1, tolerance="sd", min_linelength=2, method="python", show=False ): """**Recurrence Quantification Analysis (RQA)** A :func:`recurrence plot <recurrence_matrix>` is based on a time-delay embedding representation of a signal and is a 2D depiction of when a system revisits a state that is has been in the past. Recurrence quantification analysis (RQA) is a method of complexity analysis for the investigation of dynamical systems. It quantifies the number and duration of recurrences of a dynamical system presented by its phase space trajectory. .. figure:: ../img/douglas2022c.png :alt: Illustration of RQA (Douglas et al., 2022). Features include: * **Recurrence rate (RR)**: Proportion of points that are labelled as recurrences. Depends on the radius *r*. * **Determinism (DET)**: Proportion of recurrence points which form diagonal lines. Indicates autocorrelation. * **Divergence (DIV)**: The inverse of the longest diagonal line length (*LMax*). * **Laminarity (LAM)**: Proportion of recurrence points which form vertical lines. Indicates the amount of laminar phases (intermittency). * **Trapping Time (TT)**: Average length of vertical black lines. * **L**: Average length of diagonal black lines. Average duration that a system is staying in the same state. * **LEn**: Entropy of diagonal lines lengths. * **VMax**: Longest vertical line length. * **VEn**: Entropy of vertical lines lengths. * **W**: Average white vertical line length. * **WMax**: Longest white vertical line length. * **WEn**: Entropy of white vertical lines lengths. * **DeteRec**: The ratio of determinism / recurrence rate. * **LamiDet**: The ratio of laminarity / determinism. * **DiagRec**: Diagonal Recurrence Rates, capturing the magnitude of autocorrelation at different lags, which is related to fractal fluctuations. See Tomashin et al. (2022), approach 3. .. note:: More feature exist for RQA, such as the `trend < DynamicalSystems.jl/dev/rqa/quantification/#RecurrenceAnalysis.trend>`_. We would like to add them, but we need help. Get in touch if you're interested! 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. tolerance : float Tolerance (often denoted as *r*), distance to consider two data points as similar. If ``"sd"`` (default), will be set to :math:`0.2 * SD_{signal}`. See :func:`complexity_tolerance` to estimate the optimal value for this parameter. min_linelength : int Minimum length of diagonal and vertical lines. Default to 2. method : str Can be ``"pyrqa"`` to use the *PyRQA* package (requires to install it first). show : bool Visualise recurrence matrix. Returns ---------- rqa : DataFrame The RQA results. info : dict A dictionary containing additional information regarding the parameters used to compute RQA. Examples ---------- .. ipython:: python import neurokit2 as nk signal = nk.signal_simulate(duration=5, sampling_rate=100, frequency=[5, 6, 7], noise=0.2) # RQA @savefig p_complexity_rqa1.png scale=100% results, info = nk.complexity_rqa(signal, tolerance=1, show=True) @suppress plt.close() .. ipython:: python results # Compare to PyRQA # results1, info = nk.complexity_rqa(signal, tolerance=1, show=True, method = "pyrqa") References ---------- * Rawald, T., Sips, M., Marwan, N., & Dransch, D. (2014). Fast computation of recurrences in long time series. In Translational Recurrences (pp. 17-29). Springer, Cham. * Tomashin, A., Leonardi, G., & Wallot, S. (2022). Four Methods to Distinguish between Fractal Dimensions in Time Series through Recurrence Quantification Analysis. Entropy, 24(9), 1314. """ info = { "Tolerance": complexity_tolerance( signal, method=tolerance, delay=delay, dimension=dimension )[0] } if method == "pyrqa": data = _complexity_rqa_pyrqa( signal, delay=delay, dimension=dimension, tolerance=info["Tolerance"], linelength=min_linelength, ) rc = np.flip(data.pop("Recurrence_Matrix"), axis=0) info["Recurrence_Matrix"] = rc else: # Get recurrence matrix (rm) rc, dm = recurrence_matrix( signal, delay=delay, dimension=dimension, tolerance=info["Tolerance"], ) info["Recurrence_Matrix"] = rc info["Distance_Matrix"] = dm # Compute features data = _complexity_rqa_features(rc, min_linelength=min_linelength) data = pd.DataFrame(data, index=[0]) if show is True: try: plt.imshow(rc, cmap="Greys") # Flip the matrix to match traditional RQA representation plt.gca().invert_yaxis() plt.title("Recurrence Matrix") plt.ylabel("Time") plt.xlabel("Time") except MemoryError as e: raise MemoryError( "NeuroKit error: complexity_rqa(): the recurrence plot is too large to display. ", "You can recover the matrix from the parameters and try to display parts of it.", ) from e return data, info
def _complexity_rqa_features(rc, min_linelength=2): """Compute recurrence rate from a recurrence matrix (rc).""" width = len(rc) # Recurrence Rate (RR) # -------------------------------------------------- # Indices of the lower triangular (without the diagonal) idx = np.tril_indices(width, k=-1) # Compute percentage data = {"RecurrenceRate": (rc[idx].sum()) / len(rc[idx])} # Find diagonale lines # -------------------------------------------------- diag_lines = [] recdiag = np.zeros(width) # All diagonals except the main one (0) for i in range(1, width): diag = np.diagonal(rc, offset=i) # Get diagonal recdiag[i - 1] = np.sum(diag) / len(diag) diag = find_groups(diag) # Split into consecutives diag_lines.extend([diag[i] for i in range(len(diag)) if diag[i][0] == 1]) # Store 1s # Diagonal Recurrence Rates (Diag %REC) # Tomashin et al. (2022) distance = np.arange(1, width + 1)[recdiag > 0] recdiag = recdiag[recdiag > 0] if len(recdiag) > 2: data["DiagRec"] = np.polyfit(np.log2(distance), np.log2(recdiag), 1)[0] # plt.loglog(distance, recdiag) else: data["DiagRec"] = np.nan # Get lengths diag_lengths = np.array([len(i) for i in diag_lines]) # Exclude small diagonals (> 1) diag_lengths = diag_lengths[np.where(diag_lengths >= min_linelength)[0]] # Compute features if data["RecurrenceRate"] == 0: data["Determinism"] = np.nan data["DeteRec"] = np.nan else: data["Determinism"] = diag_lengths.sum() / rc[idx].sum() data["DeteRec"] = data["Determinism"] / data["RecurrenceRate"] data["L"] = 0 if len(diag_lengths) == 0 else np.mean(diag_lengths) data["Divergence"] = np.nan if len(diag_lengths) == 0 else 1 / np.max(diag_lengths) data["LEn"] = entropy_shannon( freq=np.unique(diag_lengths, return_counts=True)[1], base=np.e, )[0] # Find vertical lines # -------------------------------------------------- black_lines = [] white_lines = [] for i in range(width - 1): verti = rc[i, i + 1 :] verti = find_groups(verti) black_lines.extend([verti[i] for i in range(len(verti)) if verti[i][0] == 1]) white_lines.extend([verti[i] for i in range(len(verti)) if verti[i][0] == 0]) # Get lengths black_lengths = np.array([len(i) for i in black_lines]) white_lengths = np.array([len(i) for i in white_lines]) # Exclude small lines (> 1) black_lengths = black_lengths[np.where(black_lengths >= min_linelength)[0]] white_lengths = white_lengths[np.where(white_lengths >= min_linelength)[0]] # Compute features if rc[idx].sum() == 0: data["Laminarity"] = np.nan else: data["Laminarity"] = black_lengths.sum() / rc[idx].sum() if data["Determinism"] == 0 or np.isnan(data["Determinism"]): data["LamiDet"] = np.nan else: data["Laminarity"] / data["Determinism"] data["TrappingTime"] = 0 if len(black_lengths) == 0 else np.nanmean(black_lengths) data["VMax"] = 0 if len(black_lengths) == 0 else np.nanmax(black_lengths) data["VEn"] = entropy_shannon( freq=np.unique(black_lengths, return_counts=True)[1], base=np.e, )[0] data["W"] = 0 if len(white_lengths) == 0 else np.nanmean(white_lengths) data["WMax"] = 0 if len(white_lengths) == 0 else np.nanmax(white_lengths) data["WEn"] = entropy_shannon( freq=np.unique(white_lengths, return_counts=True)[1], base=np.e, )[0] return data # ============================================================================= # PyRQA # ============================================================================= def _complexity_rqa_pyrqa(signal, dimension=3, delay=1, tolerance=0.1, linelength=2): """Compute recurrence rate (imported in complexity_rqa)""" # Try loading pyrqa try: import pyrqa.analysis_type import pyrqa.computation import pyrqa.image_generator import pyrqa.metric import pyrqa.neighbourhood import pyrqa.settings import pyrqa.time_series except (ModuleNotFoundError, ImportError) as e: raise ImportError( "NeuroKit error: complexity_rqa(): the 'pyrqa' module is required for this function to run. ", "Please install it first (`pip install PyRQA`).", ) from e # Get neighbourhood r = pyrqa.neighbourhood.FixedRadius(tolerance) # Convert signal to time series signal = pyrqa.time_series.TimeSeries(signal, embedding_dimension=dimension, time_delay=delay) settings = pyrqa.settings.Settings( signal, analysis_type=pyrqa.analysis_type.Classic, neighbourhood=r, similarity_measure=pyrqa.metric.EuclideanMetric, theiler_corrector=1, ) # RQA features rqa = pyrqa.computation.RQAComputation.create(settings, verbose=False).run() # Minimum line lengths rqa.min_diagonal_line_length = linelength rqa.min_vertical_line_length = linelength rqa.min_white_vertical_line_length = linelength rp = pyrqa.computation.RPComputation.create(settings, verbose=False).run() return { "RecurrenceRate": rqa.recurrence_rate, "Determinism": rqa.determinism, "Divergence": rqa.divergence, "Laminarity": rqa.laminarity, "TrappingTime": rqa.trapping_time, "DeteRec": rqa.determinism / rqa.recurrence_rate, "LamiDet": rqa.laminarity / rqa.determinism, "L": rqa.average_diagonal_line, "LEn": rqa.entropy_diagonal_lines, "VMax": rqa.longest_vertical_line, "VEn": rqa.entropy_vertical_lines, "W": rqa.average_white_vertical_line, "WMax": rqa.longest_white_vertical_line, "W_div": rqa.longest_white_vertical_line_inverse, "WEn": rqa.entropy_white_vertical_lines, "Recurrence_Matrix": rp.recurrence_matrix_reverse, # recurrence_matrix_reverse_normalized }