viterbi 0.2.0

// Author: BolivarTech
// Version: 1.0.0
// Date: 2026-07-02
//! Shared helper for the integration test suites (Task 10).
//!
//! Each `tests/*.rs` is compiled as its own crate; this module is pulled in
//! with `mod common;` so the `codec()` fixture and the `run_ber` BER harness
//! are defined once (no shared mutable state — every caller gets fresh
//! encoder/decoder instances). `run_ber` is NOT feature-gated, so the anchor
//! gate in `ber_anchor.rs` always runs under `cargo nextest run`.

// Each `tests/*.rs` including this module uses only a subset of its items, so
// the unused ones are dead code *from that crate's view* — expected for a
// shared test fixture, not a real smell.
#![allow(dead_code)]

use rand::rngs::StdRng;
use rand::{Rng, SeedableRng};
use viterbi::{
    decoder::{CcsdsViterbiDecoder, ViterbiDecoder},
    encoder::ViterbiEncoder,
    metric::{HardHamming, SoftLlr},
    params::CodeParams,
    puncture::{DePuncturer, Puncturer},
    PunctureMatrix,
};

/// Build a fresh CCSDS K=7 R=1/2 encoder + hard-decision decoder pair.
///
/// The decoder is sized for `65_536` info bits, comfortably above every block
/// used by the integration tests.
pub fn codec() -> (ViterbiEncoder, CcsdsViterbiDecoder) {
    (
        ViterbiEncoder::new(CodeParams::ccsds_r1_2()).unwrap(),
        CcsdsViterbiDecoder::new(CodeParams::ccsds_r1_2(), 65_536).unwrap(),
    )
}

/// LLR quantization scale for [`SoftLlr::llr_from_bpsk`]. Calibrated so a typical
/// noisy BPSK sample (mean ±1, σ ≈ 0.79 at the 2 dB anchor) maps across a wide
/// span of the `i8` LLR range `[-127, 127]` without mostly saturating — that
/// preserves the reliability magnitude which is exactly the information the soft
/// path exploits and the hard slice discards. Too small ⇒ collapses to a few
/// levels (loses the gain); too large ⇒ saturates to ±127 (degenerates to hard).
const LLR_SCALE: f32 = 24.0;

/// Deterministic BER of one decode path over a seeded AWGN + BPSK channel.
///
/// Same-input determinism is the whole point: `seed` drives BOTH the random info
/// bits AND the Gaussian noise, so a hard call and a soft call with the *same*
/// `seed` decode the **identical** noisy samples. Soft's only advantage is the
/// LLR magnitude it keeps versus the hard 1-bit slice — so `soft <= hard` is
/// guaranteed by physics; if it ever fails, the bug is in this harness (σ
/// formula, LLR sign/scale, or the hard slice), never in the decoder.
///
/// # Physics
/// - BPSK: coded bit `0 → +1.0`, `1 → -1.0` (matches the `λ>0 ⇒ bit 0` LLR
///   convention, so a positive sample yields a positive LLR ⇒ bit 0).
/// - AWGN: `N(0, σ²)` with `σ = sqrt(1 / (2·R·Eb/N0))`, `R = 1/2`,
///   `Eb/N0 = 10^(eb_n0_db/10)`.
///
/// # Panics
/// Panics if `n_bits` is not a multiple of 8 (the harness works byte-aligned so
/// the recovered bytes align with the source bytes for a clean bit comparison).
///
/// # Returns
/// The measured bit-error rate `errors / n_bits` as `f64`.
pub fn run_ber(eb_n0_db: f64, n_bits: usize, seed: u64, soft: bool) -> f64 {
    run_ber_params(eb_n0_db, n_bits, seed, soft, CodeParams::ccsds_r1_2())
}

