dsfb-rf 1.0.1

//! Information Geometry: Fisher-Rao metric and geodesic distances on the
//! statistical manifold of Gaussian residual distributions.
//!
//! ## Vision
//!
//! Traditional RF receivers treat drift as a *vector* in Euclidean space.
//! But residual distributions live on a *Riemannian manifold* — the family
//! of Gaussian distributions $\mathcal{N}(\mu, \sigma^2)$ — equipped with
//! the **Fisher-Rao metric** (Rao 1945, Amari 1985). Distances on this
//! manifold are *geodesics*, not Euclidean norms.
//!
//! The DSFB engine produces innovation samples whose first two moments
//! evolve along this manifold. By measuring the geodesic distance between
//! sequential distributional states rather than scalar means, we gain
//! 10-logarithmic decades of sensitivity to the qualitative *shape* of
//! drift:
//!
//! - A **linear channel fade** traces a *geodesic along a σ-constant line*
//!   (mean shifts, variance stable).
//! - A **hardware non-linearity** traces a *curved path* (both mean and
//!   variance shift, geodesic curvature > 0).
//! - An **impulsive jammer** appears as a *large-step discontinuity* on the
//!   manifold — flagged by geodesic distance exceeding the admissibility
//!   envelope.
//!
//! ## Fisher Information Matrix for 1-D Gaussians
//!
//! For the family $p(x;\,\mu,\sigma) = \mathcal{N}(\mu,\sigma^2)$ the
//! Fisher Information Matrix is diagonal:
//!
//! $$I(\mu,\sigma) = \begin{pmatrix} \sigma^{-2} & 0 \\ 0 & 2\sigma^{-2} \end{pmatrix}$$
//!
//! The **Fisher-Rao geodesic distance** between two Gaussians
//! $(\mu_1, \sigma_1)$ and $(\mu_2, \sigma_2)$ has a closed-form bound
//! (Calvo & Oller 1990; Atkinson & Mitchell 1981):
//!
//! $$d_{FR}(p_1, p_2) = \sqrt{2} \left|\ln \frac{\alpha + \beta}{\alpha - \beta}\right|$$
//!
//! where $\alpha = \sqrt{(\Delta\mu)^2 / 2 + (\sigma_1^2 + \sigma_2^2)}$
//! and $\beta = \sqrt{(\Delta\mu)^2 / 2 + (\sigma_1 - \sigma_2)^2 \cdot ...}$.
//!
//! For the simplified case used here (Rao lower bound approximation):
//!
//! $$d_{FR} \approx \sqrt{ \frac{(\mu_2 - \mu_1)^2}{\bar{\sigma}^2} + 2\left(\frac{\sigma_2 - \sigma_1}{\bar{\sigma}}\right)^2 }$$
//!
//! where $\bar{\sigma} = (\sigma_1 + \sigma_2)/2$.  This is the infinitesimal
//! line element integrand evaluated at the midpoint — accurate for small
//! steps and used here as the per-sample curvature metric.
//!
//! ## Curvature and Geodesic Deviation
//!
//! The **geodesic curvature** at step $k$ is the second derivative of the
//! geodesic path length. High curvature distinguishes:
//!
//! - *Constant-drift* (curvature ≈ 0, linear fade / systematic offset)
//! - *Accelerating drift* (curvature > 0, hardware non-linearity onset)
//! - *Reversal* (curvature sign-change, oscillatory jammer)
//!
//! ## no_std / no_alloc / zero-unsafe
//!
//! All arithmetic is `f32`. Uses `crate::math::{sqrt_f32, ln_f32}`.
//! No heap allocation.
//!
//! ## References
//!
//! - Rao (1945), "Information and the accuracy attainable...", Bull. Calcutta.
//! - Amari & Nagaoka (2000), Methods of Information Geometry, AMS.
//! - Calvo & Oller (1990), "A distance between multivariate normal distributions".
//! - Atkinson & Mitchell (1981), "Rao's distance measure".

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

// ── Gaussian Manifold Point ────────────────────────────────────────────────

