gwseq-io 0.2.0

//! Decoding a decompressed data block.
//!
//! Wig section headers and their three item encodings,
//! zoom records, and the [`DataIntervals`] walk that turns either into one
//! stream of [`DataInterval`].
//!
//! Everything borrows. A block is inflated once into a [`Bytes`], and every
//! item is read out of it by offset, so walking a block of a thousand items
//! allocates once for the block and not at all per item.

use bytes::Bytes;

use crate::bytes::LeCursor;
use crate::error::{Error, Result};
use crate::genomic::IndexedLoc;

/// Byte size of the header a wig section starts with (Supp. Table 13).
pub const WIG_HEADER_SIZE: usize = 24;
/// Byte size of a zoom record (Supp. Table 19).
pub const ZOOM_RECORD_SIZE: usize = 32;
/// Bytes of a bed record before its tab-separated tail (Supp. Table 12).
pub const BED_RECORD_HEADER_SIZE: usize = 12;

/// The three wig item encodings, in the order the format numbers them.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum WigEncoding {
    /// Twelve bytes an item: start, end, value.
    BedGraph = 1,
    /// Eight: start, value. One span for the section.
    VarStep = 2,
    /// Four: value. One span and one step for the section.
    FixedStep = 3,
}

impl WigEncoding {
    pub fn from_u8(v: u8) -> Option<Self> {
        match v {
            1 => Some(WigEncoding::BedGraph),
            2 => Some(WigEncoding::VarStep),
            3 => Some(WigEncoding::FixedStep),
            _ => None,
        }
    }

    /// Bytes one item of this encoding occupies.
    pub fn item_size(self) -> usize {
        match self {
            WigEncoding::BedGraph => 12,
            WigEncoding::VarStep => 8,
            WigEncoding::FixedStep => 4,
        }
    }
}

#[derive(Debug, Clone, Copy)]
pub struct WigSectionHeader {
    pub chr_index: u32,
    pub chr_start: i64,
    /// Decoded but not read: the item walk is bounded by `item_count`, not by
    /// this. Kept so the struct is the on-disk header.
    #[allow(dead_code)]
    pub chr_end: i64,
    pub item_step: i64,
    pub item_span: i64,
    pub encoding: WigEncoding,
    pub item_count: u16,
}

pub fn read_wig_header(block: &[u8], path: &str) -> Result<WigSectionHeader> {
    // Checked here rather than by the caller, which needs the item count out of
    // these bytes before it can check anything else about the block.
    if block.len() < WIG_HEADER_SIZE {
        return Err(Error::corrupt(
            path,
            0,
            format!(
                "wig section header needs {WIG_HEADER_SIZE} bytes, its block holds {}",
                block.len()
            ),
        ));
    }
    let mut c = LeCursor::new(block, 0, path);
    let chr_index = c.read_u32()?;
    let chr_start = c.read_u32()? as i64;
    let chr_end = c.read_u32()? as i64;
    let item_step = c.read_u32()? as i64;
    let item_span = c.read_u32()? as i64;
    let type_byte = c.read_u8()?;
    c.skip(1)?; // reserved
    let item_count = c.read_u16()?;
    let encoding = WigEncoding::from_u8(type_byte)
        .ok_or_else(|| Error::corrupt(path, 20, format!("wig data type {type_byte} invalid")))?;
    Ok(WigSectionHeader {
        chr_index,
        chr_start,
        chr_end,
        item_step,
        item_span,
        encoding,
        item_count,
    })
}

/// One value of the file over a genomic range.
///
/// `valid_count` is the bases of the range that actually carry data: equal to
/// the span at full resolution, but a zoom record summarises a fixed-width
/// window that may be mostly empty, and weighting by the span instead would
/// inflate every summed or counted quantification as the zoom level rises.
///
/// `min_value`/`max_value`/`sum_squared` are the record's **own**, not derived
/// from `value`. At a zoom level `value` is the window mean, so taking extremes
/// and spread from it would give the largest and smallest *window* rather than
/// base — both pulled towards the middle, and further with every level — and
/// would leave `sd` with only the variance between windows.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DataInterval {
    pub chr_index: u32,
    pub start: i64,
    pub end: i64,
    pub value: f32,
    pub valid_count: i64,
    pub min_value: f32,
    pub max_value: f32,
    pub sum_squared: f64,
}