/// Same deterministic AWGN+BPSK BER sweep as [`run_ber`], but over an arbitrary rate-`1/n`
/// `params` (rate-1/2 or rate-1/3). The code rate `R = 1/n` is read from `params.n()` and drives
/// the σ formula, so the Eb/N0 → σ conversion holds *energy per info bit* constant across rates —
/// the correct basis for comparing a rate-1/2 vs a rate-1/3 code at the "same Eb/N0". Both the hard
/// and soft front-ends reuse the generic `ViterbiDecoder<64, _>::decode` path (the runtime-`n`
/// decode Task 6 adds), so this exercises the rate-1/3 decoder end-to-end.
pub fn run_ber_params(
    eb_n0_db: f64,
    n_bits: usize,
    seed: u64,
    soft: bool,
    params: CodeParams,
) -> f64 {
    assert!(
        n_bits.is_multiple_of(8),
        "run_ber expects a byte-aligned bit budget"
    );
    let code_rate = 1.0 / params.n() as f64; // R = 1/n (0.5 for rate-1/2, 1/3 for rate-1/3)
    let mut rng = StdRng::seed_from_u64(seed);

    // 1. Random info bits (byte-aligned) — first draw from the seeded stream.
    let mut data = vec![0u8; n_bits / 8];
    rng.fill(data.as_mut_slice());

    // 2. Encode: zero-tail → (n_bits + M)·n coded bits (n = params.n()).
    let enc = ViterbiEncoder::new(params.clone()).unwrap();
    let coded = enc.encode_bits(&data, n_bits).unwrap();
    let ncoded = coded.nbits;

    // 3. BPSK + AWGN. The RNG stream up to here is identical regardless of
    //    `soft`, so hard and soft observe the SAME noisy samples.
    let eb_n0_lin = 10f64.powf(eb_n0_db / 10.0);
    let sigma = (1.0 / (2.0 * code_rate * eb_n0_lin)).sqrt();
    let noise = gaussian(&mut rng, ncoded);
    let mut noisy = Vec::with_capacity(ncoded);
    for (i, &nz) in noise.iter().enumerate() {
        let symbol = if read_bit(&coded.bytes, i) == 0 {
            1.0
        } else {
            -1.0
        };
        noisy.push(symbol + sigma * nz);
    }

    // 4. Decode via the requested front-end. Both reuse the generic ACS core.
    let decoded = if soft {
        let llrs: Vec<i8> = noisy
            .iter()
            .map(|&y| SoftLlr::llr_from_bpsk(y as f32, LLR_SCALE))
            .collect();
        let mut dec = ViterbiDecoder::<64, SoftLlr>::new(params, n_bits).unwrap();
        dec.decode(&llrs, n_bits).unwrap()
    } else {
        let hard: Vec<u8> = noisy.iter().map(|&y| u8::from(y < 0.0)).collect();
        let mut dec = ViterbiDecoder::<64, HardHamming>::new(params, n_bits).unwrap();
        dec.decode(&hard, n_bits).unwrap()
    };

    // 5. Bit errors over the info bits (byte-aligned ⇒ direct byte XOR popcount).
    let errors: u32 = data
        .iter()
        .zip(decoded.bytes.iter())
        .map(|(a, b)| (a ^ b).count_ones())
        .sum();
    f64::from(errors) / n_bits as f64
}