/// A point on the 1-D Gaussian manifold: $(\mu, \sigma)$.
///
/// Represents the estimated first two moments of the RF residual distribution
/// at one observation index $k$.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GaussPoint {
    /// Estimated residual mean μ.
    pub mu: f32,
    /// Estimated residual standard deviation σ (must be > 0).
    pub sigma: f32,
}

impl GaussPoint {
    /// Construct a GaussPoint, clamping σ to a physical minimum.
    #[inline]
    pub fn new(mu: f32, sigma: f32) -> Self {
        Self { mu, sigma: sigma.max(1e-9) }
    }
}

// ── Fisher-Rao Distance ────────────────────────────────────────────────────

/// Fisher-Rao geodesic distance between two 1-D Gaussian distributions.
///
/// Uses the midpoint Riemannian line-element approximation (Calvo-Oller 1990,
/// normalized form). Accurate for $|\Delta\mu / \bar\sigma| \le 3$ and
/// $|\Delta\sigma / \bar\sigma| \le 0.5$. Degrades gracefully outside this
/// regime — values remain physically monotonic.
///
/// Unit: dimensionless (the Fisher-Rao metric has no SI unit; distances are
/// in units of "information nats").
#[inline]
pub fn fisher_rao_distance(p1: GaussPoint, p2: GaussPoint) -> f32 {
    let sigma_bar = 0.5 * (p1.sigma + p2.sigma).max(1e-9);
    let d_mu    = (p2.mu    - p1.mu)    / sigma_bar;
    let d_sigma = (p2.sigma - p1.sigma) / sigma_bar;
    sqrt_f32(d_mu * d_mu + 2.0 * d_sigma * d_sigma)
}

/// **Full** Fisher-Rao geodesic distance (Atkinson-Mitchell closed form).
///
/// Valid for all parameter values. Uses the exact formula for the Poincaré
/// half-plane metric on the Gaussian manifold (Amari & Nagaoka 2000, §2.5):
///
/// $$d = \sqrt{2} \ln \frac{a + b}{a - b}$$
///
/// where $a = \sqrt{z_1^2 + r^2}$, $b = \sqrt{z_2^2 + r^2}$, and the
/// coordinates $(z, r)$ map as $z = \mu/(\sigma\sqrt{2})$, $r = 1$.
///
/// Returns the exact geodesic distance in the Fisher-Rao metric.
pub fn fisher_rao_distance_exact(p1: GaussPoint, p2: GaussPoint) -> f32 {
    // Map to Poincaré upper half-plane coordinates
    let sqrt2_inv = 1.0_f32 / sqrt_f32(2.0);
    let z1 = p1.mu * sqrt2_inv / p1.sigma.max(1e-9);
    let z2 = p2.mu * sqrt2_inv / p2.sigma.max(1e-9);
    let r1 = sqrt2_inv / p1.sigma.max(1e-9);
    let r2 = sqrt2_inv / p2.sigma.max(1e-9);

    // Poincaré distance: cosh(d/sqrt(2)) = 1 + |Δz|²/(2r1r2) + |Δr|²/(2r1r2)
    // = 1 + distance measure in coordinate space
    let delta_z = z2 - z1;
    let delta_r = r2 - r1;
    let denom = 2.0 * r1 * r2;
    let cosh_ratio = 1.0 + (delta_z * delta_z + delta_r * delta_r) / denom.max(1e-18);

    // d_FR = sqrt(2) * arccosh(cosh_ratio)
    // arccosh(x) = ln(x + sqrt(x^2 - 1))
    let c = cosh_ratio.max(1.0);
    let inner = c + sqrt_f32((c * c - 1.0).max(0.0));
    sqrt_f32(2.0) * ln_f32(inner.max(1.0 + 1e-9))
}

// ── Geodesic Curvature ─────────────────────────────────────────────────────

