shape-runtime 0.3.2

//! Rolling standard deviation using Welford's algorithm

/// Rolling std (auto-selects Welford algorithm, always optimal)
#[inline]
pub fn rolling_std(data: &[f64], window: usize) -> Vec<f64> {
    rolling_std_welford(data, window) // Welford is always best
}

/// Welford's algorithm for numerically stable rolling standard deviation
///
/// This is significantly faster than the naive two-pass approach and
/// more numerically stable than the sum-of-squares method.
///
/// # Performance
/// - Expected speedup: 10-50x over naive O(n*w) two-pass
/// - Complexity: O(n)
/// - Numerical stability: Superior to sum-of-squares
pub fn rolling_std_welford(data: &[f64], window: usize) -> Vec<f64> {
    let n = data.len();

    if window == 0 || window > n {
        return vec![f64::NAN; n];
    }

    let mut result = vec![f64::NAN; n];

    // Initialize with first window using Welford's algorithm
    let mut mean = 0.0;
    let mut m2 = 0.0; // Sum of squared differences from mean

    for i in 0..window {
        let delta = data[i] - mean;
        mean += delta / (i + 1) as f64;
        let delta2 = data[i] - mean;
        m2 += delta * delta2;
    }

    result[window - 1] = (m2 / window as f64).sqrt();

    // Slide window: remove oldest, add newest
    for i in window..n {
        let old_val = data[i - window];
        let new_val = data[i];

        // Remove old value from running statistics
        let old_mean = mean;
        mean -= (old_val - mean) / window as f64;
        m2 -= (old_val - mean) * (old_val - old_mean);

        // Add new value
        let delta = new_val - mean;
        mean += delta / window as f64;
        let delta2 = new_val - mean;
        m2 += delta * delta2;

        // Ensure non-negative due to floating point errors
        result[i] = (m2.max(0.0) / window as f64).sqrt();
    }

    result
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_rolling_std_welford() {
        let data = vec![2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
        let result = rolling_std_welford(&data, 4);

        // First 3 are NaN
        for i in 0..3 {
            assert!(result[i].is_nan());
        }

        // Check that we get reasonable std values
        assert!(result[3] > 0.0 && result[3] < 2.0); // std of [2,4,4,4]
        assert!(result[7] > 0.0 && result[7] < 3.0); // std of [5,5,7,9]
    }
}