math-sonify 1.4.0

//! Fractal Dimension Analyzer — characterises the geometry of strange attractors.
//!
//! This module provides several complementary methods for quantifying the
//! fractal structure of attractor trajectories generated by the ODE systems in
//! [`crate::systems`]:
//!
//! | Algorithm | Struct | Notes |
//! |-----------|--------|-------|
//! | Box-counting | [`BoxCounting`] | 2-D projection; log-log slope of N(ε) |
//! | Correlation dimension | [`CorrelationDimension`] | Grassberger–Procaccia C(r) ~ r^D₂ |
//! | Lyapunov spectrum | [`LyapunovSpectrum`] | Full QR-based spectrum; Kaplan–Yorke dim |
//! | Attractor summary | [`AttractorCharacterization`] | Aggregates all metrics |
//!
//! All algorithms are designed to run on trajectory snapshots collected by the
//! simulation thread and can be polled from the MATH VIEW tab in the GUI.

#![allow(dead_code)]

// ---------------------------------------------------------------------------
// Box-counting dimension
// ---------------------------------------------------------------------------

/// Computes the box-counting (Minkowski–Bouligand) dimension of a 2-D
/// attractor projection.
///
/// The algorithm:
/// 1. Project the trajectory to a 2-D plane (any pair of state variables).
/// 2. For each of several box sizes ε, count how many ε×ε boxes the
///    trajectory visits (N(ε)).
/// 3. Fit a straight line to log(N) vs log(1/ε) — the slope is D₀.
///
/// A Lorenz attractor gives D₀ ≈ 2.05; the Hénon map gives ≈ 1.26.
pub struct BoxCounting {
    /// Box sizes to evaluate (log-spaced).  Smaller values give finer
    /// resolution but require more trajectory points to populate.
    pub box_sizes: Vec<f64>,
    /// Most recently computed dimension estimate.
    pub dimension: f64,
    /// log(1/ε) values from the most recent call to [`Self::compute`].
    pub log_inv_eps: Vec<f64>,
    /// log(N(ε)) values from the most recent call to [`Self::compute`].
    pub log_count: Vec<f64>,
}

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

impl BoxCounting {
    /// Create a box-counter with 10 log-spaced box sizes between ε=0.01 and ε=2.0.
    pub fn new() -> Self {
        let box_sizes = log_space(0.01, 2.0, 10);
        Self {
            box_sizes,
            dimension: 0.0,
            log_inv_eps: Vec::new(),
            log_count: Vec::new(),
        }
    }

    /// Compute the box-counting dimension from a 2-D trajectory slice.
    ///
    /// * `points` — `(x, y)` pairs sampled from the attractor.
    ///
    /// Returns the estimated dimension D₀ ∈ [0, 2].
    pub fn compute(&mut self, points: &[(f64, f64)]) -> f64 {
        if points.len() < 8 {
            self.dimension = 0.0;
            return 0.0;
        }

        let (x_min, x_max, y_min, y_max) = bounding_box_2d(points);
        let x_range = (x_max - x_min).max(f64::EPSILON);
        let y_range = (y_max - y_min).max(f64::EPSILON);

        let mut log_inv: Vec<f64> = Vec::with_capacity(self.box_sizes.len());
        let mut log_n: Vec<f64> = Vec::with_capacity(self.box_sizes.len());

        for &eps in &self.box_sizes {
            // Normalise box size to the data range.
            let nx = ((x_range / eps).ceil() as usize).min(4096).max(1);
            let ny = ((y_range / eps).ceil() as usize).min(4096).max(1);

            // Use a HashSet of (ix, iy) to count occupied boxes.
            let mut occupied: std::collections::HashSet<(u32, u32)> =
                std::collections::HashSet::new();
            for &(x, y) in points {
                if x.is_finite() && y.is_finite() {
                    let ix = (((x - x_min) / x_range) * nx as f64) as u32;
                    let iy = (((y - y_min) / y_range) * ny as f64) as u32;
                    occupied.insert((ix.min(nx as u32 - 1), iy.min(ny as u32 - 1)));
                }
            }

            let count = occupied.len();
            if count > 0 {
                log_inv.push((1.0 / eps).ln());
                log_n.push((count as f64).ln());
            }
        }

        self.log_inv_eps = log_inv.clone();
        self.log_count = log_n.clone();

        if log_inv.len() < 2 {
            self.dimension = 1.0;
            return 1.0;
        }

        let d = linear_regression_slope(&log_inv, &log_n).clamp(0.0, 3.0);
        self.dimension = d;
        d
    }
}

