//! Confidence and prediction intervals for LOWESS smoothing.
//!
//! ## Purpose
//!
//! This module provides tools for quantifying uncertainty in LOWESS smoothed
//! values through standard error estimation, confidence intervals, and
//! prediction intervals. These allow users to assess the reliability of
//! fitted values and make probabilistic predictions.
//!
//! ## Design notes
//!
//! * Standard errors are computed using local leverage and weighted residuals.
//! * Confidence intervals quantify uncertainty in the fitted curve (mean).
//! * Prediction intervals quantify uncertainty for new observations.
//! * Uses MAD (Median Absolute Deviation) for robust scale estimation.
//! * Z-scores are approximated using Acklam's inverse normal CDF algorithm.
//! * All computations are generic over `Float` types to support f32 and f64.
//!
//! ## Available interval types
//!
//! * **Standard errors**: Local variance estimates for each smoothed point
//! * **Confidence intervals**: Uncertainty in the fitted curve (mean function)
//! * **Prediction intervals**: Uncertainty for new observations (includes residual variance)
//!
//! ## Key concepts
//!
//! ### Standard Errors
//! Computed from local leverage and weighted residual variance. Reflects
//! the uncertainty in the fitted value at each point due to sampling
//! variability.
//!
//! ### Confidence Intervals
//! Constructed as ŷ ± z × SE, where z is the critical value for the
//! specified confidence level. These intervals contain the true mean
//! function with the specified probability.
//!
//! ### Prediction Intervals
//! Constructed as ŷ ± z × √(SE² + σ²), where σ is the residual standard
//! deviation. These intervals account for both estimation uncertainty
//! and inherent data variability.
//!
//! ### Leverage
//! The influence of each observation on its own fitted value. Computed
//! as the ratio of the observation's weight to the sum of all weights
//! in its local window.
//!
//! ### Z-Score Approximation
//! Uses Acklam's algorithm for accurate inverse normal CDF computation,
//! with fast paths for common confidence levels (90%, 95%, 99%).
//!
//! ## Invariants
//!
//! * Confidence level must be in (0, 1).
//! * Prediction intervals are always wider than confidence intervals.
//! * Standard errors are non-negative.
//! * Interval widths are non-negative and finite.
//!
//! ## Non-goals
//!
//! * This module does not perform the smoothing itself.
//! * This module does not validate input data.
//! * This module does not provide bootstrap or permutation-based intervals.
//! * This module does not handle simultaneous confidence bands.
//!
//! ## Visibility
//!
//! The [`IntervalMethod`] struct is part of the public API, allowing users
//! to configure which intervals to compute. Internal implementation details
//! may change without notice.

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

use crate::math::mad::compute_mad;
use crate::primitives::errors::LowessError;
use crate::primitives::window::Window;
use num_traits::Float;

// ============================================================================
// Interval Configuration
// ============================================================================

/// Configuration for computing confidence/prediction intervals and standard errors.
///
/// Encapsulates both the configuration (which intervals to compute) and the
/// logic for performing the computations. Users can selectively enable
/// confidence intervals, prediction intervals, and/or standard errors.
///
/// # Fields
///
/// * `level` - Confidence level (e.g., 0.95 for 95% intervals)
/// * `confidence` - Whether to compute confidence intervals for the mean
/// * `prediction` - Whether to compute prediction intervals for new observations
/// * `se` - Whether to compute standard errors
///
/// # Relationship between fields
///
/// * If `confidence` or `prediction` is true, `se` is automatically enabled
/// * Standard errors can be computed independently without intervals
/// * The `level` field is only used when computing intervals
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct IntervalMethod<T> {
    /// Confidence level (e.g., 0.95 for 95% intervals)
    pub level: T,

    /// Compute confidence intervals for the fitted curve
    pub confidence: bool,

    /// Compute prediction intervals for new observations
    pub prediction: bool,

    /// Compute standard errors
    pub se: bool,
}

impl<T: Float> Default for IntervalMethod<T> {
    fn default() -> Self {
        Self::none()
    }
}

impl<T: Float> IntervalMethod<T> {
    // ========================================================================
    // Constructors
    // ========================================================================

