///|
/// Forecasting methods supported by the production forecaster.
pub enum ForecastMethod {
  Naive
  Mean
  Median
  Drift
  Seasonal(Int)
  Ewma(Double)
  Huber
}

///|
/// A point forecast with an uncertainty band.
pub struct ForecastPoint {
  horizon : Int
  value : Double
  lower : Double
  upper : Double
}

///|
/// A report comparing predictions with observed values.
pub struct ForecastReport {
  predictions : Array[Double]
  residuals : Array[Double]
  mae : Double
  rmse : Double
  median_absolute_error : Double
  bias : Double
  coverage : Double
  interval_width : Double
}

///|
/// Stateful forecaster for repeated online predictions.
pub struct RobustForecaster {
  strategy : ForecastMethod
  horizon : Int
  confidence : Double
  mut history : Array[Double]
  mut last_forecast : Double
}

///|
fn forecast_safe_horizon(horizon : Int) -> Int {
  if horizon <= 0 {
    1
  } else {
    horizon
  }
}

///|
fn forecast_safe_confidence(confidence : Double) -> Double {
  if confidence <= 0.0 {
    0.5
  } else if confidence >= 1.0 {
    0.999
  } else {
    confidence
  }
}

///|
fn forecast_residual_scale(residuals : Array[Double]) -> Double {
  let scale = mad(residuals) * 1.4826
  if scale > 0.0 {
    scale
  } else {
    sample_stddev(residuals)
  }
}

///|
fn forecast_mean_or_zero(data : Array[Double]) -> Double {
  if data.length() == 0 {
    0.0
  } else {
    mean(data)
  }
}

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

///|
pub fn forecast_naive(data : Array[Double], horizon : Int) -> Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let value = if data.length() == 0 { 0.0 } else { data[data.length() - 1] }
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let value = forecast_mean_or_zero(data)
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let value = if data.length() == 0 { 0.0 } else { median(data) }
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let value = if data.length() == 0 { 0.0 } else { huber_location(data) }
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  if data.length() == 0 {
    return forecast_naive(data, count)
  }
  if data.length() == 1 {
    return forecast_naive(data, count)
  }
  let first = data[0]
  let last = data[data.length() - 1]
  let slope = (last - first) / (data.length() - 1).to_double()
  let result = []
  for step in 1..<=count {
    result.push(last + slope * step.to_double())
  }
  result
}

///|
pub fn forecast_robust_drift(
  data : Array[Double],
  horizon : Int,
) -> Array[Double] {
  let count = forecast_safe_horizon(horizon)
  if data.length() < 2 {
    return forecast_naive(data, count)
  }
  let x = []
  let y = []
  for index = 0; index < data.length(); index = index + 1 {
    x.push(index.to_double())
    y.push(data[index])
  }
  let fit = theil_sen_regression(x, y)
  let result = []
  let last_index = (data.length() - 1).to_double()
  for step in 1..<=count {
    result.push(fit.intercept + fit.slope * (last_index + step.to_double()))
  }
  result
}

///|
pub fn forecast_seasonal(
  data : Array[Double],
  period : Int,
  horizon : Int,
) -> Array[Double] {
  let count = forecast_safe_horizon(horizon)
  if period <= 0 || data.length() == 0 {
    return forecast_naive(data, count)
  }
  let result = []
  for step in 0..= 0 && source < data.length() {
      result.push(data[source])
    } else {
      result.push(data[data.length() - 1])
    }
  }
  result
}

