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

fdars_core/
helpers.rs
