gwseq_io/bbi/

writer.rs

//! The bigWig / bigBed writer.
//!
//! Three passes over one file:
//!
//! 1. **Data.** Values or entries accumulate into sections
//!    ([`super::section`]); a full section is compressed on the pool and
//!    appended, and its bounds recorded as an R-tree leaf item.
//! 2. **Zoom.** Coarser levels are built by reading the data section back —
//!    through the same handle, since [`LocalSink`] is open both ways.
//! 3. **Index and headers.** The R-tree is written, then the reserved header
//!    fields are patched with the offsets and counts now known.
//!
//! The file goes out in one forward pass with no temporary files. Every part of
//! a bbi file is reached through an explicit offset in its 64-byte header, so
//! the parts need not sit in the order Supp. Table 4 lists them in:
//!
//! ```text
//! 0     bbiHeader          64 bytes, patched at close, magic stamped last
//! 64    zoom header slots  24 * MAX_ZOOM_LEVELS, patched at close
//! 304   autoSql            bigBed only, known at open
//!       totalSummary       40 bytes, patched at close
//!       dataCount          8 bytes, patched at close
//!       data blocks        streamed as they fill
//!       data R-tree
//!       zoom levels        zoomCount, records, R-tree, per level
//!       chromosome B+tree
//!       trailing magic
//! ```
//!
//! The chromosome tree comes last because it is the last thing known: sizes and
//! ids are inferred from what was written. Until the magic is stamped the
//! header reads as zeroes and every reader refuses the file, so an abandoned
//! write cannot be mistaken for a short but valid one.
//!
//! Zoom levels are built by reading the data section back out of the file just
//! written, each level after the first out of the level below it. That costs
//! one inflate of what was already deflated and keeps the levels contiguous,
//! with a reduction ladder chosen from the finished data rather than guessed
//! from its first few thousand items.
//!
//! [`LocalSink`]: crate::source::LocalSink

use std::collections::BinaryHeap;

use indexmap::IndexMap;

use crate::bbi::block::{read_wig_header, read_wig_item, WigEncoding};
use crate::bbi::chr_tree::WriteEntry;
use crate::bbi::header::{
    build_auto_sql, to_bbi_f32, to_bbi_u32, write_header, write_total_summary, write_zoom_header,
    BbiHeader, BbiKind, TotalSummary, ZoomHeader, BBI_HEADER_SIZE, BBI_OUTPUT_VERSION,
    TOTAL_SUMMARY_SIZE, ZOOM_HEADER_SIZE,
};
use crate::bbi::rtree::{LeafItem, TREE_BLOCK_SIZE};
use crate::bbi::section::{Accept, CostModel, SectionPolicy, WigSection, MAX_ITEMS_PER_SECTION};
use crate::error::{Error, Result};
use crate::genomic::ChrMap;
use crate::parallel::{resolve_parallel, Executor, Promise};
use crate::source::{ByteSink, ByteSource, LocalSink};

/// Items a section or a block holds by default, matching UCSC's writers. Also
/// the unit the R-tree indexes, since each of them is compressed and fetched on
/// its own. A bed entry is bigger and more variable than a wig value, which is
/// why fewer of them go in a block.
pub const WIG_ITEMS_PER_SLOT: usize = 1024;
pub const BED_ITEMS_PER_SLOT: usize = 512;

/// zlib level for data and zoom blocks. A bbi file is written once and read
/// many times, and 6 is what UCSC's writers use.
pub const COMPRESSION_LEVEL: u32 = 6;

/// Items whose spans are averaged to anchor the zoom reduction ladder. The
/// ladder has to be fixed before the candidate levels can be counted, which is
/// why it comes from a sample rather than from the whole file.
const RESOLUTION_SAMPLE_ITEMS: u64 = 4096;

/// The eight bytes reserved after the summary for the data count.
const DATA_COUNT_SIZE: u64 = 8;

/// Slots the header reserves for zoom levels.
const MAX_ZOOM_LEVELS: usize = 10;
/// The first level summarises about ten items, and each level after it four
/// times as much.
const INITIAL_ZOOM_FACTOR: i64 = 10;
const ZOOM_INCREMENT: i64 = 4;
/// Bytes one zoom record takes (Supp. Table 19).
const ZOOM_RECORD_SIZE: u64 = 32;

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FieldType {
    String,
    Int,
    Uint,
    Float,
}

impl FieldType {
    pub fn as_str(self) -> &'static str {
        match self {
            FieldType::String => "string",
            FieldType::Int => "int",
            FieldType::Uint => "uint",
            FieldType::Float => "float",
        }
    }
}

pub struct BbiWriterOptions {
    pub kind: BbiKind,
    /// Declared sizes. Every written coordinate is checked against them, and a
    /// written id resolves against their keys the way a read one does. `None`
    /// infers each chromosome's size from what is written.
    pub chr_sizes: Option<ChrMap>,
    /// bigBed only. The first three must be the coordinates.
    pub fields: IndexMap<String, String>,
    pub items_per_slot: Option<usize>,
    pub block_size: u32,
    pub compression_level: u32,
    pub parallel: i64,
    /// How a wig section decides between widening its encoding and closing.
    /// Every setting writes a valid file; they differ in how big it is. See
    /// [`SectionPolicy`] and `examples/section_policy.rs`.
    pub section_policy: SectionPolicy,
    /// The constants behind [`SectionPolicy::Cost`]. Exposed for the same
    /// reason: so the model can be measured rather than argued about.
    pub cost_model: CostModel,
}

impl Default for BbiWriterOptions {
    fn default() -> Self {
        Self {
            kind: BbiKind::BigWig,
            chr_sizes: None,
            fields: IndexMap::new(),
            items_per_slot: None,
            block_size: TREE_BLOCK_SIZE,
            compression_level: COMPRESSION_LEVEL,
            parallel: -1,
            section_policy: SectionPolicy::default(),
            cost_model: CostModel::default(),
        }
    }
}

/// How many sections were written in each encoding. What to compare first
/// when two writes of the same values differ in size — see `bbi/section.rs`,
/// which is private and so not linked from here.
#[derive(Debug, Clone, Copy, Default)]
pub struct SectionCounts {
    pub bedgraph: u64,
    pub varstep: u64,
    pub fixedstep: u64,
}

impl SectionCounts {
    pub fn total(&self) -> u64 {
        self.bedgraph + self.varstep + self.fixedstep
    }

    fn bump(&mut self, encoding: WigEncoding) {
        match encoding {
            WigEncoding::BedGraph => self.bedgraph += 1,
            WigEncoding::VarStep => self.varstep += 1,
            WigEncoding::FixedStep => self.fixedstep += 1,
        }
    }
}

/// A chromosome as the writer has seen it.
#[derive(Debug, Clone, Copy)]
struct ChrState {
    index: u32,
    size: i64,
    declared: bool,
}

/// A block whose deflate is still running, and the index entry waiting on it.
///
/// Everything about the entry but the offset and the size is known before the
/// block is compressed, so those two are all the commit has left to fill in.
struct PendingBlock {
    leaf: LeafItem,
    /// A one-shot: the worker sends, the commit receives. `Ok` on the calling
    /// thread when there is no pool.
    result: PendingResult,
    /// Section type of a wig block, `None` for a bed or zoom one.
    encoding: Option<WigEncoding>,
}

enum PendingResult {
    InFlight(Promise<Result<Vec<u8>>>),
}

impl PendingResult {
    fn take(self) -> Result<Vec<u8>> {
        match self {
            // A worker that panicked drops its half without filling it; the
            // panic itself has already been reported by rayon, and this turns
            // the silence into an error of the call that is running.
            PendingResult::InFlight(promise) => promise
                .wait()
                .unwrap_or_else(|| Err(Error::invalid("a block failed to compress"))),
        }
    }
}

pub struct BbiWriter {
    path: String,
    kind: BbiKind,
    sink: Option<LocalSink>,
    executor: Option<Executor>,
    block_size: u32,
    items_per_slot: usize,
    compression_level: u32,

    // Offsets of the reserved prefix, fixed once at open.
    auto_sql_offset: u64,
    total_summary_offset: u64,
    full_data_offset: u64,
    full_index_offset: u64,
    chr_tree_offset: u64,

    // bigBed only: the columns of an entry, and what the header says of them.
    bed_fields: IndexMap<String, String>,
    field_count: u16,
    defined_field_count: u16,

    declared_sizes: Option<ChrMap>,
    chrs: IndexMap<String, ChrState>,
    current_chr: Option<String>,
    last_chr_id: String,
    last_chr_index: Option<u32>,
    last_chr_end: i64,
    last_entry_start: i64,

    section: WigSection,

    // bigBed only: the block being filled, its bounds, and the sweep turning
    // the entries into the coverage the summary and the zoom levels hold.
    bed_block: Vec<u8>,
    bed_block_count: usize,
    bed_bounds: Option<(u32, u32, u32, u32)>,
    coverage: CoverageSweep,

    data_items: Vec<LeafItem>,
    /// Where a placed block is indexed, and whether it counts as a section.
    ///
    /// `None` during the data pass: a block goes into `data_items` and bumps
    /// the encoding counters. `Some` during a zoom pass: it goes into that
    /// level's own items and counts towards nothing, a zoom block being an
    /// index entry rather than a section. One slot is enough because the two
    /// never overlap — the zoom pass runs inside `close()`, after the data
    /// pipeline has drained.
    zoom_items: Option<Vec<LeafItem>>,
    pending: std::collections::VecDeque<PendingBlock>,
    pending_limit: usize,
    section_count: u64,
    uncompress_buffer_size: u64,
    data_body_size: u64,
    section_counts: SectionCounts,

    summary: TotalSummary,
    item_count: u64,
    entry_count: u64,
    /// Values dropped for not being finite. A NaN or infinity leaves a gap,
    /// which is what a bigWig means by a base carrying no data.
    skipped_count: u64,
    /// Values cut back to the end of a declared chromosome they hung over. A
    /// value that *starts* at or past the end is an error instead — a span that
    /// does not divide a chromosome is ordinary, a value off the end is not.
    clipped_count: u64,

    ladder_frozen: bool,
    ladder_start_item_count: u64,
    zoom_reductions: [i64; MAX_ZOOM_LEVELS],
    zoom_res_sizes: [u64; MAX_ZOOM_LEVELS],
    zoom_res_ends: [i64; MAX_ZOOM_LEVELS],
    zoom_res_chr_index: Option<u32>,
    zoom_headers: Vec<ZoomHeader>,