///|
pub fn forecast_seasonal_median(
  data : Array[Double],
  period : Int,
  horizon : Int,
) -> Array[Double] {
  let count = forecast_safe_horizon(horizon)
  if period <= 0 || data.length() == 0 {
    return forecast_median_level(data, count)
  }
  let seasonal = []
  for slot in 0..= 0 && index < data.length() {
      values.push(data[index])
      index += period
    }
    seasonal.push(
      if values.length() == 0 {
        median(data)
      } else {
        median(values)
      },
    )
  }
  let result = []
  for step in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let weight = if alpha <= 0.0 {
    0.01
  } else if alpha > 1.0 {
    1.0
  } else {
    alpha
  }
  let level = if data.length() == 0 {
    0.0
  } else {
    let mut current = data[0]
    for index = 1; index < data.length(); index = index + 1 {
      current = weight * data[index] + (1.0 - weight) * current
    }
    current
  }
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let weight = if alpha <= 0.0 {
    0.01
  } else if alpha > 1.0 {
    1.0
  } else {
    alpha
  }
  if data.length() == 0 {
    return forecast_naive(data, count)
  }
  let scale = mad(data)
  let limit = if scale <= 0.0 { threshold } else { threshold * scale * 1.4826 }
  let mut level = data[0]
  for index = 1; index < data.length(); index = index + 1 {
    let value = data[index]
    let bounded = if abs_double(value - level) > limit {
      if value > level {
        level + limit
      } else {
        level - limit
      }
    } else {
      value
    }
    level = weight * bounded + (1.0 - weight) * level
  }
  let result = []
  for _ in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  let value = if data.length() == 0 { 0.0 } else { quantile(data, probability) }
  let result = []
  for _ in 0.. Array[Double] {
  forecast_quantile(
    data,
    if probability < 0.0 {
      0.0
    } else {
      probability
    },
    horizon,
  )
}

///|
pub fn forecast_upper_quantile(
  data : Array[Double],
  probability : Double,
  horizon : Int,
) -> Array[Double] {
  forecast_quantile(
    data,
    if probability > 1.0 {
      1.0
    } else {
      probability
    },
    horizon,
  )
}

///|
pub fn forecast_method_name(strategy : ForecastMethod) -> String {
  match strategy {
    Naive => "naive"
    Mean => "mean"
    Median => "median"
    Drift => "drift"
    Seasonal(period) => "seasonal-" + period.to_string()
    Ewma(alpha) => "ewma-" + alpha.to_string()
    Huber => "huber"
  }
}

///|
pub fn forecast_method_catalog() -> Array[ForecastMethod] {
  [Naive, Mean, Median, Drift, Seasonal(7), Ewma(0.2), Huber]
}

///|
pub fn forecast_with_method(
  data : Array[Double],
  strategy : ForecastMethod,
  horizon : Int,
) -> Array[Double] {
  match strategy {
    Naive => forecast_naive(data, horizon)
    Mean => forecast_mean_level(data, horizon)
    Median => forecast_median_level(data, horizon)
    Drift => forecast_robust_drift(data, horizon)
    Seasonal(period) => forecast_seasonal_median(data, period, horizon)
    Ewma(alpha) => forecast_robust_ewma(data, alpha, horizon)
    Huber => forecast_huber_level(data, horizon)
  }
}

///|
pub fn forecast_point(
  data : Array[Double],
  strategy : ForecastMethod,
  horizon : Int,
  confidence? : Double = 0.95,
) -> ForecastPoint {
  let safe_horizon = forecast_safe_horizon(horizon)
  let predictions = forecast_with_method(data, strategy, safe_horizon)
  let value = predictions[0]
  let residuals = []
  if data.length() > 1 {
    let fitted = forecast_with_method(
      forecast_prefix(data, data.length() - 1),
      strategy,
      1,
    )
    for index = 1; index < data.length(); index = index + 1 {
      residuals.push(data[index] - fitted[0])
      if index < data.length() - 1 {
        let next = forecast_with_method(
          forecast_prefix(data, index),
          strategy,
          1,
        )
        ignore(next)
      }
    }
  }
  let scale = forecast_residual_scale(residuals)
  let confidence_value = forecast_safe_confidence(confidence)
  let z = if confidence_value >= 0.99 {
    2.58
  } else if confidence_value >= 0.9 {
    1.96
  } else {
    1.64
  }
  {
    horizon: safe_horizon,
    value,
    lower: value - z * scale,
    upper: value + z * scale,
  }
}

///|
pub fn forecast_points(
  data : Array[Double],
  strategy : ForecastMethod,
  horizon : Int,
  confidence? : Double = 0.95,
) -> Array[ForecastPoint] {
  let predictions = forecast_with_method(data, strategy, horizon)
  let result = []
  let residuals = []
  for value in data {
    residuals.push(value - forecast_mean_or_zero(data))
  }
  let scale = forecast_residual_scale(residuals)
  let confidence_value = forecast_safe_confidence(confidence)
  let z = if confidence_value >= 0.99 {
    2.58
  } else if confidence_value >= 0.9 {
    1.96
  } else {
    1.64
  }
  for index = 0; index < predictions.length(); index = index + 1 {
    let value = predictions[index]
    result.push({
      horizon: index + 1,
      value,
      lower: value - z * scale * (index + 1).to_double().sqrt(),
      upper: value + z * scale * (index + 1).to_double().sqrt(),
    })
  }
  result
}

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

