struktura 1.8.7

#![allow(clippy::needless_range_loop, clippy::type_complexity, clippy::large_enum_variant, clippy::useless_vec, clippy::map_clone)]
//! Time-series anomaly detection with no training data.
//!
//! Struktura detects when the *structure* of a signal changes: detrended
//! fluctuation analysis (DFA) for the correlation structure, plus a streaming
//! monitor that calibrates itself on the first rows of your data.
//!
//! # Quick start
//!
//! ```
//! use struktura::{compare, is_degraded};
//!
//! # let normal_readings = vec![1.0; 256];
//! # let current_readings = vec![1.0; 256];
//! // Compare current readings against a known-good baseline
//! let result = compare(&normal_readings, &current_readings);
//! println!("{}", result); // "HEALTHY shift=+0.003" or "CRITICAL shift=-0.45"
//!
//! // Or just ask: is this signal degraded compared to baseline?
//! if is_degraded(&normal_readings, &current_readings) {
//!     trigger_alert();
//! }
//! # fn trigger_alert() {}
//! ```
//!
//! # Domains
//!
//! - [`space`]: spacecraft telemetry monitoring (reaction wheels, magnetometers, batteries)
//! - [`market`]: financial regime detection (trending / random walk / mean-reverting)
//! - [`text`]: writing rhythm analysis (human literary prose vs mechanical/AI)
//! - [`rhythm`]: event timing analysis (git commits, heartbeats, keystrokes)
//!
//! Works in `no_std` environments (`default-features = false`). 85-112x faster than Python.

#![cfg_attr(not(feature = "std"), no_std)]

#[cfg(not(feature = "std"))]
extern crate alloc;

#[cfg(not(feature = "std"))]
use alloc::vec;
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;

use core::fmt;

#[cfg(not(feature = "std"))]
fn ln(x: f64) -> f64 { libm::log(x) }
#[cfg(feature = "std")]
fn ln(x: f64) -> f64 { x.ln() }

#[cfg(not(feature = "std"))]
pub(crate) fn sqrt(x: f64) -> f64 { libm::sqrt(x) }
#[cfg(feature = "std")]
pub(crate) fn sqrt(x: f64) -> f64 { x.sqrt() }

#[cfg(not(feature = "std"))]
fn powf(x: f64, y: f64) -> f64 { libm::pow(x, y) }
#[cfg(feature = "std")]
fn powf(x: f64, y: f64) -> f64 { x.powf(y) }

#[cfg(not(feature = "std"))]
fn powi(x: f64, n: i32) -> f64 { libm::pow(x, n as f64) }
#[cfg(feature = "std")]
fn powi(x: f64, n: i32) -> f64 { x.powi(n) }

#[cfg(not(feature = "std"))]
fn sin(x: f64) -> f64 { libm::sin(x) }
#[cfg(feature = "std")]
fn sin(x: f64) -> f64 { x.sin() }

#[cfg(not(feature = "std"))]
#[allow(dead_code)]
fn cos(x: f64) -> f64 { libm::cos(x) }
#[cfg(feature = "std")]
#[allow(dead_code)]
fn cos(x: f64) -> f64 { x.cos() }

/// Result of a DFA or ACR computation.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct DfaResult {
    /// Scaling exponent (slope in log-log space).
    pub alpha: f64,
    /// Coefficient of determination of the log-log fit.
    pub r_squared: f64,
}

impl fmt::Display for DfaResult {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "alpha={:.3} R2={:.4}", self.alpha, self.r_squared)
    }
}

/// How confident the analysis is in the derived scaling exponent.
///
/// Determined by the R-squared of the log-log fit. Higher R-squared means
/// the scaling law is a better fit to the data.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[non_exhaustive]
pub enum LawQuality {
    /// R-squared > 0.95: the scaling law fits the data almost perfectly.
    Exact,
    /// R-squared > 0.85: strong confidence in the derived exponent.
    Strong,
    /// R-squared > 0.7: good enough for health monitoring.
    Good,
    /// R-squared > 0.3: approximate; use with caution.
    Approx,
    /// R-squared <= 0.3: insufficient structure; the crate abstains from diagnosis.
    Abstain,
    /// Fewer than 20 data points: not enough data to analyze.
    Insufficient,
}

/// Complete structural analysis of a time series.
///
/// Contains the DFA and ACR results plus distributional statistics.
/// Use [`analyze`] to compute this from raw data.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct StructuralLaw {
    /// Hurst exponent estimated from ACR decay: H = 1 + acr_exponent/2.
    pub hurst: f64,
    /// Detrended Fluctuation Analysis result.
    pub dfa: DfaResult,
    /// Autocorrelation decay result.
    pub acr: DfaResult,
    /// Arithmetic mean of the signal.
    pub mean: f64,
    /// Standard deviation of the signal.
    pub std_dev: f64,
    /// Kurtosis (4th moment). Values > 4 indicate heavy tails / bursty behavior.
    pub kurtosis: f64,
    /// 99th percentile value.
    pub p99: f64,
    /// Maximum observed value.
    pub max: f64,
    /// Number of samples analyzed.
    pub n: usize,
    /// Confidence classification of the analysis.
    pub quality: LawQuality,
}

/// Health verdict comparing current DFA alpha against a known baseline.
///
/// Thresholds: Healthy < 0.03, Watch < 0.08, Warning < 0.15, Critical >= 0.15.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[non_exhaustive]
pub enum HealthVerdict {
    /// Shift < 0.03 from baseline: within normal variation.
    Healthy,
    /// Shift 0.03-0.08: minor structural change, monitor closely.
    Watch,
    /// Shift 0.08-0.15: significant structural departure.
    Warning,
    /// Shift >= 0.15: major structural breakdown.
    Critical,
}

impl HealthVerdict {
    /// Classify a DFA alpha shift into a health verdict.
    ///
    /// Thresholds: |shift| < 0.03 = Healthy, < 0.08 = Watch, < 0.15 = Warning, else Critical.
    pub fn from_shift(shift: f64) -> Self {
        let s = if shift < 0.0 { -shift } else { shift };
        if s < 0.03 {
            HealthVerdict::Healthy
        } else if s < 0.08 {
            HealthVerdict::Watch
        } else if s < 0.15 {
            HealthVerdict::Warning
        } else {
            HealthVerdict::Critical
        }
    }
}

/// Compute the DFA scaling exponent of a time series.
///
/// Returns alpha (the scaling exponent) and R-squared (fit quality).
/// Requires at least 64 data points.
///
/// **Warning:** Any upstream preprocessing (filtering, detrending, artifact
/// rejection) changes the signal's correlation structure and shifts alpha.
/// Baselines calibrated on raw data are invalid after adding preprocessing.
///
/// ```
/// use struktura::dfa;
/// let noise: Vec<f64> = (0..256).map(|i| (i as f64 * 0.1).sin()).collect();
/// let result = dfa(&noise);
/// assert!(result.r_squared >= 0.0);
/// ```
///
/// Below 64 samples this returns the placeholder `alpha: 0.5, r_squared: 0.0`,
/// which is not a measurement. Check `r_squared`, or use [`dfa_short`] for
/// short series. From 64 to roughly 128 samples the box sizes (from 16 to
/// n/4, see [`dfa_box_sizes`]) span a narrow range and α is noisy (values
/// near 10 have been seen on 72-sample random walks); prefer [`dfa_short`]
/// there too.
#[must_use]
pub fn dfa(values: &[f64]) -> DfaResult {
    let n = values.len();
    if n < 64 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }
    let mut buf = Vec::with_capacity(n);
    dfa_into(values, &mut buf)
}

