fdars_core/alignment/
diagnostics.rs

1//! Registration failure detection and alignment diagnostics.
2
3use super::quality::{warp_complexity, warp_smoothness};
4use super::{AlignmentResult, KarcherMeanResult};
5use crate::error::FdarError;
6use crate::helpers::simpsons_weights;
7use crate::matrix::FdMatrix;
8
9// ─── Types ───────────────────────────────────────────────────────────────────
10
11/// Diagnostic information for a single curve's alignment.
12#[derive(Debug, Clone, PartialEq)]
13#[non_exhaustive]
14pub struct AlignmentDiagnostic {
15    /// Index of the curve in the original dataset (or 0 for pairwise).
16    pub curve_index: usize,
17    /// Geodesic distance from the warp to the identity.
18    pub warp_complexity: f64,
19    /// Bending energy of the warp.
20    pub warp_smoothness: f64,
21    /// True if the residual is barely reduced (possible under-alignment).
22    pub is_under_aligned: bool,
23    /// True if warp complexity exceeds the threshold (possible over-alignment).
24    pub is_over_aligned: bool,
25    /// True if the warp contains a non-monotone segment.
26    pub has_non_monotone: bool,
27    /// Post-alignment L2 residual (weighted).
28    pub residual: f64,
29    /// Ratio of post-alignment residual to pre-alignment distance.
30    pub distance_ratio: f64,
31    /// True if any issue was detected.
32    pub flagged: bool,
33    /// Human-readable issue descriptions.
34    pub issues: Vec<String>,
35}
36
37/// Configuration for alignment diagnostics.
38///
39/// Construct via `DiagnosticConfig::default()`, then assign the fields you need (e.g. `let mut c = DiagnosticConfig::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.
40#[non_exhaustive]
41#[derive(Debug, Clone, PartialEq)]
42pub struct DiagnosticConfig {
43    /// Warp complexity above which the curve is flagged as over-aligned.
44    pub over_alignment_threshold: f64,
45    /// Distance ratio below which the curve is flagged as under-aligned
46    /// (i.e. the alignment barely improved the fit).
47    pub under_alignment_threshold: f64,
48    /// Maximum bending energy before the warp is considered too irregular.
49    pub max_bending_energy: f64,
50    /// Minimum improvement ratio (residual / pre-distance) to avoid flagging.
51    pub min_improvement_ratio: f64,
52}
53
54impl Default for DiagnosticConfig {
55    fn default() -> Self {
56        Self {
57            over_alignment_threshold: 1.0,
58            under_alignment_threshold: 1e-6,
59            max_bending_energy: 100.0,
60            min_improvement_ratio: 0.5,
61        }
62    }
63}
64
65/// Summary of diagnostics across all curves.
66#[derive(Debug, Clone, PartialEq)]
67#[non_exhaustive]
68pub struct AlignmentDiagnosticSummary {
69    /// Per-curve diagnostics.
70    pub diagnostics: Vec<AlignmentDiagnostic>,
71    /// Indices of flagged curves.
72    pub flagged_indices: Vec<usize>,
73    /// Number of flagged curves.
74    pub n_flagged: usize,
75    /// Overall health score in [0, 1]: fraction of curves that are *not* flagged.
76    pub health_score: f64,
77}
78
79// ─── Helpers ─────────────────────────────────────────────────────────────────
80
81/// Weighted L2 distance between two slices using pre-computed Simpson weights.
82fn weighted_l2(a: &[f64], b: &[f64], weights: &[f64]) -> f64 {
83    let mut sum = 0.0;
84    for i in 0..a.len() {
85        let d = a[i] - b[i];
86        sum += d * d * weights[i];
87    }
88    sum.sqrt()
89}
90
91/// Check monotonicity of a warp: returns true if any gamma[j+1] < gamma[j].
92fn is_non_monotone(gamma: &[f64]) -> bool {
93    gamma.windows(2).any(|w| w[1] < w[0])
94}
95
96/// Build a diagnostic for one curve given its warp, pre-distance, and residual.
97fn build_diagnostic(
98    curve_index: usize,
99    gamma: &[f64],
100    argvals: &[f64],
101    pre_distance: f64,
102    residual: f64,
103    config: &DiagnosticConfig,
104) -> AlignmentDiagnostic {
105    let wc = warp_complexity(gamma, argvals);
106    let ws = warp_smoothness(gamma, argvals);
107    let non_mono = is_non_monotone(gamma);
108
109    let distance_ratio = if pre_distance > 1e-15 {
110        residual / pre_distance
111    } else {
112        0.0
113    };
114
115    let is_over = wc > config.over_alignment_threshold;
116    let is_under = distance_ratio > config.min_improvement_ratio
117        && pre_distance > config.under_alignment_threshold;
118
119    let mut issues = Vec::new();
120    if is_over {
121        issues.push(format!(
122            "warp complexity {wc:.4} exceeds threshold {}",
123            config.over_alignment_threshold
124        ));
125    }
126    if is_under {
127        issues.push(format!(
128            "distance ratio {distance_ratio:.4} exceeds improvement threshold {}",
129            config.min_improvement_ratio
130        ));
131    }
132    if non_mono {
133        issues.push("warp contains non-monotone segments".to_string());
134    }
135    if ws > config.max_bending_energy {
136        issues.push(format!(
137            "bending energy {ws:.2} exceeds threshold {}",
138            config.max_bending_energy
139        ));
140    }
141
142    let flagged = !issues.is_empty();
143
144    AlignmentDiagnostic {
145        curve_index,
146        warp_complexity: wc,
147        warp_smoothness: ws,
148        is_under_aligned: is_under,
149        is_over_aligned: is_over,
150        has_non_monotone: non_mono,
151        residual,
152        distance_ratio,
153        flagged,
154        issues,
155    }
156}
157
158// ─── Public API ──────────────────────────────────────────────────────────────
159
160/// Diagnose alignment quality for every curve after a Karcher mean computation.
161///
162/// For each curve the function computes warp complexity, smoothness, pre- and
163/// post-alignment residuals, and checks for non-monotone warps and insufficient
164/// improvement. Curves with any issue are flagged.
165///
166/// # Arguments
167/// * `data`    — Original (unaligned) functional data (n x m).
168/// * `karcher` — Result of [`super::karcher::karcher_mean`].
169/// * `argvals` — Evaluation grid (length m).
170/// * `config`  — Diagnostic thresholds.
171///
172/// # Errors
173/// Returns `FdarError::InvalidDimension` on shape mismatches.
174pub fn diagnose_alignment(
175    data: &FdMatrix,
176    karcher: &KarcherMeanResult,
177    argvals: &[f64],
178    config: &DiagnosticConfig,
179) -> Result<AlignmentDiagnosticSummary, FdarError> {
180    let (n, m) = data.shape();
181
182    if argvals.len() != m {
183        return Err(FdarError::InvalidDimension {
184            parameter: "argvals",
185            expected: format!("{m}"),
186            actual: format!("{}", argvals.len()),
187        });
188    }
189    if karcher.gammas.nrows() != n || karcher.gammas.ncols() != m {
190        return Err(FdarError::InvalidDimension {
191            parameter: "karcher.gammas",
192            expected: format!("{n} x {m}"),
193            actual: format!("{} x {}", karcher.gammas.nrows(), karcher.gammas.ncols()),
194        });
195    }
196    if karcher.aligned_data.nrows() != n || karcher.aligned_data.ncols() != m {
197        return Err(FdarError::InvalidDimension {
198            parameter: "karcher.aligned_data",
199            expected: format!("{n} x {m}"),
200            actual: format!(
201                "{} x {}",
202                karcher.aligned_data.nrows(),
203                karcher.aligned_data.ncols()
204            ),
205        });
206    }
207    if karcher.mean.len() != m {
208        return Err(FdarError::InvalidDimension {
209            parameter: "karcher.mean",
210            expected: format!("{m}"),
211            actual: format!("{}", karcher.mean.len()),
212        });
213    }
214
215    let weights = simpsons_weights(argvals);
216
217    let mut diagnostics = Vec::with_capacity(n);
218    let mut flagged_indices = Vec::new();
219
220    for i in 0..n {
221        let gamma_i: Vec<f64> = (0..m).map(|j| karcher.gammas[(i, j)]).collect();
222
223        // Pre-alignment distance: ||f_i - mean||
224        let fi = data.row(i);
225        let pre_distance = weighted_l2(&fi, &karcher.mean, &weights);
226
227        // Post-alignment residual: ||f_i_aligned - mean||
228        let fi_aligned = karcher.aligned_data.row(i);
229        let residual = weighted_l2(&fi_aligned, &karcher.mean, &weights);
230
231        let diag = build_diagnostic(i, &gamma_i, argvals, pre_distance, residual, config);
232        if diag.flagged {
233            flagged_indices.push(i);
234        }
235        diagnostics.push(diag);
236    }
237
238    let n_flagged = flagged_indices.len();
239    let health_score = if n > 0 {
240        1.0 - n_flagged as f64 / n as f64
241    } else {
242        1.0
243    };
244
245    Ok(AlignmentDiagnosticSummary {
246        diagnostics,
247        flagged_indices,
248        n_flagged,
249        health_score,
250    })
251}
252
253/// Diagnose a single pairwise alignment.
254///
255/// Examines the warp produced by [`super::pairwise::elastic_align_pair`] and
256/// checks for over-alignment, under-alignment, non-monotonicity, and excessive
257/// bending energy.
258pub fn diagnose_pairwise(
259    f1: &[f64],
260    f2: &[f64],
261    result: &AlignmentResult,
262    argvals: &[f64],
263    config: &DiagnosticConfig,
264) -> AlignmentDiagnostic {
265    let weights = simpsons_weights(argvals);
266
267    // Pre-alignment L2 distance
268    let pre_distance = weighted_l2(f1, f2, &weights);
269
270    // Post-alignment residual: ||f1 - f2_aligned||
271    let residual = weighted_l2(f1, &result.f_aligned, &weights);
272
273    build_diagnostic(0, &result.gamma, argvals, pre_distance, residual, config)
274}
275
276#[cfg(test)]
277mod tests {
278    use super::*;
279    use crate::alignment::karcher_mean;
280    use crate::alignment::pairwise::elastic_align_pair;
281    use crate::simulation::{sim_fundata, EFunType, EValType};
282    use crate::test_helpers::uniform_grid;
283
284    fn make_data(n: usize, m: usize) -> (FdMatrix, Vec<f64>) {
285        let t = uniform_grid(m);
286        let data = sim_fundata(n, &t, 3, EFunType::Fourier, EValType::Exponential, Some(99));
287        (data, t)
288    }
289
290    #[test]
291    fn diagnose_alignment_smoke() {
292        let (data, t) = make_data(8, 30);
293        let km = karcher_mean(&data, &t, 5, 1e-2, 0.0);
294        let config = DiagnosticConfig::default();
295        let summary = diagnose_alignment(&data, &km, &t, &config).unwrap();
296        assert_eq!(summary.diagnostics.len(), 8);
297        assert!(summary.health_score >= 0.0 && summary.health_score <= 1.0);
298        assert_eq!(summary.n_flagged, summary.flagged_indices.len());
299    }
300
301    #[test]
302    fn diagnose_alignment_identical_returns_low_complexity() {
303        // When data is identical, warp complexity should be small even though
304        // post-centering numerics may not yield exactly the identity warp.
305        let t = uniform_grid(30);
306        let curve: Vec<f64> = t.iter().map(|&x| x.sin()).collect();
307        let mut vals = Vec::with_capacity(5 * 30);
308        for _ in 0..5 {
309            vals.extend_from_slice(&curve);
310        }
311        let data = FdMatrix::from_column_major(vals, 5, 30).unwrap();
312        let km = karcher_mean(&data, &t, 5, 1e-3, 0.0);
313        let config = DiagnosticConfig::default();
314        let summary = diagnose_alignment(&data, &km, &t, &config).unwrap();
315        assert_eq!(summary.diagnostics.len(), 5);
316        // All warp complexities should be small (near identity)
317        for d in &summary.diagnostics {
318            assert!(
319                d.warp_complexity < 0.5,
320                "curve {} warp_complexity {} should be small for identical data",
321                d.curve_index,
322                d.warp_complexity,
323            );
324        }
325    }
326
327    #[test]
328    fn diagnose_alignment_rejects_shape_mismatch() {
329        let (data, t) = make_data(6, 30);
330        let km = karcher_mean(&data, &t, 3, 1e-2, 0.0);
331        let bad_t = uniform_grid(20);
332        let config = DiagnosticConfig::default();
333        assert!(diagnose_alignment(&data, &km, &bad_t, &config).is_err());
334    }
335
336    #[test]
337    fn diagnose_pairwise_smoke() {
338        let t = uniform_grid(30);
339        let f1: Vec<f64> = t.iter().map(|&x| (x * 6.0).sin()).collect();
340        let f2: Vec<f64> = t.iter().map(|&x| ((x + 0.15) * 6.0).sin()).collect();
341        let alignment = elastic_align_pair(&f1, &f2, &t, 0.0);
342        let config = DiagnosticConfig::default();
343        let diag = diagnose_pairwise(&f1, &f2, &alignment, &t, &config);
344        assert!(diag.warp_complexity >= 0.0);
345        assert!(diag.residual >= 0.0);
346    }
347
348    #[test]
349    fn diagnose_pairwise_identical() {
350        let t = uniform_grid(30);
351        let f: Vec<f64> = t.iter().map(|&x| x.sin()).collect();
352        let alignment = elastic_align_pair(&f, &f, &t, 0.0);
353        let config = DiagnosticConfig::default();
354        let diag = diagnose_pairwise(&f, &f, &alignment, &t, &config);
355        assert!(
356            diag.residual < 1e-3,
357            "identical curves should have near-zero residual"
358        );
359        assert!(!diag.has_non_monotone);
360    }
361}
fdars_core/alignment/diagnostics.rs

fdars_core/alignment/
diagnostics.rs
