<gh_stars>0
# -*- coding: utf-8 -*-
import itertools
from warnings import warn

import matplotlib
import matplotlib.collections
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.signal
import scipy.spatial
import scipy.stats

from ..misc import NeuroKitWarning, find_closest
from ..signal import (
    signal_autocor,
    signal_findpeaks,
    signal_psd,
    signal_surrogate,
    signal_zerocrossings,
)
from .entropy_kl import entropy_kl
from .information_mutual import mutual_information
from .utils_complexity_embedding import complexity_embedding


def complexity_delay(
    signal, delay_max=50, method="fraser1986", algorithm=None, show=False, **kwargs
):
    """**Automated selection of the optimal Delay (Tau)**

    The time delay (Tau :math:`\\tau`, also referred to as *Lag*) is one of the two critical
    parameters (the other being the :func:`Dimension <complexity_dimension>` *m*) involved in the
    construction of the time-delay embedding of a signal. It corresponds to the delay in samples
    between the original signal and its delayed version(s). In other words, how many samples do we
    consider between a given state of the signal and its closest past state.

    When :math:`\\tau` is smaller than the optimal theoretical value, consecutive coordinates of the
    system's state are correlated and the attractor is not sufficiently unfolded. Conversely, when
    :math:`\\tau` is larger than it should be, successive coordinates are almost independent,
    resulting in an uncorrelated and unstructured cloud of points.

    The selection of the parameters *delay* and :func:`*dimension* <complexity_dimension>` is a
    challenge. One approach is to select them (semi) independently (as dimension selection often
    requires the delay), using :func:`complexity_delay` and :func:`complexity_dimension`. However,
    some joint-estimation methods do exist, that attempt at finding the optimal delay and dimension
    at the same time.

    Several authors suggested different methods to guide the choice of the delay:

    * **Fraser and Swinney (1986)** suggest using the first local minimum of the mutual information
      between the delayed and non-delayed time series, effectively identifying a value of Tau for
      which they share the least information.
    * **Theiler (1990)** suggested to select Tau where the autocorrelation between the signal and
      its lagged version at Tau first crosses the value :math:`1/e`.
    * **Casdagli (1991)** suggests instead taking the first zero-crossing of the autocorrelation.
    * **Rosenstein (1993)** suggests to approximate the point where the autocorrelation function
      drops to :math:`(1 - 1/e)` of its maximum value.
    * **Rosenstein (1994)** suggests to the point close to 40% of the slope of the average
      displacement from the diagonal (ADFD).
    * **Kim (1999)** suggests estimating Tau using the correlation integral, called the C-C method,
      which has shown to agree with those obtained using the Mutual Information. This method
      makes use of a statistic within the reconstructed phase space, rather than analyzing the
      temporal evolution of the time series. However, computation times are significantly long for
      this method due to the need to compare every unique pair of pairwise vectors within the
      embedded signal per delay.
    * **Lyle (2021)** describes the "Symmetric Projection Attractor Reconstruction" (SPAR), where
      :math:`1/3` of the the dominant frequency (i.e., of the length of the average "cycle") can be
      a suitable value for approximately periodic data, and makes the attractor sensitive to
      morphological changes. See also `Aston's talk <https://youtu.be/GGrOJtcTcHA?t=730>`_. This
      method is also the fastest but might not be suitable for aperiodic signals.
      The ``algorithm`` argument (default to ``"fft"``) and will be passed as the ``method``
      argument of  ``signal_psd()``.

    **Joint-Methods for Delay and Dimension**

    * **Gautama (2003)** mentions that in practice, it is common to have a fixed time lag and to
      adjust the embedding dimension accordingly. As this can lead to large *m* values (and thus to
      embedded data of a large size) and thus, slow processing, they describe an optimisation
      method to jointly determine *m* and :math:`\\tau`, based on the **entropy ratio**.

    .. note::

      We would like to implement the joint-estimation by `Matilla-García et al. (2021)
      <https://www.mdpi.com/1099-4300/23/2/221>`_, please get in touch if you can help us!

    Parameters
    ----------
    signal : Union[list, np.array, pd.Series]
        The signal (i.e., a time series) in the form of a vector of values.
    delay_max : int
        The maximum time delay (Tau or lag) to test.
    method : str
        The method that defines what to compute for each tested value of Tau. Can be one of
        ``"fraser1986"``, ``"theiler1990"``, ``"casdagli1991"``, ``"rosenstein1993"``,
        ``"rosenstein1994"``, ``"kim1999"``, or ``"lyle2021"``.
    algorithm : str
        The method used to find the optimal value of Tau given the values computed by the method.
        If ``None`` (default), will select the algorithm according to the method. Modify only if you
        know what you are doing.
    show : bool
        If ``True``, will plot the metric values for each value of tau.
    **kwargs : optional
        Additional arguments to be passed for C-C method.

    Returns
    -------
    delay : int
        Optimal time delay.
    parameters : dict
        A dictionary containing additional information regarding the parameters used
        to compute optimal time-delay embedding.

    See Also
    ---------
    complexity, complexity_dimension, complexity_embedding, complexity_tolerance

    Examples
    ----------

    * **Example 1**: Comparison of different methods for estimating the optimal delay of an simple
      artificial signal.

    .. ipython:: python

      import neurokit2 as nk

      signal = nk.signal_simulate(duration=10, sampling_rate=100, frequency=[1, 1.5], noise=0.02)

      @savefig p_complexity_delay1.png scale=100%
      nk.signal_plot(signal)
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay2.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="fraser1986")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay3.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="theiler1990")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay4.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="casdagli1991")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay5.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="rosenstein1993")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay6.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="rosenstein1994")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay7.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="lyle2021")
      @suppress
      plt.close()

    * **Example 2**: Optimizing the delay and the dimension using joint-estimation methods.

    .. ipython:: python

      @savefig p_complexity_delay8.png scale=100%
      delay, parameters = nk.complexity_delay(
          signal,
          delay_max=np.arange(1, 30, 1),  # Can be an int or a list
          dimension_max=20,  # Can be an int or a list
          method="gautama2003",
          surrogate_n=5,  # Number of surrogate signals to generate
          surrogate_method="random",  # can be IAAFT, see nk.signal_surrogate()
          show=True)
      @suppress
      plt.close()

    .. ipython:: python

      # Optimal dimension
      dimension = parameters["Dimension"]
      dimension

    **Note**: A double-checking of that method would be appreciated! Please help us improve.

    * **Example 3**: Using a realistic signal.

    .. ipython:: python

      ecg = nk.ecg_simulate(duration=60*6, sampling_rate=200)
      signal = nk.ecg_rate(nk.ecg_peaks(ecg, sampling_rate=200),
                           sampling_rate=200,
                           desired_length=len(ecg))
      @savefig p_complexity_delay9.png scale=100%
      nk.signal_plot(signal)
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay10.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=1000, show=True)
      @suppress
      plt.close()

    References
    ------------
    * <NAME>., <NAME>., & <NAME>. (2021). Symmetric Projection Attractor Reconstruction:
      Sex Differences in the ECG. Frontiers in cardiovascular medicine, 1034.
    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.
    * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time
      series forecasting. IFAC Proceedings Volumes, 42(7), 316-320.
    * <NAME>., <NAME>., & <NAME>. (1993). A practical method for calculating
      largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65(1-2),
      117-134.
    * <NAME>., <NAME>., & <NAME>. (1994). Reconstruction expansion as a
      geometry-based framework for choosing proper delay times. Physica-Section D, 73(1), 82-98.
    * <NAME>., <NAME>., & <NAME>. (1999). Nonlinear dynamics, delay times, and embedding
      windows. Physica D: Nonlinear Phenomena, 127(1-2), 48-60.
    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.
    * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time
      series forecasting. IFAC Proceedings Volumes, 42(7), 316-320.

    """
    # Initalize vectors
    if isinstance(delay_max, int):
        tau_sequence = np.arange(1, delay_max + 1)
    else:
        tau_sequence = np.array(delay_max)

    # Method
    method = method.lower()
    if method in ["fraser", "fraser1986", "tdmi"]:
        metric = "Mutual Information"
        if algorithm is None:
            algorithm = "first local minimum"
    elif method in ["theiler", "theiler1990"]:
        metric = "Autocorrelation"
        if algorithm is None:
            algorithm = "first 1/e crossing"
    elif method in ["casdagli", "casdagli1991"]:
        metric = "Autocorrelation"
        if algorithm is None:
            algorithm = "first zero crossing"
    elif method in ["rosenstein", "rosenstein1994", "adfd"]:
        metric = "Displacement"
        if algorithm is None:
            algorithm = "closest to 40% of the slope"
    elif method in ["rosenstein1993"]:
        metric = "Autocorrelation (FFT)"
        if algorithm is None:
            algorithm = "first drop below 1-(1/e) of maximum"
    elif method in ["kim1999", "cc"]:
        metric = "Correlation Integral"
        if algorithm is None:
            algorithm = "first local minimum"
    elif method in ["aston2020", "lyle2021", "spar"]:
        return _complexity_delay_spar(signal, algorithm=algorithm, show=show, **kwargs)
    elif method in ["gautama2003", "gautama", "entropyratio"]:
        return _complexity_delay_entropyratio(signal, delay_max=delay_max, show=show, **kwargs)
    else:
        raise ValueError("NeuroKit error: complexity_delay(): 'method' not recognized.")

    # Get metric
    metric_values = _embedding_delay_metric(signal, tau_sequence, metric=metric)

    # Get optimal tau
    optimal = _embedding_delay_select(metric_values, algorithm=algorithm)
    if np.isnan(optimal):
        warn(
            "No optimal time delay is found. Nan is returned."
            " Consider using a higher `delay_max`.",
            category=NeuroKitWarning,
        )
        return optimal
    optimal = tau_sequence[optimal]

    if show is True:
        _embedding_delay_plot(
            signal,
            metric_values=metric_values,
            tau_sequence=tau_sequence,
            tau=optimal,
            metric=metric,
        )

    # Return optimal tau and info dict
    return optimal, {
        "Values": tau_sequence,
        "Scores": metric_values,
        "Algorithm": algorithm,
        "Metric": metric,
        "Method": method,
    }


