///|
pub struct LocationScaleEstimate {
  location : Double
  scale : Double
  iterations : Int
  converged : Bool
}

///|
pub fn median_of_means(data : Array[Double], groups : Int) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  if groups <= 0 {
    abort("groups must be positive")
  }
  let actual_groups = if groups > data.length() {
    data.length()
  } else {
    groups
  }
  let means = []
  for group = 0; group < actual_groups; group = group + 1 {
    let start = group * data.length() / actual_groups
    let end = (group + 1) * data.length() / actual_groups
    let block = []
    for index = start; index < end; index = index + 1 {
      block.push(data[index])
    }
    means.push(mean(block))
  }
  median(means)
}

///|
pub fn huber_loss(residual : Double, tuning : Double) -> Double {
  if tuning <= 0.0 {
    abort("tuning must be positive")
  }
  let magnitude = abs_double(residual)
  if magnitude <= tuning {
    0.5 * residual * residual
  } else {
    tuning * (magnitude - 0.5 * tuning)
  }
}

///|
pub fn huber_location_estimate(
  data : Array[Double],
  tuning? : Double = 1.345,
  max_iter? : Int = 50,
  tol? : Double = 0.0001,
) -> LocationScaleEstimate {
  if data.length() == 0 {
    return { location: 0.0, scale: 0.0, iterations: 0, converged: true }
  }
  if tuning <= 0.0 || max_iter < 1 || tol <= 0.0 {
    abort("invalid Huber configuration")
  }
  let mut location = median(data)
  let scale = mad(data)
  if scale == 0.0 {
    return { location, scale, iterations: 0, converged: true }
  }
  let cutoff = tuning * scale
  let mut iterations = 0
  let mut converged = false
  for iteration = 0; iteration < max_iter; iteration = iteration + 1 {
    let mut weighted_total = 0.0
    let mut weight_sum = 0.0
    for value in data {
      let weight = huber_weight(value - location, cutoff)
      weighted_total += weight * value
      weight_sum += weight
    }
    let next = if weight_sum == 0.0 {
      location
    } else {
      weighted_total / weight_sum
    }
    iterations = iteration + 1
    if abs_double(next - location) <= tol {
      location = next
      converged = true
      break
    }
    location = next
  }
  { location, scale, iterations, converged }
}

///|
pub fn tukey_bisquare_location(
  data : Array[Double],
  tuning? : Double = 4.685,
  max_iter? : Int = 50,
  tol? : Double = 0.0001,
) -> LocationScaleEstimate {
  if data.length() == 0 {
    return { location: 0.0, scale: 0.0, iterations: 0, converged: true }
  }
  if tuning <= 0.0 || max_iter < 1 || tol <= 0.0 {
    abort("invalid Tukey configuration")
  }
  let mut location = median(data)
  let scale = mad(data)
  if scale == 0.0 {
    return { location, scale, iterations: 0, converged: true }
  }
  let cutoff = tuning * scale
  let mut iterations = 0
  let mut converged = false
  for iteration = 0; iteration < max_iter; iteration = iteration + 1 {
    let mut weighted_total = 0.0
    let mut weight_sum = 0.0
    for value in data {
      let weight = tukey_bisquare_weight(value - location, cutoff)
      weighted_total += weight * value
      weight_sum += weight
    }
    let next = if weight_sum == 0.0 {
      location
    } else {
      weighted_total / weight_sum
    }
    iterations = iteration + 1
    if abs_double(next - location) <= tol {
      location = next
      converged = true
      break
    }
    location = next
  }
  { location, scale, iterations, converged }
}

