codec_eval/stats/
mod.rs

1//! Statistical analysis and Pareto front calculation.
2//!
3//! This module provides tools for analyzing codec comparison results,
4//! including Pareto front calculation for rate-distortion analysis.
5//!
6//! ## Core Statistics
7//!
8//! - [`Summary`]: Descriptive statistics (mean, median, std_dev, percentiles)
9//! - [`median`], [`mean`], [`std_dev`]: Basic statistical functions
10//! - [`percentile`], [`percentile_u32`]: Percentile calculation (R-7 interpolation)
11//! - [`trimmed_mean`]: Robust mean excluding outliers
12//! - [`iqr`]: Interquartile range
13//!
14//! ## Rate-Distortion Analysis
15//!
16//! - [`bd_rate`]: Bjontegaard Delta Rate calculation
17//! - [`ParetoFront`]: Pareto-optimal points on RD curve
18
19#[cfg(feature = "chart")]
20pub mod chart;
21mod pareto;
22pub mod rd_knee;
23
24#[cfg(feature = "chart")]
25pub use chart::{ChartConfig, ChartPoint, ChartSeries, generate_svg};
26pub use pareto::{ParetoFront, RDPoint};
27pub use rd_knee::{
28    AngleBin, AxisRange, BinScheme, CodecConfig, ConfiguredParetoFront, ConfiguredRDPoint,
29    CorpusAggregate, DualAngleBin, EncodeResult, FixedFrame, NormalizationContext, ParamValue,
30    QualityDirection, RDCalibration, RDKnee, RDPosition, plot_rd_svg,
31};
32
33use serde::{Deserialize, Serialize};
34
35/// Descriptive statistics for a set of measurements.
36#[derive(Debug, Clone, Serialize, Deserialize)]
37pub struct Summary {
38    /// Number of values.
39    pub count: usize,
40    /// Mean value.
41    pub mean: f64,
42    /// Median value.
43    pub median: f64,
44    /// Standard deviation.
45    pub std_dev: f64,
46    /// Minimum value.
47    pub min: f64,
48    /// Maximum value.
49    pub max: f64,
50    /// 5th percentile.
51    pub p5: f64,
52    /// 25th percentile.
53    pub p25: f64,
54    /// 75th percentile.
55    pub p75: f64,
56    /// 95th percentile.
57    pub p95: f64,
58}
59
60impl Summary {
61    /// Compute summary statistics for a slice of values.
62    ///
63    /// Returns `None` if the slice is empty.
64    #[must_use]
65    pub fn compute(values: &[f64]) -> Option<Self> {
66        if values.is_empty() {
67            return None;
68        }
69
70        let mut sorted = values.to_vec();
71        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
72
73        let count = sorted.len();
74        let sum: f64 = sorted.iter().sum();
75        let mean = sum / count as f64;
76
77        let variance: f64 = sorted.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / count as f64;
78        let std_dev = variance.sqrt();
79
80        let median = percentile_sorted(&sorted, 0.5);
81        let min = sorted[0];
82        let max = sorted[count - 1];
83
84        Some(Self {
85            count,
86            mean,
87            median,
88            std_dev,
89            min,
90            max,
91            p5: percentile_sorted(&sorted, 0.05),
92            p25: percentile_sorted(&sorted, 0.25),
93            p75: percentile_sorted(&sorted, 0.75),
94            p95: percentile_sorted(&sorted, 0.95),
95        })
96    }
97}
98
99//=============================================================================
100// Core Statistical Functions
101//=============================================================================
102
103/// Compute median of a slice.
104///
105/// For even-length slices, returns the average of the two middle values.
106///
107/// # Example
108///
109/// ```
110/// use codec_eval::stats::median;
111///
112/// assert_eq!(median(&[1.0, 2.0, 3.0, 4.0, 5.0]), 3.0);
113/// assert_eq!(median(&[1.0, 2.0, 3.0, 4.0]), 2.5);
114/// ```
115#[must_use]
116pub fn median(values: &[f64]) -> f64 {
117    if values.is_empty() {
118        return 0.0;
119    }
120    let mut sorted = values.to_vec();
121    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
122    let mid = sorted.len() / 2;
123    if sorted.len() % 2 == 0 {
124        (sorted[mid - 1] + sorted[mid]) / 2.0
125    } else {
126        sorted[mid]
127    }
128}
129
130/// Compute arithmetic mean.
131///
132/// # Example
133///
134/// ```
135/// use codec_eval::stats::mean;
136///
137/// assert!((mean(&[1.0, 2.0, 3.0, 4.0, 5.0]) - 3.0).abs() < 0.001);
138/// ```
139#[must_use]
140pub fn mean(values: &[f64]) -> f64 {
141    if values.is_empty() {
142        return 0.0;
143    }
144    values.iter().sum::<f64>() / values.len() as f64
145}
146
147/// Compute sample standard deviation.
148///
149/// Uses Bessel's correction (N-1 denominator) for sample standard deviation.
150///
151/// # Example
152///
153/// ```
154/// use codec_eval::stats::std_dev;
155///
156/// let values = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
157/// assert!((std_dev(&values) - 2.138).abs() < 0.001);
158/// ```
159#[must_use]
160pub fn std_dev(values: &[f64]) -> f64 {
161    if values.len() < 2 {
162        return 0.0;
163    }
164    let m = mean(values);
165    let variance = values.iter().map(|x| (x - m).powi(2)).sum::<f64>() / (values.len() - 1) as f64;
166    variance.sqrt()
167}
168
169/// Compute percentile using linear interpolation (R-7 method).
170///
171/// This is the default method used by R, NumPy, and Excel.
172/// The percentile `p` should be in the range 0.0 to 1.0.
173///
174/// # Example
175///
176/// ```
177/// use codec_eval::stats::percentile;
178///
179/// let values = [1.0, 2.0, 3.0, 4.0, 5.0];
180/// assert!((percentile(&values, 0.5) - 3.0).abs() < 0.001);  // median
181/// assert!((percentile(&values, 0.25) - 2.0).abs() < 0.001); // Q1
182/// assert!((percentile(&values, 0.75) - 4.0).abs() < 0.001); // Q3
183/// ```
184#[must_use]
185pub fn percentile(values: &[f64], p: f64) -> f64 {
186    if values.is_empty() {
187        return 0.0;
188    }
189    let mut sorted = values.to_vec();
190    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
191    percentile_sorted(&sorted, p)
192}
193
194/// Compute percentile for u32 values.
195///
196/// Returns the interpolated percentile rounded to the nearest integer.
197///
198/// # Example
199///
200/// ```
201/// use codec_eval::stats::percentile_u32;
202///
203/// let values = [10, 20, 30, 40, 50];
204/// assert_eq!(percentile_u32(&values, 0.5), 30);
205/// ```
206#[must_use]
207pub fn percentile_u32(values: &[u32], p: f64) -> u32 {
208    if values.is_empty() {
209        return 0;
210    }
211    let mut sorted = values.to_vec();
212    sorted.sort();
213    let pos = p.clamp(0.0, 1.0) * (sorted.len() - 1) as f64;
214    let lower = pos.floor() as usize;
215    let upper = (lower + 1).min(sorted.len() - 1);
216    let frac = pos - lower as f64;
217    let result = sorted[lower] as f64 * (1.0 - frac) + sorted[upper] as f64 * frac;
218    result.round() as u32
219}
220
221/// Compute trimmed mean (removes top and bottom `trim_pct` of values).
222///
223/// This provides a robust estimate of central tendency that is less
224/// affected by outliers than the arithmetic mean.
225///
226/// # Arguments
227///
228/// * `values` - The values to analyze
229/// * `trim_pct` - Fraction to trim from each end (e.g., 0.1 = 10% from each end)
230///
231/// # Example
232///
233/// ```
234/// use codec_eval::stats::trimmed_mean;
235///
236/// // Outliers at ends don't affect trimmed mean much
237/// let values = [1.0, 10.0, 11.0, 12.0, 13.0, 100.0];
238/// let tm = trimmed_mean(&values, 0.2);  // Trim 20% from each end
239/// assert!((tm - 11.5).abs() < 0.001);   // Average of middle values
240/// ```
241#[must_use]
242pub fn trimmed_mean(values: &[f64], trim_pct: f64) -> f64 {
243    if values.is_empty() {
244        return 0.0;
245    }
246    let mut sorted = values.to_vec();
247    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
248    let trim_count = (sorted.len() as f64 * trim_pct.clamp(0.0, 0.5)) as usize;
249    if trim_count * 2 >= sorted.len() {
250        return median(values);
251    }
252    let trimmed = &sorted[trim_count..sorted.len() - trim_count];
253    mean(trimmed)
254}
255
256/// Compute interquartile range (IQR = Q3 - Q1).
257///
258/// The IQR is a robust measure of spread that is not affected by outliers.
259///
260/// # Example
261///
262/// ```
263/// use codec_eval::stats::iqr;
264///
265/// let values = [1.0, 2.0, 3.0, 4.0, 5.0];
266/// assert!((iqr(&values) - 2.0).abs() < 0.001);  // Q3(4) - Q1(2) = 2
267/// ```
268#[must_use]
269pub fn iqr(values: &[f64]) -> f64 {
270    percentile(values, 0.75) - percentile(values, 0.25)
271}
272
273/// Internal: Calculate percentile from pre-sorted values.
274/// Accepts percentile in 0-100 range for backward compatibility with Summary.
275fn percentile_sorted(sorted: &[f64], p: f64) -> f64 {
276    if sorted.is_empty() {
277        return 0.0;
278    }
279    if sorted.len() == 1 {
280        return sorted[0];
281    }
282
283    // Normalize to 0-1 range if given as percentage
284    let p = if p > 1.0 { p / 100.0 } else { p };
285    let p = p.clamp(0.0, 1.0);
286
287    let idx = p * (sorted.len() - 1) as f64;
288    let lower = idx.floor() as usize;
289    let upper = idx.ceil() as usize;
290    let frac = idx - lower as f64;
291
292    if lower == upper {
293        sorted[lower]
294    } else {
295        sorted[lower] * (1.0 - frac) + sorted[upper] * frac
296    }
297}
298
299/// Calculate BD-Rate (Bjontegaard Delta Rate).
300///
301/// BD-Rate measures the average bitrate difference between two rate-distortion
302/// curves at the same quality level. A negative value means the test curve
303/// is more efficient (lower bitrate at same quality).
304///
305/// # Arguments
306///
307/// * `reference` - Reference curve points (bitrate, quality).
308/// * `test` - Test curve points (bitrate, quality).
309///
310/// # Returns
311///
312/// BD-Rate as a percentage. Negative = test is better.
313#[must_use]
314pub fn bd_rate(reference: &[(f64, f64)], test: &[(f64, f64)]) -> Option<f64> {
315    if reference.len() < 4 || test.len() < 4 {
316        return None;
317    }
318
319    // Sort by quality
320    let mut ref_sorted: Vec<_> = reference.to_vec();
321    let mut test_sorted: Vec<_> = test.to_vec();
322    ref_sorted.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
323    test_sorted.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
324
325    // Find overlapping quality range
326    let min_quality = ref_sorted[0].1.max(test_sorted[0].1);
327    let max_quality = ref_sorted.last()?.1.min(test_sorted.last()?.1);
328
329    if min_quality >= max_quality {
330        return None;
331    }
332
333    // Use log-rate for integration
334    let ref_log: Vec<_> = ref_sorted.iter().map(|(r, q)| (r.ln(), *q)).collect();
335    let test_log: Vec<_> = test_sorted.iter().map(|(r, q)| (r.ln(), *q)).collect();
336
337    // Integrate using trapezoidal rule (simplified)
338    let ref_area = integrate_curve(&ref_log, min_quality, max_quality);
339    let test_area = integrate_curve(&test_log, min_quality, max_quality);
340
341    let avg_ref = ref_area / (max_quality - min_quality);
342    let avg_test = test_area / (max_quality - min_quality);
343
344    // BD-Rate = (10^(avg_test - avg_ref) - 1) * 100
345    let bd = (10_f64.powf(avg_test - avg_ref) - 1.0) * 100.0;
346
347    Some(bd)
348}
349
350/// Simple trapezoidal integration of a curve.
351fn integrate_curve(points: &[(f64, f64)], min_x: f64, max_x: f64) -> f64 {
352    let mut area = 0.0;
353
354    for window in points.windows(2) {
355        let (y0, x0) = window[0];
356        let (y1, x1) = window[1];
357
358        // Skip segments outside range
359        if x1 < min_x || x0 > max_x {
360            continue;
361        }
362
363        // Clip to range
364        let x0 = x0.max(min_x);
365        let x1 = x1.min(max_x);
366
367        // Trapezoidal area
368        area += (y0 + y1) / 2.0 * (x1 - x0);
369    }
370
371    area
372}
373
374#[cfg(test)]
375mod tests {
376    use super::*;
377
378    #[test]
379    fn test_summary_compute() {
380        let values = vec![1.0, 2.0, 3.0, 4.0, 5.0];
381        let summary = Summary::compute(&values).unwrap();
382
383        assert_eq!(summary.count, 5);
384        assert!((summary.mean - 3.0).abs() < 0.001);
385        assert!((summary.median - 3.0).abs() < 0.001);
386        assert!((summary.min - 1.0).abs() < 0.001);
387        assert!((summary.max - 5.0).abs() < 0.001);
388    }
389
390    #[test]
391    fn test_summary_empty() {
392        assert!(Summary::compute(&[]).is_none());
393    }
394
395    #[test]
396    fn test_percentile() {
397        let values = vec![1.0, 2.0, 3.0, 4.0, 5.0];
398        assert!((percentile(&values, 0.0) - 1.0).abs() < 0.001);
399        assert!((percentile(&values, 0.5) - 3.0).abs() < 0.001);
400        assert!((percentile(&values, 1.0) - 5.0).abs() < 0.001);
401    }
402
403    #[test]
404    fn test_median() {
405        assert_eq!(median(&[1.0, 2.0, 3.0, 4.0, 5.0]), 3.0);
406        assert_eq!(median(&[1.0, 2.0, 3.0, 4.0]), 2.5);
407        assert_eq!(median(&[5.0]), 5.0);
408        assert_eq!(median(&[]), 0.0);
409    }
410
411    #[test]
412    fn test_trimmed_mean() {
413        // With outliers
414        let values = [1.0, 10.0, 11.0, 12.0, 13.0, 100.0];
415        let tm = trimmed_mean(&values, 0.2);
416        assert!((tm - 11.5).abs() < 0.001);
417
418        // Regular case
419        let values2 = [1.0, 2.0, 3.0, 4.0, 5.0];
420        let tm2 = trimmed_mean(&values2, 0.2);
421        assert!((tm2 - 3.0).abs() < 0.001);
422    }
423
424    #[test]
425    fn test_iqr() {
426        let values = [1.0, 2.0, 3.0, 4.0, 5.0];
427        assert!((iqr(&values) - 2.0).abs() < 0.001);
428    }
429
430    #[test]
431    fn test_percentile_u32() {
432        let values = [10, 20, 30, 40, 50];
433        assert_eq!(percentile_u32(&values, 0.0), 10);
434        assert_eq!(percentile_u32(&values, 0.5), 30);
435        assert_eq!(percentile_u32(&values, 1.0), 50);
436    }
437
438    #[test]
439    fn test_bd_rate_same_curve() {
440        let curve = vec![
441            (1000.0, 30.0),
442            (2000.0, 35.0),
443            (4000.0, 40.0),
444            (8000.0, 45.0),
445        ];
446
447        let bd = bd_rate(&curve, &curve);
448        assert!(bd.is_some());
449        assert!(bd.unwrap().abs() < 0.1); // Should be ~0 for same curve
450    }
451}
codec_eval/stats/mod.rs

codec_eval/stats/
mod.rs
