opendeviationbar-core 13.78.0

//! Bar-close-rolling features: operate over `&[f64]` of bar closes.
//!
//! Distinct from `tier2.rs` / `tier3.rs`, which operate over
//! `&[&TradeSnapshot]`: these features take N consecutive bar-close scalars as
//! input (one per prior bar), not the underlying trade-level snapshots.
//!
//! The two fractal-dimension functions (`compute_bar_petrosian_fd` /
//! `compute_bar_katz_fd`) are bit-faithful ports of the canonical `antropy` 0.2.2
//! implementations (`antropy/fractal.py`), validated by an exact-oracle unit test
//! (`assert ~1e-9` vs `antropy.petrosian_fd` / `antropy.katz_fd`) over a committed
//! real-BTCUSDT bar-close fixture. `compute_bar_dispersion_entropy` is a separate
//! Rostaghi & Azami (2016) port (NCDF dispersion entropy), validated by its own
//! 1e-9 oracle vs `entropy_candidates._dispersion_entropy`.
//!
//! Robustness contract — the two FRACTAL-DIMENSION functions (petrosian/katz):
//!   * NEVER return NaN or +/-Inf. antropy itself returns NaN on a flat window
//!     (Katz `0/0`) and +/-Inf when `d == a` (Katz `log10(1)` denominator); a
//!     leading finite-input guard plus a trailing `is_finite` guard make the
//!     output domain total.
//!   * Return the antropy value wherever antropy is finite (so a flat window
//!     yields Petrosian's natural `1.0`), and the sentinel `0.0` ONLY where the
//!     input is undefined (`len < 2`, non-finite) or where antropy would emit
//!     NaN/Inf (the Katz degenerate cases). `0.0` is disjoint from the valid
//!     output domain (a real fractal dimension is `>= 1` on bar-close data), so
//!     callers may treat `0.0` as "not applicable".
//!
//! Robustness contract — `compute_bar_dispersion_entropy` uses the INVERSE
//! convention: it returns `NaN` (NOT `0.0`) for an undefined window (`len < 2`,
//! non-finite input, flat `σ == 0`, or `< 2` embedding vectors), because `0.0` is
//! a VALID dispersion value (a single repeated ordinal pattern). The wiring layer
//! maps that `NaN → None/NULL` and never coerces `0.0`. It still NEVER returns
//! `±Inf` (proven by the `dispersion_class_*` / overflow unit tests below).
//!
//! `libm::log10` (the project-standard pure-Rust libm, also used by `tier2`/
//! `tier3`) is used rather than `f64::log10` so the column is bit-deterministic
//! across platforms; the resulting ~1 ULP difference vs numpy/glibc `log10` is
//! absorbed by the oracle test's tolerance.

use libm;

use crate::constants::entropy::{ORDINAL_M, ORDINAL_M_FACTORIAL, ORDINAL_TAU};
use crate::interbar_math::tier3::ordinal_pattern_index_m3;

// FILE-SIZE-OK: this is the contract-designated single home for ALL bar_close
// feature kernels (petrosian / katz / dispersion / cecp / recurrence-rate /
// sign-Markov-flux + the batch-5 family). It grows by design as the batch lands;
// a post-batch refactor may split it into per-feature submodules.

/// Rolling-window size for bar-close features: the number of PRIOR bar closes
/// each feature is computed over. Matches the project-wide inter-bar lookback
/// (`OPENDEVIATIONBAR_INTER_BAR_LOOKBACK_COUNT = 200`, the same window the spike
/// validated against).
pub const BAR_CLOSE_LOOKBACK_COUNT: usize = 200;

/// Petrosian fractal dimension over a bar-close lookback window.
///
/// Reversal-frequency measure: maps the number of sign changes in the
/// first-difference series onto a fractal-dimension scale.
///
/// ```text
/// N        = closes.len()
/// dx       = diff(closes)
/// N_delta  = num_zerocross(dx)   // sign-bit flips between consecutive dx values
/// PFD      = log10(N) / (log10(N) + log10(N / (N + 0.4 * N_delta)))
/// ```
///
/// `N_delta` is computed exactly as `antropy.num_zerocross(np.diff(x))`, i.e.
/// `np.diff(np.signbit(dx)).sum()`: it counts every index where the **sign bit**
/// of `dx` flips, across ALL consecutive pairs (no end-trimming). Because
/// `signbit(0.0) == false`, an exactly-flat step (`dx == 0`, i.e. two identical
/// consecutive closes — which DO occur on tick-quantized crypto data) does not
/// by itself constitute a crossing: a `…, -1, 0, +1, …` run counts as a single
/// crossing at the `-1 → 0` boundary, matching antropy.
///
/// This is the antropy / `pyrem` diff-sign variant (Goh et al. 2005), the same
/// definition `antropy.petrosian_fd` implements; the `0.4` constant is part of
/// the canonical Petrosian (1995) formula, not a tunable knob.
///
/// Empirical range on real bar-close data is `[1, ~1.07]` (the upper bound is
/// the perfectly-alternating window, `N_delta = N-2`). A flat window yields the
/// natural minimum `1.0` (antropy-faithful); `0.0` is returned only for
/// `len < 2` or non-finite input.
pub fn compute_bar_petrosian_fd(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 2 {
        return 0.0;
    }
    // Non-finite guard: a single NaN/Inf would poison the computation.
    if !closes.iter().all(|x| x.is_finite()) {
        return 0.0;
    }

    // N_delta = num_zerocross(diff(closes)): count sign-bit flips between
    // consecutive first-differences. dx[i] = closes[i+1] - closes[i].
    // (a - b with a == b yields +0.0 in IEEE-754, signbit false, so
    // is_sign_negative is a faithful np.signbit equivalent with no -0.0 hazard.)
    let mut prev_neg = (closes[1] - closes[0]).is_sign_negative();
    let mut n_delta: usize = 0;
    for i in 1..(n - 1) {
        let cur_neg = (closes[i + 1] - closes[i]).is_sign_negative();
        if cur_neg != prev_neg {
            n_delta += 1;
        }
        prev_neg = cur_neg;
    }

    let n_f = n as f64;
    let log_n = libm::log10(n_f);
    let pfd = log_n / (log_n + libm::log10(n_f / (n_f + 0.4 * (n_delta as f64))));

    if pfd.is_finite() { pfd } else { 0.0 }
}

/// Katz fractal dimension over a bar-close lookback window.
///
/// Path-length-to-displacement measure.
///
/// ```text
/// dists = |diff(closes)|
/// L     = sum(dists)                          // total path length
/// a     = mean(dists) = L / (closes.len()-1)  // mean step length
/// d     = max(|closes[i] - closes[0]|)        // max extent from the first point
/// KFD   = log10(L / a) / log10(d / a)
/// ```
///
/// This is the literal `antropy.katz_fd` 0.2.2 computation
/// (`np.log10(L / a) / np.log10(d / a)`). NOTE: antropy (following `mne-features`
/// / `pyrem`) uses the **sum-of-absolute-first-differences** path length, the
/// widely-used numerical form — NOT Katz's (1988) original Euclidean
/// `sqrt(1 + dy^2)` segment length. We intentionally match antropy (the FOSS
/// oracle the spike validated against), not the 1988 paper verbatim.
///
/// Katz FD is `>= 1` for any non-degenerate path on bar-close data (`d <= L`),
/// but is NOT universally `>= 1`: for spiky windows where the max extent from
/// the start is smaller than the mean step, antropy returns values `< 1` (even
/// negative), and so does this port — faithfully.
///
/// Robustness: antropy returns `NaN` on a flat window (`L = a = 0` → `0/0`) and
/// `+/-Inf` when `d == a` (`log10(1) == 0` denominator). Both are neutralised by
/// the leading finite guard and the trailing `is_finite` guard, so this never
/// emits NaN/Inf (it returns the `0.0` sentinel there). `0.0` is also returned
/// for `len < 2` or non-finite input.
pub fn compute_bar_katz_fd(closes: &[f64]) -> f64 {
    let n_closes = closes.len();
    if n_closes < 2 {
        return 0.0;
    }
    if !closes.iter().all(|x| x.is_finite()) {
        return 0.0;
    }

    // Total path length L and mean step length a.
    let n = n_closes - 1;
    let mut l = 0.0_f64;
    for i in 0..n {
        l += (closes[i + 1] - closes[i]).abs();
    }
    let a = l / n as f64;

    // Max extent from the first point.
    let x0 = closes[0];
    let mut d = 0.0_f64;
    for &x in closes {
        let dist = (x - x0).abs();
        if dist > d {
            d = dist;
        }
    }

    // antropy: log10(L / a) / log10(d / a)
    let kfd = libm::log10(l / a) / libm::log10(d / a);

    if kfd.is_finite() { kfd } else { 0.0 }
}

// ════════════════════════════════════════════════════════════════════════════
// Dispersion entropy (Rostaghi & Azami 2016) — bar-close lookback feature.
// ════════════════════════════════════════════════════════════════════════════

/// Number of NCDF amplitude classes. LOCKED feature-local constant (operator,
/// 2026-06-04): NOT a project-wide / tunable knob — see the ADR. Mirrors the
/// Python SSoT `entropy_candidates.DISP_C`.
const DISP_C: usize = 6;
/// Ordinal embedding dimension. LOCKED feature-local constant (`DISP_M`).
const DISP_M: usize = 2;
/// Ordinal embedding delay. LOCKED feature-local constant (`DISP_D`).
const DISP_D: usize = 1;

/// Size of the ordinal-pattern histogram for the locked (m = 2, d = 1) key encoding
/// `key = z[i] + c·z[i+1]` with classes `z ∈ [1, c]`: the maximum key is
/// `c + c·c = c² + c`, so the buffer needs `c² + c + 1` slots. Derived from
/// `DISP_C` (NOT a bare literal) so it tracks any future `c` change for this layout.
const DISPERSION_HISTOGRAM_SIZE: usize = DISP_C * DISP_C + DISP_C + 1;

// Compile-time invariant: the ordinal key `z[i] + c·z[i+d]` and the histogram
// size `c² + c + 1` are HARDCODED for the 2-element (m = 2, d = 1) embedding.
// `DISP_M` / `DISP_D` are otherwise consumed only by the (generalizable) window
// length `l = n - (m-1)·d` and the `ln(c^m)` normalizer, so flipping `DISP_M = 3`
// would compile and run WITHOUT tripping the debug_assert (a truncated 2-term key
// still fits the buffer) — yielding silently-wrong entropy. Pinning the layout
// here turns that into a build error: generalize the key encoding AND the
// histogram size before changing `DISP_M`/`DISP_D`. (No effect at m=2/d=1.)
const _: () = assert!(
    DISP_M == 2 && DISP_D == 1,
    "dispersion key `z[i] + c*z[i+d]` and DISPERSION_HISTOGRAM_SIZE = c^2+c+1 are \
     hardcoded for m=2/d=1; generalize the key encoding AND the histogram size \
     before changing DISP_M/DISP_D"
);

/// Standard-normal CDF Φ(z) via the identity `½(1 + erf(z/√2))`.
///
/// This matches `scipy.stats.norm.cdf` (= `scipy.special.ndtr`) — the Python
/// reference's NCDF map. `libm::erf` is the project-standard pure-Rust primitive
/// (the same reasoning that made #509 use `libm::log10` to match antropy). The
/// dispersion class assignment is DISCRETE, so a sub-ULP erf difference vs
/// scipy's Cephes `erf`/`erfc` cannot change the output unless a class boundary
/// is hit to within ~1e-13 — measure-zero on real bar-close data, and the F3
/// oracle (≤1e-9 vs the committed Python reference) proves it empirically.
#[inline]
fn ncdf_standard(z: f64) -> f64 {
    0.5 * (1.0 + libm::erf(z / std::f64::consts::SQRT_2))
}

/// Map a z-score to its dispersion amplitude class in `[1, c]` via the
/// Rostaghi-Azami (2016) NCDF binning `z = round(c·Φ + 0.5) ≡ floor(c·Φ) + 1`,
/// clamped to `[1, c]` — the convention used by the canonical FOSS impl
/// **EntropyHub** (`digitize(Φ, arange(0,1,1/c))`) and by the original paper.
///
/// `Φ(z) ∈ (0, 1)` ⇒ `Φ(z)·c ∈ [0, c)`, so `as usize` is the floor; `+ 1` lifts
/// it into `[1, c]`; the upper clamp catches the `Φ(z) → 1` tail saturation. The
/// mapping is TOTAL (always in `[1, c]`). Bounds + monotonicity are pinned by
/// `dispersion_class_bounded_and_monotone`.
///
/// NOTE (2026-06-06): this REPLACES the prior `int(c·Φ + 0.5)` (round-to-nearest)
/// binning. That was a hand-port slip — `.astype(int)` truncated where the paper's
/// `round(c·Φ + 0.5)` was intended — surfaced by the EntropyHub cross-validation
/// (R6) and confirmed against the paper. See the dispersion ADR.
#[inline]
fn dispersion_class(z: f64) -> usize {
    let bin = (ncdf_standard(z) * DISP_C as f64) as usize; // floor; Φ ∈ (0,1) ⇒ ≥ 0
    (bin + 1).min(DISP_C)
}

/// Dispersion entropy over a bar-close lookback window (Rostaghi & Azami 2016).
///
/// Bit-faithful port of the Python reference `entropy_candidates._dispersion_entropy`
/// with the LOCKED constants `c = 6`, `m = 2`, `d = 1`:
///
/// ```text
/// μ, σ   = mean(x), std(x)            // population std (ddof = 0, numpy default)
/// y      = Φ((x − μ) / σ)             // NCDF map into (0, 1)
/// z      = clip(floor(c·y) + 1, 1, c) // c amplitude classes, 1..=c (= paper round(c·y+0.5))
/// key_i  = z[i] + c·z[i+d]            // ordinal pattern (m = 2, d = 1)
/// p      = histogram(key) / (N − (m−1)·d)
/// DispEn = −Σ p·ln(p) / ln(c^m)       // normalized to [0, 1]
/// ```
///
/// Returns `NaN` for an undefined window — `len < 2`, a non-finite input, a flat
/// window (`σ == 0`, Python returns `nan`), or `< 2` embedding vectors. The
/// wiring layer maps a non-finite result to `None` (no spurious `0.0`, which is a
/// legitimate dispersion value: a single repeated pattern).
pub fn compute_bar_dispersion_entropy(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 2 {
        return f64::NAN;
    }
    if !closes.iter().all(|x| x.is_finite()) {
        return f64::NAN;
    }

    // Population mean and std (numpy `x.mean()` / `x.std()`, ddof = 0).
    let n_f = n as f64;
    let mean = closes.iter().sum::<f64>() / n_f;
    let var = closes
        .iter()
        .map(|&x| {
            let d = x - mean;
            d * d
        })
        .sum::<f64>()
        / n_f;
    let sigma = var.sqrt();
    // EXACT-zero check, bit-faithful to the Python reference `_dispersion_entropy`
    // (`if sigma == 0: return nan`) — required for the <=1e-9 oracle. A
    // `sigma < f64::EPSILON` threshold was considered (adversarial / SOTA review)
    // and deliberately REJECTED: it would DIVERGE from the SSoT for windows with
    // sigma in (0, EPSILON) — a divergence the oracle fixture (real closes, sigma
    // many orders of magnitude above EPSILON) cannot catch. Subnormal-but-nonzero
    // sigma is safe BY CONSTRUCTION: ncdf(z) = ½(1+erf(z/√2)) ∈ [0,1] (or NaN for
    // NaN z) ⇒ dispersion_class(z) ∈ [1,c] for ANY z incl. ±∞ (huge |z| → erf(±∞)=±1
    // → extreme classes 1/c) and NaN (→ class 1), so the histogram key stays in
    // range and the result is always NaN-or-[0,1], never ±Inf. (The
    // `dispersion_entropy_low_variance_safe` proptest exercises low-variance windows
    // but does not reach subnormal σ; the degenerate paths are pinned instead by the
    // `dispersion_class_*` and `dispersion_entropy_*_overflow*` unit tests below.)
    if sigma == 0.0 {
        return f64::NAN; // flat window → dispersion undefined (matches Python `nan`)
    }

    // NCDF map → c amplitude classes (see `dispersion_class`): clip(floor(c·Φ) + 1, 1, c).
    let classes: Vec<usize> = closes
        .iter()
        .map(|&x| dispersion_class((x - mean) / sigma))
        .collect();

    // Ordinal dispersion patterns (m = 2, d = 1): key = z[i] + c·z[i+d], z ∈ [1, c].
    // Max key = c + c·c, so a [0; c·c + c + 1] histogram fully covers the range.
    let l = n - (DISP_M - 1) * DISP_D; // number of embedding vectors
    if l <= 1 {
        return f64::NAN;
    }
    let mut counts = [0u32; DISPERSION_HISTOGRAM_SIZE];
    for i in 0..l {
        let key = classes[i] + DISP_C * classes[i + DISP_D];
        debug_assert!(
            key < DISPERSION_HISTOGRAM_SIZE,
            "dispersion ordinal key {key} exceeds histogram size {DISPERSION_HISTOGRAM_SIZE}"
        );
        counts[key] += 1;
    }

    // Shannon entropy of the pattern histogram, normalized by ln(c^m) = ln(36).
    let l_f = l as f64;
    let mut entropy = 0.0_f64;
    for &cnt in &counts {
        if cnt > 0 {
            let p = cnt as f64 / l_f;
            entropy -= p * libm::log(p);
        }
    }
    entropy / libm::log(DISP_C.pow(DISP_M as u32) as f64)
}

// ════════════════════════════════════════════════════════════════════════════
// bar_cecp_velocity — CECP (H, C) trajectory speed, bit-exact to ordpy==1.2.2
// ════════════════════════════════════════════════════════════════════════════