/// Read item `index` of a wig section.
///
/// Pulled out of the walk so the writer's encoder has exactly one decoder to be
/// the inverse of, and so the zoom pass can walk a section without standing up
/// a set of loci for it.
pub fn read_wig_item(
    block: &[u8],
    header: &WigSectionHeader,
    index: usize,
    path: &str,
) -> Result<DataInterval> {
    let item_size = header.encoding.item_size();
    let offset = WIG_HEADER_SIZE + index * item_size;
    let mut c = LeCursor::new(block, 0, path);
    c.seek(offset)?;
    let (start, end, value) = match header.encoding {
        WigEncoding::BedGraph => {
            let start = c.read_u32()? as i64;
            let end = c.read_u32()? as i64;
            (start, end, c.read_f32()?)
        }
        WigEncoding::VarStep => {
            let start = c.read_u32()? as i64;
            (start, start + header.item_span, c.read_f32()?)
        }
        WigEncoding::FixedStep => {
            let start = header.chr_start + index as i64 * header.item_step;
            (start, start + header.item_span, c.read_f32()?)
        }
    };
    let valid_count = end - start;
    Ok(DataInterval {
        chr_index: header.chr_index,
        start,
        end,
        value,
        valid_count,
        // One value over the whole range, so it is its own minimum and maximum
        // and the squares sum to it repeated over every base it covers.
        min_value: value,
        max_value: value,
        sum_squared: value as f64 * value as f64 * valid_count as f64,
    })
}

/// A zoom record: a pre-summarised window.
#[derive(Debug, Clone, Copy)]
pub struct ZoomRecord {
    pub chr_index: u32,
    pub chr_start: i64,
    pub chr_end: i64,
    pub valid_count: i64,
    pub min_value: f32,
    pub max_value: f32,
    pub sum_data: f32,
    pub sum_squared: f32,
}

pub fn read_zoom_record(block: &[u8], offset: usize, path: &str) -> Result<ZoomRecord> {
    let mut c = LeCursor::new(block, 0, path);
    c.seek(offset)?;
    Ok(ZoomRecord {
        chr_index: c.read_u32()?,
        chr_start: c.read_u32()? as i64,
        chr_end: c.read_u32()? as i64,
        valid_count: c.read_u32()? as i64,
        min_value: c.read_f32()?,
        max_value: c.read_f32()?,
        sum_data: c.read_f32()?,
        sum_squared: c.read_f32()?,
    })
}

/// A bigBed entry. Always carries its coordinates; the rest of the columns are
/// present only up to the request's `col_count`.
#[derive(Debug, Clone)]
pub struct BedEntry {
    pub chr: String,
    pub start: i64,
    pub end: i64,
    /// The remaining tab-separated columns, in file order, named by the file's
    /// autoSql. Empty when `col_count` was 3.
    pub fields: Vec<(String, String)>,
}

/// How much one block is allowed to inflate to.
///
/// A hard limit, not a hint. Without it a block whose deflate stream expands
/// far beyond what any real writer produces — a corrupt file, or a deliberate
/// one — is inflated until the process runs out of memory, and Rust aborts on
/// a failed allocation rather than returning an error.
const MAX_INFLATED_SIZE: usize = 1 << 30;

/// How much of the output buffer is reserved up front, whatever the header
/// says. The declared size is a hint from the file and a corrupt one can name
/// four gigabytes; reserving that aborts the process before a single byte is
/// inflated. Reserving less costs a few reallocations on a block that really is
/// larger, and `read_to_end` grows past it either way.
const MAX_INFLATE_RESERVE: usize = 1 << 20;

