//! Kernel (weight) functions for LOWESS smoothing.
//!
//! ## Purpose
//!
//! This module provides kernel functions that define how distance-based
//! weights are computed in local regression. The kernel controls how
//! neighboring points influence the fit at each target point, with closer
//! points typically receiving higher weights.
//!
//! ## Design notes
//!
//! * All kernels map normalized distances u = |x - x_i| / bandwidth to weights.
//! * Bounded kernels have compact support on [-1, 1] for computational efficiency.
//! * The Gaussian kernel has unbounded support but uses a cutoff for numerical stability.
//! * Kernel properties (variance, roughness, efficiency) are precomputed for performance.
//! * The tricube kernel is the default, following Cleveland's original LOWESS.
//! * All weight computations are generic over `Float` types (f32/f64).
//!
//! ## Available kernels
//!
//! * **Tricube** (default): Smooth, efficient, recommended for most use cases
//! * **Epanechnikov**: Theoretically optimal (AMISE-minimal)
//! * **Gaussian**: Infinitely smooth, unbounded support
//! * **Biweight**: Good balance between smoothness and efficiency
//! * **Cosine**: Smooth with compact support
//! * **Triangle**: Simple linear taper
//! * **Uniform**: Equal weights (simplest, least smooth)
//!
//! ## Key concepts
//!
//! ### Normalized Distance
//! Distances are normalized by bandwidth: `u = |x - x_i| / h`
//! This makes kernels scale-invariant and allows comparison across different
//! window sizes.
//!
//! ### Kernel Properties
//! Each kernel has mathematical properties that affect smoothing quality:
//! - **Variance (μ₂)**: Second moment, affects bias-variance tradeoff
//! - **Roughness (R)**: Integral of K², lower is smoother
//! - **Efficiency**: Relative AMISE performance (Epanechnikov = 1.0)
//!
//! ### Support
//! Bounded kernels return zero outside [-1, 1], enabling efficient computation.
//! The Gaussian kernel has infinite support but uses a cutoff at u = 6.
//!
//! ## Invariants
//!
//! * All kernels are non-negative: K(u) ≥ 0.
//! * Bounded kernels satisfy K(u) = 0 for |u| ≥ 1.
//! * All kernels are symmetric: K(u) = K(-u).
//! * Kernel properties are constant and precomputed.
//!
//! ## Non-goals
//!
//! * This module does not perform normalization of weights (done by caller).
//! * This module does not validate input distances or bandwidths.
//! * This module does not handle adaptive bandwidths.
//!
//! ## Visibility
//!
//! The [`WeightFunction`] enum is part of the public API. Users can select
//! different kernels to control smoothing behavior. The internal kernel
//! properties and computation details are implementation details.

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

use core::f64::consts::{PI, SQRT_2};
use num_traits::Float;

// ============================================================================
// Mathematical Constants
// ============================================================================

/// Square root of 2π, used in Gaussian kernel calculations.
const SQRT_2PI: f64 = 2.5066282746310005024157652848110452530069867406099_f64;

/// Square root of π, used in kernel property calculations.
const SQRT_PI: f64 = 1.772453850905516027298167483341145182797_f64;

/// π/2, used in cosine kernel calculations.
const PI_OVER_2: f64 = PI / 2.0;

/// Cutoff for Gaussian kernel evaluation.
///
/// Beyond this normalized distance, the Gaussian kernel value is effectively
/// zero (exp(-6²/2) ≈ 6.9e-9). This prevents numerical underflow and improves
/// performance.
const GAUSSIAN_CUTOFF: f64 = 6.0;

// ============================================================================
// Kernel Properties
// ============================================================================

/// Mathematical properties of a kernel function.
///
/// These properties are used for efficiency calculations and theoretical
/// analysis of kernel performance.
///
/// # Fields
///
/// * `integrator` - c(K) = ∫ K(u) du, the kernel integral
/// * `variance` - μ₂(K) = ∫ u² K(u) du, the unnormalized second moment
/// * `roughness` - R(K) = ∫ K(u)² du, measures kernel smoothness
struct KernelProperties {
    /// Kernel integral: c(K) = ∫ K(u) du
    integrator: f64,