    closed: bool,
    /// Set while a call is part-way through, so a `Drop` after a failure leaves
    /// the file unfinished rather than turning the caller's error into a file
    /// that looks complete.
    failed: bool,
}

impl BbiWriter {
    pub fn create(path: &str, options: BbiWriterOptions) -> Result<Self> {
        if crate::source::is_url(path) {
            return Err(Error::invalid(format!(
                "{path} is a url, which cannot be written to"
            )));
        }
        if options.kind == BbiKind::BigWig && !options.fields.is_empty() {
            return Err(Error::invalid("fields is only supported for bigbed files"));
        }
        let items_per_slot = options.items_per_slot.unwrap_or(match options.kind {
            BbiKind::BigWig => WIG_ITEMS_PER_SLOT,
            BbiKind::BigBed => BED_ITEMS_PER_SLOT,
        });
        // A wig section states its item count on 16 bits, which is the cap; a
        // bed block has no count of its own, but one rule for both is easier to
        // explain.
        if !(1..=MAX_ITEMS_PER_SECTION).contains(&items_per_slot) {
            return Err(Error::invalid(format!(
                "items_per_slot {items_per_slot} invalid (1 to {MAX_ITEMS_PER_SECTION}, \
                 or -1 for the default)"
            )));
        }
        if options.block_size < 2 {
            return Err(Error::invalid(format!(
                "block_size {} invalid (>= 2)",
                options.block_size
            )));
        }
        if options.compression_level > 9 {
            return Err(Error::invalid(format!(
                "compression_level {} invalid (0 to 9)",
                options.compression_level
            )));
        }
        if let Some(sizes) = &options.chr_sizes {
            for entry in sizes.iter() {
                if entry.size <= 0 {
                    return Err(Error::invalid(format!(
                        "size {} of chromosome {} must be positive",
                        entry.size, entry.id
                    )));
                }
                to_bbi_u32(entry.size, "chromSize")?;
            }
        }

        // Deflate is nearly the whole cost of writing a block, and blocks are
        // independent, so this is the one part of the write path worth more
        // than one thread. A single worker starts none: the block is then
        // compressed in place and the pipeline stays empty.
        let parallel = resolve_parallel(options.parallel);
        let (executor, pending_limit) = if parallel > 1 {
            // Blocks in flight, enough to keep every worker fed while the
            // oldest is waited on. A wig block runs to a few kilobytes, so the
            // whole queue is well under a megabyte, and it is what bounds the
            // memory a producer faster than the workers can run the writer into.
            (Some(Executor::new(options.parallel)?), parallel * 4)
        } else {
            (None, 0)
        };

        // The columns are known now and never change, so a bigBed's autoSql is
        // the one part of the prefix that can be written rather than reserved.
        let mut bed_fields = options.fields;
        let mut auto_sql_text = String::new();
        let (mut field_count, mut defined_field_count) = (0u16, 0u16);
        if options.kind == BbiKind::BigBed {
            if bed_fields.is_empty() {
                bed_fields = [
                    ("chr", "string"),
                    ("start", "uint"),
                    ("end", "uint"),
                    ("name", "string"),
                ]
                .into_iter()
                .map(|(a, b)| (a.to_string(), b.to_string()))
                .collect();
            }
            let described = build_auto_sql(&bed_fields)?;
            auto_sql_text = described.text;
            field_count = described.field_count;
            defined_field_count = described.defined_field_count;
        }

        let mut sink = LocalSink::create(path)?;
        // The rest of the prefix goes out as zeroes, so an abandoned file
        // carries no magic and reads as "not a bigwig or bigbed file" rather
        // than as a valid header pointing at data that was never written.
        let prefix_size = BBI_HEADER_SIZE + MAX_ZOOM_LEVELS as u64 * ZOOM_HEADER_SIZE;
        let mut prefix = vec![0u8; prefix_size as usize];
        let mut auto_sql_offset = 0;
        if !auto_sql_text.is_empty() {
            auto_sql_offset = prefix_size;
            prefix.extend_from_slice(auto_sql_text.as_bytes());
            prefix.push(0);
        }
        let total_summary_offset = prefix.len() as u64;
        let full_data_offset = total_summary_offset + TOTAL_SUMMARY_SIZE;
        prefix.resize(
            prefix.len() + (TOTAL_SUMMARY_SIZE + DATA_COUNT_SIZE) as usize,
            0,
        );
        sink.append(&prefix)?;

        Ok(Self {
            path: path.to_string(),
            kind: options.kind,
            sink: Some(sink),
            executor,
            block_size: options.block_size,
            items_per_slot,
            compression_level: options.compression_level,
            auto_sql_offset,
            total_summary_offset,
            full_data_offset,
            full_index_offset: 0,
            chr_tree_offset: 0,
            bed_fields,
            field_count,
            defined_field_count,
            declared_sizes: options.chr_sizes,
            chrs: IndexMap::new(),
            current_chr: None,
            last_chr_id: String::new(),
            last_chr_index: None,
            last_chr_end: 0,
            last_entry_start: -1,
            section: WigSection::with_policy(
                items_per_slot,
                options.section_policy,
                options.cost_model,
            ),
            bed_block: Vec::new(),
            bed_block_count: 0,
            bed_bounds: None,
            coverage: CoverageSweep::default(),
            data_items: Vec::new(),
            zoom_items: None,
            pending: std::collections::VecDeque::new(),
            pending_limit,
            section_count: 0,
            uncompress_buffer_size: 0,
            data_body_size: 0,
            section_counts: SectionCounts::default(),
            summary: TotalSummary::default(),
            item_count: 0,
            entry_count: 0,
            skipped_count: 0,
            clipped_count: 0,
            ladder_frozen: false,
            ladder_start_item_count: 0,
            zoom_reductions: [0; MAX_ZOOM_LEVELS],
            zoom_res_sizes: [0; MAX_ZOOM_LEVELS],
            zoom_res_ends: [0; MAX_ZOOM_LEVELS],
            zoom_res_chr_index: None,
            zoom_headers: Vec::new(),
            closed: false,
            failed: false,
        })
    }

    // -- byte plumbing -----------------------------------------------------

    fn sink(&mut self) -> Result<&mut LocalSink> {
        self.sink.as_mut().ok_or_else(|| Error::Closed {
            path: self.path.clone(),
        })
    }

    fn emit(&mut self, bytes: &[u8]) -> Result<()> {
        let path = self.path.clone();
        self.sink
            .as_mut()
            .ok_or(Error::Closed { path })?
            .append(bytes)?;
        Ok(())
    }

    /// Overwrite bytes already placed, for a field reserved and filled later.
    fn patch(&mut self, offset: u64, bytes: &[u8]) -> Result<()> {
        self.drain()?;
        let path = self.path.clone();
        self.sink
            .as_mut()
            .ok_or(Error::Closed { path })?
            .write_all_at(offset, bytes)
    }

    /// Offset the next emitted byte will land at, once nothing is in flight.
    ///
    /// Only meaningful with the pipeline empty: a block still compressing has
    /// room reserved for it here that it has not taken yet.
    fn sync_cursor(&mut self) -> Result<u64> {
        self.drain()?;
        Ok(self.sink()?.position())
    }

    fn is_bigwig(&self) -> bool {
        self.kind == BbiKind::BigWig
    }

    // -- block pipeline ----------------------------------------------------

    /// Compress a block, or hand it back unchanged on an uncompressed file.
    fn pack_block(body: Vec<u8>, level: u32, path: &str) -> Result<Vec<u8>> {
        if level == 0 {
            return Ok(body);
        }
        use std::io::Write as _;
        let mut encoder =
            flate2::write::ZlibEncoder::new(Vec::new(), flate2::Compression::new(level));
        encoder
            .write_all(&body)
            .and_then(|_| encoder.finish())
            .map_err(|e| Error::io(path, e))
    }

    /// Place a compressed block at the running offset and index it.
    ///
    /// The one place an offset is handed out, and so the one place the order of
    /// the file is decided.
    fn place_block(
        &mut self,
        mut leaf: LeafItem,
        block: &[u8],
        encoding: Option<WigEncoding>,
    ) -> Result<()> {
        leaf.offset = self.sink()?.position();
        leaf.size = block.len() as u64;
        let zooming = match &mut self.zoom_items {
            Some(items) => {
                items.push(leaf);
                true
            }
            None => {
                self.data_items.push(leaf);
                false
            }
        };
        self.emit(block)?;
        // A zoom block is neither a section nor an encoding of one, so the
        // counters a caller reads as `section_count` must not see it.
        if !zooming {
            if let Some(encoding) = encoding {
                self.section_counts.bump(encoding);
            }
            self.section_count += 1;
        }
        Ok(())
    }

    /// Place the oldest block in flight, waiting on its deflate if it is still
    /// running.
    ///
    /// Where the pipeline is put back in order. Blocks are compressed in
    /// whatever order the workers finish them and placed in the order they were
    /// submitted, so the offsets, the leaves of the R-tree and the bytes of the
    /// file all agree and the result is the file the serial path would have
    /// written, byte for byte.
    fn commit_one(&mut self) -> Result<()> {
        let Some(item) = self.pending.pop_front() else {
            return Ok(());
        };
        let block = item.result.take()?;
        self.place_block(item.leaf, &block, item.encoding)
    }

    fn drain(&mut self) -> Result<()> {
        while !self.pending.is_empty() {
            self.commit_one()?;
        }
        Ok(())
    }

    /// Compress a block and index it, on a worker thread if there is one.
    ///
    /// Submission blocks once enough blocks are in flight, which is what bounds
    /// the memory the pipeline holds. A producer faster than the workers then
    /// runs at their pace, which is the pace the serial path ran at anyway.
    fn submit_block(
        &mut self,
        leaf: LeafItem,
        body: Vec<u8>,
        encoding: Option<WigEncoding>,
    ) -> Result<()> {
        self.uncompress_buffer_size = self.uncompress_buffer_size.max(body.len() as u64);
        let level = self.compression_level;
        let Some(executor) = &self.executor else {
            let block = Self::pack_block(body, level, &self.path)?;
            return self.place_block(leaf, &block, encoding);
        };
        let promise: Promise<Result<Vec<u8>>> = Promise::new();
        let path = self.path.clone();
        let worker = promise.clone();
        executor.spawn(move || worker.set(Self::pack_block(body, level, &path)));
        self.pending.push_back(PendingBlock {
            leaf,
            result: PendingResult::InFlight(promise),
            encoding,
        });
        while self.pending.len() >= self.pending_limit {
            self.commit_one()?;
        }
        Ok(())
    }