///|
pub fn biweight_midvariance(
  data : Array[Double],
  tuning? : Double = 9.0,
) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  if tuning <= 0.0 {
    abort("tuning must be positive")
  }
  let center = median(data)
  let scale = mad(data)
  if scale == 0.0 {
    return 0.0
  }
  let cutoff = tuning * scale
  let mut numerator = 0.0
  let mut denominator = 0.0
  for value in data {
    let scaled = (value - center) / cutoff
    if abs_double(scaled) < 1.0 {
      let one_minus = 1.0 - scaled * scaled
      numerator += (value - center) *
        (value - center) *
        one_minus *
        one_minus *
        one_minus *
        one_minus
      denominator += one_minus * (1.0 - 5.0 * scaled * scaled)
    }
  }
  if denominator == 0.0 {
    0.0
  } else {
    (data.length().to_double() * numerator).sqrt() / abs_double(denominator)
  }
}

///|
pub fn trimmed_location_scale(
  data : Array[Double],
  trim_percent : Double,
) -> LocationScaleEstimate {
  validate_trim(trim_percent)
  if data.length() == 0 {
    return { location: 0.0, scale: 0.0, iterations: 0, converged: true }
  }
  let sorted = copy_and_sort(data)
  let cut = (sorted.length().to_double() * trim_percent).to_int()
  let retained = []
  for index = cut; index < sorted.length() - cut; index = index + 1 {
    retained.push(sorted[index])
  }
  {
    location: mean(retained),
    scale: sample_stddev(retained),
    iterations: 1,
    converged: true,
  }
}

///|
pub fn shorth_location(
  data : Array[Double],
  fraction? : Double = 0.5,
) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  if fraction <= 0.0 || fraction > 1.0 {
    abort("fraction must be in (0, 1]")
  }
  let sorted = copy_and_sort(data)
  let width = (sorted.length().to_double() * fraction).to_int()
  let window = if width < 1 { 1 } else { width }
  let mut best_start = 0
  let mut best_width = sorted[window - 1] - sorted[0]
  for start = 1; start + window <= sorted.length(); start = start + 1 {
    let current_width = sorted[start + window - 1] - sorted[start]
    if current_width < best_width {
      best_width = current_width
      best_start = start
    }
  }
  let selected = []
  for index = best_start; index < best_start + window; index = index + 1 {
    selected.push(sorted[index])
  }
  mean(selected)
}

///|
pub fn hodges_lehmann(data : Array[Double]) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let pairwise = []
  for left in data {
    for right in data {
      pairwise.push((left + right) / 2.0)
    }
  }
  median(pairwise)
}

///|
pub fn winsorized_location_scale(
  data : Array[Double],
  trim_percent : Double,
) -> LocationScaleEstimate {
  let transformed = winsorize(data, trim_percent)
  {
    location: mean(transformed),
    scale: sample_stddev(transformed),
    iterations: 1,
    converged: true,
  }
}

///|
pub fn l1_location(
  data : Array[Double],
  max_iter? : Int = 100,
  tol? : Double = 0.0001,
) -> Double {
  if max_iter < 1 || tol <= 0.0 {
    abort("invalid L1 configuration")
  }
  if data.length() == 0 {
    return 0.0
  }
  let mut location = mean(data)
  for _iteration = 0; _iteration < max_iter; _iteration = _iteration + 1 {
    let mut numerator = 0.0
    let mut denominator = 0.0
    for value in data {
      let distance = abs_double(value - location)
      let weight = if distance < tol { 1.0 / tol } else { 1.0 / distance }
      numerator += weight * value
      denominator += weight
    }
    let next = if denominator == 0.0 {
      location
    } else {
      numerator / denominator
    }
    if abs_double(next - location) <= tol {
      return next
    }
    location = next
  }
  location
}

///|
pub fn quantile_location(data : Array[Double], probability : Double) -> Double {
  quantile(data, probability)
}

///|
pub fn robust_scale_from_iqr(data : Array[Double]) -> Double {
  interquartile_range(data) / 1.3489795
}

///|
pub fn robust_scale_from_qn(data : Array[Double]) -> Double {
  if data.length() <= 1 {
    return 0.0
  }
  let distances = []
  for left = 0; left < data.length(); left = left + 1 {
    for right = left + 1; right < data.length(); right = right + 1 {
      distances.push(abs_double(data[left] - data[right]))
    }
  }
  2.2219 * quantile(distances, 0.25)
}