dataprof-metrics 0.12.0

use crate::acceleration::simd::compute_stats_auto;
use crate::types::{ColumnStats, NumericStats, Quartiles};
use std::collections::HashMap;

const SAMPLE_THRESHOLD: usize = 10_000;

/// Seed for the order-statistics reservoir. Fixed so the same values always
/// yield the same median, quartiles, mode and shape statistics: a profiler
/// that answers differently on every run for identical input is not
/// deterministic, and the streaming reservoir is seeded for the same reason.
const SAMPLE_SEED: u64 = 0xDA7A_900D_F00D_5EED;

pub fn calculate_numeric_stats(data: &[String]) -> ColumnStats {
    ColumnStats::Numeric(compute_numeric_stats(data))
}

/// Compute numeric stats and return the inner struct directly.
pub fn compute_numeric_stats(data: &[String]) -> NumericStats {
    compute_numeric_stats_with_parsed_count(data).0
}

/// Compute numeric stats plus the number of values they cover — the count of
/// cells that parsed as a finite number. Lets callers derive invalid-value
/// counts without a second parse pass over the data.
pub fn compute_numeric_stats_with_parsed_count(data: &[String]) -> (NumericStats, usize) {
    // decode-audit: no-data — cells that do not parse are non-numeric values
    // (nulls, tokens like "N/A"), excluded from numeric stats by design.
    let numbers: Vec<f64> = data
        .iter()
        .filter_map(|s| s.trim().parse::<f64>().ok())
        .filter(|x| x.is_finite())
        .collect();
    let parsed_count = numbers.len();

    if numbers.is_empty() {
        return (NumericStats::empty(), 0);
    }

    // Single stable pass for the base statistics. `NumericAccumulator` carries
    // Welford's mean and M2 alongside a compensated sum, so the variance of a
    // small spread on a large offset survives and the mean survives both
    // cancellation and an overflowing sum.
    let base = compute_stats_auto(&numbers);
    let min = base.min().expect("accumulator covers a non-empty column");
    let max = base.max().expect("accumulator covers a non-empty column");
    let mean = base.mean();
    let raw_variance = base.sample_variance();
    let raw_std_dev = base.sample_std_dev();
    let variance = raw_variance.is_finite().then_some(raw_variance);
    let std_dev = raw_std_dev.is_finite().then_some(raw_std_dev);

    // Determine if we need sampling for large datasets
    let (sample_data, is_approximate) = if numbers.len() > SAMPLE_THRESHOLD {
        (reservoir_sample(&numbers, SAMPLE_THRESHOLD), Some(true))
    } else {
        (numbers.clone(), None)
    };

    // Statistics requiring sorted/complete data
    let mut sorted_data = sample_data.clone();
    sorted_data.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let median = calculate_median(&sorted_data);
    let quartiles = calculate_quartiles(&sorted_data);
    let mode = calculate_mode(&sample_data);
    let coefficient_of_variation =
        std_dev.and_then(|value| calculate_coefficient_of_variation(value, mean));

    // Advanced metrics
    let skewness = std_dev.and_then(|value| calculate_skewness(&sample_data, mean, value));
    let kurtosis = std_dev.and_then(|value| calculate_kurtosis(&sample_data, mean, value));

    // Per-column outlier count using the same Tukey IQR rule as the global
    // accuracy.outlier_ratio dimension (k = 1.5). Returns `None` for tiny
    // samples to mirror the global metric's behaviour.
    let outlier_count = count_outliers_iqr(&numbers, quartiles.as_ref(), 4);

    (
        NumericStats {
            min,
            max,
            mean,
            std_dev,
            variance,
            median,
            quartiles,
            mode,
            coefficient_of_variation,
            skewness,
            kurtosis,
            is_approximate,
            outlier_count,
        },
        parsed_count,
    )
}

/// Count IQR outliers given pre-computed quartiles. Returns `None` when the
/// sample is too small to produce a meaningful estimate.
fn count_outliers_iqr(
    values: &[f64],
    quartiles: Option<&Quartiles>,
    min_samples: usize,
) -> Option<usize> {
    if values.len() < min_samples {
        return None;
    }
    let q = quartiles?;
    // An IQR beyond the `f64` range puts both fences beyond it too, so no
    // finite value can fall outside them.
    let Some(iqr) = q.iqr else {
        return Some(0);
    };
    let lower = q.q1 - 1.5 * iqr;
    let upper = q.q3 + 1.5 * iqr;
    Some(values.iter().filter(|&&v| v < lower || v > upper).count())
}