# =============================================================================
# Methods
# =============================================================================
def _embedding_delay_select(metric_values, algorithm="first local minimum"):

    if algorithm in ["first local minimum (corrected)"]:
        # if immediately increasing, assume it is the first is the closest
        if np.diff(metric_values)[0] > 0:
            optimal = 0
        # take last value if continuously decreasing with no inflections
        elif all(np.diff(metric_values) < 0):
            optimal = len(metric_values) - 1
        else:
            # Find reversed peaks
            optimal = signal_findpeaks(
                -1 * metric_values, relative_height_min=0.1, relative_max=True
            )["Peaks"]
    elif algorithm == "first local minimum":
        # Find reversed peaks
        try:
            optimal = signal_findpeaks(
                -1 * metric_values, relative_height_min=0.1, relative_max=True
            )["Peaks"]
        except ValueError:
            warn(
                "First local minimum detection failed. Try setting "
                + "`algorithm = 'first local minimum (corrected)'` or using another method.",
                category=NeuroKitWarning,
            )

    elif algorithm == "first 1/e crossing":
        metric_values = metric_values - 1 / np.exp(1)
        optimal = signal_zerocrossings(metric_values)
    elif algorithm == "first zero crossing":
        optimal = signal_zerocrossings(metric_values)
    elif algorithm == "closest to 40% of the slope":
        slope = np.diff(metric_values) * len(metric_values)
        slope_in_deg = np.rad2deg(np.arctan(slope))
        optimal = np.where(slope_in_deg == find_closest(40, slope_in_deg))[0]
    elif algorithm == "first drop below 1-(1/e) of maximum":
        optimal = np.where(metric_values < np.max(metric_values) * (1 - 1.0 / np.e))[0][0]

    if not isinstance(optimal, (int, float, np.integer)):
        if len(optimal) != 0:
            optimal = optimal[0]
        else:
            optimal = np.nan

    return optimal


def _embedding_delay_metric(
    signal,
    tau_sequence,
    metric="Mutual Information",
    dimensions=[2, 3, 4, 5],
    r_vals=[0.5, 1.0, 1.5, 2.0],
):
    """Iterating through dimensions and r values is relevant only if metric used is Correlation Integral.
    For this method, either first zero crossing of the statistic averages or the first local
    minimum of deviations to obtain optimal tau. This implementation takes the latter since in practice,
    they are both in close proximity.
    """

    if metric == "Autocorrelation":
        values, _ = signal_autocor(signal)
        values = values[: len(tau_sequence)]  # upper limit

    elif metric == "Autocorrelation (FFT)":
        values, _ = signal_autocor(signal, demean=False, method="fft")
        values = values[: len(tau_sequence)]

    elif metric == "Correlation Integral":
        r_vals = [i * np.std(signal) for i in r_vals]
        # initiate empty list for storing
        # averages = np.zeros(len(tau_sequence))
        values = np.zeros(len(tau_sequence))
        for i, t in enumerate(tau_sequence):
            # average = 0
            change = 0
            for m in dimensions:
                # # find average of dependence statistic
                # for r in r_vals:
                #     s = _embedding_delay_cc_statistic(signal, delay=t, dimension=m, r=r)
                #     average += s

                # find average of statistic deviations across r_vals
                deviation = _embedding_delay_cc_deviation_max(
                    signal, delay=t, dimension=m, r_vals=r_vals
                )
                change += deviation
            # averages[i] = average / 16
            values[i] = change / 4

    else:
        values = np.zeros(len(tau_sequence))

        # Loop through taus and compute all scores values
        for i, current_tau in enumerate(tau_sequence):
            embedded = complexity_embedding(signal, delay=current_tau, dimension=2)
            if metric == "Mutual Information":
                values[i] = mutual_information(
                    embedded[:, 0], embedded[:, 1], normalized=True, method="shannon"
                )
            if metric == "Displacement":
                dimension = 2

                # Reconstruct with zero time delay.
                tau0 = embedded[:, 0].repeat(dimension).reshape(len(embedded), dimension)
                dist = np.asarray(
                    [scipy.spatial.distance.euclidean(i, j) for i, j in zip(embedded, tau0)]
                )
                values[i] = np.mean(dist)

    return values


# =============================================================================
# Internals for C-C method, Kim et al. (1999)
# =============================================================================


def _embedding_delay_cc_integral_sum(signal, dimension=3, delay=10, r=0.02):
    """Correlation integral is a cumulative distribution function, which denotes
    the probability of distance between any pairs of points in phase space
    not greater than the specified `r`.
    """

    # Embed signal
    embedded = complexity_embedding(signal, delay=delay, dimension=dimension, show=False)
    M = embedded.shape[0]

    # Prepare indices for comparing all unique pairwise vectors
    combinations = list(itertools.combinations(range(0, M), r=2))
    first_index, second_index = np.transpose(combinations)[0], np.transpose(combinations)[1]

    vectorized_integral = np.vectorize(_embedding_delay_cc_integral, excluded=["embedded", "r"])
    integral = np.sum(
        vectorized_integral(
            first_index=first_index, second_index=second_index, embedded=embedded, r=r
        )
    )

    return integral


def _embedding_delay_cc_integral(first_index, second_index, embedded, r=0.02):
    M = embedded.shape[0]  # Number of embedded points
    diff = np.linalg.norm(embedded[first_index] - embedded[second_index], ord=np.inf)  # sup-norm
    h = np.heaviside(r - diff, 1)
    integral = (2 / (M * (M - 1))) * h  # find average

    return integral


