anofox-forecast 0.5.0

//! Differencing utilities for ARIMA models.

/// Apply differencing to a time series.
///
/// # Arguments
/// * `series` - The input series
/// * `d` - Differencing order (number of times to difference)
///
/// # Returns
/// The differenced series.
pub fn difference(series: &[f64], d: usize) -> Vec<f64> {
    if d == 0 || series.is_empty() {
        return series.to_vec();
    }

    let mut result = series.to_vec();
    for _ in 0..d {
        if result.len() <= 1 {
            break;
        }
        result = result.windows(2).map(|w| w[1] - w[0]).collect();
    }
    result
}

/// Apply seasonal differencing to a time series.
///
/// # Arguments
/// * `series` - The input series
/// * `d` - Seasonal differencing order
/// * `period` - Seasonal period
///
/// # Returns
/// The seasonally differenced series.
pub fn seasonal_difference(series: &[f64], d: usize, period: usize) -> Vec<f64> {
    if d == 0 || period == 0 || series.len() <= period {
        return series.to_vec();
    }

    let mut result = series.to_vec();
    for _ in 0..d {
        if result.len() <= period {
            break;
        }
        result = result
            .iter()
            .skip(period)
            .zip(result.iter())
            .map(|(curr, prev)| curr - prev)
            .collect();
    }
    result
}

/// Integrate (reverse differencing) a differenced series.
///
/// # Arguments
/// * `differenced` - The differenced series
/// * `original` - The original series (needed for initial values)
/// * `d` - Differencing order used
///
/// # Returns
/// The integrated series.
pub fn integrate(differenced: &[f64], original: &[f64], d: usize) -> Vec<f64> {
    if d == 0 || differenced.is_empty() {
        return differenced.to_vec();
    }

    let mut result = differenced.to_vec();

    // We need to reverse the differencing d times
    for level in (0..d).rev() {
        // Get the initial value at this differencing level
        let init_value = if level == 0 {
            *original.last().unwrap_or(&0.0)
        } else {
            // For higher levels, we need the last value of the intermediate difference
            let intermediate = difference(original, level);
            *intermediate.last().unwrap_or(&0.0)
        };

        // Cumulative sum starting from the initial value
        let mut integrated = Vec::with_capacity(result.len());
        let mut cumsum = init_value;
        for &diff in &result {
            cumsum += diff;
            integrated.push(cumsum);
        }
        result = integrated;
    }

    result
}

/// Integrate (reverse seasonal differencing) a seasonally differenced series.
///
/// # Arguments
/// * `differenced` - The seasonally differenced forecast values
/// * `original` - The original series (needed for the last `period * d` values)
/// * `d` - Seasonal differencing order used
/// * `period` - Seasonal period used
///
/// # Returns
/// The integrated series on the original scale.
pub fn seasonal_integrate(
    differenced: &[f64],
    original: &[f64],
    d: usize,
    period: usize,
) -> Vec<f64> {
    if d == 0 || period == 0 || differenced.is_empty() {
        return differenced.to_vec();
    }

    let mut result = differenced.to_vec();

    for level in (0..d).rev() {
        // Get the reference series at this differencing level
        let reference = if level == 0 {
            original.to_vec()
        } else {
            seasonal_difference(original, level, period)
        };

        // Each forecast value: y[t] = diff[t] + y[t - period]
        // We need the last `period` values from the reference to seed integration
        let mut integrated = Vec::with_capacity(result.len());
        let ref_len = reference.len();

        for (i, &diff_val) in result.iter().enumerate() {
            let prev = if i < period {
                // Seed from the tail of the reference series
                let ref_idx = ref_len.wrapping_sub(period).wrapping_add(i);
                if ref_idx < ref_len {
                    reference[ref_idx]
                } else {
                    0.0
                }
            } else {
                integrated[i - period]
            };
            integrated.push(diff_val + prev);
        }
        result = integrated;
    }

    result
}

/// Check if a series needs differencing using a simple variance ratio test.
///
/// # Arguments
/// * `series` - The input series
///
/// # Returns
/// Suggested differencing order (0, 1, or 2).
pub fn suggest_differencing(series: &[f64]) -> usize {
    if series.len() < 3 {
        return 0;
    }

    let var_0 = variance(series);
    let diff_1 = difference(series, 1);

    if !variance_decreased_significantly(var_0, &diff_1) {
        return 0;
    }

    // Check if second difference helps more
    let var_1 = variance(&diff_1);
    let diff_2 = difference(&diff_1, 1);
    if diff_2.len() >= 2 {
        let var_2 = variance(&diff_2);
        if var_2 / var_1 < 0.9 && var_2 < var_0 {
            return 2;
        }
    }

    1
}

