fastx/

bgzf.rs

//! BGZF: the block-compressed gzip variant used across the samtools ecosystem.
//!
//! A BGZF file is an ordinary gzip file — any gzip tool can decompress it — made
//! of independent members of at most 64 KiB each, every one carrying its own
//! compressed size in a `BC` extra field. Because the blocks are independent, a
//! reader that knows where they start can seek to any uncompressed position by
//! jumping to the enclosing block and decompressing just that. This is what makes
//! a 3 GB bgzipped reference genome randomly accessible.
//!
//! [`BgzfReader`] implements [`Read`] and [`Seek`] in *uncompressed* coordinates,
//! so anything generic over those two traits — including [`crate::IndexedFasta`] —
//! works on a bgzipped file without knowing it.
//!
//! ```no_run
//! use fastx::bgzf::{BgzfReader, BgzfWriter};
//! use std::io::{Read, Seek, SeekFrom, Write};
//!
//! // Write a seekable gzip file.
//! let mut writer = BgzfWriter::create("ref.fa.gz")?;
//! writer.write_all(b">chr1\nACGT\n")?;
//! writer.finish()?;
//!
//! // Read 4 bytes starting at uncompressed offset 6.
//! let mut reader = BgzfReader::open("ref.fa.gz")?;
//! reader.seek(SeekFrom::Start(6))?;
//! let mut buf = [0u8; 4];
//! reader.read_exact(&mut buf)?;
//! assert_eq!(&buf, b"ACGT");
//! # Ok::<(), fastx::Error>(())
//! ```

use std::fs::File;
use std::io::{self, BufReader, BufWriter, Read, Seek, SeekFrom, Write};
use std::path::{Path, PathBuf};

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

/// Largest uncompressed payload placed in one block.
///
/// The spec caps a whole block at 64 KiB; htslib uses 0xff00 for the payload so
/// that even incompressible data plus headers stays under the limit.
pub const MAX_BLOCK_PAYLOAD: usize = 0xff00;

/// Fixed part of a BGZF gzip header, up to and including `XLEN`.
const HEADER_LEN: usize = 12;
/// The `BC` extra subfield: `SI1`, `SI2`, `SLEN` and `BSIZE`.
const EXTRA_LEN: usize = 6;
/// The gzip trailer: CRC32 and ISIZE.
const TRAILER_LEN: usize = 8;

/// The 28-byte empty block that marks a complete BGZF file.
///
/// Its presence is how tools tell a truncated file from a finished one, so
/// [`BgzfWriter::finish`] always appends it.
pub const EOF_BLOCK: [u8; 28] = [
    0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x02, 0x00,
    0x1b, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
];

/// Header of one BGZF block, as read from the file.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
struct BlockHeader {
    /// Total bytes the block occupies in the file, `BSIZE + 1`.
    compressed_len: usize,
    /// Bytes of the block that precede the deflate payload.
    payload_offset: usize,
}

/// Read and validate one block header from `bytes`.
fn parse_block_header(bytes: &[u8]) -> Result<BlockHeader> {
    if bytes.len() < HEADER_LEN {
        return Err(bgzf_error("block header is truncated"));
    }
    if bytes[0] != 0x1f || bytes[1] != 0x8b {
        return Err(bgzf_error("not a gzip member"));
    }
    if bytes[2] != 8 {
        return Err(bgzf_error("unsupported compression method"));
    }
    if bytes[3] & 0x04 == 0 {
        return Err(bgzf_error(
            "gzip member has no extra field, so it is gzip but not BGZF",
        ));
    }
    let extra_len = u16::from_le_bytes([bytes[10], bytes[11]]) as usize;
    if bytes.len() < HEADER_LEN + extra_len {
        return Err(bgzf_error("extra field is truncated"));
    }
    let extra = &bytes[HEADER_LEN..HEADER_LEN + extra_len];

    // Walk the subfields looking for `BC`; other subfields are legal.
    let mut cursor = 0;
    while cursor + 4 <= extra.len() {
        let si1 = extra[cursor];
        let si2 = extra[cursor + 1];
        let slen = u16::from_le_bytes([extra[cursor + 2], extra[cursor + 3]]) as usize;
        let value = cursor + 4;
        if value + slen > extra.len() {
            return Err(bgzf_error("extra subfield runs past the extra field"));
        }
        if si1 == b'B' && si2 == b'C' {
            if slen != 2 {
                return Err(bgzf_error("BC subfield is not two bytes"));
            }
            let bsize = u16::from_le_bytes([extra[value], extra[value + 1]]) as usize;
            let compressed_len = bsize + 1;
            let overhead = HEADER_LEN + extra_len + TRAILER_LEN;
            if compressed_len <= overhead {
                return Err(bgzf_error("BSIZE is smaller than the block overhead"));
            }
            return Ok(BlockHeader {
                compressed_len,
                payload_offset: HEADER_LEN + extra_len,
            });
        }
        cursor = value + slen;
    }
    Err(bgzf_error("gzip member has no BC extra subfield"))
}

fn bgzf_error(message: &'static str) -> Error {
    Error::Other(format!("BGZF: {message}"))
}

/// Per-thread compression state.
///
/// libdeflate keeps a reusable context that is worth allocating once per worker
/// rather than once per block; flate2 needs no such thing, so without that
/// feature this is just the level.
struct Deflater {
    /// Only flate2 needs the level at deflate time; libdeflate bakes it into the
    /// compressor when that is built.
    #[cfg(not(feature = "libdeflate"))]
    level: CompressionLevel,
    #[cfg(feature = "libdeflate")]
    compressor: libdeflater::Compressor,
}

impl Deflater {
    fn new(level: CompressionLevel) -> Deflater {
        Deflater {
            #[cfg(not(feature = "libdeflate"))]
            level,
            #[cfg(feature = "libdeflate")]
            compressor: libdeflater::Compressor::new(libdeflate_level(level)),
        }
    }

    /// Deflate with no zlib or gzip wrapper, which is what a gzip member holds.
    #[cfg(feature = "libdeflate")]
    fn deflate(&mut self, data: &[u8]) -> io::Result<Vec<u8>> {
        let mut out = vec![0u8; self.compressor.deflate_compress_bound(data.len())];
        let written = self
            .compressor
            .deflate_compress(data, &mut out)
            .map_err(|e| io::Error::other(format!("libdeflate: {e}")))?;
        out.truncate(written);
        Ok(out)
    }

    #[cfg(not(feature = "libdeflate"))]
    fn deflate(&mut self, data: &[u8]) -> io::Result<Vec<u8>> {
        use flate2::write::DeflateEncoder;
        let mut encoder = DeflateEncoder::new(
            Vec::with_capacity(data.len() / 2 + 64),
            flate2::Compression::new(self.level.0.min(9)),
        );
        encoder.write_all(data)?;
        encoder.finish()
    }
}

/// Map this crate's gzip-shaped 0–9 scale onto libdeflate's 0–12.
#[cfg(feature = "libdeflate")]
fn libdeflate_level(level: CompressionLevel) -> libdeflater::CompressionLvl {
    let mapped = match level.0 {
        0 => 0,
        level => ((level.min(9) as i32 - 1) * 11 / 8 + 1).min(12),
    };
    libdeflater::CompressionLvl::new(mapped)
        .unwrap_or_else(|_| libdeflater::CompressionLvl::default())
}

