lowess/

utils.rs

//! Utility functions for LOWESS smoothing.
//!
//! This module provides lightweight, production-safe helpers used throughout
//! the LOWESS pipeline. The functions are intentionally defensive: they
//! validate inputs, document preconditions, and return safe values on
//! degenerate inputs rather than panicking in release builds.
//!
//! Global expectations
//! - Callers should prefer to pre-clean data (remove NaNs/infs) and sort `x`
//!   when windowing behavior depends on monotonic order. Several helpers
//!   assert correctness with debug-only checks but return safe outputs in
//!   production paths.
//!
//! Key parameters and their semantics
//! - x: slice of independent variable values. Many helpers assume `x` is
//!   sorted ascending; if not, use `sort_by_x` before window-based ops.
//! - y: slice of dependent variable values; must be aligned with `x`.
//! - frac: smoothing fraction in (0, 1]. Used to compute window sizes.
//!   validate_inputs will reject non-finite or out-of-range values.
//! - window_size (usize): discrete neighborhood size computed from `frac`
//!   and `n` via `calculate_window_size`. Always at least 2 and at most `n`.
//! - delta: interpolation distance threshold (`Option<T>` or T). When None
//!   `calculate_delta` returns a conservative default (~1% of x-range). A
//!   zero or negative delta disables the skip/fast-path logic.
//! - left/right/current/idx: integer indices describing the local window
//!   boundaries; many window helpers return clamped values in [0, n-1].
//!
//! Important functions (brief)
//! - validate_inputs(x, y, frac): checks lengths, finiteness, minimum points,
//!   and fraction bounds. Returns a Result for early failure in callers.
//! - validate_confidence_level(level): ensure level ∈ (0,1) for interval code.
//! - validate_delta(delta): ensures delta ≥ 0 and finite.
//! - sort_by_x: stable mapping of x/y to sorted order (returns new Vecs).
//! - calculate_delta(delta, x_sorted): resolves optional delta to a numeric
//!   value (defaults to 1% of range for None).
//! - calculate_window_size(n, frac): converts fractional span to integer
//!   window size with safe min/max clamping.
//! - initialize_window / update_window: establish and slide local windows
//!   while keeping them valid for regression computations.
//! - interpolate_gap: linear interpolation between two fitted points; used
//!   when delta indicates points can be filled in rather than re-fit.
//! - skip_close_points: delta-driven fast-path to skip/refill sequences of
//!   points close to the last fitted x. Also handles identical-x ties by
//!   copying the last fitted value for stability.
//! - normalize_weights / normalize_all_weights: numeric-safe normalization
//!   with uniform fallback when total weight is (near) zero.
//! - compute_range / find_rightmost_point: range and threshold helpers used
//!   by windowing and delta logic.
//! - compute_weighted_average / is_effectively_zero: small numeric utilities
//!   used widely across fitting and diagnostics code.
//!
//! Production recommendations
//! - Validate and deduplicate inputs upstream for large datasets. The helpers
//!   here validate and handle edge cases, but upstream cleaning improves
//!   determinism and performance.
//! - Use `delta` (default ~1% of range) for dense inputs to accelerate fitting.
//! - Choose `frac` in (0, 1]; cross-validation at the builder layer is
//!   recommended for data-driven selection.

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

use crate::{LowessError, Result};
use core::cmp::Ordering;
use num_traits::Float;

// ============================================================================
// Input Validation
// ============================================================================

/// Validate input arrays for LOWESS smoothing.
pub fn validate_inputs<T: Float>(x: &[T], y: &[T], frac: T) -> Result<()> {
    // Empty input is an error
    if x.is_empty() || y.is_empty() {
        return Err(LowessError::EmptyInput);
    }

    // Length mismatch
    if x.len() != y.len() {
        return Err(LowessError::MismatchedInputs {
            x_len: x.len(),
            y_len: y.len(),
        });
    }

    // Require at least two points for a meaningful local regression
    if x.len() < 2 {
        return Err(LowessError::TooFewPoints {
            got: x.len(),
            min: 2,
        });
    }

    // Reject NaN / infinite in inputs
    for (i, &xi) in x.iter().enumerate() {
        if !xi.is_finite() {
            return Err(LowessError::InvalidNumericValue(format!(
                "x[{}]={}",
                i,
                xi.to_f64().unwrap_or(f64::NAN)
            )));
        }
    }
    for (i, &yi) in y.iter().enumerate() {
        if !yi.is_finite() {
            return Err(LowessError::InvalidNumericValue(format!(
                "y[{}]={}",
                i,
                yi.to_f64().unwrap_or(f64::NAN)
            )));
        }
    }

    // Validate fraction
    if !frac.is_finite() || frac <= T::zero() || frac > T::one() {
        return Err(LowessError::InvalidFraction(
            frac.to_f64().unwrap_or(f64::NAN),
        ));
    }

    Ok(())
}

