///|
/// Descriptive statistics used by HRV reports and signal diagnostics.
pub(all) struct DistributionStats {
  count : Int
  sum : Double
  mean : Double
  median : Double
  variance : Double
  standard_deviation : Double
  minimum : Double
  maximum : Double
  q1 : Double
  q3 : Double
  interquartile_range : Double
  median_absolute_deviation : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// A least-squares trend fitted to an ordered sequence.
pub(all) struct RegressionSummary {
  slope : Double
  intercept : Double
  correlation : Double
  r_squared : Double
  residual_standard_error : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Sum a sequence without changing its order.
pub fn sum_values(values : Array[Double]) -> Double {
  let mut total = 0.0
  for value in values {
    total += value
  }
  total
}

///|
/// Return the arithmetic mean, or zero for an empty sequence.
pub fn mean_value(values : Array[Double]) -> Double {
  if values.length() == 0 {
    0.0
  } else {
    sum_values(values) / values.length().to_double()
  }
}

///|
/// Copy and sort a sequence. The input is never mutated.
fn sorted_values(values : Array[Double]) -> Array[Double] {
  let sorted = []
  for value in values {
    sorted.push(value)
  }
  sorted.sort()
  sorted
}

///|
/// Copy the first n values into an owned array for APIs that require Array.
fn prefix_values(values : Array[Double], n : Int) -> Array[Double] {
  let result = []
  let limit = if n < values.length() { n } else { values.length() }
  for i in 0.. Double {
  let sorted = sorted_values(values)
  let n = sorted.length()
  if n == 0 {
    0.0
  } else if n % 2 == 1 {
    sorted[n / 2]
  } else {
    (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0
  }
}

///|
/// Compute a linearly interpolated quantile in the closed interval [0, 1].
pub fn quantile_value(values : Array[Double], probability : Double) -> Double {
  let sorted = sorted_values(values)
  let n = sorted.length()
  if n == 0 {
    0.0
  } else {
    let p = probability.clamp(min=0.0, max=1.0)
    let position = p * (n - 1).to_double()
    let lower = position.floor().to_int()
    let upper = position.ceil().to_int()
    if lower == upper {
      sorted[lower]
    } else {
      let fraction = position - lower.to_double()
      sorted[lower] + fraction * (sorted[upper] - sorted[lower])
    }
  }
}

///|
/// Sample variance. A sequence with fewer than two values has zero variance.
pub fn variance_value(values : Array[Double]) -> Double {
  let n = values.length()
  if n <= 1 {
    0.0
  } else {
    let average = mean_value(values)
    let mut squared = 0.0
    for value in values {
      let delta = value - average
      squared += delta * delta
    }
    squared / (n - 1).to_double()
  }
}

///|
/// Population variance, useful for complete windows rather than samples.
pub fn population_variance(values : Array[Double]) -> Double {
  let n = values.length()
  if n == 0 {
    0.0
  } else {
    let average = mean_value(values)
    let mut squared = 0.0
    for value in values {
      let delta = value - average
      squared += delta * delta
    }
    squared / n.to_double()
  }
}

///|
/// Return the sample standard deviation.
pub fn standard_deviation(values : Array[Double]) -> Double {
  variance_value(values).sqrt()
}

///|
/// Return the median absolute deviation around the sample median.
pub fn median_absolute_deviation(values : Array[Double]) -> Double {
  let center = median_value(values)
  let deviations = []
  for value in values {
    let delta = value - center
    deviations.push(if delta < 0.0 { -delta } else { delta })
  }
  median_value(deviations)
}

///|
/// Backward-compatible mean absolute deviation around the arithmetic mean.
pub fn mean_absolute_deviation(values : Array[Double]) -> Double {
  if values.length() == 0 {
    return 0.0
  }
  let center = mean_value(values)
  let deviations = []
  for value in values {
    deviations.push(absolute_difference(value, center))
  }
  mean_value(deviations)
}

///|
/// Return the coefficient of variation as a percentage.
pub fn coefficient_of_variation(values : Array[Double]) -> Double {
  let average = mean_value(values)
  if average == 0.0 {
    0.0
  } else {
    standard_deviation(values) / average.abs() * 100.0
  }
}

///|
/// Return a full distribution summary for a sequence.
pub fn summarize_distribution(values : Array[Double]) -> DistributionStats {
  let n = values.length()
  if n == 0 {
    return {
      count: 0,
      sum: 0.0,
      mean: 0.0,
      median: 0.0,
      variance: 0.0,
      standard_deviation: 0.0,
      minimum: 0.0,
      maximum: 0.0,
      q1: 0.0,
      q3: 0.0,
      interquartile_range: 0.0,
      median_absolute_deviation: 0.0,
    }
  }
  let sorted = sorted_values(values)
  let q1 = quantile_value(sorted, 0.25)
  let q3 = quantile_value(sorted, 0.75)
  {
    count: n,
    sum: sum_values(values),
    mean: mean_value(values),
    median: median_value(sorted),
    variance: variance_value(values),
    standard_deviation: standard_deviation(values),
    minimum: sorted[0],
    maximum: sorted[n - 1],
    q1,
    q3,
    interquartile_range: q3 - q1,
    median_absolute_deviation: median_absolute_deviation(values),
  }
}

///|
/// Return covariance using the sample denominator.
pub fn covariance_value(left : Array[Double], right : Array[Double]) -> Double {
  let n = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  if n <= 1 {
    0.0
  } else {
    let left_prefix = prefix_values(left, n)
    let right_prefix = prefix_values(right, n)
    let left_mean = mean_value(left_prefix)
    let right_mean = mean_value(right_prefix)
    let mut total = 0.0
    for i in 0.. Double {
  let n = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  if n <= 1 {
    0.0
  } else {
    let covariance = covariance_value(left, right)
    let left_sd = standard_deviation(prefix_values(left, n))
    let right_sd = standard_deviation(prefix_values(right, n))
    if left_sd == 0.0 || right_sd == 0.0 {
      0.0
    } else {
      covariance / (left_sd * right_sd)
    }
  }
}

///|
/// Fit y = slope * x + intercept to an ordered sequence of y values.
pub fn fit_linear_trend(values : Array[Double]) -> RegressionSummary {
  let n = values.length()
  if n <= 1 {
    return {
      slope: 0.0,
      intercept: if n == 1 {
        values[0]
      } else {
        0.0
      },
      correlation: 0.0,
      r_squared: 0.0,
      residual_standard_error: 0.0,
    }
  }
  let x = []
  for i in 0.. Array[Double] {
  let result = []
  if window_size <= 0 {
    return result
  }
  for i in 0.. Array[Double] {
  let result = []
  if window_size <= 0 {
    return result
  }
  for i in 0.. Array[Double] {
  let result = []
  let average = mean_value(values)
  let deviation = standard_deviation(values)
  for value in values {
    result.push(
      if deviation == 0.0 {
        0.0
      } else {
        (value - average) / deviation
      },
    )
  }
  result
}

///|
/// Normalize values into [0, 1]. Constant sequences map to 0.5.
pub fn normalize_unit_interval(values : Array[Double]) -> Array[Double] {
  let result = []
  let summary = summarize_distribution(values)
  let width = summary.maximum - summary.minimum
  for value in values {
    result.push(
      if width == 0.0 {
        0.5
      } else {
        (value - summary.minimum) / width
      },
    )
  }
  result
}

///|
/// Calculate a weighted mean for paired values. Extra values are ignored.
pub fn weighted_mean(values : Array[Double], weights : Array[Double]) -> Double {
  let n = if values.length() < weights.length() {
    values.length()
  } else {
    weights.length()
  }
  let mut total = 0.0
  let mut total_weight = 0.0
  for i in 0..