///|
/// 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")
}