/// Per-thread decompression state, for the same reason as [`Deflater`].
struct Inflater {
    #[cfg(feature = "libdeflate")]
    decompressor: libdeflater::Decompressor,
}

impl Inflater {
    fn new() -> Inflater {
        Inflater {
            #[cfg(feature = "libdeflate")]
            decompressor: libdeflater::Decompressor::new(),
        }
    }

    /// Inflate a raw deflate stream whose uncompressed size is already known.
    #[cfg(feature = "libdeflate")]
    fn inflate(&mut self, payload: &[u8], expected: usize, out: &mut Vec<u8>) -> Result<()> {
        out.clear();
        if expected == 0 {
            // The EOF marker: a valid deflate stream that expands to nothing,
            // and libdeflate will not decompress into an empty buffer.
            return Ok(());
        }
        out.resize(expected, 0);
        let written = self
            .decompressor
            .deflate_decompress(payload, out)
            .map_err(|e| Error::Other(format!("BGZF: corrupt block: {e}")))?;
        if written != expected {
            return Err(bgzf_error("block size does not match its ISIZE field"));
        }
        Ok(())
    }

    #[cfg(not(feature = "libdeflate"))]
    fn inflate(&mut self, payload: &[u8], expected: usize, out: &mut Vec<u8>) -> Result<()> {
        out.clear();
        out.reserve(expected);
        let mut decoder = flate2::Decompress::new(false);
        decoder
            .decompress_vec(payload, out, flate2::FlushDecompress::Finish)
            .map_err(|e| Error::Other(format!("BGZF: corrupt block: {e}")))?;
        if out.len() != expected {
            return Err(bgzf_error("block size does not match its ISIZE field"));
        }
        Ok(())
    }
}

/// Inflate a raw deflate stream of known uncompressed size.
fn inflate(payload: &[u8], expected: usize, out: &mut Vec<u8>) -> Result<()> {
    Inflater::new().inflate(payload, expected, out)
}

/// CRC32 of a block's uncompressed bytes.
///
/// libdeflate's implementation is faster than flate2's, and this runs over every
/// byte written, so it is worth taking when available.
fn crc32(data: &[u8]) -> u32 {
    #[cfg(feature = "libdeflate")]
    {
        let mut crc = libdeflater::Crc::new();
        crc.update(data);
        crc.sum()
    }
    #[cfg(not(feature = "libdeflate"))]
    {
        let mut crc = flate2::Crc::new();
        crc.update(data);
        crc.sum()
    }
}

/// The `.gzi` index that makes a BGZF file seekable by uncompressed offset.
///
/// The on-disk layout is the one `bgzip --index` writes: a little-endian `u64`
/// count followed by that many `(compressed_offset, uncompressed_offset)` pairs.
/// The first block is implicit — it always sits at `(0, 0)` — so a file of *n*
/// blocks yields *n − 1* entries.
#[derive(Debug, Clone, Default, PartialEq, Eq)]
pub struct GziIndex {
    /// Block starts, including the implicit `(0, 0)` first entry.
    blocks: Vec<BlockOffset>,
    /// Total uncompressed size, known only when the index came from scanning a
    /// file or from the writer that produced it. A `.gzi` on disk records block
    /// starts and nothing else, so a parsed index cannot know where the data
    /// ends — and must not guess, or `SeekFrom::End` would land in the middle of
    /// the last block.
    total_uncompressed: Option<u64>,
}

/// Where one block begins, in both coordinate systems.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct BlockOffset {
    /// Byte offset of the block in the compressed file.
    pub compressed: u64,
    /// Byte offset of the block's first byte in the decompressed stream.
    pub uncompressed: u64,
}

impl GziIndex {
    /// Build an index by walking every block header in the file.
    ///
    /// Only the headers are read, not the payloads, so this is fast: it touches
    /// about 18 bytes per 64 KiB of input.
    pub fn build<R: Read + Seek>(mut reader: R) -> Result<GziIndex> {
        reader.seek(SeekFrom::Start(0))?;
        let mut reader = BufReader::with_capacity(64 * 1024, reader);
        let mut blocks = Vec::new();
        let mut compressed = 0u64;
        let mut uncompressed = 0u64;
        let mut header = [0u8; HEADER_LEN + 64];

        loop {
            if !read_exact_or_eof(&mut reader, &mut header[..HEADER_LEN])? {
                break;
            }
            let extra_len = u16::from_le_bytes([header[10], header[11]]) as usize;
            if HEADER_LEN + extra_len > header.len() {
                return Err(bgzf_error("extra field is implausibly large"));
            }
            reader.read_exact(&mut header[HEADER_LEN..HEADER_LEN + extra_len])?;
            let block = parse_block_header(&header[..HEADER_LEN + extra_len])?;

            // Skip the payload and CRC, then read ISIZE.
            let skip = block.compressed_len - block.payload_offset - 4;
            io::copy(&mut reader.by_ref().take(skip as u64), &mut io::sink())?;
            let mut isize_bytes = [0u8; 4];
            reader.read_exact(&mut isize_bytes)?;
            let payload_len = u32::from_le_bytes(isize_bytes) as u64;

            blocks.push(BlockOffset {
                compressed,
                uncompressed,
            });
            compressed += block.compressed_len as u64;
            uncompressed += payload_len;

            // A zero-length block is the EOF marker; nothing follows it.
            if payload_len == 0 {
                break;
            }
        }
        Ok(GziIndex {
            blocks,
            total_uncompressed: Some(uncompressed),
        })
    }

