gwseq-io 0.2.1

//! fqzcomp — CRAM 3.1's quality-score codec.
//!
//! Reference: `docs/cram_codecs_v3.1.md` §6.
//!
//! Every other codec here is general: hand it bytes and it compresses them.
//! This one only compresses quality strings, and it does it by *predicting* the
//! next quality from what it has already seen — the previous few qualities, how
//! far along the read it is, how many times the quality has changed so far, and
//! an arbitrary per-record selector. Those are combined into a 16-bit context,
//! and each of the `2^16` contexts has its own adaptive model. §4's range coder
//! and models do the actual coding, multiplexed into one stream.
//!
//! That makes the parameter block, not the entropy coding, the substance of
//! this module: where in the 16 bits each sub-context sits (`qloc`, `ploc`,
//! `dloc`, `sloc`), how wide the quality history is (`qbits`, `qshift`), and
//! four lookup tables that squeeze a wide input range into a narrow slice of
//! context. An encoder chooses all of it per file, so a decoder that assumes
//! any of it is wrong on the next file it meets.
//!
//! Since it stores read lengths and reverse flags of its own — CRAM keeps data
//! series strictly separate, so the codec cannot look at `RL` — a quality block
//! decodes with no reference to the rest of the slice.
//!
//! # Three corrections to the pseudocode
//!
//! §6.2's `FQZDecodeSingleParam` opens with a list of aliases that contradict
//! its own body and §6.2's flag table: it has `have_ptab <- flags AND 16`,
//! `do_rev <- flags AND 16` and `have_qmap <- flags AND 1`. The table and the
//! function body agree with each other — 16 is `have_qmap`, 32 is `have_ptab`,
//! 64 is `have_dtab`, 128 is `have_qtab` — and that is what [`pflag`] uses.
//! `do_rev` is not a per-parameter flag at all; it is global, bit 4 of
//! `gflags`.
//!
//! §6.2's `ReverseQualities` advances `i` and `rec` only inside
//! `if rev_rec != 0`, so the first record that is *not* reversed loops forever.
//! Both advance every record; only the swapping is conditional.
//!
//! §6.1 says the position context "start[s] at record length (minus 1) and
//! decrement[s]", but §6.2's main loop decrements `pos` after calling
//! `FQZUpdateContext` rather than before, which would start it at the record
//! length. The prose is right, and it is the reading that agrees with what
//! `samtools` writes.

use crate::error::{Error, Result};

use super::arith::{Model, RangeCoder};
use super::ByteReader;

/// The only version this format has had.
const VERSION: u8 = 5;

/// Global flags, §6.2.
mod gflag {
    pub const MULTI_PARAM: u8 = 1;
    pub const HAVE_STAB: u8 = 2;
    pub const DO_REV: u8 = 4;
}

/// Per-parameter-block flags, §6.2. Bit 1 is reserved.
pub mod pflag {
    pub const DO_DEDUP: u8 = 2;
    /// §6.2 calls this `do_len`, "model_len will be used for every record".
    /// It means the opposite. `htslib` sets it when every record in the block
    /// has the *same* length (`pm->fixed_len = (i == s->num_records)`) and
    /// then reads a length only for the first record
    /// (`if (!pm->fixed_len || state->first_len)`). Taken at the
    /// specification's word, every record after the first steals four bytes
    /// that belong to the next quality string.
    pub const FIXED_LEN: u8 = 4;
    pub const DO_SEL: u8 = 8;
    pub const HAVE_QMAP: u8 = 16;
    pub const HAVE_PTAB: u8 = 32;
    pub const HAVE_DTAB: u8 = 64;
    pub const HAVE_QTAB: u8 = 128;
}

/// One parameter set: how to build a context, and how to read a quality back.
struct Param {
    context: u16,
    flags: u8,
    max_sym: u8,
    qbits: u32,
    qshift: u32,
    qloc: u32,
    sloc: u32,
    ploc: u32,
    dloc: u32,
    /// Quality values, un-binned. Sized 256 rather than `max_sym` because the
    /// quality model has `max_sym + 1` symbols and so can emit `max_sym`.
    qmap: [u8; 256],
    qtab: [u8; 256],
    ptab: [u8; 1024],
    dtab: [u8; 256],
}