/// DFA for short series (from about 24 samples), or `None` when the series
/// cannot be measured (too short, constant, or fewer than 3 usable box sizes).
///
/// Box sizes are up to 10 log-spaced integers from 4 to n/4. Boxes are laid
/// from the start and again from the end of the profile (Kantelhardt et al.,
/// Physica A 316, 2002), so the tail is used; each box is linearly detrended and
/// F(s) is the root mean square over all 2 * (n / s) boxes. [`dfa`] uses
/// boxes of at least 16 samples from the start only, which leaves too few
/// sizes below a few hundred samples.
///
/// On the ESA OPS-SAT-AD test split (529 segments, median 70 samples),
/// per-channel |z| of this alpha against nominal training segments gives
/// AUC-ROC 0.943; `dfa` gives 0.770 because half the segments are under 64
/// samples. Shuffling each segment's values drops it to 0.554 overall, but
/// only to 0.743 on the segments under 64 samples, so there part of the
/// separation is not order structure (examples/opssat_eval.rs).
///
/// ```
/// use struktura::dfa_short;
/// let x: Vec<f64> = (0..70).map(|i| (i as f64 * 0.37).sin() + (i % 7) as f64 * 0.1).collect();
/// assert!(dfa_short(&x).is_some());
/// assert!(dfa_short(&[1.0; 70]).is_none());
/// assert!(dfa_short(&x[..10]).is_none());
/// ```
#[must_use]
pub fn dfa_short(values: &[f64]) -> Option<DfaResult> {
    let n = values.len();
    if n < 16 {
        return None;
    }
    let mean = values.iter().sum::<f64>() / n as f64;
    if values.iter().all(|&v| v == values[0]) {
        return None;
    }
    let mut profile = Vec::with_capacity(n);
    let mut cum = 0.0f64;
    for &v in values {
        cum += v - mean;
        profile.push(cum);
    }

    let hi = (n / 4).max(5);
    let ratio = powf(hi as f64 / 4.0, 1.0 / 9.0);
    let mut log_s = [0.0f64; 10];
    let mut log_f = [0.0f64; 10];
    let mut pts = 0usize;
    let mut prev_s = 0usize;
    for step in 0..10 {
        let s = (4.0 * powi(ratio, step) + 1e-9) as usize;
        if s == prev_s || s < 4 || s > hi {
            continue;
        }
        prev_s = s;
        let nb = n / s;
        if nb == 0 {
            continue;
        }
        let k = s as f64;
        let sx = k * (k - 1.0) / 2.0;
        let sx2 = k * (k - 1.0) * (2.0 * k - 1.0) / 6.0;
        let det = k * sx2 - sx * sx;
        let box_msr = |seg: &[f64]| {
            let (mut sy, mut sxy, mut sy2) = (0.0f64, 0.0f64, 0.0f64);
            for (j, &y) in seg.iter().enumerate() {
                sy += y;
                sxy += j as f64 * y;
                sy2 += y * y;
            }
            let a0 = (sx2 * sy - sx * sxy) / det;
            let a1 = (k * sxy - sx * sy) / det;
            (sy2 - a0 * sy - a1 * sxy).max(0.0) / k
        };
        let mut f2 = 0.0f64;
        for b in 0..nb {
            f2 += box_msr(&profile[b * s..(b + 1) * s]);
            f2 += box_msr(&profile[n - (b + 1) * s..n - b * s]);
        }
        let f = sqrt(f2 / (2 * nb) as f64);
        if f > 0.0 {
            log_s[pts] = ln(s as f64);
            log_f[pts] = ln(f);
            pts += 1;
        }
    }
    if pts < 3 {
        return None;
    }
    Some(linreg(&log_s[..pts], &log_f[..pts]))
}

/// Box sizes that [`dfa`], [`dfa_into`], [`dfa_fast_into`] and
/// [`dfa_scratch`] use for a series of `n` samples: up to 12 distinct sizes,
/// geometrically spaced from `max(16, n / 50)` to `n / 4`. Returns the sizes
/// and how many of them are valid (0 when the range is empty).
///
/// The C99 generators in [`codegen`] emit this list, so generated monitors
/// measure alpha on the same scales as the Rust code that calibrates them.
///
/// ```
/// let (sizes, count) = struktura::dfa_box_sizes(96);
/// assert_eq!(&sizes[..count], &[16, 17, 18, 19, 20, 21, 22, 23, 24]);
/// ```
#[must_use]
pub fn dfa_box_sizes(n: usize) -> ([usize; 12], usize) {
    let mut sizes = [0usize; 12];
    let s_min = 16usize.max(n / 50);
    let s_max = n / 4;
    if s_min >= s_max {
        return (sizes, 0);
    }
    let ratio = powf(s_max as f64 / s_min as f64, 1.0 / 11.0);
    let mut count = 0usize;
    let mut prev_s = 0usize;
    for step in 0..12 {
        let s = (s_min as f64 * powi(ratio, step)) as usize;
        if s == prev_s || s > s_max {
            continue;
        }
        prev_s = s;
        sizes[count] = s;
        count += 1;
    }
    (sizes, count)
}

#[must_use]
/// Prefix-sum DFA: identical boxes and mathematics to [`dfa_into`], but the
/// per-segment sums (Σy, Σj·y, Σy²) are O(1) prefix-difference lookups
/// instead of an O(s) pass per segment. One O(n) pass builds the profile
/// prefixes; each of the ≤12 box sizes then costs O(n/s) segments × O(1).
///
/// Total work: O(n + Σ n/s) versus O(n × sizes) for the naive loop.
///
/// Precision: prefix differences of Σy² cancel catastrophically only when
/// n is large enough that the prefix magnitude dwarfs a segment's sum; for
/// the streaming-monitor window sizes (≤ a few thousand samples) agreement
/// with [`dfa_into`] is at machine precision (verified to 1e-12 in tests).
/// `buf` is a scratch buffer, grown to 3(n+1) and reused across calls.
pub fn dfa_fast_into(values: &[f64], buf: &mut Vec<f64>) -> DfaResult {
    let n = values.len();
    if n < 64 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }
    let s_min = 16usize.max(n / 50);
    let s_max = n / 4;
    if s_min >= s_max {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }

    let mean = values.iter().sum::<f64>() / n as f64;

    // One pass: profile y_j = cumsum(x - mean), prefix arrays
    // Py[k] = Σ_{j<k} y_j, PJy[k] = Σ_{j<k} j·y_j, Py2[k] = Σ_{j<k} y_j².
    buf.clear();
    buf.resize(3 * (n + 1), 0.0);
    let (py, rest) = buf.split_at_mut(n + 1);
    let (pjy, py2) = rest.split_at_mut(n + 1);
    let mut cum = 0.0f64;
    let mut acc_y = 0.0f64;
    let mut acc_jy = 0.0f64;
    let mut acc_y2 = 0.0f64;
    py[0] = 0.0;
    pjy[0] = 0.0;
    py2[0] = 0.0;
    for (j, &v) in values.iter().enumerate() {
        cum += v - mean;
        acc_y += cum;
        acc_jy += j as f64 * cum;
        acc_y2 += cum * cum;
        py[j + 1] = acc_y;
        pjy[j + 1] = acc_jy;
        py2[j + 1] = acc_y2;
    }

    let (sizes, count) = dfa_box_sizes(n);
    let mut log_s = [0.0f64; 12];
    let mut log_f = [0.0f64; 12];
    let mut pts = 0usize;

    for &s in &sizes[..count] {
        let num_segs = n / s;
        if num_segs == 0 {
            continue;
        }
        let k = s as f64;
        let sx = k * (k - 1.0) / 2.0;
        let sx2 = k * (k - 1.0) * (2.0 * k - 1.0) / 6.0;
        let det = k * sx2 - sx * sx;
        if det.abs() < 1e-15 {
            continue;
        }
        let mut f2_sum = 0.0;
        for seg in 0..num_segs {
            let a = seg * s;
            let b = a + s;
            let sy = py[b] - py[a];
            // local x = j - a inside the segment
            let sxy = (pjy[b] - pjy[a]) - a as f64 * sy;
            let sy2 = py2[b] - py2[a];
            let a0 = (sx2 * sy - sx * sxy) / det;
            let a1 = (k * sxy - sx * sy) / det;
            let resid = (sy2 - a0 * sy - a1 * sxy).max(0.0);
            f2_sum += resid / k;
        }
        let f = sqrt(f2_sum / num_segs as f64);
        if f > 0.0 {
            log_s[pts] = ln(s as f64);
            log_f[pts] = ln(f);
            pts += 1;
        }
    }

    if pts < 3 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }
    linreg(&log_s[..pts], &log_f[..pts])
}

