///|
/// Poincare-plot descriptors of short-term and long-term variability.
pub(all) struct PoincareMetrics {
  sd1 : Double
  sd2 : Double
  sd1_sd2_ratio : Double
  ellipse_area : Double
  center_rr : Double
  points : Int
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Geometric HRV descriptors derived from an RR histogram.
pub(all) struct GeometricMetrics {
  triangular_index : Double
  tinn : Double
  mode_rr : Double
  mode_count : Int
  bin_width : Double
  bins : Array[Int]
} derive(FromJson, ToJson, Debug, Eq)

///|
/// A compact set of heart-rate zones for training dashboards.
pub(all) struct HeartRateSummary {
  mean_bpm : Double
  minimum_bpm : Double
  maximum_bpm : Double
  median_bpm : Double
  resting_bpm : Double
  zone1_seconds : Double
  zone2_seconds : Double
  zone3_seconds : Double
  zone4_seconds : Double
  zone5_seconds : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Calculate SD1 and SD2 from successive RR pairs.
pub fn calculate_poincare(intervals : Array[Double]) -> PoincareMetrics {
  let points = if intervals.length() > 1 { intervals.length() - 1 } else { 0 }
  if points == 0 {
    return {
      sd1: 0.0,
      sd2: 0.0,
      sd1_sd2_ratio: 0.0,
      ellipse_area: 0.0,
      center_rr: mean_value(intervals),
      points: 0,
    }
  }
  let differences = []
  let sums = []
  for i in 0.. Array[Int] {
  let bins = Array::make(if bin_count < 0 { 0 } else { bin_count }, 0)
  if bin_width <= 0.0 || bin_count <= 0 {
    return bins
  }
  for value in intervals {
    let position = ((value - minimum) / bin_width).floor().to_int()
    if position >= 0 && position < bins.length() {
      bins[position] += 1
    }
  }
  bins
}

///|
/// Return the most populated histogram bin, preferring the lower bin on ties.
pub fn histogram_mode(
  bins : Array[Int],
  minimum : Double,
  bin_width : Double,
) -> (Double, Int) {
  if bins.length() == 0 {
    return (minimum, 0)
  }
  let mut best_index = 0
  let mut best_count = bins[0]
  for i in 1.. best_count {
      best_index = i
      best_count = bins[i]
    }
  }
  (minimum + (best_index.to_double() + 0.5) * bin_width, best_count)
}

///|
/// Calculate histogram-based triangular index and an approximate TINN width.
pub fn calculate_geometric_metrics(
  intervals : Array[Double],
  bin_width : Double,
) -> GeometricMetrics {
  if intervals.length() == 0 || bin_width <= 0.0 {
    return {
      triangular_index: 0.0,
      tinn: 0.0,
      mode_rr: 0.0,
      mode_count: 0,
      bin_width,
      bins: [],
    }
  }
  let summary = summarize_distribution(intervals)
  let range = summary.maximum - summary.minimum
  let count = (range / bin_width).ceil().to_int() + 1
  let bins = rr_histogram(intervals, bin_width, summary.minimum, count)
  let (mode_rr, mode_count) = histogram_mode(bins, summary.minimum, bin_width)
  let peak = if mode_count == 0 { 1 } else { mode_count }
  let half = if peak <= 1 { 1 } else { peak / 2 }
  let mut first = 0
  while first < bins.length() && bins[first] < half {
    first += 1
  }
  let mut last = bins.length() - 1
  while last >= 0 && bins[last] < half {
    last -= 1
  }
  let tinn = if first <= last {
    (last - first + 1).to_double() * bin_width
  } else {
    0.0
  }
  {
    triangular_index: if mode_count == 0 {
      0.0
    } else {
      intervals.length().to_double() / mode_count.to_double()
    },
    tinn,
    mode_rr,
    mode_count,
    bin_width,
    bins,
  }
}

///|
/// Convenience wrapper using a 7.8125 ms histogram bin.
pub fn calculate_triangular_index(intervals : Array[Double]) -> Double {
  calculate_geometric_metrics(intervals, 7.8125).triangular_index
}

///|
/// Approximate the triangular interpolation of the NN interval histogram.
pub fn calculate_tinn(intervals : Array[Double]) -> Double {
  calculate_geometric_metrics(intervals, 7.8125).tinn
}

///|
/// Convert RR intervals in milliseconds to instantaneous heart rate.
pub fn intervals_to_heart_rate(intervals : Array[Double]) -> Array[Double] {
  let result = []
  for interval in intervals {
    result.push(if interval <= 0.0 { 0.0 } else { 60000.0 / interval })
  }
  result
}

///|
/// Calculate a bounded heart-rate summary and time in five percentage zones.
pub fn summarize_heart_rate(
  intervals : Array[Double],
  maximum_hr : Double,
) -> HeartRateSummary {
  let heart_rates = intervals_to_heart_rate(intervals)
  if heart_rates.length() == 0 {
    return {
      mean_bpm: 0.0,
      minimum_bpm: 0.0,
      maximum_bpm: 0.0,
      median_bpm: 0.0,
      resting_bpm: 0.0,
      zone1_seconds: 0.0,
      zone2_seconds: 0.0,
      zone3_seconds: 0.0,
      zone4_seconds: 0.0,
      zone5_seconds: 0.0,
    }
  }
  let summary = summarize_distribution(heart_rates)
  let peak = if maximum_hr <= 0.0 { summary.maximum } else { maximum_hr }
  let mut zone1 = 0.0
  let mut zone2 = 0.0
  let mut zone3 = 0.0
  let mut zone4 = 0.0
  let mut zone5 = 0.0
  for i in 0.. Double {
  let mean_rr = mean_value(intervals)
  let rmssd = calculate_rmssd(intervals)
  if mean_rr <= 0.0 {
    0.0
  } else {
    100.0 * rmssd / mean_rr
  }
}

///|
/// Calculate a stress index based on histogram mode and spread.
pub fn baevsky_stress_index(intervals : Array[Double]) -> Double {
  if intervals.length() == 0 {
    return 0.0
  }
  let geometric = calculate_geometric_metrics(intervals, 7.8125)
  let range = geometric.tinn
  if range <= 0.0 || geometric.mode_rr <= 0.0 {
    0.0
  } else {
    1000.0 *
    geometric.mode_count.to_double() /
    (2.0 * range * intervals.length().to_double())
  }
}

///|
/// Return a robust HRV feature vector for downstream scoring.
pub fn advanced_feature_vector(intervals : Array[Double]) -> Array[Double] {
  let poincare = calculate_poincare(intervals)
  let geometric = calculate_geometric_metrics(intervals, 7.8125)
  [
    mean_value(intervals),
    median_value(intervals),
    standard_deviation(intervals),
    median_absolute_deviation(intervals),
    coefficient_of_variation(intervals),
    calculate_sdnn(intervals),
    calculate_rmssd(intervals),
    calculate_pnn(intervals, 20.0),
    calculate_pnn(intervals, 50.0),
    poincare.sd1,
    poincare.sd2,
    poincare.sd1_sd2_ratio,
    geometric.triangular_index,
    geometric.tinn,
    cardiac_vagal_index(intervals),
    baevsky_stress_index(intervals),
  ]
}