dsfb-rf 1.0.1

//! Delay-coordinate attractor embedding and Koopman operator proxy.
//!
//! ## Theoretical Basis
//!
//! **Takens' Embedding Theorem (1981):** A scalar time-series {x(t)} generated by
//! a d-dimensional dynamical system can be reconstructed (diffeomorphically) in
//! ℝ^m using delay coordinates [x(t), x(t−τ), x(t−2τ), …, x(t−(m−1)τ)] when
//! m ≥ 2d+1 (Whitney's theorem).  For RF residual norms (d ≤ 2 in practice),
//! m = 2, τ = 2 suffices.
//!
//! **Correlation Dimension D₂ (Grassberger & Procaccia 1983):**
//! ```text
//!   C(r) = lim_{N→∞} (2 / N(N-1)) · #{pairs (i,j) : ||xᵢ − xⱼ|| < r}
//!   D₂   = lim_{r→0} log C(r) / log r
//! ```
//! Pure noise produces D₂ → m (fills the embedding space).
//! A low-dimensional attractor produces D₂ ≪ m.
//!
//! **Koopman Operator Proxy (Mezić 2005):** The spectral analysis of the linear
//! Koopman operator K acting on observables f(x) → f(F(x)) is approximated here
//! by the variance-to-mean ratio of the embedding trajectory — a computationally
//! tractable proxy that distinguishes stochastic diffusion (high VM) from
//! structured orbits (low VM, persistent modes).
//!
//! ## Design
//!
//! - `no_std`, `no_alloc`, zero `unsafe`
//! - Fixed-capacity circular buffer W for delay-coordinate history
//! - O(W²) correlation dimension per query — calibration-path only (W ≤ 64)
//!
//! ## References
//!
//! Takens, F. (1981) "Detecting strange attractors in turbulence," *Lecture Notes
//!   in Mathematics* 898:366–381. Springer. doi:10.1007/BFb0091924.
//!
//! Grassberger, P. and Procaccia, I. (1983) "Measuring the strangeness of strange
//!   attractors," *Physica D* 9:189–208. doi:10.1016/0167-2789(83)90298-1.
//!
//! Mezić, I. (2005) "Spectral properties of dynamical systems, model reduction and
//!   decompositions," *Nonlinear Dynamics* 41:309–325. doi:10.1007/s11071-005-2824-x.

use crate::math::{mean_f32, sqrt_f32};

// ── Attractor State ────────────────────────────────────────────────────────

/// Classification of the reconstructed residual attractor.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum NoiseAttractorState {
    /// Residuals fill the embedding space: purely stochastic, no hidden structure.
    /// D₂ ≈ embedding dimension.  Koopman modes absent.
    StochasticBall,
    /// Residuals form a bounded, recurring orbit: hidden determinism present.
    /// D₂ < embedding dimension.  Persistent Koopman modes detected.
    StructuredOrbit,
    /// Residuals have collapsed to a lower-dimensional manifold: systematic drift
    /// or near-periodic interference pattern dominating the residual process.
    CollapsedAttractor,
}

/// Output of an attractor analysis pass.
#[derive(Debug, Clone, Copy)]
pub struct AttractorResult {
    /// Grassberger-Procaccia correlation dimension estimate D₂.
    /// Range [0, 2] for 2-dimensional delay embedding.
    pub correlation_dimension: f32,
    /// Koopman proxy: variance-to-mean ratio of embedding norms.
    /// High (> 1.0) → stochastic; Low (< 0.3) → structured modes.
    pub koopman_proxy: f32,
    /// Inferred attractor state classification.
    pub state: NoiseAttractorState,
    /// Number of delay-coordinate pairs used for D₂ estimation.
    pub n_pairs: usize,
}

// ── Delay Embedding ────────────────────────────────────────────────────────

/// Delay-coordinate embedding estimator.
///
/// Maintains a circular buffer of W residual norms and computes the
/// 2D delay embedding [(x(t), x(t−τ))] across the window.
///
/// ## Type Parameters
/// - `W`: Embedding window size (≥ 4, ≤ 64 recommended).
pub struct DelayEmbedding<const W: usize> {
    /// Circular buffer of residual norms.
    buf: [f32; W],
    /// Write head.
    head: usize,
    /// Number of samples absorbed (saturates at W).
    count: usize,
    /// Delay lag τ.
    pub tau: usize,
}