/// DFA with a caller-provided buffer, avoiding allocation on the hot path.
///
/// `buf` is resized to `values.len()` and used for the cumulative sum.
/// Pre-allocate once and reuse across calls. For a heap-free build use
/// [`dfa_scratch`], which takes a plain slice.
#[must_use]
pub fn dfa_into(values: &[f64], buf: &mut Vec<f64>) -> DfaResult {
    let n = values.len();
    if n < 64 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }
    buf.clear();
    buf.resize(n, 0.0);
    dfa_scratch(values, buf)
}

/// Allocation-free DFA: the same mathematics as [`dfa_into`], with the
/// cumulative profile written into a caller-owned slice. Runs on `no_std`
/// targets with no heap at all; a stack array or a `static mut` buffer works.
///
/// `scratch` must hold at least `values.len()` elements; only the first
/// `values.len()` are written. Returns the neutral result
/// (`alpha = 0.5`, `r_squared = 0.0`) for fewer than 64 samples or a scratch
/// slice that is too short.
///
/// ```
/// use struktura::dfa_scratch;
/// let signal: [f64; 256] = core::array::from_fn(|i| (i as f64 * 0.1).sin());
/// let mut scratch = [0.0f64; 256];
/// let r = dfa_scratch(&signal, &mut scratch);
/// assert!(r.r_squared >= 0.0);
/// ```
#[must_use]
pub fn dfa_scratch(values: &[f64], scratch: &mut [f64]) -> DfaResult {
    let n = values.len();
    if n < 64 || scratch.len() < n {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }

    let mean = values.iter().sum::<f64>() / n as f64;

    let mut cum = 0.0;
    for (slot, &v) in scratch[..n].iter_mut().zip(values) {
        cum += v - mean;
        *slot = cum;
    }
    let buf = &scratch[..n];

    // Adaptive box sizes: geometric spacing from max(16, n/50) to n/4.
    // Gives consistent accuracy across signal lengths: short signals
    // get tighter boxes, long signals get wider coverage.
    let (sizes, count) = dfa_box_sizes(n);
    if count == 0 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }

    let mut log_s = [0.0f64; 12];
    let mut log_f = [0.0f64; 12];
    let mut pts = 0usize;

    for &s in &sizes[..count] {
        let num_segs = n / s;
        if num_segs == 0 { continue; }

        // Precompute sx, sx2, det (depend only on s, not data).
        let k = s as f64;
        let sx = k * (k - 1.0) / 2.0;
        let sx2 = k * (k - 1.0) * (2.0 * k - 1.0) / 6.0;
        let det = k * sx2 - sx * sx;
        if det.abs() < 1e-15 { continue; }

        let mut f2_sum = 0.0;
        for seg in 0..num_segs {
            let start = seg * s;
            // Single pass: accumulate sy, sxy, sy2, then use the least-squares
            // identity RSS = Σy² − a0Σy − a1Σxy (cross terms collapse via the
            // normal equations) instead of a second residual pass.
            let mut sy = 0.0;
            let mut sxy = 0.0;
            let mut sy2 = 0.0;
            for i in 0..s {
                let yi = buf[start + i];
                sy += yi;
                sxy += i as f64 * yi;
                sy2 += yi * yi;
            }
            let a0 = (sx2 * sy - sx * sxy) / det;
            let a1 = (k * sxy - sx * sy) / det;
            let resid = (sy2 - a0 * sy - a1 * sxy).max(0.0);
            f2_sum += resid / k;
        }
        let f = sqrt(f2_sum / num_segs as f64);
        if f > 0.0 {
            log_s[pts] = ln(s as f64);
            log_f[pts] = ln(f);
            pts += 1;
        }
    }

    if pts < 3 {
        return DfaResult { alpha: 0.5, r_squared: 0.0 };
    }

    linreg(&log_s[..pts], &log_f[..pts])
}

/// Compute autocorrelation decay exponent.
///
/// Measures how fast temporal correlations decay with lag.
/// Requires at least 20 data points.
#[must_use]
pub fn acr(values: &[f64]) -> DfaResult {
    let n = values.len();
    if n < 20 {
        return DfaResult { alpha: 0.0, r_squared: 0.0 };
    }

    let mean = values.iter().sum::<f64>() / n as f64;
    let var: f64 = values.iter().map(|&x| (x - mean) * (x - mean)).sum();
    if var < 1e-15 {
        return DfaResult { alpha: 0.0, r_squared: 0.0 };
    }

    const LAGS: [usize; 10] = [1, 2, 3, 5, 8, 13, 21, 34, 55, 89];
    let mut log_lag = [0.0f64; 10];
    let mut log_r = [0.0f64; 10];
    let mut pts = 0usize;

    for &lag in &LAGS {
        if lag >= n / 2 { break; }
        let mut num = 0.0;
        for i in 0..n - lag {
            num += (values[i] - mean) * (values[i + lag] - mean);
        }
        let r = num / var;
        if r > 0.001 {
            log_lag[pts] = ln(lag as f64);
            log_r[pts] = ln(r);
            pts += 1;
        }
    }

    if pts < 3 {
        return DfaResult { alpha: 0.0, r_squared: 0.0 };
    }

    linreg(&log_lag[..pts], &log_r[..pts])
}

/// Filter out NaN and Inf values from a signal.
///
/// Called automatically by [`analyze`]. You only need this if using [`dfa`] directly.
pub fn sanitize(values: &[f64]) -> Vec<f64> {
    values.iter().copied().filter(|v| v.is_finite()).collect()
}

