gwseq-io 0.2.1

//! CRAM records, rebuilt as BAM records.
//!
//! Format reference: `docs/cram_format_v3.1.md` §10.
//!
//! A CRAM record is not stored anywhere; it is *assembled* from one value taken
//! from each of twenty-odd data series, each of which lives in its own block
//! and is entropy-coded separately. This module takes one slice's worth of
//! those columns and writes out the rows — as BAM records, byte for byte what a
//! BAM file would have held, so that [`crate::bam::record::decode_block`] and
//! everything above it works unchanged. See the module doc in `cram/mod.rs` for
//! why that is the design.
//!
//! # Three things that are not local
//!
//! **Alignment starts are deltas.** `AP` is a difference from the previous
//! record's start, seeded with the slice's own start. So records must be
//! decoded in slice order and none can be skipped, however narrow the query.
//!
//! **Mates reach across records.** `NF` says how many records to skip forward
//! to this read's mate; the mate's reference, position and orientation are then
//! copied back, and `TLEN` is computed from the pair. That is a second pass
//! over the slice, after every record is decoded.
//!
//! **Sequences reach out of the file.** A mapped read stores only where it
//! *differs* from the reference. With no reference to hand, everything else —
//! flags, positions, CIGAR, names, tags, soft clips, inserted bases — still
//! decodes exactly; the matching stretches come back as `N`. See
//! [`super::reference`].

use std::sync::Arc;

use bytes::Bytes;

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

use super::compression::{tag_key, CompressionHeader};
use super::encoding::Streams;
use super::slice::Slice;

/// CRAM bit flags, §10.1's `CF` series.
mod cram_flag {
    /// Quality scores are stored as an array rather than as read features.
    pub const QUALITIES: i32 = 1;
    /// Mate information is stored verbatim rather than derived.
    pub const DETACHED: i32 = 2;
    /// The mate is further down this slice, `NF` records away.
    pub const MATE_DOWNSTREAM: i32 = 4;
    /// The sequence is unknown; any read features are there for the CIGAR.
    pub const NO_SEQUENCE: i32 = 8;
}

/// The `MF` series, §10.4 — the two SAM flags CRAM leaves out of `BF`.
mod mate_flag {
    pub const REVERSE: i32 = 1;
    pub const UNMAPPED: i32 = 2;
}

/// SAM flags this module reads or sets.
mod sam_flag {
    pub const UNMAPPED: u16 = 0x004;
    pub const MATE_UNMAPPED: u16 = 0x008;
    pub const REVERSE: u16 = 0x010;
    pub const MATE_REVERSE: u16 = 0x020;
}

/// BAM CIGAR operation codes, in the order `MIDNSHP=X` gives them.
mod cigar_op {
    pub const MATCH: u32 = 0;
    pub const INSERT: u32 = 1;
    pub const DELETE: u32 = 2;
    pub const SKIP: u32 = 3;
    pub const SOFT_CLIP: u32 = 4;
    pub const HARD_CLIP: u32 = 5;
    pub const PAD: u32 = 6;
}

/// Reference bases for the span one slice covers.
///
/// `start` is the 1-based reference position of `bases[0]`, matching the
/// coordinates the record decoder works in — CRAM's positions are 1-based
/// throughout, and only the BAM record written at the end is 0-based.
#[derive(Debug, Clone, Copy, Default)]
pub struct ReferenceBases<'a> {
    pub bases: &'a [u8],
    pub start: i64,
}

impl ReferenceBases<'_> {
    /// The base at a 1-based reference position.
    ///
    /// §11.3: "All out of range reference bases are all assumed to be 'N'." So
    /// is every base when no reference resolved at all, which is what makes a
    /// reference optional rather than required.
    #[inline]
    fn at(&self, position: i64) -> u8 {
        let index = position - self.start;
        if index < 0 {
            return b'N';
        }
        self.bases
            .get(index as usize)
            .copied()
            .map(|b| b.to_ascii_uppercase())
            .unwrap_or(b'N')
    }
}

/// Where the record decoder gets its reference bases.
///
/// The two cases are the two kinds of slice. A single-reference slice has one
/// window fetched before the decode, and every record uses it. A
/// **multi-reference** slice has a reference per record — the `RI` series says
/// which — so it carries the source itself and looks each one up as it goes.
///
/// Multi-reference slices are not a corner case to wave at: `samtools` puts one
/// wherever a container crosses a chromosome boundary, and the reads in it are
/// ordinary reads on ordinary chromosomes. Giving up on their sequences loses
/// the first several thousand reads of most chromosomes in a sorted file.
pub enum References<'a> {
    /// No reference resolved: everything that matched reads as `N`.
    None,
    /// One window, for a slice on a single reference.
    Fixed(ReferenceBases<'a>),
    /// A lookup per record, for a multi-reference slice.
    ByRefId {
        source: &'a super::reference::ReferenceSource,
        names: &'a [String],
    },
}

/// The window currently held while decoding a multi-reference slice.
struct HeldWindow {
    ref_id: i32,
    /// 0-based half-open span of `bases`.
    start: i64,
    end: i64,
    bases: Arc<Vec<u8>>,
}

/// One record, decoded but not yet serialised.
///
/// It exists because two things cannot be known until the whole slice is
/// decoded: a mate's position, and the template length computed from it.
#[derive(Debug, Default)]
struct Record {
    flags: u16,
    cram_flags: i32,
    ref_id: i32,
    /// 1-based, as CRAM stores it.
    start: i64,
    /// 1-based inclusive last reference base, or `start` when none is covered.
    end: i64,
    read_length: i32,
    read_group: i32,
    name: Vec<u8>,
    mapping_quality: u8,
    next_ref_id: i32,
    /// 1-based.
    next_pos: i64,
    template_length: i64,
    /// The record this one's mate is, once resolved. Forms a circular chain for
    /// templates of more than two segments.
    mate_line: Option<usize>,
    /// Set when the mate fields came from the file rather than from a mate.
    detached: bool,
    /// Whether the file stored an `RG` tag of its own, so the one rebuilt from
    /// the read-group series is not added on top.
    has_rg: bool,
    cigar: Vec<u32>,
    sequence: Vec<u8>,
    qualities: Vec<u8>,
    tags: Vec<u8>,
    /// Whether the file stored `MD` and `NM` itself. §10.5: a decoder that
    /// finds them verbatim must prefer those to its own.
    has_md: bool,
    has_nm: bool,
}

impl Record {
    /// Ready for the next record, keeping every buffer's memory.
    ///
    /// One record is reused for a whole slice, so every field has to be put
    /// back: some are written unconditionally by the decode and some — the
    /// mate fields, `has_md`, `has_rg` — only on the paths that find them.
    fn reset(&mut self) {
        self.name.clear();
        self.cigar.clear();
        self.sequence.clear();
        self.qualities.clear();
        self.tags.clear();
        self.flags = 0;
        self.cram_flags = 0;
        self.ref_id = 0;
        self.start = 0;
        self.end = 0;
        self.read_length = 0;
        self.read_group = 0;
        self.mapping_quality = 0;
        self.next_ref_id = -1;
        self.next_pos = 0;
        self.template_length = 0;
        self.mate_line = None;
        self.detached = false;
        self.has_rg = false;
        self.has_md = false;
        self.has_nm = false;
    }
}