/// Half-window for the CECP 50/50 split. The (H, C) "complexity-entropy causality
/// plane" coordinate is computed separately on the first and second halves of the
/// lookback window; the feature is their Euclidean displacement. The split is
/// `n // 2` — a LOCKED feature-local convention (Axis-2 PASS-with-caveat; see the
/// ADR). Derived from `BAR_CLOSE_LOOKBACK_COUNT`, never a bare literal.
const CECP_HALF: usize = BAR_CLOSE_LOOKBACK_COUNT / 2;

// Invariant 2 (even window): the 50/50 split `n // 2` only halves exactly for an
// even window; a future odd `BAR_CLOSE_LOOKBACK_COUNT` would silently bias the
// halves (first half one element short). Make that a build error here.
const _: () = assert!(
    BAR_CLOSE_LOOKBACK_COUNT.is_multiple_of(2),
    "CECP 50/50 split (n // 2) requires an even BAR_CLOSE_LOOKBACK_COUNT"
);

// Invariant 3 (consume the #513 SSoT constants) + layout pin: the (H, C) ordinal
// embedding is HARDCODED for the 3-element, τ = 1 consecutive triplet
// `(x[i], x[i+1], x[i+2])` with an `[_; ORDINAL_M_FACTORIAL]` pattern histogram.
// A different ORDINAL_M / ORDINAL_TAU would silently mis-shape the embedding;
// pin the layout so it becomes a build error until the triplet logic is
// generalized (mirrors the dispersion `DISP_M == 2 && DISP_D == 1` guard).
const _: () = assert!(
    ORDINAL_M == 3 && ORDINAL_TAU == 1,
    "compute_bar_cecp_velocity hardcodes the (x[i], x[i+1], x[i+2]) τ=1 triplet \
     and the [_; ORDINAL_M_FACTORIAL] histogram; generalize the embedding before \
     changing ORDINAL_M/ORDINAL_TAU"
);

/// `(H, C)` = (normalised permutation entropy, MPR statistical complexity) for a
/// segment, BIT-FAITHFUL to `ordpy.complexity_entropy(seg, dx=ORDINAL_M,
/// taux=ORDINAL_TAU)` (ordpy 1.2.2 — the GROUND TRUTH; read from the installed
/// source in Stage 1):
///
/// ```text
/// n          = factorial(ORDINAL_M) = ORDINAL_M_FACTORIAL            (= 6)
/// p          = histogram(ordinal triplets) / l, OCCURRING patterns only
/// H          = −Σ p·ln p / ln n                                     (natural log)
/// u          = 1/n   (uniform distribution)
/// S(½(p+u))  = −Σ_occ ½(p+u)·ln ½(p+u)  −  ½u·ln(½u)·(n − #occurring)
/// S(p)/2     = −Σ p·ln p / 2 ;   S(u)/2 = ln n / 2
/// JS         = S(½(p+u)) − S(p)/2 − S(u)/2
/// js_div_max = −½·((n+1)/n·ln(n+1) + ln n − 2·ln 2n)
/// C          = H · JS / js_div_max                 (= H · JS · Q0)
/// ```
///
/// The not-occurring ordinal patterns (probability 0) each contribute the uniform
/// term `½u·ln(½u)` to `S(½(p+u))`, batched via `(n − #occurring)` exactly as
/// ordpy does. Caller guarantees `seg.len() >= ORDINAL_M` (production half = 100).
// Single-char bindings (h, c, p, u, n, l) deliberately mirror the ordpy
// complexity_entropy notation (H, C, p, uniform, n=#patterns, l=#triplets);
// verbose names would obscure the 1:1 correspondence with the reference formula.
#[allow(clippy::many_single_char_names)]
#[inline]
fn cecp_hc(seg: &[f64]) -> (f64, f64) {
    let l = seg.len() - (ORDINAL_M - 1) * ORDINAL_TAU; // # overlapping triplets
    debug_assert!(l >= 1, "cecp_hc needs at least one ordinal triplet");

    // Ordinal-pattern histogram. `ordinal_pattern_index_m3` (the project's proven
    // m=3 primitive) classes each triplet via the `<=` tie convention, which is
    // PROVEN classing-equivalent to numpy `np.argsort` (ordpy's `ordinal_sequence`)
    // over EVERY relative-order configuration incl. ties — enforced by the
    // exhaustive `cecp_ordinal_classing_matches_argsort_over_all_triplets` test
    // (all {0,1,2}^3, machine-checked against a committed numpy reference). The
    // integer labels differ from ordpy's argsort tuples, but H and C are invariant
    // under pattern relabeling, so the (H, C) pair is bit-identical.
    let mut counts = [0u32; ORDINAL_M_FACTORIAL];
    for i in 0..l {
        let cls = ordinal_pattern_index_m3(seg[i], seg[i + 1], seg[i + 2]);
        counts[cls] += 1;
    }

    let l_f = l as f64;
    let n = ORDINAL_M_FACTORIAL as f64;
    let u = 1.0 / n;

    // One pass over the OCCURRING patterns. Sum order is pattern-index 0..n; vs
    // ordpy's unique-sorted order this differs only at ~1e-16 (float non-assoc),
    // far under the 1e-9 gate.
    let mut s_p = 0.0_f64; //   −Σ p·ln p             (Shannon entropy, nat)
    let mut s_ppu = 0.0_f64; // −Σ ½(p+u)·ln ½(p+u)   (occurring patterns only)
    let mut occurring = 0usize;
    for &cnt in &counts {
        if cnt > 0 {
            occurring += 1;
            let p = cnt as f64 / l_f;
            s_p -= p * libm::log(p);
            let ppu = 0.5 * (p + u);
            s_ppu -= ppu * libm::log(ppu);
        }
    }

    let n_not = n - occurring as f64; // patterns with probability 0
    let half_u = 0.5 * u;
    let s_of_p_plus_u_over_2 = s_ppu - half_u * libm::log(half_u) * n_not;
    let s_of_p_over_2 = 0.5 * s_p;
    let s_of_u_over_2 = 0.5 * libm::log(n);
    let js_div = s_of_p_plus_u_over_2 - s_of_p_over_2 - s_of_u_over_2;

    // Invariant 1: `js_div_max` (hence Q0 = 1/js_div_max) is derived from the
    // constant `n` — NEVER a hardcoded `ln 6`. A literal `ln(6)` is exactly the
    // two-constants-disagree trap #513 surfaced (2.72e-13 drift).
    let js_div_max =
        -0.5 * (((n + 1.0) / n) * libm::log(n + 1.0) + libm::log(n) - 2.0 * libm::log(2.0 * n));

    let h = s_p / libm::log(n);
    let c = h * js_div / js_div_max;
    (h, c)
}

/// CECP velocity over a bar-close lookback window — the Euclidean displacement of
/// the (H, C) complexity-entropy coordinate between the first and second halves of
/// the window. Bit-faithful to the Python SSoT `cecp_reference.cecp_velocity`
/// (which wraps `ordpy==1.2.2`):
///
/// ```text
/// guard  : len < 2·(ORDINAL_M+2)  OR  any non-finite  OR  ptp == 0   → NaN
/// half   = len / 2                                   (50/50 split, n // 2)
/// (H1,C1) = complexity_entropy(closes[..half])
/// (H2,C2) = complexity_entropy(closes[half..])
/// vel    = hypot(H2 − H1, C2 − C1)                   ∈ [0, √2]
/// ```
///
/// Returns `NaN` for a degenerate window (mirrors the SSoT's `float("nan")`); the
/// wiring layer maps `NaN → None/NULL`. `0.0` is a VALID velocity (the two halves
/// share an (H, C) coordinate), so it is NEVER used as a sentinel. A flat HALF
/// (when the full window is not flat) is well-defined: its single ordinal pattern
/// gives `H = 0 ⇒ C = 0`, exactly as ordpy yields — no special-casing needed.
pub fn compute_bar_cecp_velocity(closes: &[f64]) -> f64 {
    let n = closes.len();
    // SSoT guard: a window shorter than 2·(m+2) cannot form two ordinal halves.
    if n < 2 * (ORDINAL_M + 2) {
        return f64::NAN;
    }
    if !closes.iter().all(|x| x.is_finite()) {
        return f64::NAN;
    }
    // `ptp == 0` (flat full window) → undefined; matches `np.ptp(x) == 0.0`.
    let mut lo = closes[0];
    let mut hi = closes[0];
    for &x in closes {
        if x < lo {
            lo = x;
        }
        if x > hi {
            hi = x;
        }
    }
    if hi - lo == 0.0 {
        return f64::NAN;
    }

    let half = n / 2;
    debug_assert!(
        n != BAR_CLOSE_LOOKBACK_COUNT || (half == CECP_HALF && n - half == CECP_HALF),
        "at the production lookback both CECP halves must equal CECP_HALF"
    );
    let (h1, c1) = cecp_hc(&closes[..half]);
    let (h2, c2) = cecp_hc(&closes[half..]);
    let vel = libm::hypot(h2 - h1, c2 - c1);
    if vel.is_finite() { vel } else { f64::NAN }
}

/// Locked recurrence-rate threshold for `bar_categorical_recurrence_rate`: exact
/// equality (`epsilon = 0`). The RQA recurrence rate counts EXACT-tie close levels;
/// there is no tolerance knob — `epsilon = 0` IS the defined categorical estimator.
const RECURR_EPSILON: f64 = 0.0;
// Compile-time pin: an edit that introduces a non-zero epsilon (turning the exact
// categorical estimator into a binned/tolerance one) must also touch this assert.
const _: () = assert!(
    RECURR_EPSILON == 0.0,
    "categorical recurrence rate is defined at epsilon = 0 (exact equality)"
);

/// Categorical recurrence rate over a bar-close lookback window (card 66).
///
/// `RR = Σ_l m_l·(m_l − 1) / (N·(N − 1))`, where `m_l` is the number of bars in the
/// window closing at that exact (FixedPoint-quantized) close level and `N` is the
/// window size. This is the RQA recurrence rate at `epsilon = 0` — the fraction of
/// ORDERED pairs of distinct window positions whose closes are exactly equal. In
/// plain terms: how often the price lands back on an exact level it already visited
/// recently (price-level stickiness).
///
/// Range `[0, 1]`: `0` when every close is distinct, `1` when all closes are equal.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window
/// (`len < 2`, or any non-finite close) — `0.0` is a VALID recurrence rate (an
/// all-distinct window), so it must never double as the undefined marker. The
/// wiring layer maps `NaN → None/NULL`. NEVER returns `±Inf` (the sole division is
/// by `N·(N−1) ≥ 2`, with a bounded-integer numerator).
///
/// Determinism: pure integer tie-counting plus one `f64` division — no `libm`, no
/// transcendental. Bit-identical across platforms and bit-identical to the Python
/// oracle (`tie_sum` and `N·(N−1)` are exact in `f64` for the production `N = 200`,
/// both well under `2^53`).
pub fn compute_bar_categorical_recurrence_rate(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 2 || !closes.iter().all(|x| x.is_finite()) {
        return f64::NAN;
    }
    // Group exact-equal close levels by sorting a copy under a total order and
    // summing run lengths. Sorting (not float-hashing) keeps the grouping exact and
    // deterministic. After a `total_cmp` sort, equal closes are adjacent; `to_bits`
    // equality is exact equality for the finite, positive, FixedPoint-quantized
    // closes this kernel sees, coinciding with the oracle's `np.unique` grouping
    // (the only finite values where bit-equality and `==` differ are ∓0.0, which
    // never occur in close-price data).
    let mut sorted: Vec<f64> = closes.to_vec();
    sorted.sort_unstable_by(f64::total_cmp);
    let mut tie_sum: u64 = 0;
    let mut run: u64 = 1;
    for w in sorted.windows(2) {
        if w[1].to_bits() == w[0].to_bits() {
            run += 1;
        } else {
            tie_sum += run * (run - 1);
            run = 1;
        }
    }
    tie_sum += run * (run - 1);
    let denom = (n as u64) * (n as u64 - 1);
    tie_sum as f64 / denom as f64
}

/// Locked sign alphabet size + transition lag for `bar_sign_markov_flux`: the
/// structural {down=-1, flat=0, up=+1} trichotomy at lag 1. `flat` is required (a
/// 2-state up/down chain is reversible by construction), not a chosen bin; the lag
/// is the defined transition step. Neither is a tuning knob.
const SIGNFLUX_STATES: usize = 3;
const SIGNFLUX_LAG: usize = 1;
const _: () = assert!(
    SIGNFLUX_STATES == 3 && SIGNFLUX_LAG == 1,
    "sign-Markov flux is defined on the 3-state sign alphabet at lag 1"
);

/// Sign-Markov in/out cyclic-current magnitude over a bar-close lookback window (card 98).
///
/// `logret = diff(ln(close))`; `signs = sign(logret) ∈ {−1, 0, +1}`; `C` = the 3×3
/// lag-1 transition-count matrix (index map −1→0, 0→1, +1→2); `f = C / total`;
/// `J = f − fᵀ`; feature = `Σ_{i<j} |J_ij| = (|C01−C10| + |C02−C20| + |C12−C21|) /
/// total` — the L1 magnitude of the net rotational "current" of the sign sequence.
/// In plain terms: whether up/flat/down moves cycle in one direction rather than
/// reverse (a directional-flow asymmetry).
///
/// Range `[0, 1]`: `0` on a detailed-balanced (reversible) chain, `1` on a perfect
/// 3-cycle (every off-diagonal transition one-directional).
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window —
/// any `close ≤ 0` or non-finite (`ln` undefined), fewer than 4 finite signs, or
/// zero transitions. The wiring maps `NaN → None/NULL`. NEVER returns `±Inf` (the
/// sole division is by `total ≥ 1` with a bounded-integer numerator).
///
/// Determinism: the value is an EXACT rational (integer `Σ|C_ij−C_ji|` over the
/// integer `total`), so the single `f64` division is the only rounding and is
/// bit-identical to the Python oracle. `libm::log` (computed once per close) derives
/// the substrate, but ONLY the sign of each logret is used — and
/// `sign(ln c_{t+1} − ln c_t) = sign(c_{t+1} − c_t)` for positive quantized closes
/// (a 1e-8 close step is ~100+ ULP from zero), so the sign sequence — and hence every
/// transition count — is library-rounding-immune and identical to the `np.log` oracle.
pub fn compute_bar_sign_markov_flux(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Undefined: ln requires strictly-positive finite closes.
    if n < 2 || !closes.iter().all(|&c| c.is_finite() && c > 0.0) {
        return f64::NAN;
    }
    // signs has length n − 1; mirror the oracle's "< 4 finite signs → NaN" guard.
    if n - 1 < 4 {
        return f64::NAN;
    }
    // ln once per close (matches the oracle's single `np.log` per element).
    let logs: Vec<f64> = closes.iter().map(|&c| libm::log(c)).collect();
    let sign_index = |t: usize| -> usize {
        let d = logs[t + 1] - logs[t];
        if d > 0.0 {
            2 // up   (+1)
        } else if d < 0.0 {
            0 // down (−1)
        } else {
            1 // flat ( 0)
        }
    };
    // 3×3 lag-1 transition counts over the sign sequence.
    let mut counts = [[0u64; SIGNFLUX_STATES]; SIGNFLUX_STATES];
    let mut total: u64 = 0;
    let mut prev = sign_index(0);
    for t in SIGNFLUX_LAG..n - 1 {
        let cur = sign_index(t);
        counts[prev][cur] += 1;
        total += 1;
        prev = cur;
    }
    if total == 0 {
        return f64::NAN;
    }
    // Exact rational: Σ upper-triangle |C_ij − C_ji| over the integer transition total.
    let num = counts[0][1].abs_diff(counts[1][0])
        + counts[0][2].abs_diff(counts[2][0])
        + counts[1][2].abs_diff(counts[2][1]);
    num as f64 / total as f64
}

/// Locked transition lag for `bar_ramsey_rothman_bicov_lag1`: lag 1, orders (2,1),
/// population std (ddof = 0), `sigma^3` normalizer. None is a tuning knob — `ddof=0`
/// is pinned by consistency with the `1/N` moment means, and `sigma^3` is the
/// card-113 Axis-2 finite-3rd-moment normalization.
const RRBICOV_LAG: usize = 1;
const _: () = assert!(
    RRBICOV_LAG == 1,
    "Ramsey-Rothman bicovariance is defined at lag 1"
);