    /// Build an index for a file on disk.
    pub fn build_from_path<P: AsRef<Path>>(path: P) -> Result<GziIndex> {
        let path = path.as_ref();
        GziIndex::build(
            File::open(path).map_err(|e| {
                Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display())))
            })?,
        )
    }

    /// Parse a `.gzi` file.
    pub fn parse<R: Read>(mut reader: R) -> Result<GziIndex> {
        let mut count_bytes = [0u8; 8];
        reader.read_exact(&mut count_bytes)?;
        let count = u64::from_le_bytes(count_bytes);
        // Guard against a corrupt count asking for a terabyte of allocation.
        if count > 1 << 32 {
            return Err(bgzf_error("index claims an implausible number of blocks"));
        }
        // The first block is implicit.
        let mut blocks = Vec::with_capacity(count as usize + 1);
        blocks.push(BlockOffset {
            compressed: 0,
            uncompressed: 0,
        });
        let mut pair = [0u8; 16];
        for _ in 0..count {
            reader.read_exact(&mut pair)?;
            blocks.push(BlockOffset {
                compressed: u64::from_le_bytes(pair[..8].try_into().expect("8 bytes")),
                uncompressed: u64::from_le_bytes(pair[8..].try_into().expect("8 bytes")),
            });
        }
        Ok(GziIndex {
            blocks,
            // A `.gzi` does not record the total size.
            total_uncompressed: None,
        })
    }

    /// Parse a `.gzi` file from disk.
    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<GziIndex> {
        let path = path.as_ref();
        GziIndex::parse(BufReader::new(File::open(path).map_err(|e| {
            Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display())))
        })?))
    }

    /// Serialise in `bgzip --index` format, omitting the implicit first block.
    pub fn write<W: Write>(&self, out: &mut W) -> Result<()> {
        let count = self.blocks.len().saturating_sub(1) as u64;
        out.write_all(&count.to_le_bytes())?;
        for block in self.blocks.iter().skip(1) {
            out.write_all(&block.compressed.to_le_bytes())?;
            out.write_all(&block.uncompressed.to_le_bytes())?;
        }
        Ok(())
    }

    /// Write the index next to the BGZF file as `<path>.gzi`.
    pub fn write_to_path<P: AsRef<Path>>(&self, bgzf_path: P) -> Result<PathBuf> {
        let target = gzi_path(bgzf_path.as_ref());
        let mut file = BufWriter::new(File::create(&target)?);
        self.write(&mut file)?;
        file.flush()?;
        Ok(target)
    }

    /// Block starts, in file order, including the implicit first one.
    pub fn blocks(&self) -> &[BlockOffset] {
        &self.blocks
    }

    /// Number of blocks, including the EOF marker if the file has one.
    pub fn len(&self) -> usize {
        self.blocks.len()
    }

    /// True when the index describes no blocks at all.
    pub fn is_empty(&self) -> bool {
        self.blocks.is_empty()
    }

    /// Total uncompressed size.
    ///
    /// `None` for an index parsed from a `.gzi` file, which records only where
    /// blocks start. Build the index with [`GziIndex::build`] if you need this —
    /// it is a header-only scan, so it is cheap.
    pub fn uncompressed_len(&self) -> Option<u64> {
        self.total_uncompressed
    }

    /// The block containing uncompressed byte `offset`.
    fn block_for(&self, offset: u64) -> Option<BlockOffset> {
        // The last block whose uncompressed start is <= offset.
        match self
            .blocks
            .binary_search_by(|b| b.uncompressed.cmp(&offset))
        {
            Ok(i) => Some(self.blocks[i]),
            Err(0) => None,
            Err(i) => Some(self.blocks[i - 1]),
        }
    }
}

/// The conventional index path for a BGZF file: `<path>.gzi`.
pub fn gzi_path(bgzf: &Path) -> PathBuf {
    let mut name = bgzf.as_os_str().to_os_string();
    name.push(".gzi");
    PathBuf::from(name)
}

/// True when the first bytes look like a BGZF block rather than plain gzip.
///
/// ```
/// # use fastx::bgzf::{is_bgzf, EOF_BLOCK};
/// assert!(is_bgzf(&EOF_BLOCK));
/// assert!(!is_bgzf(&[0x1f, 0x8b, 0x08, 0x00]));  // gzip without FEXTRA
/// assert!(!is_bgzf(b"ACGT"));
/// ```
pub fn is_bgzf(bytes: &[u8]) -> bool {
    parse_block_header(bytes).is_ok()
}

/// A BGZF reader that seeks in uncompressed coordinates.
///
/// Sequential reads need no index. [`Seek`] does: build one with
/// [`GziIndex::build`] (fast, header-only) or load a `.gzi` file.
pub struct BgzfReader<R: Read + Seek> {
    inner: R,
    index: Option<GziIndex>,
    /// Decompressed contents of the block currently in hand.
    block: Vec<u8>,
    /// Read cursor within `block`.
    block_pos: usize,
    /// Uncompressed offset at which `block` begins.
    block_start: u64,
    /// Compressed offset of the next block to read.
    next_compressed: u64,
    eof: bool,
    /// Scratch buffer for one compressed block.
    raw: Vec<u8>,
}

impl BgzfReader<File> {
    /// Open a BGZF file, loading `<path>.gzi` if it is present.
    ///
    /// Without a `.gzi` the reader still works sequentially, and [`Seek`] will
    /// build an index on first use.
    pub fn open<P: AsRef<Path>>(path: P) -> Result<BgzfReader<File>> {
        let path = path.as_ref();
        let file = File::open(path)
            .map_err(|e| Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display()))))?;
        let index = match GziIndex::from_path(gzi_path(path)) {
            Ok(index) => Some(index),
            Err(Error::Io(e)) if e.kind() == io::ErrorKind::NotFound => None,
            Err(e) => return Err(e),
        };
        let mut reader = BgzfReader::new(file)?;
        reader.index = index;
        Ok(reader)
    }
}

impl<R: Read + Seek> BgzfReader<R> {
    /// Wrap a seekable reader, checking that it really is BGZF.
    pub fn new(mut inner: R) -> Result<BgzfReader<R>> {
        inner.seek(SeekFrom::Start(0))?;
        let mut probe = [0u8; HEADER_LEN + 64];
        let read = read_up_to(&mut inner, &mut probe)?;
        if read == 0 {
            return Err(bgzf_error("file is empty"));
        }
        parse_block_header(&probe[..read])?;
        inner.seek(SeekFrom::Start(0))?;
        Ok(BgzfReader {
            inner,
            index: None,
            block: Vec::new(),
            block_pos: 0,
            block_start: 0,
            next_compressed: 0,
            eof: false,
            raw: Vec::new(),
        })
    }

    /// Attach an index, enabling [`Seek`].
    pub fn with_index(mut self, index: GziIndex) -> Self {
        self.index = Some(index);
        self
    }

    /// The index in use, if any.
    pub fn index(&self) -> Option<&GziIndex> {
        self.index.as_ref()
    }

    /// Current uncompressed position.
    pub fn position(&self) -> u64 {
        self.block_start + self.block_pos as u64
    }

    /// Unwrap the underlying reader.
    pub fn into_inner(self) -> R {
        self.inner
    }

    /// Decompress the block at compressed offset `at`, which becomes current.
    fn load_block_at(&mut self, at: u64, uncompressed_start: u64) -> Result<()> {
        self.inner.seek(SeekFrom::Start(at))?;
        self.next_compressed = at;
        self.block_start = uncompressed_start;
        self.block_pos = 0;
        self.block.clear();
        self.eof = false;
        self.read_next_block()
    }

    /// Read and decompress the block at `self.next_compressed`.
    fn read_next_block(&mut self) -> Result<()> {
        let mut header = [0u8; HEADER_LEN + 64];
        if !read_exact_or_eof(&mut self.inner, &mut header[..HEADER_LEN])? {
            self.eof = true;
            self.block.clear();
            self.block_pos = 0;
            return Ok(());
        }
        let extra_len = u16::from_le_bytes([header[10], header[11]]) as usize;
        if HEADER_LEN + extra_len > header.len() {
            return Err(bgzf_error("extra field is implausibly large"));
        }
        self.inner
            .read_exact(&mut header[HEADER_LEN..HEADER_LEN + extra_len])?;
        let block = parse_block_header(&header[..HEADER_LEN + extra_len])?;

        let payload_len = block.compressed_len - block.payload_offset - TRAILER_LEN;
        self.raw.resize(payload_len, 0);
        self.inner.read_exact(&mut self.raw)?;
        let mut trailer = [0u8; TRAILER_LEN];
        self.inner.read_exact(&mut trailer)?;
        let expected_crc = u32::from_le_bytes(trailer[..4].try_into().expect("4 bytes"));
        let expected_len = u32::from_le_bytes(trailer[4..].try_into().expect("4 bytes")) as usize;

        let mut decompressed = std::mem::take(&mut self.block);
        let result = inflate(&self.raw, expected_len, &mut decompressed);
        self.block = decompressed;
        result?;

        if crc32(&self.block) != expected_crc {
            return Err(bgzf_error("block CRC32 does not match"));
        }

        self.block_pos = 0;
        self.next_compressed += block.compressed_len as u64;
        // The EOF marker decompresses to nothing; treat it as end of stream.
        if self.block.is_empty() {
            self.eof = true;
        }
        Ok(())
    }

