///|
pub fn huber_objective(
  data : Array[Double],
  location : Double,
  tuning : Double,
) -> Double {
  let mut total = 0.0
  for value in data {
    total += huber_loss(value - location, tuning)
  }
  total
}

///|
pub fn tukey_objective(
  data : Array[Double],
  location : Double,
  tuning : Double,
) -> Double {
  if tuning <= 0.0 {
    abort("tuning must be positive")
  }
  let mut total = 0.0
  for value in data {
    let residual = abs_double(value - location)
    let scaled = residual / tuning
    if scaled < 1.0 {
      let one_minus = 1.0 - scaled * scaled
      total += tuning * tuning / 6.0 * (1.0 - one_minus * one_minus * one_minus)
    } else {
      total += tuning * tuning / 6.0
    }
  }
  total
}

///|
pub fn grid_search_location(
  data : Array[Double],
  lower : Double,
  upper : Double,
  steps : Int,
) -> Double {
  if steps <= 0 || lower > upper {
    abort("invalid grid configuration")
  }
  if data.length() == 0 {
    return 0.0
  }
  let mut best = lower
  let mut best_loss = huber_objective(data, lower, mad(data) + 0.000001)
  for index = 1; index <= steps; index = index + 1 {
    let candidate = lower +
      (upper - lower) * index.to_double() / steps.to_double()
    let loss = huber_objective(data, candidate, mad(data) + 0.000001)
    if loss < best_loss {
      best = candidate
      best_loss = loss
    }
  }
  best
}

///|
pub fn golden_section_location(
  data : Array[Double],
  lower : Double,
  upper : Double,
  iterations : Int,
) -> Double {
  if iterations <= 0 || lower > upper {
    abort("invalid search configuration")
  }
  let tuning = if mad(data) == 0.0 { 1.0 } else { mad(data) }
  let mut left = lower
  let mut right = upper
  let ratio = 0.6180339887
  for _iteration = 0; _iteration < iterations; _iteration = _iteration + 1 {
    let first = right - ratio * (right - left)
    let second = left + ratio * (right - left)
    if huber_objective(data, first, tuning) <
      huber_objective(data, second, tuning) {
      right = second
    } else {
      left = first
    }
  }
  (left + right) / 2.0
}

///|
pub fn adaptive_huber_tuning(
  data : Array[Double],
  target_fraction : Double,
) -> Double {
  if target_fraction <= 0.0 || target_fraction >= 1.0 {
    abort("target_fraction must be in (0, 1)")
  }
  let scale = mad(data)
  if scale == 0.0 {
    1.0
  } else {
    abs_double(quantile(data, 1.0 - target_fraction) - median(data)) / scale
  }
}

///|
pub fn estimate_trim_by_loss(
  data : Array[Double],
  candidates : Array[Double],
) -> Double {
  if candidates.length() == 0 {
    return 0.0
  }
  let mut best = candidates[0]
  let mut best_loss = 0.0
  let first = winsorized_mean(data, candidates[0])
  for value in data {
    best_loss += abs_double(value - first)
  }
  for index = 1; index < candidates.length(); index = index + 1 {
    let candidate = candidates[index]
    let center = winsorized_mean(data, candidate)
    let mut loss = 0.0
    for value in data {
      loss += huber_loss(value - center, mad(data) + 0.000001)
    }
    if loss < best_loss {
      best = candidate
      best_loss = loss
    }
  }
  best
}

///|
pub fn robust_location_gradient(
  data : Array[Double],
  location : Double,
  tuning : Double,
) -> Double {
  if tuning <= 0.0 {
    abort("tuning must be positive")
  }
  let mut total = 0.0
  for value in data {
    total -= if abs_double(value - location) <= tuning {
      value - location
    } else {
      tuning * sign_double(value - location)
    }
  }
  total
}

///|
pub fn gradient_descent_location(
  data : Array[Double],
  initial : Double,
  learning_rate : Double,
  iterations : Int,
  tuning : Double,
) -> Double {
  if learning_rate <= 0.0 || iterations < 0 {
    abort("invalid gradient configuration")
  }
  let mut location = initial
  for _iteration = 0; _iteration < iterations; _iteration = _iteration + 1 {
    location -= learning_rate * robust_location_gradient(data, location, tuning)
  }
  location
}

///|
pub fn robust_scale_objective(data : Array[Double], scale : Double) -> Double {
  if scale <= 0.0 {
    abort("scale must be positive")
  }
  let center = median(data)
  let mut total = 0.0
  for value in data {
    total += huber_loss(abs_double(value - center) - scale, scale)
  }
  total
}

///|
pub fn scale_grid_search(
  data : Array[Double],
  lower : Double,
  upper : Double,
  steps : Int,
) -> Double {
  if lower <= 0.0 || upper < lower || steps <= 0 {
    abort("invalid scale grid")
  }
  let mut best = lower
  let mut best_loss = robust_scale_objective(data, lower)
  for index = 1; index <= steps; index = index + 1 {
    let candidate = lower +
      (upper - lower) * index.to_double() / steps.to_double()
    let loss = robust_scale_objective(data, candidate)
    if loss < best_loss {
      best = candidate
      best_loss = loss
    }
  }
  best
}

///|
pub fn constrained_location(
  data : Array[Double],
  lower : Double,
  upper : Double,
) -> Double {
  clamp_double(huber_location(data), lower, upper)
}

///|
pub fn robust_location_path(
  data : Array[Double],
  start : Double,
  steps : Int,
  learning_rate : Double,
  tuning : Double,
) -> Array[Double] {
  if steps < 0 {
    abort("steps must not be negative")
  }
  let result = []
  let mut location = start
  result.push(location)
  for _step = 0; _step < steps; _step = _step + 1 {
    location -= learning_rate * robust_location_gradient(data, location, tuning)
    result.push(location)
  }
  result
}