/// Full structural analysis of a time series.
///
/// Computes DFA, ACR, Hurst exponent, kurtosis, and classifies law quality.
/// Automatically filters NaN/Inf and handles constant signals.
///
/// ```
/// use struktura::{analyze, LawQuality};
/// let data: Vec<f64> = (0..256).map(|i| (i as f64 * 0.07).sin() * 3.0).collect();
/// let law = analyze(&data);
/// assert!(law.n == 256);
/// assert!(law.quality != LawQuality::Insufficient);
/// ```
#[must_use]
pub fn analyze(values: &[f64]) -> StructuralLaw {
    let values = &sanitize(values);
    let n = values.len();
    if n < 20 {
        return StructuralLaw {
            hurst: 0.5, dfa: DfaResult { alpha: 0.5, r_squared: 0.0 },
            acr: DfaResult { alpha: 0.0, r_squared: 0.0 },
            mean: 0.0, std_dev: 0.0, kurtosis: 0.0, p99: 0.0, max: 0.0,
            n, quality: LawQuality::Insufficient,
        };
    }

    let mean = values.iter().sum::<f64>() / n as f64;
    let var: f64 = values.iter().map(|&x| (x - mean) * (x - mean)).sum::<f64>() / n as f64;
    let std_dev = sqrt(var);

    if std_dev < 1e-12 {
        return StructuralLaw {
            hurst: 0.5, dfa: DfaResult { alpha: 0.5, r_squared: 0.0 },
            acr: DfaResult { alpha: 0.0, r_squared: 0.0 },
            mean, std_dev: 0.0, kurtosis: 0.0, p99: mean, max: mean,
            n, quality: LawQuality::Abstain,
        };
    }

    let sd = std_dev;
    let kurtosis = values.iter().map(|&v| {
        let z = (v - mean) / sd;
        z * z * z * z
    }).sum::<f64>() / n as f64;

    let max = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    let mut sorted = values.to_vec();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
    let p99 = sorted[((n as f64 * 0.99) as usize).min(n - 1)];

    let dfa_result = dfa(values);
    let acr_result = acr(values);
    let hurst = clamp(1.0 + acr_result.alpha / 2.0, 0.0, 1.0);

    let best_r2 = if dfa_result.r_squared > acr_result.r_squared { dfa_result.r_squared } else { acr_result.r_squared };
    let quality = if best_r2 > 0.95 { LawQuality::Exact }
        else if best_r2 > 0.85 { LawQuality::Strong }
        else if best_r2 > 0.7 { LawQuality::Good }
        else if best_r2 > 0.3 { LawQuality::Approx }
        else { LawQuality::Abstain };

    StructuralLaw { hurst, dfa: dfa_result, acr: acr_result, mean, std_dev, kurtosis, p99, max, n, quality }
}

impl StructuralLaw {
    /// Returns true if the analysis has enough structure to be meaningful (not Abstain or Insufficient).
    pub fn is_healthy(&self) -> bool {
        self.quality != LawQuality::Abstain && self.quality != LawQuality::Insufficient
    }
}

impl DfaResult {
    /// Returns true if the log-log fit is strong enough to trust (R² > 0.7).
    pub fn is_reliable(&self) -> bool {
        self.r_squared > 0.7
    }
}

/// Deterministic Fisher-Yates shuffle of a signal.
///
/// Used for control experiments: shuffle destroys sequential structure
/// while preserving the value distribution. If DFA alpha changes after
/// shuffling, the original signal had genuine long-range correlations.
pub fn shuffle(values: &[f64], seed: u64) -> Vec<f64> {
    let mut out = values.to_vec();
    let n = out.len();
    let mut state = seed;
    for i in (1..n).rev() {
        state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        let j = (state >> 33) as usize % (i + 1);
        out.swap(i, j);
    }
    out
}

#[derive(Debug, Clone)]
pub struct ShuffleProof {
    pub real_alpha: f64,
    pub real_r2: f64,
    pub shuffled_alpha: f64,
    pub shuffled_r2: f64,
    pub structure_confirmed: bool,
}

pub fn prove_structure(values: &[f64]) -> ShuffleProof {
    let real = dfa(values);
    let shuffled_values = shuffle(values, 42);
    let shuffled = dfa(&shuffled_values);
    let real_dist = (real.alpha - 0.5).abs();
    let shuf_dist = (shuffled.alpha - 0.5).abs();
    ShuffleProof {
        real_alpha: real.alpha,
        real_r2: real.r_squared,
        shuffled_alpha: shuffled.alpha,
        shuffled_r2: shuffled.r_squared,
        structure_confirmed: shuf_dist < real_dist,
    }
}

#[derive(Debug, Clone)]
pub struct BootstrapCI {
    pub alpha: f64,
    pub ci_low: f64,
    pub ci_high: f64,
    pub n_resamples: usize,
}

/// Subsampling confidence interval for α (Politis–Romano style).
///
/// α is a scaling exponent of the ORDERING, so neither an i.i.d.
/// bootstrap (resampling values, which centres the interval near the
/// shuffled ~0.5 rather than the estimate) nor a moving-block bootstrap
/// (block joins destroy scaling above the block length; DFA boxes reach
/// n/4) is valid. Both were tried on a 123K-sample rover force channel
/// with α 0.728 and produced "intervals" of [0.39, 0.63] and [0.43, 0.66].
///
/// Instead, α is re-estimated on `n_resamples` overlapping contiguous
/// windows of `n/4` samples spread across the signal, and the interval is
/// `α ± 1.96 · sd(window α)`. Measured on 1/f noise, sd(α) does not fall
/// with window length above ~1K samples (the log-log fit uses a fixed
/// number of box sizes per decade), so the window spread is used
/// unscaled. The interval is centred on the estimate by construction;
/// its width is the signal's own α variability at quarter length.
pub fn bootstrap_alpha(values: &[f64], n_resamples: usize) -> BootstrapCI {
    let n = values.len();
    let base = dfa(values);
    let m = (n / 4).max(64).min(n);
    let k = n_resamples.max(2);
    let mut alphas = Vec::with_capacity(k);
    for i in 0..k {
        let start = if n > m { (i * (n - m)) / (k - 1) } else { 0 };
        let r = dfa(&values[start..start + m]);
        if r.r_squared > 0.3 { alphas.push(r.alpha); }
    }
    if alphas.len() < 2 {
        return BootstrapCI { alpha: base.alpha, ci_low: base.alpha, ci_high: base.alpha, n_resamples: alphas.len() };
    }
    let mean = alphas.iter().sum::<f64>() / alphas.len() as f64;
    let var = alphas.iter().map(|a| (a - mean) * (a - mean)).sum::<f64>() / (alphas.len() - 1) as f64;
    let half = 1.96 * sqrt(var);
    BootstrapCI { alpha: base.alpha, ci_low: base.alpha - half, ci_high: base.alpha + half, n_resamples: alphas.len() }
}

impl fmt::Display for BootstrapCI {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{:.3} [{:.3}, {:.3}] (n={})", self.alpha, self.ci_low, self.ci_high, self.n_resamples)
    }
}

#[derive(Debug, Clone)]
pub struct SplitHalfResult {
    pub first_half_alpha: f64,
    pub second_half_alpha: f64,
    pub difference: f64,
    pub consistent: bool,
}

pub fn split_half_validate(values: &[f64]) -> SplitHalfResult {
    let mid = values.len() / 2;
    let a = dfa(&values[..mid]);
    let b = dfa(&values[mid..]);
    let diff = (a.alpha - b.alpha).abs();
    SplitHalfResult {
        first_half_alpha: a.alpha,
        second_half_alpha: b.alpha,
        difference: diff,
        consistent: diff < 0.1,
    }
}

impl fmt::Display for SplitHalfResult {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "half1={:.3} half2={:.3} delta={:.3} {}",
            self.first_half_alpha, self.second_half_alpha, self.difference,
            if self.consistent { "CONSISTENT" } else { "INCONSISTENT" })
    }
}

impl fmt::Display for ShuffleProof {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "real={:.3} shuffled={:.3} {}",
            self.real_alpha, self.shuffled_alpha,
            if self.structure_confirmed { "CONFIRMED" } else { "INCONCLUSIVE" })
    }
}