/// Inflate a data block.
///
/// `uncompress_buffer_size` of 0 means the file stores blocks uncompressed, and
/// the block is handed back as it is. Otherwise it is the size the writer
/// declared its largest block inflates to, which is what the output buffer is
/// sized from — but it is a hint from the file, so the decoder is allowed to
/// exceed it rather than truncating silently.
pub fn decompress(block: Bytes, uncompress_buffer_size: u32, path: &str) -> Result<Bytes> {
    decompress_limited(block, uncompress_buffer_size, path, MAX_INFLATED_SIZE)
}

/// The same with the cap as an argument, so the tests can reach it without
/// inflating a gibibyte to prove a limit they can prove with a kilobyte.
fn decompress_limited(
    block: Bytes,
    uncompress_buffer_size: u32,
    path: &str,
    max_size: usize,
) -> Result<Bytes> {
    if uncompress_buffer_size == 0 {
        return Ok(block);
    }
    use std::io::Read;
    let reserve = (uncompress_buffer_size as usize).min(MAX_INFLATE_RESERVE.min(max_size));
    let mut out = Vec::with_capacity(reserve);
    // `take` is what enforces the cap: the decoder stops one byte past the
    // limit and the length says whether it got there, so a stream that would
    // have gone on is refused while inflating rather than after.
    flate2::read::ZlibDecoder::new(&block[..])
        .take(max_size as u64 + 1)
        .read_to_end(&mut out)
        .map_err(|e| Error::corrupt(path, 0, format!("could not inflate data block: {e}")))?;
    if out.len() > max_size {
        return Err(Error::corrupt(
            path,
            0,
            format!("decompressed data exceeds limit ({max_size})"),
        ));
    }
    Ok(Bytes::from(out))
}

/// The bounds of the loci a block is being read for, as (chromosome, base)
/// pairs.
///
/// A block groups records by count rather than by chromosome, so comparing bare
/// coordinates would drop data at every boundary a block straddles.
#[derive(Debug, Clone, Copy)]
struct LocBounds {
    min_chr: u32,
    min_start: i64,
    max_chr: u32,
    max_end: i64,
}

impl LocBounds {
    fn of(locs: &[IndexedLoc], range: std::ops::Range<usize>) -> Self {
        let first = &locs[range.start];
        let last = &locs[range.end - 1];
        // Loci are sorted by chromosome then position, so the last one sits on
        // the highest chromosome; the furthest reach on that chromosome is the
        // widest end among the loci that share it.
        let max_chr = last.chr_index as u32;
        let mut max_end = last.binned_end;
        for loc in locs[range.start..range.end - 1].iter().rev() {
            if loc.chr_index as u32 != max_chr {
                break;
            }
            max_end = max_end.max(loc.binned_end);
        }
        Self {
            min_chr: first.chr_index as u32,
            min_start: first.binned_start,
            max_chr,
            max_end,
        }
    }
}

/// Walks a block's items as [`DataInterval`]s, skipping the ones no locus of
/// the batch can reach and stopping once the block has run past them all.
#[derive(Debug)]
pub struct DataIntervals {
    block: Bytes,
    path: String,
    bounds: LocBounds,
    /// `None` for full data, `Some(_)` for a zoom block.
    header: Option<WigSectionHeader>,
    count: usize,
    index: usize,
}

impl DataIntervals {
    /// `zoom` selects the record layout: zoom blocks are a bare run of 32-byte
    /// records, full-data blocks a wig section header followed by its items.
    pub fn new(
        block: Bytes,
        zoom: bool,
        locs: &[IndexedLoc],
        range: std::ops::Range<usize>,
        path: &str,
    ) -> Result<Self> {
        let (header, count) = if zoom {
            // Counted from the block as it stands: `data_size` is what the
            // block occupies on disk, which is the compressed size on every
            // file that sets `uncompress_buffer_size`.
            (None, block.len() / ZOOM_RECORD_SIZE)
        } else {
            let header = read_wig_header(&block, path)?;
            let count = header.item_count as usize;
            let item_size = header.encoding.item_size();
            // Checked once here rather than per item, so a block whose declared
            // item count does not fit it fails cleanly.
            if WIG_HEADER_SIZE + count * item_size > block.len() {
                return Err(Error::corrupt(
                    path,
                    0,
                    format!(
                        "wig section declares {count} items of type {:?}, which do not fit \
                         its {} byte block",
                        header.encoding,
                        block.len()
                    ),
                ));
            }
            (Some(header), count)
        };
        Ok(Self {
            block,
            path: path.to_string(),
            bounds: LocBounds::of(locs, range),
            header,
            count,
            index: 0,
        })
    }
}