/// Ramsey-Rothman bicovariance `gamma_{2,1}` at lag 1, `sigma^3`-normalised, over a
/// bar-close lookback window (card 113) — a time-irreversibility statistic.
///
/// `logret = diff(ln(close))`; `Y = logret − mean(logret)`; `sigma = sqrt(mean(Y^2))`
/// (population, ddof = 0); `gamma = mean(Y[t]^2·Y[t-1]) − mean(Y[t]·Y[t-1]^2)` over
/// the lag-1 pairs; feature = `gamma / sigma^3`. Plain terms: whether recent returns
/// "look different played backwards" (sharp-rise/slow-fall vs the reverse).
///
/// Scale-free (real-valued): scaling the closes by `k>0` leaves `logret` unchanged
/// (the constant `ln k` cancels in the difference); scaling `Y` by `a` leaves
/// `gamma/sigma^3` unchanged (3rd-order numerator over `sigma^3`). Time-reversal
/// ANTISYMMETRIC: reversing the logret series negates the value (exact — the reversed
/// `mean(Y[t]^2 Y[t-1])` equals the original `mean(Y[t] Y[t-1]^2)` over the same pairs).
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window —
/// any `close ≤ 0` or non-finite (`ln` undefined), fewer than 3 logrets, or
/// `sigma == 0` / non-finite. The wiring maps `NaN → None/NULL`. NEVER `±Inf` (the
/// `sigma == 0` guard removes the only division-by-zero; a trailing `is_finite` guard
/// absorbs any overflow).
///
/// UNLIKE the exact-rational bar_close features, this is a CONTINUOUS estimator
/// (`libm::log`, `libm::sqrt`, cubing), so the port matches the numpy oracle to
/// `<= 1e-9` (NOT bit-exact); the ~1e-13 FP divergence (libm vs glibc transcendentals
/// and summation order) sits ~4 orders under the gate. The locked conventions — plain
/// `diff(ln)`, de-mean BEFORE the moments, `ddof = 0`, `sigma^3` — are what the gate
/// pins. Means are computed as `mean(a) − mean(b)` (two `/N` divisions then subtract),
/// matching the numpy reference `np.mean(...) − np.mean(...)`.
pub fn compute_bar_ramsey_rothman_bicov_lag1(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Undefined: ln requires strictly-positive finite closes.
    if n < 2 || !closes.iter().all(|&c| c.is_finite() && c > 0.0) {
        return f64::NAN;
    }
    let m = n - 1; // logret length
    if m < 3 {
        return f64::NAN;
    }
    // logret = diff(ln(close)); ln once per close (matches the oracle's np.log).
    let logret: Vec<f64> = (0..m)
        .map(|t| libm::log(closes[t + 1]) - libm::log(closes[t]))
        .collect();
    // De-mean BEFORE the moments (locked convention).
    let mean = logret.iter().sum::<f64>() / m as f64;
    let y: Vec<f64> = logret.iter().map(|&r| r - mean).collect();
    // Population variance (ddof = 0) of the de-meaned series; sigma.
    let var = y.iter().map(|&v| v * v).sum::<f64>() / m as f64;
    let sigma = libm::sqrt(var);
    if !sigma.is_finite() || sigma == 0.0 {
        return f64::NAN;
    }
    // gamma_{2,1}(1): two third-moment means over the m-1 lag-1 pairs.
    let pairs = (m - RRBICOV_LAG) as f64;
    let mut s1 = 0.0_f64; // Σ Y[t]^2 · Y[t-1]
    let mut s2 = 0.0_f64; // Σ Y[t] · Y[t-1]^2
    for t in RRBICOV_LAG..m {
        let yt = y[t];
        let ytm1 = y[t - 1];
        s1 += yt * yt * ytm1;
        s2 += yt * ytm1 * ytm1;
    }
    let gamma = s1 / pairs - s2 / pairs;
    let result = gamma / (sigma * sigma * sigma);
    if result.is_finite() { result } else { f64::NAN }
}

/// Ehlers increment asymmetry over a bar-close lookback window (card 26) — a
/// dimensionless skewness of the bar-to-bar price changes (whether up-moves and
/// down-moves are lopsided in size).
///
/// `d_t = close_t − close_{t-1}` (RAW first differences, NOT logret);
/// `E = sum(d_t^3) / (sum(d_t^2))^1.5`.
///
/// Scale-free (real-valued): scaling the closes by `a` scales every `d_t` by `a`, and
/// `a^3·sum(d^3) / (a^2·sum(d^2))^1.5 = sum(d^3)/sum(d^2)^1.5` — units cancel. Odd:
/// negating the increments negates `E`; `E = 0` for symmetric increments.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window —
/// fewer than 3 closes, any non-finite close, or a flat window (`sum(d^2) == 0`). The
/// wiring maps `NaN → None/NULL`. NEVER `±Inf` (the `sum(d^2) == 0` guard removes the
/// only division-by-zero; a trailing `is_finite` guard absorbs any overflow). No
/// `close > 0` requirement — there is no `ln` here, only raw differences.
///
/// CONTINUOUS estimator (cubing, squaring, `libm::pow(·, 1.5)`), so the port matches
/// the numpy oracle to `<= 1e-9` (NOT bit-exact). The kernel and the oracle compute
/// `d_t` via the IDENTICAL IEEE subtraction on the same quantized closes (kernel
/// `closes[t] − closes[t-1]`, oracle `np.diff`), so `d_t` is bit-identical between
/// them regardless of Sterbenz (for consecutive same-magnitude closes it is also
/// exactly representable; a rare ≥2× jump merely rounds identically on both sides).
/// The only kernel↔oracle divergence is therefore the power sums and the `^1.5`
/// (`libm::pow(s2, 1.5)` vs numpy `s2 ** 1.5`), leaving a ~1e-17 residual ~12 orders
/// under the gate.
pub fn compute_bar_ehlers_increment_asymmetry(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 3 || !closes.iter().all(|&c| c.is_finite()) {
        return f64::NAN;
    }
    // d_t = close_t − close_{t-1} (raw first differences); power sums.
    let mut s2 = 0.0_f64; // Σ d^2
    let mut s3 = 0.0_f64; // Σ d^3
    for t in 1..n {
        let d = closes[t] - closes[t - 1];
        s2 += d * d;
        s3 += d * d * d;
    }
    if !s2.is_finite() || s2 == 0.0 {
        return f64::NAN; // flat window
    }
    let e = s3 / libm::pow(s2, 1.5);
    if e.is_finite() { e } else { f64::NAN }
}

/// Cox-Stuart trend z-score over a bar-close lookback window (card 23) — a robust
/// trend test: are later prices in the window systematically higher (or lower) than
/// earlier ones?
///
/// `m = n / 2`; pair the first half `x[0..m]` with the second half `x[n-m..n]`;
/// `S+ = #{ x[n-m+i] > x[i] }` (STRICT `>`, ties do not count);
/// `z = (S+ − m/2) / sqrt(m/4)`.
///
/// Locked convention: `c = floor(n/2)` pairing offset (`n − m`); the odd-n
/// drop-middle rule is MOOT at the production `n = 200` (even). Range: real z-score
/// (`[−10, 10]` at `n=200`, `m=100`). BIT-EXACT vs the numpy oracle: `S+` is an
/// integer count and `m/2 = 50`, `sqrt(m/4) = 5` are exact, so both compute the
/// identical `(S+ − 50)/5`.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for `n < 4` or any non-finite
/// close. The wiring maps `NaN → None/NULL`. NEVER `±Inf` (`sqrt(m/4) ≥ sqrt(0.5) > 0`
/// for `n ≥ 4`). A FLAT window is NOT undefined — no `second > first`, so `S+ = 0` and
/// `z = −(m/2)/sqrt(m/4)`, a valid "no up-trend" score. Rank/pairing invariant under
/// any positive monotone rescale (the comparison is sign-preserving and the count is
/// an integer).
pub fn compute_bar_cox_stuart_trend_z(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 4 || !closes.iter().all(|&c| c.is_finite()) {
        return f64::NAN;
    }
    let m = n / 2;
    // first half x[0..m] vs second half x[n-m..n]; pair i <-> (n-m)+i.
    let offset = n - m;
    let mut s_plus: u64 = 0;
    for i in 0..m {
        if closes[offset + i] > closes[i] {
            s_plus += 1;
        }
    }
    let mf = m as f64;
    let z = (s_plus as f64 - mf / 2.0) / libm::sqrt(mf / 4.0);
    if z.is_finite() { z } else { f64::NAN }
}

/// Groeneveld-Meeden b3 skewness over a bar-close lookback window (card 99) — a
/// robust, outlier-resistant skewness of recent returns, measured around the MEDIAN.
///
/// `x = diff(ln(close))` (logret); `med = median(x)` (numpy convention: odd-length →
/// middle order statistic, even-length → mean of the two central); `b3 = (mean(x) −
/// med) / mean(|x − med|)` — MAD about the MEDIAN (not sd, not IQR).
///
/// Range `[−1, 1]`. Scale-free (scaling closes by `k>0` leaves `logret` unchanged; a
/// logret rescale by `a` cancels in the `b3` ratio). Antisymmetric: negating the
/// logrets negates `b3` (`mean→−mean`, `med→−med`, MAD unchanged).
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window —
/// any `close ≤ 0` or non-finite (`ln`), fewer than 2 logrets, or a constant window
/// (`MAD == 0`). The wiring maps `NaN → None/NULL`. NEVER `±Inf` (the `MAD == 0` guard
/// removes the only division-by-zero; a trailing `is_finite` guard absorbs overflow).
///
/// CONTINUOUS estimator (`libm::log`, median selection, MAD), so the port matches the
/// numpy oracle to `<= 1e-9` (NOT bit-exact); the ~1e-13 FP divergence (libm vs glibc
/// `ln` + summation order) is ~4 orders under the gate. For the production `n = 200` the
/// logret length is 199 (ODD), so `med = sorted[99]` and the even-n averaging is MOOT.
/// Median stability under `ln` rounding: `ln` is STRICTLY MONOTONIC and both `libm` and
/// glibc are ~correctly-rounded (≤1 ULP), so two genuinely-distinct logrets (gap ≫ ULP)
/// keep their order in BOTH — they never reorder across Rust↔numpy. Only logrets equal to
/// within ~ULP can swap sort position, and those agree to ~ULP, so `sorted[99]` selects a
/// value identical to ~ULP either way (the ~1e-13 residual, not the inter-element gap).
pub fn compute_bar_groeneveld_meeden_b3_skewness(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Undefined: ln requires strictly-positive finite closes.
    if n < 2 || !closes.iter().all(|&c| c.is_finite() && c > 0.0) {
        return f64::NAN;
    }
    let m = n - 1; // logret length
    if m < 2 {
        return f64::NAN;
    }
    let x: Vec<f64> = (0..m)
        .map(|t| libm::log(closes[t + 1]) - libm::log(closes[t]))
        .collect();
    // median (numpy `np.median` convention) via a sorted copy.
    let mut sorted = x.clone();
    sorted.sort_unstable_by(f64::total_cmp);
    let median = if m % 2 == 1 {
        sorted[m / 2]
    } else {
        // `f64::midpoint` = exactly `(a + b) / 2.0` for |a|,|b| <= f64::MAX/2 (logrets
        // are ~1e0), so this stays bit-identical to numpy's even-n `np.median` mean.
        f64::midpoint(sorted[m / 2 - 1], sorted[m / 2])
    };
    let mean = x.iter().sum::<f64>() / m as f64;
    let mad = x.iter().map(|&v| (v - median).abs()).sum::<f64>() / m as f64;
    if !mad.is_finite() || mad == 0.0 {
        return f64::NAN; // constant window
    }
    let b3 = (mean - median) / mad;
    if b3.is_finite() { b3 } else { f64::NAN }
}

/// L-kurtosis ratio `tau4 = l4 / l2` over a bar-close lookback window (card 2) — a
/// robust, outlier-resistant 'fat-tailedness' of recent returns, built from L-moments
/// (linear combinations of order statistics) rather than central moments.
///
/// `x = sort(diff(ln(close)))` (logret, ascending). The probability-weighted moments
/// (Hosking 1990) use the PRODUCT-form order-statistic weights
/// `b_r = (1/n) Σ_i Π_{k=1..r}((i−k)/(n−k)) · x_(i)` (= `C(i−1,r)/C(n−1,r)`, i 1-indexed),
/// then `l2 = 2·b1 − b0`, `l4 = 20·b3 − 30·b2 + 12·b1 − b0`, and `tau4 = l4 / l2`.
///
/// Approx range `[−0.25, 1]` (Normal ≈ 0.1226); SAMPLE tau4 may dip a hair below the
/// population −0.25 floor (sampling). Scale-free (L-moment RATIO): scaling closes by
/// `k>0` leaves `logret` unchanged, and a logret rescale by `a` scales every `b_r` —
/// hence `l2`,`l4` — by `a`, cancelling in `l4/l2`.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window — any
/// `close ≤ 0` or non-finite (`ln`), fewer than 4 logrets (`n_closes < 5`; PWM needs
/// `b3`'s `Π((i−k)/(n−k))` up to `r=3`), or `l2 == 0` (a degenerate/constant window).
/// The wiring maps `NaN → None/NULL`. NEVER `±Inf` (the `l2 == 0` guard removes the only
/// division-by-zero; a trailing `is_finite` guard absorbs overflow).
///
/// CONTINUOUS estimator (`libm::log`, sorted PWM sums, division), so the port matches the
/// numpy oracle to `<= 1e-9` (NOT bit-exact); the ~1e-13 FP divergence (libm vs glibc
/// `ln` + summation/pairwise order) is ~4 orders under the gate. For the production
/// `n = 200` the logret length is 199 (≥ 4), so every PWM order is well-defined. Sort
/// stability under `ln` rounding: `ln` is STRICTLY MONOTONIC and ~correctly-rounded
/// (≤1 ULP) in both libm and glibc, so genuinely-distinct logrets keep their order in
/// BOTH — only equal-to-ULP logrets can swap, and they contribute identical-to-ULP terms,
/// so the PWM sums agree to the ~1e-13 residual, not the inter-element gap.
pub fn compute_bar_l_kurtosis_tau4(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Undefined: ln requires strictly-positive finite closes; PWM b3 needs >= 4 logrets
    // (n_closes >= 5).
    if n < 5 || !closes.iter().all(|&c| c.is_finite() && c > 0.0) {
        return f64::NAN;
    }
    let m = n - 1; // logret length (>= 4)
    let mut x: Vec<f64> = (0..m)
        .map(|t| libm::log(closes[t + 1]) - libm::log(closes[t]))
        .collect();
    x.sort_unstable_by(f64::total_cmp); // ascending order statistics x_(1) .. x_(n)
    let nf = m as f64;
    // PWM b_r via product-form order-statistic weights (matches the numpy oracle).
    let (mut b0, mut b1, mut b2, mut b3) = (0.0_f64, 0.0_f64, 0.0_f64, 0.0_f64);
    for (idx, &xi) in x.iter().enumerate() {
        let i = (idx + 1) as f64; // 1-indexed position
        let w1 = (i - 1.0) / (nf - 1.0);
        let w2 = (i - 1.0) * (i - 2.0) / ((nf - 1.0) * (nf - 2.0));
        let w3 = (i - 1.0) * (i - 2.0) * (i - 3.0) / ((nf - 1.0) * (nf - 2.0) * (nf - 3.0));
        b0 += xi;
        b1 += w1 * xi;
        b2 += w2 * xi;
        b3 += w3 * xi;
    }
    b0 /= nf;
    b1 /= nf;
    b2 /= nf;
    b3 /= nf;
    let l2 = 2.0 * b1 - b0;
    if !l2.is_finite() || l2 == 0.0 {
        return f64::NAN; // degenerate / constant window
    }
    let l4 = 20.0 * b3 - 30.0 * b2 + 12.0 * b1 - b0;
    let tau4 = l4 / l2;
    if tau4.is_finite() { tau4 } else { f64::NAN }
}

/// Average ranks of `x` (ties averaged) — a bit-faithful port of
/// `scipy.stats.rankdata(x, method="average")`. Sorts a copy by `total_cmp`, assigns
/// 1-indexed ranks, and gives every member of a tie group the mean of their positions
/// `(lo+1 + hi)/2`. Tie grouping uses exact float equality (`==`), matching scipy/numpy;
/// the only ties on real logret data are exact-`0.0` flat steps. Returns one rank per
/// input position (in INPUT/time order, not sorted order).
// float_cmp: exact `==` IS the scipy/numpy tie rule this port mirrors (see doc
// comment). An epsilon margin would merge near-ties into one rank group — a
// DISCRETE rank change that breaks the <=1e-9 bartels oracle. The only real ties
// on the ratio substrate are exact flat steps, seen identically by both rankers.
#[allow(clippy::float_cmp)]
fn bar_close_average_ranks(x: &[f64]) -> Vec<f64> {
    let n = x.len();
    let mut order: Vec<usize> = (0..n).collect();
    order.sort_by(|&a, &b| x[a].total_cmp(&x[b]));
    let mut ranks = vec![0.0_f64; n];
    let mut i = 0;
    while i < n {
        let mut j = i + 1;
        // Extend the tie group while values are exactly equal (scipy convention).
        while j < n && x[order[j]] == x[order[i]] {
            j += 1;
        }
        // Sorted 1-indexed positions of this group are (i+1)..=j; their mean is
        // (i+1 + j)/2 (sum of an arithmetic run / count).
        let avg = ((i + 1 + j) as f64) / 2.0;
        for &idx in &order[i..j] {
            ranks[idx] = avg;
        }
        i = j;
    }
    ranks
}

