///|
/// Non-linear HRV descriptors for rhythm complexity and recurrence.
pub(all) struct NonlinearMetrics {
  sample_entropy : Double
  approximate_entropy : Double
  dfa_alpha : Double
  turning_point_ratio : Double
  recurrence_rate : Double
  histogram_entropy : Double
  lag1_autocorrelation : Double
  complexity_index : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Return the absolute distance between two values.
fn absolute_distance(left : Double, right : Double) -> Double {
  let delta = left - right
  if delta < 0.0 {
    -delta
  } else {
    delta
  }
}

///|
/// Count matching embedded vectors of a given dimension.
fn template_match_count(
  values : Array[Double],
  dimension : Int,
  tolerance : Double,
  include_self : Bool,
) -> Int {
  let limit = values.length() - dimension + 1
  if dimension <= 0 || limit <= 0 {
    return 0
  }
  let mut count = 0
  for i in 0.. tolerance {
            matches = false
          }
        }
        if matches {
          count += 1
        }
      }
    }
  }
  count
}

///|
/// Calculate sample entropy with Chebyshev template distance.
pub fn calculate_sample_entropy(
  values : Array[Double],
  dimension : Int,
  tolerance : Double,
) -> Double {
  if values.length() <= dimension + 1 || dimension <= 0 || tolerance < 0.0 {
    return 0.0
  }
  let count_m = template_match_count(values, dimension, tolerance, false)
  let count_m1 = template_match_count(values, dimension + 1, tolerance, false)
  if count_m == 0 || count_m1 == 0 {
    0.0
  } else {
    -@math.ln(count_m1.to_double() / count_m.to_double())
  }
}

///|
/// Calculate approximate entropy with self-matches included.
pub fn calculate_approximate_entropy(
  values : Array[Double],
  dimension : Int,
  tolerance : Double,
) -> Double {
  if values.length() <= dimension + 1 || dimension <= 0 || tolerance < 0.0 {
    return 0.0
  }
  let mut phi_m = 0.0
  let mut phi_m1 = 0.0
  let limit_m = values.length() - dimension + 1
  let limit_m1 = values.length() - dimension
  for i in 0.. tolerance {
          equal = false
        }
      }
      if equal {
        matches += 1
      }
    }
    if matches > 0 {
      phi_m += @math.ln(matches.to_double() / limit_m.to_double())
    }
  }
  for i in 0.. tolerance {
          equal = false
        }
      }
      if equal {
        matches += 1
      }
    }
    if matches > 0 {
      phi_m1 += @math.ln(matches.to_double() / limit_m1.to_double())
    }
  }
  phi_m / limit_m.to_double() - phi_m1 / limit_m1.to_double()
}

///|
/// Calculate the proportion of samples that are local turning points.
pub fn turning_point_ratio(values : Array[Double]) -> Double {
  if values.length() < 3 {
    return 0.0
  }
  let mut turns = 0
  for i in 1..<(values.length() - 1) {
    let rising = values[i] > values[i - 1] && values[i] > values[i + 1]
    let falling = values[i] < values[i - 1] && values[i] < values[i + 1]
    if rising || falling {
      turns += 1
    }
  }
  turns.to_double() / (values.length() - 2).to_double()
}

