///|
/// Robust statistics for a waveform window.
pub struct RobustWaveformStatistics {
  count : Int
  median : Double
  mad : Double
  interquartile_range : Double
  trimmed_mean : Double
  trimmed_rms : Double
  robust_peak : Double
  tail_fraction : Double
  outlier_count : Int
}

///|
/// Windowed stationarity diagnostics.
pub struct StationarityResult {
  window_count : Int
  mean_drift : Double
  variance_drift : Double
  mean_difference : Double
  variance_difference : Double
  stable_fraction : Double
  is_stationary : Bool
}

///|
/// One normalized autocorrelation lag.
pub struct AutocorrelationLag {
  lag : Int
  seconds : Double
  coefficient : Double
  energy : Double
}

///|
/// A histogram bin with density and cumulative probability.
pub struct WaveformHistogramBin {
  lower : Double
  upper : Double
  count : Int
  density : Double
  cumulative_fraction : Double
}

///|
/// A detected change point with before/after statistics.
pub struct WaveformChangePoint {
  index : Int
  score : Double
  before_mean : Double
  after_mean : Double
  before_variance : Double
  after_variance : Double
}

///|
/// Return the median absolute deviation around a supplied center.
pub fn absolute_deviation_median(
  data : Array[Double],
  center : Double,
) -> Double {
  median_samples(data.map(value => (value - center).abs()))
}

///|
/// Compute trimmed mean after dropping both tails.
pub fn trimmed_mean_samples(
  data : Array[Double],
  trim_fraction : Double,
) -> Double {
  if data.is_empty() {
    return 0.0
  }
  let values = data.copy()
  values.sort()
  let fraction = trim_fraction.max(0.0).min(0.49)
  let trim = (values.length().to_double() * fraction).to_int()
  let first = trim
  let last = values.length() - trim
  mean_samples(slice_samples(values, first, last))
}

///|
/// Compute a robust summary and mark observations beyond a MAD threshold.
pub fn robust_waveform_statistics(
  data : Array[Double],
  trim_fraction : Double,
  outlier_mad : Double,
) -> RobustWaveformStatistics {
  let median = median_samples(data)
  let mad = absolute_deviation_median(data, median)
  let scale = (mad * 1.4826).max(1.0e-12)
  let threshold = outlier_mad.abs() * scale
  let mut outliers = 0
  let mut tail = 0
  for value in data {
    if (value - median).abs() > threshold {
      outliers += 1
    }
    if value.abs() > median.abs() + threshold {
      tail += 1
    }
  }
  let lower = quantile_samples(data, 0.25)
  let upper = quantile_samples(data, 0.75)
  let trimmed = data.filter(value => value >= lower && value <= upper)
  {
    count: data.length(),
    median,
    mad,
    interquartile_range: upper - lower,
    trimmed_mean: trimmed_mean_samples(data, trim_fraction),
    trimmed_rms: rms_samples(trimmed),
    robust_peak: maximum_sample(data.map(value => (value - median).abs())),
    tail_fraction: if data.is_empty() {
      0.0
    } else {
      tail.to_double() / data.length().to_double()
    },
    outlier_count: outliers,
  }
}

