///|
pub fn covariance(
  x : Array[Double],
  y : Array[Double],
  sample? : Bool = true,
) -> Double {
  if x.length() == 0 || x.length() != y.length() {
    return 0.0
  }
  let center_x = mean(x)
  let center_y = mean(y)
  let mut total = 0.0
  for index = 0; index < x.length(); index = index + 1 {
    total += (x[index] - center_x) * (y[index] - center_y)
  }
  let denominator = if sample && x.length() > 1 {
    x.length() - 1
  } else {
    x.length()
  }
  if denominator <= 0 {
    0.0
  } else {
    total / denominator.to_double()
  }
}

///|
pub fn pearson_correlation(x : Array[Double], y : Array[Double]) -> Double {
  if x.length() == 0 || x.length() != y.length() {
    return 0.0
  }
  let centered_x = []
  let centered_y = []
  let center_x = mean(x)
  let center_y = mean(y)
  for index = 0; index < x.length(); index = index + 1 {
    centered_x.push(x[index] - center_x)
    centered_y.push(y[index] - center_y)
  }
  let denominator = (sum_squared(centered_x) * sum_squared(centered_y)).sqrt()
  if denominator == 0.0 {
    0.0
  } else {
    dot_product(centered_x, centered_y) / denominator
  }
}

///|
pub fn spearman_correlation(x : Array[Double], y : Array[Double]) -> Double {
  if x.length() == 0 || x.length() != y.length() {
    0.0
  } else {
    pearson_correlation(ranks(x), ranks(y))
  }
}

///|
pub fn kendall_tau(x : Array[Double], y : Array[Double]) -> Double {
  if x.length() < 2 || x.length() != y.length() {
    return 0.0
  }
  let mut concordant = 0
  let mut discordant = 0
  let mut ties_x = 0
  let mut ties_y = 0
  for left = 0; left < x.length(); left = left + 1 {
    for right = left + 1; right < x.length(); right = right + 1 {
      let dx = x[left] - x[right]
      let dy = y[left] - y[right]
      if dx == 0.0 && dy == 0.0 {
        ties_x += 1
        ties_y += 1
      } else if dx == 0.0 {
        ties_x += 1
      } else if dy == 0.0 {
        ties_y += 1
      } else if dx * dy > 0.0 {
        concordant += 1
      } else {
        discordant += 1
      }
    }
  }
  let numerator = (concordant - discordant).to_double()
  let left_denominator = (concordant + discordant + ties_y).to_double()
  let right_denominator = (concordant + discordant + ties_x).to_double()
  let denominator = (left_denominator * right_denominator).sqrt()
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
pub fn cosine_similarity(x : Array[Double], y : Array[Double]) -> Double {
  if x.length() == 0 || x.length() != y.length() {
    return 0.0
  }
  let denominator = (sum_squared(x) * sum_squared(y)).sqrt()
  if denominator == 0.0 {
    0.0
  } else {
    dot_product(x, y) / denominator
  }
}

///|
pub fn mean_absolute_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() == 0 || actual.length() != predicted.length() {
    return 0.0
  }
  let mut total = 0.0
  for index = 0; index < actual.length(); index = index + 1 {
    total += abs_double(actual[index] - predicted[index])
  }
  total / actual.length().to_double()
}

///|
pub fn mean_squared_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() == 0 || actual.length() != predicted.length() {
    return 0.0
  }
  let mut total = 0.0
  for index = 0; index < actual.length(); index = index + 1 {
    let error = actual[index] - predicted[index]
    total += error * error
  }
  total / actual.length().to_double()
}

///|
pub fn root_mean_squared_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  mean_squared_error(actual, predicted).sqrt()
}

///|
pub fn median_absolute_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() == 0 || actual.length() != predicted.length() {
    return 0.0
  }
  let errors = []
  for index = 0; index < actual.length(); index = index + 1 {
    errors.push(abs_double(actual[index] - predicted[index]))
  }
  median(errors)
}

///|
pub fn explained_variance(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() == 0 || actual.length() != predicted.length() {
    return 0.0
  }
  let errors = []
  for index = 0; index < actual.length(); index = index + 1 {
    errors.push(actual[index] - predicted[index])
  }
  let denominator = population_variance(actual)
  if denominator == 0.0 {
    0.0
  } else {
    1.0 - population_variance(errors) / denominator
  }
}

///|
pub fn r_squared(actual : Array[Double], predicted : Array[Double]) -> Double {
  if actual.length() == 0 || actual.length() != predicted.length() {
    return 0.0
  }
  let center = mean(actual)
  let mut total = 0.0
  let mut residual = 0.0
  for index = 0; index < actual.length(); index = index + 1 {
    let total_error = actual[index] - center
    let residual_error = actual[index] - predicted[index]
    total += total_error * total_error
    residual += residual_error * residual_error
  }
  if total == 0.0 {
    0.0
  } else {
    1.0 - residual / total
  }
}

///|
pub fn robust_correlation(x : Array[Double], y : Array[Double]) -> Double {
  if x.length() == 0 || x.length() != y.length() {
    return 0.0
  }
  pearson_correlation(winsorize_by_z(x, 3.5), winsorize_by_z(y, 3.5))
}

///|
pub fn pairwise_covariance_matrix(
  data : Array[Array[Double]],
) -> Array[Array[Double]] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for left = 0; left < dimensions; left = left + 1 {
    let row = []
    for right = 0; right < dimensions; right = right + 1 {
      let x = []
      let y = []
      for observation in data {
        if observation.length() != dimensions {
          abort("inconsistent row dimensions")
        }
        x.push(observation[left])
        y.push(observation[right])
      }
      row.push(covariance(x, y))
    }
    result.push(row)
  }
  result
}

///|
pub fn correlation_matrix(data : Array[Array[Double]]) -> Array[Array[Double]] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for left = 0; left < dimensions; left = left + 1 {
    let row = []
    for right = 0; right < dimensions; right = right + 1 {
      let x = []
      let y = []
      for observation in data {
        if observation.length() != dimensions {
          abort("inconsistent row dimensions")
        }
        x.push(observation[left])
        y.push(observation[right])
      }
      row.push(pearson_correlation(x, y))
    }
    result.push(row)
  }
  result
}