///|
/// Robust additive seasonal decomposition.
pub struct SeasonalDecomposition {
  trend : Array[Double]
  seasonal : Array[Double]
  residual : Array[Double]
  period : Int
  strength : Double
}

///|
/// Configuration for seasonality detection.
pub struct SeasonalityRule {
  minimum_period : Int
  maximum_period : Int
  smoothing_window : Int
  threshold : Double
}

///|
pub fn seasonality_default_rule() -> SeasonalityRule {
  { minimum_period: 2, maximum_period: 30, smoothing_window: 5, threshold: 0.2 }
}

///|
pub fn seasonality_rule(
  minimum_period : Int,
  maximum_period : Int,
  smoothing_window : Int,
  threshold : Double,
) -> SeasonalityRule {
  let lower = if minimum_period < 2 { 2 } else { minimum_period }
  let upper = if maximum_period < lower { lower } else { maximum_period }
  {
    minimum_period: lower,
    maximum_period: upper,
    smoothing_window: if smoothing_window < 2 {
      2
    } else {
      smoothing_window
    },
    threshold: if threshold < 0.0 {
      0.0
    } else {
      threshold
    },
  }
}

///|
pub fn seasonal_trend(data : Array[Double], window : Int) -> Array[Double] {
  rolling_median(data, if window < 2 { 2 } else { window })
}

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

///|
pub fn seasonal_strength_from_components(
  residual : Array[Double],
  detrended : Array[Double],
) -> Double {
  let residual_scale = mad(residual)
  let detrended_scale = mad(detrended)
  if detrended_scale <= 1.0e-12 {
    0.0
  } else {
    1.0 - residual_scale / detrended_scale
  }
}

///|
pub fn seasonal_profile(data : Array[Double], period : Int) -> Array[Double] {
  let safe_period = if period < 2 { 2 } else { period }
  let profile = []
  for slot = 0; slot < safe_period; slot = slot + 1 {
    let values = []
    let mut index = slot
    while index < data.length() {
      values.push(data[index])
      index += safe_period
    }
    profile.push(if values.length() == 0 { 0.0 } else { median(values) })
  }
  let center = mean(profile)
  for index = 0; index < profile.length(); index = index + 1 {
    profile[index] -= center
  }
  profile
}

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

///|
pub fn seasonal_decompose(
  data : Array[Double],
  period : Int,
  window : Int,
) -> SeasonalDecomposition {
  let trend = seasonal_trend(data, window)
  let detrended = seasonal_residual(data, trend)
  let seasonal = seasonal_component(detrended, period)
  let residual = []
  let limit = if detrended.length() < seasonal.length() {
    detrended.length()
  } else {
    seasonal.length()
  }
  for index = 0; index < limit; index = index + 1 {
    residual.push(detrended[index] - seasonal[index])
  }
  {
    trend,
    seasonal,
    residual,
    period: if period < 2 {
      2
    } else {
      period
    },
    strength: seasonal_strength_from_components(residual, detrended),
  }
}

///|
pub fn seasonal_autocorrelation_at(data : Array[Double], lag : Int) -> Double {
  robust_autocorrelation(data, if lag <= 0 { 1 } else { lag })
}

///|
pub fn seasonal_candidate_periods(rule : SeasonalityRule) -> Array[Int] {
  let result = []
  for period = rule.minimum_period
      period <= rule.maximum_period
      period = period + 1 {
    result.push(period)
  }
  result
}

///|
pub fn seasonal_period_score(data : Array[Double], period : Int) -> Double {
  if period <= 1 || data.length() <= period {
    return 0.0
  }
  abs_double(seasonal_autocorrelation_at(data, period))
}

///|
pub fn seasonal_period_scores(
  data : Array[Double],
  rule : SeasonalityRule,
) -> Array[Double] {
  let result = []
  for period in seasonal_candidate_periods(rule) {
    result.push(seasonal_period_score(data, period))
  }
  result
}