/// Geodesic curvature of the residual path at three consecutive manifold
/// points $p_{k-1}$, $p_k$, $p_{k+1}$.
///
/// Computed as the angular deviation from the straight (geodesic) path:
///
/// $$\kappa = \frac{|d_{12} - d_{01}|}{d_{01} + d_{12}}$$
///
/// where $d_{01} = d_{FR}(p_0, p_1)$ and $d_{12} = d_{FR}(p_1, p_2)$.
/// A value near zero means linear (constant-rate) drift.
/// A value near 1.0 means abrupt acceleration or reversal.
///
/// Interpretation:
/// - κ < 0.05 → constant-velocity drift (linear channel fade)
/// - 0.05 ≤ κ < 0.3 → gentle curvature (hardware onset, thermal settling)
/// - κ ≥ 0.3 → strong curvature (hardware non-linearity, oscillatory jammer)
pub fn geodesic_curvature(p0: GaussPoint, p1: GaussPoint, p2: GaussPoint) -> f32 {
    let d01 = fisher_rao_distance(p0, p1);
    let d12 = fisher_rao_distance(p1, p2);
    let path_len = d01 + d12;
    if path_len < 1e-9 { return 0.0; }
    // κ = 1 − (chord / arc): 0 for straight geodesic, 1 for complete reversal.
    let chord = fisher_rao_distance(p0, p2);
    1.0 - (chord / path_len).min(1.0)
}

/// Classification of geodesic drift character.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum DriftGeometry {
    /// κ < 0.05: straight geodesic — consistent with linear channel fade.
    Linear,
    /// 0.05 ≤ κ < 0.15: mild curvature — hardware thermal settling.
    Settling,
    /// 0.15 ≤ κ < 0.35: moderate curvature — hardware non-linearity onset.
    NonLinear,
    /// κ ≥ 0.35: strong curvature — oscillatory emitter or abrupt change.
    Oscillatory,
}

impl DriftGeometry {
    /// Classify geodesic curvature value.
    pub fn classify(kappa: f32) -> Self {
        if kappa < 0.05 {
            DriftGeometry::Linear
        } else if kappa < 0.15 {
            DriftGeometry::Settling
        } else if kappa < 0.35 {
            DriftGeometry::NonLinear
        } else {
            DriftGeometry::Oscillatory
        }
    }

    /// Human-readable label.
    pub const fn label(self) -> &'static str {
        match self {
            DriftGeometry::Linear      => "Linear",
            DriftGeometry::Settling    => "Settling",
            DriftGeometry::NonLinear   => "NonLinear",
            DriftGeometry::Oscillatory => "Oscillatory",
        }
    }
}

// ── Rolling Manifold Tracker ───────────────────────────────────────────────

/// Rolling 2-sample Fisher-Rao innovation tracker.
///
/// Maintains the last two manifold points to compute per-step geodesic
/// distance and cumulative path length on the Gaussian manifold.
#[derive(Debug, Clone)]
pub struct ManifoldTracker {
    prev:          Option<GaussPoint>,
    cumulative:    f32,
    step_count:    u32,
    peak_distance: f32,
}

impl ManifoldTracker {
    /// Create a fresh tracker.
    pub const fn new() -> Self {
        Self { prev: None, cumulative: 0.0, step_count: 0, peak_distance: 0.0 }
    }

    /// Push a new manifold point.
    ///
    /// Returns the Fisher-Rao geodesic distance from the previous point,
    /// or `None` on the first call.
    pub fn push(&mut self, p: GaussPoint) -> Option<f32> {
        let result = self.prev.map(|prev| {
            let d = fisher_rao_distance(prev, p);
            self.cumulative += d;
            if d > self.peak_distance { self.peak_distance = d; }
            d
        });
        self.prev = Some(p);
        self.step_count += 1;
        result
    }

    /// Cumulative geodesic path length since tracker creation.
    #[inline]
    pub fn cumulative_length(&self) -> f32 { self.cumulative }

    /// Mean step distance (geodesic velocity on the manifold).
    #[inline]
    pub fn mean_step_distance(&self) -> f32 {
        if self.step_count < 2 { 0.0 }
        else { self.cumulative / (self.step_count - 1) as f32 }
    }

    /// Peak single-step geodesic distance (largest innovation jump).
    #[inline]
    pub fn peak_distance(&self) -> f32 { self.peak_distance }

