1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
//! Tricube-weighted moving average.
//!
//! Pure-Rust port of limma's `tricubeMovingAverage.R`
//! ([`tricube_moving_average`]): a centred moving-average smoother with tricube
//! weights, zero-padded ends, and boundary values rescaled so the padding does
//! not bias them.
/// `tricubeMovingAverage(x, span, power)`. `span` is the window width as a
/// fraction of `length(x)` (clamped to `(0, 1]`); `power` is the tricube
/// exponent (limma's default is 3, R `span = 0.5`). A non-positive `span` or a
/// window narrower than 3 returns `x` unchanged.
pub fn tricube_moving_average(x: &[f64], span: f64, power: f64) -> Vec<f64> {
let n = x.len();
let mut span = span;
if span > 1.0 {
span = 1.0;
}
if span <= 0.0 {
return x.to_vec();
}
let power = power.max(0.0);
// Window width: round span*n down to the nearest odd number, capped at n.
let mut hwidth = (span * n as f64 / 2.0).floor() as i64;
let mut width = 2 * hwidth + 1;
if width > n as i64 {
width -= 2;
hwidth -= 1;
}
if hwidth <= 0 {
return x.to_vec();
}
let width = width as usize;
let hwidth = hwidth as usize;
// Tricube weights on the centred grid, normalised to sum 1 (symmetric).
let wf = width as f64;
let mut weights = vec![0.0_f64; width];
let mut wsum = 0.0;
for (k, wk) in weights.iter_mut().enumerate() {
let u = (-1.0 + k as f64 * 2.0 / (wf - 1.0)) * wf / (wf + 1.0);
*wk = (1.0 - u.abs().powi(3)).powf(power);
wsum += *wk;
}
for wk in weights.iter_mut() {
*wk /= wsum;
}
// Zero-padded centred convolution: out[i] = sum_m weights[m] * P[i+m], where
// P is x padded with hwidth zeros on each side (so P[i+m] = x[i+m-hwidth]
// when in range, else 0). Symmetric weights make the reversal a no-op.
let mut out = vec![0.0_f64; n];
for (i, o) in out.iter_mut().enumerate() {
let mut acc = 0.0;
for (m, &wm) in weights.iter().enumerate() {
let t = i + m;
if t >= hwidth && t < hwidth + n {
acc += wm * x[t - hwidth];
}
}
*o = acc;
}
// Rescale the boundary values by the weight actually covering real data.
let mut cw = vec![0.0_f64; width];
let mut c = 0.0;
for (j, cj) in cw.iter_mut().enumerate() {
c += weights[j];
*cj = c;
}
for a in 0..hwidth {
out[a] /= cw[hwidth + a];
out[n - hwidth + a] /= cw[2 * hwidth - 1 - a];
}
out
}
#[cfg(test)]
mod tests {
use super::*;
fn close(a: &[f64], b: &[f64], tol: f64) -> bool {
a.len() == b.len()
&& a.iter()
.zip(b)
.all(|(&x, &y)| (x - y).abs() <= tol + tol * y.abs())
}
#[test]
fn ramp_span_half_matches_r() {
let x: Vec<f64> = (1..=10).map(|v| v as f64).collect();
// Reference: tricubeMovingAverage(1:10, span=0.5) in limma 3.68.3.
let want = [
1.50604198932743,
2.05598297238272,
3.0,
4.0,
5.0,
6.0,
7.0,
8.0,
8.94401702761728,
9.49395801067257,
];
assert!(close(&tricube_moving_average(&x, 0.5, 3.0), &want, 1e-9));
}
#[test]
fn noisy_span_matches_r() {
let x = [
2.1, -0.4, 1.7, 3.3, 0.9, -1.2, 2.8, 4.1, 0.2, -0.7, 1.5, 3.9,
];
// Reference: tricubeMovingAverage(x, span=0.3) in limma 3.68.3.
let want = [
1.69516075921444,
0.24108900811666,
1.63031641216123,
2.74253129728986,
0.941810152703261,
-0.349860228367038,
2.42370862567066,
3.37529068647682,
0.618101527032604,
-0.267961755399642,
1.52787343513551,
3.51135432884586,
];
assert!(close(&tricube_moving_average(&x, 0.3, 3.0), &want, 1e-9));
}
#[test]
fn span_zero_is_identity() {
let x = [1.0, 2.0, 3.0];
assert_eq!(tricube_moving_average(&x, 0.0, 3.0), x.to_vec());
}
}