/// What survives a record after it is written out.
///
/// A slice is decoded straight into its BAM buffer, one record at a time
/// through a single reused [`Record`], so nothing of a record is alive by the
/// time the next one starts — except what §10.4's second pass needs. These
/// are those fields, plus where the record's bytes begin, so the four the pass
/// changes can be written back in place.
#[derive(Debug)]
struct Placed {
    /// Where this record's BAM block starts in the slice's buffer.
    offset: usize,
    ref_id: i32,
    start: i64,
    end: i64,
    flags: u16,
    next_ref_id: i32,
    next_pos: i64,
    template_length: i64,
    mate_line: Option<usize>,
    detached: bool,
}

/// Offsets into a BAM record of the four fields mate resolution writes.
mod at {
    pub const FLAGS: usize = 18;
    pub const NEXT_REF_ID: usize = 24;
    pub const NEXT_POS: usize = 28;
    pub const TEMPLATE_LENGTH: usize = 32;
}

/// Buffers every record of a slice shares.
///
/// Each of these would otherwise be allocated and freed once per record —
/// tens of millions of times over a whole file, for buffers that are dead as
/// soon as the record is written out. Keeping them here costs one allocation
/// per slice apiece.
#[derive(Default)]
struct Scratch {
    /// Read-feature arrays: verbatim bases, qualities, inserts, soft clips.
    bytes: Vec<u8>,
    /// The `MD` string under construction.
    md: Vec<u8>,
    /// A read name generated from the record counter.
    generated: String,
}

/// The largest read this decoder will rebuild, in bases.
///
/// Sixteen mebibases is an order of magnitude past the longest read any
/// sequencer produces — an ONT ultra-long is a few megabases — and far below
/// anything that would abort the process.
const MAX_READ_LENGTH: usize = 1 << 24;

/// The largest number of records one slice may declare.
///
/// `samtools` writes ten thousand and §8.5 treats a hundred thousand as large.
/// A mebirecord is a ceiling, not a target.
const MAX_RECORDS_PER_SLICE: usize = 1 << 20;

/// What a slice's own bytes entitle it to ask for.
///
/// Every other length in this crate is bounded by the bytes behind it: a BAM
/// record's fields are in its block, and a codec's lengths are capped at
/// `MAX_CODEC_LEN`. A CRAM record is the exception, because it need not pay
/// for itself — every data series may be a one-symbol Huffman, which reads no
/// bits at all, and a mapped read with no features takes its bases from the
/// reference. So a 364-byte file can declare five million records, and an
/// 813-byte one a read of 2^30 bases, and both cost gigabytes.
///
/// This is the crate's "nothing allocated from a number a file names" rule
/// (`ARCHITECTURE.md` §12) reaching the one decoder it had not. The factor is
/// deliberately loose — a real slice carries at least a byte of quality per
/// base, so eight bases per byte is an order of magnitude of headroom — and
/// the floor keeps a legitimately tiny slice legal.
fn budget(slice: &Slice) -> usize {
    const PER_BYTE: usize = 8;
    const FLOOR: usize = 4096;
    let bytes = slice.core.len() + slice.external.iter().map(|(_, b)| b.len()).sum::<usize>();
    bytes.saturating_mul(PER_BYTE).max(FLOOR)
}

/// Decode every record of a slice into one buffer of BAM records.
///
/// The buffer is contiguous and handed out as one `Bytes`, so the records
/// sliced out of it share a single allocation — the same arrangement a
/// decompressed BGZF block already gives the BAM reader. Each record is
/// written into it as it is decoded, through one [`Record`] the whole slice
/// reuses, and [`patch_mates`] then fixes up the four fields §10.4 resolves
/// across records.
pub fn decode_slice(
    slice: &Slice,
    header: &CompressionHeader,
    references: References<'_>,
    read_groups: &[String],
    path: &str,
) -> Result<Bytes> {
    // Enough for the header of every record plus its variable parts, which is
    // the great majority of what gets written; the buffer grows for the rest.
    let mut out: Vec<u8> =
        Vec::with_capacity((slice.header.n_records.max(0) as usize * 96).min(1 << 24));
    let mut placed = decode_records(slice, header, references, read_groups, path, Some(&mut out))?;
    resolve_mates(&mut placed);
    patch_mates(&mut out, &placed);
    Ok(Bytes::from(out))
}

/// §10.4's second pass, written back into the records already serialised.
///
/// Only four fields can change, all of them fixed-width and at a fixed offset
/// from the start of a BAM record, which is what makes decoding straight into
/// the buffer possible at all.
fn patch_mates(out: &mut [u8], placed: &[Placed]) {
    for record in placed {
        let start = record.offset;
        out[start + at::FLAGS..start + at::FLAGS + 2].copy_from_slice(&record.flags.to_le_bytes());
        out[start + at::NEXT_REF_ID..start + at::NEXT_REF_ID + 4]
            .copy_from_slice(&record.next_ref_id.to_le_bytes());
        // CRAM is 1-based and BAM is 0-based; an unplaced mate is -1 either way.
        let next_pos = if record.next_ref_id < 0 && record.next_pos <= 0 {
            -1i32
        } else {
            (record.next_pos - 1) as i32
        };
        out[start + at::NEXT_POS..start + at::NEXT_POS + 4]
            .copy_from_slice(&next_pos.to_le_bytes());
        out[start + at::TEMPLATE_LENGTH..start + at::TEMPLATE_LENGTH + 4]
            .copy_from_slice(&(record.template_length as i32).to_le_bytes());
    }
}

/// Which references a slice's records actually sit on, and over what span.
///
/// Only a multi-reference slice needs this: a single-reference one says so in
/// its own header. It exists so that an index built by walking container
/// headers can file those slices per reference, as `samtools index` does,
/// rather than adding all of them to every query. §12 is explicit that this is
/// what indexing a multi-reference container costs: "the exception to this is
/// with multi-reference containers, where the `RI` data series must be read".
///
/// No reference is resolved — the spans come from `RI`, `AP` and the read
/// features, none of which need one.
pub fn slice_spans(
    slice: &Slice,
    header: &CompressionHeader,
    path: &str,
) -> Result<Vec<(i32, i64, i64)>> {
    let records = decode_records(slice, header, References::None, &[], path, None)?;
    let mut spans: Vec<(i32, i64, i64)> = Vec::new();
    for record in &records {
        if record.ref_id < 0 {
            continue;
        }
        // 0-based half-open, as an index entry is.
        let (start, end) = ((record.start - 1).max(0), record.end.max(record.start));
        match spans.iter_mut().find(|(id, ..)| *id == record.ref_id) {
            Some(span) => {
                span.1 = span.1.min(start);
                span.2 = span.2.max(end);
            }
            None => spans.push((record.ref_id, start, end)),
        }
    }
    Ok(spans)
}

