///|
/// Loss functions used by robust estimators for contaminated sensor data.
pub(all) enum RobustLossKind {
  Squared
  Huber
  Cauchy
  Tukey
  Welsch
} derive(Debug, Eq)

///|
/// Configuration shared by robust location and regression estimators.
pub struct RobustEstimatorConfig {
  loss : RobustLossKind
  tuning : Double
  iterations : Int
  tolerance : Double
  minimum_weight : Double
} derive(Debug)

///|
/// Construct a robust estimator configuration with safe defaults.
pub fn RobustEstimatorConfig::new(
  loss : RobustLossKind,
  tuning : Double,
  iterations : Int,
  tolerance : Double,
  minimum_weight : Double,
) -> RobustEstimatorConfig {
  {
    loss,
    tuning: if tuning <= 0.0 || tuning.is_nan() {
      1.345
    } else {
      tuning
    },
    iterations: if iterations < 1 {
      1
    } else {
      iterations
    },
    tolerance: if tolerance <= 0.0 || tolerance.is_nan() {
      0.000001
    } else {
      tolerance
    },
    minimum_weight: if minimum_weight < 0.0 {
      0.0
    } else if minimum_weight > 1.0 {
      1.0
    } else {
      minimum_weight
    },
  }
}

///|
/// A configuration tuned for Huber-style sensor rejection.
pub fn robust_huber_config() -> RobustEstimatorConfig {
  RobustEstimatorConfig::new(Huber, 1.345, 20, 0.000001, 0.0001)
}

///|
/// A configuration tuned for strongly contaminated streams.
pub fn robust_tukey_config() -> RobustEstimatorConfig {
  RobustEstimatorConfig::new(Tukey, 4.685, 25, 0.000001, 0.000001)
}

///|
/// Return the loss kind.
pub fn RobustEstimatorConfig::loss(
  self : RobustEstimatorConfig,
) -> RobustLossKind {
  self.loss
}

///|
/// Return the tuning constant.
pub fn RobustEstimatorConfig::tuning(self : RobustEstimatorConfig) -> Double {
  self.tuning
}

///|
/// Return the iteration budget.
pub fn RobustEstimatorConfig::iterations(self : RobustEstimatorConfig) -> Int {
  self.iterations
}

///|
/// Return the convergence tolerance.
pub fn RobustEstimatorConfig::tolerance(self : RobustEstimatorConfig) -> Double {
  self.tolerance
}

///|
/// Return the minimum retained weight.
pub fn RobustEstimatorConfig::minimum_weight(
  self : RobustEstimatorConfig,
) -> Double {
  self.minimum_weight
}

///|
/// A weighted scalar sample.
pub struct RobustSample {
  value : Double
  weight : Double
} derive(Debug)

///|
/// Construct a weighted sample.
pub fn RobustSample::new(value : Double, weight : Double) -> RobustSample {
  { value, weight: if weight < 0.0 || weight.is_nan() { 0.0 } else { weight } }
}

///|
/// Return the sample value.
pub fn RobustSample::value(self : RobustSample) -> Double {
  self.value
}

///|
/// Return the sample weight.
pub fn RobustSample::weight(self : RobustSample) -> Double {
  self.weight
}

///|
/// A scalar robust estimate with diagnostics.
pub struct RobustEstimate {
  location : Double
  scale : Double
  iterations : Int
  effective_samples : Double
  residuals : Array[Double]
  weights : Array[Double]
  converged : Bool
} derive(Debug)

///|
/// Return the estimated location.
pub fn RobustEstimate::location(self : RobustEstimate) -> Double {
  self.location
}

///|
/// Return the robust scale.
pub fn RobustEstimate::scale(self : RobustEstimate) -> Double {
  self.scale
}

///|
/// Return the number of iterations performed.
pub fn RobustEstimate::iterations(self : RobustEstimate) -> Int {
  self.iterations
}

///|
/// Return the effective sample size.
pub fn RobustEstimate::effective_samples(self : RobustEstimate) -> Double {
  self.effective_samples
}