/// Bartels rank von-Neumann ratio `z` over a bar-close lookback window (card 19) — a
/// rank-based test of whether recent returns are randomly ordered or carry serial
/// structure (momentum vs mean-reversion), robust to the marginal return distribution
/// because it uses only the RANKS of the returns, not their magnitudes.
///
/// `x = close[i+1] / close[i]` (gross-return ratio, TIME-ordered); `R = avg-rank(x)`
/// (ties averaged); `RVN = Σ diff(R)² / Σ (R − mean R)²` (the rank von-Neumann ratio,
/// mean R = (n+1)/2); `z = (RVN − 2) / sqrt(σ²)` with the Bartels (1982) exact variance
/// `σ² = 4(n−2)(5n²−2n−9) / (5n(n+1)(n−1)²)`. RVN ≈ 2 under randomness; `z < 0` =
/// positive serial correlation (momentum), `z > 0` = mean-reversion.
///
/// SUBSTRATE NOTE: the statistic ranks the RETURNS. We rank the gross-return RATIO
/// `close[i+1]/close[i]` rather than the logret `ln(close[i+1]) − ln(close[i])`. Because
/// `ln` is STRICTLY INCREASING, `rank(ratio) ≡ rank(logret)` in exact arithmetic — the
/// identical statistic — but the ratio is a SINGLE correctly-rounded IEEE division, which
/// is bit-identical between Rust and numpy/glibc. The two-`ln`-plus-subtract logret form
/// is NOT (libm vs glibc `ln` differ by ≤1 ULP), and a 1-ULP disagreement on a near-tie
/// return would FLIP an adjacent rank pair → a DISCRETE `~1e-3` jump in `RVN` (rank stats
/// are order-discontinuous, unlike the value-based continuous features). Ranking the
/// IEEE-stable ratio removes that hazard entirely → the gate holds at residual ~0.
///
/// Range: a real z-score (unbounded; a perfectly monotone series gives `z ≈ −sqrt(n)`).
/// RANK-based, so INVARIANT under any strictly-monotone transform of the closes — in
/// particular scale-free: scaling closes by `k>0` leaves the ratio `c[i+1]/c[i]` unchanged
/// (the factor cancels in the division — bit-exactly for a power-of-two `k`), so the ranks
/// and `z` are unchanged.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window — any
/// `close ≤ 0` or non-finite (`ln`), fewer than 4 logrets (`n_closes < 5`), or a zero
/// rank variance (`Σ(R − mean R)² == 0`, i.e. a constant/all-tied window). The wiring
/// maps `NaN → None/NULL`. NEVER `±Inf` (the rank-variance guard removes the only
/// division-by-zero; a trailing `is_finite` guard absorbs overflow).
///
/// EFFECTIVELY BIT-EXACT (not "continuous"): the ratio substrate is bit-identical
/// Rust↔numpy (correctly-rounded IEEE division), so the ranks match exactly; ranks are
/// half-integers, every Σ of squared rank-quantities is an exact multiple of `0.25` well
/// under `2^53` (no rounding), and the remaining ops (RVN division, the σ² integer-ratio,
/// one correctly-rounded IEEE `libm::sqrt`, the final division) are bit-identical to the
/// numpy oracle. The lone tie source — `ratio == 1.0` from flat steps (`close[i+1] ==
/// close[i]`) — is seen identically by both rankers; distinct returns are separated by
/// `≫1` ULP so the rank ORDER is stable. Residual ~0.
// Single-char bindings (n, m, x, r, d, z) deliberately mirror the Bartels (1982)
// statistic notation used in the doc comment and the numpy oracle (R, RVN, z);
// verbose names would obscure the 1:1 correspondence with the reference formula.
#[allow(clippy::many_single_char_names)]
pub fn compute_bar_bartels_rank_vn_ratio(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Undefined: ratios need strictly-positive finite closes; the rank stat needs >= 4
    // returns (n_closes >= 5).
    if n < 5 || !closes.iter().all(|&c| c.is_finite() && c > 0.0) {
        return f64::NAN;
    }
    let m = n - 1; // return count (>= 4)
    // Gross-return ratio substrate: rank-identical to the logret (ln is monotone) but a
    // single IEEE division -> bit-identical Rust↔numpy, so ranks never diverge on a near-tie.
    let x: Vec<f64> = (0..m).map(|t| closes[t + 1] / closes[t]).collect();
    let r = bar_close_average_ranks(&x);
    let mf = m as f64;
    // Exact mean rank (1+m)/2 — the midpoint of the 1-indexed rank range; holds under
    // tie-averaging. `f64::midpoint` = exactly `(a + b) / 2.0` here (values << MAX/2).
    let rbar = f64::midpoint(mf, 1.0);
    let mut denom = 0.0_f64;
    for &ri in &r {
        let d = ri - rbar;
        denom += d * d;
    }
    if !denom.is_finite() || denom == 0.0 {
        return f64::NAN; // zero rank variance (constant / all-tied window)
    }
    let mut num = 0.0_f64;
    for w in r.windows(2) {
        let d = w[1] - w[0];
        num += d * d;
    }
    let rvn = num / denom;
    let sigma2 = 4.0 * (mf - 2.0) * (5.0 * mf * mf - 2.0 * mf - 9.0)
        / (5.0 * mf * (mf + 1.0) * (mf - 1.0) * (mf - 1.0));
    if !sigma2.is_finite() || sigma2 <= 0.0 {
        return f64::NAN;
    }
    let z = (rvn - 2.0) / libm::sqrt(sigma2);
    if z.is_finite() { z } else { f64::NAN }
}

/// Ordinal (distinct, 0-indexed) ranks with STABLE tie-break by input position —
/// the exact permutation numpy's double-argsort (`argsort(argsort(x, kind="mergesort"))`)
/// produces. `Vec::sort_by` is a stable sort, so equal values keep input order (the
/// explicit index tie-break in the comparator makes the contract self-documenting).
/// Unlike [`bar_close_average_ranks`] (scipy tie-AVERAGING), every position gets a
/// DISTINCT rank `0..n-1` — the empirical-copula grid of the Hoeffding Φ² kernel
/// indexes a permutation matrix, which requires distinct coordinates.
///
/// SIGNED-ZERO CAVEAT: `total_cmp` orders `-0.0 < +0.0` while numpy's sort treats
/// them as EQUAL (stable, input order kept) — a `[-0.0, +0.0]` tie would rank
/// differently. UNREACHABLE for the hoeffding substrates (x = ratio of guarded
/// strictly-positive closes, so x > 0; y = `(close_time − open_time) as f64`,
/// an i64-difference cast that can only produce `+0.0`) — but a future caller
/// feeding computed, possibly `-0.0`, values would break the numpy-parity claim.
fn bar_close_ordinal_ranks(x: &[f64]) -> Vec<usize> {
    let n = x.len();
    let mut order: Vec<usize> = (0..n).collect();
    order.sort_by(|&a, &b| x[a].total_cmp(&x[b]).then(a.cmp(&b)));
    let mut ranks = vec![0usize; n];
    for (pos, &idx) in order.iter().enumerate() {
        ranks[idx] = pos;
    }
    ranks
}

/// Hoeffding's Φ² (phi-squared) between the bar-to-bar gross-return ratio and the
/// bar DURATION over a bar-close lookback window (card 78,
/// `bar_hoeffding_phi_squared_midreturn_duration`) — a rank-based dependence measure
/// that detects a relationship of ANY shape (not just linear/monotone) between how
/// far price moved and how long the bar took to form. Effectively `[0, 1]`: `0` iff
/// the empirical copula equals the independence copula on the rank grid; `→1` for
/// perfect (co- or counter-) monotone dependence. (Finite-m grid bias: ONLY under
/// perfect monotone dependence the grid MEAN exceeds 1 by `O(1/m)` — exactly
/// 1.0000631… at m=199; the 90· constant normalizes the INTEGRAL version, and the
/// grid mean converges to it from above. Unreachable on real market data.)
///
/// The FIRST two-substrate bar_close kernel: `closes` and `durations_us` are the
/// aligned trailing windows (`durations_us[i] = close_time − open_time` of bar `i`,
/// integer microseconds exact in f64 — only its RANKS are consumed).
///
///   m       = n − 1 pairs (returns within the window, bartels precedent)
///   x[t]    = closes[t+1] / closes[t]      (gross-return ratio, TIME-ordered)
///   y[t]    = durations_us[t+1]            (duration of the bar the return lands in
///                                           — probe convention: `ret_dur` row i pairs
///                                           `logret[i]` with `dur[i]`)
///   rk, sk  = ordinal STABLE ranks of x, y (ties keep input order — the numpy
///                                           mergesort double-argsort permutation)
///   C(i,j)  = #{k: rk[k] < i ∧ sk[k] < j} / m     (empirical copula, grid 1..m)
///   Φ²      = 90 · mean_{i,j} ( C(i,j) − (i/m)(j/m) )²
///
/// SUBSTRATE NOTE (I2, same hazard class as bartels): the statistic ranks the
/// RETURNS. We rank the gross-return RATIO `closes[t+1]/closes[t]`, not the logret —
/// `ln` is strictly increasing so `rank(ratio) ≡ rank(logret)` exactly, but the ratio
/// is a SINGLE correctly-rounded IEEE division, bit-identical Rust↔numpy, so a 1-ULP
/// `ln` disagreement can never flip an adjacent rank pair (rank statistics are
/// order-DISCONTINUOUS; one flipped pair is a discrete jump, not a 1e-16 wiggle).
///
/// EFFECTIVELY BIT-EXACT vs the committed Python oracle: the ranks match exactly
/// (identical substrate bits + identical stable tie order), the copula grid is
/// integer counts divided by `m` (same single IEEE division per cell), the
/// independence grid is `(i/m)·(j/m)` (same two roundings, same order), and the only
/// divergence source is the final 90·mean reduction (numpy pairwise vs sequential
/// summation) — bounded ~1e-13, four orders under the 1e-9 gate.
///
/// RANK-based, so scale-free on BOTH axes: closes×k (k>0) leaves every ratio's rank
/// unchanged (bit-exactly for power-of-two k); any strictly-increasing transform of
/// the durations leaves `sk` unchanged → identical Φ² bits.
///
/// Robustness contract: returns `NaN` (NOT a sentinel — `0.0` is a VALID Φ², exact
/// grid-independence) for an undefined window: `closes.len() != durations_us.len()`,
/// fewer than 4 pairs (`n < 5`), any `close ≤ 0` or non-finite, or any non-finite
/// duration. NEVER `±Inf` (all divisors are positive integers; a trailing guard
/// absorbs the impossible). The wiring maps `NaN → None/NULL`.
pub fn compute_bar_hoeffding_phi_squared_midreturn_duration(
    closes: &[f64],
    durations_us: &[f64],
) -> f64 {
    let n = closes.len();
    // Undefined: ratios need strictly-positive finite closes; the copula grid needs
    // >= 4 pairs (n_closes >= 5); the two substrates must be the SAME window.
    if n < 5 || durations_us.len() != n {
        return f64::NAN;
    }
    if !closes.iter().all(|&c| c.is_finite() && c > 0.0)
        || !durations_us.iter().all(|&d| d.is_finite())
    {
        return f64::NAN;
    }
    let m = n - 1; // pair count (>= 4)
    let x: Vec<f64> = (0..m).map(|t| closes[t + 1] / closes[t]).collect();
    let rk = bar_close_ordinal_ranks(&x);
    let sk = bar_close_ordinal_ranks(&durations_us[1..]);
    // Permutation-indicator grid + 2-D prefix sums — the numpy
    // `cumsum(cumsum(ind, axis=0), axis=1)` construction with exact integer counts.
    let mut grid = vec![0u32; m * m];
    for k in 0..m {
        grid[rk[k] * m + sk[k]] = 1;
    }
    for i in 1..m {
        for j in 0..m {
            grid[i * m + j] += grid[(i - 1) * m + j];
        }
    }
    for i in 0..m {
        for j in 1..m {
            grid[i * m + j] += grid[i * m + j - 1];
        }
    }
    let mf = m as f64;
    let mut sum = 0.0_f64;
    for i in 0..m {
        let gi = (i + 1) as f64 / mf;
        for j in 0..m {
            let cemp = f64::from(grid[i * m + j]) / mf;
            let d = cemp - gi * ((j + 1) as f64 / mf);
            sum += d * d;
        }
    }
    let phi2 = 90.0 * (sum / (mf * mf));
    if phi2.is_finite() { phi2 } else { f64::NAN }
}

/// Mean forward visibility horizon over a bar-close lookback window (card 59,
/// `bar_hvg_forward_visibility_horizon_mean`) — for each bar `t`, the first-passage
/// distance to a close at-or-above its own level,
///
///   `L+(t) = min { k >= 1 : closes[t+k] >= closes[t] }`   (within the window)
///
/// trailing non-crossings are RIGHT-CENSORED (no horizon recorded for a bar whose
/// level is never re-reached before the window ends), and the feature is the MEAN
/// of the recorded horizons. Interpretation: how many bars a price level typically
/// stays "underwater" before being reclaimed — ~1 in relentless up-drift (every
/// next close is at-or-above), large in long drawdown stretches. Range `[1, n-1]`.
///
/// HVG equivalence (why the card is named `hvg_*`): under horizontal-visibility
/// rules, node `t` sees forward exactly until the first bar at-or-above its own
/// height blocks it — and the blocker is the last forward-visible node; `L+(t)` is
/// that blocker's distance. The first-passage form is tie-free and needs NO graph
/// construction (card 60 is the one that builds the actual HVG adjacency).
///
/// IMPLEMENTATION: a single monotonic-stack "next greater-or-equal element" pass —
/// O(n) time, one `Vec<usize>` scratch — instead of the probe's naive O(n²) scan.
/// The two produce IDENTICAL integer horizons (the unit tests carry the naive
/// construction as the in-Rust independent leg); the oracle generator additionally
/// cross-checks both constructions per window, fail-closed.
///
/// BIT-EXACT vs the committed Python oracle (residual 0): closes are FixedPoint-
/// quantized identically on both sides, `>=` comparisons involve no arithmetic,
/// horizons are small integers whose sum (`< n² = 40,000`) is exact in f64, and
/// the single final division is correctly-rounded IEEE — no tolerance consumed.
///
/// Scale note: NOT scale-free in the bit-exact sense under arbitrary rescale, but
/// `>=` comparisons are preserved by any strictly-increasing transform of the
/// closes — in particular positive scaling — so the integer horizons (and hence
/// the value bits) are invariant under `closes × k, k > 0`, PROVIDED the rescale
/// neither overflows to `±Inf` nor collapses distinct closes (denormal
/// underflow); both are unreachable for real price magnitudes and the pow-2
/// factors the harness pins.
///
/// Robustness contract: returns `NaN` (NOT a sentinel) for an undefined window —
/// `n < 2`, any non-finite close, or ALL bars right-censored (a strictly-decreasing
/// window records no horizon; `1.0` is a VALID mean, so it must never double as
/// the undefined marker). NEVER `±Inf` (the mean of integers in `[1, n-1]` with a
/// non-empty count is finite by construction). The wiring maps `NaN → None/NULL`.
pub fn compute_bar_hvg_forward_visibility_horizon_mean(closes: &[f64]) -> f64 {
    let n = closes.len();
    // Exactness bound for external callers: the horizon sum is < n²/2, so the
    // u64 accumulator is exact for n < ~6e9 (production n = 200; a wrap would
    // need a ~49 GB slice). Guarded in debug builds.
    debug_assert!(n < (1 << 32), "horizon sum may wrap u64 for n >= 2^32");
    if n < 2 || !closes.iter().all(|&c| c.is_finite()) {
        return f64::NAN;
    }
    // Monotonic stack of indices with strictly-descending close values: when
    // closes[j] >= closes[stack top], the top's first-passage blocker is j.
    let mut sum: u64 = 0;
    let mut count: u64 = 0;
    let mut stack: Vec<usize> = Vec::with_capacity(n);
    for (j, &cj) in closes.iter().enumerate() {
        while let Some(&t) = stack.last() {
            if cj >= closes[t] {
                stack.pop();
                sum += (j - t) as u64;
                count += 1;
            } else {
                break;
            }
        }
        stack.push(j);
    }
    // Indices still on the stack are right-censored (their level was never
    // re-reached) — dropped, exactly like the reference construction.
    if count == 0 {
        return f64::NAN;
    }
    let v = (sum as f64) / (count as f64);
    if v.is_finite() { v } else { f64::NAN }
}

/// HVG adjacency via the tie-safe running-max LEFT-SCAN (the shared visibility-
/// graph helper, batch-6): bars `i < j` are connected iff every close strictly
/// between them is LOWER than both — `max(x[i+1..j-1]) < min(x[i], x[j])`;
/// adjacent bars always connect; an EQUAL intermediate BLOCKS (strict `<`).
///
/// For each `j` walk `i = j-1..0` keeping the running max `m` of the
/// intermediates: edge iff `m < x[i] && m < x[j]`; fold `x[i]` into `m`; stop
/// once `m >= x[j]` (everything farther is blocked). O(n · scan) with short
/// scans on real data.
///
/// TIE HAZARD (oracle-leg incident, card 60 S1): a textbook monotonic-stack
/// construction admits SPURIOUS edges on equal-height runs (the `[2,1,1,2]`
/// class — a bar "sees" through an equal-height wall) — the generator's ts2vg
/// fail-closed leg vetoed it with 3 phantom edges on the very first real
/// window. The left-scan applies the rule literally and is verified against a
/// naive O(n³) construction (in-Rust tie battery) and ts2vg (in-generator, all
/// 19,602 windows).
///
/// Returns (sorted-neighbour adjacency lists, edge bitset `i*n+j` for `i<j`).
fn hvg_adjacency_leftscan(x: &[f64]) -> (Vec<Vec<usize>>, Vec<u64>) {
    let n = x.len();
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    let mut bits = vec![0u64; (n * n).div_ceil(64)];
    let connect = |adj: &mut Vec<Vec<usize>>, bits: &mut Vec<u64>, i: usize, j: usize| {
        adj[i].push(j);
        adj[j].push(i);
        let k = i * n + j;
        bits[k / 64] |= 1u64 << (k % 64);
    };
    for j in 1..n {
        connect(&mut adj, &mut bits, j - 1, j);
        let mut m = x[j - 1];
        let mut i = j as isize - 2;
        while i >= 0 && m < x[j] {
            let iu = i as usize;
            if m < x[iu] {
                connect(&mut adj, &mut bits, iu, j);
            }
            if x[iu] > m {
                m = x[iu];
            }
            i -= 1;
        }
    }
    (adj, bits)
}

// Single-char bindings mirror graph-index notation: nodes `a`/`b` ordered into
// `i < j`, flattened to bitset key `k` over an `n × n` adjacency matrix.
#[allow(clippy::many_single_char_names)]
#[inline]
fn hvg_has_edge(bits: &[u64], n: usize, a: usize, b: usize) -> bool {
    let (i, j) = if a < b { (a, b) } else { (b, a) };
    let k = i * n + j;
    bits[k / 64] & (1u64 << (k % 64)) != 0
}