/// Decode a slice's records, writing each into `out` as it is finished.
///
/// `out` is `None` for [`slice_spans`], which wants the positions and not the
/// records; nothing else skips it.
fn decode_records(
    slice: &Slice,
    header: &CompressionHeader,
    references: References<'_>,
    read_groups: &[String],
    path: &str,
    mut out: Option<&mut Vec<u8>>,
) -> Result<Vec<Placed>> {
    let mut held: Option<HeldWindow> = None;
    let mut streams = Streams::new(&slice.core, slice.streams(), path);
    let series = &header.series;
    let preservation = &header.preservation;
    let n_records = slice.header.n_records as usize;
    let budget = budget(slice);
    if slice.header.n_records < 0 || n_records > MAX_RECORDS_PER_SLICE || n_records > budget {
        return Err(Error::corrupt(
            path,
            slice.offset,
            format!(
                "a slice of {} records, which its {} bytes of data cannot hold",
                slice.header.n_records,
                budget / 8
            ),
        ));
    }
    // What the read lengths may add up to over the whole slice. Bounding each
    // read alone is not enough: a million reads of a thousand bases is the
    // same bomb spread out.
    let mut remaining = budget;

    let mut placed: Vec<Placed> = Vec::with_capacity(n_records.min(1 << 20));
    // One record and one set of buffers, reused by every record in the slice.
    let mut record = Record::default();
    let mut scratch = Scratch::default();
    // §10.2: the first record deltas against the slice's own alignment start,
    // which is zero for a multi-reference or unmapped slice.
    let mut last_position = i64::from(slice.header.start);

    for index in 0..n_records {
        record.reset();
        record.flags = series.bf.decode_int(&mut streams)? as u16;
        record.cram_flags = series.cf.decode_int(&mut streams)?;

        // Positions, §10.2.
        record.ref_id = if slice.header.is_multi_ref() {
            series.ri.decode_int(&mut streams)?
        } else {
            slice.header.ref_id
        };
        record.read_length = series.rl.decode_int(&mut streams)?;
        if record.read_length < 0 || record.read_length as usize > MAX_READ_LENGTH {
            return Err(Error::corrupt(
                path,
                slice.offset,
                format!(
                    "record {index} of this slice is a read of {} bases",
                    record.read_length
                ),
            ));
        }
        // Charged against the slice's budget, so the sum is bounded as well as
        // each one. A mapped read with no features costs the file nothing per
        // base — the bases come from the reference — so without this a slice
        // of a few hundred bytes can ask for gigabytes of sequence.
        remaining = remaining
            .checked_sub(record.read_length as usize)
            .ok_or_else(|| {
                Error::corrupt(
                    path,
                    slice.offset,
                    format!(
                        "a slice whose reads total more than the {budget} bases its data can hold"
                    ),
                )
            })?;
        record.start = if preservation.ap_delta {
            last_position + i64::from(series.ap.decode_int(&mut streams)?)
        } else {
            i64::from(series.ap.decode_int(&mut streams)?)
        };
        last_position = record.start;
        record.read_group = series.rg.decode_int(&mut streams)?;

        // Names, §10.3.
        if preservation.read_names_included {
            series
                .rn
                .decode_array(&mut streams, None, &mut record.name)?;
        }

        // Mates, §10.4.
        record.next_ref_id = -1;
        record.next_pos = 0;
        if record.cram_flags & cram_flag::DETACHED != 0 {
            record.detached = true;
            let mate_flags = series.mf.decode_int(&mut streams)?;
            if mate_flags & mate_flag::REVERSE != 0 {
                record.flags |= sam_flag::MATE_REVERSE;
            }
            if mate_flags & mate_flag::UNMAPPED != 0 {
                record.flags |= sam_flag::MATE_UNMAPPED;
            }
            // A detached record carries its name even when the file otherwise
            // does not preserve them, because the pair spans slices and the
            // generated names would not agree.
            if !preservation.read_names_included {
                series
                    .rn
                    .decode_array(&mut streams, None, &mut record.name)?;
            }
            record.next_ref_id = series.ns.decode_int(&mut streams)?;
            record.next_pos = i64::from(series.np.decode_int(&mut streams)?);
            record.template_length = i64::from(series.ts.decode_int(&mut streams)?);
        } else if record.cram_flags & cram_flag::MATE_DOWNSTREAM != 0 {
            let skip = series.nf.decode_int(&mut streams)?;
            if skip < 0 {
                return Err(Error::corrupt(
                    path,
                    0,
                    format!("a mate {skip} records downstream"),
                ));
            }
            // "zero meaning the next fragment is also the next record".
            let mate = index + skip as usize + 1;
            if mate >= n_records {
                return Err(Error::corrupt(
                    path,
                    0,
                    format!("a mate at record {mate} of a slice holding {n_records}"),
                ));
            }
            record.mate_line = Some(mate);
        }

        // Auxiliary tags, §10.5.
        decode_tags(&mut streams, header, &mut record, path, slice.offset)?;
        // §10.5: the RG tag is rebuilt from the read group series when the
        // encoder dropped it — which `htslib` always does, but the format does
        // not require it. An encoder that kept the tag would otherwise get two.
        if record.read_group >= 0 && !record.has_rg {
            if let Some(name) = read_groups.get(record.read_group as usize) {
                record.tags.extend_from_slice(b"RGZ");
                record.tags.extend_from_slice(name.as_bytes());
                record.tags.push(0);
            }
        }

        // The reference this record needs, which for a multi-reference slice
        // is not the one the record before it needed.
        let bases = resolve_bases(&references, &mut held, &record)?;
        // §10.5: `MD` and `NM` may be rebuilt during reference-based
        // reconstruction, and `samtools` does. With no reference there is
        // nothing to rebuild them from, and inventing them from `N`s would be
        // worse than leaving them out.
        let has_reference = !bases.bases.is_empty();

        // The sequence, §10.6 or §10.7.
        if record.flags & sam_flag::UNMAPPED == 0 {
            decode_mapped(
                &mut streams,
                header,
                &mut record,
                bases,
                has_reference,
                &mut scratch,
                path,
                slice.offset,
            )?;
        } else {
            decode_unmapped(&mut streams, header, &mut record, path)?;
        }
        let offset = match out.as_deref_mut() {
            Some(out) => serialise(out, &record, index, slice, &mut scratch, path)?,
            None => 0,
        };
        placed.push(Placed {
            offset,
            ref_id: record.ref_id,
            start: record.start,
            end: record.end,
            flags: record.flags,
            next_ref_id: record.next_ref_id,
            next_pos: record.next_pos,
            template_length: record.template_length,
            mate_line: record.mate_line,
            detached: record.detached,
        });
    }

    Ok(placed)
}

