fdars_core/
helpers.rs

1//! Helper functions for numerical integration and common operations.
2
3/// Small epsilon for numerical comparisons (e.g., avoiding division by zero).
4pub const NUMERICAL_EPS: f64 = 1e-10;
5
6/// Default convergence tolerance for iterative algorithms.
7pub const DEFAULT_CONVERGENCE_TOL: f64 = 1e-6;
8
9/// Extract curves from column-major data matrix.
10///
11/// Converts a flat column-major matrix into a vector of curve vectors,
12/// where each curve contains all evaluation points for one observation.
13///
14/// # Arguments
15/// * `data` - Functional data matrix (n x m)
16///
17/// # Returns
18/// Vector of n curves, each containing m values
19pub fn extract_curves(data: &crate::matrix::FdMatrix) -> Vec<Vec<f64>> {
20    data.rows()
21}
22
23/// Compute L2 distance between two curves using integration weights.
24///
25/// # Arguments
26/// * `curve1` - First curve values
27/// * `curve2` - Second curve values
28/// * `weights` - Integration weights
29///
30/// # Returns
31/// L2 distance between the curves
32pub fn l2_distance(curve1: &[f64], curve2: &[f64], weights: &[f64]) -> f64 {
33    let mut dist_sq = 0.0;
34    for i in 0..curve1.len() {
35        let diff = curve1[i] - curve2[i];
36        dist_sq += diff * diff * weights[i];
37    }
38    dist_sq.sqrt()
39}
40
41/// Compute Simpson's rule integration weights for non-uniform grid.
42///
43/// Returns weights for trapezoidal rule integration.
44///
45/// # Arguments
46/// * `argvals` - Grid points (evaluation points)
47///
48/// # Returns
49/// Vector of integration weights
50pub fn simpsons_weights(argvals: &[f64]) -> Vec<f64> {
51    let n = argvals.len();
52    if n < 2 {
53        return vec![1.0; n];
54    }
55
56    let mut weights = vec![0.0; n];
57
58    if n == 2 {
59        // Trapezoidal rule
60        let h = argvals[1] - argvals[0];
61        weights[0] = h / 2.0;
62        weights[1] = h / 2.0;
63        return weights;
64    }
65
66    // For non-uniform spacing, use composite trapezoidal rule
67    for i in 0..n {
68        if i == 0 {
69            weights[i] = (argvals[1] - argvals[0]) / 2.0;
70        } else if i == n - 1 {
71            weights[i] = (argvals[n - 1] - argvals[n - 2]) / 2.0;
72        } else {
73            weights[i] = (argvals[i + 1] - argvals[i - 1]) / 2.0;
74        }
75    }
76
77    weights
78}
79
80/// Compute 2D integration weights using tensor product of 1D weights.
81///
82/// Returns a flattened vector of weights for an m1 x m2 grid.
83///
84/// # Arguments
85/// * `argvals_s` - Grid points in s direction
86/// * `argvals_t` - Grid points in t direction
87///
88/// # Returns
89/// Flattened vector of integration weights (row-major order)
90pub fn simpsons_weights_2d(argvals_s: &[f64], argvals_t: &[f64]) -> Vec<f64> {
91    let weights_s = simpsons_weights(argvals_s);
92    let weights_t = simpsons_weights(argvals_t);
93    let m1 = argvals_s.len();
94    let m2 = argvals_t.len();
95
96    let mut weights = vec![0.0; m1 * m2];
97    for i in 0..m1 {
98        for j in 0..m2 {
99            weights[i * m2 + j] = weights_s[i] * weights_t[j];
100        }
101    }
102    weights
103}
104
105#[cfg(test)]
106mod tests {
107    use super::*;
108
109    #[test]
110    fn test_simpsons_weights_uniform() {
111        let argvals = vec![0.0, 0.25, 0.5, 0.75, 1.0];
112        let weights = simpsons_weights(&argvals);
113        let sum: f64 = weights.iter().sum();
114        assert!((sum - 1.0).abs() < NUMERICAL_EPS);
115    }
116
117    #[test]
118    fn test_simpsons_weights_2d() {
119        let argvals_s = vec![0.0, 0.5, 1.0];
120        let argvals_t = vec![0.0, 0.5, 1.0];
121        let weights = simpsons_weights_2d(&argvals_s, &argvals_t);
122        let sum: f64 = weights.iter().sum();
123        assert!((sum - 1.0).abs() < NUMERICAL_EPS);
124    }
125
126    #[test]
127    fn test_extract_curves() {
128        // Column-major data: 2 observations, 3 points
129        // obs 0: [1, 2, 3], obs 1: [4, 5, 6]
130        let data = vec![1.0, 4.0, 2.0, 5.0, 3.0, 6.0];
131        let mat = crate::matrix::FdMatrix::from_column_major(data, 2, 3).unwrap();
132        let curves = extract_curves(&mat);
133        assert_eq!(curves.len(), 2);
134        assert_eq!(curves[0], vec![1.0, 2.0, 3.0]);
135        assert_eq!(curves[1], vec![4.0, 5.0, 6.0]);
136    }
137
138    #[test]
139    fn test_l2_distance_identical() {
140        let curve = vec![1.0, 2.0, 3.0];
141        let weights = vec![0.25, 0.5, 0.25];
142        let dist = l2_distance(&curve, &curve, &weights);
143        assert!(dist.abs() < NUMERICAL_EPS);
144    }
145
146    #[test]
147    fn test_l2_distance_different() {
148        let curve1 = vec![0.0, 0.0, 0.0];
149        let curve2 = vec![1.0, 1.0, 1.0];
150        let weights = vec![0.25, 0.5, 0.25]; // sum = 1
151        let dist = l2_distance(&curve1, &curve2, &weights);
152        // dist^2 = 0.25*1 + 0.5*1 + 0.25*1 = 1.0, so dist = 1.0
153        assert!((dist - 1.0).abs() < NUMERICAL_EPS);
154    }
155}
fdars_core/helpers.rs

fdars_core/
helpers.rs