// ---------------------------------------------------------------------------
// Correlation dimension (Grassberger–Procaccia)
// ---------------------------------------------------------------------------

/// Estimates the correlation dimension D₂ of an attractor using the
/// Grassberger–Procaccia algorithm.
///
/// The correlation integral C(r) counts the fraction of point pairs within
/// distance r:
///
///   C(r) = (2 / N(N-1)) · #{(i,j) : i<j, ‖xᵢ − xⱼ‖ < r}
///
/// For a fractal attractor C(r) ~ r^D₂ in the scaling regime, so the slope
/// of log C vs log r gives D₂.
pub struct CorrelationDimension {
    /// Number of distance thresholds evaluated.
    pub num_radii: usize,
    /// Most recently computed dimension estimate D₂.
    pub dimension: f64,
    /// log(r) values from the most recent computation.
    pub log_r: Vec<f64>,
    /// log(C(r)) values from the most recent computation.
    pub log_c: Vec<f64>,
}

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

impl CorrelationDimension {
    pub fn new(num_radii: usize) -> Self {
        Self {
            num_radii: num_radii.max(4),
            dimension: 0.0,
            log_r: Vec::new(),
            log_c: Vec::new(),
        }
    }

    /// Estimate D₂ from a set of state-space points.
    ///
    /// * `points` — each inner slice is one state-space point; all must have
    ///   the same length.
    /// * `max_pairs` — cap on the number of pair comparisons (performance).
    ///   Set to e.g. 50_000 for real-time use.
    ///
    /// Returns the estimated correlation dimension D₂.
    pub fn compute(&mut self, points: &[Vec<f64>], max_pairs: usize) -> f64 {
        if points.len() < 8 {
            self.dimension = 0.0;
            return 0.0;
        }

        // Compute pairwise distances (capped).
        let n = points.len();
        let mut distances: Vec<f64> = Vec::new();
        let mut count = 0usize;
        'outer: for i in 0..n {
            for j in (i + 1)..n {
                let d = euclidean_distance(&points[i], &points[j]);
                if d.is_finite() {
                    distances.push(d);
                }
                count += 1;
                if count >= max_pairs {
                    break 'outer;
                }
            }
        }

        if distances.is_empty() {
            self.dimension = 0.0;
            return 0.0;
        }

        distances.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let d_min = *distances.first().unwrap_or(&f64::EPSILON);
        let d_max = *distances.last().unwrap_or(&1.0);
        if d_max <= d_min {
            self.dimension = 1.0;
            return 1.0;
        }

        let radii = log_space(d_min.max(f64::EPSILON), d_max, self.num_radii);
        let total_pairs = distances.len() as f64;

        let mut log_r: Vec<f64> = Vec::with_capacity(radii.len());
        let mut log_c: Vec<f64> = Vec::with_capacity(radii.len());

        for &r in &radii {
            // Use binary search to count pairs within radius r.
            let count_r = distances.partition_point(|&d| d < r) as f64;
            let c = count_r / total_pairs;
            if c > 0.0 {
                log_r.push(r.ln());
                log_c.push(c.ln());
            }
        }

        self.log_r = log_r.clone();
        self.log_c = log_c.clone();

        if log_r.len() < 2 {
            self.dimension = 1.0;
            return 1.0;
        }

        let d = linear_regression_slope(&log_r, &log_c).clamp(0.0, 10.0);
        self.dimension = d;
        d
    }
}

// ---------------------------------------------------------------------------
// Lyapunov spectrum
// ---------------------------------------------------------------------------

