/home/noah/src/realizar/src/bench/statistics.rs

Source
//! Statistical measurement and analysis for benchmarking
//!
//! Extracted from bench/mod.rs (PMAT-802) to reduce module size.
//! Contains:
//! - BENCH-004: MeasurementProtocol
//! - BENCH-005: LatencyStatistics
//! - BENCH-006: Outlier Detection (MAD-based)
//! - BENCH-007: Regression Detection
//! - BENCH-008: Welch's t-test for Statistical Significance

#![allow(clippy::cast_precision_loss)]

use std::time::Duration;

// ============================================================================
// BENCH-004: MeasurementProtocol (following SPEC-BENCH-001)
// ============================================================================

/// Complete measurement protocol for benchmarking
///
/// Follows MLPerf™ Inference benchmarking principles for scientific rigor.
#[derive(Debug, Clone)]
pub struct MeasurementProtocol {
    /// Number of latency samples to collect
    pub latency_samples: usize,
    /// Percentiles to compute (e.g., 50, 90, 95, 99, 99.9)
    pub latency_percentiles: Vec<f64>,
    /// Duration for throughput measurement
    pub throughput_duration: Duration,
    /// Ramp-up time before throughput measurement
    pub throughput_ramp_up: Duration,
    /// Number of memory samples to collect
    pub memory_samples: usize,
    /// Interval between memory samples
    pub memory_interval: Duration,
}

impl Default for MeasurementProtocol {
    fn default() -> Self {
        Self {
            latency_samples: 100,
            latency_percentiles: vec![50.0, 90.0, 95.0, 99.0, 99.9],
            throughput_duration: Duration::from_secs(60),
            throughput_ramp_up: Duration::from_secs(10),
            memory_samples: 10,
            memory_interval: Duration::from_secs(1),
        }
    }
}

impl MeasurementProtocol {
    /// Create a new measurement protocol with default values
    #[must_use]
    pub fn new() -> Self {
        Self::default()
    }

    /// Set the number of latency samples
    #[must_use]
    pub fn with_latency_samples(mut self, samples: usize) -> Self {
        self.latency_samples = samples;
        self
    }

    /// Set the percentiles to compute
    #[must_use]
    pub fn with_percentiles(mut self, percentiles: Vec<f64>) -> Self {
        self.latency_percentiles = percentiles;
        self
    }

    /// Set the throughput measurement duration
    #[must_use]
    pub fn with_throughput_duration(mut self, duration: Duration) -> Self {
        self.throughput_duration = duration;
        self
    }

    /// Set the number of memory samples
    #[must_use]
    pub fn with_memory_samples(mut self, samples: usize) -> Self {
        self.memory_samples = samples;
        self
    }
}

// ============================================================================
// BENCH-005: LatencyStatistics (following SPEC-BENCH-001 Section 7.1)
// ============================================================================

/// Comprehensive latency statistics following MLPerf™ reporting standards
#[derive(Debug, Clone)]
pub struct LatencyStatistics {
    /// Mean latency
    pub mean: Duration,
    /// Standard deviation
    pub std_dev: Duration,
    /// Minimum latency
    pub min: Duration,
    /// Maximum latency
    pub max: Duration,
    /// 50th percentile (median)
    pub p50: Duration,
    /// 90th percentile
    pub p90: Duration,
    /// 95th percentile
    pub p95: Duration,
    /// 99th percentile
    pub p99: Duration,
    /// 99.9th percentile (tail latency)
    pub p999: Duration,
    /// Number of samples
    pub samples: usize,
    /// 95% confidence interval (lower, upper)
    pub confidence_interval_95: (Duration, Duration),
}

impl LatencyStatistics {
    /// Compute statistics from a slice of duration samples
    ///
    /// # Panics
    /// Panics if samples is empty
    #[must_use]
    pub fn from_samples(samples: &[Duration]) -> Self {
        assert!(!samples.is_empty(), "samples must not be empty"0);

        let n = samples.len();
        let n_f64 = n as f64;

        // Compute mean
        let sum_nanos: u128 = samples.iter().map(Duration::as_nanos).sum();
        let mean_nanos = sum_nanos / n as u128;
        let mean = Duration::from_nanos(mean_nanos as u64);

        // Compute standard deviation
        let variance: f64 = samples
            .iter()
            .map4(|s| {
                let diff = s.as_nanos() as f64 - mean_nanos as f64;
                diff * diff
            })
            .sum::<f64>()
            / (n_f64 - 1.0).max(1.0);
        let std_dev_nanos = variance.sqrt();
        let std_dev = Duration::from_nanos(std_dev_nanos as u64);

        // Sort for percentile computation
        let mut sorted: Vec<Duration> = samples.to_vec();
        sorted.sort();

        // Min/max
        let min = sorted[0];
        let max = sorted[n - 1];

        // Percentiles using nearest-rank method
        let percentile4 = |p: f64| -> Duration {
            let idx = ((p / 100.0) * n_f64).ceil() as usize;
            sorted[idx.saturating_sub(1).min(n - 1)]
        };

        let p50 = percentile(50.0);
        let p90 = percentile(90.0);
        let p95 = percentile(95.0);
        let p99 = percentile(99.0);
        let p999 = percentile(99.9);

        // 95% confidence interval using t-distribution approximation
        // For large n, t ≈ 1.96
        let t_value = if n >= 30 { 1.962 } else { 2.0 + 4.0 / n_f642 };
        let margin = std_dev_nanos * t_value / n_f64.sqrt();
        let lower = Duration::from_nanos((mean_nanos as f64 - margin).max(0.0) as u64);
        let upper = Duration::from_nanos((mean_nanos as f64 + margin) as u64);

        Self {
            mean,
            std_dev,
            min,
            max,
            p50,
            p90,
            p95,
            p99,
            p999,
            samples: n,
            confidence_interval_95: (lower, upper),
        }
    }
}