impl Param {
    fn read(reader: &mut ByteReader<'_>) -> Result<Self> {
        let context = u16::from(reader.u8()?) | (u16::from(reader.u8()?) << 8);
        let flags = reader.u8()?;
        let max_sym = reader.u8()?;
        let x = reader.u8()?;
        let (qbits, qshift) = (u32::from(x / 16), u32::from(x % 16));
        let x = reader.u8()?;
        let (qloc, sloc) = (u32::from(x / 16), u32::from(x % 16));
        let x = reader.u8()?;
        let (ploc, dloc) = (u32::from(x / 16), u32::from(x % 16));

        let mut param = Self {
            context,
            flags,
            max_sym,
            qbits,
            qshift,
            qloc,
            sloc,
            ploc,
            dloc,
            qmap: [0; 256],
            qtab: [0; 256],
            ptab: [0; 1024],
            dtab: [0; 256],
        };

        if flags & pflag::HAVE_QMAP != 0 {
            for i in 0..max_sym as usize {
                param.qmap[i] = reader.u8()?;
            }
        }
        // The `qbits != 0` half of this condition is `htslib`'s and not §6.2's
        // (`if (pm->qbits) { if (pm->use_qtab) ... }`). A quality context zero
        // bits wide has no use for the table, so the table is not in the
        // stream even when the flag says it is — and reading it there would
        // take its bytes out of whatever follows.
        if flags & pflag::HAVE_QTAB != 0 && qbits != 0 {
            read_array(reader, 256, &mut param.qtab)?;
        } else {
            // §6.2's default is the identity, which is what makes `qtab`
            // safe to index unconditionally below.
            for (i, slot) in param.qtab.iter_mut().enumerate() {
                *slot = i as u8;
            }
        }
        if flags & pflag::HAVE_PTAB != 0 {
            read_array(reader, 1024, &mut param.ptab)?;
        }
        if flags & pflag::HAVE_DTAB != 0 {
            read_array(reader, 256, &mut param.dtab)?;
        }
        Ok(param)
    }

    /// §6.1's `FQZUpdateContext`: the context for the *next* quality.
    ///
    /// Every field it shifts by is four bits wide, so no shift here can reach
    /// the width of a `u32` and the sum stays under `2^31`.
    fn update(&self, q: u8, pos: usize, state: &mut RecordState) -> u16 {
        // Not `self.context`. §6.1 opens `FQZUpdateContext` with
        // `ctx <- params.context`, but `htslib` opens it with
        // `unsigned int last = 0; // pm->context` — the starting context
        // applies to a record's *first* quality only, and is not re-added to
        // every one after it. Both agree while `context` is zero, which is
        // what `samtools` writes, and disagree on the first file that is not.
        let mut ctx = 0u32;
        state.qctx = (state.qctx << self.qshift) + u32::from(self.qtab[q as usize]);
        let qmask = (1u32 << self.qbits) - 1;
        ctx += (state.qctx & qmask) << self.qloc;
        // The tables are added unconditionally: a table that was not stored is
        // all zeroes, so the flags need not be re-tested here.
        ctx += u32::from(self.ptab[pos.min(1023)]) << self.ploc;
        ctx += u32::from(self.dtab[state.delta.min(255) as usize]) << self.dloc;
        ctx += u32::from(state.sel) << self.sloc;
        // Delta counts *changes*, so a run of one value does not raise it.
        state.delta += u32::from(state.prevq != q);
        state.prevq = q;
        (ctx & 0xffff) as u16
    }
}

/// The per-record state the context is built from, all of it reset at each new
/// record.
#[derive(Default)]
struct RecordState {
    qctx: u32,
    delta: u32,
    prevq: u8,
    sel: u8,
}