/// Deterministic BER of one decode path over a seeded AWGN + BPSK channel with
/// **puncturing** applied at TX and de-puncturing at RX.
///
/// Same channel/seed discipline as [`run_ber_params`], but the coded stream is
/// punctured to a higher rate before the channel and de-punctured (erasures
/// reinserted via `M::erasure()`) before the decoder — so this exercises the whole
/// puncture path end-to-end. The Eb/N0 → σ conversion uses the **punctured** code
/// rate `R = period / kept_per_period` (info bits per transmitted bit), which holds
/// energy-per-info-bit constant across puncture rates — the correct basis for a
/// fair BER comparison between rates at the "same Eb/N0". Because a higher puncture
/// rate removes more redundancy, its BER is expected to be **worse** (higher) at a
/// fixed Eb/N0; that physically-correct ordering is what the anchor gate asserts.
///
/// # Panics
/// Panics if `n_bits` is not a multiple of 8 (byte-aligned, as in [`run_ber`]).
pub fn run_ber_punctured(
    eb_n0_db: f64,
    n_bits: usize,
    seed: u64,
    soft: bool,
    params: CodeParams,
    matrix: PunctureMatrix,
) -> f64 {
    assert!(
        n_bits.is_multiple_of(8),
        "run_ber_punctured expects a byte-aligned bit budget"
    );
    // Punctured code rate R = info/coded = period / kept_per_period (e.g. 2/3 for R2_3).
    let code_rate = matrix.period() as f64 / matrix.kept_per_period() as f64;
    let mut rng = StdRng::seed_from_u64(seed);

    // 1. Random info bits (byte-aligned) — first draw from the seeded stream.
    let mut data = vec![0u8; n_bits / 8];
    rng.fill(data.as_mut_slice());

    // 2. Encode to the full rate-1/n stream, then puncture to the higher rate.
    let enc = ViterbiEncoder::new(params.clone()).unwrap();
    let coded = enc.encode_bits(&data, n_bits).unwrap();
    let punc = Puncturer::new(matrix.clone(), &params).unwrap();
    let punctured = punc.puncture(&coded, n_bits).unwrap();
    let ncoded = punctured.nbits; // on-air (transmitted) bit count

    // 3. BPSK + AWGN over the PUNCTURED (on-air) stream. The RNG stream up to here is
    //    identical regardless of `soft`, so hard and soft observe the SAME samples.
    let eb_n0_lin = 10f64.powf(eb_n0_db / 10.0);
    let sigma = (1.0 / (2.0 * code_rate * eb_n0_lin)).sqrt();
    let noise = gaussian(&mut rng, ncoded);
    let mut noisy = Vec::with_capacity(ncoded);
    for (i, &nz) in noise.iter().enumerate() {
        let symbol = if read_bit(&punctured.bytes, i) == 0 {
            1.0
        } else {
            -1.0
        };
        noisy.push(symbol + sigma * nz);
    }

    // 4. De-puncture (reinsert `M::erasure()` at punctured positions) then decode via
    //    the requested front-end — both reuse the generic ACS core unchanged.
    let decoded = if soft {
        let llrs: Vec<i8> = noisy
            .iter()
            .map(|&y| SoftLlr::llr_from_bpsk(y as f32, LLR_SCALE))
            .collect();
        let dp = DePuncturer::<SoftLlr>::new(matrix, &params).unwrap();
        let full = dp.depuncture(&llrs, n_bits).unwrap();
        let mut dec = ViterbiDecoder::<64, SoftLlr>::new(params, n_bits).unwrap();
        dec.decode(&full, n_bits).unwrap()
    } else {
        let hard: Vec<u8> = noisy.iter().map(|&y| u8::from(y < 0.0)).collect();
        let dp = DePuncturer::<HardHamming>::new(matrix, &params).unwrap();
        let full = dp.depuncture(&hard, n_bits).unwrap();
        let mut dec = ViterbiDecoder::<64, HardHamming>::new(params, n_bits).unwrap();
        dec.decode(&full, n_bits).unwrap()
    };

    // 5. Bit errors over the info bits (byte-aligned ⇒ direct byte XOR popcount).
    let errors: u32 = data
        .iter()
        .zip(decoded.bytes.iter())
        .map(|(a, b)| (a ^ b).count_ones())
        .sum();
    f64::from(errors) / n_bits as f64
}

/// Read the MSB-first bit at index `i` (bit 7 of byte 0 is index 0) — mirrors the
/// crate's internal packing convention so BPSK mapping lines up with the encoder.
fn read_bit(bytes: &[u8], i: usize) -> u8 {
    (bytes[i / 8] >> (7 - (i % 8))) & 1
}

/// `n` deterministic standard-normal `N(0, 1)` samples via Box–Muller drawn from
/// the seeded RNG (avoids pulling in `rand_distr`). Fixed seed ⇒ reproducible.
fn gaussian(rng: &mut StdRng, n: usize) -> Vec<f64> {
    let mut out = Vec::with_capacity(n);
    while out.len() < n {
        // u1 in (0, 1] so `ln` is finite (rng.gen::<f64>() is [0, 1)).
        let u1 = 1.0 - rng.gen::<f64>();
        let u2 = rng.gen::<f64>();
        let r = (-2.0 * u1.ln()).sqrt();
        let theta = 2.0 * std::f64::consts::PI * u2;
        out.push(r * theta.cos());
        if out.len() < n {
            out.push(r * theta.sin());
        }
    }
    out
}