/// Check if differencing significantly reduces variance.
#[inline]
fn variance_decreased_significantly(var_original: f64, differenced: &[f64]) -> bool {
    if differenced.len() < 2 || var_original <= 0.0 {
        return false;
    }
    variance(differenced) / var_original < 0.9
}

/// Calculate variance of a series.
fn variance(series: &[f64]) -> f64 {
    if series.len() < 2 {
        return 0.0;
    }
    crate::simd::variance_sample(series)
}

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

    #[test]
    fn difference_order_0() {
        let series = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let result = difference(&series, 0);
        assert_eq!(result, series);
    }

    #[test]
    fn difference_order_1() {
        let series = vec![1.0, 3.0, 6.0, 10.0, 15.0];
        let result = difference(&series, 1);
        assert_eq!(result, vec![2.0, 3.0, 4.0, 5.0]);
    }

    #[test]
    fn difference_order_2() {
        let series = vec![1.0, 3.0, 6.0, 10.0, 15.0];
        let result = difference(&series, 2);
        // First diff: [2, 3, 4, 5]
        // Second diff: [1, 1, 1]
        assert_eq!(result, vec![1.0, 1.0, 1.0]);
    }

    #[test]
    fn difference_constant_series() {
        let series = vec![5.0, 5.0, 5.0, 5.0];
        let result = difference(&series, 1);
        assert_eq!(result, vec![0.0, 0.0, 0.0]);
    }

    #[test]
    fn difference_empty() {
        let series: Vec<f64> = vec![];
        let result = difference(&series, 1);
        assert!(result.is_empty());
    }

    #[test]
    fn seasonal_difference_basic() {
        // Quarterly data: Q1 values increase by 10 each year
        let series = vec![
            100.0, 120.0, 80.0, 90.0, // Year 1
            110.0, 130.0, 90.0, 100.0, // Year 2
        ];
        let result = seasonal_difference(&series, 1, 4);
        // Each value minus same quarter previous year
        assert_eq!(result, vec![10.0, 10.0, 10.0, 10.0]);
    }

    #[test]
    fn seasonal_difference_no_change() {
        let series = vec![1.0, 2.0, 3.0, 1.0, 2.0, 3.0];
        let result = seasonal_difference(&series, 1, 3);
        assert_eq!(result, vec![0.0, 0.0, 0.0]);
    }

    #[test]
    fn seasonal_difference_order_0() {
        let series = vec![1.0, 2.0, 3.0, 4.0];
        let result = seasonal_difference(&series, 0, 2);
        assert_eq!(result, series);
    }

    #[test]
    fn integrate_reverses_difference() {
        let original = vec![10.0, 12.0, 15.0, 19.0, 24.0];
        let _differenced = difference(&original, 1);
        let forecast_diff = vec![6.0, 7.0]; // Forecasted differences
        let integrated = integrate(&forecast_diff, &original, 1);

        // Should continue from last value: 24 + 6 = 30, 30 + 7 = 37
        assert_relative_eq!(integrated[0], 30.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[1], 37.0, epsilon = 1e-10);
    }

    #[test]
    fn integrate_order_2() {
        let original = vec![1.0, 3.0, 6.0, 10.0, 15.0];
        let _differenced = difference(&original, 2);
        // differenced = [1, 1, 1]

        let forecast_diff2 = vec![1.0, 1.0]; // Continue the pattern
        let integrated = integrate(&forecast_diff2, &original, 2);

        // The integration should produce reasonable continuation
        assert!(integrated.len() == 2);
    }

    #[test]
    fn suggest_differencing_stationary() {
        // White noise-like stationary series
        let series = vec![1.0, 0.5, 1.2, 0.8, 1.1, 0.9, 1.0, 1.1];
        let d = suggest_differencing(&series);
        assert_eq!(d, 0);
    }

    #[test]
    fn suggest_differencing_trend() {
        // Clear upward trend
        let series: Vec<f64> = (0..20).map(|i| 10.0 + 2.0 * i as f64).collect();
        let d = suggest_differencing(&series);
        assert!(d >= 1);
    }

    #[test]
    fn suggest_differencing_quadratic() {
        // Quadratic trend
        let series: Vec<f64> = (0..20).map(|i| (i * i) as f64).collect();
        let d = suggest_differencing(&series);
        assert!(d >= 1);
    }

    #[test]
    fn seasonal_integrate_reverses_seasonal_difference() {
        // Quarterly data with year-over-year growth of 10
        let original = vec![
            100.0, 120.0, 80.0, 90.0, // Year 1
            110.0, 130.0, 90.0, 100.0, // Year 2
        ];
        let diffed = seasonal_difference(&original, 1, 4);
        assert_eq!(diffed, vec![10.0, 10.0, 10.0, 10.0]);

        // Forecast: continue the same year-over-year growth
        let forecast_diff = vec![10.0, 10.0, 10.0, 10.0];
        let integrated = seasonal_integrate(&forecast_diff, &original, 1, 4);

        // Year 3 should be year 2 + 10 for each quarter
        assert_relative_eq!(integrated[0], 120.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[1], 140.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[2], 100.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[3], 110.0, epsilon = 1e-10);
    }

    #[test]
    fn seasonal_integrate_order_0_is_identity() {
        let forecast = vec![1.0, 2.0, 3.0];
        let original = vec![10.0, 20.0, 30.0, 40.0];
        let result = seasonal_integrate(&forecast, &original, 0, 4);
        assert_eq!(result, forecast);
    }

    #[test]
    fn seasonal_integrate_empty() {
        let result = seasonal_integrate(&[], &[1.0, 2.0], 1, 2);
        assert!(result.is_empty());
    }

    // ── Extended seasonal_integrate correctness tests ────────────

    #[test]
    fn seasonal_integrate_reverses_multi_period_forecast() {
        // Weekly data (period=7) over 3 weeks with growth
        let period = 7;
        let weekly_base = [10.0, 20.0, 15.0, 25.0, 30.0, 12.0, 8.0];
        let growth = 5.0;
        let original: Vec<f64> = (0..21)
            .map(|i| weekly_base[i % period] + growth * (i / period) as f64)
            .collect();

        let diffed = seasonal_difference(&original, 1, period);
        // Each diff should be == growth (5.0) for the last two weeks
        assert_eq!(diffed.len(), 14);
        for &d in &diffed {
            assert_relative_eq!(d, growth, epsilon = 1e-10);
        }

        // Forecast 14 steps (2 more weeks) of same growth
        let forecast_diff: Vec<f64> = vec![growth; 14];
        let integrated = seasonal_integrate(&forecast_diff, &original, 1, period);

        // First 7 values should be original week 3 + 5 (= week 4)
        for i in 0..7 {
            let expected = weekly_base[i] + growth * 3.0;
            assert_relative_eq!(integrated[i], expected, epsilon = 1e-10);
        }
        // Next 7 values should be week 5 = week 4 + 5
        for i in 0..7 {
            let expected = weekly_base[i] + growth * 4.0;
            assert_relative_eq!(integrated[7 + i], expected, epsilon = 1e-10);
        }
    }

    #[test]
    fn seasonal_integrate_horizon_less_than_period() {
        // Period=7 but forecast only 3 steps — all seeds come from original
        let period = 7;
        let original = vec![
            10.0, 20.0, 15.0, 25.0, 30.0, 12.0, 8.0, // Week 1
            15.0, 25.0, 20.0, 30.0, 35.0, 17.0, 13.0, // Week 2
        ];

        // Forecast differences for 3 steps
        let forecast_diff = vec![5.0, 5.0, 5.0];
        let integrated = seasonal_integrate(&forecast_diff, &original, 1, period);

        // When horizon < period, all seeds come from the tail of the original.
        // integrated[0] = diff[0] + original[14 - 7 + 0] = 5 + original[7] = 5 + 15 = 20
        // integrated[1] = diff[1] + original[14 - 7 + 1] = 5 + original[8] = 5 + 25 = 30
        // integrated[2] = diff[2] + original[14 - 7 + 2] = 5 + original[9] = 5 + 20 = 25
        assert_eq!(integrated.len(), 3);
        assert_relative_eq!(integrated[0], 20.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[1], 30.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[2], 25.0, epsilon = 1e-10);
    }

    #[test]
    fn seasonal_integrate_horizon_greater_than_period() {
        // Period=4, forecast 6 steps — last 2 use previously integrated values
        let period = 4;
        let original = vec![
            100.0, 200.0, 150.0, 250.0, // Q1
            110.0, 210.0, 160.0, 260.0, // Q2
        ];

        // Growth of 10 for each seasonal position
        let forecast_diff = vec![10.0; 6];
        let integrated = seasonal_integrate(&forecast_diff, &original, 1, period);

        assert_eq!(integrated.len(), 6);
        // First 4: seeded from tail of original (last 4 values = Q2)
        // integrated[0] = 10 + original[8-4+0] = 10 + 110 = 120
        // integrated[1] = 10 + original[8-4+1] = 10 + 210 = 220
        // integrated[2] = 10 + original[8-4+2] = 10 + 160 = 170
        // integrated[3] = 10 + original[8-4+3] = 10 + 260 = 270
        assert_relative_eq!(integrated[0], 120.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[1], 220.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[2], 170.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[3], 270.0, epsilon = 1e-10);

        // Next 2: use previously integrated values
        // integrated[4] = 10 + integrated[4-4] = 10 + 120 = 130
        // integrated[5] = 10 + integrated[5-4] = 10 + 220 = 230
        assert_relative_eq!(integrated[4], 130.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[5], 230.0, epsilon = 1e-10);
    }

    #[test]
    fn seasonal_integrate_d2_reversal() {
        // Seasonal differencing with d=2, period=4
        let original = vec![
            100.0, 120.0, 80.0, 90.0, // Year 1
            110.0, 130.0, 90.0, 100.0, // Year 2
            125.0, 145.0, 105.0, 115.0, // Year 3
        ];

        let diffed = seasonal_difference(&original, 2, 4);
        // d=2 means: apply seasonal diff twice
        // First diff (d=1): [10, 10, 10, 10, 15, 15, 15, 15]
        // Second diff (d=1 on that): [5, 5, 5, 5]
        assert_eq!(diffed.len(), 4);
        for &d in &diffed {
            assert_relative_eq!(d, 5.0, epsilon = 1e-10);
        }

        // Forecast: same acceleration
        let forecast_diff = vec![5.0; 4];
        let integrated = seasonal_integrate(&forecast_diff, &original, 2, 4);

        // Year 4 should continue the accelerating pattern:
        // Year 3 growth was [15, 15, 15, 15], so year 4 growth should be [20, 20, 20, 20]
        assert_eq!(integrated.len(), 4);
        assert_relative_eq!(integrated[0], 145.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[1], 165.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[2], 125.0, epsilon = 1e-10);
        assert_relative_eq!(integrated[3], 135.0, epsilon = 1e-10);
    }

    #[test]
    fn integrate_reverses_difference_d1_multi_step() {
        // Verify integrate properly reverses d=1 for multi-step forecasts
        let original = vec![10.0, 13.0, 17.0, 22.0, 28.0];
        // Differences: [3, 4, 5, 6]
        // Forecast 4 more differences
        let forecast_diff = vec![7.0, 8.0, 9.0, 10.0];
        let integrated = integrate(&forecast_diff, &original, 1);

        // Should cumsum from last original value (28)
        assert_relative_eq!(integrated[0], 35.0, epsilon = 1e-10); // 28 + 7
        assert_relative_eq!(integrated[1], 43.0, epsilon = 1e-10); // 35 + 8
        assert_relative_eq!(integrated[2], 52.0, epsilon = 1e-10); // 43 + 9
        assert_relative_eq!(integrated[3], 62.0, epsilon = 1e-10); // 52 + 10
    }

    #[test]
    fn difference_then_integrate_roundtrip_d1() {
        // Verify that diff then integrate of the same data recovers it
        let original = vec![5.0, 8.0, 12.0, 17.0, 23.0, 30.0];
        let diffed = difference(&original, 1);
        // diffed = [3, 4, 5, 6, 7]

        // If we "forecast" the same differences, integration should reproduce the original
        let recovered = integrate(&diffed, &[original[0]], 1);
        // Starting from 5: 5+3=8, 8+4=12, 12+5=17, 17+6=23, 23+7=30
        for (i, &v) in recovered.iter().enumerate() {
            assert_relative_eq!(v, original[i + 1], epsilon = 1e-10);
        }
    }

    #[test]
    fn seasonal_difference_then_integrate_roundtrip() {
        // Weekly pattern: verify seasonal diff then integrate recovers scale
        let period = 7;
        let original: Vec<f64> = (0..28)
            .map(|i| {
                let base = [10.0, 20.0, 15.0, 25.0, 30.0, 12.0, 8.0];
                base[i % period] + 3.0 * (i / period) as f64
            })
            .collect();

        let diffed = seasonal_difference(&original, 1, period);
        // All diffs should be 3.0
        assert_eq!(diffed.len(), 21);
        for &d in &diffed {
            assert_relative_eq!(d, 3.0, epsilon = 1e-10);
        }

        // Forecast one more week of same growth
        let forecast_diff = vec![3.0; 7];
        let integrated = seasonal_integrate(&forecast_diff, &original, 1, period);

        // Week 5 = week 4 + 3 for each day
        let week4_base = [10.0, 20.0, 15.0, 25.0, 30.0, 12.0, 8.0];
        for i in 0..7 {
            let expected = week4_base[i] + 3.0 * 4.0; // 4th growth increment
            assert_relative_eq!(integrated[i], expected, epsilon = 1e-10);
        }
    }
}