///|
pub fn nearest_center(value : Double, centers : Array[Double]) -> Int {
  if centers.length() == 0 {
    abort("centers must not be empty")
  }
  let mut best = 0
  let mut distance = abs_double(value - centers[0])
  for index = 1; index < centers.length(); index = index + 1 {
    let current = abs_double(value - centers[index])
    if current < distance {
      distance = current
      best = index
    }
  }
  best
}

///|
pub fn one_dimensional_cluster_assignments(
  data : Array[Double],
  centers : Array[Double],
) -> Array[Int] {
  let result = []
  for value in data {
    result.push(nearest_center(value, centers))
  }
  result
}

///|
pub fn cluster_members(
  data : Array[Double],
  assignments : Array[Int],
  cluster : Int,
) -> Array[Double] {
  let result = []
  for index = 0
      index < data.length() && index < assignments.length()
      index = index + 1 {
    if assignments[index] == cluster {
      result.push(data[index])
    }
  }
  result
}

///|
pub fn cluster_medians(
  data : Array[Double],
  assignments : Array[Int],
  cluster_count : Int,
) -> Array[Double] {
  if cluster_count <= 0 {
    abort("cluster_count must be positive")
  }
  let result = []
  for cluster = 0; cluster < cluster_count; cluster = cluster + 1 {
    result.push(median(cluster_members(data, assignments, cluster)))
  }
  result
}

///|
pub fn initialize_quantile_centers(
  data : Array[Double],
  cluster_count : Int,
) -> Array[Double] {
  if cluster_count <= 0 {
    abort("cluster_count must be positive")
  }
  let result = []
  for cluster = 0; cluster < cluster_count; cluster = cluster + 1 {
    let probability = (cluster.to_double() + 0.5) / cluster_count.to_double()
    result.push(quantile(data, probability))
  }
  result
}

///|
pub fn k_medians_1d(
  data : Array[Double],
  cluster_count : Int,
  max_iter? : Int = 50,
  tol? : Double = 0.0001,
) -> Array[Double] {
  if data.length() == 0 {
    return []
  }
  if cluster_count <= 0 || max_iter <= 0 || tol <= 0.0 {
    abort("invalid clustering configuration")
  }
  let mut centers = initialize_quantile_centers(data, cluster_count)
  for _iteration = 0; _iteration < max_iter; _iteration = _iteration + 1 {
    let assignments = one_dimensional_cluster_assignments(data, centers)
    let next = cluster_medians(data, assignments, cluster_count)
    let mut movement = 0.0
    for index = 0; index < centers.length(); index = index + 1 {
      movement += abs_double(next[index] - centers[index])
    }
    centers = next
    if movement <= tol {
      break
    }
  }
  centers
}

///|
pub fn k_medians_loss(data : Array[Double], centers : Array[Double]) -> Double {
  let mut total = 0.0
  for value in data {
    total += abs_double(value - centers[nearest_center(value, centers)])
  }
  total
}

///|
pub fn k_medians_assign(
  data : Array[Double],
  cluster_count : Int,
) -> Array[Int] {
  let centers = k_medians_1d(data, cluster_count)
  one_dimensional_cluster_assignments(data, centers)
}

///|
pub fn robust_cluster_outliers(
  data : Array[Double],
  cluster_count : Int,
  threshold : Double,
) -> Array[Int] {
  let centers = k_medians_1d(data, cluster_count)
  let assignments = one_dimensional_cluster_assignments(data, centers)
  let residuals = []
  for index = 0; index < data.length(); index = index + 1 {
    residuals.push(data[index] - centers[assignments[index]])
  }
  outlier_indices_z(residuals, threshold~)
}

///|
pub fn cluster_spread(
  data : Array[Double],
  centers : Array[Double],
) -> Array[Double] {
  let assignments = one_dimensional_cluster_assignments(data, centers)
  let result = []
  for cluster = 0; cluster < centers.length(); cluster = cluster + 1 {
    let values = cluster_members(data, assignments, cluster)
    result.push(mad(values))
  }
  result
}

///|
pub fn cluster_counts(
  assignments : Array[Int],
  cluster_count : Int,
) -> Array[Int] {
  if cluster_count <= 0 {
    abort("cluster_count must be positive")
  }
  let result = []
  for cluster = 0; cluster < cluster_count; cluster = cluster + 1 {
    result.push(0)
  }
  for assignment in assignments {
    if assignment >= 0 && assignment < cluster_count {
      result[assignment] += 1
    }
  }
  result
}

///|
pub fn robust_cluster_summary(
  data : Array[Double],
  cluster_count : Int,
) -> Array[Array[Double]] {
  let centers = k_medians_1d(data, cluster_count)
  let assignments = one_dimensional_cluster_assignments(data, centers)
  let spreads = cluster_spread(data, centers)
  let counts = cluster_counts(assignments, cluster_count)
  let result = []
  for cluster = 0; cluster < cluster_count; cluster = cluster + 1 {
    result.push([
      centers[cluster],
      spreads[cluster],
      counts[cluster].to_double(),
    ])
  }
  result
}

///|
pub fn robust_partition_by_median(data : Array[Double]) -> Array[Int] {
  let center = median(data)
  let result = []
  for value in data {
    result.push(if value <= center { 0 } else { 1 })
  }
  result
}

///|
pub fn cluster_separation(
  centers : Array[Double],
  spreads : Array[Double],
) -> Double {
  if centers.length() < 2 {
    return 0.0
  }
  let mut minimum = abs_double(centers[0] - centers[1])
  for left = 0; left < centers.length(); left = left + 1 {
    for right = left + 1; right < centers.length(); right = right + 1 {
      let distance = abs_double(centers[left] - centers[right])
      let scale = spreads[left] + spreads[right]
      let normalized = if scale == 0.0 { distance } else { distance / scale }
      if normalized < minimum {
        minimum = normalized
      }
    }
  }
  minimum
}

///|
pub fn cluster_stability(data : Array[Double], cluster_count : Int) -> Double {
  if data.length() < 2 {
    return 1.0
  }
  let first = k_medians_1d(data, cluster_count)
  let reversed = []
  for index = data.length() - 1; index >= 0; index = index - 1 {
    reversed.push(data[index])
  }
  let second = k_medians_1d(reversed, cluster_count)
  1.0 / (1.0 + k_medians_loss(first, second))
}