///|
/// Return a copy of residuals.
pub fn RobustEstimate::residuals(self : RobustEstimate) -> Array[Double] {
  self.residuals.copy()
}

///|
/// Return a copy of final weights.
pub fn RobustEstimate::weights(self : RobustEstimate) -> Array[Double] {
  self.weights.copy()
}

///|
/// Return whether the location iteration converged.
pub fn RobustEstimate::converged(self : RobustEstimate) -> Bool {
  self.converged
}

///|
/// A robust straight-line fit `y = intercept + slope * x`.
pub struct RobustLineFit {
  intercept : Double
  slope : Double
  scale : Double
  residuals : Array[Double]
  weights : Array[Double]
  iterations : Int
  converged : Bool
} derive(Debug)

///|
/// Return the fitted intercept.
pub fn RobustLineFit::intercept(self : RobustLineFit) -> Double {
  self.intercept
}

///|
/// Return the fitted slope.
pub fn RobustLineFit::slope(self : RobustLineFit) -> Double {
  self.slope
}

///|
/// Return the residual scale.
pub fn RobustLineFit::scale(self : RobustLineFit) -> Double {
  self.scale
}

///|
/// Return a copy of line-fit residuals.
pub fn RobustLineFit::residuals(self : RobustLineFit) -> Array[Double] {
  self.residuals.copy()
}

///|
/// Return final line-fit weights.
pub fn RobustLineFit::weights(self : RobustLineFit) -> Array[Double] {
  self.weights.copy()
}

///|
/// Return the number of line-fit iterations.
pub fn RobustLineFit::iterations(self : RobustLineFit) -> Int {
  self.iterations
}

///|
/// Return convergence state.
pub fn RobustLineFit::converged(self : RobustLineFit) -> Bool {
  self.converged
}

///|
/// Return the loss value for a standardized residual.
pub fn robust_loss_value(
  kind : RobustLossKind,
  residual : Double,
  tuning : Double,
) -> Double {
  let c = if tuning <= 0.0 { 1.0 } else { tuning }
  let u = residual.abs() / c
  match kind {
    Squared => 0.5 * residual * residual
    Huber =>
      if u <= 1.0 {
        0.5 * residual * residual
      } else {
        c * residual.abs() - 0.5 * c * c
      }
    Cauchy => 0.5 * c * c * (u * u / (1.0 + u * u))
    Tukey =>
      if u < 1.0 {
        c * c / 6.0 * (1.0 - (1.0 - u * u) * (1.0 - u * u) * (1.0 - u * u))
      } else {
        c * c / 6.0
      }
    Welsch => 0.5 * c * c * (1.0 - 1.0 / (1.0 + u * u))
  }
}

///|
/// Return the iteratively reweighted least-squares weight.
pub fn robust_loss_weight(
  kind : RobustLossKind,
  residual : Double,
  tuning : Double,
  minimum_weight : Double,
) -> Double {
  if residual.is_nan() || residual.is_inf() {
    return minimum_weight
  }
  let c = if tuning <= 0.0 { 1.0 } else { tuning }
  let u = residual.abs() / c
  let raw = match kind {
    Squared => 1.0
    Huber => if u <= 1.0 || u == 0.0 { 1.0 } else { 1.0 / u }
    Cauchy => 1.0 / (1.0 + u * u)
    Tukey => if u < 1.0 { (1.0 - u * u) * (1.0 - u * u) } else { 0.0 }
    Welsch => 1.0 / (1.0 + u * u)
  }
  if raw < minimum_weight {
    minimum_weight
  } else if raw > 1.0 {
    1.0
  } else {
    raw
  }
}

///|
/// Compute the median of finite values.
pub fn robust_median(values : Array[Double]) -> Double {
  let finite = []
  for value in values {
    if !value.is_nan() && !value.is_inf() {
      finite.push(value)
    }
  }
  if finite.length() == 0 {
    return 0.0
  }
  for i in 1.. 0 && finite[j - 1] > current {
      finite[j] = finite[j - 1]
      j = j - 1
    }
    finite[j] = current
  }
  let middle = finite.length() / 2
  if finite.length() % 2 == 1 {
    finite[middle]
  } else {
    (finite[middle - 1] + finite[middle]) * 0.5
  }
}

