u_analytics/detection/

spectral_residual.rs

//! Spectral residual saliency for one-shot anomaly scoring of a series.
//!
//! The spectral residual transform (Hou & Zhang 2007, brought to time series
//! by Ren et al. 2019) takes the log amplitude spectrum of the series,
//! subtracts its local average, and transforms back with the original phase.
//! What is left — the saliency map — is large exactly where the series has a
//! feature its spectrum cannot explain as regular structure: a spike, a
//! step, a dropout. No model is trained and no period is assumed.
//!
//! Each point is scored relative to the saliency of the points before it,
//! `score = (S_i − mean(S over the judgement window)) / mean(...)`; a point
//! is an anomaly when its score exceeds the threshold (Ren et al. use τ = 3)
//! and, so that a bump of a few thousandths on a flat series is not one, its
//! value stands at least `min_zscore` standard deviations from the level of
//! the window before it.
//!
//! The result also carries an **expected value** — the series with the
//! anomalies replaced by their neighbours, reconstructed from its
//! low-frequency components — and a band around it whose half-width is the
//! normal quantile for `sensitivity` percent coverage times the robust
//! spread of the residuals. The band is information for a chart, not the
//! anomaly decision, which is the saliency score.
//!
//! # What the saliency does at a boundary, and at a step
//!
//! Two properties of the method are worth knowing before reading its output.
//!
//! **Near either end of a batch**, the saliency is moved by the transform's
//! own boundary handling: κ = 5 points along the trailing slope are appended
//! before the transform, and since the transform is circular that block is
//! adjacent to the first points as well as the last. Points there carry
//! [`SrPoint::near_edge`], so a flag that appears only at an end can be told
//! from one in the body of the series.
//!
//! **At a step**, the flag can fall on the point before the transition as well
//! as on the transition itself -- `[59, 60]` rather than `[60]` for a level
//! change between 59 and 60. Saliency is a measure of local spectral
//! irregularity, not an edge detector: the irregularity a step creates is
//! spread over the points around it, and which of them clear the threshold
//! depends on the noise. This is the method as published, not a defect, and it
//! is why a consumer matching detections against known transitions should
//! allow a place on either side rather than demand an exact index.
//!
//! # References
//! - Ren, H., Xu, B., Wang, Y., Yi, C., Huang, C., Kou, X., Xing, T., Yang,
//!   M., Tong, J. & Zhang, Q. (2019). "Time-Series Anomaly Detection Service
//!   at Microsoft." *KDD 2019*, 3009–3017.
//! - Hou, X. & Zhang, L. (2007). "Saliency Detection: A Spectral Residual
//!   Approach." *CVPR 2007*.

use u_numflow::fourier::{fft, ifft, Complex};
use u_numflow::special::inverse_normal_cdf;
use u_numflow::stats::median;

/// Fewest observations the transform is applied to at once.
pub const MIN_OBSERVATIONS: usize = 12;

// A refusal has to read as a sentence to whoever typed the option, so the
// `batch_size` message spells the number rather than naming the constant.
// This keeps the two from drifting apart.
const _: () = assert!(
    MIN_OBSERVATIONS == 12,
    "the batch_size requirement message spells this number"
);

/// Points appended before the transform to soften its boundary effect, and
/// the number of preceding points their slope is estimated from (κ = m = 5
/// in Ren et al.).
const EXTENSION: usize = 5;

/// Why a [`SpectralResidual`] run was refused.
///
/// One condition, not the rulebook. A consumer has to tell its user which
/// setting to change, and a refusal that lists every rule leaves it
/// re-validating the options the crate already checks -- which is how the
/// same rules end up written twice, in two places that can drift.
#[derive(Debug, Clone, Copy, PartialEq)]
#[non_exhaustive]
pub enum SpectralResidualError {
    /// An option is outside its domain.
    ///
    /// Each end of the domain is given once: `min` when the lower bound is
    /// inclusive, `exclusive_min` when it is exclusive, neither when there is
    /// none (likewise `max` / `exclusive_max`).
    OptionOutOfRange {
        /// The option's name, spelled as the builder and the wire formats
        /// spell it (`threshold`, `sensitivity`, ...).
        option: &'static str,
        /// What that option has to satisfy, for a message a consumer can show.
        requirement: &'static str,
        /// The value given (a count option's value as a float).
        got: f64,
        /// Inclusive lower bound.
        min: Option<f64>,
        /// Exclusive lower bound.
        exclusive_min: Option<f64>,
        /// Inclusive upper bound.
        max: Option<f64>,
        /// Exclusive upper bound.
        exclusive_max: Option<f64>,
    },
    /// Fewer observations than the transform is applied to at once.
    TooFewObservations {
        /// [`MIN_OBSERVATIONS`].
        needed: usize,
        /// How many were given.
        got: usize,
    },
    /// The series holds a NaN or an infinity.
    ValueNotFinite {
        /// Position of the first such value.
        index: usize,
    },
}

