anofox-forecast 0.15.9

//! Seasonal differencing as a composable seasonality component.
//!
//! Exposes ARIMA-style seasonal differencing (`y'_t = y_t - y_{t-period}`) as a
//! standalone [`SeasonalComponent`] that can be composed with other components.
//! Seasonal differencing removes stable periodic patterns of a given period
//! without imposing a parametric model on the seasonal shape.
//!
//! # Example
//!
//! ```
//! use anofox_forecast::seasonality::SeasonalDifference;
//! use anofox_forecast::seasonality::traits::SeasonalComponent;
//!
//! // Monthly data with period 12
//! let values: Vec<f64> = (0..48)
//!     .map(|i| 10.0 + (i % 12) as f64 * 2.0)
//!     .collect();
//!
//! let mut sd = SeasonalDifference::new(12).unwrap();
//! sd.fit_seasonal(&values, 12).unwrap();
//!
//! let differenced = sd.differenced();
//! assert_eq!(differenced.len(), 36); // 48 - 12
//!
//! // Inverse recovers the original series
//! let recovered = sd.inverse(differenced).unwrap();
//! assert_eq!(recovered.len(), 48);
//! ```

use super::traits::SeasonalComponent;
use crate::error::{ForecastError, Result};

/// Seasonal differencing component.
///
/// Removes periodic patterns by computing `y'_t = y_t - y_{t-period}` and
/// provides the inverse transform to reconstruct original-scale values.
#[derive(Debug, Clone)]
pub struct SeasonalDifference {
    period: usize,
    /// The first `period` values from the training data (needed for inversion).
    initial_values: Vec<f64>,
    /// The differenced series.
    differenced: Vec<f64>,
    /// Fitted seasonal values (seasonal_indices repeated over data length).
    fitted: Vec<f64>,
    /// The last cycle of seasonal indices (average of y_t - mean(y) per position).
    seasonal_indices: Vec<f64>,
}

impl SeasonalDifference {
    /// Create a new seasonal differencing component with the given period.
    ///
    /// # Errors
    ///
    /// Returns [`ForecastError::InvalidParameter`] if `period` is zero.
    pub fn new(period: usize) -> Result<Self> {
        if period == 0 {
            return Err(ForecastError::InvalidParameter(
                "seasonal differencing period must be > 0".to_string(),
            ));
        }
        Ok(Self {
            period,
            initial_values: Vec::new(),
            differenced: Vec::new(),
            fitted: Vec::new(),
            seasonal_indices: Vec::new(),
        })
    }

    /// Return the differenced series (`y'_t = y_t - y_{t-period}` for `t >= period`).
    ///
    /// Empty until [`fit_seasonal`](SeasonalComponent::fit_seasonal) is called.
    pub fn differenced(&self) -> &[f64] {
        &self.differenced
    }

    /// Reconstruct original-scale values from a differenced series.
    ///
    /// Uses the stored initial values (first `period` observations from training)
    /// to invert the differencing: `y_t = y'_t + y_{t-period}`.
    ///
    /// The returned vector has length `period + differenced.len()`.
    ///
    /// # Errors
    ///
    /// Returns [`ForecastError::FitRequired`] if the component has not been fitted.
    pub fn inverse(&self, differenced: &[f64]) -> Result<Vec<f64>> {
        if self.initial_values.is_empty() {
            return Err(ForecastError::FitRequired {
                model: Some("SeasonalDifference".to_string()),
            });
        }
        let n = self.period + differenced.len();
        let mut out = Vec::with_capacity(n);
        out.extend_from_slice(&self.initial_values);
        for (i, &d) in differenced.iter().enumerate() {
            out.push(d + out[i]); // y_t = y'_t + y_{t-period}
        }
        Ok(out)
    }

    /// The seasonal period.
    pub fn period(&self) -> usize {
        self.period
    }
}