impl Iterator for DataIntervals {
    type Item = Result<DataInterval>;

    fn next(&mut self) -> Option<Self::Item> {
        while self.index < self.count {
            let index = self.index;
            self.index += 1;

            let data = match self.header {
                Some(header) => match read_wig_item(&self.block, &header, index, &self.path) {
                    Ok(d) => d,
                    Err(e) => return Some(Err(e)),
                },
                None => {
                    let record =
                        match read_zoom_record(&self.block, index * ZOOM_RECORD_SIZE, &self.path) {
                            Ok(r) => r,
                            Err(e) => return Some(Err(e)),
                        };
                    // A window the file recorded but valued from nothing has no
                    // mean to give.
                    if record.valid_count == 0 {
                        continue;
                    }
                    DataInterval {
                        chr_index: record.chr_index,
                        start: record.chr_start,
                        end: record.chr_end,
                        value: record.sum_data / record.valid_count as f32,
                        valid_count: record.valid_count,
                        min_value: record.min_value,
                        max_value: record.max_value,
                        sum_squared: record.sum_squared as f64,
                    }
                }
            };

            let b = &self.bounds;
            if data.chr_index < b.min_chr {
                continue;
            }
            if data.chr_index == b.min_chr && data.end <= b.min_start {
                continue;
            }
            if data.chr_index > b.max_chr {
                break;
            }
            if data.chr_index == b.max_chr && data.start >= b.max_end {
                break;
            }
            return Some(Ok(data));
        }
        None
    }
}

/// Walk a bigBed block's records (Supp. Table 12), keeping the loci-relevant
/// ones and parsing only the columns asked for.
///
/// A record is `u32 chromIx, u32 start, u32 end` followed by a NUL-terminated,
/// tab-separated tail. The fields past `kept` are counted but never built,
/// which is what makes a `col_count`-limited read cheaper than a full one
/// rather than merely narrower.
#[derive(Debug)]
pub struct BedRecords<'a> {
    block: Bytes,
    path: &'a str,
    /// The file's column names, coordinates included. The tail's fields are
    /// named by the entries past the first three.
    auto_sql: &'a indexmap::IndexMap<String, String>,
    /// Tail fields the file declares, i.e. `auto_sql.len() - 3`.
    field_count: usize,
    /// Tail fields to build.
    kept: usize,
    bounds: LocBounds,
    offset: usize,
}

impl<'a> BedRecords<'a> {
    /// `col_count` counts the three coordinate columns; 0 keeps everything.
    /// Clamped rather than checked — the readers validate what they were handed
    /// against the file's own column count before they get here.
    pub fn new(
        block: Bytes,
        auto_sql: &'a indexmap::IndexMap<String, String>,
        col_count: usize,
        locs: &[IndexedLoc],
        range: std::ops::Range<usize>,
        path: &'a str,
    ) -> Result<Self> {
        if auto_sql.len() < 3 {
            return Err(Error::format(
                path,
                format!(
                    "bed entries need the 3 standard fields, autosql describes {}",
                    auto_sql.len()
                ),
            ));
        }
        let field_count = auto_sql.len() - 3;
        let kept = if col_count == 0 {
            field_count
        } else {
            field_count.min(col_count.saturating_sub(3))
        };
        Ok(Self {
            block,
            path,
            auto_sql,
            field_count,
            kept,
            bounds: LocBounds::of(locs, range),
            offset: 0,
        })
    }