///|
pub fn forecast_residuals(
  data : Array[Double],
  strategy : ForecastMethod,
) -> Array[Double] {
  let fitted = forecast_fitted_values(data, strategy)
  let result = []
  for index = 0; index < data.length(); index = index + 1 {
    result.push(data[index] - fitted[index])
  }
  result
}

///|
pub fn forecast_mae(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  mean_absolute_error(actual, predicted)
}

///|
pub fn forecast_mse(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  mean_squared_error(actual, predicted)
}

///|
pub fn forecast_rmse(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  root_mean_squared_error(actual, predicted)
}

///|
pub fn forecast_median_ae(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  median_absolute_error(actual, predicted)
}

///|
pub fn forecast_bias(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() != predicted.length() || actual.length() == 0 {
    0.0
  } else {
    let errors = []
    for index = 0; index < actual.length(); index = index + 1 {
      errors.push(predicted[index] - actual[index])
    }
    mean(errors)
  }
}

///|
pub fn forecast_mape(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() != predicted.length() || actual.length() == 0 {
    return 0.0
  }
  let errors = []
  for index = 0; index < actual.length(); index = index + 1 {
    let denominator = abs_double(actual[index])
    if denominator > 1.0e-12 {
      errors.push(abs_double(actual[index] - predicted[index]) / denominator)
    }
  }
  if errors.length() == 0 {
    0.0
  } else {
    mean(errors)
  }
}

///|
pub fn forecast_smape(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() != predicted.length() || actual.length() == 0 {
    return 0.0
  }
  let errors = []
  for index = 0; index < actual.length(); index = index + 1 {
    let denominator = abs_double(actual[index]) + abs_double(predicted[index])
    if denominator > 1.0e-12 {
      errors.push(
        2.0 * abs_double(actual[index] - predicted[index]) / denominator,
      )
    }
  }
  if errors.length() == 0 {
    0.0
  } else {
    mean(errors)
  }
}

///|
pub fn forecast_wape(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() != predicted.length() || actual.length() == 0 {
    0.0
  } else {
    let errors = []
    for index = 0; index < actual.length(); index = index + 1 {
      errors.push(abs_double(actual[index] - predicted[index]))
    }
    let denominator = sum_absolute(actual)
    if denominator == 0.0 {
      0.0
    } else {
      sum_absolute(errors) / denominator
    }
  }
}

///|
pub fn forecast_coverage(
  actual : Array[Double],
  lower : Array[Double],
  upper : Array[Double],
) -> Double {
  if actual.length() != lower.length() ||
    actual.length() != upper.length() ||
    actual.length() == 0 {
    return 0.0
  }
  let mut count = 0
  for index = 0; index < actual.length(); index = index + 1 {
    if actual[index] >= lower[index] && actual[index] <= upper[index] {
      count += 1
    }
  }
  count.to_double() / actual.length().to_double()
}

///|
pub fn forecast_interval_width(
  lower : Array[Double],
  upper : Array[Double],
) -> Double {
  if lower.length() != upper.length() || lower.length() == 0 {
    0.0
  } else {
    let widths = []
    for index = 0; index < lower.length(); index = index + 1 {
      widths.push(upper[index] - lower[index])
    }
    mean(widths)
  }
}

///|
pub fn forecast_report(
  actual : Array[Double],
  predicted : Array[Double],
  lower : Array[Double],
  upper : Array[Double],
) -> ForecastReport {
  let residuals = []
  let count = if actual.length() < predicted.length() {
    actual.length()
  } else {
    predicted.length()
  }
  for index = 0; index < count; index = index + 1 {
    residuals.push(actual[index] - predicted[index])
  }
  {
    predictions: predicted,
    residuals,
    mae: forecast_mae(actual, predicted),
    rmse: forecast_rmse(actual, predicted),
    median_absolute_error: forecast_median_ae(actual, predicted),
    bias: forecast_bias(actual, predicted),
    coverage: forecast_coverage(actual, lower, upper),
    interval_width: forecast_interval_width(lower, upper),
  }
}