/// Validate confidence level parameter.
pub fn validate_confidence_level<T: Float>(level: T) -> Result<()> {
    if !level.is_finite() || level <= T::zero() || level >= T::one() {
        return Err(LowessError::InvalidConfidenceLevel(
            level.to_f64().unwrap_or(f64::NAN),
        ));
    }
    Ok(())
}

/// Validate delta parameter.
pub fn validate_delta<T: Float>(delta: T) -> Result<()> {
    if !delta.is_finite() || delta < T::zero() {
        return Err(LowessError::InvalidDelta(
            delta.to_f64().unwrap_or(f64::NAN),
        ));
    }
    Ok(())
}

// ============================================================================
// Data Preparation
// ============================================================================

/// Sort data by x values and return sorted (x, y) pairs.
pub fn sort_by_x<T: Float>(x: &[T], y: &[T]) -> (Vec<T>, Vec<T>) {
    debug_assert_eq!(x.len(), y.len(), "Length mismatch in sort_by_x");
    let mut indices: Vec<usize> = (0..x.len()).collect();
    indices.sort_unstable_by(|&i, &j| x[i].partial_cmp(&x[j]).unwrap_or(Ordering::Equal));
    let x_sorted = indices.iter().map(|&i| x[i]).collect();
    let y_sorted = indices.iter().map(|&i| y[i]).collect();
    (x_sorted, y_sorted)
}

/// Check if x values are sorted in ascending order.
pub fn is_sorted<T: Float>(x: &[T]) -> bool {
    x.windows(2).all(|w| w[0] <= w[1])
}

// ============================================================================
// Delta Calculation
// ============================================================================

/// Calculate delta parameter for interpolation optimization.
pub fn calculate_delta<T: Float>(delta: Option<T>, x_sorted: &[T]) -> Result<T> {
    match delta {
        Some(d) => {
            validate_delta(d)?;
            Ok(d)
        }
        None => {
            if x_sorted.is_empty() {
                Ok(T::zero())
            } else {
                let range = x_sorted[x_sorted.len() - 1] - x_sorted[0];
                Ok(T::from(0.01).unwrap() * range)
            }
        }
    }
}

// ============================================================================
// Window Size Calculation
// ============================================================================

/// Calculate window size from fraction and number of points.
pub fn calculate_window_size<T: Float>(n: usize, frac: T) -> usize {
    let frac_n = frac * T::from(n).unwrap();
    let frac_n_int = frac_n.to_usize().unwrap_or(0);
    usize::max(2, usize::min(n, frac_n_int))
}

// ============================================================================
// Window Management
// ============================================================================

/// Update window boundaries to maintain optimal window around current point.
pub fn update_window<T: Float>(
    x: &[T],
    current: usize,
    left: &mut usize,
    right: &mut usize,
    n: usize,
) {
    while *right < n - 1 {
        let d_left = x[current] - x[*left];
        let d_right = x[*right + 1] - x[current];
        if d_left <= d_right {
            break;
        }
        *left += 1;
        *right += 1;
    }
}

/// Initialize window boundaries for a given point.
pub fn initialize_window(idx: usize, window_size: usize, n: usize) -> (usize, usize) {
    // If requested window covers (or exceeds) the whole array, return full range.
    if window_size >= n {
        return (0, n.saturating_sub(1));
    }

    let half = window_size / 2;
    // Start by centering the window around idx.
    let mut left = idx.saturating_sub(half);

    // Clamp left so the window fits in [0, n-1] and has the requested size.
    // Since window_size < n here, n - window_size is well-defined.
    let max_left = n - window_size;
    if left > max_left {
        left = max_left;
    }

    let right = left + window_size.saturating_sub(1);

    debug_assert_eq!(
        right - left + 1,
        window_size,
        "initialize_window produced wrong size"
    );
    (left, right)
}

// ============================================================================
// Interpolation
// ============================================================================

/// Linearly interpolate smoothed values between fitted points.
pub fn interpolate_gap<T: Float>(x: &[T], y_smooth: &mut [T], last: usize, current: usize) {
    if current <= last + 1 {
        return;
    }
    let denom = x[current] - x[last];
    if denom <= T::zero() {
        let val = y_smooth[last];
        for ys in y_smooth.iter_mut().take(current).skip(last + 1) {
            *ys = val;
        }
        return;
    }

    let y_current = y_smooth[current];
    let y_last = y_smooth[last];

    let xs = &x[(last + 1)..current];
    for (xj, ys) in xs.iter().zip(y_smooth.iter_mut().skip(last + 1)) {
        let alpha = (*xj - x[last]) / denom;
        *ys = alpha * y_current + (T::one() - alpha) * y_last;
    }
}