    // -- chromosomes and validation ----------------------------------------

    /// The index of `chr_id`, assigning it one on first sight.
    fn resolve_chr(&mut self, chr_id: &str) -> Result<u32> {
        if let Some(index) = self.last_chr_index {
            if chr_id == self.last_chr_id {
                return Ok(index);
            }
        }
        let mut name = chr_id.to_string();
        let mut declared = None;
        if let Some(sizes) = &self.declared_sizes {
            let entry = sizes.resolve(chr_id)?;
            name = entry.id.clone();
            declared = Some(entry.size);
        }
        if let Some(state) = self.chrs.get(&name) {
            if Some(state.index) != self.last_chr_index {
                return Err(Error::invalid(format!(
                    "chromosome {name} was already written, values must be pooled by chromosome"
                )));
            }
            self.last_chr_id = chr_id.to_string();
            return Ok(state.index);
        }
        // A wig section carries one chromosome id in its header, so it has to
        // end here. A bed block does not: every record of it carries its own,
        // which is what lets a genome of many small scaffolds pack them into
        // one block rather than spend a block and an index leaf on each.
        if self.is_bigwig() {
            self.flush_section()?;
        }
        let index = to_bbi_u32(self.chrs.len() as i64, "chromId")?;
        self.chrs.insert(
            name.clone(),
            ChrState {
                index,
                size: declared.unwrap_or(0),
                declared: declared.is_some(),
            },
        );
        self.current_chr = Some(name);
        self.last_chr_id = chr_id.to_string();
        self.last_chr_index = Some(index);
        self.last_chr_end = 0;
        self.last_entry_start = -1;
        Ok(index)
    }

    fn chr_state(&self) -> ChrState {
        let name = self.current_chr.as_deref().unwrap_or_default();
        self.chrs[name]
    }

    fn grow_chr(&mut self, end: i64) {
        if let Some(name) = &self.current_chr {
            let state = self.chrs.get_mut(name).expect("current chromosome exists");
            if !state.declared {
                state.size = state.size.max(end);
            }
        }
    }