///|
pub fn forecast_report_score(report : ForecastReport) -> Double {
  1.0 / (1.0 + report.mae + report.rmse + abs_double(report.bias))
}

///|
pub fn forecast_report_lines(report : ForecastReport) -> Array[String] {
  [
    "mae=" + report.mae.to_string(),
    "rmse=" + report.rmse.to_string(),
    "median_absolute_error=" + report.median_absolute_error.to_string(),
    "bias=" + report.bias.to_string(),
    "coverage=" + report.coverage.to_string(),
    "interval_width=" + report.interval_width.to_string(),
  ]
}

///|
pub fn forecast_report_string(report : ForecastReport) -> String {
  forecast_report_lines(report).join("\n")
}

///|
pub fn forecast_method_score(
  data : Array[Double],
  strategy : ForecastMethod,
) -> Double {
  let actual = data
  let predicted = forecast_fitted_values(data, strategy)
  forecast_rmse(actual, predicted)
}

///|
pub fn forecast_method_scores(
  data : Array[Double],
  methods : Array[ForecastMethod],
) -> Array[Double] {
  let result = []
  for strategy in methods {
    result.push(forecast_method_score(data, strategy))
  }
  result
}

///|
pub fn forecast_best_method(
  data : Array[Double],
  methods : Array[ForecastMethod],
) -> ForecastMethod {
  if methods.length() == 0 {
    forecast_method_catalog()[0]
  } else {
    let mut best = methods[0]
    let mut best_score = forecast_method_score(data, best)
    for index = 1; index < methods.length(); index = index + 1 {
      let score = forecast_method_score(data, methods[index])
      if score < best_score {
        best = methods[index]
        best_score = score
      }
    }
    best
  }
}

///|
pub fn forecast_cross_validate(
  data : Array[Double],
  strategy : ForecastMethod,
  folds : Int,
) -> Array[Double] {
  let result = []
  if folds <= 0 || data.length() < 3 {
    return result
  }
  for fold = 1; fold <= folds; fold = fold + 1 {
    let train_size = data.length() * fold / (folds + 1)
    if train_size > 0 && train_size < data.length() {
      let train_data = forecast_prefix(data, train_size)
      let holdout = []
      for index = train_size; index < data.length(); index = index + 1 {
        holdout.push(data[index])
      }
      let prediction = forecast_with_method(
        train_data,
        strategy,
        holdout.length(),
      )
      result.push(forecast_rmse(holdout, prediction))
    }
  }
  result
}

///|
pub fn forecast_learning_curve(
  data : Array[Double],
  strategy : ForecastMethod,
  steps : Int,
) -> Array[Double] {
  let result = []
  if steps <= 0 {
    return result
  }
  for step = 1; step <= steps; step = step + 1 {
    let size = data.length() * step / steps
    if size > 0 {
      result.push(forecast_method_score(forecast_prefix(data, size), strategy))
    } else {
      result.push(0.0)
    }
  }
  result
}

///|
pub fn forecast_residual_scale_report(
  actual : Array[Double],
  predicted : Array[Double],
) -> Array[Double] {
  let residuals = []
  let count = if actual.length() < predicted.length() {
    actual.length()
  } else {
    predicted.length()
  }
  for index = 0; index < count; index = index + 1 {
    residuals.push(actual[index] - predicted[index])
  }
  [
    mean(residuals),
    mad(residuals),
    sample_stddev(residuals),
    quantile(residuals, 0.05),
    quantile(residuals, 0.5),
    quantile(residuals, 0.95),
  ]
}

///|
pub fn forecast_direction_accuracy(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  if actual.length() < 2 || actual.length() != predicted.length() {
    return 0.0
  }
  let mut correct = 0
  let mut total = 0
  for index = 1; index < actual.length(); index = index + 1 {
    let actual_change = actual[index] - actual[index - 1]
    let predicted_change = predicted[index] - predicted[index - 1]
    if sign_double(actual_change) == sign_double(predicted_change) {
      correct += 1
    }
    total += 1
  }
  correct.to_double() / total.to_double()
}