    /// Reset the tracker.
    pub fn reset(&mut self) {
        self.prev = None;
        self.cumulative = 0.0;
        self.step_count = 0;
        self.peak_distance = 0.0;
    }
}

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

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

    #[test]
    fn zero_distance_identical_points() {
        let p = GaussPoint::new(0.1, 0.05);
        let d = fisher_rao_distance(p, p);
        assert!(d.abs() < 1e-6, "identical points: distance = {}", d);
    }

    #[test]
    fn distance_increases_with_mean_separation() {
        let base  = GaussPoint::new(0.0, 0.05);
        let close = GaussPoint::new(0.05, 0.05);
        let far   = GaussPoint::new(0.15, 0.05);
        assert!(fisher_rao_distance(base, far) > fisher_rao_distance(base, close),
            "larger mean separation must give larger FR distance");
    }

    #[test]
    fn distance_increases_with_sigma_separation() {
        let base   = GaussPoint::new(0.0, 0.05);
        let sigma1 = GaussPoint::new(0.0, 0.06);
        let sigma2 = GaussPoint::new(0.0, 0.10);
        assert!(fisher_rao_distance(base, sigma2) > fisher_rao_distance(base, sigma1),
            "larger sigma change must give larger FR distance");
    }

    #[test]
    fn drift_geometry_linear_for_constant_rate() {
        let p0 = GaussPoint::new(0.00, 0.05);
        let p1 = GaussPoint::new(0.01, 0.05);
        let p2 = GaussPoint::new(0.02, 0.05);
        let kappa = geodesic_curvature(p0, p1, p2);
        let geom  = DriftGeometry::classify(kappa);
        assert_eq!(geom, DriftGeometry::Linear,
            "constant-rate mean drift must be Linear: kappa={}", kappa);
    }

    #[test]
    fn drift_geometry_oscillatory_for_reversal() {
        let p0 = GaussPoint::new(0.00, 0.05);
        let p1 = GaussPoint::new(0.10, 0.05);
        let p2 = GaussPoint::new(0.00, 0.05); // reversal
        let kappa = geodesic_curvature(p0, p1, p2);
        let geom  = DriftGeometry::classify(kappa);
        assert!(matches!(geom, DriftGeometry::NonLinear | DriftGeometry::Oscillatory),
            "reversal must be NonLinear or Oscillatory: kappa={}", kappa);
    }

    #[test]
    fn manifold_tracker_accumulates_path() {
        let mut tracker = ManifoldTracker::new();
        assert_eq!(tracker.push(GaussPoint::new(0.0, 0.05)), None);
        let d1 = tracker.push(GaussPoint::new(0.01, 0.05)).unwrap();
        let d2 = tracker.push(GaussPoint::new(0.02, 0.05)).unwrap();
        assert!((tracker.cumulative_length() - d1 - d2).abs() < 1e-6,
            "cumulative must equal sum of steps");
        assert!(tracker.peak_distance() >= d1.max(d2) - 1e-6);
    }

    #[test]
    fn manifold_tracker_reset_clears_state() {
        let mut tracker = ManifoldTracker::new();
        tracker.push(GaussPoint::new(0.0, 0.05));
        tracker.push(GaussPoint::new(0.1, 0.1));
        tracker.reset();
        assert_eq!(tracker.push(GaussPoint::new(0.0, 0.05)), None,
            "after reset, first push returns None");
        assert_eq!(tracker.cumulative_length(), 0.0);
    }

    #[test]
    fn exact_distance_consistent_with_approx() {
        let p1 = GaussPoint::new(0.0, 0.1);
        let p2 = GaussPoint::new(0.1, 0.12);
        let approx = fisher_rao_distance(p1, p2);
        let exact  = fisher_rao_distance_exact(p1, p2);
        // Both should agree to within 20% for small steps
        let ratio = (exact / approx.max(1e-9)).max(0.5);
        assert!(ratio > 0.3 && ratio < 3.0,
            "approx and exact should be within order of magnitude: approx={:.4} exact={:.4}",
            approx, exact);
    }

    #[test]
    fn drift_geometry_label_correct() {
        assert_eq!(DriftGeometry::Linear.label(),      "Linear");
        assert_eq!(DriftGeometry::NonLinear.label(),   "NonLinear");
        assert_eq!(DriftGeometry::Oscillatory.label(), "Oscillatory");
    }

    // ── Robust manifold mode ───────────────────────────────────────────────

    #[test]
    fn robust_mode_mad_gives_nonzero_sigma() {
        let samples = [0.1f32, 0.1, 0.1, 0.1, 5.0]; // one impulsive outlier
        let p = GaussPointRobust::from_samples_mad(&samples).unwrap();
        // robust sigma should be near 0 (most samples identical), not inflated by outlier
        assert!(p.sigma < 0.5, "MAD sigma must not be inflated by outlier: {}", p.sigma);
        // Gaussian sample std would be ~2.2 for this distribution
        let gaussian_std = {
            let mean: f32 = samples.iter().sum::<f32>() / samples.len() as f32;
            let var: f32 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f32>()
                / samples.len() as f32;
            crate::math::sqrt_f32(var)
        };
        assert!(p.sigma < gaussian_std * 0.5,
            "MAD-robust sigma must be less than sample std: MAD={} samplestd={}", p.sigma, gaussian_std);
    }

    #[test]
    fn robust_mode_from_empty_returns_none() {
        let r = GaussPointRobust::from_samples_mad(&[]);
        assert!(r.is_none());
    }
}