impl<const W: usize> DelayEmbedding<W> {
    /// Create a new estimator with delay lag τ.
    pub const fn new(tau: usize) -> Self {
        Self { buf: [0.0; W], head: 0, count: 0, tau }
    }

    /// Absorb one residual norm sample.
    pub fn push(&mut self, norm: f32) {
        self.buf[self.head] = norm;
        self.head = (self.head + 1) % W;
        if self.count < W { self.count += 1; }
    }

    /// Number of valid samples in the buffer.
    pub fn len(&self) -> usize { self.count }

    /// Whether the buffer has enough samples for analysis.
    pub fn is_ready(&self) -> bool { self.count >= W }

    /// Reset the buffer.
    pub fn reset(&mut self) {
        self.buf = [0.0; W];
        self.head = 0;
        self.count = 0;
    }

    /// Access sample at position `i` relative to the most-recent sample
    /// (0 = latest, 1 = one step back, etc.).
    fn sample(&self, i: usize) -> f32 {
        let idx = (self.head + W - 1 - i) % W;
        self.buf[idx]
    }

    /// Compute attractor analysis over the current window.
    ///
    /// Returns `None` if the buffer contains fewer than `tau + 4` valid samples.
    pub fn analyse(&self, r1: f32, r2: f32) -> Option<AttractorResult> {
        if self.count < self.tau + 4 { return None; }
        let n = self.count.min(W);
        let usable = if n > self.tau { n - self.tau } else { return None; };
        if usable < 3 { return None; }
        let n_pairs = usable;

        let (c_r1, c_r2, total_pairs) = self.count_pairs_in_radii(n_pairs, r1, r2);
        let d2 = correlation_dimension_gp(c_r1, c_r2, total_pairs, r1, r2);
        let koopman_proxy = self.koopman_variance_to_mean(n);
        let state = classify_attractor_state(d2, koopman_proxy);

        Some(AttractorResult {
            correlation_dimension: d2,
            koopman_proxy,
            state,
            n_pairs,
        })
    }

    fn count_pairs_in_radii(&self, n_pairs: usize, r1: f32, r2: f32) -> (u32, u32, u32) {
        let mut c_r1: u32 = 0;
        let mut c_r2: u32 = 0;
        let mut total_pairs: u32 = 0;
        for i in 0..n_pairs {
            let xi = self.sample(i);
            let xi_d = self.sample(i + self.tau);
            for j in (i + 1)..n_pairs {
                let xj = self.sample(j);
                let xj_d = self.sample(j + self.tau);
                let dx = xi - xj;
                let dy = xi_d - xj_d;
                let dist = sqrt_f32(dx * dx + dy * dy);
                if dist < r1 { c_r1 += 1; c_r2 += 1; }
                else if dist < r2 { c_r2 += 1; }
                total_pairs += 1;
            }
        }
        (c_r1, c_r2, total_pairs)
    }

    fn koopman_variance_to_mean(&self, n: usize) -> f32 {
        let mut norms_buf = [0.0_f32; W];
        for i in 0..n {
            norms_buf[i] = self.sample(i);
        }
        let slice = &norms_buf[..n];
        let m = mean_f32(slice);
        let var = slice.iter().map(|&x| (x - m) * (x - m)).sum::<f32>() / n as f32;
        if m > 1e-6 { var / m } else { 0.0 }
    }
}

fn correlation_dimension_gp(c_r1: u32, c_r2: u32, total_pairs: u32, r1: f32, r2: f32) -> f32 {
    if total_pairs == 0 { return 2.0; }
    if c_r1 == total_pairs { return 0.0; }
    if c_r1 < 2 || c_r2 < 2 || r1 <= 0.0 || r2 <= r1 { return 2.0; }
    let fr1 = c_r1 as f32 / total_pairs as f32;
    let fr2 = c_r2 as f32 / total_pairs as f32;
    if fr1 <= 0.0 || fr2 <= fr1 { return 2.0; }
    let log_ratio_c = log2_approx(fr2 / fr1);
    let log_ratio_r = log2_approx(r2 / r1);
    (log_ratio_c / log_ratio_r).max(0.0).min(4.0)
}

