code: from warnings import warn
  
  import matplotlib.pyplot as plt
  import numpy as np
  import pandas as pd
  
  from ..misc import NeuroKitWarning, find_plateau
  
  
  def complexity_k(signal, k_max="max", show=False):
      """Automated selection of the optimal k_max parameter for Higuchi Fractal Dimension (HFD).
  
      The optimal kmax is computed based on the point at which HFD values plateau for a range of kmax 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
      ----------
      >>> import neurokit2 as nk
      >>>
      >>> signal = nk.signal_simulate(duration=2, sampling_rate=100, frequency=[5, 6], noise=0.5)
      >>> k_max, info = nk.complexity_k(signal, k_max='default', show=True)
      >>> k_max #doctest: +SKIP
  
      Reference
      ----------
      - 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://ieeexplore.ieee.org/document/7168285
      """
      # 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 = plt.cm.PuBu(np.linspace(0, 1, len(k_range)))
  
      # if single time series
      ax.plot(k_range, slope_values, color="#2196F3", zorder=1)
      for i, j 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
   # -*- coding: utf-8 -*-
  import numpy as np
  import pandas as pd
  
  
  def rsp_findpeaks(rsp_cleaned, sampling_rate=1000, method="khodadad2018", amplitude_min=0.3):
      """Extract extrema in a respiration (RSP) signal.
  
      Low-level function used by `rsp_peaks()` to identify inhalation and exhalation onsets (troughs and peaks
      respectively) in a preprocessed respiration signal using different sets of parameters.
      See `rsp_peaks()` for details.
  
      Parameters
      ----------
      rsp_cleaned : Union[list, np.array, pd.Series]
          The cleaned respiration channel as returned by `rsp_clean()`.
      sampling_rate : int
          The sampling frequency of 'rsp_cleaned' (in Hz, i.e., samples/second).
      method : str
          The processing pipeline to apply. Can be one of "khodadad2018" (default) or "biosppy".
      amplitude_min : float
          Only applies if method is "khodadad2018". Extrema that have a vertical distance smaller than
          (outlier_threshold * average vertical distance) to any direct neighbour are removed as false
          positive outliers. I.e., outlier_threshold should be a float with positive sign (the default is
          0.3). Larger values of outlier_threshold correspond to more conservative thresholds (i.e.,
          more extrema removed as outliers).
  
      Returns
      -------
      info : dict
          A dictionary containing additional information, in this case the samples at which inhalation
          onsets and exhalation onsets occur, accessible with the keys "RSP_Troughs" and "RSP_Peaks", respectively.
  
      See Also
      --------
      rsp_clean, rsp_fixpeaks, rsp_peaks, signal_rate, rsp_amplitude,
      rsp_process, rsp_plot
  
      Examples
      --------
      >>> import neurokit2 as nk
      >>>
      >>> rsp = nk.rsp_simulate(duration=30, respiratory_rate=15)
      >>> cleaned = nk.rsp_clean(rsp, sampling_rate=1000)
      >>> info = nk.rsp_findpeaks(cleaned)
      >>> fig = nk.events_plot([info["RSP_Peaks"], info["RSP_Troughs"]], cleaned)
      >>> fig #doctest: +SKIP
  
      """
      # Try retrieving correct column
      if isinstance(rsp_cleaned, pd.DataFrame):
          try:
              rsp_cleaned = rsp_cleaned["RSP_Clean"]
          except NameError:
              try:
                  rsp_cleaned = rsp_cleaned["RSP_Raw"]
              except NameError:
                  rsp_cleaned = rsp_cleaned["RSP"]
  
      cleaned = np.array(rsp_cleaned)
  
      # Find peaks
      method = method.lower()  # remove capitalised letters
      if method in ["khodadad", "khodadad2018"]:
          info = _rsp_findpeaks_khodadad(cleaned, amplitude_min)
      elif method == "biosppy":
          info = _rsp_findpeaks_biosppy(cleaned, sampling_rate=sampling_rate)
      else:
          raise ValueError("NeuroKit error: rsp_findpeaks(): 'method' should be one of 'khodadad2018' or 'biosppy'.")
  
      return info
  
  
  # =============================================================================
  # Methods
  # =============================================================================
  def _rsp_findpeaks_biosppy(rsp_cleaned, sampling_rate):
  
      extrema = _rsp_findpeaks_extrema(rsp_cleaned)
      extrema, amplitudes = _rsp_findpeaks_outliers(rsp_cleaned, extrema, amplitude_min=0)
  
      peaks, troughs = _rsp_findpeaks_sanitize(extrema, amplitudes)
  
      # Apply minimum period outlier-criterion (exclude inter-breath-intervals
      # that produce breathing rate larger than 35 breaths per minute.
      outlier_idcs = np.where((np.diff(peaks) / sampling_rate) < 1.7)[0]
  
      peaks = np.delete(peaks, outlier_idcs)
      troughs = np.delete(troughs, outlier_idcs)
  
      info = {"RSP_Peaks": peaks, "RSP_Troughs": troughs}
      return info
  
  
  def _rsp_findpeaks_khodadad(rsp_cleaned, amplitude_min=0.3):
  
      extrema = _rsp_findpeaks_extrema(rsp_cleaned)
      extrema, amplitudes = _rsp_findpeaks_outliers(rsp_cleaned, extrema, amplitude_min=amplitude_min)
      peaks, troughs = _rsp_findpeaks_sanitize(extrema, amplitudes)
  
      info = {"RSP_Peaks": peaks, "RSP_Troughs": troughs}
      return info
  
  
  # =============================================================================
  # Internals
  # =============================================================================
  def _rsp_findpeaks_extrema(rsp_cleaned):
      # Detect zero crossings (note that these are zero crossings in the raw
      # signal, not in its gradient).
      greater = rsp_cleaned > 0
      smaller = rsp_cleaned < 0
      risex = np.where(np.bitwise_and(smaller[:-1], greater[1:]))[0]
      fallx = np.where(np.bitwise_and(greater[:-1], smaller[1:]))[0]
  
      if risex[0] < fallx[0]:
          startx = "rise"
      elif fallx[0] < risex[0]:
          startx = "fall"
  
      allx = np.concatenate((risex, fallx))
      allx.sort(kind="mergesort")
  
      # Find extrema by searching minima between falling zero crossing and
      # rising zero crossing, and searching maxima between rising zero
      # crossing and falling zero crossing.
      extrema = []
      for i in range(len(allx) - 1):
  
          # Determine whether to search for minimum or maximum.
          if startx == "rise":
              if (i + 1) % 2 != 0:
                  argextreme = np.argmax
              else:
                  argextreme = np.argmin
          elif startx == "fall":
              if (i + 1) % 2 != 0:
                  argextreme = np.argmin
              else:
                  argextreme = np.argmax
  
          # Get the two zero crossings between which the extreme will be
          # searched.
          beg = allx[i]
          end = allx[i + 1]
  
          extreme = argextreme(rsp_cleaned[beg:end])
          extrema.append(beg + extreme)
  
      extrema = np.asarray(extrema)
      return extrema
  
  
  def _rsp_findpeaks_outliers(rsp_cleaned, extrema, amplitude_min=0.3):
  
      # Only consider those extrema that have a minimum vertical distance to
      # their direct neighbor, i.e., define outliers in absolute amplitude
      # difference between neighboring extrema.
      vertical_diff = np.abs(np.diff(rsp_cleaned[extrema]))
      median_diff = np.median(vertical_diff)
      min_diff = np.where(vertical_diff > (median_diff * amplitude_min))[0]
      extrema = extrema[min_diff]
  
      # Make sure that the alternation of peaks and troughs is unbroken. If
      # alternation of sign in extdiffs is broken, remove the extrema that
      # cause the breaks.
      amplitudes = rsp_cleaned[extrema]
      extdiffs = np.sign(np.diff(amplitudes))
      extdiffs = np.add(extdiffs[0:-1], extdiffs[1:])
      removeext = np.where(extdiffs != 0)[0] + 1
      extrema = np.delete(extrema, removeext)
      amplitudes = np.delete(amplitudes, removeext)
  
      return extrema, amplitudes
  
  
  def _rsp_findpeaks_sanitize(extrema, amplitudes):
      # To be able to consistently calculate breathing amplitude, make sure that
      # the extrema always start with a trough and end with a peak, since
      # breathing amplitude will be defined as vertical distance between each
      # peak and the preceding trough. Note that this also ensures that the
      # number of peaks and troughs is equal.
      if amplitudes[0] > amplitudes[1]:
          extrema = np.delete(extrema, 0)
      if amplitudes[-1] < amplitudes[-2]:
          extrema = np.delete(extrema, -1)
      peaks = extrema[1::2]
      troughs = extrema[0:-1:2]
  
      return peaks, troughs
  
apis: numpy.sum numpy.vectorize numpy.diff numpy.floor matplotlib.pyplot.subplots numpy.arange numpy.max numpy.log numpy.array numpy.linspace
apis: numpy.diff numpy.asarray numpy.median numpy.where numpy.add numpy.delete numpy.array numpy.concatenate numpy.bitwise_and