    /// Unnormalized second moment: μ₂(K) = ∫ u² K(u) du
    variance: f64,

    /// Kernel roughness: R(K) = ∫ K(u)² du
    roughness: f64,
}

/// Precomputed properties for the Cosine kernel.
const COSINE_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 4.0 / PI,
    variance: 4.0 / PI - 32.0 / (PI * PI * PI),
    roughness: 1.0,
};

/// Precomputed properties for the Epanechnikov kernel.
const EPANECHNIKOV_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 4.0 / 3.0,
    variance: 4.0 / 15.0,
    roughness: 16.0 / 15.0,
};

/// Precomputed properties for the Gaussian kernel.
const GAUSSIAN_PROPERTIES: KernelProperties = KernelProperties {
    integrator: SQRT_2PI,
    variance: SQRT_2 * SQRT_PI,
    roughness: SQRT_PI,
};

/// Precomputed properties for the Biweight kernel.
const BIWEIGHT_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 16.0 / 15.0,
    variance: 16.0 / 105.0,
    roughness: 256.0 / 315.0,
};

/// Precomputed properties for the Triangular kernel.
const TRIANGULAR_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 1.0,
    variance: 1.0 / 6.0,
    roughness: 2.0 / 3.0,
};

/// Precomputed properties for the Tricube kernel.
const TRICUBE_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 81.0 / 70.0,
    variance: 1.0 / 6.0,
    roughness: 6_561.0 / 6_916.0,
};

/// Precomputed properties for the Uniform kernel.
const UNIFORM_PROPERTIES: KernelProperties = KernelProperties {
    integrator: 2.0,
    variance: 2.0 / 3.0,
    roughness: 2.0,
};

// ============================================================================
// Weight Function Enum
// ============================================================================

