///|
pub struct OutlierReport {
  lower_fence : Double
  upper_fence : Double
  lower_count : Int
  upper_count : Int
  indices : Array[Int]
  scores : Array[Double]
}

///|
pub fn iqr_fences(
  data : Array[Double],
  multiplier? : Double = 1.5,
) -> Array[Double] {
  if multiplier < 0.0 {
    abort("multiplier must be non-negative")
  }
  let iqr = interquartile_range(data)
  [
    lower_quartile(data) - multiplier * iqr,
    upper_quartile(data) + multiplier * iqr,
  ]
}

///|
pub fn outlier_indices_iqr(
  data : Array[Double],
  multiplier? : Double = 1.5,
) -> Array[Int] {
  let fences = iqr_fences(data, multiplier~)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] < fences[0] || data[index] > fences[1] {
      result.push(index)
    }
  }
  result
}

///|
pub fn robust_z_score(
  value : Double,
  center : Double,
  scale : Double,
) -> Double {
  if scale == 0.0 {
    0.0
  } else {
    (value - center) / scale
  }
}

///|
pub fn robust_z_scores(data : Array[Double]) -> Array[Double] {
  let center = median(data)
  let scale = mad(data)
  let result = []
  for value in data {
    result.push(robust_z_score(value, center, scale))
  }
  result
}

///|
pub fn outlier_indices_z(
  data : Array[Double],
  threshold? : Double = 3.5,
) -> Array[Int] {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  let scores = robust_z_scores(data)
  let result = []
  for index = 0; index < scores.length(); index = index + 1 {
    if abs_double(scores[index]) > threshold {
      result.push(index)
    }
  }
  result
}

///|
pub fn summarize_outliers(
  data : Array[Double],
  threshold? : Double = 3.5,
) -> OutlierReport {
  let fences = iqr_fences(data)
  let scores = robust_z_scores(data)
  let indices = []
  let mut lower_count = 0
  let mut upper_count = 0
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] < fences[0] {
      indices.push(index)
      lower_count += 1
    } else if data[index] > fences[1] {
      indices.push(index)
      upper_count += 1
    } else if abs_double(scores[index]) > threshold {
      indices.push(index)
      if data[index] < median(data) {
        lower_count += 1
      } else {
        upper_count += 1
      }
    }
  }
  {
    lower_fence: fences[0],
    upper_fence: fences[1],
    lower_count,
    upper_count,
    indices,
    scores,
  }
}

///|
pub fn hampel_scores(data : Array[Double], window : Int) -> Array[Double] {
  if window < 1 {
    abort("window must be positive")
  }
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    let start = if index - window < 0 { 0 } else { index - window }
    let end = if index + window + 1 > data.length() {
      data.length()
    } else {
      index + window + 1
    }
    let neighborhood = []
    for position = start; position < end; position = position + 1 {
      neighborhood.push(data[position])
    }
    let center = median(neighborhood)
    let scale = mad(neighborhood)
    if scale == 0.0 {
      result.push(if data[index] == center { 0.0 } else { 1000000000.0 })
    } else {
      result.push(robust_z_score(data[index], center, scale))
    }
  }
  result
}

///|
pub fn hampel_filter(
  data : Array[Double],
  window : Int,
  threshold? : Double = 3.5,
) -> Array[Double] {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  let scores = hampel_scores(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    if abs_double(scores[index]) > threshold {
      let start = if index - window < 0 { 0 } else { index - window }
      let end = if index + window + 1 > data.length() {
        data.length()
      } else {
        index + window + 1
      }
      let neighborhood = []
      for position = start; position < end; position = position + 1 {
        if position != index {
          neighborhood.push(data[position])
        }
      }
      result.push(median(neighborhood))
    } else {
      result.push(data[index])
    }
  }
  result
}

///|
pub fn winsorize_by_z(
  data : Array[Double],
  threshold : Double,
) -> Array[Double] {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  let center = median(data)
  let scale = mad(data)
  let result = []
  for value in data {
    if scale == 0.0 {
      result.push(value)
    } else {
      let limit = threshold * scale
      result.push(clamp_double(value, center - limit, center + limit))
    }
  }
  result
}

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

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

///|
pub fn clip_by_quantiles(
  data : Array[Double],
  lower_probability : Double,
  upper_probability : Double,
) -> Array[Double] {
  if lower_probability > upper_probability {
    abort("lower probability must not exceed upper probability")
  }
  let lower = quantile(data, lower_probability)
  let upper = quantile(data, upper_probability)
  let result = []
  for value in data {
    result.push(clamp_double(value, lower, upper))
  }
  result
}

///|
pub fn outlier_fraction(indices : Array[Int], sample_size : Int) -> Double {
  if sample_size <= 0 {
    0.0
  } else {
    indices.length().to_double() / sample_size.to_double()
  }
}