def _embedding_delay_cc_statistic(signal, dimension=3, delay=10, r=0.02):
    """The dependence statistic as the serial correlation of a nonlinear time series."""

    # create disjoint time series
    series = [signal[i - 1 :: delay] for i in range(1, delay + 1)]

    statistic = 0
    for sub_series in series:
        diff = _embedding_delay_cc_integral_sum(
            sub_series, dimension=dimension, delay=delay, r=r
        ) - ((_embedding_delay_cc_integral_sum(signal, dimension=1, delay=delay, r=r)) ** dimension)
        statistic += diff

    return statistic / delay


def _embedding_delay_cc_deviation_max(signal, r_vals=[0.5, 1.0, 1.5, 2.0], delay=10, dimension=3):
    """A measure of the variation of the dependence statistic with r using
    several representative values of r.
    """
    vectorized_deviation = np.vectorize(
        _embedding_delay_cc_deviation, excluded=["signal", "delay", "dimension"]
    )
    deviations = vectorized_deviation(
        signal=signal, r_vals=r_vals, delay=delay, dimension=dimension
    )

    return np.max(deviations) - np.min(deviations)


def _embedding_delay_cc_deviation(signal, r_vals=[0.5, 1.0, 1.5, 2.0], delay=10, dimension=3):
    return _embedding_delay_cc_statistic(signal, delay=delay, dimension=dimension, r=r_vals)


# =============================================================================
# Plotting Generics
# =============================================================================
def _embedding_delay_plot(
    signal,
    metric_values,
    tau_sequence,
    tau=1,
    metric="Mutual Information",
    ax0=None,
    ax1=None,
    plot="2D",
):

    # Prepare figure
    if ax0 is None and ax1 is None:
        fig = plt.figure(constrained_layout=False)
        spec = matplotlib.gridspec.GridSpec(
            ncols=1, nrows=2, height_ratios=[1, 3], width_ratios=[2]
        )
        ax0 = fig.add_subplot(spec[0])
        if plot == "2D":
            ax1 = fig.add_subplot(spec[1])
        elif plot == "3D":
            ax1 = fig.add_subplot(spec[1], projection="3d")
    else:
        fig = None

    ax0.set_title("Optimization of Delay")
    ax0.set_xlabel("Time Delay")
    ax0.set_ylabel(metric)
    ax0.plot(tau_sequence, metric_values, color="#FFC107")
    ax0.axvline(x=tau, color="#E91E63", label="Optimal delay: " + str(tau))
    ax0.legend(loc="upper right")
    ax1.set_title("Attractor")
    ax1.set_xlabel("Signal [i]")
    ax1.set_ylabel("Signal [i-" + str(tau) + "]")

    # Get data points, set axis limits
    embedded = complexity_embedding(signal, delay=tau, dimension=3)
    x = embedded[:, 0]
    y = embedded[:, 1]
    z = embedded[:, 2]
    ax1.set_xlim(x.min(), x.max())
    ax1.set_ylim(x.min(), x.max())

    # Colors
    norm = plt.Normalize(z.min(), z.max())
    cmap = plt.get_cmap("plasma")
    colors = cmap(norm(x))

    # Attractor for 2D vs 3D
    if plot == "2D":
        points = np.array([x, y]).T.reshape(-1, 1, 2)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        lc = matplotlib.collections.LineCollection(segments, cmap="plasma", norm=norm)
        lc.set_array(z)
        ax1.add_collection(lc)

    elif plot == "3D":
        points = np.array([x, y, z]).T.reshape(-1, 1, 3)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        for i in range(len(x) - 1):
            seg = segments[i]
            (l,) = ax1.plot(seg[:, 0], seg[:, 1], seg[:, 2], color=colors[i])
            l.set_solid_capstyle("round")
        ax1.set_zlabel("Signal [i-" + str(2 * tau) + "]")

    return fig


# =============================================================================
# Optimal Delay via SPAR Method
# =============================================================================


def _complexity_delay_spar(signal, algorithm=None, show=False, **kwargs):
    if algorithm is None:
        algorithm = "fft"
    # Compute power in freqency domain
    psd = signal_psd(signal, sampling_rate=1000, method=algorithm, show=False, **kwargs)
    power = psd["Power"].values
    freqs = 1000 / psd["Frequency"].values  # Convert to samples

    # Get the 1/3 max frequency (in samples) (https://youtu.be/GGrOJtcTcHA?t=730)
    idx = np.argmax(power)
    optimal = int(freqs[idx] / 3)

    if show is True:
        idxs = freqs <= optimal * 6
        _embedding_delay_plot(
            signal,
            metric_values=power[idxs],
            tau_sequence=freqs[idxs],
            tau=optimal,
            metric="Power",
        )
    return optimal, {"Algorithm": algorithm, "Method": "SPAR"}


# =============================================================================
# Joint-Optimization via Entropy Ratio (Gautama, 2003)
# =============================================================================


def _complexity_delay_entropyratio(
    signal,
    delay_max=20,
    dimension_max=5,
    surrogate_n=5,
    surrogate_method="random",
    show=False,
):
    """Joint-optimization using Entropy Ratio method.

    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.

    Parameters
    ----------
    signal : Union[list, np.array, pd.Series]
        The signal (i.e., a time series) in the form of a vector of values.
    delay_max : int
        The maximum time delay (Tau) to test.
    dimension_max : int
        The maximum embedding dimension (often denoted 'm' or 'd', sometimes referred to as 'order') to test.
    surrogate_iter : int
        The maximum surrogates generated using the iAAFT method.

    """

    # Initalize vectors
    if isinstance(delay_max, int):
        delay_max = np.arange(1, delay_max + 1)
    if isinstance(dimension_max, int):
        dimension_max = np.arange(2, dimension_max + 1)

    # Prepare output
    out = pd.DataFrame(
        {
            "Dimension": np.repeat(dimension_max, len(delay_max)),
            "Delay": np.tile(delay_max, len(dimension_max)),
        }
    )
    # N is the number of delay vectors
    # out["N"] = len(signal) - (out["Dimension"] - 1) * out["Delay"]
    out["N"] = len(signal)

    klens = np.full(len(out), np.nan)
    klens_surrogates = np.full(len(out), np.nan)

    for i in range(len(out)):
        klens[i], _ = entropy_kl(
            signal,
            delay=out["Delay"][i],
            dimension=out["Dimension"][i],
        )
        klens_surrogates[i] = np.nanmean(
            [
                entropy_kl(
                    signal_surrogate(signal, method=surrogate_method),
                    delay=out["Delay"][i],
                    dimension=out["Dimension"][i],
                )[0]
                for j in range(surrogate_n)
            ]
        )

    out["KLEn_Signal"] = klens
    out["KLEn_Surrogate"] = klens_surrogates
    # Entropy Ratio (ER)
    out["KLEn_Ratio"] = klens / klens_surrogates
    # To penalise for higher embedding dimensions, the minimum
    # description length (MDL) method is superimposed, yielding
    # the "entropy ratio" (ER):
    out["Entropy_Ratio"] = out["KLEn_Ratio"] * (
        1 + (out["Dimension"] * np.log(out["N"])) / out["N"]
    )

    # optimal dimension and tau is where entropy_ratio is minimum
    idx = out["Entropy_Ratio"].argmin()
    optimal_dimension = out["Dimension"][idx]
    optimal_delay = out["Delay"][idx]

    if show is True:
        plt.figure()
        ax = plt.axes(projection="3d")
        ax.set_title("Joint-Estimation of Optimal Delay and Dimension")
        ax.plot_trisurf(
            out["Delay"],
            out["Dimension"],
            out["Entropy_Ratio"],
            cmap=plt.get_cmap("GnBu"),
            antialiased=True,
            linewidth=0,
            label="Minimum Entropy Ratio",
            zorder=1,
        )
        ax.plot(
            [optimal_delay] * 2,
            [optimal_dimension] * 2,
            [out["Entropy_Ratio"].min(), out["Entropy_Ratio"].max()],
            zorder=2,
        )
        ax.scatter(
            [optimal_delay] * 2,
            [optimal_dimension] * 2,
            [out["Entropy_Ratio"].min(), out["Entropy_Ratio"].max()],
            color="red",
            zorder=2,
        )
        ax.set_ylabel("Dimension")
        ax.set_xlabel("Delay")
        ax.set_zlabel("Entropy Ratio")
        if len(dimension_max) < 10:
            ax.set_yticks(dimension_max)
        if len(delay_max) < 10:
            ax.set_xticks(delay_max)

    return optimal_delay, {"Dimension": optimal_dimension, "Data": out}