    /// No intervals or standard errors.
    ///
    /// Carries a default confidence level (0.95) for potential future use.
    fn none() -> Self {
        Self {
            level: T::from(0.95).unwrap(),
            confidence: false,
            prediction: false,
            se: false,
        }
    }

    /// Confidence intervals only at the specified level.
    ///
    /// Standard errors are automatically enabled.
    ///
    /// # Parameters
    ///
    /// * `level` - Confidence level (e.g., 0.95 for 95% intervals)
    pub fn confidence(level: T) -> Self {
        Self {
            level,
            confidence: true,
            prediction: false,
            se: true,
        }
    }

    /// Prediction intervals only at the specified level.
    ///
    /// Standard errors are automatically enabled.
    ///
    /// # Parameters
    ///
    /// * `level` - Confidence level (e.g., 0.95 for 95% intervals)
    pub fn prediction(level: T) -> Self {
        Self {
            level,
            confidence: false,
            prediction: true,
            se: true,
        }
    }

    /// Standard errors only (no intervals).
    ///
    /// The level field defaults to 0.95 but is not used.
    pub fn se() -> Self {
        Self {
            level: T::from(0.95).unwrap(),
            confidence: false,
            prediction: false,
            se: true,
        }
    }
}

impl<T: Float> IntervalMethod<T> {
    // ========================================================================
    // Constants
    // ========================================================================

    /// Constant to convert MAD to an unbiased estimate of σ for normal data.
    ///
    /// For normally distributed data, MAD × 1.4826 ≈ standard deviation.
    const MAD_TO_STD_FACTOR: f64 = 1.4826;

    /// Minimum tuned-scale absolute epsilon to avoid division by zero.
    const MIN_TUNED_SCALE: f64 = 1e-12;

    /// Number of parameters in local linear regression (intercept + slope).
    const LINEAR_PARAMS: f64 = 2.0;

    // ========================================================================
    // Robust Scale Estimation
    // ========================================================================

    /// Estimate residual standard deviation using MAD.
    ///
    /// Uses MAD (Median Absolute Deviation) scaled by 1.4826 to provide
    /// a robust estimate of the residual standard deviation.
    ///
    /// # Parameters
    ///
    /// * `residuals` - Residual values (y - ŷ)
    ///
    /// # Returns
    ///
    /// Estimated residual standard deviation.
    fn calculate_residual_sd(residuals: &[T]) -> T {
        let n = residuals.len();
        if n == 0 {
            return T::zero();
        }

        let scale_const = T::from(Self::MAD_TO_STD_FACTOR).unwrap();

        if n == 1 {
            return residuals[0].abs() * scale_const;
        }

        let mad = compute_mad(residuals);
        if mad > T::zero() {
            mad * scale_const
        } else {
            // Apply minimum scale to avoid division by zero
            let min_eps = T::from(Self::MIN_TUNED_SCALE).unwrap();
            min_eps * scale_const
        }
    }

    // ========================================================================
    // Standard Error Computation
    // ========================================================================

    /// Core mathematical function for computing standard error.
    ///
    /// Computes the standard error for a single point using weighted
    /// residual variance and leverage.
    ///
    /// # Parameters
    ///
    /// * `sum_w` - Sum of weights in the local window
    /// * `sum_w_r2` - Sum of weighted squared residuals
    /// * `w_idx` - Weight of the current point
    ///
    /// # Returns
    ///
    /// Standard error at the current point.
    ///
    /// # Formula
    ///
    /// SE = √(variance × leverage)
    ///
    /// where:
    /// - variance = sum_w_r2 / (sum_w - 2)
    /// - leverage = w_idx / sum_w
    fn compute_se(sum_w: T, sum_w_r2: T, w_idx: T) -> T {
        // Effective degrees of freedom for weighted regression
        if sum_w <= T::zero() {
            return T::zero();
        }

        let effective_n = sum_w;
        let df = effective_n - T::from(Self::LINEAR_PARAMS).unwrap();

        if df <= T::zero() {
            return T::zero();
        }

        let variance = sum_w_r2 / df;
        let leverage = w_idx / sum_w; // Normalized leverage

        (variance * leverage).sqrt()
    }