///|
/// Compute the median absolute deviation around a location.
pub fn robust_mad(values : Array[Double], location : Double) -> Double {
  let deviations = Array::makei(values.length(), i => {
    (values[i] - location).abs()
  })
  robust_median(deviations)
}

///|
/// Convert MAD to a Gaussian-consistent standard deviation estimate.
pub fn robust_scale_from_mad(mad : Double) -> Double {
  if mad <= 0.0 || mad.is_nan() {
    0.000001
  } else {
    mad * 1.4826
  }
}

///|
/// Compute a weighted mean while ignoring invalid samples.
pub fn robust_weighted_mean(samples : Array[RobustSample]) -> Double {
  let mut numerator = 0.0
  let mut denominator = 0.0
  for sample in samples {
    if sample.weight() > 0.0 &&
      !sample.value().is_nan() &&
      !sample.value().is_inf() {
      numerator = numerator + sample.value() * sample.weight()
      denominator = denominator + sample.weight()
    }
  }
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Compute a weighted variance around a known location.
pub fn robust_weighted_variance(
  samples : Array[RobustSample],
  location : Double,
) -> Double {
  let mut numerator = 0.0
  let mut denominator = 0.0
  for sample in samples {
    if sample.weight() > 0.0 &&
      !sample.value().is_nan() &&
      !sample.value().is_inf() {
      let delta = sample.value() - location
      numerator = numerator + sample.weight() * delta * delta
      denominator = denominator + sample.weight()
    }
  }
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Estimate a robust location with iteratively reweighted means.
pub fn robust_location(
  values : Array[Double],
  config : RobustEstimatorConfig,
) -> RobustEstimate {
  let initial = robust_median(values)
  let initial_scale = robust_scale_from_mad(robust_mad(values, initial))
  let mut location = initial
  let mut scale = initial_scale
  let weights = Array::make(values.length(), 0.0)
  let residuals = Array::make(values.length(), 0.0)
  let mut converged = false
  let mut performed = 0
  for iteration in 0.. 0.0 && !value.is_nan() && !value.is_inf() {
        numerator = numerator + weight * value
        denominator = denominator + weight
      }
    }
    let next = if denominator == 0.0 {
      location
    } else {
      numerator / denominator
    }
    let change = (next - location).abs()
    location = next
    scale = robust_scale_from_mad(robust_mad(values, location))
    if change <= config.tolerance() {
      converged = true
      break
    }
  }
  let mut effective = 0.0
  for weight in weights {
    effective = effective + weight
  }
  {
    location,
    scale,
    iterations: performed,
    effective_samples: effective,
    residuals,
    weights,
    converged,
  }
}

///|
/// Estimate robust location using explicit base weights.
pub fn robust_weighted_location(
  values : Array[Double],
  base_weights : Array[Double],
  config : RobustEstimatorConfig,
) -> RobustEstimate {
  let samples = []
  let count = if values.length() < base_weights.length() {
    values.length()
  } else {
    base_weights.length()
  }
  for i in 0.. 0.0 {
        base_weights[i]
      } else {
        0.0
      }
      let residual = value - location
      residuals[i] = residual
      let weight = base *
        robust_loss_weight(
          config.loss(),
          residual / scale,
          config.tuning(),
          config.minimum_weight(),
        )
      weights[i] = weight
      numerator = numerator + weight * value
      denominator = denominator + weight
    }
    let next = if denominator == 0.0 {
      location
    } else {
      numerator / denominator
    }
    if (next - location).abs() <= config.tolerance() {
      location = next
      converged = true
      break
    }
    location = next
  }
  let mut effective = 0.0
  for weight in weights {
    effective = effective + weight
  }
  {
    location,
    scale,
    iterations: performed,
    effective_samples: effective,
    residuals,
    weights,
    converged,
  }
}