    /// The end a value may actually claim on the chromosome it sits on: its
    /// own, or the end of the chromosome when it hangs over it.
    ///
    /// Values wider than one base come off a grid, and a chromosome is rarely a
    /// whole number of bins long, so the last one of a chromosome read at a bin
    /// size of ten overshoots by up to nine bases. That is arithmetic rather
    /// than an error, and it is cut back here. Starting at or past the end is a
    /// different thing — the value belongs to no base of it — and raises.
    ///
    /// A value cut back here leaves the chromosome at its end, so nothing can
    /// follow it there: anything starting inside overlaps it and anything
    /// starting past it raises. A cut value is therefore always the last of its
    /// section, which is what lets a section carry one that no longer matches
    /// its uniform span.
    fn clip_to_chr(&mut self, start: i64, end: i64) -> Result<i64> {
        let state = self.chr_state();
        if !state.declared || end <= state.size {
            return Ok(end);
        }
        if start >= state.size {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} starts past the end of {}, which is {} bases long",
                self.last_chr_id, self.last_chr_id, state.size
            )));
        }
        self.clipped_count += 1;
        Ok(state.size)
    }

    /// Check a range against the ordering contract and the chromosome it sits
    /// on, and take it as the new high-water mark.
    ///
    /// `last_start` is the start of the last value the range covers, which is
    /// `start` for a single one. Only that value can be cut back by the end of
    /// a chromosome, so it is the one the chromosome is read against: a run
    /// reaching whole values past the end starts one of them past it and
    /// raises, where a run merely hanging over the end does not.
    fn validate_range(&mut self, start: i64, end: i64, last_start: i64) -> Result<i64> {
        if start < 0 {
            return Err(Error::invalid(format!(
                "start {start} must not be negative"
            )));
        }
        // `end == start` is allowed: BED 1.0 permits it and an insertion is
        // written that way, so refusing it made `convert_to_bigbed` fail on
        // files the format calls valid. It covers no base, which the coverage
        // sweep already handles — the run it would open closes at the position
        // it opened at — and `EntryWalk::read` already goes out of its way to
        // find such an entry when reading one back. bigWig values are a
        // different matter and stay strictly positive-width.
        if end < start {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} ends before it starts",
                self.last_chr_id
            )));
        }
        if start < self.last_chr_end {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} starts before the end {} of the previous value, values \
                 must be added in order and without overlap",
                self.last_chr_id, self.last_chr_end
            )));
        }
        let end = self.clip_to_chr(last_start, end)?;
        to_bbi_u32(end, "coordinate")?;
        self.last_chr_end = end;
        Ok(end)
    }

    /// Check a bed entry against the ordering contract and its chromosome.
    ///
    /// Only the starts have to be in order, unlike the values of a bigWig. A
    /// bed of anything real — genes, repeats, peaks called twice over — has
    /// entries that overlap and nest, and the format asks only that they be
    /// sorted by (chromosome, start).
    fn validate_entry(&mut self, start: i64, end: i64) -> Result<()> {
        if start < 0 {
            return Err(Error::invalid(format!(
                "start {start} must not be negative"
            )));
        }
        // `end == start` is allowed: BED 1.0 permits it and an insertion is
        // written that way, so refusing it made `convert_to_bigbed` fail on
        // files the format calls valid. It covers no base, which the coverage
        // sweep already handles — the run it would open closes at the position
        // it opened at — and `EntryWalk::read` already goes out of its way to
        // find such an entry when reading one back. bigWig values are a
        // different matter and stay strictly positive-width.
        if end < start {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} ends before it starts",
                self.last_chr_id
            )));
        }
        if start < self.last_entry_start {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} starts before the previous entry at {}, entries must be \
                 added in order of their start",
                self.last_chr_id, self.last_entry_start
            )));
        }
        let state = self.chr_state();
        if state.declared && end > state.size {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} runs past the end of {}, which is {} bases long",
                self.last_chr_id, self.last_chr_id, state.size
            )));
        }
        to_bbi_u32(end, "coordinate")?;
        self.last_entry_start = start;
        Ok(())
    }

    // -- statistics --------------------------------------------------------

    /// Add one value covering `span` bases to the whole-file summary.
    ///
    /// The extremes are seeded from the first value rather than compared
    /// against the NaN the struct defaults them to, since a comparison against
    /// NaN is false either way and would leave both of them NaN for ever.
    fn account_value(&mut self, value: f32, span: i64) {
        if self.summary.bases_covered == 0 {
            self.summary.min_value = value as f64;
            self.summary.max_value = value as f64;
        } else {
            if (value as f64) < self.summary.min_value {
                self.summary.min_value = value as f64;
            }
            if (value as f64) > self.summary.max_value {
                self.summary.max_value = value as f64;
            }
        }
        self.summary.bases_covered += span as u64;
        self.summary.sum_data += value as f64 * span as f64;
        self.summary.sum_squared += value as f64 * value as f64 * span as f64;
        self.item_count += 1;
        if !self.ladder_frozen && self.item_count >= RESOLUTION_SAMPLE_ITEMS {
            self.freeze_zoom_ladder();
        }
    }

    /// Fix the reduction ladder from the mean item span seen so far.
    ///
    /// UCSC anchors the first level at ten times the mean span, then
    /// quadruples, taking the mean from the whole file it wrote to a temporary
    /// one. Here the candidates are counted as the data goes by, so the mean
    /// comes from the first few thousand items. The ladder only decides how
    /// coarse the summaries are, so a poor estimate makes the file bigger or
    /// its zoomed queries wider, never wrong.
    fn freeze_zoom_ladder(&mut self) {
        let mean_span = self
            .summary
            .bases_covered
            .checked_div(self.item_count)
            .unwrap_or(1)
            .max(1) as i64;
        let mut reduction = (mean_span * INITIAL_ZOOM_FACTOR).max(1);
        for level in 0..MAX_ZOOM_LEVELS {
            self.zoom_reductions[level] = reduction;
            self.zoom_res_ends[level] = 0;
            // Capped rather than wrapped: a reduction level is a u32, and past
            // the widest chromosome a coarser one summarises nothing new.
            if reduction > 0xFFFF_FFFF / ZOOM_INCREMENT {
                continue;
            }
            reduction *= ZOOM_INCREMENT;
        }
        self.zoom_res_chr_index = None;
        self.ladder_start_item_count = self.item_count;
        self.ladder_frozen = true;
    }

    /// Count the windows each candidate reduction would open over
    /// `[start, end)`.
    ///
    /// Windows are anchored to the data rather than aligned to a grid, the way
    /// the zoom pass itself anchors them, so the counts are what that pass will
    /// actually produce. Contiguous items give the same answer whether they are
    /// counted one at a time or as a single range, which is what lets a run of
    /// values be counted in one call.
    fn count_resolutions(&mut self, chr_index: u32, start: i64, end: i64) {
        if !self.ladder_frozen {
            return;
        }
        if self.zoom_res_chr_index != Some(chr_index) {
            self.zoom_res_ends = [0; MAX_ZOOM_LEVELS];
            self.zoom_res_chr_index = Some(chr_index);
        }
        for level in 0..MAX_ZOOM_LEVELS {
            let reduction = self.zoom_reductions[level];
            if start >= self.zoom_res_ends[level] {
                self.zoom_res_sizes[level] += 1;
                self.zoom_res_ends[level] = start + reduction;
            }
            if end > self.zoom_res_ends[level] {
                // Closed form rather than a loop: one bedGraph interval can be
                // megabases wide against a reduction of a few dozen bases.
                let extra = (end - self.zoom_res_ends[level] + reduction - 1) / reduction;
                self.zoom_res_sizes[level] += extra as u64;
                self.zoom_res_ends[level] += extra * reduction;
            }
        }
    }

    // -- section accumulation ----------------------------------------------

    /// Put one validated value over `[start, end)` into the open section, and
    /// into everything the file has to say about it afterwards.
    fn place_value(
        &mut self,
        chr_index: u32,
        start: i64,
        end: i64,
        value: f32,
        batch_remaining: usize,
    ) -> Result<()> {
        if !value.is_finite() {
            self.skipped_count += 1;
            return Ok(());
        }
        self.grow_chr(end);
        self.account_value(value, end - start);
        self.count_resolutions(chr_index, start, end);
        self.add_item(chr_index, start, end - start, value, batch_remaining)
    }

    /// Add one item to the open section, opening, widening or closing it as the
    /// item requires.
    fn add_item(
        &mut self,
        chr_index: u32,
        start: i64,
        span: i64,
        value: f32,
        batch_remaining: usize,
    ) -> Result<()> {
        if self
            .section
            .offer(chr_index, start, span, value, batch_remaining)
            == Accept::Flush
        {
            self.flush_section()?;
            // Offered again, to the empty section — which cannot refuse it. The
            // retry is the whole point of `Flush`, so it is a statement and not
            // an assertion: inside a `debug_assert` it would be compiled out of
            // a release build and every item that split a section would be
            // silently dropped.
            let retried = self
                .section
                .offer(chr_index, start, span, value, batch_remaining);
            debug_assert_eq!(retried, Accept::Buffered);
        }
        if self.section.is_full() {
            self.flush_section()?;
        }
        Ok(())
    }

    /// Encode, compress and write the open section, and index the block it
    /// became.
    fn flush_section(&mut self) -> Result<()> {
        if self.section.is_empty() {
            return Ok(());
        }
        let body = self.section.encode()?;
        self.data_body_size += body.len() as u64;
        let leaf = LeafItem {
            start_chr: self.section.chr_ix,
            start_base: to_bbi_u32(self.section.first_start, "chromStart")?,
            end_chr: self.section.chr_ix,
            end_base: to_bbi_u32(self.section.last_end, "chromEnd")?,
            offset: 0,
            size: 0,
        };
        let encoding = self.section.encoding();
        self.section.clear();
        self.submit_block(leaf, body, Some(encoding))
    }

    // -- bed records -------------------------------------------------------

    /// Append one record to the block being filled (Supp. Table 12).
    ///
    /// A declared field the caller leaves out is written empty rather than
    /// dropped. Every record of a bed carries every column, so a missing one
    /// would make the record unreadable rather than merely incomplete.
    fn append_bed_record(
        &mut self,
        chr_index: u32,
        start: i64,
        end: i64,
        values: &IndexMap<String, String>,
    ) -> Result<()> {
        for name in values.keys() {
            if !self.bed_fields.contains_key(name) {
                return Err(Error::invalid(format!(
                    "field {name} is not one this file declares"
                )));
            }
            // The coordinates are the record's own three fields and are given
            // as arguments, so naming one here is a mistake worth reporting
            // rather than a value that would silently go nowhere.
            if self.bed_fields.get_index_of(name).is_some_and(|i| i < 3) {
                return Err(Error::invalid(format!(
                    "field {name} is a coordinate, which is written from start and end"
                )));
            }
        }
        let (start_u32, end_u32) = (
            to_bbi_u32(start, "chromStart")?,
            to_bbi_u32(end, "chromEnd")?,
        );
        self.bed_bounds = Some(match self.bed_bounds {
            None => (chr_index, start_u32, chr_index, end_u32),
            Some((sc, sb, ec, _)) if chr_index != ec => (sc, sb, chr_index, end_u32),
            // Entries may nest, so the furthest one reaches is not the last one
            // to start.
            Some((sc, sb, ec, eb)) => (sc, sb, ec, eb.max(end_u32)),
        });

        self.bed_block.extend_from_slice(&chr_index.to_le_bytes());
        self.bed_block.extend_from_slice(&start_u32.to_le_bytes());
        self.bed_block.extend_from_slice(&end_u32.to_le_bytes());
        for (index, (name, kind)) in self.bed_fields.iter().enumerate() {
            if index < 3 {
                continue;
            }
            if index > 3 {
                self.bed_block.push(b'\t');
            }
            match values.get(name) {
                Some(text) => {
                    if text.contains('\t') || text.contains('\0') {
                        return Err(Error::invalid(format!(
                            "field {name} value {text} contains a tab or a null byte"
                        )));
                    }
                    self.bed_block.extend_from_slice(text.as_bytes());
                }
                // A numeric column left out is written as a zero: a bed reader
                // parses it, where an empty one is not a number at all.
                None if kind != "string" => self.bed_block.push(b'0'),
                None => {}
            }
        }
        self.bed_block.push(0);
        self.bed_block_count += 1;
        self.entry_count += 1;
        Ok(())
    }

    /// Compress and write the block of bed entries being filled.
    fn flush_bed_block(&mut self) -> Result<()> {
        if self.bed_block_count == 0 {
            return Ok(());
        }
        let body = std::mem::take(&mut self.bed_block);
        // The next block starts from a fresh buffer reserved to what the last
        // one came to, rather than grown into over a hundred entries.
        self.bed_block = Vec::with_capacity(body.len());
        self.data_body_size += body.len() as u64;
        let (sc, sb, ec, eb) = self.bed_bounds.take().expect("a filled block has bounds");
        let leaf = LeafItem {
            start_chr: sc,
            start_base: sb,
            end_chr: ec,
            end_base: eb,
            offset: 0,
            size: 0,
        };
        self.bed_block_count = 0;
        self.submit_block(leaf, body, None)
    }

    /// End whatever is being filled, whichever kind of file this is.
    fn flush_pending(&mut self) -> Result<()> {
        if self.is_bigwig() {
            self.flush_section()?;
        } else {
            self.flush_bed_block()?;
        }
        self.drain()
    }

    // -- the public write API ----------------------------------------------

    /// Write one value over `[start, end)`.
    ///
    /// A value that overruns the end of a declared chromosome having started
    /// inside it is written up to that end rather than refused, since a
    /// chromosome is rarely a whole number of bins long; `clipped_count` counts
    /// them. A value that is not finite is not written at all — it leaves a
    /// gap, which is what a bigWig means by a base carrying no data.
    pub fn write_value(&mut self, chr: &str, start: i64, end: i64, value: f32) -> Result<()> {
        self.check_open()?;
        if !self.is_bigwig() {
            return Err(Error::invalid(
                "write_value is only for bigwig files, use write_entry",
            ));
        }
        self.failed = true;
        let result = (|| {
            let chr_index = self.resolve_chr(chr)?;
            let end = self.validate_range(start, end, start)?;
            self.place_value(chr_index, start, end, value, 0)
        })();
        self.failed = result.is_err();
        result
    }

    /// Write `values.len()` values of `span` bases each, the first at `start`
    /// and each one starting where the one before it ended.
    ///
    /// The fast path. Values that come in this way are a fixedStep section by
    /// construction — four bytes an item against twelve — so a run handed over
    /// in one call costs a third of what the same values cost one at a time. A
    /// run with no gaps also goes in as a single extend: the starts are implied
    /// by the step and never materialised.
    ///
    /// Only the last value of a run can overrun the end of a declared
    /// chromosome having started inside it, since the ones before it are a
    /// whole span short of it. It is cut back to that end and, no longer being
    /// a step of `span`, laid down on its own rather than in the extend.
    pub fn write_values(&mut self, chr: &str, start: i64, span: i64, values: &[f32]) -> Result<()> {
        self.check_open()?;
        if !self.is_bigwig() {
            return Err(Error::invalid(
                "write_values is only for bigwig files, use write_entry",
            ));
        }
        if values.is_empty() {
            return Ok(());
        }
        self.failed = true;
        let result = self.write_values_inner(chr, start, span, values);
        self.failed = result.is_err();
        result
    }

    fn write_values_inner(
        &mut self,
        chr: &str,
        start: i64,
        span: i64,
        values: &[f32],
    ) -> Result<()> {
        if span <= 0 {
            return Err(Error::invalid(format!("span {span} must be positive")));
        }
        let count = values.len() as i64;
        let chr_index = self.resolve_chr(chr)?;
        let last_start = start + span * (count - 1);
        let last_end = self.validate_range(start, last_start + span, last_start)?;
        // A value cut short is not a step of span, so it cannot ride the run.
        let run_count = if last_end - last_start == span {
            count
        } else {
            count - 1
        };

        let mut index = 0i64;
        while index < run_count {
            let item_start = start + span * index;
            // A section that has widened cannot narrow again, so a long
            // uniform run fed into one is written at the wider item size for
            // no reason. Close it first and let the run have a fixedStep
            // section of its own — see `WigSection::should_flush_for_run`.
            if self.section.should_flush_for_run(
                chr_index,
                item_start,
                span,
                (run_count - index) as usize,
            ) {
                self.flush_section()?;
            }
            let extends = self.section.extends_run(chr_index, item_start, span);
            if extends || self.section.is_empty() {
                // Only as far as the section can take, not as far as the run
                // goes: scanning the whole run and keeping a slotful rescans
                // the remainder on the next pass, which is quadratic in the
                // size of the call.
                let limit =
                    (run_count - index).min(self.items_per_slot as i64 - self.section.len() as i64);
                let mut run = 0i64;
                while run < limit && values[(index + run) as usize].is_finite() {
                    run += 1;
                }
                if run > 0 {
                    let slice = &values[index as usize..(index + run) as usize];
                    self.section.extend_run(chr_index, item_start, span, slice);
                    let section_end = start + span * (index + run);
                    // Once for the run, not once per value: `grow_chr` takes
                    // the maximum, so the run's own end is the only call that
                    // can change anything — and each one is a hash lookup on
                    // the chromosome name.
                    self.grow_chr(section_end);
                    for i in 0..run {
                        self.account_value(values[(index + i) as usize], span);
                    }
                    self.count_resolutions(chr_index, item_start, section_end);
                    index += run;
                    if self.section.is_full() {
                        self.flush_section()?;
                    }
                    continue;
                }
            }
            self.place_value(
                chr_index,
                item_start,
                item_start + span,
                values[index as usize],
                (count - index - 1) as usize,
            )?;
            index += 1;
        }
        if run_count != count {
            self.place_value(
                chr_index,
                last_start,
                last_end,
                values[(count - 1) as usize],
                0,
            )?;
        }
        Ok(())
    }

    /// Write one bed entry over `[start, end)`.
    ///
    /// `values` are the fields past the coordinates, already formatted as the
    /// text a bed stores, under the names the file's own fields declare. A
    /// declared field left out is written empty.
    ///
    /// Entries may overlap and nest, unlike the values of a bigWig. Only their
    /// starts have to be in order.
    pub fn write_entry(
        &mut self,
        chr: &str,
        start: i64,
        end: i64,
        values: &IndexMap<String, String>,
    ) -> Result<()> {
        self.check_open()?;
        if self.is_bigwig() {
            return Err(Error::invalid(
                "write_entry is only for bigbed files, use write_value",
            ));
        }
        self.failed = true;
        let result = (|| {
            let chr_index = self.resolve_chr(chr)?;
            self.validate_entry(start, end)?;
            self.grow_chr(end);
            self.append_bed_record(chr_index, start, end, values)?;
            // The sweep's runs are what the summary and the zoom levels of a
            // bigBed hold, so they are accumulated as the entries arrive.
            let runs = self.coverage.add(chr_index, start, end);
            for (chr, s, e, depth) in runs {
                self.account_value(depth, e - s);
                self.count_resolutions(chr, s, e);
            }
            if self.bed_block_count >= self.items_per_slot {
                self.flush_bed_block()?;
            }
            Ok(())
        })();
        self.failed = result.is_err();
        result
    }

    fn check_open(&self) -> Result<()> {
        if self.closed {
            return Err(Error::invalid(format!(
                "error writing to {} (file is closed)",
                self.path
            )));
        }
        Ok(())
    }

    // -- finalisation ------------------------------------------------------

    /// Finish the file: the index, the zoom levels, the chromosome tree and the
    /// headers, in that order.
    ///
    /// Calling it twice is harmless. Until it returns, nothing on disk is a
    /// bigWig or bigBed.
    pub fn close(&mut self) -> Result<()> {
        if self.closed {
            return Ok(());
        }
        self.failed = true;
        let result = self.finish();
        // The deflate threads stop here rather than at drop, so a closed file
        // holds none of this library's — what a caller closing one, or leaving
        // its with-block, expects. Also on the way out of a failure.
        self.executor = None;
        self.pending.clear();
        result?;
        if let Some(sink) = &mut self.sink {
            sink.close()?;
        }
        self.sink = None;
        self.closed = true;
        self.failed = false;
        Ok(())
    }

    fn finish(&mut self) -> Result<()> {
        // The sweep holds the runs its last entries left open, and they belong
        // to the summary, so it has to be finished before the block is.
        if !self.is_bigwig() {
            for (chr, s, e, depth) in self.coverage.finish() {
                self.account_value(depth, e - s);
                self.count_resolutions(chr, s, e);
            }
        }
        self.flush_pending()?;

        self.full_index_offset = self.sync_cursor()?;
        self.write_data_tree()?;
        self.write_zoom_levels()?;
        self.write_chromosome_tree()?;

        // The file's own magic repeated at the end (Supp. Table 5), as UCSC's
        // writers close a file, marking it untruncated. It has to match the
        // type: no reader here checks it, but one that does would call every
        // bigBed corrupt.
        let magic = self.magic();
        self.emit(&magic.to_le_bytes())?;
        self.sink()?.flush()?;
        self.write_headers()
    }

    fn magic(&self) -> u32 {
        match self.kind {
            BbiKind::BigWig => super::BIGWIG_MAGIC,
            BbiKind::BigBed => super::BIGBED_MAGIC,
        }
    }

    fn write_data_tree(&mut self) -> Result<()> {
        let items = std::mem::take(&mut self.data_items);
        let mut bytes = Vec::new();
        super::rtree::write_tree(
            &items,
            self.full_index_offset,
            self.block_size,
            self.items_per_slot as u32,
            self.full_index_offset,
            &mut |b| bytes.extend_from_slice(b),
        )?;
        self.data_items = items;
        self.emit(&bytes)
    }

    fn write_chromosome_tree(&mut self) -> Result<()> {
        let mut entries: Vec<WriteEntry> = self
            .chrs
            .iter()
            .map(|(id, state)| {
                Ok(WriteEntry {
                    id: id.clone(),
                    size: to_bbi_u32(state.size.max(1), "chromSize")?,
                    index: state.index,
                })
            })
            .collect::<Result<_>>()?;
        // The tree is searched by name, so its leaves are ordered by name. The
        // ids they carry follow the order the chromosomes were written in,
        // which is what keeps the data sorted by (chromosome, start), and the
        // format has never required the two orders to agree.
        entries.sort_by(|a, b| a.id.cmp(&b.id));
        self.chr_tree_offset = self.sync_cursor()?;
        let mut bytes = Vec::new();
        super::chr_tree::write_tree(&entries, self.chr_tree_offset, self.block_size, &mut |b| {
            bytes.extend_from_slice(b)
        })?;
        self.emit(&bytes)
    }

    fn write_headers(&mut self) -> Result<()> {
        let mut zoom_bytes = Vec::new();
        for header in &self.zoom_headers {
            zoom_bytes.extend_from_slice(&write_zoom_header(header));
        }
        if !zoom_bytes.is_empty() {
            self.patch(BBI_HEADER_SIZE, &zoom_bytes)?;
        }

        let summary = write_total_summary(&self.summary);
        self.patch(self.total_summary_offset, &summary)?;

        // A bigWig counts the sections it holds, a bigBed the entries.
        let count = if self.is_bigwig() {
            self.section_count
        } else {
            self.entry_count
        };
        self.patch(self.full_data_offset, &count.to_le_bytes())?;

        let header = BbiHeader {
            kind: self.kind,
            version: BBI_OUTPUT_VERSION,
            zoom_levels: self.zoom_headers.len() as u16,
            chr_tree_offset: self.chr_tree_offset,
            full_data_offset: self.full_data_offset,
            full_index_offset: self.full_index_offset,
            field_count: self.field_count,
            defined_field_count: self.defined_field_count,
            auto_sql_offset: self.auto_sql_offset,
            total_summary_offset: self.total_summary_offset,
            uncompress_buffer_size: if self.compression_level > 0 {
                to_bbi_u32(self.uncompress_buffer_size as i64, "uncompressBufSize")?
            } else {
                0
            },
        };
        let bytes = write_header(&header)?;

        // Everything but the magic, then the magic alone. The file is not a bbi
        // file until those four bytes land, so a write dying between the two
        // leaves something no reader accepts rather than a header pointing at
        // nothing.
        self.patch(4, &bytes[4..])?;
        let magic = self.magic();
        self.patch(0, &magic.to_le_bytes())
    }

    // -- accessors ---------------------------------------------------------

    pub fn path(&self) -> &str {
        &self.path
    }
    pub fn kind(&self) -> BbiKind {
        self.kind
    }
    pub fn fields(&self) -> &IndexMap<String, String> {
        &self.bed_fields
    }
    pub fn section_counts(&self) -> SectionCounts {
        self.section_counts
    }
    pub fn section_count(&self) -> u64 {
        self.section_count
    }
    pub fn entry_count(&self) -> u64 {
        self.entry_count
    }
    pub fn skipped_count(&self) -> u64 {
        self.skipped_count
    }
    pub fn clipped_count(&self) -> u64 {
        self.clipped_count
    }
    pub fn is_closed(&self) -> bool {
        self.closed
    }
    /// Chromosome sizes as written, in the order the chromosomes were written.
    pub fn chr_sizes(&self) -> Vec<(String, i64)> {
        self.chrs
            .iter()
            .map(|(id, state)| (id.clone(), state.size))
            .collect()
    }

    /// Give up on the file, stopping the workers and **removing what was
    /// written**.
    ///
    /// What a caller that will never close this writer says, so `Drop` leaves
    /// the file alone rather than finishing it. Removing it rather than leaving
    /// it is what the converters promise: `LocalSink::create` truncated the
    /// path at open, so what would otherwise stay behind is a zero-magic stub —
    /// no reader accepts it, which is right, but the caller still has to clean
    /// it up, and an interrupted conversion that says it left no file must not
    /// leave one. Harmless on a writer already closed, which owns a finished
    /// file rather than an abandoned one.
    pub fn abandon(&mut self) {
        if self.closed {
            return;
        }
        self.failed = true;
        self.pending.clear();
        self.executor = None;
        if let Some(mut sink) = self.sink.take() {
            sink.discard();
        }
    }
}