/// The reference bases for one record, refreshing the held window when the
/// record has moved off it.
fn resolve_bases<'a>(
    references: &'a References<'_>,
    held: &'a mut Option<HeldWindow>,
    record: &Record,
) -> Result<ReferenceBases<'a>> {
    match references {
        References::None => Ok(ReferenceBases::default()),
        References::Fixed(bases) => Ok(*bases),
        References::ByRefId { source, names } => {
            if record.ref_id < 0 {
                return Ok(ReferenceBases::default());
            }
            // A window wide enough for any plausible alignment from here. The
            // source reads out to a whole window in any case, so the end asked
            // for only has to be past this record.
            let start = (record.start - 1).max(0);
            let end = start + i64::from(record.read_length).max(1) + 1;
            let stale = match held.as_ref() {
                Some(window) => {
                    window.ref_id != record.ref_id || window.start > start || window.end < end
                }
                None => true,
            };
            if stale {
                *held = None;
                let Some(name) = names.get(record.ref_id as usize) else {
                    return Ok(ReferenceBases::default());
                };
                if source.has(name) {
                    let (bases, from) = source.window(name, start, end)?;
                    let length = bases.len() as i64;
                    *held = Some(HeldWindow {
                        ref_id: record.ref_id,
                        start: from,
                        end: from + length,
                        bases,
                    });
                }
            }
            Ok(match held.as_ref() {
                Some(window) => ReferenceBases {
                    bases: &window.bases[..],
                    // 1-based, as the record decoder works in.
                    start: window.start + 1,
                },
                None => ReferenceBases::default(),
            })
        }
    }
}

/// §10.5: a tag-list index, then one byte array per tag in that list.
fn decode_tags(
    streams: &mut Streams<'_>,
    header: &CompressionHeader,
    record: &mut Record,
    path: &str,
    at: u64,
) -> Result<()> {
    if header.series.tl.is_null() {
        return Ok(());
    }
    let line = header.series.tl.decode_int(streams)?;
    // Borrowed, not cloned: `header` and `record` are separate bindings, so
    // there is no conflict to allocate around — and this ran ten million
    // times on the test file.
    let list = header.tag_list(line, path)?;
    if !list.is_empty() {
        // One allocation for the whole tag block rather than one per doubling:
        // a typical `samtools` record carries four tags in about fifty bytes,
        // which a `Vec` starting empty reaches in four reallocations.
        record.tags.reserve(64);
    }
    for tag in list {
        let tag = *tag;
        let encoding = header.tags.get(&tag_key(tag)).ok_or_else(|| {
            Error::corrupt(
                path,
                at,
                format!(
                    "a record carries tag {}{}:{} and the tag encoding map does not",
                    tag[0] as char, tag[1] as char, tag[2] as char
                ),
            )
        })?;
        // The value bytes are already in BAM's own layout — including the
        // terminating nul of a `Z` or `H`, which the encoder wrote into the
        // stream before its delimiter. So the tag is its three key bytes
        // followed by whatever comes back, with nothing added.
        match &tag[..2] {
            b"MD" => record.has_md = true,
            b"NM" => record.has_nm = true,
            // Recorded here rather than by scanning the finished tag bytes,
            // which cannot tell a key from a value that looks like one.
            b"RG" => record.has_rg = true,
            _ => {}
        }
        record.tags.extend_from_slice(&tag);
        encoding.decode_array(streams, None, &mut record.tags)?;
    }
    Ok(())
}

/// Builds the `MD` string and the `NM` count as the reference is walked.
///
/// SAMtags: `MD` alternates runs of matching bases with the *reference* bases
/// that differ, and deletions as `^` followed by what was deleted. Clipping,
/// padding, reference skips and insertions do not appear in it. `NM` is the
/// edit distance: mismatches plus inserted plus deleted bases.
///
/// Numbers alternate with the other items, so two mismatches in a row are
/// separated by a `0`. Falling out of that rule is the classic way to produce
/// an `MD` that parses and means something else.
#[derive(Debug, Default)]
struct EditString {
    md: Vec<u8>,
    nm: i64,
    /// Matching bases seen since the last item was written.
    run: i64,
}

impl EditString {
    /// Over a buffer the slice lends it, emptied first. `MD` for a read that
    /// matches outright is its length in decimal, so the buffer is small and
    /// reusing it saves an allocation per record.
    fn new(md: Vec<u8>) -> Self {
        let mut md = md;
        md.clear();
        Self { md, nm: 0, run: 0 }
    }

    /// Close the current run of matches, which every non-match item must do
    /// before writing itself.
    fn flush_run(&mut self) {
        itoa(&mut self.md, self.run);
        self.run = 0;
    }

    fn matched(&mut self, count: i64) {
        self.run += count;
    }

    /// A read base that differs from the reference base under it.
    fn mismatch(&mut self, reference_base: u8) {
        self.flush_run();
        self.md.push(reference_base);
        self.nm += 1;
    }

    /// One aligned base, which is a match or a mismatch depending.
    fn aligned(&mut self, read_base: u8, reference_base: u8) {
        if read_base.to_ascii_uppercase() == reference_base {
            self.matched(1);
        } else {
            self.mismatch(reference_base);
        }
    }

    fn deleted(&mut self, bases: &[u8]) {
        self.flush_run();
        self.md.push(b'^');
        self.md.extend_from_slice(bases);
        self.nm += bases.len() as i64;
    }

    fn inserted(&mut self, count: i64) {
        self.nm += count;
    }

    /// The finished `MD`, which always ends in a number.
    fn finish(mut self) -> (Vec<u8>, i64) {
        self.flush_run();
        (self.md, self.nm)
    }
}

/// A non-negative integer as decimal bytes, appended.
///
/// Appended rather than returned because `MD` writes one of these per run of
/// matches, and a `Vec` per run is an allocation and a free per mismatch of
/// every record in the file — tens of millions of them, of a few bytes each,
/// for no purpose but to be copied and dropped.
fn itoa(out: &mut Vec<u8>, value: i64) {
    if value == 0 {
        out.push(b'0');
        return;
    }
    // Twenty digits is the widest `i64`, and this is only ever called with a
    // non-negative one.
    let mut digits = [0u8; 20];
    let mut n = value as u64;
    let mut at = digits.len();
    while n > 0 {
        at -= 1;
        digits[at] = b'0' + (n % 10) as u8;
        n /= 10;
    }
    out.extend_from_slice(&digits[at..]);
}