///|
/// Fit a robust straight line with iteratively reweighted normal equations.
pub fn robust_line_fit(
  x : Array[Double],
  y : Array[Double],
  config : RobustEstimatorConfig,
) -> RobustLineFit {
  let count = if x.length() < y.length() { x.length() } else { y.length() }
  let weights = Array::make(count, 1.0)
  let residuals = Array::make(count, 0.0)
  let mut intercept = 0.0
  let mut slope = 0.0
  let mut scale = 1.0
  let mut converged = false
  let mut performed = 0
  for iteration in 0.. Double {
  self.intercept + self.slope * x
}

///|
/// Compute a robust line's weighted residual sum of squares.
pub fn RobustLineFit::weighted_error(self : RobustLineFit) -> Double {
  let mut result = 0.0
  for i in 0.. Int {
  self.count
}

///|
/// Return finite residual count.
pub fn RobustResidualSummary::finite(self : RobustResidualSummary) -> Int {
  self.finite
}

///|
/// Return rejected/non-finite count.
pub fn RobustResidualSummary::rejected(self : RobustResidualSummary) -> Int {
  self.rejected
}

///|
/// Return residual mean.
pub fn RobustResidualSummary::mean(self : RobustResidualSummary) -> Double {
  self.mean
}

///|
/// Return residual median.
pub fn RobustResidualSummary::median(self : RobustResidualSummary) -> Double {
  self.median
}

///|
/// Return residual MAD.
pub fn RobustResidualSummary::mad(self : RobustResidualSummary) -> Double {
  self.mad
}

///|
/// Return residual RMS.
pub fn RobustResidualSummary::rms(self : RobustResidualSummary) -> Double {
  self.rms
}

///|
/// Return maximum absolute residual.
pub fn RobustResidualSummary::maximum(self : RobustResidualSummary) -> Double {
  self.maximum
}

///|
/// Return positive residual count.
pub fn RobustResidualSummary::positive(self : RobustResidualSummary) -> Int {
  self.positive
}

///|
/// Return negative residual count.
pub fn RobustResidualSummary::negative(self : RobustResidualSummary) -> Int {
  self.negative
}

///|
/// Summarize a residual array while retaining explicit invalid counts.
pub fn summarize_robust_residuals(
  values : Array[Double],
) -> RobustResidualSummary {
  let finite_values = []
  let mut sum = 0.0
  let mut square_sum = 0.0
  let mut maximum = 0.0
  let mut positive = 0
  let mut negative = 0
  let mut rejected = 0
  for value in values {
    if value.is_nan() || value.is_inf() {
      rejected = rejected + 1
    } else {
      finite_values.push(value)
      sum = sum + value
      square_sum = square_sum + value * value
      if value.abs() > maximum {
        maximum = value.abs()
      }
      if value > 0.0 {
        positive = positive + 1
      } else if value < 0.0 {
        negative = negative + 1
      }
    }
  }
  let count = finite_values.length()
  let mean = if count == 0 { 0.0 } else { sum / count.to_double() }
  let rms = if count == 0 {
    0.0
  } else {
    (square_sum / count.to_double()).sqrt()
  }
  let median = robust_median(finite_values)
  {
    count: values.length(),
    finite: count,
    rejected,
    mean,
    median,
    mad: robust_mad(finite_values, median),
    rms,
    maximum,
    positive,
    negative,
  }
}

///|
/// Return robust weights for a residual vector.
pub fn robust_weights_for_residuals(
  residuals : Array[Double],
  config : RobustEstimatorConfig,
  scale : Double,
) -> Array[Double] {
  let safe_scale = if scale <= 0.0 { 1.0 } else { scale }
  Array::makei(residuals.length(), i => {
    robust_loss_weight(
      config.loss(),
      residuals[i] / safe_scale,
      config.tuning(),
      config.minimum_weight(),
    )
  })
}

///|
/// Down-weight an observation vector using one scalar residual per row.
pub fn robust_reweight_rows(
  rows : Array[Array[Double]],
  residuals : Array[Double],
  config : RobustEstimatorConfig,
  scale : Double,
) -> Array[Array[Double]] {
  let weights = robust_weights_for_residuals(residuals, config, scale)
  Array::makei(rows.length(), i => {
    let weight = if i < weights.length() { weights[i] } else { 0.0 }
    Array::makei(rows[i].length(), j => rows[i][j] * weight)
  })
}

///|
/// Compute a weighted covariance of scalar residuals.
pub fn robust_residual_variance(
  residuals : Array[Double],
  weights : Array[Double],
) -> Double {
  let count = if residuals.length() < weights.length() {
    residuals.length()
  } else {
    weights.length()
  }
  let mut numerator = 0.0
  let mut denominator = 0.0
  for i in 0.. 0.0 {
      numerator = numerator + weights[i] * residuals[i] * residuals[i]
      denominator = denominator + weights[i]
    }
  }
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Estimate a robust scale from a weighted residual vector.
pub fn robust_weighted_scale(
  residuals : Array[Double],
  weights : Array[Double],
) -> Double {
  let variance = robust_residual_variance(residuals, weights)
  if variance <= 0.0 {
    0.000001
  } else {
    variance.sqrt()
  }
}

///|
/// Return the fraction of samples whose robust weight is above a threshold.
pub fn robust_inlier_fraction(
  weights : Array[Double],
  threshold : Double,
) -> Double {
  if weights.length() == 0 {
    return 0.0
  }
  let mut count = 0
  for weight in weights {
    if weight >= threshold {
      count = count + 1
    }
  }
  count.to_double() / weights.length().to_double()
}

///|
/// Compute a robust weighted quantile using sorted finite values.
pub fn robust_weighted_quantile(
  values : Array[Double],
  weights : Array[Double],
  probability : Double,
) -> Double {
  let entries = []
  let count = if values.length() < weights.length() {
    values.length()
  } else {
    weights.length()
  }
  for i in 0.. 0.0 {
      entries.push((values[i], weights[i]))
    }
  }
  if entries.length() == 0 {
    return 0.0
  }
  for i in 1.. 0 && entries[j - 1].0 > current.0 {
      entries[j] = entries[j - 1]
      j = j - 1
    }
    entries[j] = current
  }
  let mut total = 0.0
  for entry in entries {
    total = total + entry.1
  }
  let target = total * probability.clamp(min=0.0, max=1.0)
  let mut cumulative = 0.0
  for entry in entries {
    cumulative = cumulative + entry.1
    if cumulative >= target {
      return entry.0
    }
  }
  entries[entries.length() - 1].0
}