///|
pub fn seasonal_best_period(
  data : Array[Double],
  rule : SeasonalityRule,
) -> Int {
  let periods = seasonal_candidate_periods(rule)
  if periods.length() == 0 {
    return rule.minimum_period
  }
  let mut best = periods[0]
  let mut score = seasonal_period_score(data, best)
  for period in periods {
    let current = seasonal_period_score(data, period)
    if current > score {
      best = period
      score = current
    }
  }
  best
}

///|
pub fn seasonal_is_present(
  data : Array[Double],
  rule : SeasonalityRule,
) -> Bool {
  seasonal_period_score(data, seasonal_best_period(data, rule)) >=
  rule.threshold
}

///|
pub fn seasonal_decompose_auto(
  data : Array[Double],
  rule : SeasonalityRule,
) -> SeasonalDecomposition {
  let period = seasonal_best_period(data, rule)
  seasonal_decompose(data, period, rule.smoothing_window)
}

///|
pub fn seasonal_remove(data : Array[Double], period : Int) -> Array[Double] {
  remove_seasonal_median(data, if period < 2 { 2 } else { period })
}

///|
pub fn seasonal_add(
  base : Array[Double],
  period : Int,
  profile : Array[Double],
) -> Array[Double] {
  let result = []
  if profile.length() == 0 {
    return base.copy()
  }
  let safe_period = if period <= 0 { 1 } else { period }
  for index = 0; index < base.length(); index = index + 1 {
    result.push(base[index] + profile[index % safe_period % profile.length()])
  }
  result
}

///|
pub fn seasonal_profile_center(profile : Array[Double]) -> Double {
  mean(profile)
}

///|
pub fn seasonal_profile_scale(profile : Array[Double]) -> Double {
  mad(profile)
}

///|
pub fn seasonal_profile_range(profile : Array[Double]) -> Double {
  range(profile)
}

///|
pub fn seasonal_profile_normalized(profile : Array[Double]) -> Array[Double] {
  robust_standardize(profile)
}

///|
pub fn seasonal_profile_difference(
  left : Array[Double],
  right : Array[Double],
) -> Array[Double] {
  let result = []
  let limit = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  for index = 0; index < limit; index = index + 1 {
    result.push(left[index] - right[index])
  }
  result
}

///|
pub fn seasonal_profile_distance(
  left : Array[Double],
  right : Array[Double],
) -> Double {
  sum_absolute(seasonal_profile_difference(left, right))
}

///|
pub fn seasonal_profile_correlation(
  left : Array[Double],
  right : Array[Double],
) -> Double {
  pearson_correlation(left, right)
}

///|
pub fn seasonal_outlier_score(
  data : Array[Double],
  period : Int,
  threshold : Double,
) -> Array[Double] {
  let component = seasonal_component(data, period)
  let residual = seasonal_residual(data, component)
  let scale = mad(residual) * 1.4826
  let result = []
  for value in residual {
    result.push(
      if scale <= 1.0e-12 {
        0.0
      } else {
        abs_double(value - median(residual)) / scale / threshold
      },
    )
  }
  result
}

///|
pub fn seasonal_outlier_flags(
  data : Array[Double],
  period : Int,
  threshold : Double,
) -> Array[Bool] {
  let result = []
  for score in seasonal_outlier_score(data, period, threshold) {
    result.push(score > 1.0)
  }
  result
}

///|
pub fn seasonal_change_indices(
  data : Array[Double],
  period : Int,
  threshold : Double,
) -> Array[Int] {
  let flags = seasonal_outlier_flags(data, period, threshold)
  let result = []
  for index = 0; index < flags.length(); index = index + 1 {
    if flags[index] {
      result.push(index)
    }
  }
  result
}

///|
pub fn seasonal_forecast(
  data : Array[Double],
  period : Int,
  horizon : Int,
) -> Array[Double] {
  let safe_period = if period < 2 { 2 } else { period }
  let count = if horizon <= 0 { 1 } else { horizon }
  let profile = seasonal_profile(data, safe_period)
  let level = median(data)
  let result = []
  for index = 0; index < count; index = index + 1 {
    if profile.length() == 0 {
      result.push(level)
    } else {
      result.push(level + profile[index % profile.length()])
    }
  }
  result
}