impl Drop for BbiWriter {
    fn drop(&mut self) {
        // A writer that already threw is left alone: finishing it would turn
        // the caller's error into a file that looks complete and is not.
        if self.failed || self.closed {
            return;
        }
        let _ = self.close();
    }
}

impl std::fmt::Debug for BbiWriter {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("BbiWriter")
            .field("path", &self.path)
            .field("type", &self.kind.as_str())
            .field("closed", &self.closed)
            .finish()
    }
}

/// Rolling coverage depth over sorted, overlapping entries — what a bigBed's
/// summary and zoom levels describe.
///
/// Supp. Table 5 puts it as "the values correspond to those of a BigWig
/// constructed by the depth of coverage of bases", so a base under three
/// entries counts once towards the bases covered and three towards the sum, and
/// a base under none counts towards neither.
///
/// Entries arrive sorted by start but may overlap and nest freely, which is the
/// ordinary shape of a bed — genes overlap. A heap of the ends still open is
/// therefore what says how deep the current base is, and it holds one entry per
/// unit of depth rather than one per entry of the file.
#[derive(Debug, Default)]
pub(crate) struct CoverageSweep {
    chr_index: Option<u32>,
    position: i64,
    /// A min-heap of the ends still open, which `Reverse` is what makes of
    /// Rust's max-heap.
    ends: BinaryHeap<std::cmp::Reverse<i64>>,
}

/// One run of constant depth: chromosome, start, end, depth.
type Run = (u32, i64, i64, f32);

impl CoverageSweep {
    /// Add one entry, returning every run of constant depth it closes.
    ///
    /// A run is only emitted once nothing can change it, which is why the entry
    /// that opens it is not the entry that reports it, and why `finish` has to
    /// be called before the numbers mean anything.
    fn add(&mut self, chr: u32, start: i64, end: i64) -> Vec<Run> {
        let mut out = Vec::new();
        if self.chr_index != Some(chr) {
            out.extend(self.finish());
            self.chr_index = Some(chr);
            self.position = start;
        }
        while self.ends.peek().is_some_and(|e| e.0 <= start) {
            let expiry = self.ends.peek().expect("just peeked").0;
            if expiry > self.position {
                out.push((
                    self.chr_index.expect("a run has a chromosome"),
                    self.position,
                    expiry,
                    self.ends.len() as f32,
                ));
                self.position = expiry;
            }
            while self.ends.peek().is_some_and(|e| e.0 == expiry) {
                self.ends.pop();
            }
        }
        // Whatever is still open reaches past this entry's start, so the depth
        // over what is left before it is constant.
        if !self.ends.is_empty() && start > self.position {
            out.push((
                self.chr_index.expect("a run has a chromosome"),
                self.position,
                start,
                self.ends.len() as f32,
            ));
        }
        if start > self.position {
            self.position = start;
        }
        self.ends.push(std::cmp::Reverse(end));
        out
    }