# =============================================================================
# Joint-Optimization via SymbolicDynamic (Matilla-García, 2021)
# =============================================================================


# def _complexity_delay_symbolicdynamics(
#     signal,
#     delay_max=20,
#     dimension_max=5,
#     show=False,
# ):
#     """https://www.mdpi.com/1099-4300/23/2/221/htm"""

#     delay = 1
#     dimension = 3
#     e = 3

#     signal = [2, -7, -12, 5, -1, 9, 14]

#     embedded = nk.complexity_embedding(signal, delay=delay, dimension=dimension)
#     print(embedded)

#     # How to create the symbolic sequence?
#     symbols = np.zeros((len(embedded), dimension - 1))
#     for d in range(1, dimension):
#         # Difference with e
#         symbols[:, d - 1][np.all(np.abs(embedded[:, d - 1 : d + 1]) >= e, axis=1)] = 1

#     symbols[:, 0][np.abs(embedded[:, 0]) >= e] = 1
#     symbols[:, 1][np.all(embedded[:, 1:] < e, axis=1)] = 1

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

import matplotlib
import matplotlib.collections
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.signal
import scipy.spatial
import scipy.stats

from ..misc import NeuroKitWarning, find_closest
from ..signal import (
    signal_autocor,
    signal_findpeaks,
    signal_psd,
    signal_surrogate,
    signal_zerocrossings,
)
from .entropy_kl import entropy_kl
from .information_mutual import mutual_information
from .utils_complexity_embedding import complexity_embedding


def complexity_delay(
    signal, delay_max=50, method="fraser1986", algorithm=None, show=False, **kwargs
):
    """**Automated selection of the optimal Delay (Tau)**

    The time delay (Tau :math:`\\tau`, also referred to as *Lag*) is one of the two critical
    parameters (the other being the :func:`Dimension <complexity_dimension>` *m*) involved in the
    construction of the time-delay embedding of a signal. It corresponds to the delay in samples
    between the original signal and its delayed version(s). In other words, how many samples do we
    consider between a given state of the signal and its closest past state.

    When :math:`\\tau` is smaller than the optimal theoretical value, consecutive coordinates of the
    system's state are correlated and the attractor is not sufficiently unfolded. Conversely, when
    :math:`\\tau` is larger than it should be, successive coordinates are almost independent,
    resulting in an uncorrelated and unstructured cloud of points.

    The selection of the parameters *delay* and :func:`*dimension* <complexity_dimension>` is a
    challenge. One approach is to select them (semi) independently (as dimension selection often
    requires the delay), using :func:`complexity_delay` and :func:`complexity_dimension`. However,
    some joint-estimation methods do exist, that attempt at finding the optimal delay and dimension
    at the same time.

    Several authors suggested different methods to guide the choice of the delay:

    * **Fraser and Swinney (1986)** suggest using the first local minimum of the mutual information
      between the delayed and non-delayed time series, effectively identifying a value of Tau for
      which they share the least information.
    * **Theiler (1990)** suggested to select Tau where the autocorrelation between the signal and
      its lagged version at Tau first crosses the value :math:`1/e`.
    * **Casdagli (1991)** suggests instead taking the first zero-crossing of the autocorrelation.
    * **Rosenstein (1993)** suggests to approximate the point where the autocorrelation function
      drops to :math:`(1 - 1/e)` of its maximum value.
    * **Rosenstein (1994)** suggests to the point close to 40% of the slope of the average
      displacement from the diagonal (ADFD).
    * **Kim (1999)** suggests estimating Tau using the correlation integral, called the C-C method,
      which has shown to agree with those obtained using the Mutual Information. This method
      makes use of a statistic within the reconstructed phase space, rather than analyzing the
      temporal evolution of the time series. However, computation times are significantly long for
      this method due to the need to compare every unique pair of pairwise vectors within the
      embedded signal per delay.
    * **Lyle (2021)** describes the "Symmetric Projection Attractor Reconstruction" (SPAR), where
      :math:`1/3` of the the dominant frequency (i.e., of the length of the average "cycle") can be
      a suitable value for approximately periodic data, and makes the attractor sensitive to
      morphological changes. See also `Aston's talk <https://youtu.be/GGrOJtcTcHA?t=730>`_. This
      method is also the fastest but might not be suitable for aperiodic signals.
      The ``algorithm`` argument (default to ``"fft"``) and will be passed as the ``method``
      argument of  ``signal_psd()``.

    **Joint-Methods for Delay and Dimension**

    * **Gautama (2003)** mentions that in practice, it is common to have a fixed time lag and to
      adjust the embedding dimension accordingly. As this can lead to large *m* values (and thus to
      embedded data of a large size) and thus, slow processing, they describe an optimisation
      method to jointly determine *m* and :math:`\\tau`, based on the **entropy ratio**.

    .. note::

      We would like to implement the joint-estimation by `Matilla-García et al. (2021)
      <https://www.mdpi.com/1099-4300/23/2/221>`_, please get in touch if you can help us!

    Parameters
    ----------
    signal : Union[list, np.array, pd.Series]
        The signal (i.e., a time series) in the form of a vector of values.
    delay_max : int
        The maximum time delay (Tau or lag) to test.
    method : str
        The method that defines what to compute for each tested value of Tau. Can be one of
        ``"fraser1986"``, ``"theiler1990"``, ``"casdagli1991"``, ``"rosenstein1993"``,
        ``"rosenstein1994"``, ``"kim1999"``, or ``"lyle2021"``.
    algorithm : str
        The method used to find the optimal value of Tau given the values computed by the method.
        If ``None`` (default), will select the algorithm according to the method. Modify only if you
        know what you are doing.
    show : bool
        If ``True``, will plot the metric values for each value of tau.
    **kwargs : optional
        Additional arguments to be passed for C-C method.

    Returns
    -------
    delay : int
        Optimal time delay.
    parameters : dict
        A dictionary containing additional information regarding the parameters used
        to compute optimal time-delay embedding.

    See Also
    ---------
    complexity, complexity_dimension, complexity_embedding, complexity_tolerance

    Examples
    ----------

    * **Example 1**: Comparison of different methods for estimating the optimal delay of an simple
      artificial signal.

    .. ipython:: python

      import neurokit2 as nk

      signal = nk.signal_simulate(duration=10, sampling_rate=100, frequency=[1, 1.5], noise=0.02)

      @savefig p_complexity_delay1.png scale=100%
      nk.signal_plot(signal)
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay2.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="fraser1986")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay3.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="theiler1990")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay4.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="casdagli1991")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay5.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="rosenstein1993")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay6.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="rosenstein1994")
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay7.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,
                                              method="lyle2021")
      @suppress
      plt.close()

    * **Example 2**: Optimizing the delay and the dimension using joint-estimation methods.

    .. ipython:: python

      @savefig p_complexity_delay8.png scale=100%
      delay, parameters = nk.complexity_delay(
          signal,
          delay_max=np.arange(1, 30, 1),  # Can be an int or a list
          dimension_max=20,  # Can be an int or a list
          method="gautama2003",
          surrogate_n=5,  # Number of surrogate signals to generate
          surrogate_method="random",  # can be IAAFT, see nk.signal_surrogate()
          show=True)
      @suppress
      plt.close()

    .. ipython:: python

      # Optimal dimension
      dimension = parameters["Dimension"]
      dimension

    **Note**: A double-checking of that method would be appreciated! Please help us improve.

    * **Example 3**: Using a realistic signal.

    .. ipython:: python

      ecg = nk.ecg_simulate(duration=60*6, sampling_rate=200)
      signal = nk.ecg_rate(nk.ecg_peaks(ecg, sampling_rate=200),
                           sampling_rate=200,
                           desired_length=len(ecg))
      @savefig p_complexity_delay9.png scale=100%
      nk.signal_plot(signal)
      @suppress
      plt.close()

    .. ipython:: python

      @savefig p_complexity_delay10.png scale=100%
      delay, parameters = nk.complexity_delay(signal, delay_max=1000, show=True)
      @suppress
      plt.close()

    References
    ------------
    * <NAME>., <NAME>., & <NAME>. (2021). Symmetric Projection Attractor Reconstruction:
      Sex Differences in the ECG. Frontiers in cardiovascular medicine, 1034.
    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.
    * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time
      series forecasting. IFAC Proceedings Volumes, 42(7), 316-320.
    * <NAME>., <NAME>., & <NAME>. (1993). A practical method for calculating
      largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65(1-2),
      117-134.
    * <NAME>., <NAME>., & <NAME>. (1994). Reconstruction expansion as a
      geometry-based framework for choosing proper delay times. Physica-Section D, 73(1), 82-98.
    * <NAME>., <NAME>., & <NAME>. (1999). Nonlinear dynamics, delay times, and embedding
      windows. Physica D: Nonlinear Phenomena, 127(1-2), 48-60.
    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.
    * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time
      series forecasting. IFAC Proceedings Volumes, 42(7), 316-320.

    """
    # Initalize vectors
    if isinstance(delay_max, int):
        tau_sequence = np.arange(1, delay_max + 1)
    else:
        tau_sequence = np.array(delay_max)

    # Method
    method = method.lower()
    if method in ["fraser", "fraser1986", "tdmi"]:
        metric = "Mutual Information"
        if algorithm is None:
            algorithm = "first local minimum"
    elif method in ["theiler", "theiler1990"]:
        metric = "Autocorrelation"
        if algorithm is None:
            algorithm = "first 1/e crossing"
    elif method in ["casdagli", "casdagli1991"]:
        metric = "Autocorrelation"
        if algorithm is None:
            algorithm = "first zero crossing"
    elif method in ["rosenstein", "rosenstein1994", "adfd"]:
        metric = "Displacement"
        if algorithm is None:
            algorithm = "closest to 40% of the slope"
    elif method in ["rosenstein1993"]:
        metric = "Autocorrelation (FFT)"
        if algorithm is None:
            algorithm = "first drop below 1-(1/e) of maximum"
    elif method in ["kim1999", "cc"]:
        metric = "Correlation Integral"
        if algorithm is None:
            algorithm = "first local minimum"
    elif method in ["aston2020", "lyle2021", "spar"]:
        return _complexity_delay_spar(signal, algorithm=algorithm, show=show, **kwargs)
    elif method in ["gautama2003", "gautama", "entropyratio"]:
        return _complexity_delay_entropyratio(signal, delay_max=delay_max, show=show, **kwargs)
    else:
        raise ValueError("NeuroKit error: complexity_delay(): 'method' not recognized.")

    # Get metric
    metric_values = _embedding_delay_metric(signal, tau_sequence, metric=metric)

    # Get optimal tau
    optimal = _embedding_delay_select(metric_values, algorithm=algorithm)
    if np.isnan(optimal):
        warn(
            "No optimal time delay is found. Nan is returned."
            " Consider using a higher `delay_max`.",
            category=NeuroKitWarning,
        )
        return optimal
    optimal = tau_sequence[optimal]

    if show is True:
        _embedding_delay_plot(
            signal,
            metric_values=metric_values,
            tau_sequence=tau_sequence,
            tau=optimal,
            metric=metric,
        )

    # Return optimal tau and info dict
    return optimal, {
        "Values": tau_sequence,
        "Scores": metric_values,
        "Algorithm": algorithm,
        "Metric": metric,
        "Method": method,
    }