///|
pub fn forecast_stability(
  data : Array[Double],
  strategy : ForecastMethod,
) -> Double {
  let predictions = forecast_with_method(data, strategy, 3)
  if predictions.length() <= 1 {
    1.0
  } else {
    1.0 / (1.0 + total_variation(predictions))
  }
}

///|
pub fn forecast_robustness(
  data : Array[Double],
  strategy : ForecastMethod,
  contamination : Double,
) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let clean = forecast_with_method(data, strategy, 1)[0]
  let contaminated = data.copy()
  contaminated.push(max_value(data) + contamination)
  let changed = forecast_with_method(contaminated, strategy, 1)[0]
  1.0 / (1.0 + abs_double(changed - clean))
}

///|
pub fn forecast_residual_autocorrelation(
  actual : Array[Double],
  predicted : Array[Double],
  lag : Int,
) -> Double {
  let residuals = []
  let count = if actual.length() < predicted.length() {
    actual.length()
  } else {
    predicted.length()
  }
  for index = 0; index < count; index = index + 1 {
    residuals.push(actual[index] - predicted[index])
  }
  robust_autocorrelation(residuals, lag)
}

///|
pub fn forecast_baseline_gain(
  actual : Array[Double],
  candidate : Array[Double],
  baseline : Array[Double],
) -> Double {
  let candidate_error = forecast_rmse(actual, candidate)
  let baseline_error = forecast_rmse(actual, baseline)
  if baseline_error == 0.0 {
    0.0
  } else {
    (baseline_error - candidate_error) / baseline_error
  }
}

///|
pub fn forecast_ensemble(
  data : Array[Double],
  methods : Array[ForecastMethod],
  horizon : Int,
) -> Array[Double] {
  let predictions = []
  let count = forecast_safe_horizon(horizon)
  for step in 0.. Array[Double] {
  let count = forecast_safe_horizon(horizon)
  if methods.length() != weights.length() || methods.length() == 0 {
    return forecast_naive(data, count)
  }
  let predictions = []
  for strategy in methods {
    predictions.push(forecast_with_method(data, strategy, count))
  }
  let result = []
  for step in 0.. 0.0 {
        numerator += predictions[index][step] * weights[index]
        denominator += weights[index]
      }
    }
    result.push(
      if denominator == 0.0 {
        median(values)
      } else {
        numerator / denominator
      },
    )
  }
  result
}

///|
pub fn RobustForecaster::new(
  strategy : ForecastMethod,
  horizon : Int,
  confidence? : Double = 0.95,
) -> RobustForecaster {
  {
    strategy,
    horizon: forecast_safe_horizon(horizon),
    confidence: forecast_safe_confidence(confidence),
    history: [],
    last_forecast: 0.0,
  }
}

///|
pub fn RobustForecaster::push(
  self : RobustForecaster,
  value : Double,
) -> Double {
  self.history.push(value)
  let prediction = forecast_with_method(
    self.history,
    self.strategy,
    self.horizon,
  )
  self.last_forecast = prediction[0]
  self.last_forecast
}

///|
pub fn RobustForecaster::push_many(
  self : RobustForecaster,
  values : Array[Double],
) -> Array[Double] {
  let result = []
  for value in values {
    result.push(self.push(value))
  }
  result
}

///|
pub fn RobustForecaster::predict(
  self : RobustForecaster,
) -> Array[ForecastPoint] {
  forecast_points(
    self.history,
    self.strategy,
    self.horizon,
    confidence=self.confidence,
  )
}

///|
pub fn RobustForecaster::history(self : RobustForecaster) -> Array[Double] {
  self.history.copy()
}

///|
pub fn RobustForecaster::last(self : RobustForecaster) -> Double {
  self.last_forecast
}

///|
pub fn RobustForecaster::reset(self : RobustForecaster) -> Unit {
  self.history = []
  self.last_forecast = 0.0
}

///|
pub fn RobustForecaster::strategy(self : RobustForecaster) -> ForecastMethod {
  self.strategy
}

///|
pub fn RobustForecaster::horizon(self : RobustForecaster) -> Int {
  self.horizon
}

///|
pub fn RobustForecaster::confidence(self : RobustForecaster) -> Double {
  self.confidence
}