    /// Close every run still open, at the end of the chromosome or of the file.
    fn finish(&mut self) -> Vec<Run> {
        let mut out = Vec::new();
        while let Some(std::cmp::Reverse(expiry)) = self.ends.peek().copied() {
            if expiry > self.position {
                out.push((
                    self.chr_index.expect("a run has a chromosome"),
                    self.position,
                    expiry,
                    self.ends.len() as f32,
                ));
                self.position = expiry;
            }
            while self.ends.peek().is_some_and(|e| e.0 == expiry) {
                self.ends.pop();
            }
        }
        self.chr_index = None;
        self.position = 0;
        out
    }
}

// -- the zoom pass ---------------------------------------------------------

/// The zoom window being accumulated, holding the statistics of Supp. Table 19
/// with the two sums widened.
///
/// The sums are `f64` even though the record stores them as `f32`, because a
/// window of a few million bases passes the point where an `f32` accumulator
/// stops making progress long before it is closed.
#[derive(Debug, Default)]
struct ZoomWindow {
    open: bool,
    chr_index: u32,
    start: i64,
    end: i64,
    limit: i64,
    valid_count: u64,
    min_value: f32,
    max_value: f32,
    sum_data: f64,
    sum_squared: f64,
}

/// One zoom record (Supp. Table 19).
#[derive(Debug, Clone, Copy)]
struct ZoomRecord {
    chr_index: u32,
    chr_start: u32,
    chr_end: u32,
    valid_count: u32,
    min_value: f32,
    max_value: f32,
    sum_data: f32,
    sum_squared: f32,
}

impl ZoomRecord {
    fn write(&self, out: &mut Vec<u8>) {
        out.extend_from_slice(&self.chr_index.to_le_bytes());
        out.extend_from_slice(&self.chr_start.to_le_bytes());
        out.extend_from_slice(&self.chr_end.to_le_bytes());
        out.extend_from_slice(&self.valid_count.to_le_bytes());
        out.extend_from_slice(&self.min_value.to_le_bytes());
        out.extend_from_slice(&self.max_value.to_le_bytes());
        out.extend_from_slice(&self.sum_data.to_le_bytes());
        out.extend_from_slice(&self.sum_squared.to_le_bytes());
    }

    fn read(block: &[u8], offset: usize) -> Self {
        let u32_at = |o: usize| u32::from_le_bytes(block[o..o + 4].try_into().expect("4 bytes"));
        let f32_at = |o: usize| f32::from_le_bytes(block[o..o + 4].try_into().expect("4 bytes"));
        Self {
            chr_index: u32_at(offset),
            chr_start: u32_at(offset + 4),
            chr_end: u32_at(offset + 8),
            valid_count: u32_at(offset + 12),
            min_value: f32_at(offset + 16),
            max_value: f32_at(offset + 20),
            sum_data: f32_at(offset + 24),
            sum_squared: f32_at(offset + 28),
        }
    }
}

/// Accumulates zoom records at one reduction, closing a window whenever the
/// data leaves it, and grouping the closed records into blocks.
///
/// The blocks a level produces go out through the writer's own pipeline — the
/// same compression on the same pool, placed in submission order — but they are
/// indexed by the level rather than by the data section, and count towards
/// neither `section_count` nor the encoding counters. `BbiWriter::zoom_items`
/// is the switch.
struct ZoomLevelBuilder {
    reduction: i64,
    items_per_slot: usize,
    window: ZoomWindow,
    records: Vec<ZoomRecord>,
    /// Encoded, uncompressed blocks and their leaves, in submission order.
    ///
    /// A staging queue, not an accumulator: the writer drains it into its own
    /// pipeline after every call that can close a block, so what sits here is
    /// what one source block just produced rather than a whole level. Holding
    /// a level was the difference between a bounded `close()` and one resident
    /// gigabyte per gigabyte of data body.
    blocks: Vec<(LeafItem, Vec<u8>)>,
    record_count: u64,
}

impl ZoomLevelBuilder {
    fn new(reduction: i64, items_per_slot: usize) -> Self {
        Self {
            reduction,
            items_per_slot,
            window: ZoomWindow::default(),
            records: Vec::new(),
            blocks: Vec::new(),
            record_count: 0,
        }
    }

    fn open_window(&mut self, chr_index: u32, start: i64, min_value: f32, max_value: f32) {
        self.window = ZoomWindow {
            open: true,
            chr_index,
            start,
            end: start,
            limit: start + self.reduction,
            valid_count: 0,
            min_value,
            max_value,
            sum_data: 0.0,
            sum_squared: 0.0,
        };
    }

    fn close_window(&mut self) -> Result<()> {
        if !self.window.open {
            return Ok(());
        }
        self.window.open = false;
        if self.window.valid_count == 0 {
            return Ok(());
        }
        self.records.push(ZoomRecord {
            chr_index: self.window.chr_index,
            chr_start: to_bbi_u32(self.window.start, "chromStart")?,
            chr_end: to_bbi_u32(self.window.end, "chromEnd")?,
            valid_count: to_bbi_u32(self.window.valid_count as i64, "validCount")?,
            min_value: self.window.min_value,
            max_value: self.window.max_value,
            sum_data: to_bbi_f32(self.window.sum_data),
            sum_squared: to_bbi_f32(self.window.sum_squared),
        });
        self.record_count += 1;
        if self.records.len() >= self.items_per_slot {
            self.flush_block()?;
        }
        Ok(())
    }

    /// Add a full-resolution interval, splitting it across as many windows as
    /// it spans.
    ///
    /// An interval of a bedGraph can be megabases wide, so at the first level
    /// it has to be cut rather than summarised whole. Each piece contributes
    /// the bases it actually covers, which is what makes `validCount` the count
    /// of bases carrying data and not the width of the window — the distinction
    /// every summed or counted query depends on.
    fn add_interval(&mut self, chr_index: u32, mut start: i64, end: i64, value: f32) -> Result<()> {
        while start < end {
            if !self.window.open || self.window.chr_index != chr_index || start >= self.window.limit
            {
                self.close_window()?;
                self.open_window(chr_index, start, value, value);
            }
            let part_end = end.min(self.window.limit);
            let overlap = part_end - start;
            self.window.valid_count += overlap as u64;
            self.window.sum_data += value as f64 * overlap as f64;
            self.window.sum_squared += value as f64 * value as f64 * overlap as f64;
            if value < self.window.min_value {
                self.window.min_value = value;
            }
            if value > self.window.max_value {
                self.window.max_value = value;
            }
            self.window.end = part_end;
            start = part_end;
            if start >= self.window.limit {
                self.close_window()?;
            }
        }
        Ok(())
    }

    /// Merge a record of the level below, whole.
    ///
    /// Unlike an interval, an already summarised record is never cut. Splitting
    /// one would mean guessing how its mass sits inside it, so the coarse
    /// window is allowed to overrun instead, which is what UCSC does and what
    /// keeps the sums exactly conserved across levels.
    fn add_record(&mut self, record: &ZoomRecord) -> Result<()> {
        if record.valid_count == 0 {
            return Ok(());
        }
        if !self.window.open
            || self.window.chr_index != record.chr_index
            || record.chr_start as i64 >= self.window.limit
        {
            self.close_window()?;
            self.open_window(
                record.chr_index,
                record.chr_start as i64,
                record.min_value,
                record.max_value,
            );
        }
        self.window.valid_count += record.valid_count as u64;
        self.window.sum_data += record.sum_data as f64;
        self.window.sum_squared += record.sum_squared as f64;
        if record.min_value < self.window.min_value {
            self.window.min_value = record.min_value;
        }
        if record.max_value > self.window.max_value {
            self.window.max_value = record.max_value;
        }
        if record.chr_end as i64 > self.window.end {
            self.window.end = record.chr_end as i64;
        }
        Ok(())
    }

    /// Close the buffered records into one block.
    ///
    /// A block may straddle chromosomes, and should. A zoom record carries its
    /// own chromosome id, unlike a wig section, and a genome of thousands of
    /// scaffolds would otherwise end up with a 32-byte block per scaffold and
    /// an index larger than the data it indexes.
    fn flush_block(&mut self) -> Result<()> {
        if self.records.is_empty() {
            return Ok(());
        }
        let mut body = Vec::with_capacity(self.records.len() * ZOOM_RECORD_SIZE as usize);
        for record in &self.records {
            record.write(&mut body);
        }
        let first = self.records.first().expect("not empty");
        let last = self.records.last().expect("not empty");
        let leaf = LeafItem {
            start_chr: first.chr_index,
            start_base: first.chr_start,
            // Windows never overlap and go out in order, so the last record of
            // the block is also the furthest it reaches.
            end_chr: last.chr_index,
            end_base: last.chr_end,
            offset: 0,
            size: 0,
        };
        self.records.clear();
        self.blocks.push((leaf, body));
        Ok(())
    }

    fn finish(&mut self) -> Result<()> {
        self.close_window()?;
        self.flush_block()
    }
}

impl BbiWriter {
    /// Hand every block the builder has closed to the writer's own pipeline.
    ///
    /// Called as the level is produced rather than once at the end. Two things
    /// come of that: what the pass holds is bounded by `pending_limit` — which
    /// is what bounds the data pass too — instead of by the size of the level,
    /// and the deflate runs on the pool rather than inline on this thread,
    /// which had the workers idle through the longest part of `close()`.
    ///
    /// Appending while the pass reads is safe: it only ever reads offsets
    /// written before it started — the data section, or the level below — and
    /// only ever appends past the end of them.
    fn drain_zoom_blocks(&mut self, builder: &mut ZoomLevelBuilder) -> Result<()> {
        // Taken rather than iterated in place: `submit_block` borrows the
        // writer, and the builder is not part of it.
        for (leaf, body) in std::mem::take(&mut builder.blocks) {
            self.submit_block(leaf, body, None)?;
        }
        Ok(())
    }