/// Compare current DFA alpha against a known healthy baseline.
///
/// Returns a [`HealthVerdict`] based on how far alpha shifted from baseline.
///
/// ```
/// use struktura::{analyze, health_check, HealthVerdict};
/// let data: Vec<f64> = (0..256).map(|i| (i as f64 * 0.07).sin()).collect();
/// let law = analyze(&data);
/// let verdict = health_check(&law, 0.5);
/// // verdict is one of: Healthy, Watch, Warning, Critical
/// ```
#[must_use]
pub fn health_check(current: &StructuralLaw, baseline_alpha: f64) -> HealthVerdict {
    HealthVerdict::from_shift(current.dfa.alpha - baseline_alpha)
}

fn clamp(v: f64, lo: f64, hi: f64) -> f64 {
    if v < lo { lo } else if v > hi { hi } else { v }
}

fn linreg(x: &[f64], y: &[f64]) -> DfaResult {
    let k = x.len() as f64;
    let (mut sx, mut sy, mut sxy, mut sx2) = (0.0, 0.0, 0.0, 0.0);
    for i in 0..x.len() {
        sx += x[i]; sy += y[i]; sxy += x[i] * y[i]; sx2 += x[i] * x[i];
    }
    let slope = (k * sxy - sx * sy) / (k * sx2 - sx * sx);
    let ic = (sy - slope * sx) / k;
    let ym = sy / k;
    let mut sst = 0.0;
    let mut ssr = 0.0;
    for i in 0..x.len() {
        sst += (y[i] - ym) * (y[i] - ym);
        ssr += (y[i] - slope * x[i] - ic) * (y[i] - slope * x[i] - ic);
    }
    let r2 = 1.0 - ssr / if sst > 1e-15 { sst } else { 1e-15 };
    DfaResult { alpha: slope, r_squared: r2 }
}

pub struct SlidingWindow {
    buffer: Vec<f64>,
    capacity: usize,
    pos: usize,
    filled: bool,
}

impl SlidingWindow {
    pub fn new(capacity: usize) -> Self {
        SlidingWindow {
            buffer: vec![0.0; capacity],
            capacity,
            pos: 0,
            filled: false,
        }
    }

    pub fn push(&mut self, value: f64) {
        self.buffer[self.pos] = value;
        self.pos += 1;
        if self.pos >= self.capacity {
            self.pos = 0;
            self.filled = true;
        }
    }

    pub fn is_ready(&self) -> bool {
        self.filled
    }

    #[must_use]
pub fn analyze(&self) -> StructuralLaw {
        if !self.filled {
            return analyze(&self.buffer[..self.pos]);
        }
        let mut ordered = Vec::with_capacity(self.capacity);
        ordered.extend_from_slice(&self.buffer[self.pos..]);
        ordered.extend_from_slice(&self.buffer[..self.pos]);
        analyze(&ordered)
    }
}

pub struct BaselineTracker {
    window: SlidingWindow,
    baseline: Option<f64>,
    learning_samples: usize,
    samples_seen: usize,
}

impl BaselineTracker {
    pub fn new(window_size: usize, learning_samples: usize) -> Self {
        BaselineTracker {
            window: SlidingWindow::new(window_size),
            baseline: None,
            learning_samples,
            samples_seen: 0,
        }
    }

    pub fn push(&mut self, value: f64) -> Option<HealthVerdict> {
        self.window.push(value);
        self.samples_seen += 1;

        if !self.window.is_ready() {
            return None;
        }

        if self.samples_seen <= self.learning_samples {
            let law = self.window.analyze();
            if law.dfa.r_squared > 0.7 {
                self.baseline = Some(law.dfa.alpha);
            }
            return None;
        }

        let baseline = self.baseline?;
        let law = self.window.analyze();
        Some(health_check(&law, baseline))
    }

    pub fn baseline(&self) -> Option<f64> {
        self.baseline
    }