impl SeasonalComponent for SeasonalDifference {
    fn fit_seasonal(&mut self, values: &[f64], period: usize) -> Result<()> {
        if values.is_empty() {
            return Err(ForecastError::EmptyData);
        }
        if period == 0 {
            return Err(ForecastError::InvalidParameter(
                "seasonal differencing period must be > 0".to_string(),
            ));
        }
        let n = values.len();
        if n <= period {
            return Err(ForecastError::InsufficientData {
                needed: period + 1,
                got: n,
                hint: Some(format!(
                    "need more than {} observations for seasonal differencing with period {}",
                    period, period
                )),
            });
        }

        self.period = period;

        // Store initial values (first `period` observations).
        self.initial_values = values[..period].to_vec();

        // Compute seasonal difference: y'_t = y_t - y_{t-period} for t >= period.
        self.differenced = values[period..]
            .iter()
            .enumerate()
            .map(|(i, &v)| v - values[i])
            .collect();

        // Compute seasonal indices: for each position p in 0..period,
        // average y_t - mean(y) across all t where t % period == p.
        let mean: f64 = values.iter().sum::<f64>() / n as f64;
        let mut sums = vec![0.0; period];
        let mut counts = vec![0usize; period];
        for (t, &v) in values.iter().enumerate() {
            let p = t % period;
            sums[p] += v - mean;
            counts[p] += 1;
        }
        self.seasonal_indices = sums
            .iter()
            .zip(counts.iter())
            .map(|(&s, &c)| if c > 0 { s / c as f64 } else { 0.0 })
            .collect();

        // Fitted seasonal = seasonal_indices repeated over data length.
        self.fitted = (0..n).map(|t| self.seasonal_indices[t % period]).collect();

        Ok(())
    }

    fn fitted_seasonal(&self) -> &[f64] {
        &self.fitted
    }

    fn predict_seasonal(&self, n_ahead: usize) -> Vec<f64> {
        if self.seasonal_indices.is_empty() || self.period == 0 {
            return vec![0.0; n_ahead];
        }
        (0..n_ahead)
            .map(|i| self.seasonal_indices[i % self.period])
            .collect()
    }

    fn seasonal_features(&self) -> Vec<(&str, f64)> {
        if self.fitted.is_empty() {
            return Vec::new();
        }

        let strength = compute_strength(&self.differenced, &self.fitted);
        let variance_reduction = compute_variance_reduction(&self.differenced, &self.fitted);
        let acf1 = compute_acf1(&self.differenced);
        let amplitude = compute_amplitude(&self.seasonal_indices);

        vec![
            ("seasonal_diff_strength", strength),
            ("seasonal_diff_variance_reduction", variance_reduction),
            ("seasonal_diff_acf1", acf1),
            ("seasonal_diff_amplitude", amplitude),
        ]
    }

    fn seasonal_name(&self) -> &str {
        "seasonal_difference"
    }

    fn n_params(&self) -> usize {
        self.period
    }
}

// ── Helper functions ────────────────────────────────────────────────────────

/// Variance of a slice (population variance).
fn variance(data: &[f64]) -> f64 {
    if data.is_empty() {
        return 0.0;
    }
    let n = data.len() as f64;
    let mean = data.iter().sum::<f64>() / n;
    data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n
}

