Source code for neurokit2.complexity.optim_complexity_k

from warnings import warn

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

from ..misc import NeuroKitWarning, find_plateau

[docs] def complexity_k(signal, k_max="max", show=False): """**Automated selection of k for Higuchi Fractal Dimension (HFD)** The optimal *k-max* is computed based on the point at which HFD values plateau for a range of *k-max* values (see Vega, 2015). Parameters ---------- signal : Union[list, np.array, pd.Series] The signal (i.e., a time series) in the form of a vector of values. k_max : Union[int, str, list], optional Maximum number of interval times (should be greater than or equal to 3) to be tested. If ``max``, it selects the maximum possible value corresponding to half the length of the signal. show : bool Visualise the slope of the curve for the selected kmax value. Returns -------- k : float The optimal kmax of the time series. info : dict A dictionary containing additional information regarding the parameters used to compute optimal kmax. See Also -------- fractal_higuchi Examples ---------- .. ipython:: python import neurokit2 as nk signal = nk.signal_simulate(duration=2, sampling_rate=100, frequency=[5, 6], noise=0.5) @savefig p_complexity_k1.png scale=100% k_max, info = nk.complexity_k(signal, k_max='default', show=True) @suppress plt.close() .. ipython:: python k_max References ---------- * Higuchi, T. (1988). Approach to an irregular time series on the basis of the fractal theory. Physica D: Nonlinear Phenomena, 31(2), 277-283. * Vega, C. F., & Noel, J. (2015, June). Parameters analyzed of Higuchi's fractal dimension for EEG brain signals. In 2015 Signal Processing Symposium (SPSympo) (pp. 1-5). IEEE. https:// """ # Get the range of k-max values to be tested # ------------------------------------------ if isinstance(k_max, str): # e.g., "default" # upper limit for k value (max possible value) k_max = int(np.floor(len(signal) / 2)) # so that normalizing factor is positive if isinstance(k_max, int): kmax_range = np.arange(2, k_max + 1) elif isinstance(k_max, (list, np.ndarray, pd.Series)): kmax_range = np.array(k_max) else: warn( "k_max should be an int or a list of values of kmax to be tested.", category=NeuroKitWarning, ) # Compute the slope for each kmax value # -------------------------------------- vectorized_k_slope = np.vectorize(_complexity_k_slope, excluded=[1]) slopes, intercepts, info = vectorized_k_slope(kmax_range, signal) # k_values = [d["k_values"] for d in info] average_values = [d["average_values"] for d in info] # Find plateau (the saturation point of slope) # -------------------------------------------- optimal_point = find_plateau(slopes, show=False) if optimal_point is not None: kmax_optimal = kmax_range[optimal_point] else: kmax_optimal = np.max(kmax_range) warn( "The optimal kmax value detected is 2 or less. There may be no plateau in this case. " + f"You can inspect the plot by set `show=True`. We will return optimal k_max = {kmax_optimal} (the max).", category=NeuroKitWarning, ) # Plot if show: _complexity_k_plot(kmax_range, slopes, kmax_optimal, ax=None) # Return optimal tau and info dict return kmax_optimal, { "Values": kmax_range, "Scores": slopes, "Intercepts": intercepts, "Average_Values": average_values, }
# ============================================================================= # Utilities # ============================================================================= def _complexity_k_Lk(k, signal): n = len(signal) # Step 1: construct k number of new time series for range of k_values from 1 to kmax k_subrange = np.arange(1, k + 1) # where m = 1, 2... k idx = np.tile(np.arange(0, len(signal), k), (k, 1)).astype(float) idx += np.tile(np.arange(0, k), (idx.shape[1], 1)).T mask = idx >= len(signal) idx[mask] = 0 sig_values = signal[idx.astype(int)].astype(float) sig_values[mask] = np.nan # Step 2: Calculate length Lm(k) of each curve normalization = (n - 1) / (np.floor((n - k_subrange) / k).astype(int) * k) sets = (np.nansum(np.abs(np.diff(sig_values)), axis=1) * normalization) / k # Step 3: Compute average value over k sets of Lm(k) return np.sum(sets) / k def _complexity_k_slope(kmax, signal, k_number="max"): if k_number == "max": k_values = np.arange(1, kmax + 1) else: k_values = np.unique(np.linspace(1, kmax + 1, k_number).astype(int)) # Step 3 of Vega & Noel (2015) vectorized_Lk = np.vectorize(_complexity_k_Lk, excluded=[1]) # Compute length of the curve, Lm(k) average_values = vectorized_Lk(k_values, signal) # Slope of best-fit line through points (slope equal to FD) slope, intercept = -np.polyfit(np.log(k_values), np.log(average_values), 1) return slope, intercept, {"k_values": k_values, "average_values": average_values} # ============================================================================= # Plotting # ============================================================================= def _complexity_k_plot(k_range, slope_values, k_optimal, ax=None): # Prepare plot if ax is None: fig, ax = plt.subplots() else: fig = None ax.set_title("Optimization of $k_{max}$ parameter") ax.set_xlabel("$k_{max}$ values") ax.set_ylabel("Higuchi Fractal Dimension (HFD) values") colors =, 1, len(k_range))) # if single time series ax.plot(k_range, slope_values, color="#2196F3", zorder=1) for i, _ in enumerate(k_range): ax.scatter(k_range[i], slope_values[i], color=colors[i], marker="o", zorder=2) ax.axvline(x=k_optimal, color="#E91E63", label="Optimal $k_{max}$: " + str(k_optimal)) ax.legend(loc="upper right") return fig