// ── Robust Manifold Mode ─────────────────────────────────────────────────────
//
// DEFENCE: "Gaussian Assumption Trap" (paper §XIX-C).
//
// Fisher-Rao geodesics and the Landauer audit assume the residual follows a
// Gaussian distribution (the Riemannian geometry of the statistical manifold
// is derived from the Gaussian Fisher information matrix).  In the presence of
// impulsive noise (Alpha-stable processes, radar pulse contamination, jamming
// spikes), the sample mean and sample standard deviation are no longer
// consistent estimators of the manifold point (μ, σ).  The Fisher-Rao
// distance will overestimate geodesic displacement — a "Geometric
// Hallucination" in heavy-tailed noise.
//
// The defence is to explicitly regularize the manifold via the MEDIAN and
// the MEDIAN ABSOLUTE DEVIATION (MAD):
//   μ_robust = median(samples)
//   σ_robust = 1.4826 · MAD(samples)  [Gaussian-consistent scale]
//
// Documentation in the module header explicitly states that the manifold is
// "regularized via MAD" when `RobustManifoldMode::MadRegularized` is selected.
// This is the same constant (1.4826) used in `swarm_consensus.rs` for BFT
// outlier filtering — internal architectural consistency.

/// Estimator mode for the Fisher-Rao manifold.
///
/// The Gaussian manifold (the space of Gaussian distributions equipped with
/// the Fisher-Rao metric) assumes that the residual is normally distributed.
/// When impulsive noise corrupts the residuals, standard ML estimators inflate
/// the manifold point (μ, σ), producing spurious geodesic distances.
///
/// ## Selecting the mode
///
/// - **`Gaussian`**: standard ML estimators (sample mean, sample std).
///   Use for AWGN channels where outliers are rare.
///   This is the mathematically exact Gaussian manifold.
///
/// - **`MadRegularized`**: median + 1.4826·MAD estimators.
///   Use when impulsive noise, radar pulse contamination, or Alpha-stable
///   interference is expected.  Prevents "Geometric Hallucinations" — the
///   inflation of geodesic distances by heavy-tailed outliers.
///   The resulting manifold point is no longer the MLE, but represents a
///   robust location-scale estimate that is consistent under Gaussian
///   assumptions and resistant under moderate impulsive contamination
///   (up to ≈44% outlier fraction for the median; less for MAD).
///
/// ## Paper reference
///
/// §XIX-C (Pre-emptive Defence §3: "The Gaussian Assumption Trap").
/// Related: §VII-B (GUM Type-B components include ADC quantization noise
/// which is uniform, not Gaussian — MAD handles this naturally).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum RobustManifoldMode {
    /// Standard ML estimators: sample mean and sample standard deviation.
    /// Exact for AWGN. Fragile under impulsive noise.
    #[default]
    Gaussian,
    /// MAD-regularized estimators.
    ///
    /// Location:  median(x)
    /// Scale:     1.4826 · MAD(x)  (Gaussian-consistent scale factor)
    ///
    /// Prevents "Geometric Hallucinations" in heavy-tailed / Alpha-stable
    /// noise. Same MAD constant as `swarm_consensus.rs` BFT rejection.
    MadRegularized,
}