    /// Make sure the current block has unread bytes, or set `eof`.
    fn fill(&mut self) -> Result<()> {
        while !self.eof && self.block_pos == self.block.len() {
            let consumed = self.block.len() as u64;
            self.block_start += consumed;
            self.read_next_block()?;
        }
        Ok(())
    }

    /// The index, built on demand if it was not supplied.
    fn ensure_index(&mut self) -> Result<()> {
        if self.index.is_none() {
            let saved = self.inner.stream_position()?;
            let index = GziIndex::build(&mut self.inner)?;
            self.inner.seek(SeekFrom::Start(saved))?;
            self.index = Some(index);
        }
        Ok(())
    }
}

impl<R: Read + Seek> Read for BgzfReader<R> {
    fn read(&mut self, buf: &mut [u8]) -> io::Result<usize> {
        if buf.is_empty() {
            return Ok(0);
        }
        self.fill()?;
        if self.eof && self.block_pos == self.block.len() {
            return Ok(0);
        }
        let available = &self.block[self.block_pos..];
        let take = available.len().min(buf.len());
        buf[..take].copy_from_slice(&available[..take]);
        self.block_pos += take;
        Ok(take)
    }
}

impl<R: Read + Seek> Seek for BgzfReader<R> {
    /// Seek in *uncompressed* coordinates.
    ///
    /// `SeekFrom::End` needs the total uncompressed size, so it requires an
    /// index that covers the whole file.
    fn seek(&mut self, pos: SeekFrom) -> io::Result<u64> {
        self.ensure_index()?;

        let target = match pos {
            SeekFrom::Start(offset) => offset,
            SeekFrom::Current(delta) => add_signed(self.position(), delta)?,
            SeekFrom::End(delta) => {
                let end = self
                    .index
                    .as_ref()
                    .and_then(|index| index.uncompressed_len())
                    .ok_or_else(|| {
                        io::Error::new(
                            io::ErrorKind::InvalidInput,
                            "BGZF: this index records block starts only, so the end of the \
                             stream is unknown; rebuild it with GziIndex::build",
                        )
                    })?;
                add_signed(end, delta)?
            }
        };

        // Staying inside the current block is the common case for short hops.
        let within = target.checked_sub(self.block_start);
        if let Some(within) = within {
            if !self.block.is_empty() && within <= self.block.len() as u64 {
                self.block_pos = within as usize;
                return Ok(target);
            }
        }

        // `BlockOffset` is `Copy`, so the borrow of the index ends here and the
        // load below needs no clone of it.
        let block = self
            .index
            .as_ref()
            .and_then(|index| index.block_for(target))
            .ok_or_else(|| {
                io::Error::new(
                    io::ErrorKind::InvalidInput,
                    "BGZF: no block covers that offset",
                )
            })?;
        self.load_block_at(block.compressed, block.uncompressed)?;
        let within = (target - block.uncompressed) as usize;
        if within > self.block.len() {
            // Past the end of the data: leave the cursor at the end.
            self.block_pos = self.block.len();
            self.eof = true;
        } else {
            self.block_pos = within;
        }
        Ok(target)
    }
}

fn add_signed(base: u64, delta: i64) -> io::Result<u64> {
    let result = if delta >= 0 {
        base.checked_add(delta as u64)
    } else {
        base.checked_sub(delta.unsigned_abs())
    };
    result.ok_or_else(|| {
        io::Error::new(
            io::ErrorKind::InvalidInput,
            "BGZF: seek would leave the file",
        )
    })
}

/// One block as it sits in the file: still compressed, with its expected
/// checksum and length.
///
/// Only the parallel reader needs this: the seekable one decompresses straight
/// into its own buffer.
#[cfg(feature = "parallel")]
struct RawBlock {
    payload: Vec<u8>,
    crc: u32,
    uncompressed_len: usize,
}

/// Read one block without decompressing it. `None` at a clean end of file.
#[cfg(feature = "parallel")]
fn read_raw_block<R: Read>(reader: &mut R) -> Result<Option<RawBlock>> {
    let mut header = [0u8; HEADER_LEN + 64];
    if !read_exact_or_eof(reader, &mut header[..HEADER_LEN])? {
        return Ok(None);
    }
    let extra_len = u16::from_le_bytes([header[10], header[11]]) as usize;
    if HEADER_LEN + extra_len > header.len() {
        return Err(bgzf_error("extra field is implausibly large"));
    }
    reader.read_exact(&mut header[HEADER_LEN..HEADER_LEN + extra_len])?;
    let block = parse_block_header(&header[..HEADER_LEN + extra_len])?;

    let payload_len = block.compressed_len - block.payload_offset - TRAILER_LEN;
    let mut payload = vec![0u8; payload_len];
    reader.read_exact(&mut payload)?;
    let mut trailer = [0u8; TRAILER_LEN];
    reader.read_exact(&mut trailer)?;

    Ok(Some(RawBlock {
        payload,
        crc: u32::from_le_bytes(trailer[..4].try_into().expect("4 bytes")),
        uncompressed_len: u32::from_le_bytes(trailer[4..].try_into().expect("4 bytes")) as usize,
    }))
}

/// Decompress one raw block and check it against its trailer.
#[cfg(feature = "parallel")]
fn inflate_checked(inflater: &mut Inflater, raw: &RawBlock) -> Result<Vec<u8>> {
    let mut out = Vec::new();
    inflater.inflate(&raw.payload, raw.uncompressed_len, &mut out)?;
    if crc32(&out) != raw.crc {
        return Err(bgzf_error("block CRC32 does not match"));
    }
    Ok(out)
}

/// A BGZF reader that decompresses blocks across cores.
///
/// Blocks are independent, so inflating them is embarrassingly parallel: one
/// thread reads compressed blocks — which is cheap, mostly `memcpy` — and a
/// batch is then inflated on the rayon pool and served in order.
///
/// This trades seeking for throughput. Use it for a full pass over a compressed
/// file; use [`BgzfReader`] when you need [`Seek`], which by its nature has to
/// decompress one block at a time.
///
/// ```no_run
/// use fastx::bgzf::ParallelBgzfReader;
/// use std::io::Read;
///
/// let file = std::fs::File::open("reads.fq.gz")?;
/// let mut reader = ParallelBgzfReader::new(file);
/// let mut all = Vec::new();
/// reader.read_to_end(&mut all)?;
/// # Ok::<(), fastx::Error>(())
/// ```
#[cfg(feature = "parallel")]
pub struct ParallelBgzfReader<R: Read> {
    inner: R,
    batch: usize,
    /// The current batch, decompressed and concatenated.
    buffer: Vec<u8>,
    pos: usize,
    eof: bool,
}

