///|
/// Result of a Rauch-Tung-Striebel backward pass.
pub struct SmoothingResult {
  states : Array[Array[Double]]
  covariances : Array[Matrix]
} derive(Debug)

///|
pub fn SmoothingResult::states(self : SmoothingResult) -> Array[Array[Double]] {
  self.states.map(state => state.copy())
}

///|
pub fn SmoothingResult::covariances(self : SmoothingResult) -> Array[Matrix] {
  self.covariances.map(covariance => covariance.copy())
}

///|
pub fn SmoothingResult::length(self : SmoothingResult) -> Int {
  self.states.length()
}

///|
/// Rauch-Tung-Striebel smoother for a linear model.
pub fn rts_smooth(
  filtered_states : Array[Array[Double]],
  filtered_covariances : Array[Matrix],
  predicted_states : Array[Array[Double]],
  predicted_covariances : Array[Matrix],
  transitions : Array[Matrix],
) -> SmoothingResult {
  let count = filtered_states.length()
  if count == 0 || filtered_covariances.length() != count {
    return { states: [], covariances: [] }
  }
  let states = filtered_states.map(state => state.copy())
  let covariances = filtered_covariances.map(covariance => covariance.copy())
  if predicted_states.length() < count - 1 ||
    predicted_covariances.length() < count - 1 ||
    transitions.length() < count - 1 {
    return { states, covariances }
  }
  for index = count - 2; index >= 0; index = index - 1 {
    let inverse = predicted_covariances[index].inverse()
    match inverse {
      None => continue
      Some(predicted_inverse) => {
        let gain = filtered_covariances[index]
          .multiply(transitions[index].transpose())
          .multiply(predicted_inverse)
        let residual = vector_sub(states[index + 1], predicted_states[index])
        let corrected = vector_add(
          filtered_states[index],
          gain.multiply_vector(residual),
        )
        let covariance = filtered_covariances[index]
          .add(
            gain
            .multiply(covariances[index + 1].sub(predicted_covariances[index]))
            .multiply(gain.transpose()),
          )
          .symmetric_part()
        states[index] = corrected
        covariances[index] = covariance
      }
    }
  }
  { states, covariances }
}

///|
/// Fixed-lag history manager. It stores the data needed to smooth a bounded
/// recent window while returning the newest filtered estimate immediately.
pub struct FixedLagSmoother {
  lag : Int
  states : Array[Array[Double]]
  covariances : Array[Matrix]
  predicted_states : Array[Array[Double]]
  predicted_covariances : Array[Matrix]
  transitions : Array[Matrix]
}

///|
pub fn FixedLagSmoother::new(lag : Int) -> FixedLagSmoother {
  {
    lag: if lag < 1 {
      1
    } else {
      lag
    },
    states: [],
    covariances: [],
    predicted_states: [],
    predicted_covariances: [],
    transitions: [],
  }
}

///|
pub fn FixedLagSmoother::lag(self : FixedLagSmoother) -> Int {
  self.lag
}

///|
pub fn FixedLagSmoother::length(self : FixedLagSmoother) -> Int {
  self.states.length()
}

///|
pub fn FixedLagSmoother::push(
  self : FixedLagSmoother,
  filtered_state : Array[Double],
  filtered_covariance : Matrix,
  predicted_state : Array[Double],
  predicted_covariance : Matrix,
  transition : Matrix,
) -> Unit {
  self.states.push(filtered_state.copy())
  self.covariances.push(filtered_covariance.copy())
  self.predicted_states.push(predicted_state.copy())
  self.predicted_covariances.push(predicted_covariance.copy())
  self.transitions.push(transition.copy())
}

///|
pub fn FixedLagSmoother::smooth(self : FixedLagSmoother) -> SmoothingResult {
  let total = self.states.length()
  let start = if total > self.lag { total - self.lag } else { 0 }
  let states : Array[Array[Double]] = []
  let covariances : Array[Matrix] = []
  let predicted_states : Array[Array[Double]] = []
  let predicted_covariances : Array[Matrix] = []
  let transitions : Array[Matrix] = []
  for i in start.. Array[Double] {
  let result = self.smooth()
  match result.states().get(result.length() - 1) {
    None => []
    Some(state) => state
  }
}

///|
pub fn moving_average(values : Array[Double], radius : Int) -> Array[Double] {
  if values.length() == 0 {
    return []
  }
  let safe_radius = if radius < 0 { 0 } else { radius }
  Array::makei(values.length(), index => {
    let start = if index > safe_radius { index - safe_radius } else { 0 }
    let end = if index + safe_radius + 1 > values.length() {
      values.length()
    } else {
      index + safe_radius + 1
    }
    let mut total = 0.0
    for i in start.. Array[Double] {
  if values.length() == 0 {
    return []
  }
  let safe_radius = if radius < 0 { 0 } else { radius }
  Array::makei(values.length(), index => {
    let start = if index > safe_radius { index - safe_radius } else { 0 }
    let end = if index + safe_radius + 1 > values.length() {
      values.length()
    } else {
      index + safe_radius + 1
    }
    let window : Array[Double] = []
    for i in start.. Array[Double] {
  if values.length() == 0 {
    return []
  }
  let alpha = if factor < 0.0 {
    0.0
  } else if factor > 1.0 {
    1.0
  } else {
    factor
  }
  let result = Array::make(values.length(), 0.0)
  result[0] = values[0]
  for i in 1.. Array[Double] {
  let result : Array[Double] = []
  let mut previous = fallback
  let weight = if alpha < 0.0 { 0.0 } else if alpha > 1.0 { 1.0 } else { alpha }
  for value in values {
    match value {
      None => result.push(previous)
      Some(current) => {
        previous = if result.length() == 0 {
          current
        } else {
          weight * current + (1.0 - weight) * previous
        }
        result.push(previous)
      }
    }
  }
  result
}