//! Median Absolute Deviation (MAD) computation for robust scale estimation.
//!
//! ## Purpose
//!
//! This module provides robust scale estimation using the Median Absolute
//! Deviation (MAD), which is resistant to outliers. MAD is used in LOWESS
//! for computing robustness weights that downweight outliers in iterative
//! refinement.
//!
//! ## Design notes
//!
//! * MAD is computed as: MAD = median(|rᵢ - median(r)|)
//! * Uses in-place selection algorithms for O(n) average-case performance.
//! * Reuses allocated buffers to minimize memory allocations.
//! * All functions are generic over `Float` types to support f32 and f64.
//! * Supports both `std` and `no_std` environments.
//!
//! ## Key concepts
//!
//! ### Median Absolute Deviation
//! MAD measures the variability of a dataset by computing the median of
//! absolute deviations from the median. It is more robust to outliers
//! than standard deviation.
//!
//! ### Robustness
//! MAD has a breakdown point of 50%, meaning it remains valid even when
//! up to 50% of the data are outliers. This makes it ideal for iterative
//! outlier detection in LOWESS.
//!
//! ### Computational Efficiency
//! The implementation uses `select_nth_unstable_by` for O(n) average-case
//! median computation, avoiding the O(n log n) cost of full sorting.
//!
//! ## Invariants
//!
//! * MAD is always non-negative.
//! * MAD is zero for datasets with ≤1 elements.
//! * MAD is zero when all values are identical.
//! * The median computation handles both even and odd-length arrays correctly.
//!
//! ## Non-goals
//!
//! * This module does not validate input data for NaN or infinity.
//! * This module does not provide other robust scale estimators (e.g., IQR).
//! * This module does not handle weighted MAD computation.
//!
//! ## Visibility
//!
//! This module is an internal implementation detail used by the robustness
//! iteration logic. It is not part of the public API and may change without
//! notice.

#[cfg(not(feature = "std"))]
extern crate alloc;
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;

use core::cmp::Ordering;
use num_traits::Float;

// ============================================================================
// MAD Computation
// ============================================================================

/// Compute Median Absolute Deviation (MAD) of residuals.
///
/// MAD is a robust measure of scale, defined as:
/// ```text
/// MAD = median(|rᵢ - median(r)|)
/// ```
///
/// This is resistant to outliers and has a 50% breakdown point, making it
/// ideal for robust regression methods like LOWESS.
///
/// # Parameters
///
/// * `residuals` - Residual values (y - ŷ)
///
/// # Returns
///
/// The median absolute deviation. Returns zero for inputs with ≤1 elements
/// or when all values are identical.
///
/// # Algorithm
///
/// 1. Compute the median of the input residuals
/// 2. Compute absolute deviations from the median
/// 3. Return the median of these absolute deviations
///
/// # Complexity
///
/// * Time: O(n) average case, O(n²) worst case (due to selection algorithm)
/// * Space: O(n) for the working buffer
///
/// # Implementation notes
///
/// * Uses `select_nth_unstable_by` for efficient median computation
/// * Reuses the same buffer for both median computations to reduce allocations
/// * Handles even-length arrays by averaging the two middle values
pub fn compute_mad<T: Float>(residuals: &[T]) -> T {
    let n = residuals.len();

    // Edge case: need at least 2 points for meaningful MAD
    if n <= 1 {
        return T::zero();
    }

    // Helper function to compute median in-place
    let median_inplace = |vals: &mut [T]| -> T {
        let n = vals.len();

        if n == 0 {
            return T::zero();
        }
        if n == 1 {
            return vals[0];
        }

        let mid = n / 2;

        if n % 2 == 0 {
            // Even length: average of two middle values
            vals.select_nth_unstable_by(mid, |a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));
            let upper = vals[mid];

            // Find the largest value in the lower half
            let lower = vals[..mid].iter().copied().fold(T::neg_infinity(), T::max);

            (lower + upper) / T::from(2.0).unwrap()
        } else {
            // Odd length: middle value
            vals.select_nth_unstable_by(mid, |a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));
            vals[mid]
        }
    };

    // Step 1: Compute median of residuals
    let mut vals: Vec<T> = residuals.to_vec();
    let median: T = median_inplace(&mut vals);

    // Step 2: Compute absolute deviations from median
    // Reuse vals buffer to avoid second allocation
    for val in vals.iter_mut() {
        *val = (*val - median).abs();
    }

    // Step 3: Return median of absolute deviations
    median_inplace(&mut vals)
}