    /// Build every zoom level, from the finished data section outwards.
    fn write_zoom_levels(&mut self) -> Result<()> {
        if !self.ladder_frozen {
            self.freeze_zoom_ladder();
        }
        if self.data_items.is_empty() {
            return Ok(());
        }

        // The finest level worth keeping is the first whose records take at
        // most half the room the data does: a summary no smaller than what it
        // summarises is never worth reading in its place.
        //
        // Both sides are measured uncompressed. The two compress nothing like
        // as well as each other — a zoom record is eight numbers, four of them
        // floats, where a wig section is a column of coordinates deflate eats
        // almost entirely — so assuming they do picks a level about one rung
        // too fine, and a rung is a factor of four.
        //
        // The counters only saw items that arrived after the ladder was fixed,
        // so what they hold is scaled up to the whole file.
        let counted = (self.item_count - self.ladder_start_item_count).max(1);
        let sample_scale = self.item_count as f64 / counted as f64;
        let mut first_level = MAX_ZOOM_LEVELS - 1;
        for level in 0..MAX_ZOOM_LEVELS {
            let estimate =
                self.zoom_res_sizes[level] as f64 * sample_scale * ZOOM_RECORD_SIZE as f64;
            if estimate <= self.data_body_size as f64 / 2.0 {
                first_level = level;
                break;
            }
        }

        // Taken, not cloned: nothing reads `data_items` after this, and on a
        // whole-genome 1 bp file the list is a leaf per block — tens of
        // megabytes that were being held twice for no reason.
        let mut source = std::mem::take(&mut self.data_items);
        let mut from_data = true;
        let mut previous_count: Option<u64> = None;
        for level in first_level..MAX_ZOOM_LEVELS {
            let reduction = self.zoom_reductions[level];
            if reduction > 0xFFFF_FFFF {
                break;
            }
            let (count, header, items) = self.write_zoom_level(reduction, &source, from_data)?;
            if count == 0 {
                break;
            }
            self.zoom_headers.push(header);
            source = items;
            from_data = false;
            // A level indexed by a single node answers any query in one index
            // read and one block read, and so does every coarser level. UCSC
            // keeps going down to a handful of records; on a 3 Mb test file
            // each extra level cost 8 kB — an index node is padded to
            // block_size slots — to hold a hundred bytes, and changed no query.
            if count <= self.block_size as u64 {
                break;
            }
            // Data sparse enough that every item is its own window at every
            // resolution, so widening the windows changes nothing.
            if previous_count.is_some_and(|p| count >= p) {
                break;
            }
            previous_count = Some(count);
        }
        Ok(())
    }