#[cfg(feature = "parallel")]
impl<R: Read> ParallelBgzfReader<R> {
    /// Wrap a reader, inflating one batch of blocks per core at a time.
    pub fn new(inner: R) -> ParallelBgzfReader<R> {
        ParallelBgzfReader {
            inner,
            batch: default_batch().max(2),
            buffer: Vec::new(),
            pos: 0,
            eof: false,
        }
    }

    /// How many blocks to inflate at a time. Larger batches balance load better
    /// and cost `blocks × 64 KiB` of memory.
    pub fn blocks_per_batch(mut self, blocks: usize) -> Self {
        self.batch = blocks.max(1);
        self
    }

    /// Unwrap the underlying reader.
    pub fn into_inner(self) -> R {
        self.inner
    }

    /// Read and inflate the next batch.
    fn fill(&mut self) -> Result<()> {
        use rayon::prelude::*;

        self.buffer.clear();
        self.pos = 0;

        let mut raws: Vec<RawBlock> = Vec::with_capacity(self.batch);
        while raws.len() < self.batch {
            match read_raw_block(&mut self.inner)? {
                None => {
                    self.eof = true;
                    break;
                }
                // A zero-length block is the EOF marker; nothing follows it.
                Some(raw) if raw.uncompressed_len == 0 => {
                    self.eof = true;
                    break;
                }
                Some(raw) => raws.push(raw),
            }
        }
        if raws.is_empty() {
            return Ok(());
        }

        let blocks: Vec<Vec<u8>> = if raws.len() > 1 {
            // One decompressor per worker, not per block: see `encode_batch`.
            raws.par_iter()
                .map_init(Inflater::new, inflate_checked)
                .collect::<Result<_>>()?
        } else {
            vec![inflate_checked(&mut Inflater::new(), &raws[0])?]
        };
        self.buffer.reserve(blocks.iter().map(Vec::len).sum());
        for block in &blocks {
            self.buffer.extend_from_slice(block);
        }
        Ok(())
    }
}

#[cfg(feature = "parallel")]
impl<R: Read> Read for ParallelBgzfReader<R> {
    fn read(&mut self, buf: &mut [u8]) -> io::Result<usize> {
        if buf.is_empty() {
            return Ok(0);
        }
        while self.pos == self.buffer.len() {
            if self.eof {
                return Ok(0);
            }
            self.fill()?;
        }
        let available = &self.buffer[self.pos..];
        let take = available.len().min(buf.len());
        buf[..take].copy_from_slice(&available[..take]);
        self.pos += take;
        Ok(take)
    }
}

/// A writer that emits BGZF blocks.
///
/// Output is valid gzip, so any gzip tool can read it, and it is seekable by
/// anything that understands BGZF. Call [`BgzfWriter::finish`] to append the
/// EOF marker; without it the file looks truncated to samtools.
pub struct BgzfWriter<W: Write> {
    /// `None` once `finish` has handed the writer back. Wrapped in an `Option`
    /// only because a type with a `Drop` impl cannot give a field away.
    inner: Option<W>,
    /// The payload currently being filled.
    buffer: Vec<u8>,
    /// Full payloads waiting to be compressed as a batch.
    pending: Vec<Vec<u8>>,
    /// Payload buffers to reuse, so a long run does not allocate per block.
    recycled: Vec<Vec<u8>>,
    /// How many blocks to compress at once. Deflating one block does not depend
    /// on any other, so a batch fans out across cores; the output bytes are
    /// identical either way, which is why this needs no opt-in.
    batch: usize,
    level: CompressionLevel,
    /// Block starts, recorded so that an index can be written afterwards.
    blocks: Vec<BlockOffset>,
    compressed: u64,
    uncompressed: u64,
    finished: bool,
}

/// One block, compressed and framed, ready to be handed to the sink.
///
/// Produced by [`encode_block`], which touches no shared state and can therefore
/// run on a worker thread.
struct EncodedBlock {
    header: [u8; HEADER_LEN + EXTRA_LEN],
    payload: Vec<u8>,
    trailer: [u8; TRAILER_LEN],
    block_len: usize,
    uncompressed_len: usize,
}

/// Compress one payload and build its gzip framing.
fn encode_block(deflater: &mut Deflater, data: &[u8]) -> io::Result<EncodedBlock> {
    let payload = deflater.deflate(data)?;
    let block_len = HEADER_LEN + EXTRA_LEN + payload.len() + TRAILER_LEN;
    if block_len > u16::MAX as usize + 1 {
        // Cannot happen with MAX_BLOCK_PAYLOAD, but a wrong constant here would
        // produce files other tools silently misread.
        return Err(io::Error::new(
            io::ErrorKind::InvalidData,
            "BGZF: block would exceed 64 KiB",
        ));
    }

    let mut header = [0u8; HEADER_LEN + EXTRA_LEN];
    header[0] = 0x1f;
    header[1] = 0x8b;
    header[2] = 8; // deflate
    header[3] = 4; // FEXTRA
                   // MTIME stays zero: reproducible output matters more than a timestamp.
    header[9] = 0xff; // unknown OS
    header[10..12].copy_from_slice(&(EXTRA_LEN as u16).to_le_bytes());
    header[12] = b'B';
    header[13] = b'C';
    header[14..16].copy_from_slice(&2u16.to_le_bytes());
    header[16..18].copy_from_slice(&((block_len - 1) as u16).to_le_bytes());

    let mut trailer = [0u8; TRAILER_LEN];
    trailer[..4].copy_from_slice(&crc32(data).to_le_bytes());
    trailer[4..].copy_from_slice(&(data.len() as u32).to_le_bytes());

    Ok(EncodedBlock {
        header,
        payload,
        trailer,
        block_len,
        uncompressed_len: data.len(),
    })
}

/// Compress a batch of payloads, across cores when the `parallel` feature is on.
///
/// Order is preserved, and the result is identical to compressing them one by
/// one — deflating a BGZF block depends on nothing outside that block.
fn encode_batch(payloads: &[Vec<u8>], level: CompressionLevel) -> io::Result<Vec<EncodedBlock>> {
    #[cfg(feature = "parallel")]
    {
        if payloads.len() > 1 {
            use rayon::prelude::*;
            // `map_init` builds one compressor per worker rather than one per
            // block, which matters for libdeflate: its context is a real
            // allocation, and there are thousands of blocks in a large file.
            return payloads
                .par_iter()
                .map_init(
                    || Deflater::new(level),
                    |deflater, payload| encode_block(deflater, payload),
                )
                .collect();
        }
    }
    let mut deflater = Deflater::new(level);
    payloads
        .iter()
        .map(|payload| encode_block(&mut deflater, payload))
        .collect()
}

/// Blocks per core in a default batch.
///
/// One block per core is the obvious choice and measurably the wrong one: a
/// batch of 8 is only 512 KiB, so the fan-out happens hundreds of times over a
/// large file and every batch ends with idle cores waiting for its slowest
/// block. Several blocks per core amortise that at a cost of
/// `cores × this × 64 KiB` of memory.
#[cfg(feature = "parallel")]
const BATCHES_PER_CORE: usize = 8;