///|
/// Return a normalized autocorrelation coefficient at one lag.
pub fn autocorrelation_at_lag(data : Array[Double], lag : Int) -> Double {
  if data.is_empty() || lag < 0 || lag >= data.length() {
    return 0.0
  }
  let mean = mean_samples(data)
  let mut numerator = 0.0
  let mut left_energy = 0.0
  let mut right_energy = 0.0
  for index = lag; index < data.length(); index = index + 1 {
    let left = data[index] - mean
    let right = data[index - lag] - mean
    numerator += left * right
    left_energy += left * left
    right_energy += right * right
  }
  let denominator = (left_energy * right_energy).sqrt()
  if denominator <= 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Calculate autocorrelation lags up to a maximum lag.
pub fn autocorrelation_lags(
  data : Array[Double],
  sample_rate : Double,
  maximum_lag : Int,
) -> Array[AutocorrelationLag] {
  let result = Array::new()
  let count = maximum_lag.max(0).min((data.length() - 1).max(0))
  for lag = 0; lag <= count; lag = lag + 1 {
    result.push({
      lag,
      seconds: if sample_rate <= 0.0 {
        0.0
      } else {
        lag.to_double() / sample_rate
      },
      coefficient: autocorrelation_at_lag(data, lag),
      energy: if lag == 0 {
        signal_energy(data)
      } else {
        signal_energy(slice_samples(data, lag, data.length()))
      },
    })
  }
  result
}

///|
/// Return the first non-zero-lag local autocorrelation maximum.
pub fn first_autocorrelation_peak(
  lags : Array[AutocorrelationLag],
) -> AutocorrelationLag? {
  if lags.length() < 3 {
    return None
  }
  for index = 1; index < lags.length() - 1; index = index + 1 {
    if lags[index].coefficient >= lags[index - 1].coefficient &&
      lags[index].coefficient >= lags[index + 1].coefficient {
      return Some(lags[index])
    }
  }
  None
}

///|
/// Compute a fixed-width amplitude histogram.
pub fn waveform_histogram(
  data : Array[Double],
  bins : Int,
  minimum : Double,
  maximum : Double,
) -> Array[WaveformHistogramBin] {
  let result : Array[WaveformHistogramBin] = Array::new()
  let count = bins.max(1)
  let low = minimum.min(maximum)
  let high = maximum.max(minimum)
  let width = (high - low).max(1.0e-12) / count.to_double()
  let counts = Array::make(count, 0)
  for value in data {
    let index = ((value - low) / width).to_int().max(0).min(count - 1)
    counts[index] += 1
  }
  let mut cumulative = 0
  for index = 0; index < count; index = index + 1 {
    cumulative += counts[index]
    result.push({
      lower: low + index.to_double() * width,
      upper: low + (index + 1).to_double() * width,
      count: counts[index],
      density: if data.is_empty() {
        0.0
      } else {
        counts[index].to_double() / data.length().to_double() / width
      },
      cumulative_fraction: if data.is_empty() {
        0.0
      } else {
        cumulative.to_double() / data.length().to_double()
      },
    })
  }
  result
}

///|
/// Compare adjacent windows to estimate signal stationarity.
pub fn assess_stationarity(
  data : Array[Double],
  window_size : Int,
  step : Int,
  mean_tolerance : Double,
  variance_tolerance : Double,
) -> StationarityResult {
  let windows = split_windows(data, window_size, step, false, 0.0)
  if windows.length() < 2 {
    return {
      window_count: windows.length(),
      mean_drift: 0.0,
      variance_drift: 0.0,
      mean_difference: 0.0,
      variance_difference: 0.0,
      stable_fraction: 1.0,
      is_stationary: true,
    }
  }
  let means = windows.map(window => mean_samples(window))
  let variances = windows.map(window => variance_samples(window))
  let mut stable = 0
  let mut mean_difference = 0.0
  let mut variance_difference = 0.0
  for index = 1; index < windows.length(); index = index + 1 {
    let mean_delta = (means[index] - means[index - 1]).abs()
    let variance_delta = (variances[index] - variances[index - 1]).abs()
    mean_difference += mean_delta
    variance_difference += variance_delta
    if mean_delta <= mean_tolerance.abs() &&
      variance_delta <= variance_tolerance.abs() {
      stable += 1
    }
  }
  let pairs = (windows.length() - 1).to_double()
  let stable_fraction = stable.to_double() / pairs
  {
    window_count: windows.length(),
    mean_drift: if pairs <= 0.0 {
      0.0
    } else {
      (means[means.length() - 1] - means[0]) / pairs
    },
    variance_drift: if pairs <= 0.0 {
      0.0
    } else {
      (variances[variances.length() - 1] - variances[0]) / pairs
    },
    mean_difference: mean_difference / pairs,
    variance_difference: variance_difference / pairs,
    stable_fraction,
    is_stationary: stable_fraction >= 0.8,
  }
}

///|
/// Score every candidate split using a mean/variance contrast.
pub fn change_point_scores(
  data : Array[Double],
  minimum_window : Int,
) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  let window = minimum_window.max(2)
  if data.length() < window * 2 {
    return result
  }
  for split = window; split < data.length() - window; split = split + 1 {
    let left = slice_samples(data, split - window, split)
    let right = slice_samples(data, split, split + window)
    let mean_delta = (mean_samples(right) - mean_samples(left)).abs()
    let variance_delta = (variance_samples(right) - variance_samples(left))
      .abs()
      .sqrt()
    let scale = (standard_deviation(left) + standard_deviation(right)).max(
      1.0e-12,
    )
    result[split] = (mean_delta + variance_delta) / scale
  }
  result
}

