Source code for neurokit2.complexity.complexity_lyapunov

# -*- coding: utf-8 -*-
from warnings import warn

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sklearn.metrics.pairwise
import sklearn.neighbors

from ..misc import NeuroKitWarning, find_knee
from ..signal.signal_psd import signal_psd
from .utils_complexity_embedding import complexity_embedding

[docs] def complexity_lyapunov( signal, delay=1, dimension=2, method="rosenstein1993", separation="auto", **kwargs, ): """**(Largest) Lyapunov Exponent (LLE)** Lyapunov exponents (LE) describe the rate of exponential separation (convergence or divergence) of nearby trajectories of a dynamical system. It is a measure of sensitive dependence on initial conditions, i.e. how quickly two nearby states diverge. A system can have multiple LEs, equal to the number of the dimensionality of the phase space, and the largest LE value, "LLE" is often used to determine the overall predictability of the dynamical system. Different algorithms exist to estimate these indices: * **Rosenstein et al.'s (1993)** algorithm was designed for calculating LLEs from small datasets. The time series is first reconstructed using a delay-embedding method, and the closest neighbour of each vector is computed using the euclidean distance. These two neighbouring points are then tracked along their distance trajectories for a number of data points. The slope of the line using a least-squares fit of the mean log trajectory of the distances gives the final LLE. * **Makowski** is a custom modification of Rosenstein's algorithm, using KDTree for more efficient nearest neighbors computation. Additionally, the LLE is computed as the slope up to the changepoint of divergence rate (the point where it flattens out), making it more robust to the length trajectory parameter. * **Eckmann et al. (1986)** computes LEs by first reconstructing the time series using a delay-embedding method, and obtains the tangent that maps to the reconstructed dynamics using a least-squares fit, where the LEs are deduced from the tangent maps. .. warning:: The **Eckman (1986)** method currently does not work. Please help us fixing it by double checking the code, the paper and helping us figuring out what's wrong. Overall, we would like to improve this function to return for instance all the exponents (Lyapunov spectrum), implement newer and faster methods (e.g., Balcerzak, 2018, 2020), etc. If you're interested in helping out with this, please get in touch! 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. If method is ``"eckmann1986"``, larger values for dimension are recommended. method : str The method that defines the algorithm for computing LE. Can be one of ``"rosenstein1993"``, ``"makowski"``, or ``"eckmann1986"``. len_trajectory : int Applies when method is ``"rosenstein1993"``. The number of data points in which neighboring trajectories are followed. matrix_dim : int Applies when method is ``"eckmann1986"``. Corresponds to the number of LEs to return. min_neighbors : int, str Applies when method is ``"eckmann1986"``. Minimum number of neighbors. If ``"default"``, ``min(2 * matrix_dim, matrix_dim + 4)`` is used. **kwargs : optional Other arguments to be passed to ``signal_psd()`` for calculating the minimum temporal separation of two neighbors. Returns -------- lle : float An estimate of the largest Lyapunov exponent (LLE) if method is ``"rosenstein1993"``, and an array of LEs if ``"eckmann1986"``. info : dict A dictionary containing additional information regarding the parameters used to compute LLE. Examples ---------- .. ipython:: python import neurokit2 as nk signal = nk.signal_simulate(duration=5, sampling_rate=100, frequency=[5, 8], noise=0.1) # Rosenstein's method @savefig p_complexity_lyapunov1.png scale=100% lle, info = nk.complexity_lyapunov(signal, method="rosenstein", show=True) @suppress plt.close() lle # Makowski's change-point method @savefig p_complexity_lyapunov2.png scale=100% lle, info = nk.complexity_lyapunov(signal, method="makowski", show=True) @suppress plt.close() # Eckman's method is broken. Please help us fix-it! # lle, info = nk.complexity_lyapunov(signal, dimension=2, method="eckmann1986") References ---------- * Rosenstein, M. T., Collins, J. J., & De Luca, C. J. (1993). A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65(1-2), 117-134. * Eckmann, J. P., Kamphorst, S. O., Ruelle, D., & Ciliberto, S. (1986). Liapunov exponents from time series. Physical Review A, 34(6), 4971. """ # 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." ) # Compute Minimum temporal separation between two neighbors # ----------------------------------------------------------- # Rosenstein (1993) finds a suitable value by calculating the mean period of the data, # obtained by the reciprocal of the mean frequency of the power spectrum. # "We impose the additional constraint that nearest neighbors have a temporal separation # greater than the mean period of the time series: This allows us to consider each pair of # neighbors as nearby initial conditions for different trajectories." # "We estimated the mean period as the reciprocal of the mean frequency of the power spectrum, # although we expect any comparable estimate, e.g., using the median frequency of the magnitude # spectrum, to yield equivalent results." if separation == "auto": # Actual sampling rate does not matter psd = signal_psd( signal, sampling_rate=1000, method="fft", normalize=False, show=False ) mean_freq = np.sum(psd["Power"] * psd["Frequency"]) / np.sum(psd["Power"]) # 1 / mean_freq = seconds per cycle separation = int(np.ceil(1 / mean_freq * 1000)) else: assert isinstance(separation, int), "'separation' should be an integer." # Run algorithm # ---------------- # Method method = method.lower() if method in ["rosenstein", "rosenstein1993"]: le, parameters = _complexity_lyapunov_rosenstein( signal, delay, dimension, separation, **kwargs ) elif method in ["makowski"]: le, parameters = _complexity_lyapunov_makowski( signal, delay, dimension, separation, **kwargs ) elif method in ["eckmann", "eckmann1986", "eckmann1986"]: le, parameters = _complexity_lyapunov_eckmann( signal, dimension=dimension, separation=separation, ) else: raise ValueError( "NeuroKit error: complexity_lyapunov(): 'method' should be one of " " 'rosenstein1993', 'makowski', 'eckmann1986'." ) # Store params info = { "Dimension": dimension, "Delay": delay, "Separation": separation, "Method": method, } info.update(parameters) return le, info
# ============================================================================= # Methods # ============================================================================= def _complexity_lyapunov_makowski( signal, delay=1, dimension=2, separation=1, max_length="auto", show=False, ): # Store parameters info = { "Dimension": dimension, "Delay": delay, } # Embedding embedded = complexity_embedding(signal, delay=delay, dimension=dimension) n = len(embedded) # Set the maxiimum trajectory length to 10 times the delay if max_length == "auto": max_length = int(delay * 10) if max_length >= n / 2: max_length = n // 2 # Create KDTree and query for nearest neighbors tree = sklearn.neighbors.KDTree(embedded, metric="euclidean") # Query for nearest neighbors. To ensure we get a neighbor outside of the `separation`, # k=1 is the point itself, k=2 is the nearest neighbor, and k=3 is the second nearest neighbor. idx = tree.query(embedded, k=2 + separation, return_distance=False) # The neighbor outside the `separation` region will be the last one in the returned list. idx = idx[:, -1] # Compute the average divergence for each trajectory length trajectories = np.zeros(max_length) for k in range(1, max_length + 1): valid = np.where((np.arange(n - k) + k < n) & (idx[: n - k] + k < n))[0] if valid.size == 0: trajectories[k - 1] = -np.inf continue divergences = np.linalg.norm( embedded[valid + k] - embedded[idx[valid] + k], axis=1, ) divergences = divergences[divergences > 0] if len(divergences) == 0: trajectories[k - 1] = np.nan else: trajectories[k - 1] = np.mean(np.log(divergences)) # Change point x_axis = range(1, len(trajectories) + 1) knee = find_knee(y=trajectories, x=x_axis, show=False, verbose=False) info["Divergence_Rate"] = trajectories info["Changepoint"] = knee # Linear fit slope, intercept = np.polyfit(x_axis[0:knee], trajectories[0:knee], 1) if show is True: plt.plot(np.arange(1, len(trajectories) + 1), trajectories) plt.axvline(knee, color="red", label="Changepoint", linestyle="--") plt.axline( (0, intercept), slope=slope, color="orange", label="Least-squares Fit" ) plt.ylim(bottom=np.min(trajectories)) plt.ylabel("Divergence Rate") plt.title(f"Largest Lyapunov Exponent (slope of the line) = {slope:.3f}") plt.legend() return slope, info def _complexity_lyapunov_rosenstein( signal, delay=1, dimension=2, separation=1, len_trajectory=20, show=False, **kwargs ): # 1. Check that sufficient data points are available # Minimum length required to find single orbit vector min_len = (dimension - 1) * delay + 1 # We need len_trajectory orbit vectors to follow a complete trajectory min_len += len_trajectory - 1 # we need tolerance * 2 + 1 orbit vectors to find neighbors for each min_len += separation * 2 + 1 # Sanity check if len(signal) < min_len: warn( f"for dimension={dimension}, delay={delay}, separation={separation} and " f"len_trajectory={len_trajectory}, you need at least {min_len} datapoints in your" " time series.", category=NeuroKitWarning, ) # Embedding embedded = complexity_embedding(signal, delay=delay, dimension=dimension) m = len(embedded) # Construct matrix with pairwise distances between vectors in orbit dists = sklearn.metrics.pairwise.euclidean_distances(embedded) for i in range(m): # Exclude indices within separation dists[i, max(0, i - separation) : i + separation + 1] = np.inf # Find indices of nearest neighbours ntraj = m - len_trajectory + 1 min_dist_indices = np.argmin( dists[:ntraj, :ntraj], axis=1 ) # exclude last few indices min_dist_indices = min_dist_indices.astype(int) # Follow trajectories of neighbour pairs for len_trajectory data points trajectories = np.zeros(len_trajectory) for k in range(len_trajectory): divergence = dists[(np.arange(ntraj) + k, min_dist_indices + k)] dist_nonzero = np.where(divergence != 0)[0] if len(dist_nonzero) == 0: trajectories[k] = -np.inf else: # Get average distances of neighbour pairs along the trajectory trajectories[k] = np.mean(np.log(divergence[dist_nonzero])) divergence_rate = trajectories[np.isfinite(trajectories)] # LLE obtained by least-squares fit to average line slope, intercept = np.polyfit( np.arange(1, len(divergence_rate) + 1), divergence_rate, 1 ) # Store info parameters = { "Trajectory_Length": len_trajectory, "Divergence_Rate": divergence_rate, } if show is True: plt.plot(np.arange(1, len(divergence_rate) + 1), divergence_rate) plt.axline( (0, intercept), slope=slope, color="orange", label="Least-squares Fit" ) plt.ylabel("Divergence Rate") plt.title(f"Largest Lyapunov Exponent (slope of the line) = {slope:.3f}") plt.legend() return slope, parameters def _complexity_lyapunov_eckmann( signal, dimension=2, separation=None, matrix_dim=4, min_neighbors="default", tau=1 ): """TODO: check implementation From """ # Prepare parameters if min_neighbors == "default": min_neighbors = min(2 * matrix_dim, matrix_dim + 4) m = (dimension - 1) // (matrix_dim - 1) # minimum length required to find single orbit vector min_len = dimension # we need to follow each starting point of an orbit vector for m more steps min_len += m # we need separation * 2 + 1 orbit vectors to find neighbors for each min_len += separation * 2 # we need at least min_nb neighbors for each orbit vector min_len += min_neighbors # Sanity check if len(signal) < min_len: warn( f"for dimension={dimension}, separation={separation}, " f"matrix_dim={matrix_dim} and min_neighbors={min_neighbors}, " f"you need at least {min_len} datapoints in your time series.", category=NeuroKitWarning, ) # Storing of LEs lexp = np.zeros(matrix_dim) lexp_counts = np.zeros(matrix_dim) old_Q = np.identity(matrix_dim) # We need to be able to step m points further for the beta vector vec = signal if m == 0 else signal[:-m] # If m==0, return full signal # Reconstruction using time-delay method embedded = complexity_embedding(vec, delay=1, dimension=dimension) distances = sklearn.metrics.pairwise_distances(embedded, metric="chebyshev") for i in range(len(embedded)): # exclude difference of vector to itself and those too close in time distances[i, max(0, i - separation) : i + separation + 1] = np.inf # index of furthest nearest neighbour neighbour_furthest = np.argsort(distances[i])[min_neighbors - 1] # get neighbors within the radius r = distances[i][neighbour_furthest] neighbors = np.where(distances[i] <= r)[ 0 ] # should have length = min_neighbours # Find matrix T_i (matrix_dim * matrix_dim) that sends points from neighbourhood of x(i) to x(i+1) vec_beta = signal[neighbors + matrix_dim * m] - signal[i + matrix_dim * m] matrix = np.array([signal[j : j + dimension : m] for j in neighbors]) # x(j) matrix -= signal[i : i + dimension : m] # x(j) - x(i) # form matrix T_i t_i = np.zeros((matrix_dim, matrix_dim)) t_i[:-1, 1:] = np.identity(matrix_dim - 1) t_i[-1] = np.linalg.lstsq(matrix, vec_beta, rcond=-1)[ 0 ] # least squares solution # QR-decomposition of T * old_Q mat_Q, mat_R = np.linalg.qr(, old_Q)) # force diagonal of R to be positive sign_diag = np.sign(np.diag(mat_R)) sign_diag[np.where(sign_diag == 0)] = 1 sign_diag = np.diag(sign_diag) mat_Q =, sign_diag) mat_R =, mat_R) old_Q = mat_Q # successively build sum for Lyapunov exponents diag_R = np.diag(mat_R) # filter zeros in mat_R (would lead to -infs) positive_elements = np.where(diag_R > 0) lexp_i = np.zeros(len(diag_R)) lexp_i[positive_elements] = np.log(diag_R[positive_elements]) lexp_i[np.where(diag_R == 0)] = np.inf lexp[positive_elements] += lexp_i[positive_elements] lexp_counts[positive_elements] += 1 # normalize exponents over number of individual mat_Rs idx = np.where(lexp_counts > 0) lexp[idx] /= lexp_counts[idx] lexp[np.where(lexp_counts == 0)] = np.inf # normalize with respect to tau lexp /= tau # take m into account lexp /= m parameters = {"Minimum_Neighbors": min_neighbors} return lexp, parameters