/// Weight function (kernel) for LOWESS smoothing.
///
/// Each kernel defines a function K: ℝ → [0, ∞) that maps normalized
/// distances to weights. Bounded kernels have support on [-1, 1], while
/// the Gaussian kernel has unbounded support.
///
/// # Mathematical Properties
///
/// | Kernel       | Formula        | Efficiency† | R(K)            | μ₂(K)        |
/// |--------------|----------------|-------------|-----------------|--------------|
/// | Cosine       | `cos(πu/2)`    | 0.9995      | `π²/16`         | `1 - 8/π²`   |
/// | Epanechnikov | `(1 - u²)`     | 1.0000      | `3/5`           | `1/5`        |
/// | Gaussian     | `exp(-u²/2)`   | 0.9607      | `1/(2√π)`       | `1`          |
/// | Biweight     | `(1 - u²)²`    | 0.9951      | `5/7`           | `1/7`        |
/// | Triangular   | `(1 - |u|)`    | 0.9887      | `2/3`           | `1/6`        |
/// | Tricube      | `(1 - |u|³)³`  | 0.9983      | `1225/1729`     | `35/243`     |
/// | Uniform      | `1`            | 0.9432      | `1/2`           | `1/3`        |
///
/// † Relative AMISE efficiency (Epanechnikov = 1.0) from Wand & Jones (1995)
///
/// **Note**: The μ₂(K) column lists the unnormalized second moment.
/// Use the `integrator()` method to obtain the normalized moment when needed.
///
/// # Choosing a Kernel
///
/// * **Tricube** (default): Best all-around choice
///   - High efficiency (0.9983)
///   - Smooth derivatives
///   - Compact support (computationally efficient)
///   - Cleveland's original choice
///
/// * **Epanechnikov**: Theoretically optimal
///   - AMISE-optimal for kernel density estimation
///   - Less smooth than tricube
///   - Efficiency = 1.0 by definition
///
/// * **Gaussian**: Maximum smoothness
///   - Infinitely smooth
///   - No boundary effects
///   - More expensive to compute
///   - Good for very smooth data
///
/// * **Biweight**: Good balance
///   - High efficiency (0.9951)
///   - Smoother than Epanechnikov
///   - Compact support
///
/// * **Cosine**: Smooth and compact
///   - Good for robust smoothing contexts
///   - High efficiency (0.9995)
///
/// * **Triangle**: Simple and fast
///   - Linear taper
///   - Less smooth than other kernels
///   - Easy to understand
///
/// * **Uniform**: Simplest
///   - Equal weights within window
///   - Fastest to compute
///   - Least smooth results
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum WeightFunction {
    /// Cosine kernel: cos(πu/2) for |u| < 1
    ///
    /// Smooth, compact support, used in some robust smoothing contexts.
    /// Efficiency: 0.9995
    Cosine,

    /// Epanechnikov kernel: (1 - u²) for |u| < 1
    ///
    /// Optimal in terms of mean squared error for kernel density estimation.
    /// Faster to compute than tricube but less smooth.
    /// Efficiency: 1.0 (by definition)
    Epanechnikov,

    /// Gaussian kernel: exp(-u²/2)
    ///
    /// Infinitely smooth with unbounded support. More computationally
    /// expensive but provides very smooth results.
    /// Efficiency: 0.9607
    Gaussian,

    /// Biweight (quartic) kernel: (1 - u²)² for |u| < 1
    ///
    /// Good balance between smoothness and computational efficiency.
    /// Efficiency: 0.9951
    Biweight,

    /// Triangular (linear) kernel: (1 - |u|) for |u| < 1
    ///
    /// Simple linear taper. Fast and intuitive.
    /// Efficiency: 0.9887
    Triangle,

    /// Tricube kernel: (1 - |u|³)³ for |u| < 1
    ///
    /// This is the default and recommended weight function.
    /// It provides smooth, continuous derivatives and good local adaptation.
    /// Used in Cleveland's original LOWESS implementation.
    /// Efficiency: 0.9983
    #[default]
    Tricube,

    /// Uniform (rectangular) kernel: 1 for |u| < 1
    ///
    /// Simplest weight function, giving equal weight to all neighbors.
    /// Fastest to compute but produces less smooth results.
    /// Efficiency: 0.9432
    Uniform,
}

impl WeightFunction {
    // ========================================================================
    // Metadata Methods
    // ========================================================================