///|
/// Compute a robust confidence radius from residuals and a confidence scale.
pub fn robust_confidence_radius(
  residuals : Array[Double],
  confidence_scale : Double,
) -> Double {
  let summary = summarize_robust_residuals(residuals)
  let multiplier = if confidence_scale < 0.0 { 0.0 } else { confidence_scale }
  summary.mad() * 1.4826 * multiplier
}

///|
/// Apply a Huber-style clipping policy to a scalar residual.
pub fn huber_clip_residual(residual : Double, threshold : Double) -> Double {
  let limit = threshold.abs()
  if limit == 0.0 {
    0.0
  } else if residual > limit {
    limit
  } else if residual < -limit {
    -limit
  } else {
    residual
  }
}

///|
/// Apply symmetric clipping to a vector.
pub fn clip_residual_vector(
  values : Array[Double],
  threshold : Double,
) -> Array[Double] {
  Array::makei(values.length(), i => huber_clip_residual(values[i], threshold))
}

///|
/// Compute a robust score in `[0, 1]` from residuals and a scale limit.
pub fn robust_quality_score(
  residuals : Array[Double],
  scale_limit : Double,
) -> Double {
  let summary = summarize_robust_residuals(residuals)
  if summary.finite() == 0 {
    return 0.0
  }
  let limit = if scale_limit <= 0.0 { 1.0 } else { scale_limit }
  let score = 1.0 - summary.rms() / limit
  score.clamp(min=0.0, max=1.0)
}