    /// Build one zoom level out of `source`, which is the data section for the
    /// first level and the level below for every one after it.
    fn write_zoom_level(
        &mut self,
        reduction: i64,
        source: &[LeafItem],
        from_data: bool,
    ) -> Result<(u64, ZoomHeader, Vec<LeafItem>)> {
        let data_offset = self.sync_cursor()?;
        self.emit(&0u32.to_le_bytes())?; // record count, patched below

        let mut builder = ZoomLevelBuilder::new(reduction, self.items_per_slot);
        // A bigBed summarises the depth of coverage its entries make, so the
        // first level runs the same sweep over them the summary did, this time
        // reporting into the reducer.
        let mut sweep = CoverageSweep::default();

        // Every block placed from here until the level is finished belongs to
        // this level, not to the data section. Set before the first
        // `drain_zoom_blocks`, cleared once the pipeline has drained.
        self.zoom_items = Some(Vec::new());

        // Blocks are consumed in strict source order, which the reducer
        // requires: add_interval, add_record and the sweep all walk a position
        // that only ever moves forward.
        //
        // The result is threaded through `run` so that a failure anywhere in it
        // still clears `zoom_items` on the way out: leaving the writer in zoom
        // mode would send the *next* file's data blocks into a level's index.
        let run = (|| -> Result<()> {
            for item in source {
                let raw = {
                    let source = self.sink()?.as_source()?;
                    source.read_exact_at(item.offset, item.size as usize)?
                };
                let block = if self.compression_level > 0 {
                    super::block::decompress(raw, self.uncompress_buffer_size as u32, &self.path)?
                } else {
                    raw
                };
                if from_data && self.is_bigwig() {
                    let header = read_wig_header(&block, &self.path)?;
                    for i in 0..header.item_count as usize {
                        let item = read_wig_item(&block, &header, i, &self.path)?;
                        builder.add_interval(item.chr_index, item.start, item.end, item.value)?;
                        // Inside the loop, not only after it: one bedGraph
                        // interval can be megabases wide, and at the finest
                        // reduction that is a block of records on its own.
                        self.drain_zoom_blocks(&mut builder)?;
                    }
                } else if from_data {
                    // Collected rather than reduced inside the visitor: the
                    // sweep and the builder both need `&mut`, and the visitor
                    // already borrows the block.
                    let mut records = Vec::new();
                    super::block::visit_bed_records(&block, &self.path, |chr, start, end| {
                        records.push((chr, start, end))
                    })?;
                    for (chr, start, end) in records {
                        for (chr, s, e, depth) in sweep.add(chr, start, end) {
                            builder.add_interval(chr, s, e, depth)?;
                        }
                        self.drain_zoom_blocks(&mut builder)?;
                    }
                } else {
                    let count = block.len() / ZOOM_RECORD_SIZE as usize;
                    for i in 0..count {
                        let record = ZoomRecord::read(&block, i * ZOOM_RECORD_SIZE as usize);
                        builder.add_record(&record)?;
                    }
                    self.drain_zoom_blocks(&mut builder)?;
                }
            }
            if from_data && !self.is_bigwig() {
                for (chr, s, e, depth) in sweep.finish() {
                    builder.add_interval(chr, s, e, depth)?;
                }
            }
            builder.finish()?;
            self.drain_zoom_blocks(&mut builder)?;
            // Every block of this level is placed once the pipeline is empty,
            // which is what makes the items below complete and in file order.
            self.drain()
        })();
        let items = self.zoom_items.take().unwrap_or_default();
        run?;

        if builder.record_count == 0 {
            return Ok((
                0,
                ZoomHeader {
                    reduction_level: 0,
                    data_offset: 0,
                    index_offset: 0,
                },
                Vec::new(),
            ));
        }

        let count = to_bbi_u32(builder.record_count as i64, "zoomCount")?;
        self.patch(data_offset, &count.to_le_bytes())?;

        let index_offset = self.sync_cursor()?;
        let mut bytes = Vec::new();
        super::rtree::write_tree(
            &items,
            index_offset,
            self.block_size,
            self.items_per_slot as u32,
            index_offset,
            &mut |b| bytes.extend_from_slice(b),
        )?;
        self.emit(&bytes)?;

        Ok((
            builder.record_count,
            ZoomHeader {
                reduction_level: to_bbi_u32(reduction, "reductionLevel")?,
                data_offset,
                index_offset,
            },
            items,
        ))
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::bbi::{BbiReader, Zoom};

    fn temp(name: &str) -> std::path::PathBuf {
        let dir = std::env::temp_dir().join("gwseq_writer_tests");
        std::fs::create_dir_all(&dir).unwrap();
        dir.join(name)
    }

    fn sizes(pairs: &[(&str, i64)]) -> ChrMap {
        ChrMap::from_entries(pairs.iter().map(|(a, b)| ((*a).to_string(), *b)))
    }

    fn wig_options(chr_sizes: Option<ChrMap>, parallel: i64) -> BbiWriterOptions {
        BbiWriterOptions {
            kind: BbiKind::BigWig,
            chr_sizes,
            parallel,
            ..Default::default()
        }
    }

    /// Every read path a caller has, over the whole of one chromosome.
    fn read_back(path: &std::path::Path, chr: &str, end: i64) -> Vec<f32> {
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        let request = crate::bbi::ValuesRequest::new(
            crate::genomic::Locs::spans(&[chr.to_string()], &[0], &[end]).unwrap(),
        )
        .bin_size(1.0)
        .def_value(f32::NAN);
        reader
            .read_values(&request)
            .unwrap()
            .into_raw_vec_and_offset()
            .0
    }

    #[test]
    fn a_written_bigwig_reads_back_value_for_value() {
        let path = temp("values.bigwig");
        let values: Vec<f32> = (0..5000).map(|i| (i as f32) * 0.25).collect();
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 5000)])), 1),
            )
            .unwrap();
            w.write_values("chr1", 0, 1, &values).unwrap();
            w.close().unwrap();
        }
        let got = read_back(&path, "chr1", 5000);
        assert_eq!(got.len(), 5000);
        assert_eq!(got, values);
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn the_reader_sees_the_header_the_writer_wrote() {
        let path = temp("header.bigwig");
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 1000), ("chr2", 500)])), 1),
            )
            .unwrap();
            for i in 0..1000 {
                w.write_value("chr1", i, i + 1, i as f32).unwrap();
            }
            for i in 0..500 {
                w.write_value("chr2", i, i + 1, 1.0).unwrap();
            }
            w.close().unwrap();
        }
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        assert_eq!(reader.kind(), BbiKind::BigWig);
        assert_eq!(reader.header().version, BBI_OUTPUT_VERSION);
        assert_eq!(reader.chr_sizes().len(), 2);
        assert_eq!(reader.chr_sizes().resolve("chr2").unwrap().size, 500);
        let summary = reader.total_summary();
        assert_eq!(summary.bases_covered, 1500);
        assert_eq!(summary.min_value, 0.0);
        assert_eq!(summary.max_value, 999.0);
        // sum over 0..1000 is 499500, plus 500 ones.
        assert_eq!(summary.sum_data, 500_000.0);
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn a_written_file_carries_zoom_levels_that_read_back() {
        let path = temp("zoom.bigwig");
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 100_000)])), 1),
            )
            .unwrap();
            let values: Vec<f32> = (0..100_000).map(|i| (i % 97) as f32).collect();
            w.write_values("chr1", 0, 1, &values).unwrap();
            w.close().unwrap();
        }
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        assert!(
            !reader.zoom_headers().is_empty(),
            "no zoom levels were written"
        );
        // Reading through a zoom level has to agree with reading through the
        // data, which is what a summary is for.
        let request = |zoom| {
            crate::bbi::ValuesRequest::new(
                crate::genomic::Locs::spans(&["chr1".into()], &[0], &[100_000]).unwrap(),
            )
            .bin_size(10_000.0)
            .zoom(zoom)
        };
        let full = reader.read_values(&request(Zoom::Full)).unwrap();
        let zoomed = reader.read_values(&request(Zoom::Auto)).unwrap();
        for (a, b) in full.iter().zip(zoomed.iter()) {
            assert!((a - b).abs() < 0.5, "{a} vs {b}");
        }
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn a_written_bigbed_reads_back_entry_for_entry() {
        let path = temp("entries.bigbed");
        let fields: IndexMap<String, String> = [
            ("chr", "string"),
            ("start", "uint"),
            ("end", "uint"),
            ("name", "string"),
            ("score", "uint"),
        ]
        .into_iter()
        .map(|(a, b)| (a.to_string(), b.to_string()))
        .collect();
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                BbiWriterOptions {
                    kind: BbiKind::BigBed,
                    chr_sizes: Some(sizes(&[("chr1", 10_000)])),
                    fields: fields.clone(),
                    parallel: 1,
                    ..Default::default()
                },
            )
            .unwrap();
            for i in 0..1000i64 {
                let values: IndexMap<String, String> = [
                    ("name".to_string(), format!("item{i}")),
                    ("score".to_string(), (i % 1000).to_string()),
                ]
                .into_iter()
                .collect();
                w.write_entry("chr1", i * 5, i * 5 + 8, &values).unwrap();
            }
            w.close().unwrap();
        }
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        assert_eq!(reader.kind(), BbiKind::BigBed);
        assert_eq!(
            reader
                .auto_sql()
                .keys()
                .map(String::as_str)
                .collect::<Vec<_>>(),
            ["chrom", "chromStart", "chromEnd", "name", "score"]
        );
        let request = crate::bbi::EntriesRequest::new(
            crate::genomic::Locs::spans(&["chr1".into()], &[0], &[10_000]).unwrap(),
        );
        let per_locus = reader.read_entries(&request).unwrap();
        assert_eq!(per_locus[0].len(), 1000);
        assert_eq!(per_locus[0][0].start, 0);
        assert_eq!(per_locus[0][0].end, 8);
        assert_eq!(per_locus[0][7].fields[0].1, "item7");
        assert_eq!(per_locus[0][999].start, 4995);
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn the_deflate_pipeline_writes_the_same_bytes_as_the_serial_path() {
        // The pipeline places blocks in submission order, so the file is the
        // one a single thread would have written — byte for byte, which is the
        // whole promise.
        let values: Vec<f32> = (0..40_000).map(|i| ((i * 7) % 251) as f32).collect();
        let mut written: Vec<Vec<u8>> = Vec::new();
        for parallel in [1i64, 4] {
            let path = temp(&format!("parallel{parallel}.bigwig"));
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 40_000)])), parallel),
            )
            .unwrap();
            w.write_values("chr1", 0, 1, &values).unwrap();
            w.close().unwrap();
            written.push(std::fs::read(&path).unwrap());
            std::fs::remove_file(&path).ok();
        }
        assert_eq!(written[0].len(), written[1].len());
        assert!(written[0] == written[1], "the two files differ");
    }

    #[test]
    fn a_value_that_splits_a_section_is_still_written() {
        // The section refuses to widen once it is large enough, and the value
        // that made it refuse has to land in the *next* section rather than be
        // dropped. Worth its own test because the retry was once inside a
        // `debug_assert`, which compiles the call away in a release build: this
        // passed in debug and lost one value per split in the wheel. Run the
        // suite with `--release` as well as without.
        let path = temp("split.bigwig");
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 100_000)])), 1),
            )
            .unwrap();
            // A thousand 10 bp values, then one of a different span: too big a
            // section to widen, so it splits.
            w.write_values("chr1", 0, 10, &[1.0; 1000]).unwrap();
            w.write_value("chr1", 20_000, 20_005, 2.0).unwrap();
            w.write_value("chr1", 20_005, 20_010, 2.0).unwrap();
            w.close().unwrap();
        }
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        // 10 000 bases of the run plus 10 of the two values after it.
        assert_eq!(reader.total_summary().bases_covered, 10_010);
        let request = crate::bbi::ValuesRequest::new(
            crate::genomic::Locs::spans(&["chr1".into()], &[20_000], &[20_010]).unwrap(),
        )
        .bin_size(1.0)
        .def_value(-9.0);
        let got = reader.read_values(&request).unwrap();
        assert!(
            got.iter().all(|v| *v == 2.0),
            "the value that split the section was dropped: {got:?}"
        );
        std::fs::remove_file(&path).ok();
    }

    /// What the converters promise an interrupted run leaves behind: nothing.
    /// A zero-magic stub is no reader's idea of a bigWig, which is right, but
    /// it is still a file the caller has to clean up.
    #[test]
    fn an_abandoned_file_is_removed() {
        let path = temp("abandoned.bigwig");
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 100)])), 1),
            )
            .unwrap();
            w.write_value("chr1", 0, 10, 1.0).unwrap();
            w.abandon();
            // Twice, and after the sink has already gone.
            w.abandon();
        }
        assert!(!path.exists(), "abandon left {} behind", path.display());
    }

    #[test]
    fn out_of_order_and_overlapping_values_are_refused() {
        let path = temp("order.bigwig");
        let mut w = BbiWriter::create(
            path.to_str().unwrap(),
            wig_options(Some(sizes(&[("chr1", 1000), ("chr2", 1000)])), 1),
        )
        .unwrap();
        w.write_value("chr1", 100, 200, 1.0).unwrap();
        let err = w
            .write_value("chr1", 150, 250, 1.0)
            .unwrap_err()
            .to_string();
        assert!(err.contains("starts before the end 200"), "{err}");
        // And a chromosome cannot be returned to once left.
        let mut w2 = BbiWriter::create(
            path.to_str().unwrap(),
            wig_options(Some(sizes(&[("chr1", 1000), ("chr2", 1000)])), 1),
        )
        .unwrap();
        w2.write_value("chr1", 0, 10, 1.0).unwrap();
        w2.write_value("chr2", 0, 10, 1.0).unwrap();
        let err = w2.write_value("chr1", 20, 30, 1.0).unwrap_err().to_string();
        assert!(err.contains("was already written"), "{err}");
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn a_value_hanging_over_a_chromosome_is_clipped_and_one_past_it_is_refused() {
        let path = temp("clip.bigwig");
        let mut w = BbiWriter::create(
            path.to_str().unwrap(),
            wig_options(Some(sizes(&[("chr1", 95)])), 1),
        )
        .unwrap();
        // A 10 bp grid over a 95 bp chromosome: the last bin overshoots by 5.
        w.write_values("chr1", 0, 10, &[1.0; 10]).unwrap();
        assert_eq!(w.clipped_count(), 1);
        let err = w.write_value("chr1", 95, 105, 1.0).unwrap_err().to_string();
        assert!(err.contains("starts past the end"), "{err}");
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn a_non_finite_value_is_skipped_rather_than_written() {
        let path = temp("skip.bigwig");
        {
            let mut w = BbiWriter::create(
                path.to_str().unwrap(),
                wig_options(Some(sizes(&[("chr1", 5)])), 1),
            )
            .unwrap();
            w.write_values("chr1", 0, 1, &[1.0, f32::NAN, 3.0, f32::INFINITY, 5.0])
                .unwrap();
            assert_eq!(w.skipped_count(), 2);
            w.close().unwrap();
        }
        let got = read_back(&path, "chr1", 5);
        assert_eq!(got[0], 1.0);
        assert!(got[1].is_nan(), "the gap is a gap");
        assert_eq!(got[2], 3.0);
        assert!(got[3].is_nan());
        assert_eq!(got[4], 5.0);
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn an_undeclared_chromosome_grows_to_what_is_written_to_it() {
        let path = temp("infer.bigwig");
        {
            let mut w = BbiWriter::create(path.to_str().unwrap(), wig_options(None, 1)).unwrap();
            w.write_value("chrZ", 100, 250, 1.0).unwrap();
            w.close().unwrap();
        }
        let reader = BbiReader::open(path.to_str().unwrap(), 1, 1.0 / 3.0, None, None).unwrap();
        assert_eq!(reader.chr_sizes().resolve("chrZ").unwrap().size, 250);
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn writing_through_a_closed_writer_is_refused() {
        let path = temp("closed.bigwig");
        let mut w = BbiWriter::create(
            path.to_str().unwrap(),
            wig_options(Some(sizes(&[("chr1", 10)])), 1),
        )
        .unwrap();
        w.write_value("chr1", 0, 5, 1.0).unwrap();
        w.close().unwrap();
        w.close().unwrap(); // idempotent
        let err = w.write_value("chr1", 5, 10, 1.0).unwrap_err().to_string();
        assert!(err.contains("file is closed"), "{err}");
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn the_open_time_checks_refuse_what_the_format_cannot_hold() {
        let path = temp("bad.bigwig");
        let p = path.to_str().unwrap();
        let bad = |o: BbiWriterOptions| BbiWriter::create(p, o).unwrap_err().to_string();
        assert!(bad(BbiWriterOptions {
            items_per_slot: Some(0),
            ..Default::default()
        })
        .contains("items_per_slot 0 invalid"));
        assert!(bad(BbiWriterOptions {
            items_per_slot: Some(70_000),
            ..Default::default()
        })
        .contains("items_per_slot 70000 invalid"));
        assert!(bad(BbiWriterOptions {
            block_size: 1,
            ..Default::default()
        })
        .contains("block_size 1 invalid"));
        assert!(bad(BbiWriterOptions {
            compression_level: 10,
            ..Default::default()
        })
        .contains("compression_level 10 invalid"));
        assert!(bad(BbiWriterOptions {
            chr_sizes: Some(sizes(&[("chr1", 0)])),
            ..Default::default()
        })
        .contains("must be positive"));
        assert!(
            BbiWriter::create("https://example.org/x.bigwig", BbiWriterOptions::default())
                .unwrap_err()
                .to_string()
                .contains("cannot be written to")
        );
        std::fs::remove_file(&path).ok();
    }

    #[test]
    fn the_coverage_sweep_reports_the_depth_of_nested_entries() {
        let mut sweep = CoverageSweep::default();
        let mut runs = Vec::new();
        // [0,10) and [5,20) overlap on [5,10), where the depth is two.
        runs.extend(sweep.add(0, 0, 10));
        runs.extend(sweep.add(0, 5, 20));
        runs.extend(sweep.finish());
        assert_eq!(runs, [(0, 0, 5, 1.0), (0, 5, 10, 2.0), (0, 10, 20, 1.0)]);
    }

    #[test]
    fn the_coverage_sweep_closes_a_chromosome_before_starting_the_next() {
        let mut sweep = CoverageSweep::default();
        let mut runs = Vec::new();
        runs.extend(sweep.add(0, 0, 10));
        runs.extend(sweep.add(1, 0, 10));
        runs.extend(sweep.finish());
        assert_eq!(runs, [(0, 0, 10, 1.0), (1, 0, 10, 1.0)]);
    }
}