/// Computes the full Lyapunov exponent spectrum using the QR decomposition
/// (Benettin et al.) method.
///
/// The algorithm evolves N tangent vectors simultaneously under the linearised
/// flow (Jacobian), re-orthonormalises them via Gram–Schmidt / QR each step,
/// and accumulates the log-stretching rates.
///
/// For a 3-D system the spectrum has three exponents λ₁ ≥ λ₂ ≥ λ₃.
/// Their sum equals the divergence (trace of Jacobian = volume contraction rate).
/// The Kaplan–Yorke dimension D_KY is also computed.
pub struct LyapunovSpectrum {
    /// Accumulated log-sums for each exponent.
    log_growth: Vec<f64>,
    /// Number of steps accumulated since last reset.
    step_count: u64,
    /// Current spectrum (λ₁ ≥ … ≥ λₙ), updated each call to [`Self::update`].
    pub spectrum: Vec<f64>,
    /// Kaplan–Yorke dimension derived from the spectrum.
    pub kaplan_yorke_dim: f64,
    /// Running sum of the spectrum (= Kolmogorov entropy proxy, positive part).
    pub kolmogorov_entropy: f64,
    /// Phase-space volume contraction rate = sum of all exponents.
    pub volume_contraction: f64,
}

impl LyapunovSpectrum {
    /// Create a new spectrum estimator for a system of dimension `dim`.
    pub fn new(dim: usize) -> Self {
        Self {
            log_growth: vec![0.0; dim.max(1)],
            step_count: 0,
            spectrum: vec![0.0; dim.max(1)],
            kaplan_yorke_dim: dim as f64,
            kolmogorov_entropy: 0.0,
            volume_contraction: 0.0,
        }
    }

    /// Dimension of the system.
    pub fn dim(&self) -> usize {
        self.log_growth.len()
    }

    /// Reset all accumulators (call after a system change).
    pub fn reset(&mut self) {
        for v in &mut self.log_growth {
            *v = 0.0;
        }
        self.step_count = 0;
    }

    /// Update the spectrum from a trajectory `points` (each is a full state vector).
    ///
    /// Uses a finite-difference approximation of the Jacobian evaluated along
    /// the trajectory.  For real-time use keep `points` short (e.g. 64 steps).
    ///
    /// * `system` — must implement [`crate::systems::DynamicalSystem`].
    /// * `dt` — integration step size.
    /// * `eps` — finite-difference perturbation size.
    pub fn update_from_trajectory(
        &mut self,
        points: &[Vec<f64>],
        dt: f64,
        deriv_fn: &impl Fn(&[f64]) -> Vec<f64>,
    ) {
        let n = self.dim();
        if points.len() < 2 || n == 0 {
            return;
        }

        let eps = 1e-7;

        // We initialise an orthonormal frame.
        let mut frame: Vec<Vec<f64>> = (0..n)
            .map(|i| {
                let mut v = vec![0.0f64; n];
                if i < n { v[i] = 1.0; }
                v
            })
            .collect();

        for state in points {
            if state.len() < n {
                continue;
            }

            // Approximate action of the Jacobian on each basis vector via
            // finite differences: J·v ≈ (f(x + ε·v) − f(x − ε·v)) / (2ε).
            let f0 = deriv_fn(state);
            let mut new_frame: Vec<Vec<f64>> = Vec::with_capacity(n);

            for v in &frame {
                let x_plus: Vec<f64> = state.iter().zip(v.iter()).map(|(&s, &vi)| s + eps * vi).collect();
                let x_minus: Vec<f64> = state.iter().zip(v.iter()).map(|(&s, &vi)| s - eps * vi).collect();
                let f_plus = deriv_fn(&x_plus);
                let f_minus = deriv_fn(&x_minus);

                let jv: Vec<f64> = f_plus
                    .iter()
                    .zip(f_minus.iter())
                    .map(|(&fp, &fm)| (fp - fm) / (2.0 * eps))
                    .collect();

                // Euler step: new_v = v + dt * J·v.
                let new_v: Vec<f64> = v.iter().zip(jv.iter()).map(|(&vi, &jvi)| vi + dt * jvi).collect();
                new_frame.push(new_v);
            }

            // QR via modified Gram–Schmidt.
            let norms = gram_schmidt_inplace(&mut new_frame);

            // Accumulate log-norms.
            for (i, norm) in norms.iter().enumerate() {
                if *norm > 0.0 {
                    self.log_growth[i] += norm.ln();
                }
            }
            self.step_count += 1;
            frame = new_frame;
            let _ = f0; // suppress unused warning
        }

        self.recompute_from_accumulators(dt);
    }

