///|
pub fn seasonal_medians(data : Array[Double], period : Int) -> Array[Double] {
  if period <= 0 {
    abort("period must be positive")
  }
  let buckets = []
  for season = 0; season < period; season = season + 1 {
    buckets.push([])
  }
  for index = 0; index < data.length(); index = index + 1 {
    buckets[index % period].push(data[index])
  }
  let result = []
  for bucket in buckets {
    result.push(median(bucket))
  }
  result
}

///|
pub fn remove_seasonal_median(
  data : Array[Double],
  period : Int,
) -> Array[Double] {
  let seasonal = seasonal_medians(data, period)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(data[index] - seasonal[index % period])
  }
  result
}

///|
pub fn add_seasonal_median(
  data : Array[Double],
  period : Int,
  seasonal : Array[Double],
) -> Array[Double] {
  if period <= 0 || seasonal.length() != period {
    abort("seasonal dimensions are invalid")
  }
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(data[index] + seasonal[index % period])
  }
  result
}

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

///|
pub fn robust_moving_average(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  rolling_winsorized_mean(data, window, 0.1)
}

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

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

///|
pub fn robust_moving_signal_to_noise(
  data : Array[Double],
  window : Int,
) -> Array[Double] {
  let center = rolling_median(data, window)
  let scale = rolling_mad(data, window)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(
      if scale[index] == 0.0 {
        0.0
      } else {
        center[index] / scale[index]
      },
    )
  }
  result
}

///|
pub fn cumulative_sum(data : Array[Double]) -> Array[Double] {
  let result = []
  let mut total = 0.0
  for value in data {
    total += value
    result.push(total)
  }
  result
}

///|
pub fn cumulative_robust_sum(
  data : Array[Double],
  clip : Double,
) -> Array[Double] {
  if clip <= 0.0 {
    abort("clip must be positive")
  }
  let result = []
  let mut total = 0.0
  for value in data {
    total += clamp_double(value, -clip, clip)
    result.push(total)
  }
  result
}

///|
pub fn cumulative_median(data : Array[Double]) -> Array[Double] {
  let result = []
  let prefix = []
  for value in data {
    prefix.push(value)
    result.push(median(prefix))
  }
  result
}

///|
pub fn cumulative_mad(data : Array[Double]) -> Array[Double] {
  let result = []
  let prefix = []
  for value in data {
    prefix.push(value)
    result.push(mad(prefix))
  }
  result
}

///|
pub fn cumulative_outlier_rate(
  data : Array[Double],
  threshold? : Double = 3.5,
) -> Array[Double] {
  let result = []
  let prefix = []
  for value in data {
    prefix.push(value)
    result.push(
      outlier_fraction(outlier_indices_z(prefix, threshold~), prefix.length()),
    )
  }
  result
}

///|
pub fn robust_drift_score(data : Array[Double], window : Int) -> Array[Double] {
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    let current = extract_window(data, index, window)
    let previous_start = if index - window < 0 { 0 } else { index - window }
    let previous = []
    for position = previous_start; position < index; position = position + 1 {
      previous.push(data[position])
    }
    if previous.length() == 0 {
      result.push(0.0)
    } else {
      let scale = mad(previous)
      result.push(
        if scale == 0.0 {
          median(current) - median(previous)
        } else {
          (median(current) - median(previous)) / scale
        },
      )
    }
  }
  result
}

///|
pub fn drift_indices(
  data : Array[Double],
  window : Int,
  threshold : Double,
) -> Array[Int] {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  let scores = robust_drift_score(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_exponential_scale(
  data : Array[Double],
  alpha : Double,
) -> Array[Double] {
  let base_scale = mad(data)
  let clip = 3.5 * (if base_scale == 0.0 { 1.0 } else { base_scale })
  let centers = robust_exponentially_weighted_mean(data, alpha, clip)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(abs_double(data[index] - centers[index]))
  }
  result
}

///|
pub fn rolling_quantile_band(
  data : Array[Double],
  window : Int,
  probability : Double,
) -> Array[Array[Double]] {
  if probability <= 0.0 || probability >= 0.5 {
    abort("probability must be in (0, 0.5)")
  }
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    let values = extract_window(data, index, window)
    result.push([
      quantile(values, probability),
      quantile(values, 1.0 - probability),
    ])
  }
  result
}

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

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

///|
pub fn robust_growth_rate(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 {
    let base = abs_double(data[index - 1])
    result.push(
      if base == 0.0 {
        0.0
      } else {
        (data[index] - data[index - 1]) / base
      },
    )
  }
  result
}