impl core::fmt::Display for SpectralResidualError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            SpectralResidualError::OptionOutOfRange {
                option,
                requirement,
                ..
            } => write!(f, "{option} must be {requirement}"),
            SpectralResidualError::TooFewObservations { needed, got } => {
                write!(f, "needs at least {needed} observations, got {got}")
            }
            SpectralResidualError::ValueNotFinite { index } => {
                write!(f, "data[{index}] is not a finite number")
            }
        }
    }
}

impl std::error::Error for SpectralResidualError {}

/// The options that count something, each with its floor.
///
/// The builder takes them as `usize`, but a wire format reads a number, which
/// can be negative; checking that number here, against the same floor the
/// builder's value is checked against, keeps `-1` and `0` one refusal.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) enum CountOption {
    AveragingWindow,
    JudgementWindow,
    BatchSize,
}

impl CountOption {
    /// `got` (a whole number) against this option's floor.
    pub(crate) fn check(self, got: f64) -> Result<(), SpectralResidualError> {
        let (option, requirement, min) = match self {
            CountOption::AveragingWindow => ("averaging_window", ">= 1", 1),
            CountOption::JudgementWindow => ("judgement_window", ">= 1", 1),
            CountOption::BatchSize => ("batch_size", "unset, or at least 12", MIN_OBSERVATIONS),
        };
        if got < min as f64 {
            Err(SpectralResidualError::below(
                option,
                requirement,
                got,
                min as f64,
            ))
        } else {
            Ok(())
        }
    }
}

impl SpectralResidualError {
    /// `option` below its inclusive lower bound `min`.
    fn below(option: &'static str, requirement: &'static str, got: f64, min: f64) -> Self {
        Self::OptionOutOfRange {
            option,
            requirement,
            got,
            min: Some(min),
            exclusive_min: None,
            max: None,
            exclusive_max: None,
        }
    }

    /// `option` not above its exclusive lower bound.
    fn not_above(
        option: &'static str,
        requirement: &'static str,
        got: f64,
        exclusive_min: f64,
    ) -> Self {
        Self::OptionOutOfRange {
            option,
            requirement,
            got,
            min: None,
            exclusive_min: Some(exclusive_min),
            max: None,
            exclusive_max: None,
        }
    }

    /// `option` not strictly between `exclusive_min` and `exclusive_max`.
    fn not_between(
        option: &'static str,
        requirement: &'static str,
        got: f64,
        exclusive_min: f64,
        exclusive_max: f64,
    ) -> Self {
        Self::OptionOutOfRange {
            option,
            requirement,
            got,
            min: None,
            exclusive_min: Some(exclusive_min),
            max: None,
            exclusive_max: Some(exclusive_max),
        }
    }
}

/// Spectral residual anomaly scorer.
///
/// Built with [`SpectralResidual::new`] (the defaults of Ren et al. 2019) and
/// adjusted with the `with_*` methods; run with [`SpectralResidual::analyze`].
#[derive(Debug, Clone, PartialEq)]
pub struct SpectralResidual {
    /// Width of the moving average applied to the log amplitude spectrum
    /// (`q` in Ren et al.; default 3).
    pub averaging_window: usize,
    /// How many preceding saliency values a point is compared against
    /// (`z`; default 40).
    pub judgement_window: usize,
    /// Score above which a point is an anomaly (`τ`; default 3).
    pub threshold: f64,
    /// Minimum `|value − mean| / std` of the point against the judgement
    /// window before it, for it to be an anomaly (default 1.5), so that
    /// saliency alone does not flag noise on a flat series. A window with no
    /// spread passes any departure from its level.
    pub min_zscore: f64,
    /// Coverage of the band around the expected value, in percent (default
    /// 70): its half-width is the normal quantile for that two-sided coverage
    /// times the robust spread of the non-anomalous residuals.
    pub sensitivity: f64,
    /// Process the series in consecutive batches of this many points instead
    /// of at once (default `None`). Long series with slow drift are scored
    /// against their local context this way; the last batch absorbs a short
    /// remainder.
    pub batch_size: Option<usize>,
}

impl Default for SpectralResidual {
    fn default() -> Self {
        Self::new()
    }
}

impl SpectralResidual {
    /// The defaults of Ren et al. (2019): q = 3, z = 40, τ = 3, z-score gate
    /// 1.5, 70% band, no batching.
    pub fn new() -> Self {
        Self {
            averaging_window: 3,
            judgement_window: 40,
            threshold: 3.0,
            min_zscore: 1.5,
            sensitivity: 70.0,
            batch_size: None,
        }
    }

    /// Sets the spectral averaging window (≥ 1).
    pub fn with_averaging_window(mut self, q: usize) -> Self {
        self.averaging_window = q;
        self
    }