    /// Recompute public fields from the current accumulators.
    fn recompute_from_accumulators(&mut self, dt: f64) {
        if self.step_count == 0 || dt <= 0.0 {
            return;
        }
        let t_total = self.step_count as f64 * dt;
        for (i, &lg) in self.log_growth.iter().enumerate() {
            self.spectrum[i] = lg / t_total;
        }
        // Sort descending.
        self.spectrum.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));

        // Kaplan–Yorke dimension.
        self.kaplan_yorke_dim = kaplan_yorke(&self.spectrum);

        // Kolmogorov entropy (sum of positive exponents).
        self.kolmogorov_entropy = self.spectrum.iter().filter(|&&x| x > 0.0).sum::<f64>();

        // Volume contraction rate.
        self.volume_contraction = self.spectrum.iter().sum::<f64>();
    }

    /// Human-readable summary of the spectrum.
    pub fn summary(&self) -> String {
        let parts: Vec<String> = self
            .spectrum
            .iter()
            .map(|&x| format!("{:+.4}", x))
            .collect();
        format!(
            "λ = [{}]  D_KY={:.3}  hKS={:.4}  div={:.4}",
            parts.join(", "),
            self.kaplan_yorke_dim,
            self.kolmogorov_entropy,
            self.volume_contraction,
        )
    }
}

// ---------------------------------------------------------------------------
// Attractor characterization
// ---------------------------------------------------------------------------

/// Complete characterisation of an attractor's geometric and dynamical
/// properties.  Updated periodically by the simulation thread and displayed
/// in the MATH VIEW tab.
#[derive(Debug, Clone)]
pub struct AttractorCharacterization {
    /// Box-counting dimension D₀ ∈ [0, 3].
    pub fractal_dim: f64,
    /// Correlation dimension D₂ ∈ [0, 3].
    pub correlation_dim: f64,
    /// Full Lyapunov spectrum (λ₁ ≥ λ₂ ≥ … ≥ λₙ), all in s⁻¹.
    pub lyapunov_spectrum: Vec<f64>,
    /// Kaplan–Yorke dimension derived from the spectrum.
    pub kaplan_yorke_dim: f64,
    /// Kolmogorov–Sinai entropy proxy = sum of positive Lyapunov exponents.
    pub kolmogorov_entropy: f64,
    /// Phase-space volume contraction rate = sum of all Lyapunov exponents.
    pub phase_space_volume_contraction: f64,
    /// Qualitative attractor type label.
    pub attractor_type: AttractorType,
    /// Number of trajectory points used for the most recent analysis.
    pub sample_size: usize,
    /// Timestamp (monotonic ticks) of the last update.
    pub last_updated_ticks: u64,
}

/// Qualitative attractor type determined from the Lyapunov spectrum.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum AttractorType {
    /// All Lyapunov exponents negative — stable fixed point.
    FixedPoint,
    /// One zero exponent, all others negative — limit cycle.
    LimitCycle,
    /// Two zero exponents — quasi-periodic torus (T²).
    QuasiPeriodicTorus,
    /// Exactly one positive exponent — strange attractor / deterministic chaos.
    StrangeAttractor,
    /// Two or more positive exponents — hyperchaos.
    Hyperchaos,
    /// Spectrum not yet computed.
    Unknown,
}

impl AttractorType {
    /// Derive the attractor type from a Lyapunov spectrum (must be sorted descending).
    pub fn from_spectrum(spectrum: &[f64]) -> Self {
        if spectrum.is_empty() {
            return Self::Unknown;
        }
        let tol = 0.01;
        let positive = spectrum.iter().filter(|&&x| x > tol).count();
        let zero = spectrum.iter().filter(|&&x| x.abs() <= tol).count();

        match (positive, zero) {
            (0, 0) => Self::FixedPoint,
            (0, 1) => Self::LimitCycle,
            (0, 2) => Self::QuasiPeriodicTorus,
            (1, _) => Self::StrangeAttractor,
            _ => Self::Hyperchaos,
        }
    }