# =============================================================================
# Methods
# =============================================================================
def _embedding_delay_select(metric_values, algorithm="first local minimum"):

    if algorithm in ["first local minimum (corrected)"]:
        # if immediately increasing, assume it is the first is the closest
        if np.diff(metric_values)[0] > 0:
            optimal = 0
        # take last value if continuously decreasing with no inflections
        elif all(np.diff(metric_values) < 0):
            optimal = len(metric_values) - 1
        else:
            # Find reversed peaks
            optimal = signal_findpeaks(
                -1 * metric_values, relative_height_min=0.1, relative_max=True
            )["Peaks"]
    elif algorithm == "first local minimum":
        # Find reversed peaks
        try:
            optimal = signal_findpeaks(
                -1 * metric_values, relative_height_min=0.1, relative_max=True
            )["Peaks"]
        except ValueError:
            warn(
                "First local minimum detection failed. Try setting "
                + "`algorithm = 'first local minimum (corrected)'` or using another method.",
                category=NeuroKitWarning,
            )

    elif algorithm == "first 1/e crossing":
        metric_values = metric_values - 1 / np.exp(1)
        optimal = signal_zerocrossings(metric_values)
    elif algorithm == "first zero crossing":
        optimal = signal_zerocrossings(metric_values)
    elif algorithm == "closest to 40% of the slope":
        slope = np.diff(metric_values) * len(metric_values)
        slope_in_deg = np.rad2deg(np.arctan(slope))
        optimal = np.where(slope_in_deg == find_closest(40, slope_in_deg))[0]
    elif algorithm == "first drop below 1-(1/e) of maximum":
        optimal = np.where(metric_values < np.max(metric_values) * (1 - 1.0 / np.e))[0][0]

    if not isinstance(optimal, (int, float, np.integer)):
        if len(optimal) != 0:
            optimal = optimal[0]
        else:
            optimal = np.nan

    return optimal


def _embedding_delay_metric(
    signal,
    tau_sequence,
    metric="Mutual Information",
    dimensions=[2, 3, 4, 5],
    r_vals=[0.5, 1.0, 1.5, 2.0],
):
    """Iterating through dimensions and r values is relevant only if metric used is Correlation Integral.
    For this method, either first zero crossing of the statistic averages or the first local
    minimum of deviations to obtain optimal tau. This implementation takes the latter since in practice,
    they are both in close proximity.
    """

    if metric == "Autocorrelation":
        values, _ = signal_autocor(signal)
        values = values[: len(tau_sequence)]  # upper limit

    elif metric == "Autocorrelation (FFT)":
        values, _ = signal_autocor(signal, demean=False, method="fft")
        values = values[: len(tau_sequence)]

    elif metric == "Correlation Integral":
        r_vals = [i * np.std(signal) for i in r_vals]
        # initiate empty list for storing
        # averages = np.zeros(len(tau_sequence))
        values = np.zeros(len(tau_sequence))
        for i, t in enumerate(tau_sequence):
            # average = 0
            change = 0
            for m in dimensions:
                # # find average of dependence statistic
                # for r in r_vals:
                #     s = _embedding_delay_cc_statistic(signal, delay=t, dimension=m, r=r)
                #     average += s

                # find average of statistic deviations across r_vals
                deviation = _embedding_delay_cc_deviation_max(
                    signal, delay=t, dimension=m, r_vals=r_vals
                )
                change += deviation
            # averages[i] = average / 16
            values[i] = change / 4

    else:
        values = np.zeros(len(tau_sequence))

        # Loop through taus and compute all scores values
        for i, current_tau in enumerate(tau_sequence):
            embedded = complexity_embedding(signal, delay=current_tau, dimension=2)
            if metric == "Mutual Information":
                values[i] = mutual_information(
                    embedded[:, 0], embedded[:, 1], normalized=True, method="shannon"
                )
            if metric == "Displacement":
                dimension = 2

                # Reconstruct with zero time delay.
                tau0 = embedded[:, 0].repeat(dimension).reshape(len(embedded), dimension)
                dist = np.asarray(
                    [scipy.spatial.distance.euclidean(i, j) for i, j in zip(embedded, tau0)]
                )
                values[i] = np.mean(dist)

    return values


