///|
pub fn median_filter(data : Array[Double], window : Int) -> Array[Double] {
  rolling_median(data, window)
}

///|
pub fn winsorized_filter(
  data : Array[Double],
  window : Int,
  trim_percent : Double,
) -> Array[Double] {
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    let values = extract_window(data, index, window)
    result.push(median(winsorize(values, trim_percent)))
  }
  result
}

///|
pub fn hampel_filter_series(
  data : Array[Double],
  window : Int,
  threshold? : Double = 3.5,
) -> Array[Double] {
  hampel_filter(data, window, threshold~)
}

///|
pub fn rolling_hampel_scores(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  hampel_scores(data, window)
}

///|
pub fn robust_smooth(data : Array[Double], window : Int) -> Array[Double] {
  rolling_median(hampel_filter(data, window), window)
}

///|
pub fn robust_residuals(data : Array[Double], window : Int) -> Array[Double] {
  let baseline = robust_smooth(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(data[index] - baseline[index])
  }
  result
}

///|
pub fn total_variation(data : Array[Double]) -> Double {
  if data.length() <= 1 {
    return 0.0
  }
  let mut total = 0.0
  for index = 1; index < data.length(); index = index + 1 {
    total += abs_double(data[index] - data[index - 1])
  }
  total
}

///|
pub fn robust_total_variation(data : Array[Double]) -> Double {
  total_variation(hampel_filter(data, 3))
}

///|
pub fn first_difference(data : Array[Double]) -> Array[Double] {
  let result = []
  if data.length() == 0 {
    return result
  }
  result.push(0.0)
  for index = 1; index < data.length(); index = index + 1 {
    result.push(data[index] - data[index - 1])
  }
  result
}

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

///|
pub fn robust_change_scores(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  rolling_z_scores(first_difference(data), window)
}

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

///|
pub fn robust_envelope(
  data : Array[Double],
  window : Int,
  multiplier? : Double = 3.5,
) -> Array[Array[Double]] {
  if multiplier <= 0.0 {
    abort("multiplier must be positive")
  }
  let centers = rolling_median(data, window)
  let scales = rolling_mad(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push([
      centers[index] - multiplier * scales[index],
      centers[index] + multiplier * scales[index],
    ])
  }
  result
}

///|
pub fn envelope_violations(
  data : Array[Double],
  window : Int,
  multiplier? : Double = 3.5,
) -> Array[Int] {
  let envelope = robust_envelope(data, window, multiplier~)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] < envelope[index][0] || data[index] > envelope[index][1] {
      result.push(index)
    }
  }
  result
}

///|
pub fn winsorized_residuals(
  data : Array[Double],
  trim_percent : Double,
) -> Array[Double] {
  let transformed = winsorize(data, trim_percent)
  let center = mean(transformed)
  let result = []
  for value in data {
    result.push(value - center)
  }
  result
}

///|
pub fn robust_autocorrelation(data : Array[Double], lag : Int) -> Double {
  if lag < 0 {
    abort("lag must not be negative")
  }
  if lag >= data.length() || data.length() == 0 {
    return 0.0
  }
  let first = []
  let second = []
  for index = lag; index < data.length(); index = index + 1 {
    first.push(data[index])
    second.push(data[index - lag])
  }
  pearson_correlation(first, second)
}

///|
pub fn autocorrelation_series(
  data : Array[Double],
  max_lag : Int,
) -> Array[Double] {
  if max_lag < 0 {
    abort("max_lag must not be negative")
  }
  let result = []
  let limit = if max_lag >= data.length() { data.length() - 1 } else { max_lag }
  for lag = 0; lag <= limit; lag = lag + 1 {
    result.push(robust_autocorrelation(data, lag))
  }
  result
}

///|
pub fn robust_forecast_next(data : Array[Double], window : Int) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let values = extract_window(data, data.length() - 1, window)
  let model_x = []
  for index = 0; index < values.length(); index = index + 1 {
    model_x.push(index.to_double())
  }
  let model = theil_sen_regression(model_x, values)
  model.intercept + model.slope * values.length().to_double()
}

///|
pub fn rolling_forecast_errors(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    if index < window {
      result.push(0.0)
    } else {
      let history = []
      for position = index - window; position < index; position = position + 1 {
        history.push(data[position])
      }
      let model_x = []
      for position = 0; position < history.length(); position = position + 1 {
        model_x.push(position.to_double())
      }
      let model = theil_sen_regression(model_x, history)
      let forecast = model.intercept +
        model.slope * history.length().to_double()
      result.push(data[index] - forecast)
    }
  }
  result
}

///|
pub fn robust_forecast_outliers(
  data : Array[Double],
  window : Int,
  threshold? : Double = 3.5,
) -> Array[Int] {
  outlier_indices_z(rolling_forecast_errors(data, window), threshold~)
}

///|
pub fn rolling_median_residuals(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  let centers = rolling_median(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(data[index] - centers[index])
  }
  result
}

///|
pub fn bounded_slope(data : Array[Double], maximum : Double) -> Array[Double] {
  if maximum <= 0.0 {
    abort("maximum must be positive")
  }
  let slopes = centered_difference(data)
  let result = []
  for slope in slopes {
    result.push(clamp_double(slope, -maximum, maximum))
  }
  result
}

///|
pub fn robust_interpolate(data : Array[Double], window : Int) -> Array[Double] {
  let filtered = hampel_filter(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] == filtered[index] {
      result.push(data[index])
    } else {
      result.push(filtered[index])
    }
  }
  result
}