    pub fn label(&self) -> &'static str {
        match self {
            Self::FixedPoint => "Fixed Point",
            Self::LimitCycle => "Limit Cycle",
            Self::QuasiPeriodicTorus => "Quasi-Periodic (T²)",
            Self::StrangeAttractor => "Strange Attractor",
            Self::Hyperchaos => "Hyperchaos",
            Self::Unknown => "Unknown",
        }
    }
}

impl Default for AttractorCharacterization {
    fn default() -> Self {
        Self {
            fractal_dim: 0.0,
            correlation_dim: 0.0,
            lyapunov_spectrum: Vec::new(),
            kaplan_yorke_dim: 0.0,
            kolmogorov_entropy: 0.0,
            phase_space_volume_contraction: 0.0,
            attractor_type: AttractorType::Unknown,
            sample_size: 0,
            last_updated_ticks: 0,
        }
    }
}

/// Convenience builder that runs all analyses on a trajectory snapshot.
pub struct FractalAnalyzer {
    box_counting: BoxCounting,
    correlation: CorrelationDimension,
    lyapunov: LyapunovSpectrum,
}

impl FractalAnalyzer {
    /// Create an analyzer for an N-dimensional system.
    pub fn new(dim: usize) -> Self {
        Self {
            box_counting: BoxCounting::new(),
            correlation: CorrelationDimension::new(16),
            lyapunov: LyapunovSpectrum::new(dim),
        }
    }

    /// Run all analyses on `trajectory` (each entry is a full state vector).
    ///
    /// * `deriv_fn` — derivative function for the Lyapunov computation.
    /// * `dt` — ODE step size.
    /// * `current_tick` — monotonic tick counter for the timestamp field.
    pub fn analyze(
        &mut self,
        trajectory: &[Vec<f64>],
        deriv_fn: &impl Fn(&[f64]) -> Vec<f64>,
        dt: f64,
        current_tick: u64,
    ) -> AttractorCharacterization {
        if trajectory.is_empty() {
            return AttractorCharacterization::default();
        }

        // 2-D projection (XY) for box-counting.
        let pts_2d: Vec<(f64, f64)> = trajectory
            .iter()
            .map(|s| (s.first().copied().unwrap_or(0.0), s.get(1).copied().unwrap_or(0.0)))
            .collect();
        let fractal_dim = self.box_counting.compute(&pts_2d);

        // Correlation dimension (cap pairs to keep it real-time friendly).
        let correlation_dim = self.correlation.compute(trajectory, 20_000);

        // Lyapunov spectrum.
        self.lyapunov.update_from_trajectory(trajectory, dt, deriv_fn);
        let spectrum = self.lyapunov.spectrum.clone();
        let kaplan_yorke_dim = self.lyapunov.kaplan_yorke_dim;
        let kolmogorov_entropy = self.lyapunov.kolmogorov_entropy;
        let volume_contraction = self.lyapunov.volume_contraction;
        let attractor_type = AttractorType::from_spectrum(&spectrum);

        AttractorCharacterization {
            fractal_dim,
            correlation_dim,
            lyapunov_spectrum: spectrum,
            kaplan_yorke_dim,
            kolmogorov_entropy,
            phase_space_volume_contraction: volume_contraction,
            attractor_type,
            sample_size: trajectory.len(),
            last_updated_ticks: current_tick,
        }
    }

    /// Access the Lyapunov spectrum directly (e.g. for streaming display).
    pub fn lyapunov_spectrum(&self) -> &LyapunovSpectrum {
        &self.lyapunov
    }

    /// Reset the Lyapunov accumulator (e.g. after a system change).
    pub fn reset_lyapunov(&mut self) {
        self.lyapunov.reset();
    }
}

// ---------------------------------------------------------------------------
// Internal helpers
// ---------------------------------------------------------------------------