    /// Get the name of the weight function.
    ///
    /// # Returns
    ///
    /// A static string with the kernel name.
    #[inline]
    pub const fn name(&self) -> &'static str {
        match self {
            WeightFunction::Cosine => "Cosine",
            WeightFunction::Epanechnikov => "Epanechnikov",
            WeightFunction::Gaussian => "Gaussian",
            WeightFunction::Biweight => "Biweight",
            WeightFunction::Triangle => "Triangle",
            WeightFunction::Tricube => "Tricube",
            WeightFunction::Uniform => "Uniform",
        }
    }

    /// Get the kernel properties.
    ///
    /// # Returns
    ///
    /// Reference to the precomputed kernel properties.
    const fn properties(&self) -> &'static KernelProperties {
        match self {
            WeightFunction::Cosine => &COSINE_PROPERTIES,
            WeightFunction::Epanechnikov => &EPANECHNIKOV_PROPERTIES,
            WeightFunction::Gaussian => &GAUSSIAN_PROPERTIES,
            WeightFunction::Biweight => &BIWEIGHT_PROPERTIES,
            WeightFunction::Triangle => &TRIANGULAR_PROPERTIES,
            WeightFunction::Tricube => &TRICUBE_PROPERTIES,
            WeightFunction::Uniform => &UNIFORM_PROPERTIES,
        }
    }

    // ========================================================================
    // Kernel Property Accessors
    // ========================================================================

    /// Get the unnormalized variance (second moment) of the kernel.
    ///
    /// This is μ₂(K) = ∫ u² K(u) du.
    ///
    /// # Returns
    ///
    /// The unnormalized second moment.
    #[inline]
    pub fn variance(&self) -> f64 {
        self.properties().variance
    }

    /// Get the roughness of the kernel.
    ///
    /// This is R(K) = ∫ K(u)² du.
    ///
    /// Lower roughness generally indicates a smoother kernel.
    ///
    /// # Returns
    ///
    /// The kernel roughness.
    #[inline]
    pub fn roughness(&self) -> f64 {
        self.properties().roughness
    }

    /// Get the kernel integrator.
    ///
    /// This is c(K) = ∫ K(u) du.
    ///
    /// For normalized kernels that integrate to 1, this equals 1.
    /// For unnormalized kernels, this is the normalization constant.
    ///
    /// # Returns
    ///
    /// The kernel integral.
    #[inline]
    pub fn integrator(&self) -> f64 {
        self.properties().integrator
    }

    /// AMISE relative efficiency (Epanechnikov = 1.0).
    ///
    /// This measures how efficient the kernel is for density estimation
    /// in terms of asymptotic mean integrated squared error (AMISE).
    /// Higher values are better.
    ///
    /// # Formula
    ///
    /// Calculated as:
    /// ```text
    /// (R(K_E)/R(K))^(4/5) * (μ₂(K_E)/μ₂(K))^(2/5) * (c(K)/c(K_E))²
    /// ```
    ///
    /// where K_E is the Epanechnikov kernel (the AMISE-optimal kernel).
    ///
    /// # Returns
    ///
    /// The relative efficiency (Epanechnikov = 1.0).
    #[inline]
    pub fn efficiency(&self) -> f64 {
        let stats = self.properties();
        let ep_stats = &EPANECHNIKOV_PROPERTIES;

        let r_ratio = ep_stats.roughness / stats.roughness;
        let v_ratio = ep_stats.variance / stats.variance;
        let c_ratio = stats.integrator / ep_stats.integrator;

        r_ratio.powf(4.0 / 5.0) * v_ratio.powf(2.0 / 5.0) * c_ratio.powf(2.0)
    }

    // ========================================================================
    // Support Methods
    // ========================================================================

    /// Returns the support interval for bounded kernels.
    ///
    /// # Returns
    ///
    /// * `None` for unbounded kernels (Gaussian)
    /// * `Some((-1.0, 1.0))` for all bounded kernels
    #[inline]
    pub fn support(&self) -> Option<(f64, f64)> {
        match self {
            WeightFunction::Gaussian => None, // Unbounded
            _ => Some((-1.0, 1.0)),           // All others are bounded on [-1, 1]
        }
    }

    /// Returns `true` if the kernel has bounded support.
    ///
    /// # Returns
    ///
    /// `true` for all kernels except Gaussian, `false` for Gaussian.
    #[inline]
    fn is_bounded(&self) -> bool {
        self.support().is_some()
    }

    // ========================================================================
    // Weight Computation
    // ========================================================================

    /// Compute the unnormalized weight for a given normalized distance.
    ///
    /// # Parameters
    ///
    /// * `u` - Normalized distance: u = |x - x_i| / bandwidth
    ///
    /// # Returns
    ///
    /// Weight value K(u). For bounded kernels, returns 0 if |u| ≥ 1.
    ///
    /// # Algorithm
    ///
    /// 1. Compute |u| once
    /// 2. For bounded kernels, early-return 0 if |u| ≥ 1
    /// 3. Evaluate the kernel formula
    #[inline]
    pub fn compute_weight<T: Float>(&self, u: T) -> T {
        let abs_u = u.abs();

        // Fast path for bounded kernels: return 0 if outside support
        if self.is_bounded() && abs_u >= T::one() {
            return T::zero();
        }

        match self {
            WeightFunction::Cosine => {
                let pi_over_2 = T::from(PI_OVER_2).unwrap();
                (pi_over_2 * abs_u).cos()
            }

            WeightFunction::Epanechnikov => T::one() - abs_u * abs_u,

            WeightFunction::Gaussian => {
                // Convert to f64 for exponential calculation
                let u_f64 = abs_u.to_f64().unwrap_or(f64::INFINITY);

                // Use cutoff to avoid underflow to zero
                if u_f64 > GAUSSIAN_CUTOFF {
                    T::from(f64::MIN_POSITIVE).unwrap_or_else(T::zero)
                } else {
                    let val = (-0.5 * u_f64 * u_f64).exp().max(f64::MIN_POSITIVE);
                    T::from(val).unwrap_or_else(T::zero)
                }
            }

            WeightFunction::Biweight => {
                let tmp = T::one() - abs_u * abs_u;
                tmp * tmp
            }

            WeightFunction::Triangle => T::one() - abs_u,

            WeightFunction::Tricube => {
                let tmp = T::one() - abs_u * abs_u * abs_u;
                tmp * tmp * tmp
            }

            WeightFunction::Uniform => T::one(),
        }
    }

    /// Apply the kernel to a window of points.
    ///
    /// This method computes unnormalized weights for all points in a window.
    /// It is the core kernel weighting logic used by regression fitters.
    ///
    /// # Parameters
    ///
    /// * `x` - Sorted array of x-coordinates
    /// * `left` - Start index of window (inclusive)
    /// * `right` - End index of window (inclusive)
    /// * `x_current` - Center x-coordinate for the current fit
    /// * `bandwidth` - Bandwidth for distance normalization
    /// * `h1` - Distance threshold for weight 1.0 (very close points)
    /// * `h9` - Distance threshold for inclusion (usually ~1.0 * bandwidth)
    /// * `weights` - Output buffer for computed weights (will be modified)
    ///
    /// # Returns
    ///
    /// Tuple of `(unnormalized_weight_sum, rightmost_included_index)`:
    /// - `unnormalized_weight_sum`: Sum of all computed weights
    /// - `rightmost_included_index`: Index of the rightmost point with non-zero weight
    ///
    /// # Algorithm
    ///
    /// 1. Handle degenerate bandwidth (≤ 0) by zeroing all weights
    /// 2. Skip points to the left of `x_current - h9` for efficiency
    /// 3. Compute weights for points in range `[x_current - h9, x_current + h9]`
    /// 4. Early exit when exceeding `h9` on the right (data is sorted)
    /// 5. Track sum of weights and rightmost included index
    #[allow(clippy::too_many_arguments)]
    pub fn compute_window_weights<T: Float>(
        &self,
        x: &[T],
        left: usize,
        right: usize,
        x_current: T,
        bandwidth: T,
        h1: T,
        h9: T,
        weights: &mut [T],
    ) -> (T, usize) {
        let n = x.len();

        // Degenerate bandwidth: zero all weights in window
        if bandwidth <= T::zero() {
            for w in weights.iter_mut().take(n).skip(left) {
                *w = T::zero();
            }
            return (T::zero(), left);
        }

        let mut sum = T::zero();
        let mut rightmost = left;

        // Skip points to the left of (x_current - h9) for efficiency
        let lower_bound = x_current - h9;
        let mut start = left;
        while start < n && x[start] < lower_bound {
            start += 1;
        }

        // Zero the skipped region [left..start)
        for w in weights.iter_mut().take(start).skip(left) {
            *w = T::zero();
        }

        // Compute weights from start onwards
        for j in start..=right {
            let xj = x[j];
            let distance = (xj - x_current).abs();

            if distance > h9 {
                if xj >= x_current {
                    // Beyond h9 on right side (x is sorted): zero remaining in window and break
                    for w in weights.iter_mut().take(right + 1).skip(j) {
                        *w = T::zero();
                    }
                    break;
                }
                // Beyond h9 on left side (defensive)
                weights[j] = T::zero();
                continue;
            }

            // Compute weight: use 1.0 for very close points, otherwise evaluate kernel
            let w_k = if distance <= h1 {
                T::one()
            } else {
                let u = distance / bandwidth;
                self.compute_weight(u)
            };

            weights[j] = w_k;
            sum = sum + w_k;
            rightmost = j;
        }

        (sum, rightmost)
    }
}