# =============================================================================
# Internals for C-C method, Kim et al. (1999)
# =============================================================================


def _embedding_delay_cc_integral_sum(signal, dimension=3, delay=10, r=0.02):
    """Correlation integral is a cumulative distribution function, which denotes
    the probability of distance between any pairs of points in phase space
    not greater than the specified `r`.
    """

    # Embed signal
    embedded = complexity_embedding(signal, delay=delay, dimension=dimension, show=False)
    M = embedded.shape[0]

    # Prepare indices for comparing all unique pairwise vectors
    combinations = list(itertools.combinations(range(0, M), r=2))
    first_index, second_index = np.transpose(combinations)[0], np.transpose(combinations)[1]

    vectorized_integral = np.vectorize(_embedding_delay_cc_integral, excluded=["embedded", "r"])
    integral = np.sum(
        vectorized_integral(
            first_index=first_index, second_index=second_index, embedded=embedded, r=r
        )
    )

    return integral


def _embedding_delay_cc_integral(first_index, second_index, embedded, r=0.02):
    M = embedded.shape[0]  # Number of embedded points
    diff = np.linalg.norm(embedded[first_index] - embedded[second_index], ord=np.inf)  # sup-norm
    h = np.heaviside(r - diff, 1)
    integral = (2 / (M * (M - 1))) * h  # find average

    return integral


def _embedding_delay_cc_statistic(signal, dimension=3, delay=10, r=0.02):
    """The dependence statistic as the serial correlation of a nonlinear time series."""

    # create disjoint time series
    series = [signal[i - 1 :: delay] for i in range(1, delay + 1)]

    statistic = 0
    for sub_series in series:
        diff = _embedding_delay_cc_integral_sum(
            sub_series, dimension=dimension, delay=delay, r=r
        ) - ((_embedding_delay_cc_integral_sum(signal, dimension=1, delay=delay, r=r)) ** dimension)
        statistic += diff

    return statistic / delay


def _embedding_delay_cc_deviation_max(signal, r_vals=[0.5, 1.0, 1.5, 2.0], delay=10, dimension=3):
    """A measure of the variation of the dependence statistic with r using
    several representative values of r.
    """
    vectorized_deviation = np.vectorize(
        _embedding_delay_cc_deviation, excluded=["signal", "delay", "dimension"]
    )
    deviations = vectorized_deviation(
        signal=signal, r_vals=r_vals, delay=delay, dimension=dimension
    )

    return np.max(deviations) - np.min(deviations)


def _embedding_delay_cc_deviation(signal, r_vals=[0.5, 1.0, 1.5, 2.0], delay=10, dimension=3):
    return _embedding_delay_cc_statistic(signal, delay=delay, dimension=dimension, r=r_vals)


# =============================================================================
# Plotting Generics
# =============================================================================
def _embedding_delay_plot(
    signal,
    metric_values,
    tau_sequence,
    tau=1,
    metric="Mutual Information",
    ax0=None,
    ax1=None,
    plot="2D",
):

    # Prepare figure
    if ax0 is None and ax1 is None:
        fig = plt.figure(constrained_layout=False)
        spec = matplotlib.gridspec.GridSpec(
            ncols=1, nrows=2, height_ratios=[1, 3], width_ratios=[2]
        )
        ax0 = fig.add_subplot(spec[0])
        if plot == "2D":
            ax1 = fig.add_subplot(spec[1])
        elif plot == "3D":
            ax1 = fig.add_subplot(spec[1], projection="3d")
    else:
        fig = None

    ax0.set_title("Optimization of Delay")
    ax0.set_xlabel("Time Delay")
    ax0.set_ylabel(metric)
    ax0.plot(tau_sequence, metric_values, color="#FFC107")
    ax0.axvline(x=tau, color="#E91E63", label="Optimal delay: " + str(tau))
    ax0.legend(loc="upper right")
    ax1.set_title("Attractor")
    ax1.set_xlabel("Signal [i]")
    ax1.set_ylabel("Signal [i-" + str(tau) + "]")

    # Get data points, set axis limits
    embedded = complexity_embedding(signal, delay=tau, dimension=3)
    x = embedded[:, 0]
    y = embedded[:, 1]
    z = embedded[:, 2]
    ax1.set_xlim(x.min(), x.max())
    ax1.set_ylim(x.min(), x.max())

    # Colors
    norm = plt.Normalize(z.min(), z.max())
    cmap = plt.get_cmap("plasma")
    colors = cmap(norm(x))

    # Attractor for 2D vs 3D
    if plot == "2D":
        points = np.array([x, y]).T.reshape(-1, 1, 2)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        lc = matplotlib.collections.LineCollection(segments, cmap="plasma", norm=norm)
        lc.set_array(z)
        ax1.add_collection(lc)

    elif plot == "3D":
        points = np.array([x, y, z]).T.reshape(-1, 1, 3)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        for i in range(len(x) - 1):
            seg = segments[i]
            (l,) = ax1.plot(seg[:, 0], seg[:, 1], seg[:, 2], color=colors[i])
            l.set_solid_capstyle("round")
        ax1.set_zlabel("Signal [i-" + str(2 * tau) + "]")

    return fig


# =============================================================================
# Optimal Delay via SPAR Method
# =============================================================================


def _complexity_delay_spar(signal, algorithm=None, show=False, **kwargs):
    if algorithm is None:
        algorithm = "fft"
    # Compute power in freqency domain
    psd = signal_psd(signal, sampling_rate=1000, method=algorithm, show=False, **kwargs)
    power = psd["Power"].values
    freqs = 1000 / psd["Frequency"].values  # Convert to samples

    # Get the 1/3 max frequency (in samples) (https://youtu.be/GGrOJtcTcHA?t=730)
    idx = np.argmax(power)
    optimal = int(freqs[idx] / 3)

    if show is True:
        idxs = freqs <= optimal * 6
        _embedding_delay_plot(
            signal,
            metric_values=power[idxs],
            tau_sequence=freqs[idxs],
            tau=optimal,
            metric="Power",
        )
    return optimal, {"Algorithm": algorithm, "Method": "SPAR"}


# =============================================================================
# Joint-Optimization via Entropy Ratio (Gautama, 2003)
# =============================================================================