    /// Sets the judgement window (≥ 1).
    pub fn with_judgement_window(mut self, z: usize) -> Self {
        self.judgement_window = z;
        self
    }

    /// Sets the anomaly threshold on the score (> 0).
    pub fn with_threshold(mut self, threshold: f64) -> Self {
        self.threshold = threshold;
        self
    }

    /// Sets the z-score gate (≥ 0; 0 disables it).
    pub fn with_min_zscore(mut self, min_zscore: f64) -> Self {
        self.min_zscore = min_zscore;
        self
    }

    /// Sets the band coverage in percent (0 < s < 100).
    pub fn with_sensitivity(mut self, sensitivity: f64) -> Self {
        self.sensitivity = sensitivity;
        self
    }

    /// Sets the batch size (≥ [`MIN_OBSERVATIONS`]), or `None` for one batch.
    pub fn with_batch_size(mut self, batch_size: Option<usize>) -> Self {
        self.batch_size = batch_size;
        self
    }

    /// The first option that is outside its domain, in a fixed order.
    fn check_options(&self) -> Result<(), SpectralResidualError> {
        use SpectralResidualError as E;
        CountOption::AveragingWindow.check(self.averaging_window as f64)?;
        CountOption::JudgementWindow.check(self.judgement_window as f64)?;
        if !self.threshold.is_finite() || self.threshold <= 0.0 {
            return Err(E::not_above(
                "threshold",
                "a finite number > 0",
                self.threshold,
                0.0,
            ));
        }
        if !self.min_zscore.is_finite() || self.min_zscore < 0.0 {
            return Err(E::below(
                "min_zscore",
                "a finite number >= 0",
                self.min_zscore,
                0.0,
            ));
        }
        if !self.sensitivity.is_finite() || self.sensitivity <= 0.0 || self.sensitivity >= 100.0 {
            return Err(E::not_between(
                "sensitivity",
                "a finite number strictly between 0 and 100",
                self.sensitivity,
                0.0,
                100.0,
            ));
        }
        if let Some(b) = self.batch_size {
            CountOption::BatchSize.check(b as f64)?;
        }
        Ok(())
    }

    /// Scores every point of the series.
    ///
    /// The result has one [`SrPoint`] per input point, in order.
    ///
    /// # Errors
    ///
    /// Returns the **one** condition that was not met — the option that is out
    /// of its domain, the shortfall in observations, or the position of the
    /// first non-finite value — rather than a refusal a caller has to match
    /// against the whole rulebook to interpret.
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::detection::{SpectralResidual, SpectralResidualError};
    ///
    /// let mut series: Vec<f64> = (0..60).map(|t| (t as f64 * 0.3).sin()).collect();
    /// series[40] += 4.0; // a spike
    /// let points = SpectralResidual::new().analyze(&series).unwrap();
    /// assert!(points[40].is_anomaly);
    /// assert!(points[40].score > points[20].score);
    ///
    /// // A refusal names the option it is about, the value and its domain.
    /// let refused = SpectralResidual::new().with_threshold(0.0).analyze(&series);
    /// assert_eq!(
    ///     refused,
    ///     Err(SpectralResidualError::OptionOutOfRange {
    ///         option: "threshold",
    ///         requirement: "a finite number > 0",
    ///         got: 0.0,
    ///         min: None,
    ///         exclusive_min: Some(0.0),
    ///         max: None,
    ///         exclusive_max: None,
    ///     })
    /// );
    /// ```
    pub fn analyze(&self, series: &[f64]) -> Result<Vec<SrPoint>, SpectralResidualError> {
        self.check_options()?;
        if series.len() < MIN_OBSERVATIONS {
            return Err(SpectralResidualError::TooFewObservations {
                needed: MIN_OBSERVATIONS,
                got: series.len(),
            });
        }
        if let Some(index) = series.iter().position(|v| !v.is_finite()) {
            return Err(SpectralResidualError::ValueNotFinite { index });
        }
        let n = series.len();
        let batch = self.batch_size.unwrap_or(n).min(n);
        let mut points = Vec::with_capacity(n);
        let mut start = 0;
        while start < n {
            // The last batch absorbs a remainder shorter than a batch.
            let end = if n - start < 2 * batch {
                n
            } else {
                start + batch
            };
            points.extend(self.analyze_batch(&series[start..end], start));
            start = end;
        }
        Ok(points)
    }