    pub fn is_learning(&self) -> bool {
        self.samples_seen <= self.learning_samples
    }
}

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

    fn white_noise(n: usize, seed: u64) -> Vec<f64> {
        let mut state = seed;
        (0..n).map(|_| {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            (state >> 33) as f64 / (1u64 << 31) as f64 - 0.5
        }).collect()
    }

    fn brownian(n: usize, seed: u64) -> Vec<f64> {
        let noise = white_noise(n, seed);
        let mut walk = Vec::with_capacity(n);
        let mut sum = 0.0;
        for v in noise {
            sum += v;
            walk.push(sum);
        }
        walk
    }

    #[test]
    fn dfa_short_separates_white_noise_from_random_walk_at_70_samples() {
        // Mean over seeds: white noise ~0.5, random walk ~1.5. The ranges are
        // wide because 70 samples is short; the point is that short series
        // get a real estimate where `dfa` returns the 0.5 placeholder.
        let (mut w, mut b) = (0.0, 0.0);
        for seed in 0..50 {
            w += dfa_short(&white_noise(70, seed)).expect("white").alpha;
            b += dfa_short(&brownian(70, seed)).expect("walk").alpha;
        }
        let (w, b) = (w / 50.0, b / 50.0);
        assert!((0.35..0.7).contains(&w), "white mean alpha {w}");
        assert!((1.25..1.75).contains(&b), "walk mean alpha {b}");
    }

    #[test]
    fn dfa_short_refuses_what_it_cannot_measure() {
        assert!(dfa_short(&white_noise(15, 1)).is_none(), "too short");
        assert!(dfa_short(&[3.0; 200]).is_none(), "constant");
        assert!(dfa_short(&white_noise(24, 1)).is_some(), "24 samples gives 3 box sizes");
    }

    /// Analytic standard deviation of the DFA α estimator at window n.
    ///
    /// Model: F²(s) averages W/s per-segment residual variances, each with
    /// s−2 degrees of freedom (linear detrend), so F² is chi-square-like
    /// with dof ≈ (W/s)(s−2) and var(ln F(s)) = ¼·var(ln F²) ≈
    /// 1 / (2·(W/s)·(s−2)). α is the OLS slope of ln F on ln s, hence
    /// var(α̂) = Σ (x_j − x̄)² v_j / (Σ (x_j − x̄)²)², x_j = ln s_j.
    /// The model treats box scales as independent, but they share the same
    /// profile, so this is a LOWER bound on the true variance.
    fn analytic_alpha_sd(n: usize) -> f64 {
        let s_min = 16usize.max(n / 50);
        let s_max = n / 4;
        let ratio = powf(s_max as f64 / s_min as f64, 1.0 / 11.0);
        let mut xs = Vec::new();
        let mut vs = Vec::new();
        let mut prev_s = 0usize;
        for step in 0..12 {
            let s = (s_min as f64 * powi(ratio, step)) as usize;
            if s == prev_s || s > s_max {
                continue;
            }
            prev_s = s;
            let num_segs = (n / s) as f64;
            xs.push(ln(s as f64));
            vs.push(1.0 / (2.0 * num_segs * (s as f64 - 2.0)));
        }
        let xbar = xs.iter().sum::<f64>() / xs.len() as f64;
        let sxx: f64 = xs.iter().map(|x| (x - xbar) * (x - xbar)).sum();
        let num: f64 = xs
            .iter()
            .zip(vs.iter())
            .map(|(x, v)| (x - xbar) * (x - xbar) * v)
            .sum();
        sqrt(num / (sxx * sxx))
    }

    #[test]
    fn analytic_alpha_sd_bounds_measured_scatter() {
        // Measured: sd of α over 800 independent white-noise windows.
        for &n in &[96usize, 192, 384] {
            let mut alphas = Vec::new();
            let mut buf = Vec::new();
            for seed in 0..800u64 {
                let w = white_noise(n, seed * 13 + 7);
                alphas.push(dfa_into(&w, &mut buf).alpha);
            }
            let mean = alphas.iter().sum::<f64>() / alphas.len() as f64;
            let var = alphas.iter().map(|a| (a - mean).powi(2)).sum::<f64>()
                / alphas.len() as f64;
            let measured = var.sqrt();
            let derived = analytic_alpha_sd(n);
            // The independence model is a LOWER bound: the box scales share
            // one profile, and that correlation inflates the true variance
            // by an n-dependent factor (measured: ~1.3x at n=96 rising to
            // ~4x at n=384, as more scales means more shared structure).
            // Assert the bound direction, and that the inflation stays below 6x
            // over the monitor's window range.
            let ratio = measured / derived;
            assert!(
                derived <= measured * 1.25,
                "n={}: derived {:.4} should not exceed measured {:.4}",
                n, derived, measured
            );
            assert!(
                ratio < 6.0,
                "n={}: inflation {:.1}x (derived {:.4}, measured {:.4})",
                n, ratio, derived, measured
            );
        }
    }

    #[test]
    fn dfa_fast_matches_dfa_into_exactly() {
        // 1000 random windows across lengths and signal classes:
        // prefix-sum DFA must agree with the reference at 1e-12.
        let mut buf_a = Vec::new();
        let mut buf_b = Vec::new();
        let mut worst = 0.0f64;
        for trial in 0..1000u64 {
            let n = 96 + (trial as usize * 37) % 417; // 96..512
            let data = if trial % 2 == 0 {
                white_noise(n, trial + 1)
            } else {
                brownian(n, trial + 1)
            };
            let a = dfa_into(&data, &mut buf_a);
            let b = dfa_fast_into(&data, &mut buf_b);
            let d = (a.alpha - b.alpha).abs();
            if d > worst {
                worst = d;
            }
            // 1e-9: prefix-difference Σy² reassociates floating-point ops;
            // Brownian-class signals (double-integrated profiles) cost a few
            // ulps. Far below any physically meaningful alpha difference.
            assert!(d < 1e-9, "trial {} n {} diff {}", trial, n, d);
            assert!((a.r_squared - b.r_squared).abs() < 1e-9);
        }
        assert!(worst < 1e-9, "worst diff {}", worst);
    }

    #[test]
    fn white_noise_alpha_near_half() {
        let data = white_noise(4096, 42);
        let result = dfa(&data);
        assert!(result.alpha > 0.35 && result.alpha < 0.65,
            "white noise DFA alpha should be near 0.5, got {}", result.alpha);
        assert!(result.r_squared > 0.8, "R2 should be high, got {}", result.r_squared);
    }

    #[test]
    fn brownian_alpha_above_one() {
        let data = brownian(4096, 42);
        let result = dfa(&data);
        assert!(result.alpha > 1.2 && result.alpha < 1.8,
            "brownian DFA alpha should be near 1.5, got {}", result.alpha);
    }

    #[test]
    fn deterministic() {
        let data = white_noise(1024, 7);
        let r1 = dfa(&data);
        let r2 = dfa(&data);
        assert!((r1.alpha - r2.alpha).abs() < 1e-10);
    }

    #[test]
    fn too_short_returns_half() {
        let data = [1.0; 10];
        let result = dfa(&data);
        assert_eq!(result.alpha, 0.5);
        assert_eq!(result.r_squared, 0.0);
    }

    #[test]
    fn analyze_produces_quality() {
        let data = white_noise(2048, 7);
        let law = analyze(&data);
        assert_eq!(law.n, 2048);
        assert!(law.quality != LawQuality::Insufficient);
    }

    #[test]
    fn health_verdict_thresholds() {
        assert_eq!(HealthVerdict::from_shift(0.01), HealthVerdict::Healthy);
        assert_eq!(HealthVerdict::from_shift(0.05), HealthVerdict::Watch);
        assert_eq!(HealthVerdict::from_shift(0.10), HealthVerdict::Warning);
        assert_eq!(HealthVerdict::from_shift(0.20), HealthVerdict::Critical);
        assert_eq!(HealthVerdict::from_shift(-0.20), HealthVerdict::Critical);
    }

    #[test]
    fn acr_detects_correlation() {
        let data = brownian(2048, 99);
        let result = acr(&data);
        assert!(result.alpha < -0.05, "brownian ACR exponent should be negative, got {}", result.alpha);
    }

    #[test]
    fn sliding_window_detects_after_fill() {
        let mut sw = SlidingWindow::new(256);
        assert!(!sw.is_ready());
        let noise = white_noise(256, 77);
        for v in &noise { sw.push(*v); }
        assert!(sw.is_ready());
        let law = sw.analyze();
        assert!(law.n == 256);
        assert!(law.dfa.alpha > 0.3);
    }

    #[test]
    fn baseline_tracker_learns_then_verdicts() {
        let mut bt = BaselineTracker::new(256, 500);
        let normal = brownian(600, 88);
        for (i, v) in normal.iter().enumerate() {
            let result = bt.push(*v);
            if i < 500 {
                assert!(result.is_none(), "should be learning at sample {}", i);
            }
        }
        assert!(!bt.is_learning());
    }

    #[test]
    fn sliding_window_before_fill_still_works() {
        let mut sw = SlidingWindow::new(512);
        for i in 0..100 {
            sw.push(i as f64 * 0.1);
        }
        assert!(!sw.is_ready());
        let law = sw.analyze();
        assert!(law.n == 100);
    }

    #[test]
    fn builtin_demo_data_detects_fault() {
        let normal: Vec<f64> = include_str!("../data/normal_sample.csv")
            .lines().filter_map(|l| l.trim().parse().ok()).collect();
        let fault: Vec<f64> = include_str!("../data/fault_sample.csv")
            .lines().filter_map(|l| l.trim().parse().ok()).collect();
        let law_n = analyze(&normal);
        let law_f = analyze(&fault);
        let verdict = health_check(&law_f, law_n.dfa.alpha);
        assert_eq!(verdict, HealthVerdict::Critical);
        assert!(law_n.dfa.r_squared > 0.9);
        assert!(law_f.dfa.r_squared > 0.9);
    }

    #[test]
    fn empty_input_does_not_panic() {
        let empty: Vec<f64> = vec![];
        let law = analyze(&empty);
        assert_eq!(law.quality, LawQuality::Insufficient);
        let result = dfa(&empty);
        assert_eq!(result.alpha, 0.5);
    }

    #[test]
    fn single_value_does_not_panic() {
        let law = analyze(&[42.0]);
        assert_eq!(law.quality, LawQuality::Insufficient);
    }

    #[test]
    fn all_nan_produces_abstain() {
        let nans = vec![f64::NAN; 100];
        let law = analyze(&nans);
        assert_eq!(law.quality, LawQuality::Insufficient);
    }

    #[test]
    fn inf_values_filtered() {
        let mut data = white_noise(256, 55);
        data[50] = f64::INFINITY;
        data[100] = f64::NEG_INFINITY;
        let law = analyze(&data);
        assert!(law.n < 256, "inf values should be filtered out");
    }

    #[test]
    fn constant_signal_abstains() {
        let constant = vec![core::f64::consts::PI; 200];
        let law = analyze(&constant);
        assert_eq!(law.quality, LawQuality::Abstain);
    }

    #[test]
    fn compare_identical_signals_healthy() {
        let data = white_noise(1024, 42);
        let result = compare(&data, &data);
        assert_eq!(result.verdict, HealthVerdict::Healthy);
        assert!(result.shift.abs() < 1e-10);
    }

    #[test]
    fn is_degraded_catches_structural_change() {
        let normal = white_noise(1024, 42);
        let brownian = brownian(1024, 42);
        assert!(is_degraded(&normal, &brownian));
    }

    #[test]
    fn has_changed_more_sensitive_than_is_degraded() {
        let data1 = white_noise(1024, 42);
        let data2 = white_noise(1024, 99);
        // Two different white noise samples should have similar alpha
        // but has_changed might catch tiny differences
        let _ = has_changed(&data1, &data2); // just verify no panic
    }

    #[test]
    fn white_noise_alpha_exact() {
        let data = white_noise(4096, 42);
        let result = dfa(&data);
        assert!((result.alpha - 0.5246244706).abs() < 1e-9,
            "expected alpha=0.5246244706, got {}", result.alpha);
        assert!((result.r_squared - 0.9853682683).abs() < 1e-9,
            "expected r_squared=0.9853682683, got {}", result.r_squared);
    }

    #[test]
    fn brownian_alpha_exact() {
        let data = brownian(4096, 42);
        let result = dfa(&data);
        assert!((result.alpha - 1.4180278792).abs() < 1e-9,
            "expected alpha=1.4180278792, got {}", result.alpha);
        assert!((result.r_squared - 0.9894368223).abs() < 1e-9,
            "expected r_squared=0.9894368223, got {}", result.r_squared);
    }

    #[test]
    fn health_verdict_exact_boundaries() {
        assert_eq!(HealthVerdict::from_shift(0.029999), HealthVerdict::Healthy);
        assert_eq!(HealthVerdict::from_shift(0.03), HealthVerdict::Watch);
        assert_eq!(HealthVerdict::from_shift(0.079999), HealthVerdict::Watch);
        assert_eq!(HealthVerdict::from_shift(0.08), HealthVerdict::Warning);
        assert_eq!(HealthVerdict::from_shift(0.149999), HealthVerdict::Warning);
        assert_eq!(HealthVerdict::from_shift(0.15), HealthVerdict::Critical);
        assert_eq!(HealthVerdict::from_shift(-0.03), HealthVerdict::Watch);
        assert_eq!(HealthVerdict::from_shift(-0.08), HealthVerdict::Warning);
        assert_eq!(HealthVerdict::from_shift(-0.15), HealthVerdict::Critical);
    }

    #[test]
    fn anomaly_scores_detects_shift() {
        let mut signal = white_noise(2048, 5);
        signal.extend(brownian(2048, 5));
        let scores = anomaly_scores(&signal, 256, 128, 0.05);
        assert!(!scores.is_empty(), "should produce per-window scores");
        let n = scores.len();
        let baseline_region = &scores[..n / 3];
        let shifted_region = &scores[2 * n / 3..];
        assert!(
            baseline_region.iter().all(|&s| s < 1.0),
            "baseline scores should stay under 1.0, got {baseline_region:?}"
        );
        assert!(
            shifted_region.iter().all(|&s| s > 1.0),
            "shifted-region scores should exceed 1.0, got {shifted_region:?}"
        );
    }

    #[test]
    fn anomaly_scores_too_short_returns_empty() {
        assert_eq!(anomaly_scores(&[1.0; 10], 256, 128, 0.05), Vec::<f64>::new());
        assert_eq!(anomaly_scores(&[1.0; 300], 32, 16, 0.05), Vec::<f64>::new());
    }

    #[test]
    fn anomaly_scores_window_boundary() {
        let data = white_noise(256, 99);
        let scores = anomaly_scores(&data, 64, 32, 0.05);
        assert!(!scores.is_empty(), "window=64 on 256 pts should produce scores");
        let scores_63 = anomaly_scores(&data, 63, 32, 0.05);
        assert!(scores_63.is_empty(), "window<64 should return empty");
    }

    #[test]
    fn anomaly_scores_baseline_uses_first_third() {
        let mut signal = white_noise(4096, 7);
        signal.extend(brownian(4096, 7));
        let scores = anomaly_scores(&signal, 256, 128, 0.05);
        let n = scores.len();
        let first_third_max = scores[..n/3].iter().cloned().fold(0.0f64, f64::max);
        let last_third_min = scores[2*n/3..].iter().cloned().fold(f64::MAX, f64::min);
        assert!(last_third_min > first_third_max,
            "shifted region should score higher than baseline: min={} vs max={}",
            last_third_min, first_third_max);
    }

    #[test]
    fn dfa_box_fit_operator_check() {
        let data = white_noise(512, 1);
        let r1 = dfa(&data);
        let data2 = brownian(512, 1);
        let r2 = dfa(&data2);
        assert!(r1.alpha < r2.alpha,
            "white noise alpha ({}) should be less than brownian ({})", r1.alpha, r2.alpha);
        assert!(r1.alpha > 0.3 && r1.alpha < 0.7, "white noise alpha out of range: {}", r1.alpha);
        assert!(r2.alpha > 1.0 && r2.alpha < 2.0, "brownian alpha out of range: {}", r2.alpha);
    }

    #[test]
    fn analyze_kurtosis_and_p99_computed() {
        let data = white_noise(1024, 42);
        let law = analyze(&data);
        assert!(law.kurtosis > 0.0, "kurtosis should be positive");
        assert!(law.p99 > law.mean, "p99 should exceed mean for noise");
        assert!(law.std_dev > 0.0, "std_dev should be positive for noise");
        assert!(law.max >= law.p99, "max should be >= p99");
    }
}