/// Compute seasonal differencing strength: `1 - var(differenced) / var(original_fitted)`.
///
/// `original_fitted` here is the fitted seasonal values (same length as the original series),
/// and we need the original series variance. We reconstruct it from the fitted seasonal
/// and the knowledge that the original = fitted_seasonal + remainder.
///
/// More directly: strength = 1 - var(differenced) / var(original).
/// Since we store the fitted seasonal (same length as original), we can compute
/// var(original) from the initial_values and differenced series via inversion,
/// but the simplest correct approach uses the stored fitted array length to know n,
/// and we compute var over the original via the seasonal indices.
///
/// Actually, the cleanest formulation: we need var(original). But we don't store
/// the original directly. We do have initial_values + differenced which lets us
/// reconstruct it. However, for the standalone functions we get the original directly.
///
/// For the trait impl, we compute from the fitted seasonal array (which has the same
/// length as the original) and use it to reconstruct the original variance.
/// But simpler: we just store enough info. Let's compute it from differenced and fitted.
///
/// The fitted array IS the seasonal component. The original = seasonal + remainder.
/// var(original) = var(seasonal + remainder). We need the actual original variance.
///
/// Since we have initial_values and differenced, we can reconstruct original.
/// But that's wasteful. Instead, let's compute from the differenced series and the
/// fitted seasonal: the original's variance can be computed from the reconstructed series.
///
/// For simplicity and correctness: strength = 1 - var(differenced) / var(reconstructed_original).
fn compute_strength(differenced: &[f64], fitted: &[f64]) -> f64 {
    // We need variance of the original. The fitted array has the same length as the original,
    // and we can use the initial_values + differenced to get original. But we don't have
    // initial_values here. Instead, strength is defined as measuring how much variance the
    // seasonal differencing removed. We use:
    //   strength = 1 - var(differenced) / var(original)
    // where original is represented by the fitted seasonal pattern.
    // Actually: "original" means the input series. The fitted is the seasonal component.
    // For the feature, we use: 1 - var(remainder) / var(data)
    // where remainder = data - fitted_seasonal and data is the original.
    // But we don't store the original directly.
    //
    // A pragmatic approach: compute var(remainder) and var(data) from fitted.
    // remainder = original - fitted_seasonal. We don't have the original here.
    //
    // The cleanest correct measure using what we have:
    //   strength = max(0, 1 - var(differenced) / var_of_original)
    // Since we don't store the original but DO have fitted (the seasonal indices repeated),
    // we can compute: var_of_original ~ var(fitted) + var(differenced) in the ideal case.
    //
    // Instead: use the simple ratio approach with the original reconstructed from fitted.
    // The fitted seasonal has variance = var(seasonal_indices repeated).
    // Actually the simplest correct thing: fitted = seasonal_indices repeated.
    // original variance ≈ fitted variance + residual variance (if uncorrelated).
    //
    // Let's just use: strength = 1 - var(differenced) / (var(fitted) + var(differenced))
    // This gives 0 when var(fitted)=0, 1 when var(differenced)=0.
    // But this isn't the standard definition.
    //
    // Standard: strength = max(0, 1 - var(R_t) / var(Y_t - T_t))
    // where R_t is the remainder and Y_t - T_t is the detrended series.
    //
    // For seasonal differencing without an explicit trend, we use:
    //   strength = max(0, 1 - var(differenced) / (var(fitted) + var(differenced)))
    // when var(fitted) + var(differenced) > 0.

    let var_diff = variance(differenced);
    let var_fitted = variance(fitted);
    let var_total = var_fitted + var_diff;

    if var_total <= 0.0 {
        return 0.0;
    }
    (1.0 - var_diff / var_total).clamp(0.0, 1.0)
}

/// Compute variance reduction ratio: var(original) / var(differenced).
/// A ratio > 1 means differencing reduced variance (seasonality was present).
fn compute_variance_reduction(differenced: &[f64], fitted: &[f64]) -> f64 {
    let var_diff = variance(differenced);
    let var_fitted = variance(fitted);
    let var_original = var_fitted + var_diff;

    if var_diff <= 0.0 {
        return f64::INFINITY;
    }
    var_original / var_diff
}

/// Compute the first autocorrelation (ACF at lag 1) of the differenced series.
///
/// acf(1) = sum((y_t - mean)(y_{t-1} - mean)) / sum((y_t - mean)^2)
fn compute_acf1(data: &[f64]) -> f64 {
    if data.len() < 2 {
        return 0.0;
    }
    let n = data.len() as f64;
    let mean = data.iter().sum::<f64>() / n;

    let denom: f64 = data.iter().map(|&x| (x - mean).powi(2)).sum();
    if denom == 0.0 {
        return 0.0;
    }

    let numer: f64 = data.windows(2).map(|w| (w[1] - mean) * (w[0] - mean)).sum();

    numer / denom
}