    fn analyze_batch(&self, data: &[f64], offset: usize) -> Vec<SrPoint> {
        let n = data.len();
        let saliency = saliency_map(data, self.averaging_window);

        // Score against the preceding judgement window (the point included,
        // as in Ren et al.'s implementation), relative to that average.
        let scores: Vec<f64> = (0..n)
            .map(|i| {
                let from = i.saturating_sub(self.judgement_window - 1);
                let avg = saliency[from..=i].iter().sum::<f64>() / (i + 1 - from) as f64;
                let avg = avg.max(1e-8);
                ((saliency[i] - avg) / avg).max(0.0)
            })
            .collect();

        // The z-score gate is local: the point against the mean and spread of
        // the judgement window before it, so a step is judged against the
        // level it departs from, not against a mean that already straddles it.
        let flagged: Vec<bool> = (0..n)
            .map(|i| {
                if scores[i] <= self.threshold {
                    return false;
                }
                if self.min_zscore == 0.0 || i == 0 {
                    return true;
                }
                let from = i.saturating_sub(self.judgement_window);
                let window = &data[from..i];
                let mean = window.iter().sum::<f64>() / window.len() as f64;
                let std = (window.iter().map(|v| (v - mean).powi(2)).sum::<f64>()
                    / window.len() as f64)
                    .sqrt();
                if std > 0.0 {
                    (data[i] - mean).abs() / std >= self.min_zscore
                } else {
                    data[i] != mean
                }
            })
            .collect();

        let expected = expected_values(data, &flagged);
        let residuals: Vec<f64> = (0..n)
            .filter(|&i| !flagged[i])
            .map(|i| (data[i] - expected[i]).abs())
            .collect();
        // Robust spread: 1.4826 · median absolute residual (the MAD of a
        // centred residual), then the two-sided normal quantile for the
        // requested coverage.
        let spread = median(&residuals).unwrap_or(0.0) * 1.4826;
        let z = inverse_normal_cdf(0.5 + self.sensitivity / 200.0);
        let margin = z * spread;

        (0..n)
            .map(|i| SrPoint {
                index: offset + i,
                value: data[i],
                saliency: saliency[i],
                score: scores[i],
                expected: expected[i],
                lower: expected[i] - margin,
                upper: expected[i] + margin,
                is_anomaly: flagged[i],
                near_edge: i < EXTENSION || i + EXTENSION >= n,
            })
            .collect()
    }
}

/// One scored point.
#[derive(Debug, Clone, PartialEq)]
pub struct SrPoint {
    /// Position in the input series.
    pub index: usize,
    /// The observed value.
    pub value: f64,
    /// Spectral residual saliency of the point (≥ 0).
    pub saliency: f64,
    /// Saliency relative to the preceding judgement window (≥ 0).
    pub score: f64,
    /// Low-frequency reconstruction of the series with anomalies removed.
    pub expected: f64,
    /// `expected − margin`.
    pub lower: f64,
    /// `expected + margin`.
    pub upper: f64,
    /// `score > threshold` and the value clears the z-score gate.
    pub is_anomaly: bool,
    /// The point lies within κ = 5 places of an end of its batch, where the
    /// transform's own boundary handling moves the saliency most.
    ///
    /// The method extends the series by κ points before transforming, and the
    /// transform is circular, so that appended block sits next to the last
    /// points *and* the first ones. Measured on a noisy sawtooth of 200
    /// points, replacing the predicted extension with the series' real
    /// continuation moves the saliency of the five points at either end about
    /// eight times as much as it moves the middle -- and the two ends by the
    /// same amount, which is why this is not a tail-only flag.
    ///
    /// It is a property of the position, not a verdict: a real anomaly near an
    /// end is still reported as one. What it says is that a lone flag here is
    /// the one worth a second look, and an anomaly count is worth reporting
    /// with and without these.
    pub near_edge: bool,
}

/// The saliency map of Ren et al. (2019): extend the series by κ = 5 points
/// along the slope of its last m = 5, take the FFT, average the log
/// amplitude over `q` bins, keep the residual with the original phase,
/// transform back, and drop the extension.
fn saliency_map(data: &[f64], q: usize) -> Vec<f64> {
    let n = data.len();
    let mut extended = data.to_vec();
    let m = EXTENSION.min(n - 1);
    let last = data[n - 1];
    let slope = (1..=m)
        .map(|i| (last - data[n - 1 - i]) / i as f64)
        .sum::<f64>()
        / m as f64;
    // x_{n+1} = x_{n-m+1} + g * m (Ren et al., eq. 9): the point m places before
    // the one being predicted, so `data[n - m]` with the last point `data[n - 1]`.
    let predicted = data[n - m] + slope * m as f64;
    extended.extend(std::iter::repeat_n(predicted, EXTENSION));

    let spectrum = fft(&extended
        .iter()
        .map(|&x| Complex::real(x))
        .collect::<Vec<_>>());
    let len = spectrum.len();
    // The floor that keeps the log finite at a zero bin is relative to the
    // spectrum's own size: a fixed one is nothing for a series in units of 1
    // and everything for the same series in nanoamperes (Ren et al.'s 1e-8 is
    // for data of order 1).
    let peak = spectrum.iter().map(|z| z.norm()).fold(0.0_f64, f64::max);
    let floor = if peak > 0.0 { 1e-8 * peak } else { 1e-8 };
    let log_amp: Vec<f64> = spectrum.iter().map(|z| (z.norm() + floor).ln()).collect();
    // Moving average of width q over the preceding bins (cumulative form).
    let mut avg = vec![0.0; len];
    let mut cum = 0.0;
    let mut cums = Vec::with_capacity(len + 1);
    cums.push(0.0);
    for &v in &log_amp {
        cum += v;
        cums.push(cum);
    }
    for i in 0..len {
        let from = (i + 1).saturating_sub(q);
        avg[i] = (cums[i + 1] - cums[from]) / (i + 1 - from) as f64;
    }
    let residual_spectrum: Vec<Complex> = spectrum
        .iter()
        .enumerate()
        .map(|(i, z)| {
            let r = (log_amp[i] - avg[i]).exp();
            let phase = z.arg();
            Complex::new(r * phase.cos(), r * phase.sin())
        })
        .collect();
    ifft(&residual_spectrum)
        .into_iter()
        .take(n)
        .map(|z| z.norm())
        .collect()
}