fn classify_attractor_state(d2: f32, koopman_proxy: f32) -> NoiseAttractorState {
    if d2 < 0.6 {
        NoiseAttractorState::CollapsedAttractor
    } else if d2 < 1.4 && koopman_proxy < 0.5 {
        NoiseAttractorState::StructuredOrbit
    } else {
        NoiseAttractorState::StochasticBall
    }
}

impl<const W: usize> Default for DelayEmbedding<W> {
    fn default() -> Self { Self::new(2) }
}

// ── Private helpers ────────────────────────────────────────────────────────

/// log₂(x) approximation via identity log₂(x) = log(x) / log(2).
///
/// Uses the identity log₂(x) = ln(x) / ln(2) with ln(x) approximated by
/// a 4-term Padé expansion around x=1: ln(x) ≈ 2·(u + u³/3 + u⁵/5) where u=(x-1)/(x+1).
/// Accurate to < 0.5% for x ∈ [0.01, 100].
fn log2_approx(x: f32) -> f32 {
    if x <= 0.0 { return -100.0; }
    // Find integer power of 2 to reduce to [1, 2). f32 exponent range
    // is bounded by [−149, 127]; 300 iterations covers the full range.
    let mut y = x;
    let mut exp = 0i32;
    for _ in 0..300 {
        if y < 2.0 { break; }
        y *= 0.5; exp += 1;
    }
    for _ in 0..300 {
        if y >= 1.0 { break; }
        y *= 2.0; exp -= 1;
    }
    // y ∈ [1, 2); compute log₂(y) via ln Padé + convert
    let u = (y - 1.0) / (y + 1.0);
    let u2 = u * u;
    let ln_y = 2.0 * u * (1.0 + u2 / 3.0 + u2 * u2 / 5.0);
    let ln2 = 0.693_147_f32;
    exp as f32 + ln_y / ln2
}

// ── Tests ──────────────────────────────────────────────────────────────────

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

    #[test]
    fn push_fills_buffer() {
        let mut emb = DelayEmbedding::<8>::new(2);
        assert!(!emb.is_ready());
        for i in 0..8 { emb.push(i as f32 * 0.01); }
        assert!(emb.is_ready());
    }

    #[test]
    fn reset_clears_buffer() {
        let mut emb = DelayEmbedding::<8>::new(2);
        for _ in 0..8 { emb.push(1.0); }
        emb.reset();
        assert_eq!(emb.len(), 0);
        assert!(!emb.is_ready());
    }

    #[test]
    fn constant_signal_collapsed_attractor() {
        let mut emb = DelayEmbedding::<32>::new(2);
        for _ in 0..32 { emb.push(0.1_f32); }
        // Constant → all pairs at distance 0 → very high C(r1) → D₂ → 0 → Collapsed
        let res = emb.analyse(0.05, 0.2).unwrap();
        assert_eq!(res.state, NoiseAttractorState::CollapsedAttractor,
            "constant signal: D₂={:.3}", res.correlation_dimension);
    }

    #[test]
    fn stochastic_signal_fills_space() {
        // Widely varying norms span the embedding space: D₂ should be higher
        let mut emb = DelayEmbedding::<32>::new(2);
        // LCG pseudo-random
        let mut s: u32 = 12345;
        for _ in 0..32 {
            s = s.wrapping_mul(1664525).wrapping_add(1013904223);
            let v = (s >> 16) as f32 / 65535.0;
            emb.push(v);
        }
        let res = emb.analyse(0.05, 0.3).unwrap();
        // Stochastic → D₂ tends toward 2.0; at minimum, state should not be Collapsed
        assert_ne!(res.state, NoiseAttractorState::CollapsedAttractor,
            "random norms: D₂={:.3}", res.correlation_dimension);
    }

    #[test]
    fn not_enough_samples_returns_none() {
        let emb = DelayEmbedding::<32>::new(2);
        assert!(emb.analyse(0.05, 0.2).is_none());
    }

    #[test]
    fn log2_approx_accuracy() {
        assert!((log2_approx(1.0) - 0.0).abs() < 0.01, "log₂(1)=0");
        assert!((log2_approx(2.0) - 1.0).abs() < 0.02, "log₂(2)=1");
        assert!((log2_approx(4.0) - 2.0).abs() < 0.02, "log₂(4)=2");
        assert!((log2_approx(0.5) + 1.0).abs() < 0.02, "log₂(0.5)=-1");
    }
}