/// Compute the amplitude of seasonal indices: max - min.
fn compute_amplitude(indices: &[f64]) -> f64 {
    if indices.is_empty() {
        return 0.0;
    }
    let min = indices.iter().cloned().fold(f64::INFINITY, f64::min);
    let max = indices.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    max - min
}

// ── Standalone feature functions ────────────────────────────────────────────

/// Compute seasonal differencing strength for the given values and period.
///
/// Returns `1 - var(differenced) / (var(seasonal_component) + var(differenced))`,
/// clamped to `[0, 1]`. Values near 1 indicate strong seasonality; near 0 indicate
/// no seasonal pattern at the given period.
///
/// Returns 0.0 if `values.len() <= period` or `period == 0`.
pub fn seasonal_diff_strength(values: &[f64], period: usize) -> f64 {
    if period == 0 || values.len() <= period {
        return 0.0;
    }
    let mut sd = match SeasonalDifference::new(period) {
        Ok(sd) => sd,
        Err(_) => return 0.0,
    };
    if sd.fit_seasonal(values, period).is_err() {
        return 0.0;
    }
    let features = sd.seasonal_features();
    features
        .iter()
        .find(|(name, _)| *name == "seasonal_diff_strength")
        .map(|(_, v)| *v)
        .unwrap_or(0.0)
}