/// The series with each flagged point replaced by the mean of its nearest
/// unflagged neighbours, reconstructed from its low-frequency components:
/// the middle band of the spectrum (bins from 3/8 to 5/8 of the length) is
/// zeroed, as in the Microsoft service's expected value.
fn expected_values(data: &[f64], flagged: &[bool]) -> Vec<f64> {
    let n = data.len();
    let mut cleaned = data.to_vec();
    for i in 0..n {
        if !flagged[i] {
            continue;
        }
        let left = (0..i).rev().find(|&j| !flagged[j]).map(|j| data[j]);
        let right = (i + 1..n).find(|&j| !flagged[j]).map(|j| data[j]);
        cleaned[i] = match (left, right) {
            (Some(a), Some(b)) => (a + b) / 2.0,
            (Some(a), None) | (None, Some(a)) => a,
            (None, None) => data[i],
        };
    }
    let mut spectrum = fft(&cleaned
        .iter()
        .map(|&x| Complex::real(x))
        .collect::<Vec<_>>());
    let lo = n * 3 / 8;
    let hi = n * 5 / 8;
    for (i, z) in spectrum.iter_mut().enumerate() {
        if i > lo && i < hi {
            *z = Complex::default();
        }
    }
    ifft(&spectrum).into_iter().map(|z| z.re).collect()
}

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

    fn noise(seed: u64, n: usize) -> Vec<f64> {
        let mut state = seed;
        (0..n)
            .map(|_| {
                state = state
                    .wrapping_mul(6364136223846793005)
                    .wrapping_add(1442695040888963407);
                ((state >> 11) as f64 / (1u64 << 53) as f64) * 2.0 - 1.0
            })
            .collect()
    }

    fn check_shape(points: &[SrPoint], n: usize) {
        assert_eq!(points.len(), n);
        for (i, p) in points.iter().enumerate() {
            assert_eq!(p.index, i);
            assert!(p.saliency >= 0.0 && p.saliency.is_finite(), "{p:?}");
            assert!(p.score >= 0.0 && p.score.is_finite(), "{p:?}");
            assert!(
                p.expected.is_finite() && p.lower <= p.expected && p.expected <= p.upper,
                "{p:?}"
            );
        }
    }

    /// The same series in another unit is the same series: saliency, score and
    /// the anomalies are dimensionless, the expected value and band scale with
    /// the data.
    #[test]
    fn the_unit_of_the_series_does_not_change_the_answer() {
        let mut series: Vec<f64> = noise(3, 120)
            .iter()
            .enumerate()
            .map(|(t, e)| (t as f64 * 0.25).sin() + 0.1 * e)
            .collect();
        series[80] += 5.0;
        let base = SpectralResidual::new().analyze(&series).unwrap();
        for k in [1e-15, 1e-12, 1e-9, 1e-6, 1e6, 1e12] {
            let scaled: Vec<f64> = series.iter().map(|v| v * k).collect();
            let got = SpectralResidual::new().analyze(&scaled).unwrap();
            for (g, b) in got.iter().zip(&base) {
                assert_eq!(g.is_anomaly, b.is_anomaly, "scale {k:e}, point {}", b.index);
                assert!(
                    (g.saliency - b.saliency).abs() <= 1e-6 * b.saliency.max(1e-3),
                    "scale {k:e}, point {}: saliency {} vs {}",
                    b.index,
                    g.saliency,
                    b.saliency
                );
                assert!(
                    (g.expected / k - b.expected).abs() <= 1e-6 * b.expected.abs().max(1.0),
                    "scale {k:e}, point {}: expected",
                    b.index
                );
            }
        }
    }

    #[test]
    fn a_spike_on_a_sine_is_the_anomaly() {
        let mut series: Vec<f64> = noise(3, 120)
            .iter()
            .enumerate()
            .map(|(t, e)| (t as f64 * 0.25).sin() + 0.1 * e)
            .collect();
        series[80] += 5.0;
        let points = SpectralResidual::new().analyze(&series).unwrap();
        check_shape(&points, 120);
        assert!(points[80].is_anomaly, "{:?}", points[80]);
        let flagged: Vec<usize> = points
            .iter()
            .filter(|p| p.is_anomaly)
            .map(|p| p.index)
            .collect();
        assert!(flagged.len() <= 3, "flagged {flagged:?}");
        assert!(points[80].value > points[80].upper);
        // The expected value is not dragged up by the spike.
        assert!(
            (points[80].expected - (80.0f64 * 0.25).sin()).abs() < 0.6,
            "{:?}",
            points[80]
        );
    }

    #[test]
    fn a_step_is_flagged_where_it_happens() {
        let series: Vec<f64> = (0..100).map(|t| if t < 60 { 1.0 } else { 6.0 }).collect();
        let points = SpectralResidual::new().analyze(&series).unwrap();
        check_shape(&points, 100);
        let first = points.iter().position(|p| p.is_anomaly);
        assert_eq!(
            first,
            Some(60),
            "{:?}",
            points
                .iter()
                .filter(|p| p.is_anomaly)
                .map(|p| p.index)
                .collect::<Vec<_>>()
        );
    }

    #[test]
    fn a_constant_and_plain_noise_raise_no_anomaly() {
        let constant = vec![3.0; 50];
        let points = SpectralResidual::new().analyze(&constant).unwrap();
        check_shape(&points, 50);
        assert!(points.iter().all(|p| !p.is_anomaly));
        assert!(points.iter().all(|p| (p.expected - 3.0).abs() < 1e-9));

        let series = noise(9, 200);
        let points = SpectralResidual::new().analyze(&series).unwrap();
        check_shape(&points, 200);
        let flagged = points.iter().filter(|p| p.is_anomaly).count();
        assert!(flagged <= 4, "{flagged} of 200 noise points flagged");
    }

    #[test]
    fn scoring_is_deterministic_and_batching_keeps_the_shape() {
        let mut series: Vec<f64> = noise(5, 300).iter().map(|e| 10.0 + e).collect();
        series[150] = 20.0;
        series[260] = -5.0;
        let sr = SpectralResidual::new();
        let a = sr.analyze(&series).unwrap();
        let b = sr.analyze(&series).unwrap();
        assert_eq!(a, b);
        let batched = sr.with_batch_size(Some(100)).analyze(&series).unwrap();
        check_shape(&batched, 300);
        assert!(batched[150].is_anomaly && batched[260].is_anomaly);
        // Each batch is scored on its own: the first point of a batch has no
        // history but the shape is unchanged and indices are global.
        assert_eq!(batched[200].index, 200);
    }

    #[test]
    fn the_band_widens_with_sensitivity() {
        let series: Vec<f64> = noise(1, 80).iter().map(|e| 5.0 + e).collect();
        let narrow = SpectralResidual::new()
            .with_sensitivity(50.0)
            .analyze(&series)
            .unwrap();
        let wide = SpectralResidual::new()
            .with_sensitivity(99.0)
            .analyze(&series)
            .unwrap();
        for (a, b) in narrow.iter().zip(&wide) {
            assert!(b.upper - b.lower > a.upper - a.lower);
            assert_eq!(a.expected, b.expected);
        }
    }

    /// The transform's boundary handling reaches both ends, so the marker
    /// does too -- a tail-only flag would leave half the affected points
    /// unmarked.
    #[test]
    fn points_near_either_end_of_a_batch_are_marked() {
        let series: Vec<f64> = (0..200).map(|i| (i % 7) as f64).collect();
        let points = SpectralResidual::new().analyze(&series).expect("ok");
        let marked: Vec<usize> = points
            .iter()
            .filter(|p| p.near_edge)
            .map(|p| p.index)
            .collect();
        assert_eq!(marked, vec![0, 1, 2, 3, 4, 195, 196, 197, 198, 199]);

        // Batching moves the boundaries, so it moves the marks with them.
        let batched = SpectralResidual::new()
            .with_batch_size(Some(50))
            .analyze(&series)
            .expect("ok");
        let marked: Vec<usize> = batched
            .iter()
            .filter(|p| p.near_edge)
            .map(|p| p.index)
            .collect();
        // 49 and 50 straddle a batch boundary, so both are near an edge; 75
        // sits in the body of a batch and is not.
        assert!(marked.contains(&49) && marked.contains(&50), "{marked:?}");
        assert!(!marked.contains(&75), "{marked:?}");
        assert_eq!(
            marked.len(),
            4 * 10,
            "four batches, ten marks each: {marked:?}"
        );

        // It marks a position, it does not suppress a detection: a spike
        // planted in the last few points is still an anomaly.
        let mut spiked: Vec<f64> = (0..60).map(|i| (i % 7) as f64).collect();
        spiked[57] += 30.0;
        let points = SpectralResidual::new().analyze(&spiked).expect("ok");
        assert!(points[57].is_anomaly);
        assert!(points[57].near_edge);
    }

    /// Documented, not corrected (see the module docs): the flag at a level
    /// change can land on the point before it as well.
    #[test]
    fn a_step_can_flag_the_point_before_the_transition() {
        let step: Vec<f64> = (0..120)
            .map(|i| {
                let jitter = ((i as f64 * 7.3891).sin() * 10000.0).fract().abs() * 0.2;
                (if i < 60 { 1.0 } else { 6.0 }) + jitter
            })
            .collect();
        let points = SpectralResidual::new().analyze(&step).expect("ok");
        let flagged: Vec<usize> = points
            .iter()
            .filter(|p| p.is_anomaly && !p.near_edge)
            .map(|p| p.index)
            .collect();
        assert!(!flagged.is_empty(), "the transition is detected");
        assert!(
            flagged.iter().all(|&i| (59..=61).contains(&i)),
            "detections sit at the transition or its neighbour: {flagged:?}"
        );
    }

    /// A refusal names the one condition that failed, so a consumer can point
    /// at the setting to change instead of restating every rule.
    #[test]
    fn a_refusal_names_the_condition_that_failed() {
        let ok = vec![1.0; 20];
        use SpectralResidualError as E;

        assert_eq!(
            SpectralResidual::new().analyze(&[1.0; 11]),
            Err(SpectralResidualError::TooFewObservations {
                needed: MIN_OBSERVATIONS,
                got: 11
            })
        );

        let mut series = vec![1.0; 20];
        series[7] = f64::INFINITY;
        assert_eq!(
            SpectralResidual::new().analyze(&series),
            Err(SpectralResidualError::ValueNotFinite { index: 7 })
        );

        // Every option, each naming itself, with the value given and the
        // bounds of its domain -- inclusive as `min`, exclusive as
        // `exclusive_min` / `exclusive_max`.
        assert_eq!(
            SpectralResidual::new().with_threshold(0.0).analyze(&ok),
            Err(E::not_above("threshold", "a finite number > 0", 0.0, 0.0))
        );
        let sensitivity = |got| {
            Err(E::not_between(
                "sensitivity",
                "a finite number strictly between 0 and 100",
                got,
                0.0,
                100.0,
            ))
        };
        assert_eq!(
            SpectralResidual::new().with_sensitivity(100.0).analyze(&ok),
            sensitivity(100.0)
        );
        assert_eq!(
            SpectralResidual::new().with_sensitivity(0.0).analyze(&ok),
            sensitivity(0.0)
        );
        assert_eq!(
            SpectralResidual::new()
                .with_batch_size(Some(5))
                .analyze(&ok),
            Err(E::below("batch_size", "unset, or at least 12", 5.0, 12.0))
        );
        assert_eq!(
            SpectralResidual::new()
                .with_averaging_window(0)
                .analyze(&ok),
            Err(E::below("averaging_window", ">= 1", 0.0, 1.0))
        );
        assert_eq!(
            SpectralResidual::new()
                .with_judgement_window(0)
                .analyze(&ok),
            Err(E::below("judgement_window", ">= 1", 0.0, 1.0))
        );
        assert_eq!(
            SpectralResidual::new().with_min_zscore(-1.0).analyze(&ok),
            Err(E::below("min_zscore", "a finite number >= 0", -1.0, 0.0))
        );
        // NaN is refused too; `got` is the NaN itself, so match the rest.
        assert!(matches!(
            SpectralResidual::new()
                .with_threshold(f64::NAN)
                .analyze(&ok),
            Err(E::OptionOutOfRange {
                option: "threshold",
                min: None,
                exclusive_min,
                max: None,
                exclusive_max: None,
                got,
                ..
            }) if got.is_nan() && exclusive_min == Some(0.0)
        ));

        // The message a transport shows repeats the option, not the rulebook.
        let message = SpectralResidual::new()
            .with_threshold(0.0)
            .analyze(&ok)
            .expect_err("refused")
            .to_string();
        assert_eq!(message, "threshold must be a finite number > 0");
        assert!(!message.contains("sensitivity"));

        // `min_zscore = 0` disables the gate; it is not out of range.
        assert!(SpectralResidual::new()
            .with_min_zscore(0.0)
            .analyze(&ok)
            .is_ok());
    }
}