// ============================================================================
// Delta Optimization
// ============================================================================

/// Skip points that are within delta distance of the last fitted point.
pub fn skip_close_points<T: Float>(
    x: &[T],
    y_smooth: &mut [T],
    delta: T,
    last_fitted: &mut usize,
    n: usize,
) -> usize {
    let mut last = *last_fitted;

    // No delta optimization -> advance one step (cap to last index).
    if delta <= T::zero() {
        return usize::min(last.saturating_add(1), n.saturating_sub(1));
    }

    // If at or beyond the end, return final index.
    if last + 1 >= n {
        return n.saturating_sub(1);
    }

    // 1) Handle contiguous identical-x ties immediately after `last`.
    // Copy the last fitted value into tied neighbors and advance `last_fitted`
    // to the end of that tie-block so they are treated as processed.
    let mut tie_end = last;
    while tie_end + 1 < n && x[tie_end + 1] == x[last] {
        tie_end += 1;
    }
    if tie_end > last {
        for i in (last + 1)..=tie_end {
            y_smooth[i] = y_smooth[last];
        }
        *last_fitted = tie_end;
        last = *last_fitted;

        if last + 1 >= n {
            return n.saturating_sub(1);
        }
    }

    // 2) If any later point already has a fitted value, jump to it.
    for (i, ys) in y_smooth.iter().enumerate().skip(last + 1) {
        if *ys != T::zero() {
            return i;
        }
    }

    // 3) Apply delta-based skip logic from the updated `last`.
    let cut = x[last] + delta;
    let slice = &x[(last + 1)..n];
    let pos = slice.partition_point(|&xj| xj <= cut);

    // Next index should be the first index beyond those <= cut:
    let mut next = last.saturating_add(pos).saturating_add(1);

    // Cap to last valid index
    if next > n.saturating_sub(1) {
        next = n.saturating_sub(1);
    }

    // Ensure we advance at least one step and never go out of bounds.
    let candidate = usize::max(last.saturating_add(1), next);
    usize::min(candidate, n.saturating_sub(1))
}

// ============================================================================
// Weight Normalization
// ============================================================================

/// Normalize weights to sum to 1 over a specified range.
pub fn normalize_weights<T: Float>(weights: &mut [T], left: usize, right: usize, sum: T) {
    debug_assert!(sum > T::zero(), "Sum must be positive for normalization");
    let inv = T::one() / sum;
    for w in weights.iter_mut().take(right + 1).skip(left) {
        *w = *w * inv;
    }
}

/// Normalize entire weight array to sum to 1.
pub fn normalize_all_weights<T: Float>(weights: &mut [T]) {
    let sum: T = weights.iter().copied().fold(T::zero(), |acc, v| acc + v);
    if sum > T::zero() {
        let inv = T::one() / sum;
        for w in weights.iter_mut() {
            *w = *w * inv;
        }
    } else if !weights.is_empty() {
        let uniform = T::one() / T::from(weights.len()).unwrap();
        weights.fill(uniform);
    }
}

// ============================================================================
// Range and Distance Utilities
// ============================================================================

/// Compute the range of x values.
pub fn compute_range<T: Float>(x: &[T]) -> T {
    if x.is_empty() {
        T::zero()
    } else {
        x[x.len() - 1] - x[0]
    }
}

/// Find index of rightmost point within a distance threshold.
pub fn find_rightmost_point<T: Float>(
    x: &[T],
    x_current: T,
    left: usize,
    n: usize,
    threshold: T,
) -> usize {
    let target = x_current + threshold;
    x[left..n]
        .partition_point(|&xj| xj <= target)
        .saturating_sub(1)
        + left
}

// ============================================================================
// Statistical Utilities
// ============================================================================

/// Compute weighted average over a range.
pub fn compute_weighted_average<T: Float>(
    values: &[T],
    weights: &[T],
    left: usize,
    right: usize,
) -> T {
    let mut sum = T::zero();
    for j in left..=right {
        sum = sum + weights[j] * values[j];
    }
    sum
}

/// Check if a value is effectively zero.
pub fn is_effectively_zero<T: Float>(value: T, scale: T) -> bool {
    if scale <= T::zero() {
        return value == T::zero();
    }
    // Machine-epsilon based term (keeps behavior proportional to precision)
    let multiplier = T::from(8.0).unwrap();
    let eps_term = T::epsilon() * multiplier;

    // Practical absolute floor:
    // ~1e-9 relative to scale for f64, while f32 will use its larger eps_term.
    let abs_floor = T::from(1e-9).unwrap();

    let base = if eps_term > abs_floor {
        eps_term
    } else {
        abs_floor
    };
    let threshold = base * scale;
    value.abs() <= threshold
}