///|
/// Calculate recurrence rate using a fixed radius around the median.
pub fn recurrence_rate(
  values : Array[Double],
  radius : Double,
  minimum_lag : Int,
) -> Double {
  let n = values.length()
  if n <= 1 || radius < 0.0 {
    return 0.0
  }
  let lag = if minimum_lag < 1 { 1 } else { minimum_lag }
  let mut recurrent = 0
  let mut pairs = 0
  for i in 0.. Double {
  if values.length() == 0 || bin_count <= 0 {
    return 0.0
  }
  let summary = summarize_distribution(values)
  let width = summary.maximum - summary.minimum
  let bins = Array::make(bin_count, 0)
  if width == 0.0 {
    bins[0] = values.length()
  } else {
    for value in values {
      let index = ((value - summary.minimum) / width * bin_count.to_double())
        .floor()
        .to_int()
      let bounded = if index < 0 {
        0
      } else if index >= bin_count {
        bin_count - 1
      } else {
        index
      }
      bins[bounded] += 1
    }
  }
  let mut entropy = 0.0
  for count in bins {
    if count > 0 {
      let p = count.to_double() / values.length().to_double()
      entropy -= p * @math.ln(p)
    }
  }
  if bin_count <= 1 {
    0.0
  } else {
    entropy / @math.ln(bin_count.to_double())
  }
}

///|
/// Build a cumulative sum of deviations from the mean for DFA.
fn cumulative_deviation(values : Array[Double]) -> Array[Double] {
  let result = []
  let average = mean_value(values)
  let mut current = 0.0
  for value in values {
    current += value - average
    result.push(current)
  }
  result
}

///|
/// Calculate a compact detrended fluctuation analysis exponent.
pub fn calculate_dfa_alpha(
  values : Array[Double],
  minimum_scale : Int,
  maximum_scale : Int,
) -> Double {
  if values.length() < 8 {
    return 0.0
  }
  let profile = cumulative_deviation(values)
  let low = if minimum_scale < 2 { 2 } else { minimum_scale }
  let high = if maximum_scale >= values.length() / 2 {
    values.length() / 2
  } else {
    maximum_scale
  }
  let log_scales = []
  let log_fluctuations = []
  let mut scale = low
  while scale <= high {
    let windows = values.length() / scale
    let fluctuations = []
    for w in 0.. 0.0 {
      log_scales.push(@math.ln(scale.to_double()))
      log_fluctuations.push(@math.ln(fluctuation))
    }
    scale += if scale < 16 { 2 } else { scale / 4 }
  }
  if log_scales.length() <= 1 {
    0.0
  } else {
    fit_linear_trend(log_fluctuations).slope
  }
}

///|
/// Estimate sample entropy tolerance from a standard deviation multiplier.
pub fn entropy_tolerance(values : Array[Double], multiplier : Double) -> Double {
  standard_deviation(values) * (if multiplier < 0.0 { 0.0 } else { multiplier })
}

///|
/// Calculate all non-linear metrics with robust default parameters.
pub fn calculate_nonlinear_metrics(values : Array[Double]) -> NonlinearMetrics {
  let tolerance = entropy_tolerance(values, 0.2)
  let sample = calculate_sample_entropy(values, 2, tolerance)
  let approximate = calculate_approximate_entropy(values, 2, tolerance)
  let radius = median_absolute_deviation(values) * 1.5
  let recurrence = recurrence_rate(values, radius, 1)
  let dfa = calculate_dfa_alpha(values, 4, 64)
  let turning = turning_point_ratio(values)
  let histogram = histogram_entropy(values, 16)
  let lag1 = autocorrelation(values, 1)
  let complexity = (sample + approximate + histogram + turning) / 4.0
  {
    sample_entropy: sample,
    approximate_entropy: approximate,
    dfa_alpha: dfa,
    turning_point_ratio: turning,
    recurrence_rate: recurrence,
    histogram_entropy: histogram,
    lag1_autocorrelation: lag1,
    complexity_index: complexity,
  }
}

///|
/// Return the main non-linear features in stable column order.
pub fn nonlinear_feature_vector(values : Array[Double]) -> Array[Double] {
  let metrics = calculate_nonlinear_metrics(values)
  [
    metrics.sample_entropy,
    metrics.approximate_entropy,
    metrics.dfa_alpha,
    metrics.turning_point_ratio,
    metrics.recurrence_rate,
    metrics.histogram_entropy,
    metrics.lag1_autocorrelation,
    metrics.complexity_index,
  ]
}

///|
/// Compare two aligned sequences with a symmetric normalized distance.
pub fn normalized_sequence_distance(
  left : Array[Double],
  right : Array[Double],
) -> Double {
  let n = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  if n == 0 {
    0.0
  } else {
    let mut total = 0.0
    for i in 0..