///|
/// Return the strongest separated change points.
pub fn strongest_change_points(
  data : Array[Double],
  minimum_window : Int,
  minimum_separation : Int,
  minimum_score : Double,
) -> Array[WaveformChangePoint] {
  let scores = change_point_scores(data, minimum_window)
  let candidates = Array::new()
  for index = 0; index < scores.length(); index = index + 1 {
    if scores[index] >= minimum_score {
      candidates.push(index)
    }
  }
  candidates.sort_by((left, right) => {
    if scores[left] > scores[right] {
      -1
    } else if scores[left] < scores[right] {
      1
    } else {
      0
    }
  })
  let result : Array[WaveformChangePoint] = Array::new()
  for index in candidates {
    let too_close = result.any(point => {
      (point.index - index).abs() < minimum_separation.max(1)
    })
    if !too_close {
      let left = slice_samples(data, index - minimum_window.max(2), index)
      let right = slice_samples(data, index, index + minimum_window.max(2))
      result.push({
        index,
        score: scores[index],
        before_mean: mean_samples(left),
        after_mean: mean_samples(right),
        before_variance: variance_samples(left),
        after_variance: variance_samples(right),
      })
    }
  }
  result.sort_by((left, right) => left.index - right.index)
  result
}

///|
/// Compute an RMS envelope over fixed windows.
pub fn rms_envelope(
  data : Array[Double],
  window_size : Int,
  step : Int,
) -> Array[Double] {
  let result = Array::new()
  if window_size <= 0 || step <= 0 {
    return result
  }
  let mut start = 0
  while start + window_size <= data.length() {
    result.push(rms_samples(slice_samples(data, start, start + window_size)))
    start += step
  }
  result
}

///|
/// Compute an amplitude percentile envelope over windows.
pub fn percentile_envelope(
  data : Array[Double],
  window_size : Int,
  step : Int,
  probability : Double,
) -> Array[Double] {
  let result = Array::new()
  if window_size <= 0 || step <= 0 {
    return result
  }
  let mut start = 0
  while start + window_size <= data.length() {
    result.push(
      quantile_samples(
        slice_samples(data, start, start + window_size),
        probability,
      ),
    )
    start += step
  }
  result
}

///|
/// Return a deterministic bootstrap mean confidence interval.
pub fn bootstrap_mean_interval(
  data : Array[Double],
  resamples : Int,
  seed : Int,
) -> (Double, Double) {
  if data.is_empty() || resamples <= 0 {
    return (0.0, 0.0)
  }
  let means = Array::new()
  let mut state = if seed == 0 { 1 } else { seed.abs() }
  for _ in 0.. String {
  let mut output = "index,score,before_mean,after_mean,before_variance,after_variance\n"
  for point in points {
    output = output +
      point.index.to_string() +
      "," +
      point.score.to_string() +
      "," +
      point.before_mean.to_string() +
      "," +
      point.after_mean.to_string() +
      "," +
      point.before_variance.to_string() +
      "," +
      point.after_variance.to_string() +
      "\n"
  }
  output
}

///|
/// Format robust statistics for quality reports.
pub fn robust_statistics_text(statistics : RobustWaveformStatistics) -> String {
  "count=" +
  statistics.count.to_string() +
  " median=" +
  statistics.median.to_string() +
  " mad=" +
  statistics.mad.to_string() +
  " iqr=" +
  statistics.interquartile_range.to_string() +
  " trimmed_mean=" +
  statistics.trimmed_mean.to_string() +
  " robust_peak=" +
  statistics.robust_peak.to_string() +
  " tail_fraction=" +
  statistics.tail_fraction.to_string() +
  " outliers=" +
  statistics.outlier_count.to_string()
}