///|
/// Adaptive reweighted penalized least squares with sample deviation of
/// negative residuals and a bounded logistic update. At least two distinct
/// negative residuals are required for another statistically defined update.
pub fn arpls(
  input : Spectrum,
  options~ : IterationOptions = iteration_options(),
) -> BaselineResult raise SpectralError {
  validate_iteration_options(options)
  let y = input.values
  let mut weights = Array::make(y.length(), 1.0)
  for iteration in 1..<=options.max_iterations {
    let baseline = whittaker(y, weights, options.lambda)
    let negatives : Array[Double] = []
    let residual = Array::makei(y.length(), fn(i) { y[i] - baseline[i] })
    for value in residual {
      if value < 0.0 {
        negatives.push(value)
      }
    }
    let residual_scale = root_mean_square(residual)
    let signal_scale = root_mean_square(y)
    if residual_scale <= 1.0e-10 * (1.0 + signal_scale) {
      return baseline_result(y, baseline, weights, iteration, true, 0.0)
    }
    if negatives.length() < 2 {
      return baseline_result(y, baseline, weights, iteration, false, 1.0)
    }
    let deviation = sample_deviation(negatives)
    if deviation <= 1.0e-14 * (1.0 + signal_scale) {
      return baseline_result(y, baseline, weights, iteration, false, 1.0)
    }
    let threshold = 2.0 * deviation - mean(negatives)
    let next = Array::makei(y.length(), fn(i) {
      let exponent = 2.0 * (residual[i] - threshold) / deviation
      if exponent >= 40.0 {
        1.0e-8
      } else if exponent <= -40.0 {
        1.0
      } else {
        let weight = 1.0 / (1.0 + @math.exp(exponent))
        if weight < 1.0e-8 {
          1.0e-8
        } else {
          weight
        }
      }
    })
    let change = relative_difference(next, weights)
    let converged = change <= options.tolerance
    if converged || iteration == options.max_iterations {
      return baseline_result(y, baseline, weights, iteration, converged, change)
    }
    weights = next
  }
  raise Numerical("unreachable arPLS iteration state")
}