/// Gram–Schmidt orthonormalisation in place.  Returns the norms before
/// normalisation (used for Lyapunov exponent accumulation).
fn gram_schmidt_inplace(vecs: &mut Vec<Vec<f64>>) -> Vec<f64> {
    let n = vecs.len();
    let mut norms = vec![1.0f64; n];
    for i in 0..n {
        // Subtract projections onto already orthonormal vectors.
        for j in 0..i {
            let proj = dot(&vecs[i], &vecs[j]);
            let vj = vecs[j].clone();
            let vi = &mut vecs[i];
            for (a, b) in vi.iter_mut().zip(vj.iter()) {
                *a -= proj * b;
            }
        }
        let norm = vec_norm(&vecs[i]);
        norms[i] = norm;
        if norm > f64::EPSILON {
            for x in &mut vecs[i] {
                *x /= norm;
            }
        }
    }
    norms
}

fn dot(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum()
}

fn vec_norm(v: &[f64]) -> f64 {
    v.iter().map(|x| x * x).sum::<f64>().sqrt()
}

fn euclidean_distance(a: &[f64], b: &[f64]) -> f64 {
    a.iter()
        .zip(b.iter())
        .map(|(&x, &y)| (x - y).powi(2))
        .sum::<f64>()
        .sqrt()
}

/// Generate `n` log-evenly-spaced values from `lo` to `hi`.
fn log_space(lo: f64, hi: f64, n: usize) -> Vec<f64> {
    if n < 2 {
        return vec![lo];
    }
    let log_lo = lo.max(f64::EPSILON).ln();
    let log_hi = hi.max(f64::EPSILON).ln();
    (0..n)
        .map(|i| (log_lo + (log_hi - log_lo) * i as f64 / (n - 1) as f64).exp())
        .collect()
}

/// Ordinary least-squares slope of y ~ slope * x.
fn linear_regression_slope(x: &[f64], y: &[f64]) -> f64 {
    let n = x.len().min(y.len());
    if n < 2 {
        return 0.0;
    }
    let xm: f64 = x[..n].iter().sum::<f64>() / n as f64;
    let ym: f64 = y[..n].iter().sum::<f64>() / n as f64;
    let num: f64 = x[..n].iter().zip(y[..n].iter()).map(|(&xi, &yi)| (xi - xm) * (yi - ym)).sum();
    let den: f64 = x[..n].iter().map(|&xi| (xi - xm).powi(2)).sum();
    if den.abs() < f64::EPSILON { 0.0 } else { num / den }
}

/// Kaplan–Yorke dimension from a sorted (descending) Lyapunov spectrum.
fn kaplan_yorke(spectrum: &[f64]) -> f64 {
    if spectrum.is_empty() {
        return 0.0;
    }
    let mut sum = 0.0f64;
    let mut j = 0usize;
    for (i, &l) in spectrum.iter().enumerate() {
        sum += l;
        if sum < 0.0 {
            j = i;
            break;
        }
        j = i + 1;
    }
    if j == 0 {
        return 0.0;
    }
    // D_KY = j + (sum_{i=0}^{j-1} λᵢ) / |λⱼ|
    let partial_sum: f64 = spectrum[..j].iter().sum();
    let lambda_j = spectrum.get(j).copied().unwrap_or(f64::NEG_INFINITY);
    if lambda_j >= 0.0 {
        // All exponents non-negative: KY = full dimension.
        return spectrum.len() as f64;
    }
    let d = j as f64 + partial_sum / (-lambda_j).max(f64::EPSILON);
    d.clamp(0.0, spectrum.len() as f64)
}