    /// Compute standard errors for all smoothed values.
    ///
    /// For each point, constructs a local window, computes weighted residuals,
    /// and estimates the standard error using local leverage.
    ///
    /// # Parameters
    ///
    /// * `x` - Sorted x-coordinates
    /// * `y` - Original y-values
    /// * `y_smooth` - Smoothed y-values
    /// * `window_size` - Size of the local window
    /// * `robustness_weights` - Robustness weights from iterative refinement
    /// * `std_errors` - Output buffer for standard errors (will be modified)
    /// * `weight_fn` - Kernel weight function
    ///
    /// # Algorithm
    ///
    /// For each point i:
    /// 1. Initialize and center a local window
    /// 2. Compute bandwidth from window extent
    /// 3. Compute weights for all points in window
    /// 4. Compute weighted residual variance
    /// 5. Compute leverage for point i
    /// 6. Compute SE using `compute_se`
    #[allow(clippy::too_many_arguments)]
    pub fn compute_window_se<F>(
        &self,
        x: &[T],
        y: &[T],
        y_smooth: &[T],
        window_size: usize,
        robustness_weights: &[T],
        std_errors: &mut [T],
        weight_fn: &F,
    ) where
        F: Fn(T) -> T,
    {
        // Early exit if no intervals or SE requested
        if !self.se && !self.confidence && !self.prediction {
            return;
        }

        let n = x.len();

        for (i, se) in std_errors.iter_mut().enumerate().take(n) {
            // Initialize and center window
            let mut window = Window::initialize(i, window_size, n);
            window.recenter(x, i, n);

            let idx = i;
            let left = window.left;
            let right = window.right;

            // Compute bandwidth
            let x_current = x[idx];
            let bandwidth_left = x_current - x[left];
            let bandwidth_right = x[right] - x_current;
            let bandwidth = T::max(bandwidth_left, bandwidth_right);

            if bandwidth <= T::zero() {
                *se = T::zero();
                continue;
            }

            // Compute weight for current point (distance = 0)
            let u_idx = T::zero();
            let kernel_val = weight_fn(u_idx);
            let w_idx = kernel_val * robustness_weights[idx];

            // Accumulate weighted residual variance
            let mut sum_w_r2 = T::zero();
            let mut sum_w = T::zero();

            for j in left..=right {
                let dist = (x[j] - x_current).abs();
                let u = dist / bandwidth;
                let w = if j == idx {
                    w_idx
                } else {
                    weight_fn(u) * robustness_weights[j]
                };

                let r = y[j] - y_smooth[j];
                sum_w_r2 = sum_w_r2 + w * r * r;
                sum_w = sum_w + w;
            }

            *se = Self::compute_se(sum_w, sum_w_r2, w_idx);
        }
    }

    // ========================================================================
    // Interval Computation
    // ========================================================================

    /// Compute confidence and/or prediction intervals.
    ///
    /// Based on the configuration, computes the requested interval types
    /// and validates that they are well-formed.
    ///
    /// # Parameters
    ///
    /// * `y_smooth` - Smoothed y-values
    /// * `std_errors` - Standard errors for each point
    /// * `residuals` - Residuals (y - ŷ)
    ///
    /// # Returns
    ///
    /// Tuple of `(conf_lower, conf_upper, pred_lower, pred_upper)`:
    /// - Each component is `Some(Vec<T>)` if requested, `None` otherwise
    ///
    /// # Errors
    ///
    /// Returns `LowessError::InvalidIntervals` if the confidence level is
    /// invalid or if interval computation fails.
    #[allow(clippy::type_complexity)]
    pub fn compute_intervals(
        &self,
        y_smooth: &[T],
        std_errors: &[T],
        residuals: &[T],
    ) -> Result<
        (
            Option<Vec<T>>, // confidence lower
            Option<Vec<T>>, // confidence upper
            Option<Vec<T>>, // prediction lower
            Option<Vec<T>>, // prediction upper
        ),
        LowessError,
    > {
        // Compute confidence intervals if requested
        let (mut conf_lower, mut conf_upper) = if self.confidence {
            let (lower, upper) = self
                .compute_confidence_intervals_impl(y_smooth, std_errors)
                .map_err(|_| LowessError::InvalidIntervals(self.level.to_f64().unwrap_or(0.0)))?;
            (Some(lower), Some(upper))
        } else {
            (None, None)
        };

        // Compute prediction intervals if requested
        let (mut pred_lower, mut pred_upper) = if self.prediction {
            let residual_sd = Self::calculate_residual_sd(residuals);
            let (lower, upper) = self
                .compute_prediction_intervals_impl(y_smooth, std_errors, residual_sd)
                .map_err(|_| LowessError::InvalidIntervals(self.level.to_f64().unwrap_or(0.0)))?;
            (Some(lower), Some(upper))
        } else {
            (None, None)
        };

        // Guard against degenerate intervals
        let residual_sd = Self::calculate_residual_sd(residuals);
        let any_std_nonzero = std_errors.iter().any(|&s| s > T::zero());

        if residual_sd > T::zero() || any_std_nonzero {
            let eps = T::from(1e-12).unwrap_or_else(|| T::from(1e-6).unwrap());

            // Fix degenerate confidence intervals
            if let (Some(lo), Some(hi)) = (&mut conf_lower, &mut conf_upper) {
                for (l, h) in lo.iter_mut().zip(hi.iter_mut()) {
                    let width = *h - *l;
                    if !width.is_finite() || width <= T::zero() {
                        *h = *l + eps;
                    }
                }
            }

            // Fix degenerate prediction intervals
            if let (Some(lo), Some(hi)) = (&mut pred_lower, &mut pred_upper) {
                for (l, h) in lo.iter_mut().zip(hi.iter_mut()) {
                    let width = *h - *l;
                    if !width.is_finite() || width <= T::zero() {
                        *h = *l + eps;
                    }
                }
            }
        }

        Ok((conf_lower, conf_upper, pred_lower, pred_upper))
    }