///|
pub fn seasonal_forecast_with_trend(
  data : Array[Double],
  period : Int,
  horizon : Int,
  window : Int,
) -> Array[Double] {
  let decomposition = seasonal_decompose(data, period, window)
  let trend_forecast = forecast_robust_drift(decomposition.trend, horizon)
  let result = []
  for index = 0; index < trend_forecast.length(); index = index + 1 {
    let seasonal_value = if decomposition.seasonal.length() == 0 {
      0.0
    } else {
      decomposition.seasonal[(data.length() + index) %
      decomposition.seasonal.length()]
    }
    result.push(trend_forecast[index] + seasonal_value)
  }
  result
}

///|
pub fn seasonal_forecast_interval(
  data : Array[Double],
  period : Int,
  horizon : Int,
  confidence : Double,
) -> Array[Array[Double]] {
  let predictions = seasonal_forecast(data, period, horizon)
  let residual = seasonal_residual(data, seasonal_component(data, period))
  let scale = mad(residual) * 1.4826
  let z = if confidence >= 0.99 {
    2.58
  } else if confidence >= 0.9 {
    1.96
  } else {
    1.64
  }
  let result = []
  for index = 0; index < predictions.length(); index = index + 1 {
    let width = z * scale * (index + 1).to_double().sqrt()
    result.push([
      predictions[index] - width,
      predictions[index],
      predictions[index] + width,
    ])
  }
  result
}

///|
pub fn seasonal_decomposition_score(
  decomposition : SeasonalDecomposition,
) -> Double {
  decomposition.strength / (1.0 + mad(decomposition.residual))
}

///|
pub fn seasonal_decomposition_vector(
  decomposition : SeasonalDecomposition,
) -> Array[Double] {
  [
    decomposition.period.to_double(),
    decomposition.strength,
    mean(decomposition.trend),
    mad(decomposition.trend),
    mean(decomposition.seasonal),
    mad(decomposition.seasonal),
    mean(decomposition.residual),
    mad(decomposition.residual),
  ]
}

///|
pub fn seasonal_decomposition_lines(
  decomposition : SeasonalDecomposition,
) -> Array[String] {
  [
    "period=" + decomposition.period.to_string(),
    "strength=" + decomposition.strength.to_string(),
    "trend_scale=" + mad(decomposition.trend).to_string(),
    "seasonal_scale=" + mad(decomposition.seasonal).to_string(),
    "residual_scale=" + mad(decomposition.residual).to_string(),
  ]
}

///|
pub fn seasonal_decomposition_string(
  decomposition : SeasonalDecomposition,
) -> String {
  seasonal_decomposition_lines(decomposition).join("\n")
}

///|
pub fn seasonal_residual_quality(
  decomposition : SeasonalDecomposition,
) -> Double {
  robust_signal_quality(decomposition.residual)
}

///|
pub fn seasonal_stability(data : Array[Double], period : Int) -> Double {
  let first = seasonal_profile(data, period)
  let shifted = seasonal_profile(seasonal_remove(data, period), period)
  1.0 / (1.0 + seasonal_profile_distance(first, shifted))
}

///|
pub fn seasonal_period_consensus(
  data : Array[Double],
  periods : Array[Int],
) -> Int {
  if periods.length() == 0 {
    return 0
  }
  let scores = []
  for period in periods {
    scores.push(seasonal_period_score(data, period))
  }
  let mut best = 0
  let mut value = scores[0]
  for index = 1; index < scores.length(); index = index + 1 {
    if scores[index] > value {
      best = index
      value = scores[index]
    }
  }
  periods[best]
}

///|
pub fn seasonal_segment_strength(
  data : Array[Double],
  segments : Int,
  period : Int,
) -> Array[Double] {
  let result = []
  for values in segment_values(data, segments) {
    result.push(seasonal_decompose(values, period, period).strength)
  }
  result
}

///|
pub fn seasonal_strength_change(
  data : Array[Double],
  segments : Int,
  period : Int,
) -> Array[Double] {
  let strengths = seasonal_segment_strength(data, segments, period)
  let result = []
  for index = 1; index < strengths.length(); index = index + 1 {
    result.push(strengths[index] - strengths[index - 1])
  }
  result
}