fn bounding_box_2d(pts: &[(f64, f64)]) -> (f64, f64, f64, f64) {
    let mut xmin = f64::MAX;
    let mut xmax = f64::MIN;
    let mut ymin = f64::MAX;
    let mut ymax = f64::MIN;
    for &(x, y) in pts {
        if x < xmin { xmin = x; }
        if x > xmax { xmax = x; }
        if y < ymin { ymin = y; }
        if y > ymax { ymax = y; }
    }
    (xmin, xmax, ymin, ymax)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    // ---- log_space ----------------------------------------------------------

    #[test]
    fn log_space_length() {
        let v = log_space(0.01, 2.0, 10);
        assert_eq!(v.len(), 10);
    }

    #[test]
    fn log_space_bounds() {
        let v = log_space(0.1, 10.0, 8);
        assert!((v[0] - 0.1).abs() < 1e-6);
        assert!((v[7] - 10.0).abs() < 1e-6);
    }

    // ---- linear_regression_slope --------------------------------------------

    #[test]
    fn linear_regression_unit_slope() {
        let x: Vec<f64> = (0..10).map(|i| i as f64).collect();
        let y: Vec<f64> = x.clone();
        let s = linear_regression_slope(&x, &y);
        assert!((s - 1.0).abs() < 1e-10, "slope should be 1, got {s}");
    }

    #[test]
    fn linear_regression_zero_slope() {
        let x: Vec<f64> = (0..10).map(|i| i as f64).collect();
        let y: Vec<f64> = vec![3.0; 10];
        let s = linear_regression_slope(&x, &y);
        assert!(s.abs() < 1e-10, "slope should be 0, got {s}");
    }

    // ---- kaplan_yorke -------------------------------------------------------

    #[test]
    fn kaplan_yorke_lorenz_like() {
        // Lorenz exponents ≈ +0.9, 0, -14.5 → D_KY ≈ 2.06
        let spectrum = vec![0.9, 0.0, -14.5];
        let dky = kaplan_yorke(&spectrum);
        assert!(dky > 2.0 && dky < 2.2, "D_KY={dky} should be ≈2.06");
    }

    #[test]
    fn kaplan_yorke_all_negative() {
        let spectrum = vec![-1.0, -2.0, -3.0];
        let dky = kaplan_yorke(&spectrum);
        assert_eq!(dky, 0.0);
    }

    // ---- AttractorType ------------------------------------------------------

    #[test]
    fn attractor_type_fixed_point() {
        assert_eq!(
            AttractorType::from_spectrum(&[-1.0, -2.0, -3.0]),
            AttractorType::FixedPoint
        );
    }

    #[test]
    fn attractor_type_limit_cycle() {
        assert_eq!(
            AttractorType::from_spectrum(&[-0.001, -1.0, -2.0]),
            AttractorType::LimitCycle
        );
    }

    #[test]
    fn attractor_type_strange_attractor() {
        assert_eq!(
            AttractorType::from_spectrum(&[0.9, 0.0, -14.5]),
            AttractorType::StrangeAttractor
        );
    }

    #[test]
    fn attractor_type_hyperchaos() {
        assert_eq!(
            AttractorType::from_spectrum(&[0.5, 0.3, 0.0, -5.0]),
            AttractorType::Hyperchaos
        );
    }

    // ---- BoxCounting --------------------------------------------------------

    #[test]
    fn box_counting_circle_dim_near_one() {
        // A circle in 2D has box-counting dimension 1.
        let n = 400usize;
        let pts: Vec<(f64, f64)> = (0..n)
            .map(|i| {
                let theta = 2.0 * std::f64::consts::PI * i as f64 / n as f64;
                (theta.cos(), theta.sin())
            })
            .collect();
        let mut bc = BoxCounting::new();
        let d = bc.compute(&pts);
        // Box-counting of a circle ≈ 1; allow generous tolerance for finite samples.
        assert!(d > 0.8 && d < 1.3, "circle dimension should be ≈1, got {d}");
    }

    #[test]
    fn box_counting_empty_returns_zero() {
        let mut bc = BoxCounting::new();
        assert_eq!(bc.compute(&[]), 0.0);
    }

    // ---- CorrelationDimension -----------------------------------------------

    #[test]
    fn correlation_dimension_line_near_one() {
        // Points on a line → D₂ ≈ 1.
        let pts: Vec<Vec<f64>> = (0..200)
            .map(|i| vec![i as f64 * 0.05, 0.0, 0.0])
            .collect();
        let mut cd = CorrelationDimension::new(12);
        let d = cd.compute(&pts, 5000);
        assert!(d > 0.7 && d < 1.4, "line D₂ should be ≈1, got {d}");
    }

    #[test]
    fn correlation_dimension_too_few_points() {
        let pts: Vec<Vec<f64>> = vec![vec![1.0, 0.0], vec![2.0, 0.0]];
        let mut cd = CorrelationDimension::default();
        assert_eq!(cd.compute(&pts, 1000), 0.0);
    }

    // ---- LyapunovSpectrum ---------------------------------------------------

    #[test]
    fn lyapunov_spectrum_lorenz_positive_largest() {
        // Lorenz linearized: check that the largest exponent is positive.
        let sigma = 10.0f64;
        let rho = 28.0f64;
        let beta = 8.0 / 3.0;

        let deriv = |s: &[f64]| {
            if s.len() < 3 {
                return vec![0.0; s.len()];
            }
            vec![
                sigma * (s[1] - s[0]),
                s[0] * (rho - s[2]) - s[1],
                s[0] * s[1] - beta * s[2],
            ]
        };

        // Build a short trajectory.
        let mut state = vec![1.0f64, 0.0, 0.0];
        let dt = 0.01;
        let mut traj = Vec::new();
        for _ in 0..100 {
            traj.push(state.clone());
            let d = deriv(&state);
            for (s, ds) in state.iter_mut().zip(d.iter()) {
                *s += dt * ds;
            }
        }

        let mut spec = LyapunovSpectrum::new(3);
        spec.update_from_trajectory(&traj, dt, &deriv);

        // After a short run the signs may not be settled, but at minimum it
        // should not panic and the spectrum should have 3 entries.
        assert_eq!(spec.spectrum.len(), 3);
        for &l in &spec.spectrum {
            assert!(l.is_finite(), "spectrum entry {l} is not finite");
        }
    }

    #[test]
    fn lyapunov_spectrum_reset_clears_accumulators() {
        let mut spec = LyapunovSpectrum::new(3);
        spec.log_growth[0] = 999.0;
        spec.step_count = 100;
        spec.reset();
        assert_eq!(spec.step_count, 0);
        assert!(spec.log_growth.iter().all(|&v| v == 0.0));
    }

    // ---- gram_schmidt -------------------------------------------------------

    #[test]
    fn gram_schmidt_orthonormal() {
        let mut vecs = vec![
            vec![3.0f64, 1.0, 0.0],
            vec![2.0, 2.0, 0.0],
            vec![0.0, 0.0, 5.0],
        ];
        gram_schmidt_inplace(&mut vecs);
        // Check orthogonality.
        let tol = 1e-10;
        for i in 0..3 {
            for j in 0..3 {
                let d = dot(&vecs[i], &vecs[j]);
                if i == j {
                    assert!((d - 1.0).abs() < tol, "v{i}·v{i}={d}, expected 1");
                } else {
                    assert!(d.abs() < tol, "v{i}·v{j}={d}, expected 0");
                }
            }
        }
    }

    // ---- FractalAnalyzer ----------------------------------------------------

    #[test]
    fn fractal_analyzer_returns_characterization() {
        let lorenz_deriv = |s: &[f64]| {
            if s.len() < 3 { return vec![0.0f64; s.len()]; }
            let (sigma, rho, beta) = (10.0f64, 28.0, 8.0 / 3.0);
            vec![sigma * (s[1] - s[0]), s[0] * (rho - s[2]) - s[1], s[0] * s[1] - beta * s[2]]
        };

        let mut state = vec![1.0f64, 0.0, 0.0];
        let dt = 0.01;
        let mut traj = Vec::new();
        for _ in 0..200 {
            traj.push(state.clone());
            let d = lorenz_deriv(&state);
            for (s, ds) in state.iter_mut().zip(d.iter()) {
                *s += dt * ds;
            }
        }

        let mut analyzer = FractalAnalyzer::new(3);
        let ch = analyzer.analyze(&traj, &lorenz_deriv, dt, 1000);

        assert_eq!(ch.sample_size, 200);
        assert_eq!(ch.last_updated_ticks, 1000);
        assert!(ch.fractal_dim >= 0.0);
        assert!(ch.correlation_dim >= 0.0);
        assert!(ch.kaplan_yorke_dim >= 0.0);
    }
}