/// Compute seasonal differencing variance reduction for the given values and period.
///
/// Returns `var(original) / var(differenced)`. A ratio greater than 1 means the
/// seasonal differencing reduced variance, indicating seasonality was present.
///
/// Returns 1.0 if `values.len() <= period` or `period == 0`.
pub fn seasonal_diff_variance_reduction(values: &[f64], period: usize) -> f64 {
    if period == 0 || values.len() <= period {
        return 1.0;
    }
    let mut sd = match SeasonalDifference::new(period) {
        Ok(sd) => sd,
        Err(_) => return 1.0,
    };
    if sd.fit_seasonal(values, period).is_err() {
        return 1.0;
    }
    let features = sd.seasonal_features();
    features
        .iter()
        .find(|(name, _)| *name == "seasonal_diff_variance_reduction")
        .map(|(_, v)| *v)
        .unwrap_or(1.0)
}

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

    /// Generate a pure seasonal pattern: each position p has value `p as f64`.
    fn pure_seasonal(period: usize, n_cycles: usize) -> Vec<f64> {
        (0..period * n_cycles)
            .map(|i| (i % period) as f64)
            .collect()
    }

    /// Generate random-ish noise (deterministic, no external crate).
    fn pseudo_noise(n: usize, seed: u64) -> Vec<f64> {
        let mut state = seed;
        (0..n)
            .map(|_| {
                // Simple LCG
                state = state.wrapping_mul(6364136223846793005).wrapping_add(1);
                // Map to [-1, 1]
                (state >> 33) as f64 / (1u64 << 31) as f64 - 1.0
            })
            .collect()
    }

    // ── Perfect seasonal pattern ────────────────────────────────────────

    #[test]
    fn perfect_seasonal_strength_near_one() {
        let period = 7;
        let values = pure_seasonal(period, 10);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let features = sd.seasonal_features();
        let strength = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_strength")
            .unwrap()
            .1;

        // Differenced series should be all zeros for a perfect repeating pattern.
        for &d in sd.differenced() {
            assert_abs_diff_eq!(d, 0.0, epsilon = 1e-12);
        }
        assert_abs_diff_eq!(strength, 1.0, epsilon = 1e-10);
    }

    // ── No seasonal pattern ────────────────────────────────────────────

    #[test]
    fn no_seasonal_strength_near_zero() {
        // Constant series: no seasonal pattern.
        let values = vec![5.0; 100];

        let mut sd = SeasonalDifference::new(7).unwrap();
        sd.fit_seasonal(&values, 7).unwrap();

        let features = sd.seasonal_features();
        let strength = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_strength")
            .unwrap()
            .1;

        assert_abs_diff_eq!(strength, 0.0, epsilon = 1e-10);
    }

    #[test]
    fn noise_has_low_strength() {
        let values = pseudo_noise(200, 42);

        let mut sd = SeasonalDifference::new(7).unwrap();
        sd.fit_seasonal(&values, 7).unwrap();

        let features = sd.seasonal_features();
        let strength = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_strength")
            .unwrap()
            .1;

        // Pure noise should have low seasonal strength (not necessarily zero,
        // but far from 1).
        assert!(
            strength < 0.5,
            "noise strength {} should be < 0.5",
            strength
        );
    }

    // ── Inverse recovers original ──────────────────────────────────────

    #[test]
    fn inverse_recovers_original() {
        let period = 4;
        let values: Vec<f64> = (0..20)
            .map(|i| (i as f64).sin() * 10.0 + i as f64)
            .collect();

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let recovered = sd.inverse(sd.differenced()).unwrap();
        assert_eq!(recovered.len(), values.len());

        for (&orig, &rec) in values.iter().zip(recovered.iter()) {
            assert_abs_diff_eq!(orig, rec, epsilon = 1e-10);
        }
    }

    #[test]
    fn inverse_with_known_values() {
        let values = vec![1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 11.0];
        let period = 4;

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        // differenced = [5-1, 7-2, 9-3, 11-4] = [4, 5, 6, 7]
        let diff = sd.differenced();
        assert_eq!(diff.len(), 4);
        assert_abs_diff_eq!(diff[0], 4.0, epsilon = 1e-12);
        assert_abs_diff_eq!(diff[1], 5.0, epsilon = 1e-12);
        assert_abs_diff_eq!(diff[2], 6.0, epsilon = 1e-12);
        assert_abs_diff_eq!(diff[3], 7.0, epsilon = 1e-12);

        let recovered = sd.inverse(diff).unwrap();
        for (&o, &r) in values.iter().zip(recovered.iter()) {
            assert_abs_diff_eq!(o, r, epsilon = 1e-12);
        }
    }

    // ── Predict continues pattern ──────────────────────────────────────

    #[test]
    fn predict_continues_seasonal_pattern() {
        let period = 4;
        let values = pure_seasonal(period, 5); // [0,1,2,3, 0,1,2,3, ...]

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let prediction = sd.predict_seasonal(8);
        assert_eq!(prediction.len(), 8);

        // Seasonal indices should capture the repeating pattern.
        // For a 0-mean-centred pattern: indices = [-1.5, -0.5, 0.5, 1.5].
        // Prediction repeats these.
        for i in 0..8 {
            assert_abs_diff_eq!(prediction[i], prediction[i % period], epsilon = 1e-10,);
        }
        // The pattern should repeat with period.
        assert_abs_diff_eq!(prediction[0], prediction[4], epsilon = 1e-10);
        assert_abs_diff_eq!(prediction[1], prediction[5], epsilon = 1e-10);
    }

    #[test]
    fn predict_zero_ahead() {
        let values = pure_seasonal(4, 3);
        let mut sd = SeasonalDifference::new(4).unwrap();
        sd.fit_seasonal(&values, 4).unwrap();

        let prediction = sd.predict_seasonal(0);
        assert!(prediction.is_empty());
    }

    // ── Features extraction ────────────────────────────────────────────

    #[test]
    fn features_all_present() {
        let values = pure_seasonal(7, 5);
        let mut sd = SeasonalDifference::new(7).unwrap();
        sd.fit_seasonal(&values, 7).unwrap();

        let features = sd.seasonal_features();
        assert_eq!(features.len(), 4);

        let names: Vec<&str> = features.iter().map(|(n, _)| *n).collect();
        assert!(names.contains(&"seasonal_diff_strength"));
        assert!(names.contains(&"seasonal_diff_variance_reduction"));
        assert!(names.contains(&"seasonal_diff_acf1"));
        assert!(names.contains(&"seasonal_diff_amplitude"));
    }

    #[test]
    fn features_empty_before_fit() {
        let sd = SeasonalDifference::new(7).unwrap();
        let features = sd.seasonal_features();
        assert!(features.is_empty());
    }

    #[test]
    fn variance_reduction_greater_than_one_for_seasonal() {
        let period = 12;
        let values: Vec<f64> = (0..120)
            .map(|i| (i % period) as f64 * 3.0 + 100.0)
            .collect();

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let features = sd.seasonal_features();
        let vr = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_variance_reduction")
            .unwrap()
            .1;

        assert!(
            vr > 1.0,
            "variance reduction {} should be > 1 for seasonal data",
            vr
        );
    }

    #[test]
    fn acf1_near_zero_for_perfect_seasonal() {
        let period = 4;
        let values = pure_seasonal(period, 20);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let features = sd.seasonal_features();
        let acf1 = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_acf1")
            .unwrap()
            .1;

        // Differenced series is all zeros, so acf1 should be 0 (or NaN handled as 0).
        assert_abs_diff_eq!(acf1, 0.0, epsilon = 1e-10);
    }

    #[test]
    fn amplitude_captures_range() {
        let period = 4;
        // seasonal indices for pattern [0,1,2,3] repeated are centred:
        // [-1.5, -0.5, 0.5, 1.5], amplitude = 3.0
        let values = pure_seasonal(period, 5);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let features = sd.seasonal_features();
        let amp = features
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_amplitude")
            .unwrap()
            .1;

        assert_abs_diff_eq!(amp, 3.0, epsilon = 1e-10);
    }

    // ── Standalone functions ───────────────────────────────────────────

    #[test]
    fn standalone_strength_matches_trait() {
        let period = 7;
        let values = pure_seasonal(period, 8);

        let standalone = seasonal_diff_strength(&values, period);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();
        let trait_strength = sd
            .seasonal_features()
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_strength")
            .unwrap()
            .1;

        assert_abs_diff_eq!(standalone, trait_strength, epsilon = 1e-12);
    }

    #[test]
    fn standalone_variance_reduction_matches_trait() {
        let period = 12;
        // Add a small trend so differenced is not all-zero (avoids inf).
        let values: Vec<f64> = (0..120)
            .map(|i| (i % period) as f64 * 2.0 + 50.0 + i as f64 * 0.1)
            .collect();

        let standalone = seasonal_diff_variance_reduction(&values, period);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();
        let trait_vr = sd
            .seasonal_features()
            .iter()
            .find(|(n, _)| *n == "seasonal_diff_variance_reduction")
            .unwrap()
            .1;

        assert_abs_diff_eq!(standalone, trait_vr, epsilon = 1e-12);
    }

    #[test]
    fn standalone_edge_cases() {
        // period == 0
        assert_abs_diff_eq!(seasonal_diff_strength(&[1.0, 2.0], 0), 0.0, epsilon = 1e-12);
        assert_abs_diff_eq!(
            seasonal_diff_variance_reduction(&[1.0, 2.0], 0),
            1.0,
            epsilon = 1e-12
        );

        // values.len() <= period
        assert_abs_diff_eq!(
            seasonal_diff_strength(&[1.0, 2.0, 3.0], 5),
            0.0,
            epsilon = 1e-12
        );
        assert_abs_diff_eq!(
            seasonal_diff_variance_reduction(&[1.0, 2.0, 3.0], 5),
            1.0,
            epsilon = 1e-12
        );
    }

    // ── Error cases ────────────────────────────────────────────────────

    #[test]
    fn new_period_zero_errors() {
        let result = SeasonalDifference::new(0);
        assert!(matches!(result, Err(ForecastError::InvalidParameter(_))));
    }

    #[test]
    fn fit_empty_data_errors() {
        let mut sd = SeasonalDifference::new(4).unwrap();
        let result = sd.fit_seasonal(&[], 4);
        assert!(matches!(result, Err(ForecastError::EmptyData)));
    }

    #[test]
    fn fit_insufficient_data_errors() {
        let mut sd = SeasonalDifference::new(10).unwrap();
        let values = vec![1.0; 10]; // need > period observations
        let result = sd.fit_seasonal(&values, 10);
        assert!(matches!(
            result,
            Err(ForecastError::InsufficientData { .. })
        ));
    }

    #[test]
    fn fit_period_zero_errors() {
        let mut sd = SeasonalDifference::new(1).unwrap();
        let result = sd.fit_seasonal(&[1.0, 2.0, 3.0], 0);
        assert!(matches!(result, Err(ForecastError::InvalidParameter(_))));
    }

    #[test]
    fn inverse_before_fit_errors() {
        let sd = SeasonalDifference::new(4).unwrap();
        let result = sd.inverse(&[1.0, 2.0]);
        assert!(matches!(result, Err(ForecastError::FitRequired { .. })));
    }

    // ── Fitted seasonal values ─────────────────────────────────────────

    #[test]
    fn fitted_seasonal_same_length_as_input() {
        let period = 5;
        let values: Vec<f64> = (0..30).map(|i| (i % period) as f64 * 2.0).collect();

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        assert_eq!(sd.fitted_seasonal().len(), values.len());
    }

    #[test]
    fn fitted_seasonal_is_periodic() {
        let period = 3;
        let values: Vec<f64> = (0..21).map(|i| (i % period) as f64 + 10.0).collect();

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let fitted = sd.fitted_seasonal();
        for i in period..fitted.len() {
            assert_abs_diff_eq!(fitted[i], fitted[i % period], epsilon = 1e-12);
        }
    }

    // ── Name ───────────────────────────────────────────────────────────

    #[test]
    fn name_is_correct() {
        let sd = SeasonalDifference::new(4).unwrap();
        assert_eq!(sd.seasonal_name(), "seasonal_difference");
    }

    // ── Differenced length ─────────────────────────────────────────────

    #[test]
    fn differenced_length() {
        let period = 6;
        let n = 30;
        let values: Vec<f64> = (0..n).map(|i| i as f64).collect();

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        assert_eq!(sd.differenced().len(), n - period);
    }

    // ── Refit resets state ──────────────────────────────────────────────

    #[test]
    fn refit_resets_state() {
        let mut sd = SeasonalDifference::new(4).unwrap();

        let values1: Vec<f64> = (0..20).map(|i| (i % 4) as f64).collect();
        sd.fit_seasonal(&values1, 4).unwrap();
        let fitted1 = sd.fitted_seasonal().to_vec();

        let values2: Vec<f64> = (0..20).map(|i| (i % 4) as f64 * 10.0).collect();
        sd.fit_seasonal(&values2, 4).unwrap();
        let fitted2 = sd.fitted_seasonal().to_vec();

        // Second fit should produce different results.
        let diff: f64 = fitted1
            .iter()
            .zip(fitted2.iter())
            .map(|(a, b)| (a - b).abs())
            .sum();
        assert!(diff > 1.0, "refitting should change fitted values");
    }

    // ── Clone preserves state ──────────────────────────────────────────

    #[test]
    fn clone_preserves_state() {
        let period = 5;
        let values = pure_seasonal(period, 6);

        let mut sd = SeasonalDifference::new(period).unwrap();
        sd.fit_seasonal(&values, period).unwrap();

        let cloned = sd.clone();
        assert_eq!(sd.differenced(), cloned.differenced());
        assert_eq!(sd.fitted_seasonal(), cloned.fitted_seasonal());
        assert_eq!(sd.predict_seasonal(10), cloned.predict_seasonal(10));
    }
}