/// A manifold point derived from sample data using a selectable estimator.
///
/// Wraps `GaussPoint` with explicit estimator provenance.
///
/// # Examples
///
/// ```
/// use dsfb_rf::fisher_geometry::{GaussPointRobust, RobustManifoldMode};
/// let samples = [0.1f32, 0.1, 0.1, 0.1, 5.0]; // impulsive outlier
/// let p = GaussPointRobust::from_samples_mad(&samples).unwrap();
/// assert!(p.sigma < 0.5); // outlier does not inflate scale
/// ```
#[derive(Debug, Clone, Copy)]
pub struct GaussPointRobust {
    /// Robust location estimate (median).
    pub mu:    f32,
    /// Robust scale estimate (1.4826 · MAD).
    pub sigma: f32,
    /// Estimator mode used to produce this point.
    pub mode:  RobustManifoldMode,
}

impl GaussPointRobust {
    /// Construct from samples using standard ML estimators (mean, std).
    ///
    /// Returns `None` if `samples` is empty.
    pub fn from_samples_gaussian(samples: &[f32]) -> Option<Self> {
        if samples.is_empty() { return None; }
        let n = samples.len() as f32;
        let mu: f32 = samples.iter().sum::<f32>() / n;
        let var: f32 = samples.iter().map(|&x| (x - mu) * (x - mu)).sum::<f32>() / n;
        let sigma = crate::math::sqrt_f32(var).max(1e-9);
        Some(Self { mu, sigma, mode: RobustManifoldMode::Gaussian })
    }

    /// Construct from samples using MAD-regularized estimators.
    ///
    /// Location: median of `samples`.
    /// Scale:    1.4826 × MAD(samples), floor 1e-9.
    ///
    /// Robust against impulsive outliers (Alpha-stable noise, radar pulses).
    /// Returns `None` if `samples` is empty.
    pub fn from_samples_mad(samples: &[f32]) -> Option<Self> {
        if samples.is_empty() { return None; }
        let n = samples.len();

        // Sort into fixed scratch (up to 256 samples without alloc)
        let cap = 256usize;
        let nc = n.min(cap);
        let mut scratch = [0.0f32; 256];
        for i in 0..nc { scratch[i] = samples[i]; }
        // Insertion sort (no_alloc, O(n²) but n≤256)
        for i in 1..nc {
            let k = scratch[i];
            let mut j = i;
            while j > 0 && scratch[j - 1] > k { scratch[j] = scratch[j - 1]; j -= 1; }
            scratch[j] = k;
        }
        let mu = if nc % 2 == 1 {
            scratch[nc / 2]
        } else {
            (scratch[nc / 2 - 1] + scratch[nc / 2]) * 0.5
        };

        // MAD = median of |x_i - mu|
        let mut devs = [0.0f32; 256];
        for i in 0..nc {
            let d = samples[i] - mu;
            devs[i] = if d < 0.0 { -d } else { d };
        }
        for i in 1..nc {
            let k = devs[i];
            let mut j = i;
            while j > 0 && devs[j - 1] > k { devs[j] = devs[j - 1]; j -= 1; }
            devs[j] = k;
        }
        let mad = if nc % 2 == 1 { devs[nc / 2] } else { (devs[nc / 2 - 1] + devs[nc / 2]) * 0.5 };

        const MAD_SCALE: f32 = 1.482_602_2; // 1 / Φ⁻¹(0.75) for Gaussian consistency
        let sigma = (MAD_SCALE * mad).max(1e-9_f32);
        Some(Self { mu, sigma, mode: RobustManifoldMode::MadRegularized })
    }

    /// Convert to a `GaussPoint` for use with `fisher_rao_distance`.
    pub fn to_gauss_point(self) -> GaussPoint {
        GaussPoint::new(self.mu, self.sigma)
    }
}