/// Calculate median from sorted data
pub fn calculate_median(sorted_data: &[f64]) -> Option<f64> {
    if sorted_data.is_empty() {
        return None;
    }

    let len = sorted_data.len();
    #[allow(clippy::manual_is_multiple_of)]
    if len % 2 == 0 {
        // Even number of elements - average of two middle values
        let mid1 = sorted_data[len / 2 - 1];
        let mid2 = sorted_data[len / 2];
        Some(mid1.midpoint(mid2))
    } else {
        // Odd number of elements - middle value
        Some(sorted_data[len / 2])
    }
}

/// Calculate quartiles using Type 7 interpolation (R/Excel default)
pub fn calculate_quartiles(sorted_data: &[f64]) -> Option<Quartiles> {
    if sorted_data.len() < 4 {
        return None;
    }

    let n = sorted_data.len() as f64;

    let q1 = calculate_percentile(sorted_data, 0.25, n);
    let q2 = calculate_percentile(sorted_data, 0.50, n);
    let q3 = calculate_percentile(sorted_data, 0.75, n);

    let iqr = q3 - q1;

    Some(Quartiles {
        q1,
        q2,
        q3,
        iqr: iqr.is_finite().then_some(iqr),
    })
}

/// Calculate percentile using Type 7 linear interpolation (R/Excel/pandas default)
fn calculate_percentile(sorted_data: &[f64], p: f64, n: f64) -> f64 {
    // Type 7 (0-indexed): pos = p * (n - 1)
    let pos = p * (n - 1.0);

    if pos < 0.0 {
        return sorted_data[0];
    }
    if pos >= n - 1.0 {
        return sorted_data[(n as usize) - 1];
    }

    let lower_idx = pos.floor() as usize;
    let upper_idx = pos.ceil() as usize;

    if lower_idx == upper_idx {
        sorted_data[lower_idx]
    } else {
        let weight = pos - lower_idx as f64;
        sorted_data[lower_idx] * (1.0 - weight) + sorted_data[upper_idx] * weight
    }
}

/// Calculate mode (most frequent value)
/// For multimodal distributions, returns the smallest value with maximum frequency
pub fn calculate_mode(data: &[f64]) -> Option<f64> {
    if data.is_empty() {
        return None;
    }

    let mut freq_map: HashMap<u64, usize> = HashMap::new();

    // Count the actual values, without rounding distinct small numbers into
    // the same bucket. Signed zeros compare equal numerically.
    for &value in data {
        if !value.is_finite() {
            continue;
        }
        let key = if value == 0.0 { 0 } else { value.to_bits() };
        *freq_map.entry(key).or_insert(0) += 1;
    }

    // Find maximum frequency
    let max_freq = *freq_map.values().max()?;

    // If all values are unique, return None
    if max_freq == 1 {
        return None;
    }

    // Find all values with maximum frequency and return the smallest (deterministic)
    freq_map
        .iter()
        .filter(|(_, count)| **count == max_freq)
        .map(|(&bits, _)| f64::from_bits(bits))
        .min_by(f64::total_cmp)
}

/// Calculate coefficient of variation (CV)
pub fn calculate_coefficient_of_variation(std_dev: f64, mean: f64) -> Option<f64> {
    // A near-zero mean would divide by zero, and an infinite one would divide
    // any spread down to 0.0. A non-finite spread fails the result check.
    if !mean.is_finite() || mean.abs() < 0.001 {
        None
    } else {
        let value = (std_dev / mean.abs()) * 100.0;
        value.is_finite().then_some(value)
    }
}

/// Calculate skewness (Pearson's moment coefficient)
pub fn calculate_skewness(data: &[f64], mean: f64, std_dev: f64) -> Option<f64> {
    // An infinite std_dev would standardize every value to 0.0; a non-finite
    // mean fails the result check below.
    if data.len() < 3 || !std_dev.is_finite() || std_dev < 1e-10 {
        return None;
    }

    let n = data.len() as f64;
    let sum_cubed: f64 = data
        .iter()
        .map(|&x| {
            let z = (x - mean) / std_dev;
            z * z * z
        })
        .sum();

    let value = sum_cubed / n;
    value.is_finite().then_some(value)
}