/// Default number of blocks compressed or decompressed together.
fn default_batch() -> usize {
    #[cfg(feature = "parallel")]
    {
        std::thread::available_parallelism().map_or(1, |n| n.get() * BATCHES_PER_CORE)
    }
    #[cfg(not(feature = "parallel"))]
    {
        1
    }
}

impl BgzfWriter<BufWriter<File>> {
    /// Create a BGZF file.
    pub fn create<P: AsRef<Path>>(path: P) -> Result<BgzfWriter<BufWriter<File>>> {
        let path = path.as_ref();
        let file = File::create(path)
            .map_err(|e| Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display()))))?;
        Ok(BgzfWriter::new(BufWriter::with_capacity(128 * 1024, file)))
    }
}

impl<W: Write> BgzfWriter<W> {
    /// Wrap a writer, compressing at the default level.
    pub fn new(inner: W) -> BgzfWriter<W> {
        BgzfWriter::with_level(inner, CompressionLevel::default())
    }

    /// Wrap a writer, compressing at `level`.
    pub fn with_level(inner: W, level: CompressionLevel) -> BgzfWriter<W> {
        BgzfWriter {
            inner: Some(inner),
            buffer: Vec::with_capacity(MAX_BLOCK_PAYLOAD),
            pending: Vec::new(),
            recycled: Vec::new(),
            batch: default_batch(),
            level,
            blocks: vec![BlockOffset {
                compressed: 0,
                uncompressed: 0,
            }],
            compressed: 0,
            uncompressed: 0,
            finished: false,
        }
    }

    /// Borrow the sink, or `None` once `finish` has taken it.
    pub fn get_ref(&self) -> Option<&W> {
        self.inner.as_ref()
    }

    /// How many blocks to compress at a time.
    ///
    /// Defaults to the number of available cores with the `parallel` feature and
    /// to 1 without it. Compressing a batch produces byte-for-byte the same file
    /// as compressing one at a time, since blocks are independent — this only
    /// trades memory (`blocks × 64 KiB`) for cores.
    ///
    /// Setting 1 forces the single-threaded path, which is what you want if the
    /// surrounding program is already saturating every core itself.
    pub fn blocks_per_batch(mut self, blocks: usize) -> Self {
        self.batch = blocks.max(1);
        self
    }

    /// The index describing the blocks written *so far*.
    ///
    /// Buffered data that has not been compressed yet is not in it, so calling
    /// this before [`BgzfWriter::finish_with_index`] gives an index that stops
    /// short of the end of the file. Prefer `finish_with_index`, which cannot be
    /// wrong.
    pub fn index(&self) -> GziIndex {
        GziIndex {
            blocks: self.blocks.clone(),
            total_uncompressed: Some(self.uncompressed),
        }
    }

    /// The underlying writer, or an error once `finish` has taken it.
    fn sink(&mut self) -> io::Result<&mut W> {
        self.inner.as_mut().ok_or_else(|| {
            io::Error::new(
                io::ErrorKind::BrokenPipe,
                "BGZF: writer was already finished",
            )
        })
    }

    /// Compress and emit whatever is buffered as one block.
    /// Move the payload being filled into the pending batch, compressing the
    /// batch once it is full.
    fn seal_buffer(&mut self) -> io::Result<()> {
        if self.buffer.is_empty() {
            return Ok(());
        }
        let mut fresh = self.recycled.pop().unwrap_or_else(|| {
            let mut buffer = Vec::new();
            buffer.reserve_exact(MAX_BLOCK_PAYLOAD);
            buffer
        });
        std::mem::swap(&mut self.buffer, &mut fresh);
        self.pending.push(fresh);
        if self.pending.len() >= self.batch {
            self.compress_pending()?;
        }
        Ok(())
    }

    /// Compress every pending payload and write the blocks out in order.
    fn compress_pending(&mut self) -> io::Result<()> {
        if self.pending.is_empty() {
            return Ok(());
        }
        let level = self.level;
        let encoded = encode_batch(&self.pending, level)?;
        for block in encoded {
            self.emit(block)?;
        }
        // Hand the payload buffers back for reuse rather than dropping them.
        for mut buffer in self.pending.drain(..) {
            buffer.clear();
            if self.recycled.len() < self.batch {
                self.recycled.push(buffer);
            }
        }
        Ok(())
    }

    /// Write one already-compressed block and record where it landed.
    fn emit(&mut self, block: EncodedBlock) -> io::Result<()> {
        let sink = self.sink()?;
        sink.write_all(&block.header)?;
        sink.write_all(&block.payload)?;
        sink.write_all(&block.trailer)?;

        self.compressed += block.block_len as u64;
        self.uncompressed += block.uncompressed_len as u64;
        self.blocks.push(BlockOffset {
            compressed: self.compressed,
            uncompressed: self.uncompressed,
        });
        Ok(())
    }

    /// Seal whatever is buffered and compress everything outstanding.
    fn flush_block(&mut self) -> io::Result<()> {
        self.seal_buffer()?;
        self.compress_pending()
    }

    /// Flush the pending block, append the EOF marker and hand the writer back.
    ///
    /// Dropping the writer does the same on a best-effort basis, but only
    /// `finish` reports a failure.
    pub fn finish(self) -> Result<W> {
        self.finish_with_index().map(|(inner, _)| inner)
    }

    /// Finish, and return the complete index alongside the writer.
    ///
    /// This is the safe way to obtain a `.gzi`: every block has been emitted by
    /// the time the index is taken, so it cannot be missing the tail.
    ///
    /// ```no_run
    /// use fastx::bgzf::BgzfWriter;
    /// use std::io::Write;
    ///
    /// let mut writer = BgzfWriter::create("reads.fq.gz")?;
    /// writer.write_all(b"@r\nACGT\n+\nIIII\n")?;
    /// let (_file, index) = writer.finish_with_index()?;
    /// index.write_to_path("reads.fq.gz")?;   // reads.fq.gz.gzi
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn finish_with_index(mut self) -> Result<(W, GziIndex)> {
        self.finish_in_place()?;
        let index = self.index();
        let inner = self
            .inner
            .take()
            .ok_or_else(|| Error::Other("BGZF: writer was already finished".to_string()))?;
        Ok((inner, index))
    }

    fn finish_in_place(&mut self) -> Result<()> {
        if self.finished || self.inner.is_none() {
            return Ok(());
        }
        self.flush_block()?;
        let sink = self.sink()?;
        sink.write_all(&EOF_BLOCK)?;
        sink.flush()?;
        self.compressed += EOF_BLOCK.len() as u64;
        self.finished = true;
        Ok(())
    }
}

impl<W: Write> Write for BgzfWriter<W> {
    fn write(&mut self, buf: &[u8]) -> io::Result<usize> {
        let room = MAX_BLOCK_PAYLOAD - self.buffer.len();
        let take = room.min(buf.len());
        self.buffer.extend_from_slice(&buf[..take]);
        if self.buffer.len() == MAX_BLOCK_PAYLOAD {
            // Only seal it: compression waits until a whole batch has gathered.
            self.seal_buffer()?;
        }
        Ok(take)
    }