def _complexity_delay_entropyratio(
    signal,
    delay_max=20,
    dimension_max=5,
    surrogate_n=5,
    surrogate_method="random",
    show=False,
):
    """Joint-optimization using Entropy Ratio method.

    * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based
      method for determining the optimal embedding parameters of a signal. In 2003 IEEE
      International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.
      (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.

    Parameters
    ----------
    signal : Union[list, np.array, pd.Series]
        The signal (i.e., a time series) in the form of a vector of values.
    delay_max : int
        The maximum time delay (Tau) to test.
    dimension_max : int
        The maximum embedding dimension (often denoted 'm' or 'd', sometimes referred to as 'order') to test.
    surrogate_iter : int
        The maximum surrogates generated using the iAAFT method.

    """

    # Initalize vectors
    if isinstance(delay_max, int):
        delay_max = np.arange(1, delay_max + 1)
    if isinstance(dimension_max, int):
        dimension_max = np.arange(2, dimension_max + 1)

    # Prepare output
    out = pd.DataFrame(
        {
            "Dimension": np.repeat(dimension_max, len(delay_max)),
            "Delay": np.tile(delay_max, len(dimension_max)),
        }
    )
    # N is the number of delay vectors
    # out["N"] = len(signal) - (out["Dimension"] - 1) * out["Delay"]
    out["N"] = len(signal)

    klens = np.full(len(out), np.nan)
    klens_surrogates = np.full(len(out), np.nan)

    for i in range(len(out)):
        klens[i], _ = entropy_kl(
            signal,
            delay=out["Delay"][i],
            dimension=out["Dimension"][i],
        )
        klens_surrogates[i] = np.nanmean(
            [
                entropy_kl(
                    signal_surrogate(signal, method=surrogate_method),
                    delay=out["Delay"][i],
                    dimension=out["Dimension"][i],
                )[0]
                for j in range(surrogate_n)
            ]
        )

    out["KLEn_Signal"] = klens
    out["KLEn_Surrogate"] = klens_surrogates
    # Entropy Ratio (ER)
    out["KLEn_Ratio"] = klens / klens_surrogates
    # To penalise for higher embedding dimensions, the minimum
    # description length (MDL) method is superimposed, yielding
    # the "entropy ratio" (ER):
    out["Entropy_Ratio"] = out["KLEn_Ratio"] * (
        1 + (out["Dimension"] * np.log(out["N"])) / out["N"]
    )

    # optimal dimension and tau is where entropy_ratio is minimum
    idx = out["Entropy_Ratio"].argmin()
    optimal_dimension = out["Dimension"][idx]
    optimal_delay = out["Delay"][idx]

    if show is True:
        plt.figure()
        ax = plt.axes(projection="3d")
        ax.set_title("Joint-Estimation of Optimal Delay and Dimension")
        ax.plot_trisurf(
            out["Delay"],
            out["Dimension"],
            out["Entropy_Ratio"],
            cmap=plt.get_cmap("GnBu"),
            antialiased=True,
            linewidth=0,
            label="Minimum Entropy Ratio",
            zorder=1,
        )
        ax.plot(
            [optimal_delay] * 2,
            [optimal_dimension] * 2,
            [out["Entropy_Ratio"].min(), out["Entropy_Ratio"].max()],
            zorder=2,
        )
        ax.scatter(
            [optimal_delay] * 2,
            [optimal_dimension] * 2,
            [out["Entropy_Ratio"].min(), out["Entropy_Ratio"].max()],
            color="red",
            zorder=2,
        )
        ax.set_ylabel("Dimension")
        ax.set_xlabel("Delay")
        ax.set_zlabel("Entropy Ratio")
        if len(dimension_max) < 10:
            ax.set_yticks(dimension_max)
        if len(delay_max) < 10:
            ax.set_xticks(delay_max)

    return optimal_delay, {"Dimension": optimal_dimension, "Data": out}


# =============================================================================
# Joint-Optimization via SymbolicDynamic (Matilla-García, 2021)
# =============================================================================


# def _complexity_delay_symbolicdynamics(
#     signal,
#     delay_max=20,
#     dimension_max=5,
#     show=False,
# ):
#     """https://www.mdpi.com/1099-4300/23/2/221/htm"""

#     delay = 1
#     dimension = 3
#     e = 3

#     signal = [2, -7, -12, 5, -1, 9, 14]

#     embedded = nk.complexity_embedding(signal, delay=delay, dimension=dimension)
#     print(embedded)

#     # How to create the symbolic sequence?
#     symbols = np.zeros((len(embedded), dimension - 1))
#     for d in range(1, dimension):
#         # Difference with e
#         symbols[:, d - 1][np.all(np.abs(embedded[:, d - 1 : d + 1]) >= e, axis=1)] = 1

