fdars_core/alignment/srsf.rs
1//! SRSF (Square-Root Slope Function) transforms and warping utilities.
2
3use crate::fdata::deriv_1d;
4use crate::helpers::{cumulative_trapz, linear_interp};
5use crate::matrix::FdMatrix;
6
7// ─── SRSF Transform and Inverse ─────────────────────────────────────────────
8
9/// Compute the Square-Root Slope Function (SRSF) transform.
10///
11/// For each curve f, the SRSF is: `q(t) = sign(f'(t)) * sqrt(|f'(t)|)`
12///
13/// # Arguments
14/// * `data` — Functional data matrix (n × m)
15/// * `argvals` — Evaluation points (length m)
16///
17/// # Returns
18/// FdMatrix of SRSFs with the same shape as input.
19#[must_use = "expensive computation whose result should not be discarded"]
20pub fn srsf_transform(data: &FdMatrix, argvals: &[f64]) -> FdMatrix {
21 let (n, m) = data.shape();
22 if n == 0 || m == 0 || argvals.len() != m {
23 return FdMatrix::zeros(n, m);
24 }
25
26 let deriv = deriv_1d(data, argvals, 1);
27
28 let mut result = FdMatrix::zeros(n, m);
29 for i in 0..n {
30 for j in 0..m {
31 let d = deriv[(i, j)];
32 result[(i, j)] = d.signum() * d.abs().sqrt();
33 }
34 }
35 result
36}
37
38/// Reconstruct a curve from its SRSF representation.
39///
40/// Given SRSF q and initial value f0, reconstructs: `f(t) = f0 + ∫₀ᵗ q(s)|q(s)| ds`
41///
42/// # Arguments
43/// * `q` — SRSF values (length m)
44/// * `argvals` — Evaluation points (length m)
45/// * `f0` — Initial value f(argvals\[0\])
46///
47/// # Returns
48/// Reconstructed curve values.
49pub fn srsf_inverse(q: &[f64], argvals: &[f64], f0: f64) -> Vec<f64> {
50 let m = q.len();
51 if m == 0 {
52 return Vec::new();
53 }
54
55 // Integrand: q(s) * |q(s)|
56 let integrand: Vec<f64> = q.iter().map(|&qi| qi * qi.abs()).collect();
57 let integral = cumulative_trapz(&integrand, argvals);
58
59 integral.iter().map(|&v| f0 + v).collect()
60}
61
62// ─── Reparameterization ─────────────────────────────────────────────────────
63
64/// Reparameterize a curve by a warping function.
65///
66/// Computes `f(γ(t))` via linear interpolation.
67///
68/// # Arguments
69/// * `f` — Curve values (length m)
70/// * `argvals` — Evaluation points (length m)
71/// * `gamma` — Warping function values (length m)
72pub fn reparameterize_curve(f: &[f64], argvals: &[f64], gamma: &[f64]) -> Vec<f64> {
73 gamma
74 .iter()
75 .map(|&g| linear_interp(argvals, f, g))
76 .collect()
77}
78
79/// Compose two warping functions: `(γ₁ ∘ γ₂)(t) = γ₁(γ₂(t))`.
80///
81/// # Arguments
82/// * `gamma1` — Outer warping function (length m)
83/// * `gamma2` — Inner warping function (length m)
84/// * `argvals` — Evaluation points (length m)
85pub fn compose_warps(gamma1: &[f64], gamma2: &[f64], argvals: &[f64]) -> Vec<f64> {
86 gamma2
87 .iter()
88 .map(|&g| linear_interp(argvals, gamma1, g))
89 .collect()
90}
91
92/// Compute a single SRSF from a slice (single-row convenience).
93pub(crate) fn srsf_single(f: &[f64], argvals: &[f64]) -> Vec<f64> {
94 let m = f.len();
95 let mat = FdMatrix::from_slice(f, 1, m).expect("dimension invariant: data.len() == n * m");
96 let q_mat = srsf_transform(&mat, argvals);
97 q_mat.row(0)
98}