    /// Ends the current block, so the next byte starts a fresh one.
    ///
    /// This is what makes a flush point seekable, and it costs a little
    /// compression, so flush at record boundaries rather than per record.
    fn flush(&mut self) -> io::Result<()> {
        self.flush_block()?;
        self.sink()?.flush()
    }
}

impl<W: Write> Drop for BgzfWriter<W> {
    fn drop(&mut self) {
        // Best effort: a caller who wants to see errors uses finish().
        let _ = self.finish_in_place();
    }
}

/// Read into `buf` fully, or report `false` if the reader was already at EOF.
fn read_exact_or_eof<R: Read>(reader: &mut R, buf: &mut [u8]) -> Result<bool> {
    let mut filled = 0;
    while filled < buf.len() {
        match reader.read(&mut buf[filled..]) {
            Ok(0) if filled == 0 => return Ok(false),
            Ok(0) => return Err(bgzf_error("file ends in the middle of a block")),
            Ok(n) => filled += n,
            Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
            Err(e) => return Err(Error::Io(e)),
        }
    }
    Ok(true)
}

/// Read as much as is available, up to `buf.len()`.
fn read_up_to<R: Read>(reader: &mut R, buf: &mut [u8]) -> Result<usize> {
    let mut filled = 0;
    while filled < buf.len() {
        match reader.read(&mut buf[filled..]) {
            Ok(0) => break,
            Ok(n) => filled += n,
            Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
            Err(e) => return Err(Error::Io(e)),
        }
    }
    Ok(filled)
}

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

    fn compress(data: &[u8]) -> Vec<u8> {
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(data).unwrap();
        writer.finish().unwrap()
    }

    #[test]
    fn round_trips_through_our_own_reader() {
        for size in [
            0usize,
            1,
            100,
            MAX_BLOCK_PAYLOAD - 1,
            MAX_BLOCK_PAYLOAD,
            MAX_BLOCK_PAYLOAD + 1,
            300_000,
        ] {
            let data: Vec<u8> = (0..size).map(|i| b"ACGTN"[i % 5]).collect();
            let compressed = compress(&data);
            assert!(is_bgzf(&compressed), "size {size} did not produce BGZF");

            let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
            let mut out = Vec::new();
            reader.read_to_end(&mut out).unwrap();
            assert_eq!(out, data, "size {size}");
        }
    }

    #[test]
    fn output_is_plain_gzip_too() {
        // The whole point of BGZF: ordinary gzip tools must still read it.
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let mut out = Vec::new();
        flate2::read::MultiGzDecoder::new(&compressed[..])
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, data);
    }

    #[test]
    fn ends_with_the_eof_marker() {
        let compressed = compress(b"ACGT");
        assert_eq!(
            &compressed[compressed.len() - EOF_BLOCK.len()..],
            &EOF_BLOCK
        );
        // An empty file is just the marker.
        assert_eq!(compress(b""), EOF_BLOCK.to_vec());
    }

    #[test]
    fn blocks_stay_within_the_size_limit() {
        // Incompressible data is the case that could overflow a block.
        let data: Vec<u8> = (0..500_000)
            .map(|i| ((i * 2_654_435_761u64 as usize) >> 7) as u8)
            .collect();
        let compressed = compress(&data);
        let index = GziIndex::build(Cursor::new(&compressed)).unwrap();
        for pair in index.blocks().windows(2) {
            let block_len = pair[1].compressed - pair[0].compressed;
            assert!(block_len <= 65_536, "block of {block_len} bytes");
        }
        let mut out = Vec::new();
        BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, data);
    }

    #[test]
    fn index_from_the_writer_matches_a_rescan() {
        let data: Vec<u8> = (0..250_000).map(|i| b"ACGTN"[i % 5]).collect();
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(&data).unwrap();
        let from_writer = writer.index();
        let compressed = writer.finish().unwrap();

        let from_scan = GziIndex::build(Cursor::new(&compressed)).unwrap();
        // The rescan also sees the EOF block, which the writer's index predates.
        assert_eq!(
            &from_scan.blocks()[..from_writer.len()],
            from_writer.blocks()
        );
        assert_eq!(from_scan.uncompressed_len(), Some(data.len() as u64));
    }

    #[test]
    fn gzi_round_trips_and_omits_the_first_block() {
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let index = GziIndex::build(Cursor::new(&compressed)).unwrap();

        let mut text = Vec::new();
        index.write(&mut text).unwrap();
        // 8-byte count plus 16 bytes per entry, first block implicit.
        assert_eq!(text.len(), 8 + 16 * (index.len() - 1));
        assert_eq!(
            u64::from_le_bytes(text[..8].try_into().unwrap()),
            index.len() as u64 - 1
        );

        let reparsed = GziIndex::parse(&text[..]).unwrap();
        assert_eq!(reparsed.blocks(), index.blocks());
        // A parsed index cannot know where the data ends, and must not pretend.
        assert_eq!(reparsed.uncompressed_len(), None);
        assert_eq!(index.uncompressed_len(), Some(200_000));
    }

    #[test]
    fn seek_from_end_needs_a_scanned_index() {
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let scanned = GziIndex::build(Cursor::new(&compressed)).unwrap();
        let mut text = Vec::new();
        scanned.write(&mut text).unwrap();
        let parsed = GziIndex::parse(&text[..]).unwrap();

        // With block starts only, End-relative seeks must fail loudly rather
        // than landing at the start of the last block.
        let mut reader = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .with_index(parsed);
        assert!(reader.seek(SeekFrom::End(0)).is_err());
        // Absolute seeks still work.
        reader.seek(SeekFrom::Start(199_998)).unwrap();
        let mut tail = Vec::new();
        reader.read_to_end(&mut tail).unwrap();
        assert_eq!(tail, &data[199_998..]);

        let mut reader = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .with_index(scanned);
        assert_eq!(reader.seek(SeekFrom::End(0)).unwrap(), 200_000);
    }

    #[test]
    fn seeks_to_any_offset() {
        let data: Vec<u8> = (0..300_000).map(|i| (i % 251) as u8).collect();
        let compressed = compress(&data);
        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();

        // Offsets that land in the first block, deep inside, and on boundaries.
        for target in [
            0usize,
            1,
            MAX_BLOCK_PAYLOAD - 1,
            MAX_BLOCK_PAYLOAD,
            MAX_BLOCK_PAYLOAD + 1,
            2 * MAX_BLOCK_PAYLOAD,
            299_999,
        ] {
            reader.seek(SeekFrom::Start(target as u64)).unwrap();
            assert_eq!(reader.position(), target as u64);
            let mut buf = [0u8; 8];
            let want = (data.len() - target).min(buf.len());
            reader.read_exact(&mut buf[..want]).unwrap();
            assert_eq!(&buf[..want], &data[target..target + want], "at {target}");
        }

        // Relative and end-relative seeks.
        reader.seek(SeekFrom::Start(10)).unwrap();
        reader.seek(SeekFrom::Current(5)).unwrap();
        assert_eq!(reader.position(), 15);
        assert_eq!(reader.seek(SeekFrom::End(0)).unwrap(), data.len() as u64);
        let mut rest = Vec::new();
        reader.read_to_end(&mut rest).unwrap();
        assert!(rest.is_empty());

        // Seeking backwards must work as well as forwards.
        reader.seek(SeekFrom::Start(7)).unwrap();
        let mut buf = [0u8; 4];
        reader.read_exact(&mut buf).unwrap();
        assert_eq!(&buf, &data[7..11]);
    }

    #[test]
    fn seek_before_the_start_is_an_error() {
        let compressed = compress(b"ACGT");
        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
        assert!(reader.seek(SeekFrom::Current(-1)).is_err());
        assert!(reader.seek(SeekFrom::End(-100)).is_err());
    }

    #[test]
    fn rejects_plain_gzip_and_garbage() {
        // Valid gzip, but no BC extra field: not BGZF.
        let mut plain = Vec::new();
        {
            let mut encoder =
                flate2::write::GzEncoder::new(&mut plain, flate2::Compression::default());
            encoder.write_all(b"ACGT").unwrap();
            encoder.finish().unwrap();
        }
        assert!(!is_bgzf(&plain));
        assert!(BgzfReader::new(Cursor::new(&plain)).is_err());

        assert!(BgzfReader::new(Cursor::new(b"not gzip at all".to_vec())).is_err());
        assert!(BgzfReader::new(Cursor::new(Vec::new())).is_err());
    }

    #[test]
    fn detects_a_corrupt_block() {
        let mut compressed = compress(&vec![b'A'; 5_000]);
        // Flip a byte in the deflate payload; the CRC or the inflate must catch it.
        let victim = HEADER_LEN + EXTRA_LEN + 5;
        compressed[victim] ^= 0xff;
        let mut out = Vec::new();
        let result = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .read_to_end(&mut out);
        assert!(result.is_err(), "corruption went unnoticed");
    }

    #[test]
    fn truncated_file_is_an_error_not_silent_truncation() {
        let compressed = compress(&vec![b'A'; 200_000]);
        let cut = compressed.len() / 2;
        let mut out = Vec::new();
        let result = BgzfReader::new(Cursor::new(compressed[..cut].to_vec()))
            .unwrap()
            .read_to_end(&mut out);
        assert!(result.is_err(), "truncation went unnoticed");
    }

    #[test]
    fn parse_rejects_an_implausible_index() {
        let mut bad = u64::MAX.to_le_bytes().to_vec();
        bad.extend_from_slice(&[0u8; 16]);
        assert!(GziIndex::parse(&bad[..]).is_err());
        assert!(GziIndex::parse(&[0u8; 3][..]).is_err());
    }

    #[test]
    fn batch_size_never_changes_the_bytes() {
        // Deflating a BGZF block depends on nothing outside that block, so the
        // file must be byte-identical however the work is divided. If this ever
        // fails, parallel output has stopped being reproducible.
        let data: Vec<u8> = (0..400_000).map(|i| b"ACGTN"[i % 5]).collect();
        let mut reference = None;
        for blocks in [1usize, 2, 3, 7, 64, 1024] {
            let mut writer = BgzfWriter::new(Vec::new()).blocks_per_batch(blocks);
            writer.write_all(&data).unwrap();
            let (bytes, index) = writer.finish_with_index().unwrap();

            match &reference {
                None => reference = Some((bytes, index)),
                Some((expected_bytes, expected_index)) => {
                    assert_eq!(&bytes, expected_bytes, "batch of {blocks} differs");
                    assert_eq!(index.blocks(), expected_index.blocks(), "index differs");
                }
            }
        }
        // And the bytes still decompress to the input.
        let (bytes, _) = reference.unwrap();
        let mut out = Vec::new();
        BgzfReader::new(Cursor::new(&bytes))
            .unwrap()
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, data);
    }

    #[test]
    fn write_flush_boundaries_survive_batching() {
        // A flush mid-batch has to seal the current block *and* drain everything
        // pending, or the flushed bytes would sit in memory unwritten.
        let mut writer = BgzfWriter::new(Vec::new()).blocks_per_batch(16);
        writer.write_all(b"first").unwrap();
        writer.flush().unwrap();
        assert!(
            !writer.get_ref().expect("still open").is_empty(),
            "flush left the batch uncompressed"
        );
        writer.write_all(b"second").unwrap();
        let (bytes, index) = writer.finish_with_index().unwrap();

        assert_eq!(index.len(), 3); // implicit start plus two blocks
        let mut out = Vec::new();
        BgzfReader::new(Cursor::new(&bytes))
            .unwrap()
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, b"firstsecond");
    }

    #[cfg(feature = "parallel")]
    #[test]
    fn parallel_reader_agrees_with_the_serial_one() {
        for size in [
            0usize,
            1,
            MAX_BLOCK_PAYLOAD - 1,
            MAX_BLOCK_PAYLOAD,
            MAX_BLOCK_PAYLOAD + 1,
            500_000,
        ] {
            let data: Vec<u8> = (0..size).map(|i| (i % 251) as u8).collect();
            let compressed = compress(&data);

            for blocks in [1usize, 2, 5, 64] {
                let mut parallel =
                    ParallelBgzfReader::new(Cursor::new(&compressed)).blocks_per_batch(blocks);
                let mut out = Vec::new();
                parallel.read_to_end(&mut out).unwrap();
                assert_eq!(out, data, "size {size}, batch {blocks}");
            }

            // Byte-at-a-time reads must work too, since `Read` allows any size.
            let mut parallel = ParallelBgzfReader::new(Cursor::new(&compressed));
            let mut out = Vec::new();
            let mut byte = [0u8; 1];
            while parallel.read(&mut byte).unwrap() == 1 {
                out.push(byte[0]);
            }
            assert_eq!(out, data, "size {size}, byte at a time");
        }
    }

    #[cfg(feature = "parallel")]
    #[test]
    fn parallel_reader_still_catches_corruption() {
        let mut compressed = compress(&vec![b'A'; 300_000]);
        // Damage a block that is not the first, so it is corrupted mid-batch.
        let victim = compressed.len() / 2;
        compressed[victim] ^= 0xff;
        let mut out = Vec::new();
        let result = ParallelBgzfReader::new(Cursor::new(&compressed)).read_to_end(&mut out);
        assert!(result.is_err(), "corruption went unnoticed");

        // Truncation as well.
        let whole = compress(&vec![b'C'; 300_000]);
        let mut out = Vec::new();
        let result = ParallelBgzfReader::new(Cursor::new(whole[..whole.len() / 2].to_vec()))
            .read_to_end(&mut out);
        assert!(result.is_err(), "truncation went unnoticed");
    }

    #[test]
    fn flush_starts_a_new_block() {
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(b"first").unwrap();
        writer.flush().unwrap();
        writer.write_all(b"second").unwrap();
        // index() before finishing sees only the flushed block, which is exactly
        // why finish_with_index exists.
        assert_eq!(writer.index().len(), 2);
        let (compressed, index) = writer.finish_with_index().unwrap();

        // Two data blocks, each a seek point.
        assert_eq!(index.len(), 3); // implicit start + two blocks
        assert_eq!(index.blocks()[1].uncompressed, 5);
        assert_eq!(index.uncompressed_len(), Some(11));

        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
        reader.seek(SeekFrom::Start(5)).unwrap();
        let mut out = Vec::new();
        reader.read_to_end(&mut out).unwrap();
        assert_eq!(out, b"second");
    }
}