/// §10.6: read features, then the mapping quality and the qualities.
#[allow(clippy::too_many_arguments)]
fn decode_mapped(
    streams: &mut Streams<'_>,
    header: &CompressionHeader,
    record: &mut Record,
    reference: ReferenceBases<'_>,
    has_reference: bool,
    scratch: &mut Scratch,
    path: &str,
    at: u64,
) -> Result<()> {
    let series = &header.series;
    let matrix = &header.preservation.substitution_matrix;
    let read_length = record.read_length as usize;

    let n_features = series.fn_.decode_int(streams)?;
    // Two per base, not one. A `Q` feature sets the quality of a position and
    // moves neither cursor, so a file preserving mismatch qualities writes a
    // `Q` *and* an `X` at the same position — both legal, and `read_length + 1`
    // refuses it. The bound is here only to stop a record claiming more
    // features than the bases it has; the streams bound it in any case unless
    // `FC`/`FP` are constants, which the slice budget now covers.
    if n_features < 0 || n_features as usize > 2 * read_length + 2 {
        return Err(Error::corrupt(
            path,
            at,
            format!("a read of {read_length} bases carrying {n_features} read features"),
        ));
    }

    // Taken from the record and put back at the end: the buffer, and its
    // memory, are the slice's rather than this record's.
    let mut sequence = std::mem::take(&mut record.sequence);
    sequence.clear();
    sequence.reserve(read_length);
    // A quality of 0xFF is BAM's "not stored"; the `Q` and `B` features
    // overwrite individual positions when the file kept only some.
    let mut qualities = std::mem::take(&mut record.qualities);
    qualities.clear();
    qualities.resize(read_length, 0xFF);
    let mut cigar: Vec<u32> = std::mem::take(&mut record.cigar);
    cigar.clear();
    let mut open_op = u32::MAX;
    let mut open_len = 0u32;

    /// Close the CIGAR operation being accumulated and start another.
    fn push_op(cigar: &mut Vec<u32>, open_op: &mut u32, open_len: &mut u32, op: u32, len: u32) {
        if len == 0 {
            return;
        }
        if *open_op == op {
            *open_len += len;
            return;
        }
        if *open_len > 0 {
            cigar.push((*open_len << 4) | *open_op);
        }
        *open_op = op;
        *open_len = len;
    }

    // §10.5's optional rebuild, done only when there is a reference to do it
    // from and the file did not store the tags itself.
    // Not for a record the file says has no bases: walking its features
    // against the reference produces an `MD:Z:50 NM:i:0` describing a sequence
    // that is not there, and `serialise` then drops the sequence and keeps the
    // tags. Nothing to reconstruct means nothing to describe.
    let no_sequence = record.cram_flags & cram_flag::NO_SEQUENCE != 0;
    let mut edits = (has_reference && !no_sequence && !(record.has_md && record.has_nm))
        .then(|| EditString::new(std::mem::take(&mut scratch.md)));

    // Both 1-based: where we are in the read, and where on the reference.
    let mut read_pos: i64 = 1;
    let mut ref_pos: i64 = record.start;
    let mut feature_pos: i64 = 0;

    for _ in 0..n_features {
        let code = series.fc.decode_byte(streams)?;
        feature_pos += i64::from(series.fp.decode_int(streams)?);
        let position = feature_pos;
        // `Q` and `q` set the quality of a position without consuming a base,
        // so they may name the position the feature before them just consumed:
        // after an `X` at p, `read_pos` is p+1 and "the quality of the base
        // just substituted" is a `Q` at p. Every other feature must be at or
        // past the cursor.
        // Never below 1: the two handlers that use `position` index
        // `position - 1`, so a position of zero would underflow the index
        // rather than miss it.
        let earliest = if matches!(code, b'Q' | b'q') {
            (read_pos - 1).max(1)
        } else {
            read_pos
        };
        if position < earliest || position > read_length as i64 + 1 {
            return Err(Error::corrupt(
                path,
                at,
                format!(
                    "a read feature at read position {position}, which is outside a \
                     {read_length}-base read read forwards from {read_pos}"
                ),
            ));
        }
        // Everything between the last feature and this one matched the
        // reference, so it is copied from there.
        let gap = position - read_pos;
        if gap > 0 {
            for offset in 0..gap {
                sequence.push(reference.at(ref_pos + offset));
            }
            push_op(
                &mut cigar,
                &mut open_op,
                &mut open_len,
                cigar_op::MATCH,
                gap as u32,
            );
            if let Some(edits) = &mut edits {
                edits.matched(gap);
            }
            read_pos += gap;
            ref_pos += gap;
        }

        match code {
            // A base substitution: the reference base, changed by the code.
            b'X' => {
                let code = series.bs.decode_byte(streams)?;
                // A substitution names a base *relative to the reference*, so
                // with no reference there is no base to name. Substituting
                // anyway picks row `N` of the matrix, whose four codes map to
                // `A`, `C`, `G` and `T` — a concrete base with no relation to
                // the read's, presented as data. The README and ARCHITECTURE
                // §14.2 both promise Ns here, and a caller counting mismatches
                // has no way to tell a fabricated base from a real one.
                //
                // §11.3's "all out of range reference bases are assumed to be
                // 'N'" is a different case and still substitutes: there the
                // reference genuinely reads `N` and the matrix row is the
                // specified answer. `has_reference` separates the two.
                sequence.push(if has_reference {
                    matrix.substitute(reference.at(ref_pos), code)
                } else {
                    b'N'
                });
                if let Some(edits) = &mut edits {
                    edits.mismatch(reference.at(ref_pos));
                }
                push_op(&mut cigar, &mut open_op, &mut open_len, cigar_op::MATCH, 1);
                read_pos += 1;
                ref_pos += 1;
            }
            // A base and its quality, stored verbatim.
            b'B' => {
                let base = series.ba.decode_byte(streams)?;
                let quality = series.qs.decode_byte(streams)?;
                sequence.push(base);
                if let Some(slot) = qualities.get_mut(read_pos as usize - 1) {
                    *slot = quality;
                }
                if let Some(edits) = &mut edits {
                    edits.aligned(base, reference.at(ref_pos));
                }
                push_op(&mut cigar, &mut open_op, &mut open_len, cigar_op::MATCH, 1);
                read_pos += 1;
                ref_pos += 1;
            }
            // A stretch of bases, stored verbatim.
            b'b' => {
                scratch.bytes.clear();
                series.bb.decode_array(streams, None, &mut scratch.bytes)?;
                let len = scratch.bytes.len() as i64;
                sequence.extend_from_slice(&scratch.bytes);
                if let Some(edits) = &mut edits {
                    for (offset, base) in scratch.bytes.iter().enumerate() {
                        edits.aligned(*base, reference.at(ref_pos + offset as i64));
                    }
                }
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::MATCH,
                    len as u32,
                );
                read_pos += len;
                ref_pos += len;
            }
            // A stretch of qualities. Moves neither cursor.
            b'q' => {
                scratch.bytes.clear();
                series.qq.decode_array(streams, None, &mut scratch.bytes)?;
                for (offset, quality) in scratch.bytes.iter().enumerate() {
                    if let Some(slot) = qualities.get_mut(position as usize - 1 + offset) {
                        *slot = *quality;
                    }
                }
            }
            // One quality. Moves neither cursor.
            b'Q' => {
                let quality = series.qs.decode_byte(streams)?;
                if let Some(slot) = qualities.get_mut(position as usize - 1) {
                    *slot = quality;
                }
            }
            // Inserted bases: the read advances, the reference does not.
            b'I' => {
                scratch.bytes.clear();
                series.in_.decode_array(streams, None, &mut scratch.bytes)?;
                let len = scratch.bytes.len() as i64;
                sequence.extend_from_slice(&scratch.bytes);
                if let Some(edits) = &mut edits {
                    edits.inserted(len);
                }
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::INSERT,
                    len as u32,
                );
                read_pos += len;
            }
            b'i' => {
                let base = series.ba.decode_byte(streams)?;
                sequence.push(base);
                if let Some(edits) = &mut edits {
                    edits.inserted(1);
                }
                push_op(&mut cigar, &mut open_op, &mut open_len, cigar_op::INSERT, 1);
                read_pos += 1;
            }
            // Soft-clipped bases: in the read, not aligned to the reference.
            b'S' => {
                scratch.bytes.clear();
                series.sc.decode_array(streams, None, &mut scratch.bytes)?;
                let len = scratch.bytes.len() as i64;
                sequence.extend_from_slice(&scratch.bytes);
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::SOFT_CLIP,
                    len as u32,
                );
                read_pos += len;
            }
            // Deletions and skips: the reference advances, the read does not.
            b'D' => {
                let len = i64::from(series.dl.decode_int(streams)?.max(0));
                if let Some(edits) = &mut edits {
                    let deleted: Vec<u8> = (0..len)
                        .map(|offset| reference.at(ref_pos + offset))
                        .collect();
                    edits.deleted(&deleted);
                }
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::DELETE,
                    len as u32,
                );
                ref_pos += len;
            }
            b'N' => {
                let len = series.rs.decode_int(streams)?;
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::SKIP,
                    len.max(0) as u32,
                );
                ref_pos += i64::from(len.max(0));
            }
            // Padding and hard clips: neither cursor moves. Hard-clipped bases
            // are not in the read at all, which is what distinguishes them from
            // soft clips.
            b'P' => {
                let len = series.pd.decode_int(streams)?;
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::PAD,
                    len.max(0) as u32,
                );
            }
            b'H' => {
                let len = series.hc.decode_int(streams)?;
                push_op(
                    &mut cigar,
                    &mut open_op,
                    &mut open_len,
                    cigar_op::HARD_CLIP,
                    len.max(0) as u32,
                );
            }
            other => {
                return Err(Error::corrupt(
                    path,
                    0,
                    format!(
                        "read feature code {:?} ({other:#04x}), which §10.6 does not define",
                        other as char
                    ),
                ))
            }
        }
    }

    // Whatever is left of the read matched the reference.
    let remaining = read_length as i64 - (read_pos - 1);
    if remaining > 0 {
        for offset in 0..remaining {
            sequence.push(reference.at(ref_pos + offset));
        }
        push_op(
            &mut cigar,
            &mut open_op,
            &mut open_len,
            cigar_op::MATCH,
            remaining as u32,
        );
        if let Some(edits) = &mut edits {
            edits.matched(remaining);
        }
        ref_pos += remaining;
    }
    if open_len > 0 {
        cigar.push((open_len << 4) | open_op);
    }

    record.mapping_quality = series.mq.decode_int(streams)? as u8;
    if record.cram_flags & cram_flag::QUALITIES != 0 {
        qualities.clear();
        series
            .qs
            .decode_array(streams, Some(read_length), &mut qualities)?;
    }

    if let Some(edits) = edits {
        let (mut md, nm) = edits.finish();
        if !record.has_md {
            record.tags.extend_from_slice(b"MDZ");
            record.tags.extend_from_slice(&md);
            record.tags.push(0);
        }
        if !record.has_nm {
            // BAM types an integer by width; `samtools` writes NM as the
            // narrowest unsigned type that holds it, which for an edit distance
            // is nearly always one byte.
            if (0..=255).contains(&nm) {
                record.tags.extend_from_slice(b"NMC");
                record.tags.push(nm as u8);
            } else {
                record.tags.extend_from_slice(b"NMi");
                record.tags.extend_from_slice(&(nm as i32).to_le_bytes());
            }
        }
        // Back to the slice, with its memory.
        md.clear();
        scratch.md = md;
    }

    record.cigar = cigar;
    record.sequence = sequence;
    record.qualities = qualities;
    // The last reference base the alignment covers, 1-based and inclusive. A
    // read consuming no reference reports its own start, which is what the BAM
    // reader's filter already expects of a placed unmapped read.
    record.end = (ref_pos - 1).max(record.start);
    Ok(())
}

