///|
/// Population variance. Returns zero for an empty or one-element sample.
pub fn population_variance(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let center = mean(data)
  let mut total = 0.0
  for value in data {
    let deviation = value - center
    total += deviation * deviation
  }
  total / data.length().to_double()
}

///|
pub fn sample_variance(data : Array[Double]) -> Double {
  if data.length() <= 1 {
    return 0.0
  }
  let center = mean(data)
  let mut total = 0.0
  for value in data {
    let deviation = value - center
    total += deviation * deviation
  }
  total / (data.length() - 1).to_double()
}

///|
pub fn population_stddev(data : Array[Double]) -> Double {
  population_variance(data).sqrt()
}

///|
pub fn sample_stddev(data : Array[Double]) -> Double {
  sample_variance(data).sqrt()
}

///|
pub fn mean_absolute_deviation(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let center = mean(data)
  let mut total = 0.0
  for value in data {
    total += abs_double(value - center)
  }
  total / data.length().to_double()
}

///|
pub fn central_moment(data : Array[Double], order : Int) -> Double {
  if order < 0 {
    abort("moment order must be non-negative")
  }
  if data.length() == 0 {
    return 0.0
  }
  if order == 0 {
    return 1.0
  }
  let center = mean(data)
  let mut total = 0.0
  for value in data {
    let mut power = 1.0
    let deviation = value - center
    for power_index = 0; power_index < order; power_index = power_index + 1 {
      power *= deviation
    }
    total += power
  }
  total / data.length().to_double()
}

///|
pub fn skewness(data : Array[Double]) -> Double {
  let deviation = population_stddev(data)
  if deviation == 0.0 {
    0.0
  } else {
    central_moment(data, 3) / (deviation * deviation * deviation)
  }
}

///|
pub fn excess_kurtosis(data : Array[Double]) -> Double {
  let deviation = population_stddev(data)
  if deviation == 0.0 {
    0.0
  } else {
    central_moment(data, 4) / (deviation * deviation * deviation * deviation) -
    3.0
  }
}

///|
pub fn weighted_mean(data : Array[Double], weights : Array[Double]) -> Double {
  if data.length() == 0 || data.length() != weights.length() {
    return 0.0
  }
  let mut numerator = 0.0
  let mut denominator = 0.0
  for index = 0; index < data.length(); index = index + 1 {
    if weights[index] < 0.0 {
      abort("weights must be non-negative")
    }
    numerator += data[index] * weights[index]
    denominator += weights[index]
  }
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
pub fn weighted_variance(
  data : Array[Double],
  weights : Array[Double],
  sample? : Bool = false,
) -> Double {
  if data.length() == 0 || data.length() != weights.length() {
    return 0.0
  }
  let center = weighted_mean(data, weights)
  let mut total = 0.0
  let mut weight_total = 0.0
  let mut squared_weight_total = 0.0
  for index = 0; index < data.length(); index = index + 1 {
    if weights[index] < 0.0 {
      abort("weights must be non-negative")
    }
    let deviation = data[index] - center
    total += weights[index] * deviation * deviation
    weight_total += weights[index]
    squared_weight_total += weights[index] * weights[index]
  }
  if weight_total == 0.0 {
    return 0.0
  }
  if sample {
    let correction = weight_total - squared_weight_total / weight_total
    if correction <= 0.0 {
      0.0
    } else {
      total / correction
    }
  } else {
    total / weight_total
  }
}

///|
pub fn coefficient_of_variation(data : Array[Double]) -> Double {
  let center = mean(data)
  if center == 0.0 {
    0.0
  } else {
    population_stddev(data) / abs_double(center)
  }
}

///|
pub fn robust_scale(data : Array[Double]) -> Double {
  mad(data)
}

///|
pub fn median_squared_deviation(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let center = median(data)
  let deviations = []
  for value in data {
    let delta = value - center
    deviations.push(delta * delta)
  }
  median(deviations)
}

///|
pub fn gini_mean_difference(data : Array[Double]) -> Double {
  if data.length() <= 1 {
    return 0.0
  }
  let mut total = 0.0
  for left in data {
    for right in data {
      total += abs_double(left - right)
    }
  }
  total / (data.length() * (data.length() - 1)).to_double()
}

///|
pub fn root_mean_square(data : Array[Double]) -> Double {
  if data.length() == 0 {
    0.0
  } else {
    (sum_squared(data) / data.length().to_double()).sqrt()
  }
}

///|
pub fn harmonic_mean(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let mut reciprocal_sum = 0.0
  for value in data {
    if value <= 0.0 {
      abort("harmonic mean requires positive values")
    }
    reciprocal_sum += 1.0 / value
  }
  data.length().to_double() / reciprocal_sum
}

///|
pub fn geometric_mean(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let mut product = 1.0
  for value in data {
    if value <= 0.0 {
      abort("geometric mean requires positive values")
    }
    product *= value
  }
  // Newton iteration avoids requiring a logarithm dependency and is stable for
  // the moderate-sized monitoring batches this package targets.
  let mut result = if product > 1.0 { product } else { 1.0 }
  let n = data.length().to_double()
  for iteration = 0; iteration < 32; iteration = iteration + 1 {
    let denominator = result
    let mut power = 1.0
    for power_index = 1
        power_index < data.length()
        power_index = power_index + 1 {
      power *= denominator
    }
    if power == 0.0 {
      break
    }
    result = ((n - 1.0) * result + product / power) / n
  }
  result
}