/// §6.2's `ReadArray`: a monotonically increasing map, run-length encoded
/// twice over.
///
/// The first level turns the map into run lengths — how many inputs share each
/// output value, zeros included — with runs over 255 split into several bytes.
/// The second level collapses repeats of *those*: whenever a length equals the
/// one before it, a count of extra copies follows. §6.2's worked example is
/// `A = 0,1,3,4,5,6,7,7,7,7` becoming `R = 1,1,0,1,1,1,1,4` becoming
/// `R2 = 1,1,+0,0,1,1,+2,4`.
fn read_array(reader: &mut ByteReader<'_>, n: usize, out: &mut [u8]) -> Result<()> {
    debug_assert!(out.len() >= n);
    // One entry per distinct output value, plus one per 255-byte split. The
    // ceiling is generous — the real bound is 256 + n/255 — and it is here
    // because a stream of zero runs would otherwise grow this without ever
    // advancing `total`.
    let limit = n + 512;
    let mut runs: Vec<u8> = Vec::new();
    let mut total = 0usize;
    let mut last: i32 = -1;
    while total < n {
        let run = reader.u8()?;
        runs.push(run);
        total += run as usize;
        if i32::from(run) == last {
            let copies = reader.u8()? as usize;
            for _ in 0..copies {
                runs.push(run);
            }
            total += run as usize * copies;
        }
        last = i32::from(run);
        if runs.len() > limit {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                "an fqzcomp lookup table whose run lengths never reach its size",
            ));
        }
    }

    let mut value = 0usize;
    let mut written = 0usize;
    let mut j = 0usize;
    while written < n {
        // A run of exactly 255 means "and the next byte continues it".
        let mut run_len = 0usize;
        loop {
            let part = *runs.get(j).ok_or_else(|| {
                Error::corrupt(
                    reader.path(),
                    reader.offset(),
                    "an fqzcomp lookup table that ends inside a run",
                )
            })?;
            j += 1;
            run_len += part as usize;
            if part != 255 {
                break;
            }
        }
        if value > 255 {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                "an fqzcomp lookup table mapping past the 255 values a byte holds",
            ));
        }
        for _ in 0..run_len {
            if written >= n {
                break;
            }
            out[written] = value as u8;
            written += 1;
        }
        value += 1;
    }
    Ok(())
}

/// Everything §6.2's parameter block holds.
struct Params {
    params: Vec<Param>,
    /// Selector to parameter index. The identity unless `have_stab`.
    stab: [u8; 256],
    max_sel: usize,
    max_sym: u8,
    do_rev: bool,
}

impl Params {
    /// §6.2's `FQZDecodeParams`.
    fn read(reader: &mut ByteReader<'_>) -> Result<Self> {
        let version = reader.u8()?;
        if version != VERSION {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                format!("an fqzcomp block of version {version}, where only {VERSION} is defined"),
            ));
        }
        let gflags = reader.u8()?;
        let (n_param, mut max_sel) = if gflags & gflag::MULTI_PARAM != 0 {
            let n = reader.u8()? as usize;
            (n, n)
        } else {
            (1, 0)
        };
        if n_param == 0 {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                "an fqzcomp block with no parameter sets",
            ));
        }

        let mut stab = [0u8; 256];
        for (i, slot) in stab.iter_mut().enumerate() {
            *slot = i as u8;
        }
        if gflags & gflag::HAVE_STAB != 0 {
            max_sel = reader.u8()? as usize;
            read_array(reader, 256, &mut stab)?;
        }

        let mut params = Vec::with_capacity(n_param.min(256));
        let mut max_sym = 0u8;
        for _ in 0..n_param {
            let param = Param::read(reader)?;
            max_sym = max_sym.max(param.max_sym);
            params.push(param);
        }
        Ok(Self {
            params,
            stab,
            max_sel,
            max_sym,
            do_rev: gflags & gflag::DO_REV != 0,
        })
    }
}

/// The quality models: one per 16-bit context, built when first reached.
///
/// There are 65536 of them and a file may use a handful. Built eagerly at
/// `max_sym` of 255 they would be tens of megabytes per block, on every worker
/// thread at once; boxed and lazy they cost half a megabyte of pointers plus
/// what the data actually touches.
struct QualityModels {
    models: Vec<Option<Box<Model>>>,
    n_symbols: usize,
}

impl QualityModels {
    fn new(max_sym: u8) -> Self {
        Self {
            models: (0..1 << 16).map(|_| None).collect(),
            n_symbols: max_sym as usize + 1,
        }
    }

    fn get(&mut self, context: u16) -> &mut Model {
        let n_symbols = self.n_symbols;
        self.models[context as usize].get_or_insert_with(|| Box::new(Model::new(n_symbols)))
    }
}