    /// Compute confidence intervals for the fitted curve.
    ///
    /// Confidence intervals quantify uncertainty in the mean function.
    ///
    /// # Parameters
    ///
    /// * `y_smooth` - Smoothed y-values
    /// * `std_errors` - Standard errors for each point
    ///
    /// # Returns
    ///
    /// Tuple of `(lower, upper)` bounds.
    ///
    /// # Formula
    ///
    /// CI = ŷ ± z × SE
    ///
    /// where z is the critical value for the specified confidence level.
    fn compute_confidence_intervals_impl(
        &self,
        y_smooth: &[T],
        std_errors: &[T],
    ) -> Result<(Vec<T>, Vec<T>), &'static str> {
        let z = Self::approximate_z_score(self.level)?;

        let lower: Vec<T> = y_smooth
            .iter()
            .zip(std_errors.iter())
            .map(|(&ys, &se)| ys - z * se)
            .collect();

        let upper: Vec<T> = y_smooth
            .iter()
            .zip(std_errors.iter())
            .map(|(&ys, &se)| ys + z * se)
            .collect();

        Ok((lower, upper))
    }

    /// Compute prediction intervals for future observations.
    ///
    /// Prediction intervals account for both estimation uncertainty and
    /// inherent data variability.
    ///
    /// # Parameters
    ///
    /// * `y_smooth` - Smoothed y-values
    /// * `std_errors` - Standard errors for each point
    /// * `residual_sd` - Residual standard deviation
    ///
    /// # Returns
    ///
    /// Tuple of `(lower, upper)` bounds.
    ///
    /// # Formula
    ///
    /// PI = ŷ ± z × √(SE² + σ²)
    ///
    /// where:
    /// - z is the critical value for the specified confidence level
    /// - σ is the residual standard deviation
    fn compute_prediction_intervals_impl(
        &self,
        y_smooth: &[T],
        std_errors: &[T],
        residual_sd: T,
    ) -> Result<(Vec<T>, Vec<T>), &'static str> {
        let z = Self::approximate_z_score(self.level)?;
        let rsd_sq = residual_sd * residual_sd;

        let lower: Vec<T> = y_smooth
            .iter()
            .zip(std_errors.iter())
            .map(|(&ys, &se)| {
                let pred_se = (se * se + rsd_sq).sqrt();
                ys - z * pred_se
            })
            .collect();

        let upper: Vec<T> = y_smooth
            .iter()
            .zip(std_errors.iter())
            .map(|(&ys, &se)| {
                let pred_se = (se * se + rsd_sq).sqrt();
                ys + z * pred_se
            })
            .collect();

        Ok((lower, upper))
    }

    // ========================================================================
    // Z-Score Approximation
    // ========================================================================

    /// Approximate the z-score (critical value) for a given confidence level.
    ///
    /// Uses fast paths for common confidence levels (90%, 95%, 99%) and
    /// Acklam's algorithm for other values.
    ///
    /// # Parameters
    ///
    /// * `confidence_level` - Confidence level in (0, 1)
    ///
    /// # Returns
    ///
    /// The z-score (critical value from standard normal distribution).
    ///
    /// # Errors
    ///
    /// Returns an error if the confidence level is not in (0, 1).
    ///
    /// # Common values
    ///
    /// * 0.90 → z ≈ 1.645
    /// * 0.95 → z ≈ 1.960
    /// * 0.99 → z ≈ 2.576
    pub fn approximate_z_score(confidence_level: T) -> Result<T, &'static str> {
        let cl_f = confidence_level.to_f64().unwrap_or(0.95);

        if !(cl_f > 0.0 && cl_f < 1.0) {
            return Err("confidence_level must be in (0, 1) (exclusive)");
        }

        // Convert confidence level to cumulative probability
        let p = (1.0 + cl_f) / 2.0;

        // Fast paths for common confidence levels
        let z = if (cl_f - 0.99).abs() < 1e-6 {
            2.576
        } else if (cl_f - 0.95).abs() < 1e-6 {
            1.960
        } else if (cl_f - 0.90).abs() < 1e-6 {
            1.645
        } else {
            // Use Acklam's algorithm for other values
            Self::acklam_inverse_cdf(p)
        };

        Ok(T::from(z).unwrap_or_else(|| T::one()))
    }

    /// Acklam's inverse normal CDF algorithm.
    ///
    /// Provides accurate approximation of the inverse cumulative distribution
    /// function of the standard normal distribution.
    ///
    /// # Parameters
    ///
    /// * `p` - Cumulative probability in (0, 1)
    ///
    /// # Returns
    ///
    /// The z-score corresponding to the cumulative probability.
    ///
    /// # Algorithm
    ///
    /// Uses rational approximations with different coefficients for three
    /// regions:
    /// - p < 0.02425: Lower tail
    /// - p > 0.97575: Upper tail
    /// - Otherwise: Central region
    ///
    /// # Reference
    ///
    /// Acklam, P. J. (2010). An algorithm for computing the inverse normal
    /// cumulative distribution function.
    fn acklam_inverse_cdf(p: f64) -> f64 {
        if p <= 0.0 || p >= 1.0 {
            return 0.0;
        }

        // Coefficients for central region
        const A: [f64; 6] = [
            -3.969_683_028_665_376e1,
            2.209_460_984_245_205e2,
            -2.759_285_104_469_687e2,
            1.383_577_518_672_69e2,
            -3.066_479_806_614_716e1,
            2.506_628_277_459_239e0,
        ];
        const B: [f64; 5] = [
            -5.447_609_879_822_406e1,
            1.615_858_368_580_409e2,
            -1.556_989_798_598_866e2,
            6.680_131_188_771_972e1,
            -1.328_068_155_288_572e1,
        ];

        // Coefficients for tail regions
        const C: [f64; 6] = [
            -7.784_894_002_430_293e-3,
            -3.223_964_580_411_365e-1,
            -2.400_758_277_161_838e0,
            -2.549_732_539_343_734e0,
            4.374_664_141_464_968e0,
            2.938_163_982_698_783e0,
        ];
        const D: [f64; 4] = [
            7.784_695_709_041_462e-3,
            3.224_671_290_700_398e-1,
            2.445_134_137_142_996e0,
            3.754_408_661_907_416e0,
        ];

        const P_LOW: f64 = 0.02425;
        const P_HIGH: f64 = 0.97575;

        if p < P_LOW {
            // Lower tail
            let q = (-2.0 * p.ln()).sqrt();
            let num = ((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5];
            let den = (((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0;
            -(num / den)
        } else if p > P_HIGH {
            // Upper tail
            let q = (-2.0 * (1.0 - p).ln()).sqrt();
            let num = ((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5];
            let den = (((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0;
            num / den
        } else {
            // Central region
            let q = p - 0.5;
            let r = q * q;
            let num = (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q;
            let den = ((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0;
            num / den
        }
    }
}
