///|
/// A pairwise waveform comparison result.
pub struct WaveformComparison {
  sample_count : Int
  mean_absolute_error : Double
  root_mean_square_error : Double
  maximum_absolute_error : Double
  correlation : Double
  normalized_distance : Double
  lag_samples : Int
}

///|
/// Compute absolute, squared, and signed error summaries.
pub fn compare_waveforms(
  left : Array[Double],
  right : Array[Double],
  lag_samples : Int,
) -> WaveformComparison {
  let count = left.length().min(right.length())
  if count == 0 {
    return {
      sample_count: 0,
      mean_absolute_error: 0.0,
      root_mean_square_error: 0.0,
      maximum_absolute_error: 0.0,
      correlation: 0.0,
      normalized_distance: 0.0,
      lag_samples,
    }
  }
  let mut absolute = 0.0
  let mut squared = 0.0
  let mut maximum = 0.0
  let mut dot = 0.0
  let mut left_energy = 0.0
  let mut right_energy = 0.0
  for index = 0; index < count; index = index + 1 {
    let error = left[index] - right[index]
    absolute += error.abs()
    squared += error * error
    maximum = maximum.max(error.abs())
    dot += left[index] * right[index]
    left_energy += left[index] * left[index]
    right_energy += right[index] * right[index]
  }
  let denominator = (left_energy * right_energy).sqrt()
  let distance = (squared / (left_energy + right_energy).max(1.0e-20)).sqrt()
  {
    sample_count: count,
    mean_absolute_error: absolute / count.to_double(),
    root_mean_square_error: (squared / count.to_double()).sqrt(),
    maximum_absolute_error: maximum,
    correlation: if denominator <= 0.0 {
      0.0
    } else {
      dot / denominator
    },
    normalized_distance: distance,
    lag_samples,
  }
}

///|
/// Return a lagged copy with zero padding.
pub fn lagged_samples(data : Array[Double], lag : Int) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  for index = 0; index < data.length(); index = index + 1 {
    let source = index - lag
    if source >= 0 && source < data.length() {
      result[index] = data[source]
    }
  }
  result
}

///|
/// Find the integer lag maximizing normalized correlation.
pub fn best_integer_lag(
  left : Array[Double],
  right : Array[Double],
  maximum_lag : Int,
) -> Int {
  let limit = maximum_lag.abs()
  let mut best_lag = 0
  let mut best_score = -1.0e300
  for lag = -limit; lag <= limit; lag = lag + 1 {
    let shifted = lagged_samples(right, lag)
    let comparison = compare_waveforms(left, shifted, lag)
    if comparison.correlation > best_score {
      best_score = comparison.correlation
      best_lag = lag
    }
  }
  best_lag
}

///|
/// Compute a normalized cross-spectrum magnitude between two arrays.
pub fn cross_spectrum_magnitude(
  left : Array[Double],
  right : Array[Double],
) -> Array[Double] {
  let count = left.length().min(right.length())
  let left_spectrum = real_dft(slice_samples(left, 0, count))
  let right_spectrum = real_dft(slice_samples(right, 0, count))
  let result = Array::make(
    left_spectrum.length().min(right_spectrum.length()),
    0.0,
  )
  for index = 0; index < result.length(); index = index + 1 {
    result[index] = left_spectrum[index].magnitude() *
      right_spectrum[index].magnitude()
  }
  result
}

///|
/// Calculate a band-limited coherence proxy.
pub fn waveform_band_coherence(
  left : Array[Double],
  right : Array[Double],
  sample_rate : Double,
  band : SpectralBand,
) -> Double {
  let count = left.length().min(right.length())
  if count == 0 {
    return 0.0
  }
  let left_spectrum = hann_periodogram(
    slice_samples(left, 0, count),
    sample_rate,
  )
  let right_spectrum = hann_periodogram(
    slice_samples(right, 0, count),
    sample_rate,
  )
  let distance = advanced_spectral_distance(left_spectrum, right_spectrum)
  let left_power = advanced_band_power(left_spectrum, band)
  let right_power = advanced_band_power(right_spectrum, band)
  if left_power <= 0.0 || right_power <= 0.0 {
    0.0
  } else {
    (1.0 - distance).max(0.0).min(1.0)
  }
}

///|
/// Apply a scalar gain and offset to a copied waveform.
pub fn affine_waveform(
  data : Array[Double],
  gain : Double,
  offset : Double,
) -> Array[Double] {
  data.map(value => value * gain + offset)
}

///|
/// Mix two waveforms with a clipped blend factor.
pub fn blend_waveforms(
  left : Array[Double],
  right : Array[Double],
  right_weight : Double,
) -> Array[Double] {
  let count = left.length().min(right.length())
  let result = Array::make(count, 0.0)
  let weight = right_weight.max(0.0).min(1.0)
  for index = 0; index < count; index = index + 1 {
    result[index] = left[index] * (1.0 - weight) + right[index] * weight
  }
  result
}

///|
/// Return a stable comparison summary string.
pub fn waveform_comparison_text(comparison : WaveformComparison) -> String {
  "samples=" +
  comparison.sample_count.to_string() +
  " mae=" +
  comparison.mean_absolute_error.to_string() +
  " rmse=" +
  comparison.root_mean_square_error.to_string() +
  " max_error=" +
  comparison.maximum_absolute_error.to_string() +
  " correlation=" +
  comparison.correlation.to_string() +
  " distance=" +
  comparison.normalized_distance.to_string() +
  " lag=" +
  comparison.lag_samples.to_string()
}