//! Interpolation and delta optimization for LOWESS smoothing.
//!
//! ## Purpose
//!
//! This module provides utilities for optimizing LOWESS performance through
//! delta-based point skipping and linear interpolation. When data points are
//! densely sampled, fitting every point is computationally expensive and
//! often unnecessary. This module implements the "delta optimization" that
//! fits anchor points and interpolates between them.
//!
//! ## Design notes
//!
//! * Delta controls the distance threshold for point skipping during fitting.
//! * Points closer than delta may be interpolated rather than explicitly fitted.
//! * Uses linear interpolation to fill gaps between fitted anchor points.
//! * Handles tied x-values (duplicates) by copying fitted values.
//! * Delta defaults to 1% of the x-range when not specified.
//! * All operations are generic over `Float` types to support f32 and f64.
//!
//! ## Key concepts
//!
//! ### Delta Optimization
//! Instead of fitting every point, the algorithm:
//! 1. Selects "anchor" points spaced at least delta apart
//! 2. Explicitly fits these anchor points
//! 3. Linearly interpolates between anchors for intermediate points
//!
//! This provides significant speedup on dense data with minimal accuracy loss.
//!
//! ### Linear Interpolation
//! For points between two fitted anchors at (x₀, y₀) and (x₁, y₁), the
//! interpolated value at x is:
//! ```text
//! y = y₀ + α(y₁ - y₀)  where  α = (x - x₀)/(x₁ - x₀)
//! ```
//!
//! ### Tied Values
//! When multiple points have identical x-values, they all receive the same
//! fitted value (the value computed for that x-coordinate).
//!
//! ## Invariants
//!
//! * Input x-values must be sorted in ascending order.
//! * Delta must be non-negative and finite.
//! * Interpolation preserves monotonicity between anchor points.
//! * At least one point is always fitted (no infinite loops).
//!
//! ## Non-goals
//!
//! * This module does not perform the actual smoothing/fitting.
//! * This module does not sort the input data.
//! * This module does not provide higher-order interpolation (only linear).
//! * This module does not validate that x-values are sorted.
//!
//! ## Visibility
//!
//! This module is an internal implementation detail used by the LOWESS
//! engine. It is not part of the public API and may change without notice.

#[cfg(not(feature = "std"))]
extern crate alloc;

use crate::primitives::errors::LowessError;
use core::result::Result;
use num_traits::Float;

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

/// Calculate delta parameter for interpolation optimization.
///
/// Delta controls the distance threshold for point skipping during fitting.
/// Points closer than delta may be interpolated rather than explicitly fitted,
/// improving performance on dense data.
///
/// # Parameters
///
/// * `delta` - Optional delta value (must be non-negative and finite if provided)
/// * `x_sorted` - Sorted x-values
///
/// # Returns
///
/// Validated or computed delta value.
///
/// # Errors
///
/// Returns `LowessError::InvalidDelta` if the provided delta is negative,
/// infinite, or NaN.
///
/// # Default behavior
///
/// If delta is `None`, computes a conservative default as 1% of the x-range:
/// ```text
/// delta = 0.01 × (max(x) - min(x))
/// ```
///
/// This default provides a good balance between performance and accuracy
/// for most datasets.
///
/// # Algorithm
///
/// 1. If delta is provided:
///    - Validate that it is finite and non-negative
///    - Return the validated value
/// 2. If delta is `None`:
///    - Compute range as max(x) - min(x)
///    - Return 1% of the range
pub fn calculate_delta<T: Float>(delta: Option<T>, x_sorted: &[T]) -> Result<T, LowessError> {
    match delta {
        Some(d) => {
            // Validate provided delta
            if !d.is_finite() || d < T::zero() {
                return Err(LowessError::InvalidDelta(d.to_f64().unwrap_or(f64::NAN)));
            }
            Ok(d)
        }
        None => {
            // Compute default delta as 1% of x-range
            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)
            }
        }
    }
}

// ============================================================================
// Linear Interpolation
// ============================================================================

/// Interpolate gap between two fitted anchor points.
///
/// Uses linear interpolation to fill in smoothed values between `last_fitted`
/// and `current`. This is part of the delta optimization: when points are
/// close together, we fit anchor points and interpolate between them.
///
/// # Parameters
///
/// * `x` - Independent variable values (must be sorted)
/// * `y_smooth` - Smoothed values (will be modified in-place)
/// * `last_fitted` - Index of last explicitly fitted point
/// * `current` - Index of current fitted point
///
/// # Algorithm
///
/// For each index k ∈ (last_fitted, current):
/// 1. Compute interpolation weight: α = (x[k] - x₀) / (x₁ - x₀)
/// 2. Compute interpolated value: y[k] = y₀ + α(y₁ - y₀)
///
/// # Special cases
///
/// * **No gap**: If current ≤ last_fitted + 1, no interpolation is needed
/// * **Tied x-values**: If x₁ = x₀, uses simple average of y-values
/// * **Decreasing x**: Treated same as tied values (uses average)
///
/// # Preconditions
///
/// * Both `last_fitted` and `current` must have been explicitly fitted
/// * `x` must be sorted in ascending order
/// * `last_fitted < current`
///
/// # Complexity
///
/// * Time: O(current - last_fitted)
/// * Space: O(1)
pub fn interpolate_gap<T: Float>(x: &[T], y_smooth: &mut [T], last_fitted: usize, current: usize) {
    // No gap to interpolate
    if current <= last_fitted + 1 {
        return;
    }

    let x0 = x[last_fitted];
    let x1 = x[current];
    let y0 = y_smooth[last_fitted];
    let y1 = y_smooth[current];

    let denom = x1 - x0;

    if denom <= T::zero() {
        // Duplicate or decreasing x-values: use simple average
        let avg = (y0 + y1) / T::from(2.0).unwrap();
        y_smooth
            .iter_mut()
            .take(current)
            .skip(last_fitted + 1)
            .for_each(|ys| *ys = avg);
        return;
    }

    // Linear interpolation
    for k in (last_fitted + 1)..current {
        let xi = x[k];
        let alpha = (xi - x0) / denom;
        y_smooth[k] = y0 + alpha * (y1 - y0);
    }
}

