///|
/// Agreement statistics for two synchronized measurement channels.
pub struct SignalAgreement {
  sample_count : Int
  mean_absolute_error : Double
  mean_squared_error : Double
  root_mean_squared_error : Double
  mean_bias : Double
  correlation : Double
  maximum_absolute_error : Double
} derive(Debug, Eq)

///|
fn valid_pair(left : Double, right : Double) -> Bool {
  !left.is_nan() && !left.is_inf() && !right.is_nan() && !right.is_inf()
}

///|
/// Compares finite pairs with an unbiased Pearson correlation estimate.
pub fn compare_signals(
  reference : Array[Double],
  measured : Array[Double],
) -> SignalAgreement {
  let n = if reference.length() < measured.length() {
    reference.length()
  } else {
    measured.length()
  }
  let left : Array[Double] = []
  let right : Array[Double] = []
  for i = 0; i < n; i = i + 1 {
    if valid_pair(reference[i], measured[i]) {
      left.push(reference[i])
      right.push(measured[i])
    }
  }
  if left.length() == 0 {
    return {
      sample_count: 0,
      mean_absolute_error: 0.0,
      mean_squared_error: 0.0,
      root_mean_squared_error: 0.0,
      mean_bias: 0.0,
      correlation: 0.0,
      maximum_absolute_error: 0.0,
    }
  }
  let mut abs_error = 0.0
  let mut squared_error = 0.0
  let mut bias = 0.0
  let mut maximum = 0.0
  for i = 0; i < left.length(); i = i + 1 {
    let difference = measured[i] - reference[i]
    abs_error = abs_error + difference.abs()
    squared_error = squared_error + difference * difference
    bias = bias + difference
    if difference.abs() > maximum {
      maximum = difference.abs()
    }
  }
  let mean_left = compute_mean(left)
  let mean_right = compute_mean(right)
  let mut covariance = 0.0
  let mut variance_left = 0.0
  let mut variance_right = 0.0
  for i = 0; i < left.length(); i = i + 1 {
    let dl = left[i] - mean_left
    let dr = right[i] - mean_right
    covariance = covariance + dl * dr
    variance_left = variance_left + dl * dl
    variance_right = variance_right + dr * dr
  }
  let denominator = (variance_left * variance_right).sqrt()
  {
    sample_count: left.length(),
    mean_absolute_error: abs_error / left.length().to_double(),
    mean_squared_error: squared_error / left.length().to_double(),
    root_mean_squared_error: (squared_error / left.length().to_double()).sqrt(),
    mean_bias: bias / left.length().to_double(),
    correlation: if denominator <= 0.000000001 {
      0.0
    } else {
      covariance / denominator
    },
    maximum_absolute_error: maximum,
  }
}

///|
/// Maps mean absolute error to a bounded 0--100 agreement score.
pub fn compute_signal_agreement_score(
  reference : Array[Double],
  measured : Array[Double],
  tolerance : Double,
) -> Double {
  if tolerance <= 0.0 || tolerance.is_nan() || tolerance.is_inf() {
    return 0.0
  }
  let agreement = compare_signals(reference, measured)
  if agreement.sample_count == 0 {
    0.0
  } else {
    (100.0 * (1.0 - agreement.mean_absolute_error / tolerance)).clamp(
      min=0.0,
      max=100.0,
    )
  }
}

///|
/// Correlates two channels with a non-negative integer sample lag.
pub fn lagged_correlation(
  reference : Array[Double],
  measured : Array[Double],
  lag : Int,
) -> Double {
  let safe_lag = if lag < 0 { 0 } else { lag }
  if safe_lag >= reference.length() || safe_lag >= measured.length() {
    return 0.0
  }
  let left : Array[Double] = []
  let right : Array[Double] = []
  let n = if reference.length() < measured.length() - safe_lag {
    reference.length()
  } else {
    measured.length() - safe_lag
  }
  for i = 0; i < n; i = i + 1 {
    if valid_pair(reference[i], measured[i + safe_lag]) {
      left.push(reference[i])
      right.push(measured[i + safe_lag])
    }
  }
  compare_signals(left, right).correlation
}