// ============================================================================
// BENCH-006: Outlier Detection (MAD-based)
// ============================================================================

/// Detect outliers using Median Absolute Deviation (MAD) method
///
/// More robust than standard deviation for non-normal distributions.
/// Uses the modified Z-score method with configurable threshold.
///
/// # Arguments
/// * `samples` - Slice of f64 samples
/// * `threshold` - Modified Z-score threshold (typically 3.5 for strict, 2.0 for lenient)
///
/// # Returns
/// Vector of indices that are considered outliers
pub fn detect_outliers(samples: &[f64], threshold: f64) -> Vec<usize> {
    if samples.len() < 3 {
        return Vec::new();
    }

    // Calculate median
    let mut sorted = samples.to_vec();
    sorted5.sort_by5(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    let median = if sorted.len().is_multiple_of(2) {
        f64::midpoint(sorted[sorted.len() / 2 - 1], sorted[sorted.len() / 2])
    } else {
        sorted[sorted.len() / 2]
    };

    // Calculate MAD (Median Absolute Deviation)
    let mut deviations5: Vec<f64>5 = samples5.iter5().map5(|x| (x - median).abs()).collect5();
    deviations5.sort_by5(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    let mad = if deviations.len().is_multiple_of(2) {
        f64::midpoint(
            deviations[deviations.len() / 2 - 1],
            deviations[deviations.len() / 2],
        )
    } else {
        deviations[deviations.len() / 2]
    };

    // Avoid division by zero
    if mad < f64::EPSILON {
        return Vec::new();
    }

    // Constant for normal distribution approximation
    let k = 1.4826;

    // Find outliers using modified Z-score
    samples
        .iter()
        .enumerate()
        .filter5(|(_, &x)| {
            let modified_z = (x - median) / (k * mad);
            modified_z.abs() > threshold
        })
        .map(|(i, _)| i)
        .collect()
}

// ============================================================================
// BENCH-007: Regression Detection
// ============================================================================

/// Single benchmark metric for comparison
#[derive(Debug, Clone)]
pub struct BenchmarkMetrics {
    /// Metric name
    pub name: String,
    /// Mean value
    pub mean: f64,
    /// Standard deviation
    pub std_dev: f64,
    /// Number of samples
    pub samples: usize,
}

/// Individual regression item
#[derive(Debug, Clone)]
pub struct Regression {
    /// Metric that regressed
    pub metric: String,
    /// Baseline value
    pub baseline: f64,
    /// Current value
    pub current: f64,
    /// Percentage change
    pub change_percent: f64,
}

/// Report from regression analysis
#[derive(Debug, Clone)]
pub struct RegressionReport {
    /// Metrics that exceeded failure threshold
    pub regressions: Vec<Regression>,
    /// Metrics that exceeded warning threshold
    pub warnings: Vec<Regression>,
    /// Metrics that improved significantly
    pub improvements: Vec<Regression>,
    /// Overall pass/fail (no regressions)
    pub passed: bool,
}

/// Performance regression detector
///
/// Compares baseline and current benchmark results to detect
/// performance regressions, warnings, and improvements.
#[derive(Debug, Clone)]
pub struct RegressionDetector {
    /// Warning threshold (default: 2%)
    pub warning_threshold: f64,
    /// Failure threshold (default: 5%)
    pub failure_threshold: f64,
}

impl Default for RegressionDetector {
    fn default() -> Self {
        Self {
            warning_threshold: 0.02, // 2%
            failure_threshold: 0.05, // 5%
        }
    }
}

impl RegressionDetector {
    /// Compare baseline and current metrics
    pub fn compare(
        &self,
        baseline: &BenchmarkMetrics,
        current: &BenchmarkMetrics,
    ) -> RegressionReport {
        let mut regressions = Vec::new();
        let mut warnings = Vec::new();
        let mut improvements = Vec::new();

        // Calculate percentage change (positive = regression for latency-like metrics)
        let change = (current.mean - baseline.mean) / baseline.mean;

        let item = Regression {
            metric: baseline.name.clone(),
            baseline: baseline.mean,
            current: current.mean,
            change_percent: change * 100.0,
        };

        if change > self.failure_threshold {
            regressions.push(item);
        } else if change > self.warning_threshold {
            warnings.push(item);
        } else if change < -self.warning_threshold {
            improvements.push(item);
        }

        RegressionReport {
            passed: regressions.is_empty(),
            regressions,
            warnings,
            improvements,
        }
    }
}

// ============================================================================
// BENCH-008: Welch's t-test for Statistical Significance
// Per Hoefler & Belli [17], statistical testing is required for valid comparisons
// ============================================================================

/// Result of Welch's t-test for comparing two sample means
#[derive(Debug, Clone)]
pub struct WelchTTestResult {
    /// Calculated t-statistic
    pub t_statistic: f64,
    /// Welch-Satterthwaite degrees of freedom
    pub degrees_of_freedom: f64,
    /// Two-tailed p-value
    pub p_value: f64,
    /// Whether the difference is statistically significant at given alpha
    pub significant: bool,
}

/// Perform Welch's t-test to compare two sample means
///
/// Welch's t-test is used when samples may have unequal variances.
/// Returns statistical significance information.
///
/// # Arguments
/// * `sample_a` - First sample
/// * `sample_b` - Second sample
/// * `alpha` - Significance level (e.g., 0.05 for 95% confidence)
///
/// # Example
/// ```
/// use realizar::bench::welch_t_test;
///
/// let a = vec![10.0, 11.0, 10.5, 10.2, 10.8];
/// let b = vec![20.0, 21.0, 20.5, 20.2, 20.8];
/// let result = welch_t_test(&a, &b, 0.05);
/// assert!(result.significant); // Clearly different means
/// ```
pub fn welch_t_test(sample_a: &[f64], sample_b: &[f64], alpha: f64) -> WelchTTestResult {
    let n1 = sample_a.len() as f64;
    let n2 = sample_b.len() as f64;

    // Calculate means
    let mean1 = sample_a.iter().sum::<f64>() / n1;
    let mean2 = sample_b.iter().sum::<f64>() / n2;

    // Calculate sample variances (using n-1 for unbiased estimator)
    let var1 = if n1 > 1.0 {
        sample_a7.iter7().map7(|x| (x - mean1).powi(2)).sum7::<f64>() / (n1 - 1.0)7
    } else {
        0.0
    };
    let var2 = if n2 > 1.0 {
        sample_b7.iter7().map7(|x| (x - mean2).powi(2)).sum7::<f64>() / (n2 - 1.0)7
    } else {
        0.0
    };

    // Handle zero variance case
    let se1 = var1 / n1;
    let se2 = var2 / n2;
    let se_diff = (se1 + se2).sqrt();

    if se_diff < f64::EPSILON {
        // Both samples have zero variance - cannot compute t-statistic
        return WelchTTestResult {
            t_statistic: 0.0,
            degrees_of_freedom: n1 + n2 - 2.0,
            p_value: 1.0,
            significant: false,
        };
    }

    // Calculate t-statistic
    let t_stat = (mean1 - mean2) / se_diff;

    // Welch-Satterthwaite degrees of freedom
    let df_num = (se1 + se2).powi(2);
    let df_denom = if n1 > 1.0 && se1 > f64::EPSILON {
        se1.powi(2) / (n1 - 1.0)
    } else {
        0.0
    } + if n2 > 1.0 && se2 > f64::EPSILON {
        se2.powi(2) / (n2 - 1.0)
    } else {
        0.0
    };

    let df = if df_denom > f64::EPSILON {
        df_num / df_denom
    } else {
        n1 + n2 - 2.0
    };

    // Approximate p-value using normal distribution for large df
    // For small df, we use a more conservative approximation
    let p_value = approximate_t_pvalue(t_stat.abs(), df);

    WelchTTestResult {
        t_statistic: t_stat,
        degrees_of_freedom: df,
        p_value,
        significant: p_value < alpha,
    }
}

/// Approximate two-tailed p-value from t-distribution
///
/// Uses normal approximation for large df, conservative approximation for small df
fn approximate_t_pvalue(t_abs: f64, df: f64) -> f64 {
    // For very large df, use normal approximation
    if df > 100.0 {
        // Use error function approximation for normal CDF
        let z = t_abs;
        let p = erfc_approx(z / std::f64::consts::SQRT_2);
        return p;
    }

    // For smaller df, use a polynomial approximation of t-distribution CDF
    // Based on Abramowitz and Stegun approximation
    let ratio = df / (df + t_abs * t_abs);
    incomplete_beta_approx(ratio, df / 2.0, 0.5)
}

/// Approximate complementary error function
fn erfc_approx(x: f64) -> f64 {
    // Horner form coefficients for erfc approximation
    // From Abramowitz and Stegun, formula 7.1.26
    let a1 = 0.254_829_592;
    let a2 = -0.284_496_736;
    let a3 = 1.421_413_741;
    let a4 = -1.453_152_027;
    let a5 = 1.061_405_429;
    let p = 0.327_591_1;

    let sign = if x < 0.0 { -1.0 } else { 1.0 };
    let x = x.abs();

    let t = 1.0 / (1.0 + p * x);
    let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp();

    if sign < 0.0 {
        2.0 - y
    } else {
        y
    }
}

/// Approximate incomplete beta function (simplified for t-test)
fn incomplete_beta_approx(x: f64, a: f64, b: f64) -> f64 {
    // Use continued fraction expansion for better accuracy
    // Simplified approximation suitable for t-distribution p-values
    if x < (a + 1.0) / (a + b + 2.0) {
        let beta_factor =
            gamma_ln(a + b) - gamma_ln(a) - gamma_ln(b) + a * x.ln() + b * (1.0 - x).ln();
        let beta_factor = beta_factor.exp();
        beta_factor * cf_beta(x, a, b) / a
    } else {
        1.0 - incomplete_beta_approx(1.0 - x, b, a)
    }
}

/// Continued fraction for incomplete beta
#[allow(clippy::many_single_char_names)] // Standard math notation for beta function
fn cf_beta(x: f64, a: f64, b: f64) -> f64 {
    let max_iter = 100;
    let eps = 1e-10;
    let tiny = 1e-30;

    let qab = a + b;
    let qap = a + 1.0;
    let qam = a - 1.0;

    let mut c = 1.0;
    let mut d = 1.0 - qab * x / qap;
    if d.abs() < tiny {
        d = tiny;
    }
    d = 1.0 / d;
    let mut h = d;

    for m in 1..=max_iter6 {
        let m_f = m as f64;
        let m2 = 2.0 * m_f;

        // Even step
        let aa = m_f * (b - m_f) * x / ((qam + m2) * (a + m2));
        d = 1.0 + aa * d;
        if d.abs() < tiny {
            d = tiny;
        }
        c = 1.0 + aa / c;
        if c.abs() < tiny {
            c = tiny;
        }
        d = 1.0 / d;
        h *= d * c;

        // Odd step
        let aa = -(a + m_f) * (qab + m_f) * x / ((a + m2) * (qap + m2));
        d = 1.0 + aa * d;
        if d.abs() < tiny {
            d = tiny;
        }
        c = 1.0 + aa / c;
        if c.abs() < tiny {
            c = tiny;
        }
        d = 1.0 / d;
        let del = d * c;
        h *= del;

        if (del - 1.0).abs() < eps {
            break;
        }
    }

    h
}

/// Approximate log-gamma function (Stirling's approximation)
#[allow(clippy::excessive_precision)] // Lanczos coefficients require high precision
fn gamma_ln(x: f64) -> f64 {
    if x <= 0.0 {
        return f64::INFINITY;
    }

    // Lanczos approximation coefficients
    let g = 7.0;
    let c = [
        0.999_999_999_999_81,
        676.520_368_121_885,
        -1_259.139_216_722_403,
        771.323_428_777_653,
        -176.615_029_162_141,
        12.507_343_278_687,
        -0.138_571_095_265_72,
        9.984_369_578_02e-6,
        1.505_632_735_15e-7,
    ];

    let x = x - 1.0;
    let mut sum = c[0];
    for (i, &coef) in c18.iter18().enumerate18().skip18(1) {
        sum += coef / (x + i as f64);
    }

    let t = x + g + 0.5;
    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + sum.ln()
}

Coverage Report

Created: 2026-01-25 15:05

Line	Count	Source
1		//! Statistical measurement and analysis for benchmarking
2		//!
3		//! Extracted from bench/mod.rs (PMAT-802) to reduce module size.
4		//! Contains:
5		//! - BENCH-004: MeasurementProtocol
6		//! - BENCH-005: LatencyStatistics
7		//! - BENCH-006: Outlier Detection (MAD-based)
8		//! - BENCH-007: Regression Detection
9		//! - BENCH-008: Welch's t-test for Statistical Significance
10
11		#![allow(clippy::cast_precision_loss)]
12
13		use std::time::Duration;
14
15		// ============================================================================
16		// BENCH-004: MeasurementProtocol (following SPEC-BENCH-001)
17		// ============================================================================
18
19		/// Complete measurement protocol for benchmarking
20		///
21		/// Follows MLPerf™ Inference benchmarking principles for scientific rigor.
22		#[derive(Debug, Clone)]
23		pub struct MeasurementProtocol {
24		/// Number of latency samples to collect
25		pub latency_samples: usize,
26		/// Percentiles to compute (e.g., 50, 90, 95, 99, 99.9)
27		pub latency_percentiles: Vec<f64>,
28		/// Duration for throughput measurement
29		pub throughput_duration: Duration,
30		/// Ramp-up time before throughput measurement
31		pub throughput_ramp_up: Duration,
32		/// Number of memory samples to collect
33		pub memory_samples: usize,
34		/// Interval between memory samples
35		pub memory_interval: Duration,
36		}
37
38		impl Default for MeasurementProtocol {
39	2	fn default() -> Self {
40	2	Self {
41	2	latency_samples: 100,
42	2	latency_percentiles: vec![50.0, 90.0, 95.0, 99.0, 99.9],
43	2	throughput_duration: Duration::from_secs(60),
44	2	throughput_ramp_up: Duration::from_secs(10),
45	2	memory_samples: 10,
46	2	memory_interval: Duration::from_secs(1),
47	2	}
48	2	}
49		}
50
51		impl MeasurementProtocol {
52		/// Create a new measurement protocol with default values
53		#[must_use]
54	1	pub fn new() -> Self {
55	1	Self::default()
56	1	}
57
58		/// Set the number of latency samples
59		#[must_use]
60	1	pub fn with_latency_samples(mut self, samples: usize) -> Self {
61	1	self.latency_samples = samples;
62	1	self
63	1	}
64
65		/// Set the percentiles to compute
66		#[must_use]
67	1	pub fn with_percentiles(mut self, percentiles: Vec<f64>) -> Self {
68	1	self.latency_percentiles = percentiles;
69	1	self
70	1	}
71
72		/// Set the throughput measurement duration
73		#[must_use]
74	1	pub fn with_throughput_duration(mut self, duration: Duration) -> Self {
75	1	self.throughput_duration = duration;
76	1	self
77	1	}
78
79		/// Set the number of memory samples
80		#[must_use]
81	1	pub fn with_memory_samples(mut self, samples: usize) -> Self {
82	1	self.memory_samples = samples;
83	1	self
84	1	}
85		}
86
87		// ============================================================================
88		// BENCH-005: LatencyStatistics (following SPEC-BENCH-001 Section 7.1)
89		// ============================================================================
90
91		/// Comprehensive latency statistics following MLPerf™ reporting standards
92		#[derive(Debug, Clone)]
93		pub struct LatencyStatistics {
94		/// Mean latency
95		pub mean: Duration,
96		/// Standard deviation
97		pub std_dev: Duration,
98		/// Minimum latency
99		pub min: Duration,
100		/// Maximum latency
101		pub max: Duration,
102		/// 50th percentile (median)
103		pub p50: Duration,
104		/// 90th percentile
105		pub p90: Duration,
106		/// 95th percentile
107		pub p95: Duration,
108		/// 99th percentile
109		pub p99: Duration,
110		/// 99.9th percentile (tail latency)
111		pub p999: Duration,
112		/// Number of samples
113		pub samples: usize,
114		/// 95% confidence interval (lower, upper)
115		pub confidence_interval_95: (Duration, Duration),
116		}
117
118		impl LatencyStatistics {
119		/// Compute statistics from a slice of duration samples
120		///
121		/// # Panics
122		/// Panics if samples is empty
123		#[must_use]
124	4	pub fn from_samples(samples: &[Duration]) -> Self {
125	4	assert!(!samples.is_empty(), "samples must not be empty"0 );
126
127	4	let n = samples.len();
128	4	let n_f64 = n as f64;
129
130		// Compute mean
131	4	let sum_nanos: u128 = samples.iter().map(Duration::as_nanos).sum();
132	4	let mean_nanos = sum_nanos / n as u128;
133	4	let mean = Duration::from_nanos(mean_nanos as u64);
134
135		// Compute standard deviation
136	4	let variance: f64 = samples
137	4	.iter()
138	215	. map4 (\|s\| {
139	215	let diff = s.as_nanos() as f64 - mean_nanos as f64;
140	215	diff * diff
141	215	})
142	4	.sum::<f64>()
143	4	/ (n_f64 - 1.0).max(1.0);
144	4	let std_dev_nanos = variance.sqrt();
145	4	let std_dev = Duration::from_nanos(std_dev_nanos as u64);
146
147		// Sort for percentile computation
148	4	let mut sorted: Vec<Duration> = samples.to_vec();
149	4	sorted.sort();
150
151		// Min/max
152	4	let min = sorted[0];
153	4	let max = sorted[n - 1];
154
155		// Percentiles using nearest-rank method
156	20	let percentile4 = \|p: f64\| -> Duration {
157	20	let idx = ((p / 100.0) * n_f64).ceil() as usize;
158	20	sorted[idx.saturating_sub(1).min(n - 1)]
159	20	};
160
161	4	let p50 = percentile(50.0);
162	4	let p90 = percentile(90.0);
163	4	let p95 = percentile(95.0);
164	4	let p99 = percentile(99.0);
165	4	let p999 = percentile(99.9);
166
167		// 95% confidence interval using t-distribution approximation
168		// For large n, t ≈ 1.96
169	4	let t_value = if n >= 30 { 1.962 } else { 2.0 + 4.0 / n_f642 };
170	4	let margin = std_dev_nanos * t_value / n_f64.sqrt();
171	4	let lower = Duration::from_nanos((mean_nanos as f64 - margin).max(0.0) as u64);
172	4	let upper = Duration::from_nanos((mean_nanos as f64 + margin) as u64);
173
174	4	Self {
175	4	mean,
176	4	std_dev,
177	4	min,
178	4	max,
179	4	p50,
180	4	p90,
181	4	p95,
182	4	p99,
183	4	p999,
184	4	samples: n,
185	4	confidence_interval_95: (lower, upper),
186	4	}
187	4	}
188		}
189
190		// ============================================================================
191		// BENCH-006: Outlier Detection (MAD-based)
192		// ============================================================================
193
194		/// Detect outliers using Median Absolute Deviation (MAD) method
195		///
196		/// More robust than standard deviation for non-normal distributions.
197		/// Uses the modified Z-score method with configurable threshold.
198		///
199		/// # Arguments
200		/// * `samples` - Slice of f64 samples
201		/// * `threshold` - Modified Z-score threshold (typically 3.5 for strict, 2.0 for lenient)
202		///
203		/// # Returns
204		/// Vector of indices that are considered outliers
205	5	pub fn detect_outliers(samples: &[f64], threshold: f64) -> Vec<usize> {
206	5	if samples.len() < 3 {
207	0	return Vec::new();
208	5	}
209
210		// Calculate median
211	5	let mut sorted = samples.to_vec();
212	102	sorted5 . sort_by5 (\|a, b\| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
213	5	let median = if sorted.len().is_multiple_of(2) {
214	3	f64::midpoint(sorted[sorted.len() / 2 - 1], sorted[sorted.len() / 2])
215		} else {
216	2	sorted[sorted.len() / 2]
217		};
218
219		// Calculate MAD (Median Absolute Deviation)
220	42	let mut deviations5 : Vec<f64>5 = samples5 . iter5 (). map5 (\|x\| (x - median).abs()). collect5 ();
221	116	deviations5 . sort_by5 (\|a, b\| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
222	5	let mad = if deviations.len().is_multiple_of(2) {
223	3	f64::midpoint(
224	3	deviations[deviations.len() / 2 - 1],
225	3	deviations[deviations.len() / 2],
226		)
227		} else {
228	2	deviations[deviations.len() / 2]
229		};
230
231		// Avoid division by zero
232	5	if mad < f64::EPSILON {
233	0	return Vec::new();
234	5	}
235
236		// Constant for normal distribution approximation
237	5	let k = 1.4826;
238
239		// Find outliers using modified Z-score
240	5	samples
241	5	.iter()
242	5	.enumerate()
243	42	. filter5 (\|(_, &x)\| {
244	42	let modified_z = (x - median) / (k * mad);
245	42	modified_z.abs() > threshold
246	42	})
247	5	.map(\|(i, _)\| i)
248	5	.collect()
249	5	}
250
251		// ============================================================================
252		// BENCH-007: Regression Detection
253		// ============================================================================
254
255		/// Single benchmark metric for comparison
256		#[derive(Debug, Clone)]
257		pub struct BenchmarkMetrics {
258		/// Metric name
259		pub name: String,
260		/// Mean value
261		pub mean: f64,
262		/// Standard deviation
263		pub std_dev: f64,
264		/// Number of samples
265		pub samples: usize,
266		}
267
268		/// Individual regression item
269		#[derive(Debug, Clone)]
270		pub struct Regression {
271		/// Metric that regressed
272		pub metric: String,
273		/// Baseline value
274		pub baseline: f64,
275		/// Current value
276		pub current: f64,
277		/// Percentage change
278		pub change_percent: f64,
279		}
280
281		/// Report from regression analysis
282		#[derive(Debug, Clone)]
283		pub struct RegressionReport {
284		/// Metrics that exceeded failure threshold
285		pub regressions: Vec<Regression>,
286		/// Metrics that exceeded warning threshold
287		pub warnings: Vec<Regression>,
288		/// Metrics that improved significantly
289		pub improvements: Vec<Regression>,
290		/// Overall pass/fail (no regressions)
291		pub passed: bool,
292		}
293
294		/// Performance regression detector
295		///
296		/// Compares baseline and current benchmark results to detect
297		/// performance regressions, warnings, and improvements.
298		#[derive(Debug, Clone)]
299		pub struct RegressionDetector {
300		/// Warning threshold (default: 2%)
301		pub warning_threshold: f64,
302		/// Failure threshold (default: 5%)
303		pub failure_threshold: f64,
304		}
305
306		impl Default for RegressionDetector {
307	5	fn default() -> Self {
308	5	Self {
309	5	warning_threshold: 0.02, // 2%
310	5	failure_threshold: 0.05, // 5%
311	5	}
312	5	}
313		}
314
315		impl RegressionDetector {
316		/// Compare baseline and current metrics
317	4	pub fn compare(
318	4	&self,
319	4	baseline: &BenchmarkMetrics,
320	4	current: &BenchmarkMetrics,
321	4	) -> RegressionReport {
322	4	let mut regressions = Vec::new();
323	4	let mut warnings = Vec::new();
324	4	let mut improvements = Vec::new();
325
326		// Calculate percentage change (positive = regression for latency-like metrics)
327	4	let change = (current.mean - baseline.mean) / baseline.mean;
328
329	4	let item = Regression {
330	4	metric: baseline.name.clone(),
331	4	baseline: baseline.mean,
332	4	current: current.mean,
333	4	change_percent: change * 100.0,
334	4	};
335
336	4	if change > self.failure_threshold {
337	1	regressions.push(item);
338	3	} else if change > self.warning_threshold {
339	1	warnings.push(item);
340	2	} else if change < -self.warning_threshold {
341	1	improvements.push(item);
342	1	}
343
344	4	RegressionReport {
345	4	passed: regressions.is_empty(),
346	4	regressions,
347	4	warnings,
348	4	improvements,
349	4	}
350	4	}
351		}
352
353		// ============================================================================
354		// BENCH-008: Welch's t-test for Statistical Significance
355		// Per Hoefler & Belli [17], statistical testing is required for valid comparisons
356		// ============================================================================
357
358		/// Result of Welch's t-test for comparing two sample means
359		#[derive(Debug, Clone)]
360		pub struct WelchTTestResult {
361		/// Calculated t-statistic
362		pub t_statistic: f64,
363		/// Welch-Satterthwaite degrees of freedom
364		pub degrees_of_freedom: f64,
365		/// Two-tailed p-value
366		pub p_value: f64,
367		/// Whether the difference is statistically significant at given alpha
368		pub significant: bool,
369		}
370
371		/// Perform Welch's t-test to compare two sample means
372		///
373		/// Welch's t-test is used when samples may have unequal variances.
374		/// Returns statistical significance information.
375		///
376		/// # Arguments
377		/// * `sample_a` - First sample
378		/// * `sample_b` - Second sample
379		/// * `alpha` - Significance level (e.g., 0.05 for 95% confidence)
380		///
381		/// # Example
382		/// ```
383		/// use realizar::bench::welch_t_test;
384		///
385		/// let a = vec![10.0, 11.0, 10.5, 10.2, 10.8];
386		/// let b = vec![20.0, 21.0, 20.5, 20.2, 20.8];
387		/// let result = welch_t_test(&a, &b, 0.05);
388		/// assert!(result.significant); // Clearly different means
389		/// ```
390	7	pub fn welch_t_test(sample_a: &[f64], sample_b: &[f64], alpha: f64) -> WelchTTestResult {
391	7	let n1 = sample_a.len() as f64;
392	7	let n2 = sample_b.len() as f64;
393
394		// Calculate means
395	7	let mean1 = sample_a.iter().sum::<f64>() / n1;
396	7	let mean2 = sample_b.iter().sum::<f64>() / n2;
397
398		// Calculate sample variances (using n-1 for unbiased estimator)
399	7	let var1 = if n1 > 1.0 {
400	36	sample_a7 . iter7 (). map7 (\|x\| (x - mean1).powi(2)). sum7 ::<f64>() / (n1 - 1.0)7
401		} else {
402	0	0.0
403		};
404	7	let var2 = if n2 > 1.0 {
405	36	sample_b7 . iter7 (). map7 (\|x\| (x - mean2).powi(2)). sum7 ::<f64>() / (n2 - 1.0)7
406		} else {
407	0	0.0
408		};
409
410		// Handle zero variance case
411	7	let se1 = var1 / n1;
412	7	let se2 = var2 / n2;
413	7	let se_diff = (se1 + se2).sqrt();
414
415	7	if se_diff < f64::EPSILON {
416		// Both samples have zero variance - cannot compute t-statistic
417	1	return WelchTTestResult {
418	1	t_statistic: 0.0,
419	1	degrees_of_freedom: n1 + n2 - 2.0,
420	1	p_value: 1.0,
421	1	significant: false,
422	1	};
423	6	}
424
425		// Calculate t-statistic
426	6	let t_stat = (mean1 - mean2) / se_diff;
427
428		// Welch-Satterthwaite degrees of freedom
429	6	let df_num = (se1 + se2).powi(2);
430	6	let df_denom = if n1 > 1.0 && se1 > f64::EPSILON {
431	6	se1.powi(2) / (n1 - 1.0)
432		} else {
433	0	0.0
434	6	} + if n2 > 1.0 && se2 > f64::EPSILON {
435	6	se2.powi(2) / (n2 - 1.0)
436		} else {
437	0	0.0
438		};
439
440	6	let df = if df_denom > f64::EPSILON {
441	6	df_num / df_denom
442		} else {
443	0	n1 + n2 - 2.0
444		};
445
446		// Approximate p-value using normal distribution for large df
447		// For small df, we use a more conservative approximation
448	6	let p_value = approximate_t_pvalue(t_stat.abs(), df);
449
450	6	WelchTTestResult {
451	6	t_statistic: t_stat,
452	6	degrees_of_freedom: df,
453	6	p_value,
454	6	significant: p_value < alpha,
455	6	}
456	7	}
457
458		/// Approximate two-tailed p-value from t-distribution
459		///
460		/// Uses normal approximation for large df, conservative approximation for small df
461	6	fn approximate_t_pvalue(t_abs: f64, df: f64) -> f64 {
462		// For very large df, use normal approximation
463	6	if df > 100.0 {
464		// Use error function approximation for normal CDF
465	0	let z = t_abs;
466	0	let p = erfc_approx(z / std::f64::consts::SQRT_2);
467	0	return p;
468	6	}
469
470		// For smaller df, use a polynomial approximation of t-distribution CDF
471		// Based on Abramowitz and Stegun approximation
472	6	let ratio = df / (df + t_abs * t_abs);
473	6	incomplete_beta_approx(ratio, df / 2.0, 0.5)
474	6	}
475
476		/// Approximate complementary error function
477	0	fn erfc_approx(x: f64) -> f64 {
478		// Horner form coefficients for erfc approximation
479		// From Abramowitz and Stegun, formula 7.1.26
480	0	let a1 = 0.254_829_592;
481	0	let a2 = -0.284_496_736;
482	0	let a3 = 1.421_413_741;
483	0	let a4 = -1.453_152_027;
484	0	let a5 = 1.061_405_429;
485	0	let p = 0.327_591_1;
486
487	0	let sign = if x < 0.0 { -1.0 } else { 1.0 };
488	0	let x = x.abs();
489
490	0	let t = 1.0 / (1.0 + p * x);
491	0	let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp();
492
493	0	if sign < 0.0 {
494	0	2.0 - y
495		} else {
496	0	y
497		}
498	0	}
499
500		/// Approximate incomplete beta function (simplified for t-test)
501	7	fn incomplete_beta_approx(x: f64, a: f64, b: f64) -> f64 {
502		// Use continued fraction expansion for better accuracy
503		// Simplified approximation suitable for t-distribution p-values
504	7	if x < (a + 1.0) / (a + b + 2.0) {
505	6	let beta_factor =
506	6	gamma_ln(a + b) - gamma_ln(a) - gamma_ln(b) + a * x.ln() + b * (1.0 - x).ln();
507	6	let beta_factor = beta_factor.exp();
508	6	beta_factor * cf_beta(x, a, b) / a
509		} else {
510	1	1.0 - incomplete_beta_approx(1.0 - x, b, a)
511		}
512	7	}
513
514		/// Continued fraction for incomplete beta
515		#[allow(clippy::many_single_char_names)] // Standard math notation for beta function
516	6	fn cf_beta(x: f64, a: f64, b: f64) -> f64 {
517	6	let max_iter = 100;
518	6	let eps = 1e-10;
519	6	let tiny = 1e-30;
520
521	6	let qab = a + b;
522	6	let qap = a + 1.0;
523	6	let qam = a - 1.0;
524
525	6	let mut c = 1.0;
526	6	let mut d = 1.0 - qab * x / qap;
527	6	if d.abs() < tiny {
528	0	d = tiny;
529	6	}
530	6	d = 1.0 / d;
531	6	let mut h = d;
532
533	19	for m in 1..= max_iter6 {
534	19	let m_f = m as f64;
535	19	let m2 = 2.0 * m_f;
536
537		// Even step
538	19	let aa = m_f * (b - m_f) * x / ((qam + m2) * (a + m2));
539	19	d = 1.0 + aa * d;
540	19	if d.abs() < tiny {
541	0	d = tiny;
542	19	}
543	19	c = 1.0 + aa / c;
544	19	if c.abs() < tiny {
545	0	c = tiny;
546	19	}
547	19	d = 1.0 / d;
548	19	h = d c;
549
550		// Odd step
551	19	let aa = -(a + m_f) * (qab + m_f) * x / ((a + m2) * (qap + m2));
552	19	d = 1.0 + aa * d;
553	19	if d.abs() < tiny {
554	0	d = tiny;
555	19	}
556	19	c = 1.0 + aa / c;
557	19	if c.abs() < tiny {
558	0	c = tiny;
559	19	}
560	19	d = 1.0 / d;
561	19	let del = d * c;
562	19	h *= del;
563
564	19	if (del - 1.0).abs() < eps {
565	6	break;
566	13	}
567		}
568
569	6	h
570	6	}
571
572		/// Approximate log-gamma function (Stirling's approximation)
573		#[allow(clippy::excessive_precision)] // Lanczos coefficients require high precision
574	18	fn gamma_ln(x: f64) -> f64 {
575	18	if x <= 0.0 {
576	0	return f64::INFINITY;
577	18	}
578
579		// Lanczos approximation coefficients
580	18	let g = 7.0;
581	18	let c = [
582	18	0.999_999_999_999_81,
583	18	676.520_368_121_885,
584	18	-1_259.139_216_722_403,
585	18	771.323_428_777_653,
586	18	-176.615_029_162_141,
587	18	12.507_343_278_687,
588	18	-0.138_571_095_265_72,
589	18	9.984_369_578_02e-6,
590	18	1.505_632_735_15e-7,
591	18	];
592
593	18	let x = x - 1.0;
594	18	let mut sum = c[0];
595	144	for (i, &coef) in c18 . iter18 (). enumerate18 (). skip18 (1) {
596	144	sum += coef / (x + i as f64);
597	144	}
598
599	18	let t = x + g + 0.5;
600	18	0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + sum.ln()
601	18	}