fdars_core/alignment/
outlier.rs

1//! SRVF-based outlier detection using elastic distances.
2//!
3//! Detects outlier curves by computing elastic distances from a reference
4//! (Karcher mean or median) and applying the Tukey fence rule.
5
6use super::karcher::karcher_mean;
7use super::pairwise::{amplitude_distance, elastic_distance, phase_distance_pair};
8use super::robust_karcher::{karcher_median, RobustKarcherConfig};
9use crate::error::FdarError;
10use crate::matrix::FdMatrix;
11
12/// Configuration for elastic outlier detection.
13///
14/// Construct via `ElasticOutlierConfig::default()`, then assign the fields you need (e.g. `let mut c = ElasticOutlierConfig::default(); c.field = …;`). This struct is `#[non_exhaustive]`, so external crates cannot build it with a struct literal — not even functional-update `..Default::default()` form.
15#[non_exhaustive]
16#[derive(Debug, Clone, PartialEq)]
17pub struct ElasticOutlierConfig {
18    /// Roughness penalty for elastic alignment (0.0 = no penalty).
19    pub lambda: f64,
20    /// Significance level (controls threshold sensitivity; currently used
21    /// to document intent — the actual threshold uses the Tukey fence).
22    pub alpha: f64,
23    /// If `true`, use the Karcher median as reference (more robust).
24    /// If `false`, use the Karcher mean.
25    pub use_median: bool,
26}
27
28impl Default for ElasticOutlierConfig {
29    fn default() -> Self {
30        Self {
31            lambda: 0.0,
32            alpha: 0.05,
33            use_median: true,
34        }
35    }
36}
37
38/// Result of elastic outlier detection.
39#[derive(Debug, Clone, PartialEq)]
40#[non_exhaustive]
41pub struct ElasticOutlierResult {
42    /// Indices of detected outlier curves.
43    pub outlier_indices: Vec<usize>,
44    /// Total elastic distance from each curve to the reference (length n).
45    pub distances: Vec<f64>,
46    /// Cutoff distance (Tukey fence: Q3 + 1.5 * IQR).
47    pub threshold: f64,
48    /// Amplitude component of elastic distance for each curve (length n).
49    pub amplitude_distances: Vec<f64>,
50    /// Phase component of elastic distance for each curve (length n).
51    pub phase_distances: Vec<f64>,
52}
53
54/// Detect outlier curves using elastic distances and the Tukey fence.
55///
56/// Computes a reference curve (Karcher mean or median), then measures the
57/// elastic distance from each curve to the reference. Curves exceeding the
58/// Tukey fence threshold (Q3 + 1.5 * IQR) are flagged as outliers.
59///
60/// # Arguments
61/// * `data`    — Functional data matrix (n x m).
62/// * `argvals` — Evaluation points (length m).
63/// * `config`  — Configuration parameters.
64///
65/// # Errors
66/// Returns [`FdarError::InvalidDimension`] if `argvals` length does not match `m`
67/// or `n < 2`.
68#[must_use = "expensive computation whose result should not be discarded"]
69pub fn elastic_outlier_detection(
70    data: &FdMatrix,
71    argvals: &[f64],
72    config: &ElasticOutlierConfig,
73) -> Result<ElasticOutlierResult, FdarError> {
74    let (n, m) = data.shape();
75
76    if argvals.len() != m {
77        return Err(FdarError::InvalidDimension {
78            parameter: "argvals",
79            expected: format!("{m}"),
80            actual: format!("{}", argvals.len()),
81        });
82    }
83    if n < 2 {
84        return Err(FdarError::InvalidDimension {
85            parameter: "data",
86            expected: "at least 2 rows".to_string(),
87            actual: format!("{n} rows"),
88        });
89    }
90
91    // Step 1: Compute reference curve.
92    let reference = if config.use_median {
93        let median_config = RobustKarcherConfig {
94            max_iter: 15,
95            tol: 1e-3,
96            lambda: config.lambda,
97            trim_fraction: 0.1,
98        };
99        let result = karcher_median(data, argvals, &median_config)?;
100        result.mean
101    } else {
102        let result = karcher_mean(data, argvals, 15, 1e-3, config.lambda);
103        result.mean
104    };
105
106    // Step 2: Compute elastic distances from reference.
107    let distances: Vec<f64> = (0..n)
108        .map(|i| {
109            let fi = data.row(i);
110            elastic_distance(&reference, &fi, argvals, config.lambda)
111        })
112        .collect();
113
114    // Step 3: Compute amplitude and phase distances.
115    let amplitude_distances: Vec<f64> = (0..n)
116        .map(|i| {
117            let fi = data.row(i);
118            amplitude_distance(&reference, &fi, argvals, config.lambda)
119        })
120        .collect();
121
122    let phase_distances: Vec<f64> = (0..n)
123        .map(|i| {
124            let fi = data.row(i);
125            phase_distance_pair(&reference, &fi, argvals, config.lambda)
126        })
127        .collect();
128
129    // Step 4: Tukey fence on total distances.
130    let threshold = tukey_fence(&distances);
131
132    // Step 5: Identify outliers.
133    let outlier_indices: Vec<usize> = (0..n).filter(|&i| distances[i] > threshold).collect();
134
135    Ok(ElasticOutlierResult {
136        outlier_indices,
137        distances,
138        threshold,
139        amplitude_distances,
140        phase_distances,
141    })
142}
143
144/// Compute the Tukey fence threshold: Q3 + 1.5 * IQR.
145fn tukey_fence(values: &[f64]) -> f64 {
146    let n = values.len();
147    if n < 4 {
148        // With very few values, use max + epsilon as a permissive threshold.
149        return values.iter().copied().fold(f64::NEG_INFINITY, f64::max) + 1.0;
150    }
151
152    let mut sorted = values.to_vec();
153    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
154
155    let q1 = percentile_sorted(&sorted, 25.0);
156    let q3 = percentile_sorted(&sorted, 75.0);
157    let iqr = q3 - q1;
158
159    q3 + 1.5 * iqr
160}
161
162/// Compute a percentile from a sorted slice using linear interpolation.
163fn percentile_sorted(sorted: &[f64], pct: f64) -> f64 {
164    let n = sorted.len();
165    if n == 0 {
166        return 0.0;
167    }
168    if n == 1 {
169        return sorted[0];
170    }
171
172    let rank = pct / 100.0 * (n - 1) as f64;
173    let lo = rank.floor() as usize;
174    let hi = rank.ceil() as usize;
175    let frac = rank - lo as f64;
176
177    if lo >= n || hi >= n {
178        sorted[n - 1]
179    } else {
180        sorted[lo] * (1.0 - frac) + sorted[hi] * frac
181    }
182}
183
184#[cfg(test)]
185mod tests {
186    use super::*;
187    use crate::test_helpers::uniform_grid;
188
189    fn make_clean_data(n: usize, m: usize) -> (FdMatrix, Vec<f64>) {
190        let t = uniform_grid(m);
191        let mut data_vec = vec![0.0; n * m];
192        for i in 0..n {
193            let phase = 0.02 * i as f64;
194            for j in 0..m {
195                data_vec[i + j * n] = ((t[j] + phase) * 4.0).sin();
196            }
197        }
198        let data = FdMatrix::from_column_major(data_vec, n, m).unwrap();
199        (data, t)
200    }
201
202    #[test]
203    fn outlier_detection_no_outliers() {
204        let (data, t) = make_clean_data(6, 20);
205        let config = ElasticOutlierConfig::default();
206        let result = elastic_outlier_detection(&data, &t, &config).unwrap();
207
208        assert_eq!(result.distances.len(), 6);
209        assert_eq!(result.amplitude_distances.len(), 6);
210        assert_eq!(result.phase_distances.len(), 6);
211        assert!(result.threshold > 0.0);
212
213        // With clean homogeneous data, expect no (or very few) outliers.
214        assert!(
215            result.outlier_indices.len() <= 1,
216            "clean data should have at most 1 outlier, got {}",
217            result.outlier_indices.len()
218        );
219    }
220
221    #[test]
222    fn outlier_detection_finds_extreme() {
223        let m = 20;
224        let t = uniform_grid(m);
225        let n = 8;
226        let mut data_vec = vec![0.0; n * m];
227
228        // 7 clean curves.
229        for i in 0..7 {
230            let phase = 0.02 * i as f64;
231            for j in 0..m {
232                data_vec[i + j * n] = ((t[j] + phase) * 4.0).sin();
233            }
234        }
235        // 1 extreme outlier (curve 7).
236        for j in 0..m {
237            data_vec[7 + j * n] = (t[j] * 20.0).cos() * 10.0;
238        }
239        let data = FdMatrix::from_column_major(data_vec, n, m).unwrap();
240
241        let config = ElasticOutlierConfig::default();
242        let result = elastic_outlier_detection(&data, &t, &config).unwrap();
243
244        // The outlier (index 7) should be detected.
245        assert!(
246            result.outlier_indices.contains(&7),
247            "should detect curve 7 as outlier, detected: {:?}",
248            result.outlier_indices
249        );
250
251        // The outlier's distance should exceed the threshold.
252        assert!(
253            result.distances[7] > result.threshold,
254            "outlier distance ({}) should exceed threshold ({})",
255            result.distances[7],
256            result.threshold
257        );
258    }
259
260    #[test]
261    fn outlier_detection_config_default() {
262        let cfg = ElasticOutlierConfig::default();
263        assert!((cfg.lambda - 0.0).abs() < f64::EPSILON);
264        assert!((cfg.alpha - 0.05).abs() < f64::EPSILON);
265        assert!(cfg.use_median);
266    }
267}
fdars_core/alignment/outlier.rs

fdars_core/alignment/
outlier.rs
