opendeviationbar-core 13.76.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;

/// 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 }
}

#[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]
                );
            }
        }
    }
}