    /// Read the tab-separated tail into named fields.
    ///
    /// A trailing delimiter closes a field of its own, so `"name\t0\t"` is three
    /// fields, the last blank — dropping it would fail the record for a field
    /// count mismatch it does not have. By the same count an empty tail is one
    /// blank field, since n fields carry n-1 delimiters; that is also what a
    /// bed3 record looks like, and the two cannot be told apart from the bytes.
    /// So the autoSql decides, and only a file declaring no tail fields reads an
    /// empty tail as none.
    fn read_fields(&self, tail: &[u8]) -> Result<Vec<(String, String)>> {
        let mut fields = Vec::with_capacity(self.kept);
        let mut found = 0usize;
        if !tail.is_empty() || self.field_count > 0 {
            for part in tail.split(|b| *b == b'\t') {
                if found < self.kept {
                    let name = self
                        .auto_sql
                        .get_index(3 + found)
                        .map(|(k, _)| k.clone())
                        .unwrap_or_else(|| format!("field{}", 4 + found));
                    fields.push((name, String::from_utf8_lossy(part).into_owned()));
                }
                found += 1;
            }
        }
        if found != self.field_count {
            return Err(Error::corrupt(
                self.path,
                self.offset as u64,
                format!(
                    "invalid bed entry (found {found} fields past the first 3, \
                     autosql declares {})",
                    self.field_count
                ),
            ));
        }
        Ok(fields)
    }
}

impl Iterator for BedRecords<'_> {
    type Item = Result<(u32, i64, i64, Vec<(String, String)>)>;

    fn next(&mut self) -> Option<Self::Item> {
        while self.offset < self.block.len() {
            if self.offset + BED_RECORD_HEADER_SIZE > self.block.len() {
                return Some(Err(Error::corrupt(
                    self.path,
                    self.offset as u64,
                    format!(
                        "truncated bed record at {} ({} bytes left in its block)",
                        self.offset,
                        self.block.len() - self.offset
                    ),
                )));
            }
            let head = &self.block[self.offset..self.offset + BED_RECORD_HEADER_SIZE];
            let chr_index = u32::from_le_bytes([head[0], head[1], head[2], head[3]]);
            let start = u32::from_le_bytes([head[4], head[5], head[6], head[7]]) as i64;
            let end = u32::from_le_bytes([head[8], head[9], head[10], head[11]]) as i64;

            let tail_start = self.offset + BED_RECORD_HEADER_SIZE;
            let Some(nul) = memchr::memchr(0, &self.block[tail_start..]) else {
                return Some(Err(Error::corrupt(
                    self.path,
                    tail_start as u64,
                    "invalid bed entry (null terminator not found)",
                )));
            };
            let tail_end = tail_start + nul;

            // Bounds first: a record no locus of this batch can reach costs one
            // comparison rather than a field split.
            //
            // `reach` rather than `end`: BED allows `start == end` — an
            // insertion — and a half-open interval of no width overlaps
            // nothing at all, so such an entry would be dropped by every test
            // it met. It occupies the single base it names instead, which is
            // the rule `advance_cursor` and `walk_bed_batch` also apply.
            let b = &self.bounds;
            let reach = end.max(start + 1);
            let reachable =
                !(chr_index < b.min_chr || (chr_index == b.min_chr && reach <= b.min_start));
            let past = chr_index > b.max_chr || (chr_index == b.max_chr && start >= b.max_end);

            let fields = if reachable && !past {
                match self.read_fields(&self.block[tail_start..tail_end]) {
                    Ok(f) => Some(f),
                    Err(e) => return Some(Err(e)),
                }
            } else {
                None
            };
            self.offset = tail_end + 1;

            if past {
                break;
            }
            if let Some(fields) = fields {
                return Some(Ok((chr_index, start, end, fields)));
            }
        }
        None
    }
}