impl fmt::Display for LawQuality {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            LawQuality::Exact => write!(f, "EXACT"),
            LawQuality::Strong => write!(f, "STRONG"),
            LawQuality::Good => write!(f, "GOOD"),
            LawQuality::Approx => write!(f, "APPROX"),
            LawQuality::Abstain => write!(f, "ABSTAIN"),
            LawQuality::Insufficient => write!(f, "INSUFFICIENT"),
        }
    }
}

impl fmt::Display for HealthVerdict {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            HealthVerdict::Healthy => write!(f, "HEALTHY"),
            HealthVerdict::Watch => write!(f, "WATCH"),
            HealthVerdict::Warning => write!(f, "WARNING"),
            HealthVerdict::Critical => write!(f, "CRITICAL"),
        }
    }
}

impl fmt::Display for StructuralLaw {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "alpha={:.3} R2={:.4} H={:.3} quality={}", self.dfa.alpha, self.dfa.r_squared, self.hurst, self.quality)
    }
}

impl From<&[f64]> for SlidingWindow {
    fn from(data: &[f64]) -> Self {
        let mut sw = SlidingWindow::new(data.len().max(64));
        for &v in data { sw.push(v); }
        sw
    }
}

impl From<Vec<f64>> for SlidingWindow {
    fn from(data: Vec<f64>) -> Self {
        SlidingWindow::from(data.as_slice())
    }
}

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

impl Default for BaselineTracker {
    fn default() -> Self {
        BaselineTracker::new(256, 1000)
    }
}

impl HealthVerdict {
    pub fn from_shift_threshold(shift: f64, threshold: f64) -> Self {
        let s = if shift < 0.0 { -shift } else { shift };
        if s < threshold * 0.375 {
            HealthVerdict::Healthy
        } else if s < threshold {
            HealthVerdict::Watch
        } else if s < threshold * 1.875 {
            HealthVerdict::Warning
        } else {
            HealthVerdict::Critical
        }
    }
}