#     symbols[:, 0][np.abs(embedded[:, 0]) >= e] = 1
#     symbols[:, 1][np.all(embedded[:, 1:] < e, axis=1)] = 1
# -*- coding: utf-8 -*- **Automated selection of the optimal Delay (Tau)**      The time delay (Tau :math:`\\tau`, also referred to as *Lag*) is one of the two critical     parameters (the other being the :func:`Dimension <complexity_dimension>` *m*) involved in the     construction of the time-delay embedding of a signal. It corresponds to the delay in samples     between the original signal and its delayed version(s). In other words, how many samples do we     consider between a given state of the signal and its closest past state.      When :math:`\\tau` is smaller than the optimal theoretical value, consecutive coordinates of the     system's state are correlated and the attractor is not sufficiently unfolded. Conversely, when     :math:`\\tau` is larger than it should be, successive coordinates are almost independent,     resulting in an uncorrelated and unstructured cloud of points.      The selection of the parameters *delay* and :func:`*dimension* <complexity_dimension>` is a     challenge. One approach is to select them (semi) independently (as dimension selection often     requires the delay), using :func:`complexity_delay` and :func:`complexity_dimension`. However,     some joint-estimation methods do exist, that attempt at finding the optimal delay and dimension     at the same time.      Several authors suggested different methods to guide the choice of the delay:      * **Fraser and Swinney (1986)** suggest using the first local minimum of the mutual information       between the delayed and non-delayed time series, effectively identifying a value of Tau for       which they share the least information.     * **Theiler (1990)** suggested to select Tau where the autocorrelation between the signal and       its lagged version at Tau first crosses the value :math:`1/e`.     * **Casdagli (1991)** suggests instead taking the first zero-crossing of the autocorrelation.     * **Rosenstein (1993)** suggests to approximate the point where the autocorrelation function       drops to :math:`(1 - 1/e)` of its maximum value.     * **Rosenstein (1994)** suggests to the point close to 40% of the slope of the average       displacement from the diagonal (ADFD).     * **Kim (1999)** suggests estimating Tau using the correlation integral, called the C-C method,       which has shown to agree with those obtained using the Mutual Information. This method       makes use of a statistic within the reconstructed phase space, rather than analyzing the       temporal evolution of the time series. However, computation times are significantly long for       this method due to the need to compare every unique pair of pairwise vectors within the       embedded signal per delay.     * **Lyle (2021)** describes the "Symmetric Projection Attractor Reconstruction" (SPAR), where       :math:`1/3` of the the dominant frequency (i.e., of the length of the average "cycle") can be       a suitable value for approximately periodic data, and makes the attractor sensitive to       morphological changes. See also `Aston's talk <https://youtu.be/GGrOJtcTcHA?t=730>`_. This       method is also the fastest but might not be suitable for aperiodic signals.       The ``algorithm`` argument (default to ``"fft"``) and will be passed as the ``method``       argument of  ``signal_psd()``.      **Joint-Methods for Delay and Dimension**      * **Gautama (2003)** mentions that in practice, it is common to have a fixed time lag and to       adjust the embedding dimension accordingly. As this can lead to large *m* values (and thus to       embedded data of a large size) and thus, slow processing, they describe an optimisation       method to jointly determine *m* and :math:`\\tau`, based on the **entropy ratio**.      .. note::        We would like to implement the joint-estimation by `Matilla-García et al. (2021)       <https://www.mdpi.com/1099-4300/23/2/221>`_, please get in touch if you can help us!      Parameters     ----------     signal : Union[list, np.array, pd.Series]         The signal (i.e., a time series) in the form of a vector of values.     delay_max : int         The maximum time delay (Tau or lag) to test.     method : str         The method that defines what to compute for each tested value of Tau. Can be one of         ``"fraser1986"``, ``"theiler1990"``, ``"casdagli1991"``, ``"rosenstein1993"``,         ``"rosenstein1994"``, ``"kim1999"``, or ``"lyle2021"``.     algorithm : str         The method used to find the optimal value of Tau given the values computed by the method.         If ``None`` (default), will select the algorithm according to the method. Modify only if you         know what you are doing.     show : bool         If ``True``, will plot the metric values for each value of tau.     **kwargs : optional         Additional arguments to be passed for C-C method.      Returns     -------     delay : int         Optimal time delay.     parameters : dict         A dictionary containing additional information regarding the parameters used         to compute optimal time-delay embedding.      See Also     ---------     complexity, complexity_dimension, complexity_embedding, complexity_tolerance      Examples     ----------      * **Example 1**: Comparison of different methods for estimating the optimal delay of an simple       artificial signal.      .. ipython:: python        import neurokit2 as nk        signal = nk.signal_simulate(duration=10, sampling_rate=100, frequency=[1, 1.5], noise=0.02)        @savefig p_complexity_delay1.png scale=100%       nk.signal_plot(signal)       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay2.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="fraser1986")       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay3.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="theiler1990")       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay4.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="casdagli1991")       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay5.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="rosenstein1993")       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay6.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="rosenstein1994")       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay7.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=100, show=True,                                               method="lyle2021")       @suppress       plt.close()      * **Example 2**: Optimizing the delay and the dimension using joint-estimation methods.      .. ipython:: python        @savefig p_complexity_delay8.png scale=100%       delay, parameters = nk.complexity_delay(           signal,           delay_max=np.arange(1, 30, 1),  # Can be an int or a list           dimension_max=20,  # Can be an int or a list           method="gautama2003",           surrogate_n=5,  # Number of surrogate signals to generate           surrogate_method="random",  # can be IAAFT, see nk.signal_surrogate()           show=True)       @suppress       plt.close()      .. ipython:: python        # Optimal dimension       dimension = parameters["Dimension"]       dimension      **Note**: A double-checking of that method would be appreciated! Please help us improve.      * **Example 3**: Using a realistic signal.      .. ipython:: python        ecg = nk.ecg_simulate(duration=60*6, sampling_rate=200)       signal = nk.ecg_rate(nk.ecg_peaks(ecg, sampling_rate=200),                            sampling_rate=200,                            desired_length=len(ecg))       @savefig p_complexity_delay9.png scale=100%       nk.signal_plot(signal)       @suppress       plt.close()      .. ipython:: python        @savefig p_complexity_delay10.png scale=100%       delay, parameters = nk.complexity_delay(signal, delay_max=1000, show=True)       @suppress       plt.close()      References     ------------     * <NAME>., <NAME>., & <NAME>. (2021). Symmetric Projection Attractor Reconstruction:       Sex Differences in the ECG. Frontiers in cardiovascular medicine, 1034.     * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based       method for determining the optimal embedding parameters of a signal. In 2003 IEEE       International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.       (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.     * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time       series forecasting. IFAC Proceedings Volumes, 42(7), 316-320.     * <NAME>., <NAME>., & <NAME>. (1993). A practical method for calculating       largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65(1-2),       117-134.     * <NAME>., <NAME>., & <NAME>. (1994). Reconstruction expansion as a       geometry-based framework for choosing proper delay times. Physica-Section D, 73(1), 82-98.     * <NAME>., <NAME>., & <NAME>. (1999). Nonlinear dynamics, delay times, and embedding       windows. Physica D: Nonlinear Phenomena, 127(1-2), 48-60.     * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based       method for determining the optimal embedding parameters of a signal. In 2003 IEEE       International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.       (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.     * <NAME>., & <NAME>. (2009). The role of the embedding dimension and time delay in time       series forecasting. IFAC Proceedings Volumes, 42(7), 316-320. # Initalize vectors # Method # Get metric # Get optimal tau # Return optimal tau and info dict # ============================================================================= # Methods # ============================================================================= # if immediately increasing, assume it is the first is the closest # take last value if continuously decreasing with no inflections # Find reversed peaks # Find reversed peaks Iterating through dimensions and r values is relevant only if metric used is Correlation Integral.     For this method, either first zero crossing of the statistic averages or the first local     minimum of deviations to obtain optimal tau. This implementation takes the latter since in practice,     they are both in close proximity. # upper limit # initiate empty list for storing # averages = np.zeros(len(tau_sequence)) # average = 0 # # find average of dependence statistic # for r in r_vals: #     s = _embedding_delay_cc_statistic(signal, delay=t, dimension=m, r=r) #     average += s # find average of statistic deviations across r_vals # averages[i] = average / 16 # Loop through taus and compute all scores values # Reconstruct with zero time delay. # ============================================================================= # Internals for C-C method, Kim et al. (1999) # ============================================================================= Correlation integral is a cumulative distribution function, which denotes     the probability of distance between any pairs of points in phase space     not greater than the specified `r`. # Embed signal # Prepare indices for comparing all unique pairwise vectors # Number of embedded points # sup-norm # find average The dependence statistic as the serial correlation of a nonlinear time series. # create disjoint time series A measure of the variation of the dependence statistic with r using     several representative values of r. # ============================================================================= # Plotting Generics # ============================================================================= # Prepare figure # Get data points, set axis limits # Colors # Attractor for 2D vs 3D # ============================================================================= # Optimal Delay via SPAR Method # ============================================================================= # Compute power in freqency domain # Convert to samples # Get the 1/3 max frequency (in samples) (https://youtu.be/GGrOJtcTcHA?t=730) # ============================================================================= # Joint-Optimization via Entropy Ratio (Gautama, 2003) # ============================================================================= Joint-optimization using Entropy Ratio method.      * <NAME>., <NAME>., & <NAME>. (2003, April). A differential entropy based       method for determining the optimal embedding parameters of a signal. In 2003 IEEE       International Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings.       (ICASSP'03). (Vol. 6, pp. VI-29). IEEE.      Parameters     ----------     signal : Union[list, np.array, pd.Series]         The signal (i.e., a time series) in the form of a vector of values.     delay_max : int         The maximum time delay (Tau) to test.     dimension_max : int         The maximum embedding dimension (often denoted 'm' or 'd', sometimes referred to as 'order') to test.     surrogate_iter : int         The maximum surrogates generated using the iAAFT method. # Initalize vectors # Prepare output # N is the number of delay vectors # out["N"] = len(signal) - (out["Dimension"] - 1) * out["Delay"] # Entropy Ratio (ER) # To penalise for higher embedding dimensions, the minimum # description length (MDL) method is superimposed, yielding # the "entropy ratio" (ER): # optimal dimension and tau is where entropy_ratio is minimum # ============================================================================= # Joint-Optimization via SymbolicDynamic (Matilla-García, 2021) # ============================================================================= # def _complexity_delay_symbolicdynamics( #     signal, #     delay_max=20, #     dimension_max=5, #     show=False, # ): #     """https://www.mdpi.com/1099-4300/23/2/221/htm""" #     delay = 1 #     dimension = 3 #     e = 3 #     signal = [2, -7, -12, 5, -1, 9, 14] #     embedded = nk.complexity_embedding(signal, delay=delay, dimension=dimension) #     print(embedded) #     # How to create the symbolic sequence? #     symbols = np.zeros((len(embedded), dimension - 1)) #     for d in range(1, dimension): #         # Difference with e #         symbols[:, d - 1][np.all(np.abs(embedded[:, d - 1 : d + 1]) >= e, axis=1)] = 1 #     symbols[:, 0][np.abs(embedded[:, 0]) >= e] = 1 #     symbols[:, 1][np.all(embedded[:, 1:] < e, axis=1)] = 1