/// HVG time-directed clustering asymmetry — signed MEAN-GAP of past-vs-future
/// local clustering — over a bar-close lookback window (card 60,
/// `bar_vg_time_directed_clustering_meangap`; the card slug's `_kld` is RENAMED:
/// a true D_KL needs an unpinned density estimator, and what was evaluated and
/// ships is this mean gap). A detector for the ARROW OF TIME in the close
/// series: build the horizontal visibility graph, split each node's neighbours
/// by time index into past (`u < v`) and future (`u > v`), take the local
/// clustering coefficient of each set (fraction of connected pairs within it),
/// and average:
///
///   C_v^past = links(past)/(kp·(kp−1)/2)  (defined when kp ≥ 2; future alike)
///   value    = mean_v(C_v^past) − mean_v(C_v^future)   (means over defined nodes)
///
/// `0` = clustering-symmetric (time-reversible); real markets sit slightly
/// positive (~+0.05..0.06 on the committed fixtures). Range `[-1, 1]`.
/// Time-reversal flips the sign EXACTLY in real arithmetic — but only to
/// ~1e-12 in IEEE (reversal permutes the mean's summation order; oracle-leg
/// incident P2) — property tests must NOT assert `to_bits` negation.
///
/// NEAR-BIT-EXACT vs the committed Python oracle: identical FixedPoint inputs,
/// identical integer edge sets and link counts, identical per-node divisions;
/// the only divergence is the final mean reduction order (numpy pairwise vs
/// sequential) — residual ~1e-15, four+ orders under the 1e-9 gate.
///
/// Invariant under any strictly-increasing close transform (visibility
/// comparisons preserved) — bit-exactly for a power-of-two rescale.
///
/// Robustness contract: returns `NaN` (NOT a sentinel — `0.0` is a VALID,
/// meaningful value: perfect time-symmetry) for an undefined window — `n < 5`,
/// any non-finite close, or either defined-node list empty (e.g. a strictly
/// monotone window whose HVG is a bare path). NEVER `±Inf`. The wiring maps
/// `NaN → None/NULL`.
pub fn compute_bar_vg_time_directed_clustering_meangap(closes: &[f64]) -> f64 {
    let n = closes.len();
    if n < 5 || !closes.iter().all(|&c| c.is_finite()) {
        return f64::NAN;
    }
    let (adj, bits) = hvg_adjacency_leftscan(closes);
    let mut past_sum = 0.0_f64;
    let mut past_cnt = 0u32;
    let mut fut_sum = 0.0_f64;
    let mut fut_cnt = 0u32;
    let mut past: Vec<usize> = Vec::with_capacity(8);
    let mut fut: Vec<usize> = Vec::with_capacity(8);
    for (v, nbrs) in adj.iter().enumerate() {
        past.clear();
        fut.clear();
        for &u in nbrs {
            if u < v {
                past.push(u);
            } else {
                fut.push(u);
            }
        }
        for (set, sum, cnt) in [
            (&past, &mut past_sum, &mut past_cnt),
            (&fut, &mut fut_sum, &mut fut_cnt),
        ] {
            let k = set.len();
            if k >= 2 {
                let mut links = 0u32;
                for a in 0..k {
                    for b in a + 1..k {
                        if hvg_has_edge(&bits, n, set[a], set[b]) {
                            links += 1;
                        }
                    }
                }
                // Mirror the oracle exactly: 2.0*links / (k*(k-1)) — one f64 division.
                *sum += 2.0 * f64::from(links) / ((k * (k - 1)) as f64);
                *cnt += 1;
            }
        }
    }
    if past_cnt == 0 || fut_cnt == 0 {
        return f64::NAN;
    }
    let val = past_sum / f64::from(past_cnt) - fut_sum / f64::from(fut_cnt);
    if val.is_finite() { val } else { f64::NAN }
}

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

    /// `[1, 2, 3]` z-maps (floor(c·Φ)+1) to classes `[1, 4, 6]` (all safely mid-bin),
    /// giving two distinct ordinal pairs with equal probability → entropy =
    /// ln(2)/ln(36) EXACTLY, independent of erf sub-ULP precision (classes robust).
    #[test]
    fn dispersion_entropy_three_distinct_classes_is_ln2_over_ln36() {
        let v = compute_bar_dispersion_entropy(&[1.0, 2.0, 3.0]);
        let expected = 2f64.ln() / 36f64.ln();
        assert!((v - expected).abs() < 1e-12, "got {v}, want {expected}");
    }

    /// A flat window has σ == 0 → dispersion undefined → NaN (Python parity).
    #[test]
    fn dispersion_entropy_flat_window_is_nan() {
        assert!(compute_bar_dispersion_entropy(&[42.0; 50]).is_nan());
    }

    /// Non-finite input and too-short input are undefined → NaN.
    #[test]
    fn dispersion_entropy_degenerate_inputs_are_nan() {
        assert!(compute_bar_dispersion_entropy(&[1.0]).is_nan());
        assert!(compute_bar_dispersion_entropy(&[1.0, f64::NAN, 3.0]).is_nan());
    }

    /// On a varied window the result is finite and normalized to [0, 1].
    #[test]
    fn dispersion_entropy_bounded_on_varied_window() {
        let w: Vec<f64> = (0..BAR_CLOSE_LOOKBACK_COUNT)
            .map(|i| 100.0 + (i as f64 * 0.7).sin() * 3.0 + i as f64 * 0.01)
            .collect();
        let v = compute_bar_dispersion_entropy(&w);
        assert!(
            v.is_finite() && (0.0..=1.0).contains(&v),
            "out of [0,1]: {v}"
        );
    }

    /// Φ(z) = ½(1 + erf(z/√2)) must be centered (Φ(0) == 0.5 exactly), symmetric
    /// (Φ(z) + Φ(−z) == 1), non-decreasing, bounded in [0, 1], and saturate in the
    /// tails — the NCDF correctness the discrete class assignment relies on.
    #[test]
    fn ncdf_standard_is_centered_and_symmetric() {
        assert_eq!(ncdf_standard(0.0), 0.5, "Φ(0) must be exactly 0.5");
        for &z in &[0.1, 0.5, 1.0, 2.0, 3.0, 7.5] {
            let s = ncdf_standard(z) + ncdf_standard(-z);
            assert!(
                (s - 1.0).abs() < 1e-12,
                "Φ(z)+Φ(-z) must be 1, got {s} at z={z}"
            );
        }
        // Non-decreasing + bounded on a fine grid (tails saturate to 1.0, hence ≥).
        let mut prev = ncdf_standard(-8.0);
        for k in -79..=80 {
            let y = ncdf_standard(k as f64 * 0.1);
            assert!(
                y >= prev - 1e-15 && (0.0..=1.0).contains(&y),
                "Φ non-monotone/out-of-range at k={k}: {y}"
            );
            prev = y;
        }
        assert!(
            ncdf_standard(40.0) >= 1.0 - 1e-12,
            "right tail saturates to 1"
        );
        assert!(ncdf_standard(-40.0) <= 1e-12, "left tail saturates to 0");
    }

    /// `dispersion_class` is TOTAL: always in [1, c] for ANY z (incl. ±extremes,
    /// pathological f64::MIN/MAX, and exact NCDF boundaries), and non-decreasing in
    /// z. Proves the clamp + rounding can never produce an out-of-range class.
    #[test]
    fn dispersion_class_bounded_and_monotone() {
        // Incl. NaN and ±∞: a totally undefined / saturated z must still yield a
        // valid class (NaN→1 via NaN-as-usize=0; +∞→c; -∞→1), never an OOB key.
        for &z in &[
            f64::NAN,
            f64::NEG_INFINITY,
            f64::MIN,
            -1e12,
            -40.0,
            -1.0,
            0.0,
            1.0,
            40.0,
            1e12,
            f64::MAX,
            f64::INFINITY,
        ] {
            let c = dispersion_class(z);
            assert!(
                (1..=DISP_C).contains(&c),
                "class {c} out of [1,{DISP_C}] at z={z}"
            );
        }
        assert_eq!(
            dispersion_class(f64::NAN),
            1,
            "NaN z → class 1 (no OOB key)"
        );
        assert_eq!(dispersion_class(f64::INFINITY), DISP_C, "+∞ z → class c");
        assert_eq!(dispersion_class(f64::NEG_INFINITY), 1, "-∞ z → class 1");
        // Sweep across the NCDF domain (crosses every interior class boundary):
        // classes stay in range and never decrease as z increases.
        let mut prev = 0usize;
        let mut z = -6.0_f64;
        while z <= 6.0 {
            let c = dispersion_class(z);
            assert!((1..=DISP_C).contains(&c), "class {c} out of range at z={z}");
            assert!(
                c >= prev,
                "class must be non-decreasing in z (got {c} after {prev})"
            );
            prev = c;
            z += 0.05;
        }
        // The full class range is reachable at the tails.
        assert_eq!(dispersion_class(-40.0), 1, "deep left tail → class 1");
        assert_eq!(dispersion_class(40.0), DISP_C, "deep right tail → class c");
    }

    /// Mean-overflow / NaN-z path: a window of `f64::MAX` makes the mean overflow
    /// to +∞, so every `z = (x - mean)/σ` is NaN. dispersion_class(NaN) = 1 for all
    /// elements → a single ordinal pattern → entropy 0.0 (finite). The contract is
    /// only "never ±Inf"; assert finiteness, not a specific value.
    #[test]
    fn dispersion_entropy_mean_overflow_is_finite_never_inf() {
        let w = [f64::MAX; BAR_CLOSE_LOOKBACK_COUNT];
        let v = compute_bar_dispersion_entropy(&w);
        assert!(!v.is_infinite(), "must never be ±Inf, got {v}");
        // NaN-z everywhere collapses to one class → degenerate (zero) entropy.
        assert!(
            v == 0.0 || v.is_nan(),
            "expected degenerate 0.0 or NaN, got {v}"
        );
    }

    /// Variance-overflow path: alternating ±1e200 makes σ overflow to +∞ while the
    /// mean stays finite, so every `z = (x-mean)/∞ = ±0` → class is well-defined and
    /// the result is finite in [0, 1] (never ±Inf).
    #[test]
    fn dispersion_entropy_variance_overflow_is_finite() {
        let w: Vec<f64> = (0..BAR_CLOSE_LOOKBACK_COUNT)
            .map(|i| if i % 2 == 0 { 1e200 } else { -1e200 })
            .collect();
        let v = compute_bar_dispersion_entropy(&w);
        assert!(
            !v.is_infinite() && (v.is_nan() || (0.0..=1.0).contains(&v)),
            "must be NaN-or-[0,1], never ±Inf, got {v}"
        );
    }

    // ── bar_cecp_velocity ────────────────────────────────────────────────────

    /// A flat full window has `ptp == 0` → CECP velocity undefined → NaN (mirrors
    /// the SSoT `np.ptp(x) == 0.0` guard).
    #[test]
    fn cecp_velocity_flat_window_is_nan() {
        assert!(compute_bar_cecp_velocity(&[42.0; BAR_CLOSE_LOOKBACK_COUNT]).is_nan());
    }

    /// Degenerate inputs are undefined → NaN: too short (`< 2·(m+2) = 10`) and
    /// non-finite values both hit the leading guards.
    #[test]
    fn cecp_velocity_degenerate_inputs_are_nan() {
        assert!(
            compute_bar_cecp_velocity(&[1.0, 2.0, 3.0]).is_nan(),
            "len < 10"
        );
        let mut w = vec![1.0_f64; BAR_CLOSE_LOOKBACK_COUNT];
        w[100] = f64::NAN;
        assert!(compute_bar_cecp_velocity(&w).is_nan(), "non-finite input");
    }

    /// On a varied window the velocity is finite and within `[0, √2]` (both H and
    /// C live in `[0, 1]`, so their displacement cannot exceed √2). Stage-3 adds
    /// the dedicated invariant-5 bounds test; this is the Stage-1 sanity guard.
    #[test]
    fn cecp_velocity_bounded_on_varied_window() {
        let w: Vec<f64> = (0..BAR_CLOSE_LOOKBACK_COUNT)
            .map(|i| 100.0 + (i as f64 * 0.7).sin() * 3.0 + i as f64 * 0.01)
            .collect();
        let v = compute_bar_cecp_velocity(&w);
        assert!(
            v.is_finite() && (0.0..=std::f64::consts::SQRT_2).contains(&v),
            "out of [0, √2]: {v}"
        );
    }

    /// A window whose two halves are IDENTICAL must yield velocity 0 exactly: the
    /// (H, C) coordinate is the same on both halves, so the displacement is 0.
    /// (`0.0` is a valid value, never the NaN sentinel.)
    #[test]
    fn cecp_velocity_identical_halves_is_zero() {
        // Per-half varied (so each half is non-degenerate) but the two halves are
        // bit-identical; the full window is not flat, so the guard passes.
        let mut w = Vec::with_capacity(BAR_CLOSE_LOOKBACK_COUNT);
        let half: Vec<f64> = (0..CECP_HALF)
            .map(|i| 100.0 + (i as f64 * 0.5).sin())
            .collect();
        w.extend_from_slice(&half);
        w.extend_from_slice(&half);
        assert_eq!(w.len(), BAR_CLOSE_LOOKBACK_COUNT);
        let v = compute_bar_cecp_velocity(&w);
        assert_eq!(v, 0.0, "identical halves → zero displacement, got {v}");
    }

    /// Invariant 5 (the deeper bounds, beyond the velocity range): each (H, C)
    /// coordinate ordpy produces must itself lie in the unit square — H =
    /// normalised permutation entropy ∈ [0, 1] and C = MPR statistical complexity
    /// ∈ [0, 1]. This is what makes the velocity bound √2 sound, and it catches a
    /// Q0 / js_div_max / JS-summation error the specific oracle sample might miss
    /// (e.g. a wrong normalizer would push C outside [0, 1] long before the
    /// velocity left [0, √2]). Exercised over diverse segment shapes incl. the
    /// degenerate constant half (→ (0, 0)) and a near-uniform-ordinal half
    /// (H → 1, C → 0).
    #[test]
    fn cecp_hc_coordinates_in_unit_square() {
        let constant: Vec<f64> = vec![42.0; CECP_HALF];
        let trend: Vec<f64> = (0..CECP_HALF).map(|i| i as f64).collect();
        let wave: Vec<f64> = (0..CECP_HALF)
            .map(|i| 100.0 + (i as f64 * 0.7).sin() * 3.0 + i as f64 * 0.01)
            .collect();
        let noisy: Vec<f64> = (0..CECP_HALF)
            .map(|i| ((i * 2_654_435_761usize) % 1000) as f64) // cheap deterministic pseudo-noise
            .collect();
        for seg in [&constant, &trend, &wave, &noisy] {
            let (h, c) = cecp_hc(seg);
            assert!(
                h.is_finite() && (0.0..=1.0).contains(&h),
                "H out of [0,1]: {h}"
            );
            assert!(
                c.is_finite() && (0.0..=1.0).contains(&c),
                "C out of [0,1]: {c}"
            );
        }
        // The constant half is the exact corner (H, C) == (0, 0).
        assert_eq!(cecp_hc(&constant), (0.0, 0.0));
    }

    /// Machine-checked proof that the reused `ordinal_pattern_index_m3` (`<=` tie
    /// convention) classes ordinal triplets into the SAME equivalence partition as
    /// numpy `np.argsort` — i.e. ordpy's `ordinal_sequence` symbol. Two triplets
    /// share an ordpy symbol IFF they share a Rust class. Since H and C are
    /// invariant under pattern relabeling, this partition-equivalence is exactly
    /// what makes the `(H, C)` pair bit-identical to ordpy, INCLUDING ties. The
    /// `ARGSORT_REF` table is `np.argsort([a,b,c])` for every `(a,b,c)` in
    /// `{0,1,2}^3` (indexed `9a+3b+c`), generated once from numpy 1.x and committed
    /// here; `{0,1,2}^3` covers every possible relative-order configuration of three
    /// values (all distinct, every single tie, and the all-equal case), so it is an
    /// EXHAUSTIVE proof, not a sample. (Replaces the prose "PROVEN" claim with an
    /// enforced gate.)
    #[test]
    fn cecp_ordinal_classing_matches_argsort_over_all_triplets() {
        // np.argsort([a,b,c]) for (a,b,c) in {0,1,2}^3, idx = 9a + 3b + c.
        #[rustfmt::skip]
        const ARGSORT_REF: [[u8; 3]; 27] = [
            [0,1,2],[0,1,2],[0,1,2],[0,2,1],[0,1,2],[0,1,2],[0,2,1],[0,2,1],[0,1,2],
            [1,2,0],[1,0,2],[1,0,2],[2,0,1],[0,1,2],[0,1,2],[2,0,1],[0,2,1],[0,1,2],
            [1,2,0],[1,2,0],[1,0,2],[2,1,0],[1,2,0],[1,0,2],[2,0,1],[2,0,1],[0,1,2],
        ];
        let mut triplets = Vec::with_capacity(27);
        for a in 0..3u8 {
            for b in 0..3u8 {
                for c in 0..3u8 {
                    triplets.push((a, b, c));
                }
            }
        }
        assert_eq!(triplets.len(), ARGSORT_REF.len());
        // For EVERY ordered pair of triplets: the Rust classing puts them in the
        // same class IFF numpy argsort gives them the same symbol.
        for i in 0..triplets.len() {
            for j in 0..triplets.len() {
                let (ai, bi, ci) = triplets[i];
                let (aj, bj, cj) = triplets[j];
                let rust_same = ordinal_pattern_index_m3(ai as f64, bi as f64, ci as f64)
                    == ordinal_pattern_index_m3(aj as f64, bj as f64, cj as f64);
                let argsort_same = ARGSORT_REF[i] == ARGSORT_REF[j];
                assert_eq!(
                    rust_same, argsort_same,
                    "classing partition mismatch: {:?} vs {:?} (rust_same={rust_same}, argsort_same={argsort_same})",
                    triplets[i], triplets[j]
                );
            }
        }
    }

    // ── bar_categorical_recurrence_rate (recurr, card 66) ─────────────────────
    // The property battery below IS this feature's independent leg: the estimator
    // is an exact closed form on integer tie-counts, so a separate brute-force
    // would be the identical arithmetic. These hand cases + invariants pin it.

    /// Exact RR == num/den computed the SAME way the kernel divides (bit-identical).
    fn rr_eq(got: f64, num: u64, den: u64) -> bool {
        got.to_bits() == (num as f64 / den as f64).to_bits()
    }

    #[test]
    fn recurr_all_distinct_is_zero() {
        let w: Vec<f64> = (1..=BAR_CLOSE_LOOKBACK_COUNT).map(|i| i as f64).collect();
        assert!(rr_eq(compute_bar_categorical_recurrence_rate(&w), 0, 39800));
    }

    #[test]
    fn recurr_all_equal_is_one() {
        let w = [42.0_f64; BAR_CLOSE_LOOKBACK_COUNT];
        // tie_sum = N(N-1) = denom → RR = 1.0 exactly.
        assert!(rr_eq(
            compute_bar_categorical_recurrence_rate(&w),
            39800,
            39800
        ));
    }

    #[test]
    fn recurr_hand_three_tie() {
        // One level appears 3× (m=3 → 3·2=6), two singletons (0). N=5, denom=20.
        let w = [1.0, 1.0, 1.0, 2.0, 3.0];
        assert!(rr_eq(compute_bar_categorical_recurrence_rate(&w), 6, 20));
    }

    #[test]
    fn recurr_two_pairs() {
        // [1,1,2,2]: two levels ×2 → tie_sum = 2+2 = 4. N=4, denom=12.
        let w = [1.0, 1.0, 2.0, 2.0];
        assert!(rr_eq(compute_bar_categorical_recurrence_rate(&w), 4, 12));
    }

    #[test]
    fn recurr_mixed_multiplicities_in_unit_range() {
        // [1,2,2,3,3,3,4]: m = 1,2,3,1 → tie_sum = 0+2+6+0 = 8. N=7, denom=42.
        let w = [1.0, 2.0, 2.0, 3.0, 3.0, 3.0, 4.0];
        let got = compute_bar_categorical_recurrence_rate(&w);
        assert!((0.0..=1.0).contains(&got));
        assert!(rr_eq(got, 8, 42));
    }

    #[test]
    fn recurr_nan_on_short_or_empty() {
        assert!(compute_bar_categorical_recurrence_rate(&[]).is_nan());
        assert!(compute_bar_categorical_recurrence_rate(&[5.0]).is_nan());
    }

    #[test]
    fn recurr_nan_on_nonfinite() {
        assert!(compute_bar_categorical_recurrence_rate(&[1.0, f64::NAN, 2.0]).is_nan());
        assert!(compute_bar_categorical_recurrence_rate(&[1.0, f64::INFINITY, 2.0]).is_nan());
        assert!(compute_bar_categorical_recurrence_rate(&[f64::NEG_INFINITY, 1.0]).is_nan());
    }

    #[test]
    fn recurr_never_inf_on_degenerate() {
        for w in [vec![0.0_f64; 200], vec![-3.5_f64; 200], vec![7.0_f64; 2]] {
            assert!(!compute_bar_categorical_recurrence_rate(&w).is_infinite());
        }
    }

    #[test]
    fn recurr_tie_preserving_under_constant_rescale() {
        // scale_check: a constant rescale preserves exact ties AND distinctness, so
        // RR is bit-identical (rank/tie feature — invariant, not just within tol).
        let base = [1.0, 1.0, 2.0, 3.0, 3.0, 3.0];
        for factor in [2.0_f64, 0.5, 100.0] {
            let scaled: Vec<f64> = base.iter().map(|&x| x * factor).collect();
            assert_eq!(
                compute_bar_categorical_recurrence_rate(&base).to_bits(),
                compute_bar_categorical_recurrence_rate(&scaled).to_bits(),
                "RR not tie-preserving under ×{factor}"
            );
        }
    }

    // ── bar_sign_markov_flux (signflux, card 98) ──────────────────────────────
    // The property battery + the in-Rust brute-force transition double-loop (in
    // bar_close_robustness_test.rs) are this feature's independent leg.

    #[test]
    fn signflux_zero_on_monotone_series() {
        // Strictly increasing → every sign is `up` → only the 2→2 self-transition
        // occurs → no off-diagonal current → flux = 0 (detailed-balance limit).
        let w = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        assert_eq!(
            compute_bar_sign_markov_flux(&w).to_bits(),
            0.0_f64.to_bits()
        );
    }

    #[test]
    fn signflux_one_on_perfect_three_cycle() {
        // up,flat,down repeating: every transition is one-directional around the
        // 3-cycle (up→flat→down→up), so each off-diagonal pair has a zero reverse
        // leg → Σ|C_ij−C_ji| == total → flux = 1.
        let w = [1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.0];
        assert_eq!(
            compute_bar_sign_markov_flux(&w).to_bits(),
            1.0_f64.to_bits()
        );
    }

    #[test]
    fn signflux_in_unit_range_on_mixed_window() {
        let w = [1.0, 2.0, 1.5, 1.5, 3.0, 2.0, 2.5, 2.5, 1.0, 4.0];
        let got = compute_bar_sign_markov_flux(&w);
        assert!((0.0..=1.0).contains(&got), "flux {got} out of [0,1]");
    }

    #[test]
    fn signflux_nan_on_short_or_undefined() {
        // < 4 signs (n < 5) → NaN.
        assert!(compute_bar_sign_markov_flux(&[1.0, 2.0, 3.0]).is_nan());
        assert!(compute_bar_sign_markov_flux(&[]).is_nan());
        // any close ≤ 0 (ln undefined) → NaN.
        assert!(compute_bar_sign_markov_flux(&[1.0, 2.0, 0.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_sign_markov_flux(&[1.0, -2.0, 3.0, 4.0, 5.0]).is_nan());
        // non-finite close → NaN.
        assert!(compute_bar_sign_markov_flux(&[1.0, f64::NAN, 3.0, 4.0, 5.0]).is_nan());
        assert!(compute_bar_sign_markov_flux(&[1.0, f64::INFINITY, 3.0, 4.0, 5.0]).is_nan());
    }

    #[test]
    fn signflux_never_inf_on_degenerate() {
        // Flat full window → all flat signs → only 1→1 self-transition → flux 0.
        let flat = vec![5.0_f64; 200];
        let v = compute_bar_sign_markov_flux(&flat);
        assert!(!v.is_infinite());
        assert_eq!(v.to_bits(), 0.0_f64.to_bits());
        for w in [vec![7.0_f64; 2], vec![1.0, 2.0, 3.0, 4.0]] {
            assert!(!compute_bar_sign_markov_flux(&w).is_infinite());
        }
    }

    #[test]
    fn signflux_invariant_under_positive_rescale() {
        // scale_check: a positive constant rescale preserves the sign sequence
        // exactly (sign(ln(k·a) − ln(k·b)) = sign(ln a − ln b)), so flux is
        // bit-identical — sign-sequence invariant, not merely within tol.
        let base = [1.0, 2.0, 2.0, 1.0, 3.0, 2.5, 2.5, 4.0];
        for factor in [2.0_f64, 0.5, 1000.0] {
            let scaled: Vec<f64> = base.iter().map(|&x| x * factor).collect();
            assert_eq!(
                compute_bar_sign_markov_flux(&base).to_bits(),
                compute_bar_sign_markov_flux(&scaled).to_bits(),
                "sign-flux not invariant under ×{factor}"
            );
        }
    }

    // ── bar_ramsey_rothman_bicov_lag1 (rrbicov, card 113) ─────────────────────
    // CONTINUOUS estimator (ln + std + cubing): the gate is <=1e-9, NOT bit-exact.
    // The property battery (time-reversal antisymmetry, scale-freedom, NaN paths) is
    // this feature's lightweight check; the in-Rust brute-force independent leg lives
    // in bar_close_robustness_test.rs.

    /// Build a strictly-positive close series whose `diff(ln(close))` equals `logret`
    /// (`c[0]=100; c[i+1]=c[i]·exp(logret[i])`), for the antisymmetry property test.
    fn closes_from_logret(logret: &[f64]) -> Vec<f64> {
        let mut c = Vec::with_capacity(logret.len() + 1);
        c.push(100.0_f64);
        for &r in logret {
            let prev = *c.last().unwrap();
            c.push(prev * libm::exp(r));
        }
        c
    }

    #[test]
    fn rrbicov_time_reversal_antisymmetry() {
        // Reversing the logret series negates the bicovariance (time-irreversibility).
        // Continuous estimator → assert the sum is ~0 within a tight FP tolerance.
        let logret = [
            0.10, 0.04, -0.18, 0.06, -0.02, 0.09, -0.15, 0.03, 0.07, -0.05, 0.12, -0.08,
        ];
        let fwd = closes_from_logret(&logret);
        let mut rev_lr = logret;
        rev_lr.reverse();
        let rev = closes_from_logret(&rev_lr);
        let v_fwd = compute_bar_ramsey_rothman_bicov_lag1(&fwd);
        let v_rev = compute_bar_ramsey_rothman_bicov_lag1(&rev);
        assert!(v_fwd.is_finite() && v_rev.is_finite());
        assert!(
            v_fwd.abs() > 1e-6,
            "asymmetric series must be non-trivially irreversible (got {v_fwd})"
        );
        assert!(
            (v_fwd + v_rev).abs() < 1e-12,
            "time-reversal antisymmetry breach: {v_fwd} + {v_rev} != 0"
        );
    }

    #[test]
    fn rrbicov_nan_on_short_or_undefined() {
        // < 3 logrets (n < 4) → NaN.
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[1.0, 2.0, 3.0]).is_nan());
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[]).is_nan());
        // any close ≤ 0 (ln undefined) → NaN.
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[1.0, 2.0, 0.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[1.0, -2.0, 3.0, 4.0, 5.0]).is_nan());
        // non-finite close → NaN.
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[1.0, f64::NAN, 3.0, 4.0, 5.0]).is_nan());
        // flat window → sigma == 0 → NaN.
        assert!(compute_bar_ramsey_rothman_bicov_lag1(&[5.0; 200]).is_nan());
    }

    #[test]
    fn rrbicov_never_inf_on_degenerate() {
        for w in [
            vec![5.0_f64; 200],
            vec![7.0_f64; 4],
            vec![1.0, 2.0, 3.0],
            vec![1.0, 1.0, 1.0, 1.0, 2.0],
        ] {
            assert!(!compute_bar_ramsey_rothman_bicov_lag1(&w).is_infinite());
        }
    }

    // ── bar_ehlers_increment_asymmetry (ehlers, card 26) ──────────────────────
    // CONTINUOUS estimator (cube/square/pow): the gate is <=1e-9, NOT bit-exact.
    // The property battery (scale-invariance, odd-negation, =0 symmetric, NaN flat)
    // IS this feature's independent leg — a brute-force would be the identical
    // closed-form power-sum arithmetic.

    /// Build a close series whose raw first differences equal `incr`
    /// (`c[0]=100; c[i+1]=c[i]+incr[i]`), for the property tests.
    fn closes_from_increments(incr: &[f64]) -> Vec<f64> {
        let mut c = Vec::with_capacity(incr.len() + 1);
        c.push(100.0_f64);
        for &d in incr {
            let prev = *c.last().unwrap();
            c.push(prev + d);
        }
        c
    }

    #[test]
    fn ehlers_zero_on_symmetric_increments() {
        // A symmetric increment multiset (every +d has a matching −d) → Σd^3 = 0 → E = 0.
        let closes = closes_from_increments(&[1.0, -1.0, 2.0, -2.0, 3.0, -3.0]);
        assert_eq!(
            compute_bar_ehlers_increment_asymmetry(&closes).to_bits(),
            0.0_f64.to_bits()
        );
    }

    #[test]
    fn ehlers_odd_under_increment_negation() {
        // Negating the increments negates E (Σ(−d)^3 = −Σd^3, Σ(−d)^2 = Σd^2) — exact
        // in FP (negation distributes over the sums and the quotient).
        let incr = [0.7, -0.2, 0.5, 0.9, -0.4, 0.3, -0.8, 0.6];
        let fwd = closes_from_increments(&incr);
        let neg_incr: Vec<f64> = incr.iter().map(|&d| -d).collect();
        let neg = closes_from_increments(&neg_incr);
        let e = compute_bar_ehlers_increment_asymmetry(&fwd);
        let e_neg = compute_bar_ehlers_increment_asymmetry(&neg);
        assert!(
            e.is_finite() && e.abs() > 1e-6,
            "need a non-trivially skewed window"
        );
        assert_eq!(e_neg.to_bits(), (-e).to_bits(), "E not odd under negation");
    }

    #[test]
    fn ehlers_nan_on_short_or_flat() {
        // < 3 closes → NaN.
        assert!(compute_bar_ehlers_increment_asymmetry(&[1.0, 2.0]).is_nan());
        assert!(compute_bar_ehlers_increment_asymmetry(&[]).is_nan());
        // flat window → sum(d^2) == 0 → NaN.
        assert!(compute_bar_ehlers_increment_asymmetry(&[5.0; 200]).is_nan());
        // non-finite close → NaN.
        assert!(compute_bar_ehlers_increment_asymmetry(&[1.0, f64::NAN, 3.0]).is_nan());
        assert!(compute_bar_ehlers_increment_asymmetry(&[1.0, 2.0, f64::INFINITY]).is_nan());
    }

    #[test]
    fn ehlers_never_inf_on_degenerate() {
        for w in [
            vec![5.0_f64; 200],
            vec![1.0, 2.0, 3.0],
            vec![1.0, 1.0, 1.0, 1.0, 5.0],
        ] {
            assert!(!compute_bar_ehlers_increment_asymmetry(&w).is_infinite());
        }
    }

    // ── bar_cox_stuart_trend_z (coxstuart, card 23) ───────────────────────────
    // BIT-EXACT (integer S+ count -> (S+-m/2)/sqrt(m/4)). The property battery +
    // the in-Rust brute-force (in bar_close_robustness_test.rs) are the independent
    // leg: the FOSS randtests cox.stuart.test is non-hermetic (no R) AND may use a
    // different sign-test normalization than the card's pinned formula.

    #[test]
    fn coxstuart_zero_on_trendless_pairing() {
        // n=4, m=2: pairs (x0,x2),(x1,x3). One pair up, one not → S+=1=m/2 → z=0.
        let w = [1.0, 2.0, 3.0, 1.0]; // (1<3 up), (2>1 not) → S+=1
        assert_eq!(
            compute_bar_cox_stuart_trend_z(&w).to_bits(),
            0.0_f64.to_bits()
        );
    }

    #[test]
    fn coxstuart_max_positive_on_monotone_up() {
        // Strictly increasing 200-window → every second>first → S+=100 → z=(100-50)/5=10.
        let w: Vec<f64> = (1..=BAR_CLOSE_LOOKBACK_COUNT).map(|i| i as f64).collect();
        assert_eq!(
            compute_bar_cox_stuart_trend_z(&w).to_bits(),
            10.0_f64.to_bits()
        );
    }

    #[test]
    fn coxstuart_max_negative_on_monotone_down() {
        // Strictly decreasing → no second>first → S+=0 → z=(0-50)/5=-10.
        let w: Vec<f64> = (1..=BAR_CLOSE_LOOKBACK_COUNT)
            .rev()
            .map(|i| i as f64)
            .collect();
        assert_eq!(
            compute_bar_cox_stuart_trend_z(&w).to_bits(),
            (-10.0_f64).to_bits()
        );
    }

    #[test]
    fn coxstuart_sign_tracks_trend_direction() {
        let up = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let down = [6.0, 5.0, 4.0, 3.0, 2.0, 1.0];
        assert!(compute_bar_cox_stuart_trend_z(&up) > 0.0);
        assert!(compute_bar_cox_stuart_trend_z(&down) < 0.0);
    }

    #[test]
    fn coxstuart_invariant_under_positive_rescale() {
        // The pairwise comparison is sign-preserving under a positive constant, and
        // the count is an integer → z is bit-identical (not just within tol).
        let base = [1.0, 5.0, 2.0, 9.0, 3.0, 1.0, 4.0, 8.0, 2.5, 7.0];
        for factor in [2.0_f64, 0.5, 1000.0, 0.001] {
            let scaled: Vec<f64> = base.iter().map(|&x| x * factor).collect();
            assert_eq!(
                compute_bar_cox_stuart_trend_z(&base).to_bits(),
                compute_bar_cox_stuart_trend_z(&scaled).to_bits(),
                "cox-stuart z not invariant under ×{factor}"
            );
        }
    }

    #[test]
    fn coxstuart_nan_on_short_or_nonfinite() {
        assert!(compute_bar_cox_stuart_trend_z(&[1.0, 2.0, 3.0]).is_nan()); // n<4
        assert!(compute_bar_cox_stuart_trend_z(&[]).is_nan());
        assert!(compute_bar_cox_stuart_trend_z(&[1.0, f64::NAN, 3.0, 4.0]).is_nan());
        assert!(compute_bar_cox_stuart_trend_z(&[1.0, 2.0, f64::INFINITY, 4.0]).is_nan());
        // A flat window is DEFINED (z = -(m/2)/sqrt(m/4)), not NaN.
        assert!(compute_bar_cox_stuart_trend_z(&[5.0; 200]).is_finite());
    }

    #[test]
    fn coxstuart_never_inf_on_degenerate() {
        for w in [
            vec![5.0_f64; 200],
            vec![1.0, 2.0, 3.0, 4.0],
            vec![1.0, 2.0, 3.0],
        ] {
            assert!(!compute_bar_cox_stuart_trend_z(&w).is_infinite());
        }
    }

    // ── bar_groeneveld_meeden_b3_skewness (gmskew, card 99) ───────────────────
    // CONTINUOUS estimator (ln + median + MAD): the gate is <=1e-9, NOT bit-exact.
    // The property battery + the in-Rust brute-force (sort+median+MAD, in
    // bar_close_robustness_test.rs) are the independent leg (no faithful FOSS).

    #[test]
    fn gmskew_zero_on_symmetric_logrets() {
        // Logrets symmetric about their median → mean == median → b3 = 0 (within the
        // exp/ln round-trip + FP noise of a continuous estimator).
        let logret = [-0.02, -0.01, 0.01, 0.02]; // sorted; median=0, mean=0
        let closes = closes_from_logret(&logret);
        let b3 = compute_bar_groeneveld_meeden_b3_skewness(&closes);
        assert!(
            b3.is_finite() && b3.abs() <= 1e-9,
            "symmetric b3 should be ~0, got {b3}"
        );
    }

    #[test]
    fn gmskew_antisymmetric_under_logret_negation() {
        // Negating the logrets negates b3 (mean→−mean, med→−med, MAD unchanged).
        let logret = [
            0.03, -0.01, 0.05, 0.12, -0.02, 0.04, -0.08, 0.01, 0.06, -0.03,
        ];
        let fwd = closes_from_logret(&logret);
        let neg: Vec<f64> = logret.iter().map(|&r| -r).collect();
        let rev = closes_from_logret(&neg);
        let b = compute_bar_groeneveld_meeden_b3_skewness(&fwd);
        let bn = compute_bar_groeneveld_meeden_b3_skewness(&rev);
        assert!(
            b.is_finite() && b.abs() > 1e-6,
            "need a non-trivially skewed window"
        );
        assert!(
            (b + bn).abs() <= 1e-9,
            "b3 not antisymmetric: {b} + {bn} != 0"
        );
    }

    #[test]
    fn gmskew_in_unit_range_on_mixed_window() {
        // b3 ∈ [-1,1] is a mathematical property of the GM estimator for ANY data.
        let logret = [
            0.03, -0.01, 0.05, 0.12, -0.02, 0.04, -0.08, 0.01, 0.06, -0.03, 0.09, -0.04, 0.02,
            -0.07, 0.11,
        ];
        let closes = closes_from_logret(&logret);
        let b3 = compute_bar_groeneveld_meeden_b3_skewness(&closes);
        assert!(
            b3.is_finite() && (-1.0..=1.0).contains(&b3),
            "b3 {b3} out of [-1,1]"
        );
    }

    #[test]
    fn gmskew_nan_on_short_or_undefined() {
        // < 2 logrets (n < 3) → NaN.
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[1.0, 2.0]).is_nan());
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[]).is_nan());
        // close ≤ 0 (ln undefined) → NaN.
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[1.0, 0.0, 2.0, 3.0]).is_nan());
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[1.0, -2.0, 3.0]).is_nan());
        // non-finite close → NaN.
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[1.0, f64::NAN, 3.0]).is_nan());
        // constant window → all logret 0 → MAD == 0 → NaN.
        assert!(compute_bar_groeneveld_meeden_b3_skewness(&[5.0; 200]).is_nan());
    }

    #[test]
    fn gmskew_never_inf_on_degenerate() {
        for w in [
            vec![5.0_f64; 200],
            vec![1.0, 2.0, 3.0],
            vec![1.0, 1.0, 1.0, 1.0, 2.0],
        ] {
            assert!(!compute_bar_groeneveld_meeden_b3_skewness(&w).is_infinite());
        }
    }

    #[test]
    fn gmskew_even_n_median_averaging_path() {
        // Explicit coverage of the EVEN-n median path (mean of the 2 central order
        // statistics). 5 closes → 4 logrets (EVEN). Hand-built skewed logrets:
        // sorted [0.01,0.02,0.03,0.10] → median=(0.02+0.03)/2=0.025; mean=0.04;
        // MAD=(0.075+0.015+0.005+0.005)/4=0.025; b3=(0.04-0.025)/0.025=0.6.
        let closes = closes_from_logret(&[0.10, 0.01, 0.02, 0.03]);
        assert_eq!(closes.len(), 5, "5 closes → 4 logrets (even median path)");
        let b3 = compute_bar_groeneveld_meeden_b3_skewness(&closes);
        assert!(b3.is_finite() && (-1.0..=1.0).contains(&b3));
        assert!(
            (b3 - 0.6).abs() <= 1e-9,
            "even-n median-averaged b3 should be ~0.6, got {b3}"
        );
    }

    // ── bar_l_kurtosis_tau4 (lkurt, card 2) ───────────────────────────────────
    // L-kurtosis ratio tau4 = l4/l2 via probability-weighted moments (Hosking 1990).

    /// Deterministic ~N(0,1) draw via an LCG + Irwin–Hall CLT (sum of 12 uniforms − 6).
    /// Fixed seed → fully reproducible; no flakiness, no `rand` dependency.
    fn normalish_logrets(count: usize, seed: u64) -> Vec<f64> {
        let mut state = seed;
        let mut next_u = || -> f64 {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            // top 53 bits → [0, 1)
            (state >> 11) as f64 / (1u64 << 53) as f64
        };
        (0..count)
            .map(|_| {
                let s: f64 = (0..12).map(|_| next_u()).sum();
                0.01 * (s - 6.0) // ~N(0, 0.01^2): realistic logret magnitude
            })
            .collect()
    }

    #[test]
    fn lkurt_near_normal_tau4_approx_0p1226() {
        // For a (near-)Gaussian return window the population L-kurtosis is tau4 ≈ 0.1226.
        // 200 logrets → 201 closes; deterministic CLT-normal sample. Finite-sample tau4
        // scatters around 0.1226, so assert a generous bracket (NOT the 1e-9 oracle gate —
        // that lives in the parity test against the frozen fixture).
        let lr = normalish_logrets(200, 0x9E3779B97F4A7C15);
        let closes = closes_from_logret(&lr);
        assert_eq!(closes.len(), 201);
        let t4 = compute_bar_l_kurtosis_tau4(&closes);
        assert!(t4.is_finite(), "near-normal tau4 must be finite, got {t4}");
        assert!(
            (0.05..=0.20).contains(&t4),
            "near-normal tau4 should sit near the Gaussian 0.1226, got {t4}"
        );
    }

    #[test]
    fn lkurt_nan_on_short_or_undefined() {
        // < 4 logrets (n_closes < 5): PWM b3 undefined.
        assert!(compute_bar_l_kurtosis_tau4(&[1.0, 2.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_l_kurtosis_tau4(&[]).is_nan());
        // close <= 0 / non-finite → ln undefined.
        assert!(compute_bar_l_kurtosis_tau4(&[1.0, 0.0, 2.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_l_kurtosis_tau4(&[1.0, -2.0, 3.0, 4.0, 5.0]).is_nan());
        assert!(compute_bar_l_kurtosis_tau4(&[1.0, f64::NAN, 3.0, 4.0, 5.0]).is_nan());
        // constant window → all logret 0 → l2 == 0.
        assert!(compute_bar_l_kurtosis_tau4(&[5.0; 200]).is_nan());
    }

    #[test]
    fn lkurt_never_inf_on_degenerate() {
        // Near-constant + extreme-ratio windows must yield NaN or finite, NEVER ±Inf.
        let cases: Vec<Vec<f64>> = vec![
            vec![1.0, 1.0, 1.0, 1.0, 1.000_000_001],
            vec![1e-12, 1.0, 1e12, 1.0, 1e-12, 1.0],
            vec![1.0; 5],
        ];
        for w in cases {
            assert!(!compute_bar_l_kurtosis_tau4(&w).is_infinite());
        }
    }

    #[test]
    fn lkurt_finite_in_sample_range_on_mixed_window() {
        // A real-ish mixed window stays in the valid L-kurtosis band. SAMPLE tau4 can dip
        // a touch below the population −0.25 floor, so use a generous lower bound.
        let lr = normalish_logrets(120, 0x123456789ABCDEF0);
        let closes = closes_from_logret(&lr);
        let t4 = compute_bar_l_kurtosis_tau4(&closes);
        assert!(t4.is_finite(), "tau4 must be finite, got {t4}");
        assert!(
            (-0.30..=1.0).contains(&t4),
            "tau4 out of valid L-kurtosis band: {t4}"
        );
    }

    // ── bar_bartels_rank_vn_ratio (bartels, card 19) ──────────────────────────
    // Rank von-Neumann ratio z = (RVN-2)/sqrt(sigma2); RVN = Σdiff(R)²/Σ(R-meanR)²,
    // R = avg-rank(close[i+1]/close[i]); sigma2 = Bartels (1982) exact variance.

    #[test]
    fn bartels_average_ranks_match_scipy_convention() {
        // ties averaged, 1-indexed: [3.0, 1.0, 1.0, 5.0] (2nd & 3rd tie) →
        // sorted positions: the two 1.0s occupy positions 1,2 → avg 1.5; 3.0 → 3; 5.0 → 4.
        let r = bar_close_average_ranks(&[3.0, 1.0, 1.0, 5.0]);
        assert_eq!(r, vec![3.0, 1.5, 1.5, 4.0]);
        // all-distinct, ascending input → identity ranks.
        let r2 = bar_close_average_ranks(&[10.0, 20.0, 30.0]);
        assert_eq!(r2, vec![1.0, 2.0, 3.0]);
        // all-equal → all get the mean rank (n+1)/2.
        let r3 = bar_close_average_ranks(&[7.0, 7.0, 7.0, 7.0]);
        assert_eq!(r3, vec![2.5, 2.5, 2.5, 2.5]);
    }

    #[test]
    fn bartels_hand_pinned_small_window() {
        // 5 closes → 4 logrets. closes=[100,101,103,106,102]:
        // logret ≈ [+,+,+,−]; ascending order → ranks (time order) R=[2,3,4,1].
        // diff(R)=[1,1,-3] → num=1+1+9=11; meanR=2.5 → denom=0.25+0.25+2.25+2.25=5;
        // RVN=11/5=2.2; sigma2=4·2·(80-8-9)/(5·4·5·9)=504/900=0.56;
        // z=(2.2-2)/sqrt(0.56)=0.2/0.748331...=0.2672612419124244.
        let closes = [100.0_f64, 101.0, 103.0, 106.0, 102.0];
        let z = compute_bar_bartels_rank_vn_ratio(&closes);
        assert!(
            (z - 0.267_261_241_912_424_4).abs() <= 1e-12,
            "hand-pinned bartels z mismatch: got {z}"
        );
    }

    #[test]
    fn bartels_strong_negative_z_on_monotone_series() {
        // A strictly-increasing logret series → ranks 1..n in time order → RVN→0 →
        // z ≈ −sqrt(n) (perfect positive serial correlation / momentum). Strongly negative.
        let lr: Vec<f64> = (0..150).map(|i| 0.001 + 0.0001 * i as f64).collect();
        let closes = closes_from_logret(&lr);
        let z = compute_bar_bartels_rank_vn_ratio(&closes);
        assert!(
            z.is_finite() && z < -3.0,
            "monotone series should give strong negative z, got {z}"
        );
    }

    #[test]
    fn bartels_scale_invariant_rank_feature() {
        // Pure rank statistic: scaling closes by k>0 preserves the logret ORDER, so the
        // ranks — and z — are unchanged (within FP; ranks identical → z bit-identical).
        let base = closes_from_logret(&normalish_logrets(80, 0xA5A5A5A55A5A5A5A));
        let z0 = compute_bar_bartels_rank_vn_ratio(&base);
        assert!(z0.is_finite());
        for factor in [0.5_f64, 2.0, 13.0, 1e-4, 9999.0] {
            let scaled: Vec<f64> = base.iter().map(|&c| c * factor).collect();
            let z = compute_bar_bartels_rank_vn_ratio(&scaled);
            assert!(
                (z0 - z).abs() <= 1e-9,
                "bartels not rank/scale-invariant under ×{factor}: {z0} vs {z}"
            );
        }
    }

    #[test]
    fn bartels_random_series_is_finite_and_sane() {
        // A pseudo-random logret series should yield a finite z of N(0,1) scale.
        let lr = normalish_logrets(199, 0x0123456789ABCDEF);
        let closes = closes_from_logret(&lr);
        let z = compute_bar_bartels_rank_vn_ratio(&closes);
        assert!(
            z.is_finite() && z.abs() < 8.0,
            "random-series bartels z out of sane band: {z}"
        );
    }

    #[test]
    fn bartels_nan_on_short_or_undefined() {
        // < 4 logrets (n_closes < 5).
        assert!(compute_bar_bartels_rank_vn_ratio(&[1.0, 2.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_bartels_rank_vn_ratio(&[]).is_nan());
        // close <= 0 / non-finite → ln undefined.
        assert!(compute_bar_bartels_rank_vn_ratio(&[1.0, 0.0, 2.0, 3.0, 4.0]).is_nan());
        assert!(compute_bar_bartels_rank_vn_ratio(&[1.0, -2.0, 3.0, 4.0, 5.0]).is_nan());
        assert!(compute_bar_bartels_rank_vn_ratio(&[1.0, f64::NAN, 3.0, 4.0, 5.0]).is_nan());
        // constant window → all logret 0 → all ranks tie → rank variance 0.
        assert!(compute_bar_bartels_rank_vn_ratio(&[5.0; 200]).is_nan());
    }

    #[test]
    fn bartels_never_inf_on_degenerate() {
        let cases: Vec<Vec<f64>> = vec![
            vec![1.0, 1.0, 1.0, 1.0, 1.000_000_001],
            vec![1e-9, 1.0, 1e9, 1.0, 1e-9, 1.0],
            vec![1.0; 5],
        ];
        for w in cases {
            assert!(!compute_bar_bartels_rank_vn_ratio(&w).is_infinite());
        }
    }

    // ── bar_hoeffding_phi_squared_midreturn_duration (hoeffding, card 78) ────────

    /// INDEPENDENT brute-force leg: the copula grid by DIRECT O(m³) counting —
    /// no indicator matrix, no prefix sums — so a bookkeeping bug in the kernel's
    /// cumsum construction cannot be self-confirming.
    fn hoeffding_phi2_brute(closes: &[f64], durations_us: &[f64]) -> f64 {
        let n = closes.len();
        assert!(n >= 5 && durations_us.len() == n);
        let m = n - 1;
        let x: Vec<f64> = (0..m).map(|t| closes[t + 1] / closes[t]).collect();
        let rk = bar_close_ordinal_ranks(&x);
        let sk = bar_close_ordinal_ranks(&durations_us[1..]);
        let mf = m as f64;
        let mut sum = 0.0_f64;
        for i in 1..=m {
            for j in 1..=m {
                let mut cnt = 0usize;
                for k in 0..m {
                    if rk[k] < i && sk[k] < j {
                        cnt += 1;
                    }
                }
                let d = (cnt as f64) / mf - ((i as f64) / mf) * ((j as f64) / mf);
                sum += d * d;
            }
        }
        90.0 * (sum / (mf * mf))
    }

    /// Deterministic pseudo-random-ish (but fixed) fixture for the brute-force leg:
    /// a multiplicative close walk + an unrelated integer duration pattern, both
    /// generated from a tiny LCG so the test needs no external data and never drifts.
    fn hoeffding_lcg_fixture(n: usize) -> (Vec<f64>, Vec<f64>) {
        let mut state: u64 = 0x00c0_ffee_d00d_5eed;
        let mut next = || {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1_442_695_040_888_963_407);
            (state >> 33) as f64 / f64::from(1u32 << 31)
        };
        let mut closes = Vec::with_capacity(n);
        let mut durs = Vec::with_capacity(n);
        let mut c = 40_000.0_f64;
        for _ in 0..n {
            c *= 1.0 + (next() - 0.5) * 0.004;
            closes.push((c * 1e8).round() / 1e8); // FixedPoint 8dp, like production
            durs.push((1.0 + (next() * 3_000_000.0)).floor()); // integer µs
        }
        (closes, durs)
    }

    /// Kernel (cumsum grid) == brute force (direct counting) on a 61-close window —
    /// the in-Rust independent leg of the 3-tier cross-check (S2).
    #[test]
    fn hoeffding_kernel_matches_brute_force_independent_leg() {
        let (closes, durs) = hoeffding_lcg_fixture(61);
        let fast = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &durs);
        let brute = hoeffding_phi2_brute(&closes, &durs);
        assert!(
            (fast - brute).abs() <= 1e-12,
            "cumsum-grid {fast} vs direct-count {brute} diverge"
        );
    }

    /// Perfect co-monotone dependence (both substrates strictly increasing) → Φ²
    /// near 1; perfect counter-monotone → also near 1 (Φ² is direction-blind);
    /// the exact grid value approaches 1 from below as m grows.
    #[test]
    fn hoeffding_monotone_dependence_near_one() {
        // closes rising 0.1%/bar → ratios all distinct & increasing? No — constant
        // ratio ties everywhere; use an ACCELERATING walk so ratios strictly rise.
        let n = 200usize;
        let mut closes = Vec::with_capacity(n);
        let mut c = 1000.0_f64;
        for t in 0..n {
            c *= 1.0 + 0.0001 * (t as f64 + 1.0);
            closes.push(c);
        }
        let up: Vec<f64> = (0..n).map(|t| 1_000.0 + 1_000.0 * t as f64).collect();
        let down: Vec<f64> = (0..n).map(|t| 1_000_000.0 - 1_000.0 * t as f64).collect();
        let co = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &up);
        let counter = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &down);
        assert!(co > 0.9, "co-monotone Φ² should approach 1, got {co}");
        assert!(
            counter > 0.9,
            "counter-monotone Φ² should approach 1, got {counter}"
        );
        // Finite-m grid bias: the exact perfect-monotone grid value at m=199 is
        // 1.0000631… (the 90· normalization makes the INTEGRAL version 1; the
        // grid MEAN exceeds it by O(1/m)). Bound with a small explicit ceiling.
        assert!(
            co <= 1.001 && counter <= 1.001,
            "grid Φ² exceeded 1 + O(1/m): {co} / {counter}"
        );
    }

    /// Φ² stays in [0,1] and is deterministic (same input → identical bits).
    #[test]
    fn hoeffding_bounded_and_deterministic() {
        let (closes, durs) = hoeffding_lcg_fixture(200);
        let a = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &durs);
        let b = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &durs);
        assert!((0.0..=1.0).contains(&a), "Φ² out of [0,1]: {a}");
        assert_eq!(a.to_bits(), b.to_bits(), "non-deterministic Φ²");
    }

    /// RANK invariance: power-of-two close scaling is bit-exact (the factor cancels
    /// in each ratio division); any strictly-increasing duration transform (here
    /// affine) leaves the duration ranks — and the Φ² bits — unchanged.
    #[test]
    fn hoeffding_scale_and_monotone_transform_invariant() {
        let (closes, durs) = hoeffding_lcg_fixture(200);
        let base = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &durs);
        let scaled: Vec<f64> = closes.iter().map(|&c| c * 4.0).collect();
        let stretched: Vec<f64> = durs.iter().map(|&d| d * 7.0 + 3.0).collect();
        let v_scaled = compute_bar_hoeffding_phi_squared_midreturn_duration(&scaled, &durs);
        let v_stretched = compute_bar_hoeffding_phi_squared_midreturn_duration(&closes, &stretched);
        assert_eq!(
            base.to_bits(),
            v_scaled.to_bits(),
            "closes×4 changed Φ² bits"
        );
        assert_eq!(
            base.to_bits(),
            v_stretched.to_bits(),
            "monotone duration transform changed Φ² bits"
        );
    }

    /// NaN contract: short window, mismatched substrate lengths, close <= 0,
    /// non-finite close, non-finite duration → NaN. All-tied durations are DEFINED
    /// (stable ordinal ranks fall back to time order), not NaN.
    #[test]
    fn hoeffding_nan_policy() {
        let f = compute_bar_hoeffding_phi_squared_midreturn_duration;
        assert!(f(&[1.0, 2.0, 3.0, 4.0], &[1.0, 2.0, 3.0, 4.0]).is_nan()); // n < 5
        assert!(f(&[1.0, 2.0, 3.0, 4.0, 5.0], &[1.0, 2.0, 3.0, 4.0]).is_nan()); // len mismatch
        assert!(f(&[1.0, 0.0, 3.0, 4.0, 5.0], &[1.0; 5]).is_nan()); // close == 0
        assert!(f(&[1.0, -2.0, 3.0, 4.0, 5.0], &[1.0; 5]).is_nan()); // close < 0
        assert!(f(&[1.0, f64::NAN, 3.0, 4.0, 5.0], &[1.0; 5]).is_nan()); // NaN close
        assert!(f(&[1.0, 2.0, 3.0, 4.0, 5.0], &[1.0, 2.0, f64::NAN, 4.0, 5.0]).is_nan()); // NaN dur
        let tied = f(&[1.0, 2.0, 1.5, 3.0, 2.5, 4.0], &[7.0; 6]);
        assert!(
            tied.is_finite() && (0.0..=1.0).contains(&tied),
            "all-tied durations must be DEFINED"
        );
    }

    /// Never ±Inf on degenerate-but-defined windows.
    #[test]
    fn hoeffding_never_inf_on_degenerate() {
        let cases: Vec<(Vec<f64>, Vec<f64>)> = vec![
            (vec![1.0, 1.0, 1.0, 1.0, 1.000_000_001], vec![1.0; 5]),
            (
                vec![1e-9, 1.0, 1e9, 1.0, 1e-9, 1.0],
                vec![1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
            ),
            (vec![1.0; 5], vec![5.0, 4.0, 3.0, 2.0, 1.0]),
        ];
        for (c, d) in cases {
            assert!(!compute_bar_hoeffding_phi_squared_midreturn_duration(&c, &d).is_infinite());
        }
    }

    /// Stable-tie ordinal ranks: equal values keep input order (the numpy
    /// mergesort double-argsort permutation), distinct values rank by magnitude.
    #[test]
    fn hoeffding_ordinal_ranks_stable_ties() {
        assert_eq!(
            bar_close_ordinal_ranks(&[2.0, 1.0, 2.0, 0.5]),
            vec![2, 1, 3, 0]
        );
        assert_eq!(bar_close_ordinal_ranks(&[7.0; 4]), vec![0, 1, 2, 3]);
    }

    // ── bar_hvg_forward_visibility_horizon_mean (hvghorizon, card 59) ────────────

    /// INDEPENDENT brute-force leg: the probe's naive per-bar forward first-passage
    /// scan (O(n²), no stack) — a structurally different construction, so a stack
    /// bookkeeping bug in the kernel cannot be self-confirming. Returns NaN under
    /// the identical undefined-window contract.
    fn hvghorizon_naive(closes: &[f64]) -> f64 {
        let n = closes.len();
        if n < 2 || !closes.iter().all(|&c| c.is_finite()) {
            return f64::NAN;
        }
        let mut sum: u64 = 0;
        let mut count: u64 = 0;
        for t in 0..n - 1 {
            for k in 1..n - t {
                if closes[t + k] >= closes[t] {
                    sum += k as u64;
                    count += 1;
                    break;
                }
            }
        }
        if count == 0 {
            f64::NAN
        } else {
            (sum as f64) / (count as f64)
        }
    }

    /// Deterministic LCG close walk (same generator family as the hoeffding tests).
    fn hvghorizon_lcg_closes(n: usize) -> Vec<f64> {
        let mut state: u64 = 0x5eed_cafe_f00d_d00d;
        let mut next = || {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1_442_695_040_888_963_407);
            (state >> 33) as f64 / f64::from(1u32 << 31)
        };
        let mut closes = Vec::with_capacity(n);
        let mut c = 40_000.0_f64;
        for _ in 0..n {
            c *= 1.0 + (next() - 0.5) * 0.004;
            closes.push((c * 1e8).round() / 1e8); // FixedPoint 8dp, like production
        }
        closes
    }

    /// Kernel (monotonic stack) == brute force (naive scan), bit-for-bit, on
    /// deterministic pseudo-random windows of several sizes — the in-Rust
    /// independent leg of the 3-tier cross-check (S2).
    #[test]
    fn hvghorizon_kernel_matches_naive_brute_force_leg() {
        for n in [5usize, 23, 61, 200] {
            let closes = hvghorizon_lcg_closes(n);
            let fast = compute_bar_hvg_forward_visibility_horizon_mean(&closes);
            let brute = hvghorizon_naive(&closes);
            assert_eq!(
                fast.to_bits(),
                brute.to_bits(),
                "stack {fast} vs naive {brute} diverge at n={n}"
            );
        }
    }

    /// Closed-form anchors: strictly-increasing and constant windows give mean
    /// EXACTLY 1.0 (every horizon is 1; ties count via `>=`); strictly-decreasing
    /// gives NaN (every bar right-censored); hand-pinned [3,1,2,3] gives 5/3.
    #[test]
    fn hvghorizon_closed_form_anchors() {
        let inc: Vec<f64> = (1..=200).map(|i| i as f64).collect();
        assert_eq!(compute_bar_hvg_forward_visibility_horizon_mean(&inc), 1.0);
        assert_eq!(
            compute_bar_hvg_forward_visibility_horizon_mean(&[7.0; 200]),
            1.0
        );
        let dec: Vec<f64> = (1..=200).rev().map(|i| i as f64).collect();
        assert!(compute_bar_hvg_forward_visibility_horizon_mean(&dec).is_nan());
        // L+(0)=3 (first >= 3.0 at index 3), L+(1)=1, L+(2)=1 -> mean 5/3.
        let v = compute_bar_hvg_forward_visibility_horizon_mean(&[3.0, 1.0, 2.0, 3.0]);
        assert_eq!(v, 5.0 / 3.0, "hand-pinned [3,1,2,3] must give exactly 5/3");
    }

    /// Positive scaling preserves every `>=` comparison → identical integer
    /// horizons → identical value bits (exactly, not to a tolerance).
    #[test]
    fn hvghorizon_positive_scale_invariant_bits() {
        let closes = hvghorizon_lcg_closes(200);
        let base = compute_bar_hvg_forward_visibility_horizon_mean(&closes);
        for &k in &[2.0_f64, 0.5, 4.0, 1000.0] {
            let scaled: Vec<f64> = closes.iter().map(|&c| c * k).collect();
            let got = compute_bar_hvg_forward_visibility_horizon_mean(&scaled);
            assert_eq!(base.to_bits(), got.to_bits(), "scale x{k} changed the bits");
        }
    }

    /// Determinism + bounds: same input → identical bits; defined values lie in
    /// `[1, n-1]`.
    #[test]
    fn hvghorizon_deterministic_and_bounded() {
        let closes = hvghorizon_lcg_closes(200);
        let a = compute_bar_hvg_forward_visibility_horizon_mean(&closes);
        let b = compute_bar_hvg_forward_visibility_horizon_mean(&closes);
        assert_eq!(a.to_bits(), b.to_bits(), "non-deterministic value");
        assert!(
            (1.0..=199.0).contains(&a),
            "mean horizon out of [1, n-1]: {a}"
        );
    }

    /// NaN contract: short window, non-finite close, all-censored window → NaN.
    /// Never ±Inf on degenerate-but-defined windows.
    #[test]
    fn hvghorizon_nan_policy_and_never_inf() {
        let f = compute_bar_hvg_forward_visibility_horizon_mean;
        assert!(f(&[]).is_nan());
        assert!(f(&[42.0]).is_nan()); // n < 2
        assert!(f(&[1.0, f64::NAN, 3.0]).is_nan());
        assert!(f(&[1.0, f64::INFINITY, 3.0]).is_nan());
        assert!(f(&[3.0, 2.0, 1.0]).is_nan()); // strictly decreasing: all censored
        // Degenerate-but-defined: negative and zero closes are comparable values.
        let cases: Vec<Vec<f64>> = vec![
            vec![0.0, 0.0, 0.0, 0.0],
            vec![-3.0, -1.0, -2.0, -0.5],
            vec![1e-12, 1e12, 1e-12, 1e12],
        ];
        for w in cases {
            let v = f(&w);
            assert!(!v.is_infinite(), "±Inf on {w:?}");
        }
    }
    // ── bar_vg_time_directed_clustering_meangap (vgclust, card 60) ──────────────

    /// INDEPENDENT brute-force leg: naive O(n³) pairwise-visibility HVG straight
    /// from the definition + an independent set-based clustering recompute — no
    /// left-scan, no bitset — so a construction bug cannot self-confirm.
    fn vgclust_naive(closes: &[f64]) -> f64 {
        let n = closes.len();
        if n < 5 || !closes.iter().all(|&c| c.is_finite()) {
            return f64::NAN;
        }
        let mut edge = vec![false; n * n];
        for i in 0..n - 1 {
            for j in i + 1..n {
                let lo = closes[i].min(closes[j]);
                if (i + 1..j).all(|k| closes[k] < lo) {
                    edge[i * n + j] = true;
                }
            }
        }
        let has = |a: usize, b: usize| {
            let (i, j) = if a < b { (a, b) } else { (b, a) };
            edge[i * n + j]
        };
        let (mut ps, mut pc, mut fs, mut fc) = (0.0_f64, 0u32, 0.0_f64, 0u32);
        for v in 0..n {
            let past: Vec<usize> = (0..v).filter(|&u| has(u, v)).collect();
            let fut: Vec<usize> = (v + 1..n).filter(|&u| has(v, u)).collect();
            for (set, sum, cnt) in [(&past, &mut ps, &mut pc), (&fut, &mut fs, &mut fc)] {
                let k = set.len();
                if k >= 2 {
                    let mut links = 0u32;
                    for a in 0..k {
                        for b in a + 1..k {
                            if has(set[a], set[b]) {
                                links += 1;
                            }
                        }
                    }
                    *sum += 2.0 * f64::from(links) / ((k * (k - 1)) as f64);
                    *cnt += 1;
                }
            }
        }
        if pc == 0 || fc == 0 {
            return f64::NAN;
        }
        ps / f64::from(pc) - fs / f64::from(fc)
    }

    fn vgclust_lcg_closes(n: usize) -> Vec<f64> {
        let mut state: u64 = 0xdead_beef_cafe_1234;
        let mut next = || {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1_442_695_040_888_963_407);
            (state >> 33) as f64 / f64::from(1u32 << 31)
        };
        let mut closes = Vec::with_capacity(n);
        let mut c = 40_000.0_f64;
        for _ in 0..n {
            c *= 1.0 + (next() - 0.5) * 0.004;
            closes.push((c * 1e8).round() / 1e8);
        }
        closes
    }

    /// Kernel (left-scan + bitset) == brute force (naive pairwise + sets),
    /// bit-for-bit, on pseudo-random windows AND the tie battery that killed the
    /// monotonic-stack variant in the oracle generator (incident P1).
    #[test]
    fn vgclust_kernel_matches_naive_brute_force_leg() {
        for n in [5usize, 23, 61, 200] {
            let closes = vgclust_lcg_closes(n);
            let fast = compute_bar_vg_time_directed_clustering_meangap(&closes);
            let brute = vgclust_naive(&closes);
            assert_eq!(
                fast.to_bits(),
                brute.to_bits(),
                "left-scan {fast} vs naive {brute} diverge at n={n}"
            );
        }
        // Tie battery (P1): equal-height runs must BLOCK visibility identically.
        let ties: Vec<Vec<f64>> = vec![
            vec![2.0, 1.0, 1.0, 2.0, 1.5],
            vec![3.0, 3.0, 1.0, 3.0, 2.0],
            vec![5.0, 5.0, 5.0, 5.0, 5.0],
            vec![1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.0],
            vec![2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
        ];
        for w in ties {
            let fast = compute_bar_vg_time_directed_clustering_meangap(&w);
            let brute = vgclust_naive(&w);
            assert!(
                fast.to_bits() == brute.to_bits() || (fast.is_nan() && brute.is_nan()),
                "tie case {w:?}: left-scan {fast} vs naive {brute}"
            );
        }
    }

    /// Time-reversal flips the sign — EXACT in real arithmetic, ~1e-12 in IEEE
    /// (incident P2: reversal permutes the mean summation order). NOT to_bits.
    #[test]
    fn vgclust_time_reversal_flips_sign_within_1e12() {
        for n in [23usize, 61, 200] {
            let closes = vgclust_lcg_closes(n);
            let v = compute_bar_vg_time_directed_clustering_meangap(&closes);
            let rev: Vec<f64> = closes.iter().rev().copied().collect();
            let vr = compute_bar_vg_time_directed_clustering_meangap(&rev);
            assert!(
                (vr + v).abs() <= 1e-12,
                "reversal flip violated at n={n}: {v} vs {vr}"
            );
        }
    }

    /// Anchors: the symmetric sawtooth [2,1,2,1,2] gives EXACTLY 0.0 (hand-worked:
    /// past_c = fut_c = [1.0, 1.0]); strictly monotone windows give NaN (the HVG
    /// is a bare path — no node has 2 past or 2 future neighbours).
    #[test]
    fn vgclust_closed_form_anchors() {
        let v = compute_bar_vg_time_directed_clustering_meangap(&[2.0, 1.0, 2.0, 1.0, 2.0]);
        assert_eq!(v, 0.0, "symmetric sawtooth must give exactly 0.0, got {v}");
        let inc: Vec<f64> = (1..=200).map(|i| i as f64).collect();
        assert!(compute_bar_vg_time_directed_clustering_meangap(&inc).is_nan());
        let dec: Vec<f64> = (1..=200).rev().map(|i| i as f64).collect();
        assert!(compute_bar_vg_time_directed_clustering_meangap(&dec).is_nan());
    }

    /// Positive pow-2 rescale preserves every visibility comparison → identical
    /// edge set → identical value bits. Determinism + [-1,1] bounds.
    #[test]
    fn vgclust_scale_invariant_deterministic_bounded() {
        let closes = vgclust_lcg_closes(200);
        let base = compute_bar_vg_time_directed_clustering_meangap(&closes);
        for &k in &[2.0_f64, 0.5, 4.0, 0.25] {
            let scaled: Vec<f64> = closes.iter().map(|&c| c * k).collect();
            let got = compute_bar_vg_time_directed_clustering_meangap(&scaled);
            assert_eq!(
                base.to_bits(),
                got.to_bits(),
                "pow2 rescale x{k} changed bits"
            );
        }
        let again = compute_bar_vg_time_directed_clustering_meangap(&closes);
        assert_eq!(base.to_bits(), again.to_bits(), "non-deterministic");
        assert!((-1.0..=1.0).contains(&base), "out of [-1,1]: {base}");
    }

    /// NaN contract + never ±Inf on degenerate windows.
    #[test]
    fn vgclust_nan_policy_and_never_inf() {
        let f = compute_bar_vg_time_directed_clustering_meangap;
        assert!(f(&[]).is_nan());
        assert!(f(&[1.0, 2.0, 3.0, 4.0]).is_nan()); // n < 5
        assert!(f(&[1.0, f64::NAN, 3.0, 4.0, 5.0]).is_nan());
        assert!(f(&[1.0, f64::INFINITY, 3.0, 4.0, 5.0]).is_nan());
        let cases: Vec<Vec<f64>> = vec![
            vec![7.0; 200],
            vec![0.0, 1.0, 0.0, 1.0, 0.0],
            vec![-3.0, -1.0, -2.0, -0.5, -4.0],
        ];
        for w in cases {
            assert!(!f(&w).is_infinite(), "±Inf on {w:?}");
        }
    }
}