/// Call `visit(chr_index, start, end)` for every record of a bed block,
/// reading the coordinates alone.
///
/// What the zoom pass wants: the tail is skipped rather than split, which is all
/// the coverage a bigBed's summaries need.
pub fn visit_bed_records(
    block: &[u8],
    path: &str,
    mut visit: impl FnMut(u32, i64, i64),
) -> Result<()> {
    let mut offset = 0usize;
    while offset < block.len() {
        if offset + BED_RECORD_HEADER_SIZE > block.len() {
            return Err(Error::corrupt(
                path,
                offset as u64,
                format!(
                    "truncated bed record at {offset} ({} bytes left in its block)",
                    block.len() - offset
                ),
            ));
        }
        let head = &block[offset..offset + BED_RECORD_HEADER_SIZE];
        let chr_index = u32::from_le_bytes([head[0], head[1], head[2], head[3]]);
        let start = u32::from_le_bytes([head[4], head[5], head[6], head[7]]) as i64;
        let end = u32::from_le_bytes([head[8], head[9], head[10], head[11]]) as i64;
        let tail_start = offset + BED_RECORD_HEADER_SIZE;
        let Some(nul) = memchr::memchr(0, &block[tail_start..]) else {
            return Err(Error::corrupt(
                path,
                tail_start as u64,
                "invalid bed entry (null terminator not found)",
            ));
        };
        offset = tail_start + nul + 1;
        visit(chr_index, start, end);
    }
    Ok(())
}

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

    /// Deflate `bytes` the way a bbi writer would, so `decompress` can read it.
    fn deflated(bytes: &[u8]) -> Vec<u8> {
        use std::io::Write as _;
        let mut e = flate2::write::ZlibEncoder::new(Vec::new(), flate2::Compression::best());
        e.write_all(bytes).unwrap();
        e.finish().unwrap()
    }

    #[test]
    fn a_declared_buffer_size_is_a_hint_and_not_an_allocation() {
        // Four gigabytes, which is what a mutated header can name and what
        // reserving would abort the process for. The block is 32 bytes.
        let block = bytes::Bytes::from(deflated(&[7u8; 32]));
        let out = decompress(block, u32::MAX, "corrupt.bigwig").unwrap();
        assert_eq!(out.len(), 32);
        assert!(out.iter().all(|b| *b == 7));
    }

    #[test]
    fn a_block_that_inflates_past_the_limit_is_refused_rather_than_read() {
        // The cheap half of a decompression bomb: a run of zeros deflates to
        // almost nothing, so refusing it has to happen while inflating rather
        // than after. The real cap is a gibibyte; proving the rule does not
        // need one.
        let block = bytes::Bytes::from(deflated(&vec![0u8; 4097]));
        let err = decompress_limited(block, 4096, "bomb.bigwig", 4096)
            .unwrap_err()
            .to_string();
        assert!(err.contains("exceeds limit (4096)"), "{err}");
    }

    #[test]
    fn a_block_exactly_at_the_limit_is_still_read() {
        let block = bytes::Bytes::from(deflated(&vec![0u8; 4096]));
        assert_eq!(
            decompress_limited(block, 4096, "big.bigwig", 4096)
                .unwrap()
                .len(),
            4096
        );
    }

    fn loc(chr: usize, start: i64, end: i64) -> IndexedLoc {
        IndexedLoc {
            chr_index: chr,
            start,
            end,
            binned_start: start,
            binned_end: end,
            bin_size: 1.0,
            reverse: false,
            output_start: 0,
            output_end: (end - start) as usize,
        }
    }

    fn wig_block(
        encoding: WigEncoding,
        items: &[(u32, u32, f32)],
        span: u32,
        step: u32,
    ) -> Vec<u8> {
        let mut b = Vec::new();
        b.extend_from_slice(&7u32.to_le_bytes()); // chromId
        b.extend_from_slice(&items[0].0.to_le_bytes()); // chromStart
        b.extend_from_slice(&items[items.len() - 1].1.to_le_bytes()); // chromEnd
        b.extend_from_slice(&step.to_le_bytes());
        b.extend_from_slice(&span.to_le_bytes());
        b.push(encoding as u8);
        b.push(0);
        b.extend_from_slice(&(items.len() as u16).to_le_bytes());
        assert_eq!(b.len(), WIG_HEADER_SIZE);
        for (start, end, value) in items {
            match encoding {
                WigEncoding::BedGraph => {
                    b.extend_from_slice(&start.to_le_bytes());
                    b.extend_from_slice(&end.to_le_bytes());
                    b.extend_from_slice(&value.to_le_bytes());
                }
                WigEncoding::VarStep => {
                    b.extend_from_slice(&start.to_le_bytes());
                    b.extend_from_slice(&value.to_le_bytes());
                }
                WigEncoding::FixedStep => b.extend_from_slice(&value.to_le_bytes()),
            }
        }
        b
    }

    fn collect(block: Vec<u8>, zoom: bool, locs: &[IndexedLoc]) -> Vec<DataInterval> {
        DataIntervals::new(Bytes::from(block), zoom, locs, 0..locs.len(), "test")
            .unwrap()
            .map(|r| r.unwrap())
            .collect()
    }

    #[test]
    fn the_three_encodings_decode_to_the_same_intervals() {
        let items = [(100u32, 110u32, 1.5f32), (110, 120, 2.5), (120, 130, 3.5)];
        let locs = [loc(7, 0, 1000)];

        let bg = collect(
            wig_block(WigEncoding::BedGraph, &items, 10, 10),
            false,
            &locs,
        );
        let vs = collect(
            wig_block(WigEncoding::VarStep, &items, 10, 10),
            false,
            &locs,
        );
        let fs = collect(
            wig_block(WigEncoding::FixedStep, &items, 10, 10),
            false,
            &locs,
        );

        assert_eq!(bg.len(), 3);
        assert_eq!(bg, vs);
        assert_eq!(bg, fs);
        assert_eq!(bg[0].start, 100);
        assert_eq!(bg[0].end, 110);
        assert_eq!(bg[0].value, 1.5);
        assert_eq!(bg[0].valid_count, 10);
        assert_eq!(bg[0].min_value, 1.5);
        assert_eq!(bg[0].sum_squared, 1.5f64 * 1.5 * 10.0);
        assert_eq!(bg[0].chr_index, 7);
    }

    #[test]
    fn items_outside_the_batchs_bounds_are_skipped_and_the_walk_stops() {
        let items = [
            (0u32, 10u32, 1.0f32),
            (10, 20, 2.0),
            (100, 110, 3.0),
            (200, 210, 4.0),
            (300, 310, 5.0),
        ];
        // Only [100, 210) is wanted.
        let locs = [loc(7, 100, 210)];
        let got = collect(
            wig_block(WigEncoding::BedGraph, &items, 10, 10),
            false,
            &locs,
        );
        assert_eq!(got.iter().map(|d| d.start).collect::<Vec<_>>(), [100, 200]);
    }

    #[test]
    fn a_block_straddling_two_chromosomes_keeps_only_the_reachable_side() {
        // Items on chr 7 then chr 8; the batch only wants chr 8.
        let mut b = wig_block(WigEncoding::BedGraph, &[(0, 10, 1.0)], 10, 10);
        b[0..4].copy_from_slice(&8u32.to_le_bytes());
        let locs = [loc(8, 0, 10)];
        assert_eq!(collect(b.clone(), false, &locs).len(), 1);
        let locs = [loc(9, 0, 10)];
        assert!(collect(b, false, &locs).is_empty());
    }

    #[test]
    fn zoom_records_become_intervals_and_zero_count_ones_are_dropped() {
        let mut b = Vec::new();
        for (start, end, valid, min, max, sum, sq) in [
            (0u32, 100u32, 100u32, 1.0f32, 5.0f32, 300.0f32, 1000.0f32),
            (100, 200, 0, 0.0, 0.0, 0.0, 0.0), // recorded but valued from nothing
            (200, 300, 50, 2.0, 4.0, 150.0, 500.0),
        ] {
            b.extend_from_slice(&7u32.to_le_bytes());
            b.extend_from_slice(&start.to_le_bytes());
            b.extend_from_slice(&end.to_le_bytes());
            b.extend_from_slice(&valid.to_le_bytes());
            b.extend_from_slice(&min.to_le_bytes());
            b.extend_from_slice(&max.to_le_bytes());
            b.extend_from_slice(&sum.to_le_bytes());
            b.extend_from_slice(&sq.to_le_bytes());
        }
        let locs = [loc(7, 0, 1000)];
        let got = collect(b, true, &locs);
        assert_eq!(got.len(), 2);
        // value is the window mean, not the stored sum.
        assert_eq!(got[0].value, 3.0);
        assert_eq!(got[0].valid_count, 100);
        assert_eq!(got[0].min_value, 1.0);
        assert_eq!(got[0].max_value, 5.0);
        assert_eq!(got[1].value, 3.0);
        assert_eq!(got[1].valid_count, 50);
    }

    #[test]
    fn a_bad_wig_type_is_corrupt_not_a_panic() {
        let mut b = wig_block(WigEncoding::BedGraph, &[(0, 10, 1.0)], 10, 10);
        b[20] = 9;
        let locs = [loc(7, 0, 10)];
        let err = DataIntervals::new(Bytes::from(b), false, &locs, 0..1, "test").unwrap_err();
        assert!(err.to_string().contains("wig data type 9 invalid"), "{err}");
    }

    #[test]
    fn a_block_declaring_more_items_than_it_holds_is_refused_up_front() {
        let mut b = wig_block(WigEncoding::BedGraph, &[(0, 10, 1.0)], 10, 10);
        b[22..24].copy_from_slice(&500u16.to_le_bytes());
        let locs = [loc(7, 0, 10)];
        let err = DataIntervals::new(Bytes::from(b), false, &locs, 0..1, "test").unwrap_err();
        assert!(err.to_string().contains("do not fit"), "{err}");
    }

    #[test]
    fn a_block_too_short_for_a_header_is_refused() {
        let locs = [loc(7, 0, 10)];
        let err =
            DataIntervals::new(Bytes::from(vec![0u8; 8]), false, &locs, 0..1, "test").unwrap_err();
        assert!(err.to_string().contains("header needs 24 bytes"), "{err}");
    }

    #[test]
    fn an_uncompressed_file_hands_its_block_back_untouched() {
        let block = Bytes::from(vec![1u8, 2, 3]);
        assert_eq!(decompress(block.clone(), 0, "test").unwrap(), block);
    }

    #[test]
    fn a_compressed_block_round_trips() {
        use std::io::Write;
        let raw: Vec<u8> = (0..5000).map(|i| (i % 251) as u8).collect();
        let mut encoder = flate2::write::ZlibEncoder::new(Vec::new(), flate2::Compression::new(6));
        encoder.write_all(&raw).unwrap();
        let compressed = Bytes::from(encoder.finish().unwrap());
        assert_eq!(&decompress(compressed, 8192, "test").unwrap()[..], &raw[..]);
    }

    #[test]
    fn garbage_where_a_deflate_stream_should_be_is_corrupt() {
        let err = decompress(Bytes::from(vec![9u8; 64]), 4096, "test").unwrap_err();
        assert!(err.to_string().contains("could not inflate"), "{err}");
    }

    #[test]
    fn bounds_take_the_widest_end_on_the_last_chromosome() {
        // Sorted by (chr, start): the last locus is not the one reaching
        // furthest, which is the case a naive `last.binned_end` gets wrong.
        let locs = [loc(1, 0, 10), loc(2, 0, 500), loc(2, 100, 200)];
        let bounds = LocBounds::of(&locs, 0..3);
        assert_eq!(bounds.min_chr, 1);
        assert_eq!(bounds.min_start, 0);
        assert_eq!(bounds.max_chr, 2);
        assert_eq!(bounds.max_end, 500);
    }
}