/// §10.7: bases and qualities, with no reference and no features.
fn decode_unmapped(
    streams: &mut Streams<'_>,
    header: &CompressionHeader,
    record: &mut Record,
    _path: &str,
) -> Result<()> {
    let series = &header.series;
    let read_length = record.read_length as usize;
    // The record's own buffers, emptied and refilled — the slice reuses one
    // record, so these are the same allocations every record here uses.
    let mut sequence = std::mem::take(&mut record.sequence);
    sequence.clear();
    sequence.reserve(read_length.min(1 << 16));
    series
        .ba
        .decode_array(streams, Some(read_length), &mut sequence)?;
    let mut qualities = std::mem::take(&mut record.qualities);
    qualities.clear();
    if record.cram_flags & cram_flag::QUALITIES != 0 {
        series
            .qs
            .decode_array(streams, Some(read_length), &mut qualities)?;
    } else {
        qualities.resize(read_length, 0xFF);
    }
    record.sequence = sequence;
    record.qualities = qualities;
    record.end = record.start;
    Ok(())
}

/// §10.4's second pass: close each mate chain, copy mate fields across it, and
/// compute the template length from its extremes.
///
/// A template may hold more than two segments, in which case the mate links
/// form a circle. The walk below is bounded by the record count, so a malformed
/// chain gives up rather than spinning.
fn resolve_mates(records: &mut [Placed]) {
    let n = records.len();
    // Every record a walk has already reached. Two jobs: a well-formed chain
    // is walked once rather than once per member, and a malformed one — two
    // records naming the same mate, which nothing in the format prevents —
    // cannot be walked again from each of them. Without it every start costs
    // up to `n` steps, so a slice of a million records shaped that way is
    // 10^12 steps; with it each step reaches a record no step has reached, so
    // the whole pass is linear.
    let mut visited = vec![false; n];
    for start in 0..n {
        if visited[start] || records[start].mate_line.is_none() || records[start].detached {
            continue;
        }
        // Walk the chain, noting its extremes and closing it back on itself.
        let reference_id = records[start].ref_id;
        let mut leftmost = records[start].start;
        let mut rightmost = records[start].end;
        let mut left_count = 0usize;
        let mut chain = vec![start];
        let mut at = start;
        visited[start] = true;
        for _ in 0..n {
            let record = &records[at];
            if record.start < leftmost {
                leftmost = record.start;
                left_count = 1;
            } else if record.start == leftmost {
                left_count += 1;
            }
            if record.end > rightmost {
                rightmost = record.end;
            }
            match records[at].mate_line {
                // `next == start` closes a well-formed circle, and is the
                // commonest way this stops; any other visited record means the
                // chain runs into one already resolved, which only a malformed
                // slice does.
                Some(next) if !visited[next] => {
                    visited[next] = true;
                    at = next;
                    chain.push(at);
                }
                Some(_) => break,
                None => {
                    // The last segment closes the circle back to the first.
                    records[at].mate_line = Some(start);
                    break;
                }
            }
        }

        // Mate fields come from the next record in the chain — except for a
        // record that stored them verbatim, which is the whole point of being
        // detached. §10.4: an encoder marks a record detached even inside one
        // slice "to preserve inconsistent data so that it round-trips through
        // the CRAM format with full fidelity", and a pair where one read
        // contains the other is exactly that: the derived template length is
        // zero and the original was not.
        for &index in &chain {
            if records[index].detached {
                continue;
            }
            let Some(mate) = records[index].mate_line else {
                continue;
            };
            let (mate_ref, mate_pos, mate_flags) = {
                let m = &records[mate];
                (m.ref_id, m.start, m.flags)
            };
            let record = &mut records[index];
            record.next_ref_id = mate_ref;
            record.next_pos = mate_pos;
            if mate_flags & sam_flag::REVERSE != 0 {
                record.flags |= sam_flag::MATE_REVERSE;
            }
            if mate_flags & sam_flag::UNMAPPED != 0 {
                record.flags |= sam_flag::MATE_UNMAPPED;
            }
        }

        // §10.4, and the SAM specification behind it: the leftmost segment
        // gets a positive template length and the rightmost a negative one.
        // An unmapped chain has no length to speak of, and neither has a
        // template of one segment — which is what a record naming a mate some
        // other record already claimed comes out as.
        if reference_id < 0 || chain.len() < 2 {
            continue;
        }
        // Nor has one whose segments sit on different references: the extremes
        // below would subtract positions on two sequences, which is a number
        // with no meaning. SAM's convention is zero and `htslib` writes zero.
        // An encoder marks such a pair detached precisely because the derived
        // value would not round-trip, but a reader may not depend on that.
        if chain
            .iter()
            .any(|&index| records[index].ref_id != reference_id)
        {
            continue;
        }
        let length = rightmost - leftmost + 1;
        // "The sign of segments in the middle is undefined" — but a template of
        // *two* segments has no middle, and that distinction is load-bearing.
        // A pair where one read wholly contains the other has a single segment
        // that is both leftmost and rightmost, and the other is neither; read
        // as "leftmost, rightmost, or else zero" that second read comes out
        // zero, where every other tool gives it the negative. Six reads in ten
        // million of a real file are shaped like this, all of them a read with
        // an indel against a mate that swallows it.
        let is_pair = chain.len() == 2;
        let mut left_remaining = left_count;
        for &index in &chain {
            let record = &mut records[index];
            if record.detached {
                // Its template length came from the file. A detached record
                // still counts towards the extremes above — the records that do
                // derive theirs need it — but nothing here may overwrite it.
                continue;
            }
            record.template_length = if record.start == leftmost && left_remaining > 0 {
                left_remaining -= 1;
                length
            } else if record.end == rightmost || is_pair {
                -length
            } else {
                0
            };
        }
    }
}