impl PartialEq for StructuralLaw {
    fn eq(&self, other: &Self) -> bool {
        self.quality == other.quality
            && (self.dfa.alpha - other.dfa.alpha).abs() < 1e-10
            && self.n == other.n
    }
}
// Simple API (start here)

/// Compare two signals and get a verdict: is the structure the same?
///
/// `baseline` is the known-good signal. `current` is what you're checking.
/// Returns a [`CompareResult`] with the verdict and the structural shift.
///
/// ```
/// use struktura::compare;
/// # let baseline = vec![1.0; 256];
/// # let current = vec![1.0; 256];
/// let result = compare(&baseline, &current);
/// println!("{}", result.verdict); // HEALTHY, WATCH, WARNING, or CRITICAL
/// ```
#[must_use]
pub fn compare(baseline: &[f64], current: &[f64]) -> CompareResult {
    let law_b = analyze(baseline);
    let law_c = analyze(current);
    let shift = law_c.dfa.alpha - law_b.dfa.alpha;
    let verdict = health_check(&law_c, law_b.dfa.alpha);
    CompareResult {
        baseline_alpha: law_b.dfa.alpha,
        current_alpha: law_c.dfa.alpha,
        shift,
        verdict,
        confidence: law_c.dfa.r_squared.min(law_b.dfa.r_squared),
    }
}

/// Is the current signal structurally degraded compared to baseline?
///
/// Returns `true` if the structural shift exceeds the Watch threshold (0.03).
/// For more detail, use [`compare`].
#[must_use]
pub fn is_degraded(baseline: &[f64], current: &[f64]) -> bool {
    let result = compare(baseline, current);
    result.verdict != HealthVerdict::Healthy
}

/// Has the signal's structure changed at all?
///
/// More sensitive than [`is_degraded`]: returns `true` on any measurable
/// shift (> 0.01), even if below the Watch threshold.
#[must_use]
pub fn has_changed(baseline: &[f64], current: &[f64]) -> bool {
    let result = compare(baseline, current);
    result.shift.abs() > 0.01 && result.confidence > 0.5
}

/// Result of comparing two signals.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct CompareResult {
    pub baseline_alpha: f64,
    pub current_alpha: f64,
    pub shift: f64,
    pub verdict: HealthVerdict,
    pub confidence: f64,
}

impl fmt::Display for CompareResult {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{} shift={:+.3} (baseline={:.3} current={:.3} R²={:.3})",
            self.verdict, self.shift, self.baseline_alpha, self.current_alpha, self.confidence)
    }
}

/// Per-window anomaly scores from sliding DFA.
///
/// Learns baseline α from the first `learn_windows` windows, then
/// scores each subsequent window as `|α_current - α_baseline| / threshold`.
/// Score > 1.0 means anomaly detected.
///
/// This matches the detector interface used in telemetry assurance
/// benchmarks: input signal → per-window anomaly score.
#[must_use]
pub fn anomaly_scores(values: &[f64], window: usize, step: usize, threshold: f64) -> Vec<f64> {
    if values.len() < window || window < 64 { return vec![]; }
    let mut alphas = Vec::new();
    let mut i = 0;
    while i + window <= values.len() {
        let w = &values[i..i + window];
        let d = dfa(w);
        alphas.push(d.alpha);
        i += step;
    }
    if alphas.is_empty() { return vec![]; }
    let learn_n = alphas.len() / 3; // use first third as baseline
    let learn_n = learn_n.max(3).min(alphas.len());
    let baseline: f64 = alphas[..learn_n].iter().sum::<f64>() / learn_n as f64;
    let var: f64 = alphas[..learn_n].iter().map(|a| powi(a - baseline, 2)).sum::<f64>() / learn_n as f64;
    let std = sqrt(var).max(threshold * 0.1);
    alphas.iter().map(|a| (a - baseline).abs() / (std + threshold)).collect()
}

// Domain modules

#[cfg(test)]
mod scratch_tests {
    use super::{dfa_into, dfa_scratch};

    /// `dfa_scratch` must reproduce `dfa_into` bit for bit: same profile,
    /// same boxes, same fit. A tolerance here would hide a real divergence.
    #[test]
    fn scratch_matches_into_bitwise() {
        let signals: [Vec<f64>; 3] = [
            (0..1024)
                .map(|i| ((i as f64 * 1103515245.0 + 12345.0) % 65536.0) / 65536.0 - 0.5)
                .collect(),
            (0..2048).map(|i| (i as f64 * 0.013).sin()).collect(),
            (0..4096)
                .map(|i| (i as f64 * 0.007).sin() + (i as f64 * 0.0003))
                .collect(),
        ];
        let mut buf = Vec::new();
        for sig in &signals {
            let mut scratch = vec![0.0f64; sig.len() + 7];
            let a = dfa_into(sig, &mut buf);
            let b = dfa_scratch(sig, &mut scratch);
            assert_eq!(a.alpha.to_bits(), b.alpha.to_bits(), "alpha {} vs {}", a.alpha, b.alpha);
            assert_eq!(a.r_squared.to_bits(), b.r_squared.to_bits());
        }
    }

    #[test]
    fn short_scratch_is_neutral() {
        let sig: Vec<f64> = (0..256).map(|i| (i as f64 * 0.1).sin()).collect();
        let mut short = [0.0f64; 255];
        let r = dfa_scratch(&sig, &mut short);
        assert_eq!(r.alpha, 0.5);
        assert_eq!(r.r_squared, 0.0);
    }
}

pub mod ffi;
pub mod space;
pub mod text;
pub mod market;
pub mod rhythm;
pub mod genome;
#[cfg(feature = "std")]
pub mod telemetry_bench;
pub mod monitor;
pub mod prognosis;
pub mod autopilot;
#[cfg(feature = "std")]
pub mod redblue;
#[cfg(feature = "std")]
pub mod evolve_real;
#[cfg(feature = "std")]
pub mod smap_eval;
#[cfg(feature = "std")]
pub mod rover;
pub mod rover_flight;
pub mod conformal;

/// Solve A x = b in place (Gauss–Jordan, partial pivoting). Shared by the
/// reconstruction and autoregression fitters. Returns false if singular.
pub(crate) fn solve_ridge(a: &mut [f64], b: &mut [f64], n: usize) -> bool {
    for col in 0..n {
        let mut pivot = col;
        for row in col + 1..n {
            if a[row * n + col].abs() > a[pivot * n + col].abs() {
                pivot = row;
            }
        }
        if a[pivot * n + col].abs() < 1e-12 {
            return false;
        }
        if pivot != col {
            for k in 0..n {
                a.swap(col * n + k, pivot * n + k);
            }
            b.swap(col, pivot);
        }
        let d = a[col * n + col];
        for k in 0..n {
            a[col * n + k] /= d;
        }
        b[col] /= d;
        for row in 0..n {
            if row != col {
                let f = a[row * n + col];
                if f != 0.0 {
                    for k in 0..n {
                        a[row * n + k] -= f * a[col * n + k];
                    }
                    b[row] -= f * b[col];
                }
            }
        }
    }
    true
}
pub mod mfdfa;
pub mod trend;
pub mod classify;
pub mod changepoint;
pub mod fingerprint;
#[cfg(feature = "std")]
pub mod codegen;

// Debugger modules (v1.8 telemetry investigation)
pub mod context;
pub mod incident;
#[cfg(feature = "std")]
pub mod case;
#[cfg(feature = "std")]
pub mod replay;
#[cfg(feature = "std")]
pub mod report;