///|
/// A fixed-width histogram for telemetry and residual distributions.
pub struct HistogramBucket {
  lower : Double
  upper : Double
  count : Int
} derive(Debug)

///|
pub fn HistogramBucket::lower(self : HistogramBucket) -> Double {
  self.lower
}

///|
pub fn HistogramBucket::upper(self : HistogramBucket) -> Double {
  self.upper
}

///|
pub fn HistogramBucket::count(self : HistogramBucket) -> Int {
  self.count
}

///|
pub struct Histogram {
  minimum : Double
  maximum : Double
  buckets : Array[HistogramBucket]
  mut underflow : Int
  mut overflow : Int
  mut samples : Int
}

///|
pub fn Histogram::new(
  minimum : Double,
  maximum : Double,
  bucket_count : Int,
) -> Histogram {
  let safe_minimum = if minimum.is_nan() { 0.0 } else { minimum }
  let safe_maximum = if maximum <= safe_minimum || maximum.is_nan() {
    safe_minimum + 1.0
  } else {
    maximum
  }
  let count = if bucket_count < 1 { 1 } else { bucket_count }
  let width = (safe_maximum - safe_minimum) / count.to_double()
  let buckets : Array[HistogramBucket] = []
  for index in 0.. Bool {
  if value.is_nan() || value.is_inf() {
    return false
  }
  self.samples = self.samples + 1
  if value < self.minimum {
    self.underflow = self.underflow + 1
    return true
  }
  if value > self.maximum {
    self.overflow = self.overflow + 1
    return true
  }
  let width = (self.maximum - self.minimum) / self.buckets.length().to_double()
  let raw = ((value - self.minimum) / width).to_int()
  let index = if raw < 0 {
    0
  } else if raw >= self.buckets.length() {
    self.buckets.length() - 1
  } else {
    raw
  }
  let bucket = self.buckets[index]
  self.buckets[index] = {
    lower: bucket.lower,
    upper: bucket.upper,
    count: bucket.count + 1,
  }
  true
}

///|
pub fn Histogram::add_many(self : Histogram, values : Array[Double]) -> Int {
  let mut accepted = 0
  for value in values {
    if self.add(value) {
      accepted = accepted + 1
    }
  }
  accepted
}

///|
pub fn Histogram::minimum(self : Histogram) -> Double {
  self.minimum
}

///|
pub fn Histogram::maximum(self : Histogram) -> Double {
  self.maximum
}

///|
pub fn Histogram::bucket_count(self : Histogram) -> Int {
  self.buckets.length()
}

///|
pub fn Histogram::buckets(self : Histogram) -> Array[HistogramBucket] {
  self.buckets.copy()
}

///|
pub fn Histogram::underflow(self : Histogram) -> Int {
  self.underflow
}

///|
pub fn Histogram::overflow(self : Histogram) -> Int {
  self.overflow
}

///|
pub fn Histogram::samples(self : Histogram) -> Int {
  self.samples
}

///|
pub fn Histogram::in_range(self : Histogram) -> Int {
  let mut total = 0
  for bucket in self.buckets {
    total = total + bucket.count
  }
  total
}

///|
pub fn Histogram::percentile(self : Histogram, probability : Double) -> Double {
  if self.samples == 0 {
    return 0.0
  }
  let p = if probability < 0.0 {
    0.0
  } else if probability > 1.0 {
    1.0
  } else {
    probability
  }
  let target = (p * (self.samples - 1).to_double()).to_int()
  let mut cumulative = 0
  for bucket in self.buckets {
    cumulative = cumulative + bucket.count
    if cumulative > target {
      return (bucket.lower + bucket.upper) * 0.5
    }
  }
  self.maximum
}

///|
pub fn Histogram::density(self : Histogram) -> Array[Double] {
  let total = self.samples.to_double()
  if total <= 0.0 {
    return Array::make(self.buckets.length(), 0.0)
  }
  self.buckets.map(bucket => bucket.count.to_double() / total)
}

///|
pub fn Histogram::reset(self : Histogram) -> Unit {
  for i in 0.. ExponentialStats {
  let safe_alpha = if alpha <= 0.0 {
    0.01
  } else if alpha > 1.0 {
    1.0
  } else {
    alpha
  }
  {
    alpha: safe_alpha,
    count: 0,
    mean: 0.0,
    variance: 0.0,
    minimum: 0.0,
    maximum: 0.0,
  }
}

///|
pub fn ExponentialStats::add(self : ExponentialStats, value : Double) -> Bool {
  if value.is_nan() || value.is_inf() {
    return false
  }
  if self.count == 0 {
    self.mean = value
    self.minimum = value
    self.maximum = value
  } else {
    let difference = value - self.mean
    self.mean = self.mean + self.alpha * difference
    self.variance = (1.0 - self.alpha) *
      (self.variance + self.alpha * difference * difference)
    if value < self.minimum {
      self.minimum = value
    }
    if value > self.maximum {
      self.maximum = value
    }
  }
  self.count = self.count + 1
  true
}

///|
pub fn ExponentialStats::alpha(self : ExponentialStats) -> Double {
  self.alpha
}

///|
pub fn ExponentialStats::count(self : ExponentialStats) -> Int {
  self.count
}

///|
pub fn ExponentialStats::mean(self : ExponentialStats) -> Double {
  self.mean
}

///|
pub fn ExponentialStats::variance(self : ExponentialStats) -> Double {
  self.variance
}

///|
pub fn ExponentialStats::standard_deviation(self : ExponentialStats) -> Double {
  self.variance.sqrt()
}

///|
pub fn ExponentialStats::minimum(self : ExponentialStats) -> Double {
  self.minimum
}

///|
pub fn ExponentialStats::maximum(self : ExponentialStats) -> Double {
  self.maximum
}

///|
pub fn ExponentialStats::reset(self : ExponentialStats) -> Unit {
  self.count = 0
  self.mean = 0.0
  self.variance = 0.0
  self.minimum = 0.0
  self.maximum = 0.0
}

///|
/// Online vector statistics using a numerically stable rank-one update.
pub struct VectorAccumulator {
  dimension : Int
  mut count : Int
  mut mean : Array[Double]
  mut scatter : Matrix
}

///|
pub fn VectorAccumulator::new(dimension : Int) -> VectorAccumulator {
  let safe_dimension = if dimension < 0 { 0 } else { dimension }
  {
    dimension: safe_dimension,
    count: 0,
    mean: Array::make(safe_dimension, 0.0),
    scatter: Matrix::zeros(safe_dimension, safe_dimension),
  }
}

///|
pub fn VectorAccumulator::add(
  self : VectorAccumulator,
  sample : Array[Double],
) -> Bool {
  if sample.length() != self.dimension || !vector_is_finite(sample) {
    return false
  }
  self.count = self.count + 1
  let delta = vector_sub(sample, self.mean)
  self.mean = vector_axpy(1.0 / self.count.to_double(), delta, self.mean)
  let corrected = vector_sub(sample, self.mean)
  self.scatter = self.scatter.add(Matrix::outer(delta, corrected))
  true
}

///|
pub fn VectorAccumulator::dimension(self : VectorAccumulator) -> Int {
  self.dimension
}

///|
pub fn VectorAccumulator::count(self : VectorAccumulator) -> Int {
  self.count
}

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

///|
pub fn VectorAccumulator::covariance(self : VectorAccumulator) -> Matrix {
  if self.count < 2 {
    Matrix::zeros(self.dimension, self.dimension)
  } else {
    self.scatter.scale(1.0 / (self.count - 1).to_double())
  }
}

///|
pub fn VectorAccumulator::standard_deviation(
  self : VectorAccumulator,
) -> Array[Double] {
  let covariance = self.covariance()
  Array::makei(self.dimension, index => covariance.get(index, index).sqrt())
}

///|
pub fn VectorAccumulator::reset(self : VectorAccumulator) -> Unit {
  self.count = 0
  self.mean = Array::make(self.dimension, 0.0)
  self.scatter = Matrix::zeros(self.dimension, self.dimension)
}

///|
/// A weighted mean useful for blending calibrated sensor channels.
pub struct WeightedAccumulator {
  dimension : Int
  mut total_weight : Double
  mut weighted_sum : Array[Double]
}

///|
pub fn WeightedAccumulator::new(dimension : Int) -> WeightedAccumulator {
  let safe_dimension = if dimension < 0 { 0 } else { dimension }
  {
    dimension: safe_dimension,
    total_weight: 0.0,
    weighted_sum: Array::make(safe_dimension, 0.0),
  }
}

///|
pub fn WeightedAccumulator::add(
  self : WeightedAccumulator,
  value : Array[Double],
  weight : Double,
) -> Bool {
  if value.length() != self.dimension ||
    !vector_is_finite(value) ||
    weight <= 0.0 ||
    weight.is_nan() {
    return false
  }
  self.total_weight = self.total_weight + weight
  self.weighted_sum = vector_axpy(weight, value, self.weighted_sum)
  true
}

///|
pub fn WeightedAccumulator::mean(self : WeightedAccumulator) -> Array[Double] {
  if self.total_weight <= 0.0 {
    Array::make(self.dimension, 0.0)
  } else {
    vector_scale(self.weighted_sum, 1.0 / self.total_weight)
  }
}

///|
pub fn WeightedAccumulator::weight(self : WeightedAccumulator) -> Double {
  self.total_weight
}

///|
pub fn WeightedAccumulator::dimension(self : WeightedAccumulator) -> Int {
  self.dimension
}

///|
pub fn WeightedAccumulator::reset(self : WeightedAccumulator) -> Unit {
  self.total_weight = 0.0
  self.weighted_sum = Array::make(self.dimension, 0.0)
}