/// BAM's 4-bit sequence alphabet, `=ACMGRSVTWYHKDBN`, indexed by ASCII.
static BASE_CODES: [u8; 256] = {
    let mut table = [15u8; 256]; // anything unrecognised is N
    table[b'=' as usize] = 0;
    table[b'A' as usize] = 1;
    table[b'C' as usize] = 2;
    table[b'M' as usize] = 3;
    table[b'G' as usize] = 4;
    table[b'R' as usize] = 5;
    table[b'S' as usize] = 6;
    table[b'V' as usize] = 7;
    table[b'T' as usize] = 8;
    table[b'W' as usize] = 9;
    table[b'Y' as usize] = 10;
    table[b'H' as usize] = 11;
    table[b'K' as usize] = 12;
    table[b'D' as usize] = 13;
    table[b'B' as usize] = 14;
    table[b'N' as usize] = 15;
    table[b'a' as usize] = 1;
    table[b'c' as usize] = 2;
    table[b'm' as usize] = 3;
    table[b'g' as usize] = 4;
    table[b'r' as usize] = 5;
    table[b's' as usize] = 6;
    table[b'v' as usize] = 7;
    table[b't' as usize] = 8;
    table[b'w' as usize] = 9;
    table[b'y' as usize] = 10;
    table[b'h' as usize] = 11;
    table[b'k' as usize] = 12;
    table[b'd' as usize] = 13;
    table[b'b' as usize] = 14;
    table
};

/// The BAI bin a record declares itself in, from the SAM specification's own
/// `reg2bin`.
///
/// The specification writes each level's first bin as `((1 << k) - 1) / 7`;
/// those are worked out here, because the expression is five constants in a
/// trench coat and one of them is `7 / 7`.
fn reg2bin(start: i64, end: i64) -> u16 {
    /// `(shift, first bin at that level)`, finest level first.
    const LEVELS: [(u32, i64); 5] = [(14, 4681), (17, 585), (20, 73), (23, 9), (26, 1)];
    let beg = start;
    // A record covering no reference sits in the bin its own position does.
    let end = (end - 1).max(start);
    for (shift, first) in LEVELS {
        if beg >> shift == end >> shift {
            // Wrapped rather than checked: a position past the 2^29 the bin
            // scheme covers is a file's to name, and this is a reported field
            // rather than one anything is looked up by.
            return (first + (beg >> shift)) as u16;
        }
    }
    0
}