/// Calculate kurtosis (excess kurtosis)
pub fn calculate_kurtosis(data: &[f64], mean: f64, std_dev: f64) -> Option<f64> {
    // Same guards as `calculate_skewness`.
    if data.len() < 4 || !std_dev.is_finite() || std_dev < 1e-10 {
        return None;
    }

    let n = data.len() as f64;
    let sum_fourth: f64 = data
        .iter()
        .map(|&x| {
            let z = (x - mean) / std_dev;
            let z2 = z * z;
            z2 * z2
        })
        .sum();

    // Excess kurtosis (subtract 3 for normal distribution)
    let value = (sum_fourth / n) - 3.0;
    value.is_finite().then_some(value)
}

/// Reservoir sampling for large datasets, seeded so the sample is the same
/// on every run for the same values.
fn reservoir_sample(data: &[f64], k: usize) -> Vec<f64> {
    use rand::rngs::SmallRng;
    use rand::{Rng, SeedableRng};

    if data.len() <= k {
        return data.to_vec();
    }

    let mut rng = SmallRng::seed_from_u64(SAMPLE_SEED);
    let mut reservoir: Vec<f64> = data[0..k].to_vec();

    for (i, &value) in data.iter().enumerate().skip(k) {
        let j = rng.random_range(0..=i);
        if j < k {
            reservoir[j] = value;
        }
    }

    reservoir
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn sampled_statistics_are_the_same_on_every_run() {
        // Above the threshold the order statistics come from a reservoir; an
        // unseeded one gave a different median on each call for the same data.
        let data: Vec<String> = (0..50_000)
            .map(|n| ((n * 7919) % 50_000).to_string())
            .collect();
        let first = compute_numeric_stats(&data);
        let second = compute_numeric_stats(&data);
        assert_eq!(first.is_approximate, Some(true));
        assert_eq!(first.median, second.median);
        assert_eq!(first.quartiles, second.quartiles);
        assert_eq!(first.mode, second.mode);
        assert_eq!(first.skewness, second.skewness);
        assert_eq!(first.kurtosis, second.kurtosis);
        assert_eq!(first.outlier_count, second.outlier_count);
    }

    #[test]
    fn test_median_odd_count() {
        let sorted = vec![1.0, 2.0, 3.0];
        assert_eq!(calculate_median(&sorted), Some(2.0));
    }

    #[test]
    fn test_median_even_count() {
        let sorted = vec![1.0, 2.0, 3.0, 4.0];
        assert_eq!(calculate_median(&sorted), Some(2.5));
    }

    #[test]
    fn test_median_empty() {
        let sorted: Vec<f64> = vec![];
        assert_eq!(calculate_median(&sorted), None);
    }

    #[test]
    fn median_preserves_finite_extremes_and_subnormals() {
        assert_eq!(calculate_median(&[f64::MAX, f64::MAX]), Some(f64::MAX));
        assert_eq!(calculate_median(&[-f64::MAX, f64::MAX]), Some(0.0));
        let smallest = f64::from_bits(1);
        assert_eq!(calculate_median(&[smallest, smallest]), Some(smallest));
    }

    #[test]
    fn test_quartiles() {
        let sorted = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let q = calculate_quartiles(&sorted).unwrap();
        assert!((q.q1 - 2.0).abs() < 0.01);
        assert_eq!(q.q2, 3.0);
        assert!((q.q3 - 4.0).abs() < 0.01);
        assert!((q.iqr.unwrap() - 2.0).abs() < 0.01);
    }

    #[test]
    fn test_quartiles_small_dataset() {
        let sorted = vec![1.0, 2.0];
        assert_eq!(calculate_quartiles(&sorted), None);
    }

    #[test]
    fn test_mode_unique_values() {
        let data = vec![1.0, 2.0, 3.0];
        assert_eq!(calculate_mode(&data), None);
    }

    #[test]
    fn test_mode_with_duplicates() {
        let data = vec![1.0, 2.0, 2.0, 3.0];
        assert_eq!(calculate_mode(&data), Some(2.0));
    }

    #[test]
    fn test_mode_empty() {
        let data: Vec<f64> = vec![];
        assert_eq!(calculate_mode(&data), None);
    }

    #[test]
    fn mode_counts_exact_values_and_breaks_ties_numerically() {
        assert_eq!(calculate_mode(&[1e-12, 2e-12, 3e-12]), None);
        assert_eq!(calculate_mode(&[1e-12, 1e-12, 2e-12]), Some(1e-12));
        assert_eq!(calculate_mode(&[-0.0, 0.0, 1.0, 1.0]), Some(0.0));
        assert_eq!(calculate_mode(&[2.0, -3.0, 2.0, -3.0]), Some(-3.0));
        assert_eq!(calculate_mode(&[f64::NAN, f64::INFINITY]), None);
    }

    #[test]
    fn test_skewness_symmetric() {
        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let mean = 3.0;
        let std_dev = (2.0_f64).sqrt();
        let skew = calculate_skewness(&data, mean, std_dev).unwrap();
        assert!(skew.abs() < 0.1); // Symmetric distribution
    }

    #[test]
    fn test_coefficient_of_variation() {
        let cv = calculate_coefficient_of_variation(10.0, 100.0);
        assert_eq!(cv, Some(10.0));
    }

    #[test]
    fn test_cv_zero_mean() {
        let cv = calculate_coefficient_of_variation(1.0, 0.0);
        assert_eq!(cv, None);
    }

    #[test]
    fn test_kurtosis() {
        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let mean = 3.0;
        let std_dev = (2.0_f64).sqrt();
        let kurt = calculate_kurtosis(&data, mean, std_dev);
        assert!(kurt.is_some());
    }

    /// The batch path is what an in-memory source (a dict, a DataFrame) is
    /// profiled with. No engine supplies exact stream aggregates to override
    /// these numbers there, so this is the whole answer (#670, #671).
    #[test]
    fn batch_stats_are_numerically_stable() {
        fn column(values: &[f64]) -> NumericStats {
            let cells: Vec<String> = values.iter().map(|value| value.to_string()).collect();
            compute_numeric_stats(&cells)
        }

        // A small spread on a large offset: `sum_squares - n * mean²` cancels
        // it away entirely, reporting a varying column as constant.
        for base in [1e6, 1e8, 1e9, 1e12] {
            let stats = column(&[base, base + 1.0, base + 2.0, base + 3.0]);
            // Equality, not a tolerance: four values accumulate to the
            // correctly rounded answer at every offset here.
            assert_eq!(stats.variance, Some(5.0 / 3.0), "offset {base}");
            assert_eq!(stats.std_dev, Some((5.0f64 / 3.0).sqrt()), "offset {base}");
            assert_eq!(stats.mean, base + 1.5, "offset {base}");
        }

        // Past the SIMD threshold, where the lanes accumulate independently.
        let many: Vec<f64> = (0..500).map(|i| 1e9 + (i % 4) as f64).collect();
        let expected_variance = (0..4)
            .map(|i| 125.0 * (i as f64 - 1.5).powi(2))
            .sum::<f64>()
            / (many.len() as f64 - 1.0);
        let stats = column(&many);
        assert!(
            (stats.variance.unwrap() - expected_variance).abs() < 1e-6,
            "500 values reported variance {:?}",
            stats.variance
        );

        // Large values that cancel: the unit contribution is the whole mean.
        for values in [[1e16, 1.0, -1e16], [1e16, -1e16, 1.0]] {
            let stats = column(&values);
            assert!(
                (stats.mean - 1.0 / 3.0).abs() < 1e-9,
                "{values:?} reported mean {}",
                stats.mean
            );
        }

        // A representable mean whose naive sum is not representable.
        assert_eq!(column(&[1e308, 1e308]).mean, 1e308);

        // And a genuinely constant column still reports no spread at all.
        let constant = column(&[1e9; 4]);
        assert_eq!(constant.variance, Some(0.0));
        assert_eq!(constant.std_dev, Some(0.0));
        assert_eq!(constant.mean, 1e9);
    }

    /// The batch path has no exact aggregates to fall back on, so a spread
    /// that overflows `f64` must be absent here, not a number (#784).
    #[test]
    fn batch_stats_report_an_overflowing_spread_as_absent() {
        let cells: Vec<String> = [1e308, -1e308, 5.0, 0.5, 5.0, 7.0]
            .iter()
            .map(|value| value.to_string())
            .collect();
        let stats = compute_numeric_stats(&cells);
        assert_eq!(stats.variance, None);
        assert_eq!(stats.std_dev, None);
        assert_eq!(stats.coefficient_of_variation, None);
        assert_eq!(stats.skewness, None);
        assert_eq!(stats.kurtosis, None);
        assert!(stats.mean.is_finite());
    }

    /// `q3 - q1` can overflow while both quartiles are representable (#802).
    #[test]
    fn quartiles_report_an_overflowing_iqr_as_absent() {
        let q = calculate_quartiles(&[-1e308, -1e308, 1e308, 1e308]).unwrap();
        assert_eq!((q.q1, q.q2, q.q3, q.iqr), (-1e308, 0.0, 1e308, None));
        assert_eq!(calculate_quartiles(&[2.0; 4]).unwrap().iqr, Some(0.0));

        // The extremes lie outside [q1, q3], but not outside fences that are
        // themselves beyond the `f64` range: no value is an outlier.
        let cells: Vec<String> = [-1.7e308, -1e308, -1e308, 1e308, 1e308, 1.7e308]
            .iter()
            .map(|value| value.to_string())
            .collect();
        let stats = compute_numeric_stats(&cells);
        let quartiles = stats.quartiles.unwrap();
        assert_eq!((quartiles.q1, quartiles.q3), (-1e308, 1e308));
        assert_eq!(quartiles.iqr, None);
        assert_eq!(stats.outlier_count, Some(0));
    }

    #[test]
    fn shape_helpers_reject_a_non_finite_spread() {
        let data = [1.0, 2.0, 3.0, 4.0];
        // An infinite std_dev would standardize every value to 0.0.
        assert_eq!(calculate_skewness(&data, 2.5, f64::INFINITY), None);
        assert_eq!(calculate_kurtosis(&data, 2.5, f64::INFINITY), None);
        assert_eq!(calculate_coefficient_of_variation(f64::NAN, 2.5), None);
        // An infinite mean would divide any spread down to 0.0.
        assert_eq!(calculate_coefficient_of_variation(1.0, f64::INFINITY), None);
        // Finite arguments whose result overflows.
        assert_eq!(calculate_skewness(&[1e200, -1e200, 0.0], 0.0, 1e-10), None);
        assert_eq!(
            calculate_kurtosis(&[1e200, -1e200, 0.0, 0.0], 0.0, 1e-10),
            None
        );
        assert_eq!(calculate_coefficient_of_variation(1e307, 0.01), None);
    }

    #[test]
    fn test_calculate_numeric_stats_basic() {
        let data = vec!["10".to_string(), "20".to_string(), "30".to_string()];
        let stats = calculate_numeric_stats(&data);

        match stats {
            ColumnStats::Numeric(n) => {
                assert_eq!(n.min, 10.0);
                assert_eq!(n.max, 30.0);
                assert_eq!(n.mean, 20.0);
                assert!(n.std_dev.unwrap() > 0.0);
                assert_eq!(n.median, Some(20.0));
            }
            _ => panic!("Expected Numeric stats"),
        }
    }

    #[test]
    fn test_calculate_numeric_stats_empty() {
        let data: Vec<String> = vec![];
        let stats = calculate_numeric_stats(&data);

        match stats {
            ColumnStats::Numeric(n) => {
                assert_eq!(n.min, 0.0);
                assert_eq!(n.max, 0.0);
                assert_eq!(n.mean, 0.0);
                assert_eq!(n.median, None);
            }
            _ => panic!("Expected Numeric stats"),
        }
    }

    #[test]
    fn test_numeric_stats_filters_nan() {
        let data = vec!["1.0".to_string(), "NaN".to_string(), "3.0".to_string()];
        let stats = calculate_numeric_stats(&data);

        match stats {
            ColumnStats::Numeric(n) => {
                assert_eq!(n.min, 1.0);
                assert_eq!(n.max, 3.0);
                assert_eq!(n.mean, 2.0);
            }
            _ => panic!("Expected Numeric stats"),
        }
    }

    #[test]
    fn test_numeric_stats_filters_infinity() {
        let data = vec![
            "1.0".to_string(),
            "inf".to_string(),
            "3.0".to_string(),
            "-inf".to_string(),
        ];
        let stats = calculate_numeric_stats(&data);

        match stats {
            ColumnStats::Numeric(n) => {
                assert_eq!(n.min, 1.0);
                assert_eq!(n.max, 3.0);
                assert_eq!(n.mean, 2.0);
            }
            _ => panic!("Expected Numeric stats"),
        }
    }

    #[test]
    fn test_numeric_stats_all_nan() {
        let data = vec!["NaN".to_string(), "inf".to_string(), "-inf".to_string()];
        let stats = calculate_numeric_stats(&data);

        match stats {
            ColumnStats::Numeric(n) => {
                assert_eq!(n.min, 0.0);
                assert_eq!(n.max, 0.0);
                assert_eq!(n.mean, 0.0);
                assert_eq!(n.median, None);
            }
            _ => panic!("Expected Numeric stats"),
        }
    }
}