fdars_core/detrend/
stl.rs

1use crate::iter_maybe_parallel;
2use crate::matrix::FdMatrix;
3#[cfg(feature = "parallel")]
4use rayon::iter::ParallelIterator;
5
6// ============================================================================
7// STL Decomposition (Cleveland et al., 1990)
8// ============================================================================
9
10/// Result of STL decomposition including robustness weights.
11#[derive(Debug, Clone)]
12#[non_exhaustive]
13pub struct StlResult {
14    /// Trend component (n x m)
15    pub trend: FdMatrix,
16    /// Seasonal component (n x m)
17    pub seasonal: FdMatrix,
18    /// Remainder/residual component (n x m)
19    pub remainder: FdMatrix,
20    /// Robustness weights per point (n x m)
21    pub weights: FdMatrix,
22    /// Period used for decomposition
23    pub period: usize,
24    /// Seasonal smoothing window
25    pub s_window: usize,
26    /// Trend smoothing window
27    pub t_window: usize,
28    /// Number of inner loop iterations performed
29    pub inner_iterations: usize,
30    /// Number of outer loop iterations performed
31    pub outer_iterations: usize,
32}
33
34/// Configuration for STL decomposition.
35///
36/// Collects all tuning parameters for [`stl_decompose_with_config`], with sensible
37/// defaults obtained via [`StlConfig::default()`].
38///
39/// # Example
40/// ```no_run
41/// use fdars_core::detrend::stl::StlConfig;
42///
43/// let config = StlConfig {
44///     robust: true,
45///     s_window: Some(13),
46///     ..StlConfig::default()
47/// };
48/// ```
49#[derive(Debug, Clone, Default)]
50pub struct StlConfig {
51    /// Seasonal smoothing window (default: `None` for auto = 7).
52    pub s_window: Option<usize>,
53    /// Trend smoothing window (default: `None` for auto).
54    pub t_window: Option<usize>,
55    /// Low-pass filter window (default: `None` for auto = period).
56    pub l_window: Option<usize>,
57    /// Whether to use robust (bisquare) weights (default: false).
58    pub robust: bool,
59    /// Number of inner loop iterations (default: `None` for auto = 2).
60    pub inner_iterations: Option<usize>,
61    /// Number of outer loop iterations (default: `None` for auto = 1 or 15 if robust).
62    pub outer_iterations: Option<usize>,
63}
64
65/// STL Decomposition using a [`StlConfig`] struct.
66///
67/// This is the config-based alternative to [`stl_decompose`]. It takes data
68/// and period directly, and reads all tuning parameters from the config.
69///
70/// # Arguments
71/// * `data` — Functional data matrix (n x m)
72/// * `period` — Seasonal period length
73/// * `config` — Tuning parameters
74pub fn stl_decompose_with_config(data: &FdMatrix, period: usize, config: &StlConfig) -> StlResult {
75    stl_decompose(
76        data,
77        period,
78        config.s_window,
79        config.t_window,
80        config.l_window,
81        config.robust,
82        config.inner_iterations,
83        config.outer_iterations,
84    )
85}
86
87/// STL Decomposition: Seasonal and Trend decomposition using LOESS.
88///
89/// # Arguments
90/// * `data` - Functional data matrix (n x m)
91/// * `period` - Seasonal period length
92/// * `s_window` - Seasonal smoothing window (None for auto)
93/// * `t_window` - Trend smoothing window (None for auto)
94/// * `l_window` - Low-pass filter window (None for auto)
95/// * `robust` - Whether to use robust weights
96/// * `inner_iterations` - Number of inner loop iterations (None for auto)
97/// * `outer_iterations` - Number of outer loop iterations (None for auto)
98///
99/// # Examples
100///
101/// ```
102/// use fdars_core::matrix::FdMatrix;
103/// use fdars_core::detrend::stl::stl_decompose;
104///
105/// let n = 3;
106/// let m = 40; // must be >= 2 * period
107/// let data = FdMatrix::from_column_major(
108///     (0..n * m).map(|i| {
109///         let t = (i % m) as f64;
110///         (t * std::f64::consts::PI / 5.0).sin() + t * 0.01
111///     }).collect(),
112///     n, m,
113/// ).unwrap();
114/// let result = stl_decompose(&data, 10, None, None, None, false, None, None);
115/// assert_eq!(result.trend.shape(), (n, m));
116/// assert_eq!(result.seasonal.shape(), (n, m));
117/// assert_eq!(result.remainder.shape(), (n, m));
118/// ```
119pub fn stl_decompose(
120    data: &FdMatrix,
121    period: usize,
122    s_window: Option<usize>,
123    t_window: Option<usize>,
124    l_window: Option<usize>,
125    robust: bool,
126    inner_iterations: Option<usize>,
127    outer_iterations: Option<usize>,
128) -> StlResult {
129    let (n, m) = data.shape();
130    if n == 0 || m < 2 * period || period < 2 {
131        return StlResult {
132            trend: FdMatrix::zeros(n, m),
133            seasonal: FdMatrix::zeros(n, m),
134            remainder: FdMatrix::from_slice(data.as_slice(), n, m)
135                .unwrap_or_else(|_| FdMatrix::zeros(n, m)),
136            weights: FdMatrix::from_column_major(vec![1.0; n * m], n, m)
137                .unwrap_or_else(|_| FdMatrix::zeros(n, m)),
138            period,
139            s_window: 0,
140            t_window: 0,
141            inner_iterations: 0,
142            outer_iterations: 0,
143        };
144    }
145    let s_win = s_window.unwrap_or(7).max(3) | 1;
146    let t_win = t_window.unwrap_or_else(|| {
147        let ratio = 1.5 * period as f64 / (1.0 - 1.5 / s_win as f64);
148        let val = ratio.ceil() as usize;
149        val.max(3) | 1
150    });
151    let l_win = l_window.unwrap_or(period) | 1;
152    let n_inner = inner_iterations.unwrap_or(2);
153    let n_outer = outer_iterations.unwrap_or(if robust { 15 } else { 1 });
154    let results: Vec<(Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>)> = iter_maybe_parallel!(0..n)
155        .map(|i| {
156            let curve: Vec<f64> = (0..m).map(|j| data[(i, j)]).collect();
157            stl_single_series(
158                &curve, period, s_win, t_win, l_win, robust, n_inner, n_outer,
159            )
160        })
161        .collect();
162    let mut trend = FdMatrix::zeros(n, m);
163    let mut seasonal = FdMatrix::zeros(n, m);
164    let mut remainder = FdMatrix::zeros(n, m);
165    let mut weights = FdMatrix::from_column_major(vec![1.0; n * m], n, m)
166        .expect("dimension invariant: data.len() == n * m");
167    for (i, (t, s, r, w)) in results.into_iter().enumerate() {
168        for j in 0..m {
169            trend[(i, j)] = t[j];
170            seasonal[(i, j)] = s[j];
171            remainder[(i, j)] = r[j];
172            weights[(i, j)] = w[j];
173        }
174    }
175    StlResult {
176        trend,
177        seasonal,
178        remainder,
179        weights,
180        period,
181        s_window: s_win,
182        t_window: t_win,
183        inner_iterations: n_inner,
184        outer_iterations: n_outer,
185    }
186}
187
188fn stl_single_series(
189    data: &[f64],
190    period: usize,
191    s_window: usize,
192    t_window: usize,
193    l_window: usize,
194    robust: bool,
195    n_inner: usize,
196    n_outer: usize,
197) -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
198    let m = data.len();
199    let mut trend = vec![0.0; m];
200    let mut seasonal = vec![0.0; m];
201    let mut weights = vec![1.0; m];
202    for outer in 0..n_outer {
203        for _inner in 0..n_inner {
204            let detrended: Vec<f64> = data
205                .iter()
206                .zip(trend.iter())
207                .map(|(&y, &t)| y - t)
208                .collect();
209            let cycle_smoothed = smooth_cycle_subseries(&detrended, period, s_window, &weights);
210            let low_pass = stl_lowpass_filter(&cycle_smoothed, period, l_window);
211            seasonal = cycle_smoothed
212                .iter()
213                .zip(low_pass.iter())
214                .map(|(&c, &l)| c - l)
215                .collect();
216            let deseasonalized: Vec<f64> = data
217                .iter()
218                .zip(seasonal.iter())
219                .map(|(&y, &s)| y - s)
220                .collect();
221            trend = weighted_loess(&deseasonalized, t_window, &weights);
222        }
223        if robust && outer < n_outer - 1 {
224            let remainder: Vec<f64> = data
225                .iter()
226                .zip(trend.iter())
227                .zip(seasonal.iter())
228                .map(|((&y, &t), &s)| y - t - s)
229                .collect();
230            weights = compute_robustness_weights(&remainder);
231        }
232    }
233    let remainder: Vec<f64> = data
234        .iter()
235        .zip(trend.iter())
236        .zip(seasonal.iter())
237        .map(|((&y, &t), &s)| y - t - s)
238        .collect();
239    (trend, seasonal, remainder, weights)
240}
241
242fn smooth_cycle_subseries(
243    data: &[f64],
244    period: usize,
245    s_window: usize,
246    weights: &[f64],
247) -> Vec<f64> {
248    let m = data.len();
249    let n_cycles = m.div_ceil(period);
250    let mut result = vec![0.0; m];
251    for pos in 0..period {
252        let mut subseries_idx: Vec<usize> = Vec::new();
253        let mut subseries_vals: Vec<f64> = Vec::new();
254        let mut subseries_weights: Vec<f64> = Vec::new();
255        for cycle in 0..n_cycles {
256            let idx = cycle * period + pos;
257            if idx < m {
258                subseries_idx.push(idx);
259                subseries_vals.push(data[idx]);
260                subseries_weights.push(weights[idx]);
261            }
262        }
263        if subseries_vals.is_empty() {
264            continue;
265        }
266        let smoothed = weighted_loess(&subseries_vals, s_window, &subseries_weights);
267        for (i, &idx) in subseries_idx.iter().enumerate() {
268            result[idx] = smoothed[i];
269        }
270    }
271    result
272}
273
274fn stl_lowpass_filter(data: &[f64], period: usize, _l_window: usize) -> Vec<f64> {
275    let ma1 = moving_average(data, period);
276    let ma2 = moving_average(&ma1, period);
277    moving_average(&ma2, 3)
278}
279
280fn moving_average(data: &[f64], window: usize) -> Vec<f64> {
281    let m = data.len();
282    if m == 0 || window == 0 {
283        return data.to_vec();
284    }
285    let half = window / 2;
286    let mut result = vec![0.0; m];
287    for i in 0..m {
288        let start = i.saturating_sub(half);
289        let end = (i + half + 1).min(m);
290        let sum: f64 = data[start..end].iter().sum();
291        let count = (end - start) as f64;
292        result[i] = sum / count;
293    }
294    result
295}
296
297fn weighted_loess(data: &[f64], window: usize, weights: &[f64]) -> Vec<f64> {
298    let m = data.len();
299    if m == 0 {
300        return vec![];
301    }
302    let half = window / 2;
303    let mut result = vec![0.0; m];
304    for i in 0..m {
305        let start = i.saturating_sub(half);
306        let end = (i + half + 1).min(m);
307        let mut sum_w = 0.0;
308        let mut sum_wx = 0.0;
309        let mut sum_wy = 0.0;
310        let mut sum_wxx = 0.0;
311        let mut sum_wxy = 0.0;
312        for j in start..end {
313            let dist = (j as f64 - i as f64).abs() / (half.max(1) as f64);
314            let tricube = if dist < 1.0 {
315                (1.0 - dist.powi(3)).powi(3)
316            } else {
317                0.0
318            };
319            let w = tricube * weights[j];
320            let x = j as f64;
321            let y = data[j];
322            sum_w += w;
323            sum_wx += w * x;
324            sum_wy += w * y;
325            sum_wxx += w * x * x;
326            sum_wxy += w * x * y;
327        }
328        if sum_w > 1e-10 {
329            let denom = sum_w * sum_wxx - sum_wx * sum_wx;
330            if denom.abs() > 1e-10 {
331                let intercept = (sum_wxx * sum_wy - sum_wx * sum_wxy) / denom;
332                let slope = (sum_w * sum_wxy - sum_wx * sum_wy) / denom;
333                result[i] = intercept + slope * i as f64;
334            } else {
335                result[i] = sum_wy / sum_w;
336            }
337        } else {
338            result[i] = data[i];
339        }
340    }
341    result
342}
343
344fn compute_robustness_weights(residuals: &[f64]) -> Vec<f64> {
345    let m = residuals.len();
346    if m == 0 {
347        return vec![];
348    }
349    let mut abs_residuals: Vec<f64> = residuals.iter().map(|&r| r.abs()).collect();
350    crate::helpers::sort_nan_safe(&mut abs_residuals);
351    let median_idx = m / 2;
352    let mad = if m % 2 == 0 {
353        (abs_residuals[median_idx - 1] + abs_residuals[median_idx]) / 2.0
354    } else {
355        abs_residuals[median_idx]
356    };
357    let h = 6.0 * mad.max(1e-10);
358    residuals
359        .iter()
360        .map(|&r| {
361            let u = r.abs() / h;
362            if u < 1.0 {
363                (1.0 - u * u).powi(2)
364            } else {
365                0.0
366            }
367        })
368        .collect()
369}
370
371/// Wrapper function for functional data STL decomposition.
372pub fn stl_fdata(
373    data: &FdMatrix,
374    _argvals: &[f64],
375    period: usize,
376    s_window: Option<usize>,
377    t_window: Option<usize>,
378    robust: bool,
379) -> StlResult {
380    stl_decompose(data, period, s_window, t_window, None, robust, None, None)
381}
fdars_core/detrend/stl.rs

fdars_core/detrend/
stl.rs