/// The saliency map and the scores against a NumPy implementation of the
/// published method (Ren et al. 2019, section 4.1): the series extended by
/// kappa = 5 copies of x_{n+1} = x_{n-m+1} + g * m (equations 8 and 9, with
/// m = 5), the log amplitude spectrum less its trailing average over q = 3
/// bins, the original phase, the inverse transform. One sine of period 16 with
/// seeded noise and a spike of 4 at point 30. The reference is a restatement
/// of the paper's equations, not Microsoft's code, whose extension leaves out
/// the last point.
#[cfg(test)]
mod numpy_reference {
    use super::*;

    const SERIES: [f64; 64] = [
        0.061, 0.585, 0.52, 1.134, 1.216, 0.847, 0.806, 0.369, -0.109, -0.265, -0.526, -1.008,
        -1.139, -0.838, -0.869, -0.672, -0.068, 0.575, 0.594, 0.806, 0.924, 0.716, 0.839, 0.411,
        0.217, -0.158, -0.583, -1.146, -1.238, -0.975, 3.168, -0.527, 0.231, 0.418, 0.492, 1.129,
        0.898, 1.212, 0.745, 0.17, -0.134, -0.308, -0.865, -0.867, -0.867, -0.913, -0.746, -0.503,
        0.119, 0.423, 0.811, 0.793, 0.878, 1.089, 0.49, 0.568, 0.154, -0.484, -0.677, -1.11,
        -0.991, -0.84, -0.504, -0.44,
    ];
    const SALIENCY: [f64; 64] = [
        0.0369614318108603,
        0.0796307697333925,
        0.0271892659280146,
        0.152493866186463,
        0.166141948303089,
        0.0816656918320778,
        0.107174248120332,
        0.0556332303417036,
        0.0251300576938508,
        0.00299138253005822,
        0.00594392958654813,
        0.106106695034459,
        0.119126565861323,
        0.0424769851751049,
        0.0893964151003082,
        0.0920927965046129,
        0.0198418341809379,
        0.0882438287553024,
        0.0549823475501431,
        0.067869463097445,
        0.113994918150788,
        0.105398733553807,
        0.119888725456503,
        0.0746552451898384,
        0.109038966340197,
        0.0173152365850345,
        0.0229006048210958,
        0.132971301432821,
        0.205056469142638,
        0.150938526823712,
        0.830517907514111,
        0.128639067122104,
        0.0093039215691647,
        0.0611436432409572,
        0.0458141563355618,
        0.143893509514701,
        0.129454085999137,
        0.179180669344178,
        0.105280578308881,
        0.0679545755367448,
        0.0275280318011067,
        0.0125531569333055,
        0.0908993020495842,
        0.0747621046819059,
        0.0497022369555018,
        0.0789105027673606,
        0.0665323274213843,
        0.0500154087434401,
        0.0336638535386336,
        0.0524451689724362,
        0.0975047615363696,
        0.0534062948099107,
        0.0728496986560913,
        0.133492464750329,
        0.0136337226253535,
        0.068520339745976,
        0.0280483371477404,
        0.102625433672842,
        0.106450145288677,
        0.183534057756558,
        0.135468059382387,
        0.0877555778478272,
        0.0278059231072768,
        0.0277019054815613,
    ];
    const SCORE: [f64; 64] = [
        0.0,
        0.365970771264122,
        0.0,
        1.05881285226553,
        0.796450465672768,
        0.0,
        0.15195611134601,
        0.0,
        0.0,
        0.0,
        0.0,
        0.503171625125652,
        0.602838806095091,
        0.0,
        0.22119300766101,
        0.238060922935183,
        0.0,
        0.223493163232391,
        0.0,
        0.0,
        0.559450528883167,
        0.413466232376651,
        0.566393797404165,
        0.0,
        0.402200065392971,
        0.0,
        0.0,
        0.758500430509727,
        1.5606529881759,
        0.830851346829295,
        6.79293349706241,
        0.199289129287858,
        0.0,
        0.0,
        0.0,
        0.402864454850786,
        0.253212692235468,
        0.701706509860138,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.315186058686334,
        0.0,
        0.0,
        0.0,
        0.0305605858088266,
        0.0553323044375994,
        0.76882472342237,
        0.298864697731819,
        0.0,
        0.0,
        0.0,
    ];

    /// The published method floors the amplitude at an absolute 1e-8; this
    /// crate at 1e-8 of the spectrum's peak (the same for data of order one,
    /// independent of the unit otherwise), which moves this series' values by
    /// less than 1e-7.
    fn close(got: f64, want: f64) -> bool {
        (got - want).abs() <= 1e-6 * want.abs().max(1.0)
    }

    #[test]
    fn saliency_and_scores_reproduce_the_published_method() {
        let points = SpectralResidual::new()
            .analyze(&SERIES)
            .expect("a valid series");
        for (i, p) in points.iter().enumerate() {
            assert!(
                close(p.saliency, SALIENCY[i]),
                "saliency[{i}]: got {}, reference {}",
                p.saliency,
                SALIENCY[i]
            );
            assert!(
                close(p.score, SCORE[i]),
                "score[{i}]: got {}, reference {}",
                p.score,
                SCORE[i]
            );
        }
        let flagged: Vec<usize> = points
            .iter()
            .filter(|p| p.is_anomaly)
            .map(|p| p.index)
            .collect();
        assert_eq!(flagged, vec![30]);
    }
}