// ============================================================================
// Point Skipping (Delta Optimization)
// ============================================================================

/// Skip points that are close together based on delta threshold.
///
/// This function implements the delta optimization: instead of fitting every
/// point, we "anchor" on sparse points and interpolate between them when
/// points are densely sampled.
///
/// # Parameters
///
/// * `x` - Sorted array of x-values
/// * `y_smooth` - Smoothed values (partially filled for fitted points)
/// * `delta` - Distance threshold for skipping (0 disables optimization)
/// * `last_fitted` - Mutable reference to last fitted index
/// * `n` - Total number of points
/// * `fitted_mask` - Optional boolean mask indicating which points are fitted
///
/// # Returns
///
/// Index of the next point to fit.
///
/// # Algorithm
///
/// 1. **Handle tied x-values**: If points immediately after `last_fitted` have
///    identical x-values, copy the fitted value to them and advance `last_fitted`.
///
/// 2. **Check for already-fitted points**: If `fitted_mask` is provided, find
///    the next point that's already marked as fitted. Otherwise, infer from
///    non-zero `y_smooth` values.
///
/// 3. **Apply delta-based skip**: Find the first point beyond `x[last] + delta`.
///    This skips over densely-sampled regions.
///
/// 4. **Ensure progress**: Always advance at least one step to prevent infinite loops.
///
/// # Fitted mask rationale
///
/// The `fitted_mask` parameter is necessary because we cannot reliably infer
/// whether a point has been fitted from `y_smooth` alone—a legitimate fitted
/// value might be exactly zero. When `fitted_mask` is `None`, the function
/// falls back to a heuristic (non-zero check) that may be unreliable.
///
/// # Special cases
///
/// * **Delta = 0**: Disables optimization, returns `last_fitted + 1`
/// * **At end**: Returns `n - 1` if already at or beyond the last point
/// * **Tied values**: Copies fitted value to all tied neighbors
///
/// # Complexity
///
/// * Time: O(log n) for delta-based skip (binary search), O(k) for tied values
/// * Space: O(1)
pub fn skip_close_points<T: Float>(
    x: &[T],
    y_smooth: &mut [T],
    delta: T,
    last_fitted: &mut usize,
    n: usize,
    fitted_mask: Option<&[bool]>,
) -> usize {
    // No delta optimization: advance one step
    if delta <= T::zero() {
        let idx = last_fitted.saturating_add(1);
        return if n == 0 { 0 } else { idx.min(n - 1) };
    }

    let mut last = *last_fitted;

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

    // Step 1: Handle contiguous tied x-values immediately after `last`
    // Copy the fitted value to all tied neighbors and advance last_fitted
    let x_last = x[last];
    let mut tie_end = last;

    while tie_end + 1 < n && x[tie_end + 1] == x_last {
        tie_end += 1;
        y_smooth[tie_end] = y_smooth[last];
    }

    if tie_end > last {
        *last_fitted = tie_end;
        last = tie_end;

        if last + 1 >= n {
            return n - 1;
        }
    }

    // Step 2: Find next already-fitted point (if any)
    if let Some(mask) = fitted_mask {
        debug_assert!(
            mask.len() >= n,
            "fitted_mask length must be at least n when provided"
        );

        // Search for next fitted point using mask
        let start = last + 1;
        for (i, &m) in mask.iter().enumerate().skip(start).take(n - start) {
            if m {
                return i;
            }
        }
    } else {
        // Infer fitted points from non-zero y_smooth values
        // NOTE: This is a fallback heuristic and may be unreliable if
        // legitimate fitted values can be zero
        let start = last + 1;
        for (i, &ys) in y_smooth.iter().enumerate().skip(start).take(n - start) {
            if ys != T::zero() {
                return i;
            }
        }
    }

    // Step 3: Apply delta-based skip logic
    // Find first point with x > x[last] + delta
    let cut = x[last] + delta;
    let slice = &x[(last + 1)..n];
    let pos = slice.partition_point(|&xj| xj <= cut);

    // next_candidate is the first index beyond those <= cut
    let next_candidate = last + 1 + pos;

    // Step 4: Ensure we advance at least one step and stay in bounds
    usize::max(last + 1, next_candidate).min(n - 1)
}