/// Write one record as BAM, appended, and say where it starts.
///
/// Called as each record finishes rather than over a finished `Vec<Record>`:
/// the only fields §10.4's second pass can still change are four of fixed
/// width at a fixed offset, so the rest can be written once and left alone.
fn serialise(
    out: &mut Vec<u8>,
    record: &Record,
    index: usize,
    slice: &Slice,
    scratch: &mut Scratch,
    path: &str,
) -> Result<usize> {
    let offset = out.len();
    {
        // §10.3: names may be dropped, and are then rebuilt from the slice's
        // record counter — which is what makes them stable across a re-read.
        let name: &[u8] = if record.name.is_empty() {
            scratch.generated.clear();
            use std::fmt::Write as _;
            let _ = write!(
                scratch.generated,
                "{}",
                slice.header.record_counter + index as i64 + 1
            );
            scratch.generated.as_bytes()
        } else {
            &record.name
        };
        // BAM's read name is nul-terminated and its length byte counts the nul.
        let name_len = name.len() + 1;
        if name_len > 255 {
            return Err(Error::corrupt(
                path,
                slice.offset,
                format!("a read name of {} bytes, which bam cannot hold", name.len()),
            ));
        }
        if record.cigar.len() > u16::MAX as usize {
            return Err(Error::corrupt(
                path,
                slice.offset,
                format!(
                    "a cigar of {} operations, which this reader does not yet spill \
                     into a CG tag",
                    record.cigar.len()
                ),
            ));
        }

        // §10.1: a record may say its sequence is unknown, the read features
        // being there only to rebuild the cigar. BAM spells that `l_seq = 0`.
        let no_sequence = record.cram_flags & cram_flag::NO_SEQUENCE != 0;
        let sequence: &[u8] = if no_sequence { &[] } else { &record.sequence };
        let qualities: &[u8] = if no_sequence { &[] } else { &record.qualities };
        let l_seq = sequence.len();

        let variable =
            name_len + record.cigar.len() * 4 + l_seq.div_ceil(2) + l_seq + record.tags.len();
        let block_size = super::super::bam::record::RECORD_HEADER_SIZE - 4 + variable;

        out.extend_from_slice(&(block_size as u32).to_le_bytes());
        out.extend_from_slice(&record.ref_id.to_le_bytes());
        // CRAM is 1-based and BAM is 0-based; an unplaced read is -1 either way.
        let pos = if record.ref_id < 0 && record.start <= 0 {
            -1i32
        } else {
            (record.start - 1) as i32
        };
        out.extend_from_slice(&pos.to_le_bytes());
        out.push(name_len as u8);
        out.push(record.mapping_quality);
        let bin = if record.ref_id < 0 {
            0
        } else {
            reg2bin(record.start - 1, record.end)
        };
        out.extend_from_slice(&bin.to_le_bytes());
        out.extend_from_slice(&(record.cigar.len() as u16).to_le_bytes());
        out.extend_from_slice(&record.flags.to_le_bytes());
        out.extend_from_slice(&(l_seq as i32).to_le_bytes());
        out.extend_from_slice(&record.next_ref_id.to_le_bytes());
        let next_pos = if record.next_ref_id < 0 && record.next_pos <= 0 {
            -1i32
        } else {
            (record.next_pos - 1) as i32
        };
        out.extend_from_slice(&next_pos.to_le_bytes());
        out.extend_from_slice(&(record.template_length as i32).to_le_bytes());

        out.extend_from_slice(name);
        out.push(0);
        for op in &record.cigar {
            out.extend_from_slice(&op.to_le_bytes());
        }
        // Two bases to a byte, the first in the high nibble.
        for pair in sequence.chunks(2) {
            let high = BASE_CODES[pair[0] as usize];
            let low = pair.get(1).map(|b| BASE_CODES[*b as usize]).unwrap_or(0);
            out.push((high << 4) | low);
        }
        if qualities.len() == l_seq {
            out.extend_from_slice(qualities);
        } else {
            // A file that kept no qualities, or fewer than it has bases.
            out.extend_from_slice(&qualities[..qualities.len().min(l_seq)]);
            out.resize(out.len() + l_seq.saturating_sub(qualities.len()), 0xFF);
        }
        out.extend_from_slice(&record.tags);
    }
    Ok(offset)
}

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

    use crate::fuzz::{cram_compression_header, cram_slice};

    fn header_for_tests() -> CompressionHeader {
        CompressionHeader::parse(&cram_compression_header(), "test").expect("a valid header")
    }

    fn decode(slice: &Slice) -> Result<bytes::Bytes> {
        decode_slice(slice, &header_for_tests(), References::None, &[], "test")
    }

    /// A slice may declare any record count it likes, and a record need not
    /// cost the file anything: every series may be a one-symbol Huffman, which
    /// reads no bits, and a mapped read with no features takes its bases from
    /// the reference. A 364-byte file declaring five million records used to
    /// return five million records and 4.9 GB of resident memory.
    #[test]
    fn a_slice_cannot_declare_more_records_than_its_bytes_could_hold() {
        let mut slice = cram_slice(&[0u8; 16], 1);
        slice.header.n_records = 5_000_000;
        let error = decode(&slice).expect_err("five million records from sixteen bytes");
        assert!(error.to_string().contains("cannot hold"), "{error}");

        // The ceiling holds even when the data is large enough to justify it.
        let mut big = cram_slice(&[0u8; 4096], 1);
        big.header.n_records = (MAX_RECORDS_PER_SLICE + 1) as i32;
        assert!(decode(&big).is_err(), "past the per-slice ceiling");

        // And an ordinary count is still accepted.
        let ordinary = cram_slice(&[0u8; 4096], 8);
        let _ = decode(&ordinary);
    }

    /// `RL` is an `i32` the file chooses and a mapped read with no features
    /// costs nothing per base, so an 813-byte file used to ask for a gibibyte
    /// of sequence and get it.
    #[test]
    fn a_read_longer_than_its_slice_could_hold_is_refused() {
        // ITF8 for 2^30, which every series in this header reads from its own
        // block — including `RL`.
        let huge = [0xf4u8, 0x00, 0x00, 0x00, 0x00];
        let slice = cram_slice(&huge, 1);
        let error = decode(&slice).expect_err("a read of 2^30 bases");
        assert!(error.to_string().contains("a read of"), "{error}");
    }

    /// The sum, not just each one: a million reads of a thousand bases is the
    /// same bomb spread out.
    #[test]
    fn the_read_lengths_of_a_slice_are_bounded_in_total() {
        // ITF8 for 100 000, repeated: each record reads one and the budget for
        // these few bytes is far below their sum.
        let mut data = Vec::new();
        for _ in 0..64 {
            data.extend_from_slice(&[0xe0, 0x01, 0x86, 0xa0]);
        }
        let slice = cram_slice(&data, 64);
        let error = decode(&slice).expect_err("sixty-four hundred-kilobase reads");
        let message = error.to_string();
        assert!(
            message.contains("total more than") || message.contains("cannot hold"),
            "{message}"
        );
    }

    #[test]
    fn reference_bases_outside_the_window_read_as_n() {
        let reference = ReferenceBases {
            bases: b"ACGT",
            start: 101,
        };
        assert_eq!(reference.at(101), b'A');
        assert_eq!(reference.at(104), b'T');
        assert_eq!(reference.at(100), b'N');
        assert_eq!(reference.at(105), b'N');
        // And a reader with no reference at all answers N everywhere.
        let none = ReferenceBases::default();
        assert_eq!(none.at(1), b'N');
    }

    #[test]
    fn reference_bases_are_upper_cased() {
        let reference = ReferenceBases {
            bases: b"acgt",
            start: 1,
        };
        assert_eq!(reference.at(1), b'A');
        assert_eq!(reference.at(3), b'G');
    }

    /// The SAM specification's own `reg2bin`, checked at the boundaries where
    /// the bin level changes.
    #[test]
    fn reg2bin_matches_the_sam_specification() {
        assert_eq!(reg2bin(0, 1), 4681);
        assert_eq!(reg2bin(0, 16384), 4681);
        assert_eq!(reg2bin(16384, 16385), 4682);
        assert_eq!(reg2bin(0, 16385), 585);
        assert_eq!(reg2bin(0, 131073), 73);
        assert_eq!(reg2bin(0, 1 << 20), 73);
        assert_eq!(reg2bin(0, 1 << 26), 1);
        // Past the last level, everything falls into bin 0.
        assert_eq!(reg2bin(0, 1 << 29), 0);
        assert_eq!(reg2bin(1 << 28, (1 << 28) + 1), 21065);
    }

    #[test]
    fn base_codes_cover_the_iupac_alphabet_in_both_cases() {
        for (base, code) in b"=ACMGRSVTWYHKDBN".iter().zip(0u8..) {
            assert_eq!(BASE_CODES[*base as usize], code, "{}", *base as char);
        }
        assert_eq!(BASE_CODES[b'a' as usize], BASE_CODES[b'A' as usize]);
        assert_eq!(BASE_CODES[b'n' as usize], BASE_CODES[b'N' as usize]);
        // Anything else is an N rather than an index off the end of the table.
        assert_eq!(BASE_CODES[b'*' as usize], 15);
        assert_eq!(BASE_CODES[0], 15);
    }
}