/// §6.2's `DecodeLength`: four bytes, each with its own model.
fn decode_length(
    models: &mut [Model],
    rc: &mut RangeCoder,
    reader: &mut ByteReader<'_>,
) -> Result<usize> {
    let mut len = 0usize;
    for (i, model) in models.iter_mut().enumerate() {
        len |= (model.decode(rc, reader)? as usize) << (i * 8);
    }
    Ok(len)
}

/// Decode a complete fqzcomp block.
pub fn decode(data: &[u8], path: &str, offset: u64) -> Result<Vec<u8>> {
    let mut reader = ByteReader::new(data, path, offset);
    let buf_len = reader.length()?;
    let params = Params::read(&mut reader)?;

    let mut rc = RangeCoder::new(&mut reader)?;
    let mut len_models: Vec<Model> = (0..4).map(|_| Model::new(256)).collect();
    let mut qual_models = QualityModels::new(params.max_sym);
    let mut dup_model = Model::new(2);
    let mut rev_model = Model::new(2);
    let mut sel_model = Model::new(params.max_sel + 1);

    let mut out = vec![0u8; buf_len];
    // Only kept when `do_rev` is set, and then one entry per record.
    let mut reversed: Vec<(usize, bool)> = Vec::new();

    let mut state = RecordState::default();
    let mut first_len = true;
    let mut last_len = 0usize;
    let mut i = 0usize;
    let mut pos = 0usize;
    let mut which = 0usize;
    let mut ctx = 0u16;

    while i < buf_len {
        if pos == 0 {
            // ---- §6.2's FQZNewRecord ------------------------------------
            state = RecordState::default();
            // Gated on the parameter set that is *current* when the record
            // opens — the previous record's, or the first set to begin with —
            // because the selector is what chooses the next one. `htslib`
            // reads `pm->do_sel` before reassigning `pm`.
            if params.params[which].flags & pflag::DO_SEL != 0 {
                state.sel = sel_model.decode(&mut rc, &mut reader)?;
            }
            which = params.stab[state.sel as usize] as usize;
            if which >= params.params.len() {
                return Err(Error::corrupt(
                    path,
                    reader.offset(),
                    format!(
                        "an fqzcomp selector choosing parameter set {which} of {}",
                        params.params.len()
                    ),
                ));
            }

            // A length is stored when the block's records vary in length, and
            // for the first record either way. Both flags are the decoder's,
            // not the parameter set's — `htslib` keeps them in `fqz_state`,
            // so switching parameter sets mid-block does not reset them.
            let rec_len = if params.params[which].flags & pflag::FIXED_LEN == 0 || first_len {
                first_len = false;
                last_len = decode_length(&mut len_models, &mut rc, &mut reader)?;
                last_len
            } else {
                last_len
            };
            if rec_len == 0 || rec_len > buf_len - i {
                return Err(Error::corrupt(
                    path,
                    reader.offset(),
                    format!(
                        "an fqzcomp record of {rec_len} qualities with {} left in the block",
                        buf_len - i
                    ),
                ));
            }
            pos = rec_len;

            if params.do_rev {
                let flag = rev_model.decode(&mut rc, &mut reader)? != 0;
                reversed.push((rec_len, flag));
            }

            let is_dup = params.params[which].flags & pflag::DO_DEDUP != 0
                && dup_model.decode(&mut rc, &mut reader)? > 0;
            if is_dup {
                // The whole record repeats the one before it, which must
                // therefore exist and be at least as long.
                if rec_len > i || i + rec_len > buf_len {
                    return Err(Error::corrupt(
                        path,
                        reader.offset(),
                        "an fqzcomp duplicate record with nothing before it to copy",
                    ));
                }
                out.copy_within(i - rec_len..i, i);
                i += rec_len;
                pos = 0;
                continue;
            }
            ctx = params.params[which].context;
        }

        let param = &params.params[which];
        let q = qual_models.get(ctx).decode(&mut rc, &mut reader)?;
        out[i] = if param.flags & pflag::HAVE_QMAP != 0 {
            param.qmap[q as usize]
        } else {
            q
        };
        // §6.1: the position sub-context is the number of qualities *left*, so
        // it is decremented before the context for the next one is built.
        pos -= 1;
        ctx = param.update(q, pos, &mut state);
        i += 1;
    }

    if params.do_rev {
        reverse_qualities(&mut out, &reversed);
    }
    Ok(out)
}

