///|
/// Descriptive statistics for one finite signal. Non-finite samples are
/// excluded from all numeric estimates while `count` records usable samples.
pub struct SignalSummary {
  count : Int
  mean : Double
  standard_deviation : Double
  minimum : Double
  maximum : Double
  median : Double
  first_quartile : Double
  third_quartile : Double
  interquartile_range : Double
  rms : Double
  finite_fraction : Double
} derive(Debug, Eq)

///|
fn finite_values(values : Array[Double]) -> Array[Double] {
  let result : Array[Double] = Array::new(capacity=values.length())
  for value in values {
    if !value.is_nan() && !value.is_inf() {
      result.push(value)
    }
  }
  result
}

///|
/// Computes a numerically stable two-pass standard deviation for finite data.
fn finite_standard_deviation(values : Array[Double], mean : Double) -> Double {
  if values.length() <= 1 {
    return 0.0
  }
  let mut squared = 0.0
  for value in values {
    let delta = value - mean
    squared = squared + delta * delta
  }
  (squared / (values.length() - 1).to_double()).sqrt()
}

///|
/// Summarizes a signal with robust quartiles and a population-style RMS.
pub fn summarize_signal(values : Array[Double]) -> SignalSummary {
  let finite = finite_values(values)
  let count = finite.length()
  if count == 0 {
    return {
      count: 0,
      mean: 0.0,
      standard_deviation: 0.0,
      minimum: 0.0,
      maximum: 0.0,
      median: 0.0,
      first_quartile: 0.0,
      third_quartile: 0.0,
      interquartile_range: 0.0,
      rms: 0.0,
      finite_fraction: 0.0,
    }
  }
  let mut sum = 0.0
  let mut square_sum = 0.0
  let mut minimum = finite[0]
  let mut maximum = finite[0]
  for value in finite {
    sum = sum + value
    square_sum = square_sum + value * value
    if value < minimum {
      minimum = value
    }
    if value > maximum {
      maximum = value
    }
  }
  let mean = sum / count.to_double()
  let q1 = quantile(finite, 0.25)
  let q3 = quantile(finite, 0.75)
  {
    count,
    mean,
    standard_deviation: finite_standard_deviation(finite, mean),
    minimum,
    maximum,
    median: quantile(finite, 0.5),
    first_quartile: q1,
    third_quartile: q3,
    interquartile_range: q3 - q1,
    rms: (square_sum / count.to_double()).sqrt(),
    finite_fraction: count.to_double() / values.length().to_double(),
  }
}

///|
/// Returns a clamped half-open window. Reversed windows become empty.
pub fn window_indices(length : Int, start : Int, end : Int) -> SampleRange {
  clamp_range(length, start, end)
}

///|
/// Copies the vertical force samples in a safe, half-open window.
pub fn extract_force_window(
  trial : ForceTrial,
  start : Int,
  end : Int,
) -> Array[Double] {
  let force = trial.get_fz()
  copy_range(force, start, end)
}

///|
/// Copies COP samples in a safe, half-open window.
pub fn extract_cop_window(
  trial : ForceTrial,
  start : Int,
  end : Int,
) -> Array[CoPPoint] {
  let points = trial.get_cop_points()
  let range = clamp_range(points.length(), start, end)
  let result : Array[CoPPoint] = Array::new(capacity=range.end - range.start)
  for i = range.start; i < range.end; i = i + 1 {
    result.push(points[i])
  }
  result
}

///|
/// Feature vector shared by quality dashboards and downstream pipelines.
pub struct TrialFeatureSummary {
  sample_count : Int
  duration_s : Double
  sampling_rate_hz : Double
  body_mass_kg : Double
  force : SignalSummary
  cop_path_length : Double
  cop_mean_x : Double
  cop_mean_y : Double
  cop_range_x : Double
  cop_range_y : Double
  missing_cop_fraction : Double
} derive(Debug, Eq)

///|
fn cop_finite(point : CoPPoint) -> Bool {
  !point.x.is_nan() &&
  !point.x.is_inf() &&
  !point.y.is_nan() &&
  !point.y.is_inf()
}

///|
fn cop_path_length_finite(points : Array[CoPPoint]) -> Double {
  let mut result = 0.0
  let mut previous : CoPPoint? = None
  for point in points {
    if cop_finite(point) {
      match previous {
        Some(last) => result = result + last.distance_to(point)
        None => ()
      }
      previous = Some(point)
    } else {
      previous = None
    }
  }
  result
}

///|
/// Produces a compact finite-aware feature summary of a force plate trial.
pub fn summarize_trial(trial : ForceTrial) -> TrialFeatureSummary {
  let forces = trial.get_fz()
  let points = trial.get_cop_points()
  let timestamps = trial.get_timestamps()
  let force = summarize_signal(forces)
  let timebase = summarize_timebase(timestamps)
  let mut sum_x = 0.0
  let mut sum_y = 0.0
  let mut count = 0
  let mut min_x = 0.0
  let mut max_x = 0.0
  let mut min_y = 0.0
  let mut max_y = 0.0
  for point in points {
    if cop_finite(point) {
      if count == 0 {
        min_x = point.x
        max_x = point.x
        min_y = point.y
        max_y = point.y
      } else {
        if point.x < min_x {
          min_x = point.x
        }
        if point.x > max_x {
          max_x = point.x
        }
        if point.y < min_y {
          min_y = point.y
        }
        if point.y > max_y {
          max_y = point.y
        }
      }
      sum_x = sum_x + point.x
      sum_y = sum_y + point.y
      count = count + 1
    }
  }
  let mean_x = if count == 0 { 0.0 } else { sum_x / count.to_double() }
  let mean_y = if count == 0 { 0.0 } else { sum_y / count.to_double() }
  {
    sample_count: trial.length(),
    duration_s: trial.duration(),
    sampling_rate_hz: timebase.nominal_rate_hz,
    body_mass_kg: trial.metadata.nominal_body_mass,
    force,
    cop_path_length: cop_path_length_finite(points),
    cop_mean_x: mean_x,
    cop_mean_y: mean_y,
    cop_range_x: if count == 0 {
      0.0
    } else {
      max_x - min_x
    },
    cop_range_y: if count == 0 {
      0.0
    } else {
      max_y - min_y
    },
    missing_cop_fraction: if points.length() == 0 {
      0.0
    } else {
      1.0 - count.to_double() / points.length().to_double()
    },
  }
}