/// §6.2's `ReverseQualities`, with its loop fixed: every record advances the
/// cursor, and only the reversal is conditional.
fn reverse_qualities(out: &mut [u8], records: &[(usize, bool)]) {
    let mut i = 0usize;
    for &(len, reverse) in records {
        if i + len > out.len() {
            return;
        }
        if reverse {
            out[i..i + len].reverse();
        }
        i += len;
    }
}

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

    /// §6.2's own worked example for `ReadArray`, which is the one part of
    /// this format that can be tested against the specification's numbers
    /// rather than against an encoder written here.
    #[test]
    fn the_specifications_read_array_example_decodes_as_documented() {
        // R2 = 1, 1, +0, 0, 1, 1, +2, 4
        let bytes = [1u8, 1, 0, 0, 1, 1, 2, 4];
        let mut reader = ByteReader::new(&bytes, "spec", 0);
        let mut out = [0u8; 10];
        read_array(&mut reader, 10, &mut out).expect("the spec's own array");
        assert_eq!(out, [0, 1, 3, 4, 5, 6, 7, 7, 7, 7]);
        assert!(reader.is_empty());
    }

    /// §6.2: "a run of length 600 becomes 255 255 90", and the second level
    /// then makes that "255 255 0 90".
    #[test]
    fn a_run_longer_than_255_is_split_and_rejoined() {
        let bytes = [255u8, 255, 0, 90];
        let mut reader = ByteReader::new(&bytes, "spec", 0);
        let mut out = vec![0u8; 600];
        read_array(&mut reader, 600, &mut out).expect("split run");
        assert!(out.iter().all(|v| *v == 0), "every input maps to value 0");
    }

    #[test]
    fn a_table_whose_runs_never_reach_its_size_is_refused() {
        // Nothing but zero runs: `total` never advances.
        let bytes = [0u8; 4096];
        let mut reader = ByteReader::new(&bytes, "test", 0);
        let mut out = [0u8; 256];
        let error = read_array(&mut reader, 256, &mut out).expect_err("never terminates");
        assert!(
            error.to_string().contains("never reach its size"),
            "{error}"
        );
    }

    // ---- a minimal encoder, for the round trips ------------------------
    //
    // It writes one parameter set with no lookup tables and no quality map,
    // which is enough to exercise the range coder, the length models, the
    // quality models and the context update. The tables and the multi-set
    // paths are exercised by the files `samtools` writes; see
    // `cram/local_checks.rs`.

    use super::super::arith::testing::{ModelEncoder, RangeEncoder};

    struct Encoder {
        qbits: u32,
        qshift: u32,
        qloc: u32,
        ploc: u32,
        dloc: u32,
        flags: u8,
        ptab: [u8; 1024],
        dtab: [u8; 256],
    }

    impl Encoder {
        fn simple() -> Self {
            Self {
                qbits: 8,
                qshift: 2,
                qloc: 0,
                ploc: 8,
                dloc: 14,
                // No flags at all: variable-length records, so every one
                // carries its own length.
                flags: 0,
                ptab: [0; 1024],
                dtab: [0; 256],
            }
        }

        fn encode(&self, records: &[Vec<u8>]) -> Vec<u8> {
            let max_sym = records
                .iter()
                .flatten()
                .map(|q| u16::from(*q) + 1)
                .max()
                .unwrap_or(1) as u8;

            let mut header = Vec::new();
            let total: usize = records.iter().map(Vec::len).sum();
            write_uint7(&mut header, total as u32);
            header.push(VERSION);
            header.push(0); // gflags: one parameter set, no stab, no reversal
            header.extend_from_slice(&[0, 0]); // context
            header.push(self.flags);
            header.push(max_sym);
            header.push(((self.qbits as u8) << 4) | self.qshift as u8);
            header.push((self.qloc as u8) << 4); // sloc = 0
            header.push(((self.ploc as u8) << 4) | self.dloc as u8);
            if self.flags & pflag::HAVE_PTAB != 0 {
                write_array(&mut header, &self.ptab);
            }
            if self.flags & pflag::HAVE_DTAB != 0 {
                write_array(&mut header, &self.dtab);
            }

            let mut rc = RangeEncoder::new();
            let mut len_models: Vec<ModelEncoder> =
                (0..4).map(|_| ModelEncoder::new(256)).collect();
            let mut qual: Vec<Option<ModelEncoder>> = (0..1 << 16).map(|_| None).collect();

            for (r, record) in records.iter().enumerate() {
                // `FIXED_LEN` means the length is stored once, for the first
                // record, and reused — see [`pflag::FIXED_LEN`].
                if self.flags & pflag::FIXED_LEN == 0 || r == 0 {
                    for (i, model) in len_models.iter_mut().enumerate() {
                        model.encode(&mut rc, ((record.len() >> (i * 8)) & 0xff) as u8);
                    }
                }
                let mut state = State::default();
                let mut ctx = 0u16;
                for (n, q) in record.iter().enumerate() {
                    qual[ctx as usize]
                        .get_or_insert_with(|| ModelEncoder::new(max_sym as usize + 1))
                        .encode(&mut rc, *q);
                    let pos = record.len() - n - 1;
                    ctx = self.update(*q, pos, &mut state);
                }
            }
            header.extend_from_slice(&rc.finish());
            header
        }

        /// The same arithmetic as [`Param::update`], over the encoder's own
        /// parameters — written out rather than shared, so that a mistake in
        /// one does not cancel a mistake in the other.
        fn update(&self, q: u8, pos: usize, state: &mut State) -> u16 {
            let mut ctx = 0u32;
            state.qctx = (state.qctx << self.qshift) + u32::from(q);
            ctx += (state.qctx & ((1 << self.qbits) - 1)) << self.qloc;
            if self.flags & pflag::HAVE_PTAB != 0 {
                ctx += u32::from(self.ptab[pos.min(1023)]) << self.ploc;
            }
            if self.flags & pflag::HAVE_DTAB != 0 {
                ctx += u32::from(self.dtab[state.delta.min(255) as usize]) << self.dloc;
                if state.prevq != q {
                    state.delta += 1;
                }
                state.prevq = q;
            }
            (ctx & 0xffff) as u16
        }
    }

    #[derive(Default)]
    struct State {
        qctx: u32,
        delta: u32,
        prevq: u8,
    }

    /// The inverse of [`read_array`], both levels of it.
    ///
    /// The one trap is that trailing zero runs must *not* be written: the
    /// decoder stops as soon as the run lengths account for `n` values, so
    /// anything after that is never consumed and every later field in the
    /// stream shifts by however many bytes were left behind.
    fn write_array(out: &mut Vec<u8>, table: &[u8]) {
        // Level one: how many inputs share each output value, long runs split.
        let mut runs: Vec<u8> = Vec::new();
        let highest = table.iter().copied().max().unwrap_or(0);
        for value in 0..=highest {
            let mut count = table.iter().filter(|v| **v == value).count();
            while count >= 255 {
                runs.push(255);
                count -= 255;
            }
            runs.push(count as u8);
        }
        // Level two: a repeat of a run length becomes a count of extra copies.
        let mut last: i32 = -1;
        let mut i = 0usize;
        while i < runs.len() {
            let run = runs[i];
            out.push(run);
            i += 1;
            if i32::from(run) == last {
                let mut copies = 0u8;
                while i < runs.len() && runs[i] == run && copies < 255 {
                    copies += 1;
                    i += 1;
                }
                out.push(copies);
            }
            last = i32::from(run);
        }
    }

    fn write_uint7(out: &mut Vec<u8>, value: u32) {
        let mut groups = Vec::new();
        let mut value = value;
        loop {
            groups.push((value & 0x7f) as u8);
            value >>= 7;
            if value == 0 {
                break;
            }
        }
        for (i, byte) in groups.iter().enumerate().rev() {
            out.push(if i == 0 { *byte } else { byte | 0x80 });
        }
    }

    fn roundtrip(records: &[Vec<u8>], encoder: &Encoder) {
        let expected: Vec<u8> = records.iter().flatten().copied().collect();
        let stream = encoder.encode(records);
        let decoded = decode(&stream, "test", 0).expect("decode");
        assert_eq!(decoded, expected);
    }

    #[test]
    fn every_lookup_table_survives_its_own_encoding() {
        for table in [
            (0..1024)
                .map(|i| (i / 16).min(15) as u8)
                .collect::<Vec<u8>>(),
            (0..256).map(|i| (i / 8).min(3) as u8).collect(),
            (0..256).map(|i| i as u8).collect(),
            vec![0u8; 1024],
            (0..1024).map(|i| u8::from(i > 900)).collect(),
        ] {
            let mut bytes = Vec::new();
            write_array(&mut bytes, &table);
            let mut reader = ByteReader::new(&bytes, "test", 0);
            let mut out = vec![0u8; table.len()];
            read_array(&mut reader, table.len(), &mut out).expect("read back");
            assert_eq!(out, table, "table of {} entries", table.len());
            // Nothing may be left over: the decoder reads the fields that
            // follow straight after this one.
            assert!(reader.is_empty(), "{} bytes unread", reader.remaining());
        }
    }

    #[test]
    fn a_block_of_equal_length_records_round_trips() {
        let records: Vec<Vec<u8>> = (0..20)
            .map(|r| (0..50).map(|i| ((i + r) % 40) as u8).collect())
            .collect();
        roundtrip(&records, &Encoder::simple());
    }

    #[test]
    fn records_of_different_lengths_round_trip() {
        let records: Vec<Vec<u8>> = (1..30)
            .map(|r| (0..r).map(|i| (i % 8) as u8).collect())
            .collect();
        roundtrip(&records, &Encoder::simple());
    }

    /// The flag §6.2 misnames `do_len`.
    ///
    /// Read as the specification describes it, this block decodes its first
    /// record and then takes the next record's four length bytes out of the
    /// quality stream — which is exactly what a `samtools` file of
    /// fixed-length reads does.
    #[test]
    fn a_fixed_length_block_stores_its_length_once() {
        let records: Vec<Vec<u8>> = (0..25)
            .map(|r| (0..40).map(|i| ((i + r) % 30) as u8).collect())
            .collect();

        let mut fixed = Encoder::simple();
        fixed.flags |= pflag::FIXED_LEN;
        roundtrip(&records, &fixed);

        // And it really did store one length rather than twenty-five: the
        // same records without the flag cost more bytes.
        let varying = Encoder::simple().encode(&records);
        assert!(
            fixed.encode(&records).len() < varying.len(),
            "the fixed-length stream should be the shorter of the two"
        );
    }

    #[test]
    fn the_position_and_delta_contexts_round_trip() {
        // Turning on ptab and dtab puts `pos` and `delta` into the context,
        // which is where the off-by-one in §6.2's main loop would show.
        let mut encoder = Encoder::simple();
        encoder.flags |= pflag::HAVE_PTAB | pflag::HAVE_DTAB;
        for (i, slot) in encoder.ptab.iter_mut().enumerate() {
            *slot = (i / 16).min(15) as u8;
        }
        for (i, slot) in encoder.dtab.iter_mut().enumerate() {
            *slot = (i / 8).min(3) as u8;
        }
        let records: Vec<Vec<u8>> = (0..15)
            .map(|r| (0..60).map(|i| ((i * 7 + r) % 12) as u8).collect())
            .collect();
        roundtrip(&records, &encoder);
    }

    #[test]
    fn a_single_quality_value_throughout_round_trips() {
        let records: Vec<Vec<u8>> = (0..10).map(|_| vec![30u8; 40]).collect();
        roundtrip(&records, &Encoder::simple());
    }

    #[test]
    fn a_block_of_a_wrong_version_is_refused() {
        let mut stream = Vec::new();
        write_uint7(&mut stream, 8);
        stream.push(4); // version 4 does not exist
        stream.extend_from_slice(&[0u8; 16]);
        let error = decode(&stream, "test", 0).expect_err("version 4");
        assert!(error.to_string().contains("version 4"), "{error}");
    }

    #[test]
    fn every_prefix_of_a_real_stream_fails_without_panicking() {
        let records: Vec<Vec<u8>> = (0..8)
            .map(|r| (0..30).map(|i| ((i + r) % 20) as u8).collect())
            .collect();
        let stream = Encoder::simple().encode(&records);
        for cut in 0..stream.len() {
            let _ = decode(&stream[..cut], "test", 0);
        }
    }
}