gwseq_io/bbi/

reader.rs

//! The bigWig / bigBed reader and its request builders.
//!
//! The plumbing only; the extraction kernels are in [`super::extract`].

use std::fmt::Write as _;
use std::io::Write as _;
use std::sync::Arc;

use ndarray::{Array1, Array2};

use crate::bbi::extract::Extraction;
use crate::bbi::header::{BbiHeader, BbiKind, TotalSummary, ZoomHeader};
use crate::error::{Error, Result};
use crate::genomic::{BinMode, ChrMap, IndexedLocs, LocBatch, Locs, Reduce};
use crate::parallel::Executor;
use crate::progress::{ProgressFn, ProgressTracker};
use crate::source::ByteSource;

/// Which zoom level a read uses.
#[derive(Debug, Clone, Copy, Default)]
pub enum Zoom {
    /// Full data. The default, and the only option for a bigBed — which carries
    /// no zoom data, so naming a level there is refused rather than ignored.
    #[default]
    Full,
    /// The coarsest level whose bin size is under `bin_size * zoom_correction`,
    /// which may be the full data.
    Auto,
    Level(usize),
}

/// Everything a reader holds that `close()` gives back.
///
/// Split out so that `close()` is `self.inner.take()`: dropping it drops the
/// executor, which joins its threads, and the source, which releases the
/// handle. The headers live outside it and stay readable after close, as the
/// API documents.
#[derive(Debug)]
struct Inner {
    source: Arc<dyn ByteSource>,
    executor: Executor,
}

#[derive(Debug)]
pub struct BbiReader {
    inner: Option<Inner>,
    path: String,
    zoom_correction: f64,

    // Read once at open, and readable after close.
    pub(crate) header: BbiHeader,
    pub(crate) zoom_headers: Vec<ZoomHeader>,
    pub(crate) total_summary: TotalSummary,
    pub(crate) chr_map: ChrMap,
    /// Chromosome names indexed by the file's own chromosome index.
    pub(crate) chr_names: Vec<String>,
    /// bigBed only: the entry columns the file declares.
    pub(crate) auto_sql: indexmap::IndexMap<String, String>,
}

/// The three resolutions a read pulls apart, which are not always the request's
/// own `bin_size`.
///
/// One struct rather than three positional arguments because two of them are an
/// `Option<f64>` and swapping those would still compile.
#[derive(Default, Clone, Copy)]
struct Grid {
    /// Output bins per locus. `None` lets the locus span and the bin size
    /// decide between them.
    bin_count: Option<usize>,
    /// What the loci snap to, in base pairs. `None` is the request's own
    /// `bin_size`. The entry paths pass 1.0: an entry is not a bin, so snapping
    /// a whole-chromosome locus down to a bin boundary would only lose the
    /// entries past the last whole one.
    snap: Option<f64>,
    /// What a zoom level has to be finer than, when the output bins do not say.
    /// `quantify` pins `bin_count` to 1 and still reads at `bin_size`, so its
    /// single output bin says nothing about the resolution asked for.
    zoom: Option<f64>,
}

impl Grid {
    /// What `read_entries`, `read_all_entries` and `to_bed` read on: one bin
    /// per locus, and the loci on the base grid rather than the request's.
    fn entries() -> Self {
        Self {
            bin_count: Some(1),
            snap: Some(1.0),
            zoom: None,
        }
    }
}

impl BbiReader {
    pub fn open(
        path: &str,
        parallel: i64,
        zoom_correction: f64,
        block_size: Option<u64>,
        max_blocks: Option<usize>,
    ) -> Result<Self> {
        let source = crate::source::open(path, block_size, max_blocks)?;
        Self::from_source(source, path, parallel, zoom_correction)
    }

    /// Open over an already-built source. What `gwseq_io::open` calls once it
    /// has sniffed the magic, so the file is not opened twice.
    pub(crate) fn from_source(
        source: Arc<dyn ByteSource>,
        path: &str,
        parallel: i64,
        zoom_correction: f64,
    ) -> Result<Self> {
        let header = super::header::read_header(source.as_ref())?;
        let zoom_headers = super::header::read_zoom_headers(source.as_ref(), header.zoom_levels)?;
        let total_summary =
            super::header::read_total_summary(source.as_ref(), header.total_summary_offset)?;
        let (chr_map, _tree) = super::chr_tree::read(source.as_ref(), header.chr_tree_offset)?;
        // Only a bigBed has columns; a bigWig's field_count is 0 and its
        // autoSql offset points at nothing.
        let auto_sql = if header.kind.is_bigbed() {
            super::header::read_auto_sql(
                source.as_ref(),
                header.auto_sql_offset,
                header.field_count,
            )?
        } else {
            indexmap::IndexMap::new()
        };
        // Chromosome names by the index the file gives them, which is what data
        // records and R-tree items carry and what an entry has to be named from.
        let mut chr_names =
            vec![String::new(); chr_map.iter().map(|e| e.index + 1).max().unwrap_or(0)];
        for entry in chr_map.iter() {
            chr_names[entry.index] = entry.id.clone();
        }
        let executor = Executor::new(parallel)?;
        Ok(Self {
            inner: Some(Inner { source, executor }),
            path: path.to_string(),
            zoom_correction,
            header,
            zoom_headers,
            total_summary,
            chr_map,
            chr_names,
            auto_sql,
        })
    }

    pub fn kind(&self) -> BbiKind {
        self.header.kind
    }
    pub fn path(&self) -> &str {
        &self.path
    }
    pub fn chr_sizes(&self) -> &ChrMap {
        &self.chr_map
    }
    pub fn header(&self) -> &BbiHeader {
        &self.header
    }
    pub fn zoom_headers(&self) -> &[ZoomHeader] {
        &self.zoom_headers
    }
    pub fn total_summary(&self) -> &TotalSummary {
        &self.total_summary
    }
    /// bigBed only: the entry columns the file declares, in file order, read at
    /// open from the autoSql block. Empty for a bigWig, and empty for a bigBed
    /// whose header names no autoSql offset.
    pub fn auto_sql(&self) -> &indexmap::IndexMap<String, String> {
        &self.auto_sql
    }
    pub fn is_closed(&self) -> bool {
        self.inner.is_none()
    }
    pub fn parallel(&self) -> usize {
        self.inner.as_ref().map_or(0, |i| i.executor.parallel())
    }

    /// Give back the threads and the file handle. Idempotent.
    pub fn close(&mut self) {
        if let Some(inner) = self.inner.take() {
            inner.source.close();
        }
    }

    /// At the top of every call that reads, not deep inside one, so the error
    /// names the file rather than a handle.
    fn inner(&self) -> Result<&Inner> {
        self.inner.as_ref().ok_or_else(|| Error::Closed {
            path: self.path.clone(),
        })
    }

    /// Everything the read paths do before they diverge: resolve the loci,
    /// split them into batches, pick a zoom level, and find the root of the
    /// index that level's data hangs off.
    #[allow(clippy::type_complexity)]
    fn prepare<'a>(
        &'a self,
        inner: &'a Inner,
        locs: &Locs,
        common: &ReadCommon,
        grid: Grid,
    ) -> Result<(
        IndexedLocs,
        Vec<LocBatch>,
        Option<usize>,
        u64,
        ProgressTracker,
    )> {
        // bigBed files carry no zoom data: the entries are piled up from the
        // full index whatever is asked for, so an explicit request is refused
        // rather than silently ignored.
        if self.header.kind.is_bigbed() && !matches!(common.zoom, Zoom::Full) {
            return Err(Error::invalid("zoom is only supported for bigwig files"));
        }
        let indexed = IndexedLocs::build(
            &self.chr_map,
            locs,
            grid.snap.unwrap_or(common.bin_size),
            grid.bin_count,
            common.full_bin,
        )?;
        let (batches, coverage) = indexed.batches(inner.executor.parallel());
        let level = self.select_zoom(
            grid.zoom
                .unwrap_or_else(|| indexed.effective_bin_size(common.bin_size)),
            common.zoom,
        )?;
        let index_offset = match level {
            Some(i) => self.zoom_headers[i].index_offset,
            None => self.header.full_index_offset,
        };
        super::header::check_data_tree_magic(inner.source.as_ref(), index_offset)?;
        let tree_root = index_offset + super::header::DATA_TREE_HEADER_SIZE;
        Ok((
            indexed,
            batches,
            level,
            tree_root,
            ProgressTracker::with_callback(coverage, common.progress.clone()),
        ))
    }

    fn extraction<'a>(
        &'a self,
        inner: &'a Inner,
        indexed: &'a IndexedLocs,
        batches: &'a [LocBatch],
        level: Option<usize>,
        tree_root: u64,
        tracker: &'a ProgressTracker,
    ) -> Extraction<'a> {
        Extraction {
            source: inner.source.as_ref(),
            locs: indexed,
            batches,
            tree_root,
            zoom: level.is_some(),
            uncompress_buffer_size: self.header.uncompress_buffer_size,
            tracker,
        }
    }

    /// (loci, bins) of `f32`.
    pub fn read_values(&self, req: &ValuesRequest) -> Result<Array2<f32>> {
        let inner = self.inner()?;
        let (indexed, batches, level, root, tracker) = self.prepare(
            inner,
            &req.common.locs,
            &req.common,
            Grid {
                bin_count: req.common.bin_count,
                ..Grid::default()
            },
        )?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        // For a bigBed the "values" are the depth of coverage its entries make,
        // which has no bin_mode to choose between — the three coincide.
        let flat = if self.header.kind.is_bigbed() {
            super::extract::entries_pileup(
                &ex,
                &inner.executor,
                &self.auto_sql,
                req.common.def_value,
            )?
        } else {
            super::extract::values(&ex, &inner.executor, req.bin_mode, req.common.def_value)?
        };
        tracker.done_report();
        let rows = indexed.locs.len();
        let cols = indexed.bin_count;
        Array2::from_shape_vec((rows, cols), flat)
            .map_err(|e| Error::invalid(format!("output shape {rows}x{cols}: {e}")))
    }

    /// One value per locus.
    pub fn quantify(&self, req: &QuantifyRequest) -> Result<Array1<f32>> {
        let inner = self.inner()?;
        let is_bigbed = self.header.kind.is_bigbed();
        // A bigWig reads a single bin per locus, the extraction weighting each
        // value by the bases it covers, so the reduction runs over the bases
        // without the bins having to exist. `zoom_bin_size` carries the
        // resolution the caller asked for past that pinned bin_count.
        //
        // A bigBed has no values to weight — its signal is the depth its
        // entries make, which exists only once they are piled up — so it is
        // binned at the size the caller asked for and reduced over those bins,
        // as read_values bins it.
        let (indexed, batches, level, root, tracker) = self.prepare(
            inner,
            &req.common.locs,
            &req.common,
            Grid {
                bin_count: if is_bigbed {
                    req.common.bin_count
                } else {
                    Some(1)
                },
                zoom: Some(req.common.bin_size),
                ..Grid::default()
            },
        )?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        let mut stats = if is_bigbed {
            let pileup = super::extract::entries_pileup(
                &ex,
                &inner.executor,
                &self.auto_sql,
                req.common.def_value,
            )?;
            super::extract::pileup_stats(&indexed, &pileup)
        } else {
            super::extract::values_stats(&ex, &inner.executor)?
        };
        tracker.done_report();

        // The uncovered part of a locus still takes part, with def_value
        // standing in for it, unless def_value is NaN — which is how a caller
        // asks for the covered part alone. Filling it in rather than only
        // stretching the denominator keeps every reduce mode agreeing on what
        // those bases hold: min and max see def_value, and a non-zero one
        // reaches the sum instead of reading as zero.
        let def = req.common.def_value;
        if !def.is_nan() {
            for loc in &indexed.locs {
                let s = &mut stats[loc.row(indexed.bin_count)];
                // Counted in whatever unit the extraction accumulates: bases of
                // the locus for the values, bins of it for the bed pileup —
                // where nothing is ever missing, the pileup having already put
                // def_value in the bins no entry reached.
                let total = if is_bigbed {
                    (loc.output_end - loc.output_start) as i64
                } else {
                    loc.binned_end - loc.binned_start
                };
                let missing = total - s.count;
                if missing <= 0 {
                    continue;
                }
                s.add_repeated(def, missing);
            }
        }
        Ok(Array1::from_vec(
            stats.iter().map(|s| s.reduce(req.reduce, def)).collect(),
        ))
    }

    /// One value per bin, reduced across loci.
    pub fn profile(&self, req: &ProfileRequest) -> Result<Array1<f32>> {
        let inner = self.inner()?;
        let (indexed, batches, level, root, tracker) = self.prepare(
            inner,
            &req.common.locs,
            &req.common,
            Grid {
                bin_count: req.common.bin_count,
                ..Grid::default()
            },
        )?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        // `def_value` reaches the pileup, so a bin no entry covered holds what
        // the caller asked for rather than a bare 0. Dropping it would make
        // `profile(def_value=NaN)` report 0 over an uncovered region where
        // read_values reports NaN.
        let mut stats = if self.header.kind.is_bigbed() {
            let pileup = super::extract::entries_pileup(
                &ex,
                &inner.executor,
                &self.auto_sql,
                req.common.def_value,
            )?;
            super::extract::pileup_profile(&indexed, &pileup)
        } else {
            super::extract::values_profile(&ex, &inner.executor, req.bin_mode)?
        };
        tracker.done_report();

        // A locus whose bin held no data still takes part in the profile, with
        // def_value standing in for it, unless def_value is NaN — which is how
        // a caller asks for those loci to be left out of the reduction instead.
        let def = req.common.def_value;
        if !def.is_nan() {
            let loc_count = indexed.locs.len() as i64;
            for s in &mut stats {
                let missing = loc_count - s.count;
                if missing <= 0 {
                    continue;
                }
                s.add_repeated(def, missing);
            }
        }
        Ok(Array1::from_vec(
            stats.iter().map(|s| s.reduce(req.reduce, def)).collect(),
        ))
    }

    /// bigBed entries, per locus, in the order the loci were given.
    ///
    /// Each locus gets its own full list, so two overlapping loci both report
    /// the entries they share.
    pub fn read_entries(&self, req: &EntriesRequest) -> Result<Vec<Vec<super::BedEntry>>> {
        let inner = self.inner()?;
        self.require_bigbed("read_entries")?;
        self.check_col_count(req.col_count, 3)?;
        let (indexed, batches, level, root, tracker) =
            self.prepare(inner, &req.common.locs, &req.common, Grid::entries())?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        let out = super::extract::entries(
            &ex,
            &inner.executor,
            &self.auto_sql,
            &self.chr_names,
            req.col_count,
        )?;
        tracker.done_report();
        Ok(out)
    }

    /// Every entry on the named chromosomes, in chromosome then coordinate
    /// order.
    pub fn read_all_entries(&self, req: &EntriesRequest) -> Result<Vec<super::BedEntry>> {
        let inner = self.inner()?;
        self.require_bigbed("read_all_entries")?;
        self.check_col_count(req.col_count, 3)?;
        let locs = Locs::whole_chromosomes(&self.chr_map, &req.common.locs.chr_ids)?;
        let (indexed, batches, level, root, tracker) =
            self.prepare(inner, &locs, &req.common, Grid::entries())?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        let by_chr = super::extract::entries(
            &ex,
            &inner.executor,
            &self.auto_sql,
            &self.chr_names,
            req.col_count,
        )?;
        tracker.done_report();
        Ok(by_chr.into_iter().flatten().collect())
    }

    /// A method that only exists for one of the two formats says which.
    fn require_bigbed(&self, what: &str) -> Result<()> {
        if self.header.kind.is_bigbed() {
            Ok(())
        } else {
            Err(Error::invalid(format!("{what} only for bigbed")))
        }
    }

    fn require_bigwig(&self, what: &str) -> Result<()> {
        if self.header.kind.is_bigbed() {
            Err(Error::invalid(format!("{what} only for bigwig")))
        } else {
            Ok(())
        }
    }

    /// Check a `col_count` against the columns the file holds.
    ///
    /// `min` is 3 for the entry readers, whose entries are built from the
    /// coordinates, and 1 for `to_bed`, which writes as many columns as it is
    /// told and has no such floor.
    pub(crate) fn check_col_count(&self, col_count: usize, min: usize) -> Result<()> {
        if col_count == 0 {
            return Ok(());
        }
        if col_count < min {
            return Err(Error::invalid(format!(
                "col_count {col_count} must be 0 or at least {min}"
            )));
        }
        if col_count > self.header.field_count as usize {
            return Err(Error::invalid(format!(
                "col_count {col_count} exceeds number of fields {}",
                self.header.field_count
            )));
        }
        Ok(())
    }

    /// Bytes of data the file holds at a given resolution.
    ///
    /// What a whole-file walk sizes its pieces from: a zoom level holds a
    /// fraction of what the full data does, so a walk reading one is worth
    /// splitting far less finely.
    pub fn data_size(&self, zoom: Option<usize>) -> u64 {
        match zoom.and_then(|i| self.zoom_headers.get(i)) {
            Some(z) => z.index_offset.saturating_sub(z.data_offset),
            None => self
                .header
                .full_index_offset
                .saturating_sub(self.header.full_data_offset),
        }
    }

    /// Base pairs the file's chromosomes hold between them.
    pub fn genome_size(&self) -> i64 {
        self.chr_map.genome_size()
    }

    /// Successive windows of values over whole chromosomes.
    ///
    /// A plain [`Iterator`], exhausted after one pass. `len()` and `locs()`
    /// are on the iterator, so
    /// `iter.locs().iter().zip(iter)` works once, as documented.
    pub fn iter_all_values(
        &self,
        req: &ValuesRequest,
        window: i64,
    ) -> Result<super::ValuesWindows<'_>> {
        super::extract::ValuesWindows::plan(self, req, window)
    }

    pub fn iter_all_entries(
        &self,
        req: &EntriesRequest,
        window: i64,
    ) -> Result<super::EntryWindows<'_>> {
        super::extract::EntryWindows::plan(self, req, window)
    }

    /// Every bin of the named chromosomes, as `chr\tstart\tend\tvalue`.
    ///
    /// The values are the ones [`Self::read_values`] gives for the same
    /// chromosomes at the same `bin_size`, `bin_mode`, `full_bin`, `def_value`
    /// and `zoom` — the export bins, it does not copy the file's own intervals.
    /// A default request therefore writes one interval per base, and a
    /// `bin_size` of 10,000 writes one per 10,000 bases.
    ///
    /// With `merge_bins`, adjacent bins carrying the same value become one
    /// interval — the shape a bedGraph is usually in, and what keeps a
    /// whole-genome export from being one line per base wherever the file says
    /// nothing. Without it every bin is its own interval, so the line count is
    /// the bin count and the grid is visible in the file.
    ///
    /// A `def_value` of NaN is how a caller asks for the covered part alone: a
    /// bin holding NaN is not written at all, so the gaps stay gaps.
    ///
    /// `bin_count` is ignored — the bins are the chromosome's, not a fixed
    /// number of them.
    pub fn to_bedgraph(
        &self,
        out: &std::path::Path,
        req: &ValuesRequest,
        merge_bins: bool,
    ) -> Result<()> {
        self.require_bigwig("to_bedgraph")?;
        self.export_bins(out, req, BedGraphSink::new(merge_bins))
    }

    /// The same bins as fixedStep WIG sections.
    ///
    /// A fixedStep section walks a fixed step from its own start, so anything
    /// that breaks the run of bins has to open a new one: a change of
    /// chromosome, a bin left out because its value is NaN, and the shorter
    /// last bin a `full_bin` export ends a chromosome on.
    ///
    /// Bins are written one per line whatever their values, since a fixedStep
    /// section has no way to say "and again" — the `merge_bins` of
    /// [`Self::to_bedgraph`] has no equivalent here.
    pub fn to_wig(&self, out: &std::path::Path, req: &ValuesRequest) -> Result<()> {
        self.require_bigwig("to_wig")?;
        self.export_bins(out, req, WigSink::default())
    }

    /// Every entry as a tab-separated BED line.
    ///
    /// `col_count` here has a floor of 1, not 3: this writes as many columns as
    /// it is told, where the entry readers build an entry from its coordinates
    /// and cannot go below them.
    pub fn to_bed(&self, out: &std::path::Path, req: &EntriesRequest) -> Result<()> {
        let inner = self.inner()?;
        self.require_bigbed("to_bed")?;
        self.check_col_count(req.col_count, 1)?;
        let col_count = if req.col_count == 0 {
            self.header.field_count as usize
        } else {
            req.col_count
        };

        let locs = Locs::whole_chromosomes(&self.chr_map, &req.common.locs.chr_ids)?;
        let (indexed, batches, level, root, tracker) =
            self.prepare(inner, &locs, &req.common, Grid::entries())?;
        let ex = self.extraction(inner, &indexed, &batches, level, root, &tracker);
        let wanted = self.walked_chrs(&indexed);

        let mut writer = std::io::BufWriter::new(
            std::fs::File::create(out).map_err(|e| Error::io(out.to_string_lossy(), e))?,
        );
        let mut line = String::new();
        for batch in &batches {
            ex.walk_bed_batch(*batch, &self.auto_sql, col_count, |entry, _| {
                if !wanted.contains(&entry.chr_index) {
                    return Ok(());
                }
                line.clear();
                line.push_str(self.chr_name(entry.chr_index));
                if col_count >= 2 {
                    let _ = write!(line, "\t{}", entry.start);
                }
                if col_count >= 3 {
                    let _ = write!(line, "\t{}", entry.end);
                }
                for (_, value) in entry.fields.iter().take(col_count.saturating_sub(3)) {
                    line.push('\t');
                    line.push_str(value);
                }
                line.push('\n');
                write_line(&mut writer, &line, out)
            })?;
        }
        flush(&mut writer, out)?;
        tracker.done_report();
        Ok(())
    }

    /// The shared body of `to_bedgraph` and `to_wig`: bin the named
    /// chromosomes exactly as `read_values` bins them and let `sink` turn each
    /// bin into text.
    ///
    /// The read runs on the pool a window at a time, through the very walk
    /// `iter_all_values` hands to a caller, which is what makes an export and a
    /// read of the same chromosomes agree by construction rather than by two
    /// implementations of the same arithmetic. The formatting and the writing
    /// are single-threaded on purpose: the output is one file in chromosome
    /// then coordinate order, so formatted text produced out of order would
    /// have to be held until its turn came, and a window of it is the memory
    /// this avoids.
    fn export_bins(
        &self,
        out: &std::path::Path,
        req: &ValuesRequest,
        mut sink: impl BinSink,
    ) -> Result<()> {
        // Named before the file is created: an export through a closed reader
        // should not leave an empty file behind.
        self.inner()?;
        let walk = super::extract::ValuesWindows::plan(self, req, Self::export_window(req)?)?;
        let bin_size = walk.bin_size();
        // Cloned up front: the walk is consumed by the loop below, and each
        // window needs the chromosome and start coordinate of its own region.
        let locs = walk.locs().to_vec();

        let mut writer = std::io::BufWriter::new(
            std::fs::File::create(out).map_err(|e| Error::io(out.to_string_lossy(), e))?,
        );
        let mut line = String::new();
        for (index, values) in walk.enumerate() {
            let values = values?;
            let (chr, window_start, window_end) = &locs[index];
            for (i, &value) in values.iter().enumerate() {
                // A NaN bin is left out of the file rather than written as the
                // text "NaN", which no reader of either format accepts. That is
                // what makes `def_value = NaN` an export of the covered part
                // alone, and the coordinate discontinuity it leaves is what
                // both sinks break a run on.
                if value.is_nan() {
                    continue;
                }
                let start = window_start + i as i64 * bin_size;
                // The last bin of a `full_bin` chromosome is the short one.
                let end = (start + bin_size).min(*window_end);
                line.clear();
                sink.bin(&mut line, chr, start, end, value);
                if !line.is_empty() {
                    write_line(&mut writer, &line, out)?;
                }
            }
        }
        line.clear();
        sink.finish(&mut line);
        if !line.is_empty() {
            write_line(&mut writer, &line, out)?;
        }
        flush(&mut writer, out)
    }

    /// How much of a chromosome an export reads at a time, in base pairs.
    ///
    /// Two caps, and the narrower one wins. [`EXPORT_WINDOW_BINS`] bounds what
    /// a window costs to hold — 4 MiB of `f32` for a 1 bp export of a human
    /// chromosome as much as for a 10 kb one — and [`EXPORT_WINDOW_BASES`]
    /// bounds how much of a chromosome goes by between two progress reports,
    /// which a bin-only cap would leave at "the whole thing" for any bin wide
    /// enough. Never below one bin, since a bin may be wider than either cap.
    ///
    /// `saturating_mul` because a bin size is any positive whole number of
    /// bases, including one larger than the genome.
    fn export_window(req: &ValuesRequest) -> Result<i64> {
        let bin_size =
            crate::genomic::BinPlan::new(req.common.bin_size, None, req.common.full_bin)?
                .whole_bin_size();
        Ok(bin_size
            .saturating_mul(EXPORT_WINDOW_BINS)
            .min(EXPORT_WINDOW_BASES)
            .max(bin_size))
    }

    /// The chromosomes a walk was asked for.
    ///
    /// A tree item may straddle a chromosome boundary, so a whole-file walk can
    /// reach data from chromosomes that were never asked for. The exporters
    /// filter their output on this.
    fn walked_chrs(&self, locs: &IndexedLocs) -> std::collections::HashSet<u32> {
        locs.locs.iter().map(|l| l.chr_index as u32).collect()
    }

    pub(crate) fn chr_name(&self, index: u32) -> &str {
        self.chr_names
            .get(index as usize)
            .map(String::as_str)
            .unwrap_or("")
    }

    /// Which zoom level to read, or `None` for the full data.
    ///
    /// The automatic choice is the **coarsest level still finer than** the
    /// output bins, so that each bin averages several summaries rather than
    /// inheriting the edges of one. `zoom_correction` is how much finer. It may
    /// well settle on the full data, which is what a file whose coarsest level
    /// is still wider than the output bins has to be read at.
    pub(crate) fn select_zoom(&self, bin_size: f64, zoom: Zoom) -> Result<Option<usize>> {
        let count = self.zoom_headers.len();
        match zoom {
            Zoom::Full => Ok(None),
            Zoom::Level(level) => {
                if level < count {
                    Ok(Some(level))
                } else if count == 0 {
                    Err(Error::invalid("file has no zoom level"))
                } else {
                    Err(Error::invalid(format!(
                        "requested zoom level {level} exceeds max zoom level {}",
                        count - 1
                    )))
                }
            }
            Zoom::Auto => {
                let threshold = (bin_size * self.zoom_correction).round() as i64;
                let mut best: Option<usize> = None;
                let mut best_reduction = 0i64;
                for (i, zoom) in self.zoom_headers.iter().enumerate() {
                    let reduction = zoom.reduction_level as i64;
                    if reduction <= threshold && reduction > best_reduction {
                        best_reduction = reduction;
                        best = Some(i);
                    }
                }
                Ok(best)
            }
        }
    }
}

// ---------------------------------------------------------------------------
// Requests
// ---------------------------------------------------------------------------
//
// A read takes up to fourteen optional parameters, and Rust has no defaulted
// parameters, so each request shape gets a builder — shared by the
// Python layer, the CLI and Rust callers, so there is one place a default is
// written down.

/// Shared by every read: which loci, at what resolution, from which zoom.
pub struct ReadCommon {
    pub locs: Locs,
    pub bin_size: f64,
    pub bin_count: Option<usize>,
    pub full_bin: bool,
    pub def_value: f32,
    pub zoom: Zoom,
    pub progress: Option<ProgressFn>,
}

impl ReadCommon {
    pub fn new(locs: Locs) -> Self {
        Self {
            locs,
            bin_size: 1.0,
            bin_count: None,
            full_bin: false,
            def_value: 0.0,
            zoom: Zoom::Full,
            progress: None,
        }
    }
}

macro_rules! read_common_builders {
    ($t:ty) => {
        impl $t {
            pub fn bin_size(mut self, v: f64) -> Self {
                self.common.bin_size = v;
                self
            }
            pub fn bin_count(mut self, v: usize) -> Self {
                self.common.bin_count = Some(v);
                self
            }
            pub fn full_bin(mut self, v: bool) -> Self {
                self.common.full_bin = v;
                self
            }
            pub fn def_value(mut self, v: f32) -> Self {
                self.common.def_value = v;
                self
            }
            pub fn zoom(mut self, v: Zoom) -> Self {
                self.common.zoom = v;
                self
            }
            pub fn progress(mut self, f: ProgressFn) -> Self {
                self.common.progress = Some(f);
                self
            }
        }
    };
}

pub struct ValuesRequest {
    pub common: ReadCommon,
    pub bin_mode: BinMode,
}

pub struct QuantifyRequest {
    pub common: ReadCommon,
    pub reduce: Reduce,
}

pub struct ProfileRequest {
    pub common: ReadCommon,
    pub bin_mode: BinMode,
    pub reduce: Reduce,
}

pub struct EntriesRequest {
    pub common: ReadCommon,
    /// 0 for all, otherwise at least 3 — or at least 1 for [`BbiReader::to_bed`],
    /// which writes the columns it is told to rather than building an entry out
    /// of them. Columns left out are never parsed, so a narrower read is a
    /// cheaper one.
    pub col_count: usize,
}

read_common_builders!(ValuesRequest);
read_common_builders!(QuantifyRequest);
read_common_builders!(ProfileRequest);
read_common_builders!(EntriesRequest);

impl ValuesRequest {
    pub fn new(locs: Locs) -> Self {
        Self {
            common: ReadCommon::new(locs),
            bin_mode: BinMode::Mean,
        }
    }
    pub fn bin_mode(mut self, v: BinMode) -> Self {
        self.bin_mode = v;
        self
    }
}

impl QuantifyRequest {
    pub fn new(locs: Locs) -> Self {
        Self {
            common: ReadCommon::new(locs),
            reduce: Reduce::Mean,
        }
    }
    pub fn reduce(mut self, v: Reduce) -> Self {
        self.reduce = v;
        self
    }
}

impl ProfileRequest {
    pub fn new(locs: Locs) -> Self {
        Self {
            common: ReadCommon::new(locs),
            bin_mode: BinMode::Mean,
            reduce: Reduce::Mean,
        }
    }
    pub fn bin_mode(mut self, v: BinMode) -> Self {
        self.bin_mode = v;
        self
    }
    pub fn reduce(mut self, v: Reduce) -> Self {
        self.reduce = v;
        self
    }
}

impl EntriesRequest {
    pub fn new(locs: Locs) -> Self {
        Self {
            common: ReadCommon::new(locs),
            col_count: 0,
        }
    }
    pub fn col_count(mut self, v: usize) -> Self {
        self.col_count = v;
        self
    }
}

// ---------------------------------------------------------------------------
// Export sinks
// ---------------------------------------------------------------------------

/// Bins per window of the walk the exporters stream through. A window of
/// `f32`s is 4 MiB, which is what a whole-genome export holds however wide its
/// bins are.
const EXPORT_WINDOW_BINS: i64 = 1 << 20;

/// Base pairs per window, the other half of the cap. Progress is reported once
/// per window, so this is the resolution of a progress bar over an export —
/// a few hundred steps across a human chromosome — and it is what keeps a
/// wide-binned export from reading a whole chromosome before saying anything.
const EXPORT_WINDOW_BASES: i64 = 16 << 20;

/// How an export turns a bin into text.
///
/// Both implementations carry a run across bins — an interval being extended,
/// a section being filled — so a bin often appends nothing and the pending
/// line comes out one bin later. That is why this is a trait with a `finish`
/// rather than a closure: the last run has to be written after the last bin.
trait BinSink {
    /// Append whatever lines this bin completes. Bins arrive in chromosome
    /// then coordinate order, and a bin whose value is NaN never arrives.
    fn bin(&mut self, line: &mut String, chr: &str, start: i64, end: i64, value: f32);
    /// Append whatever is still pending once the last bin has gone by.
    fn finish(&mut self, line: &mut String);
}

/// `chr\tstart\tend\tvalue`, one line per run of equal-valued bins — or per
/// bin, without `merge`.
#[derive(Default)]
struct BedGraphSink {
    chr: String,
    start: i64,
    end: i64,
    value: f32,
    /// Whether the three fields above stand for an interval yet. A bare `end`
    /// of 0 cannot say so: a first bin at `0-1` is a real interval.
    open: bool,
    merge: bool,
}

impl BedGraphSink {
    fn new(merge: bool) -> Self {
        Self {
            merge,
            ..Self::default()
        }
    }

    fn flush(&mut self, line: &mut String) {
        if !self.open {
            return;
        }
        let _ = write!(line, "{}\t{}\t{}\t", self.chr, self.start, self.end);
        super::text::push_float(line, self.value);
        line.push('\n');
        self.open = false;
    }
}

impl BinSink for BedGraphSink {
    fn bin(&mut self, line: &mut String, chr: &str, start: i64, end: i64, value: f32) {
        // `start == self.end` is what a skipped NaN bin breaks, and the chromosome
        // test is what a new chromosome starting at the old one's end would
        // otherwise slip past.
        if self.merge && self.open && self.end == start && self.value == value && self.chr == chr {
            self.end = end;
            return;
        }
        self.flush(line);
        self.chr.clear();
        self.chr.push_str(chr);
        self.start = start;
        self.end = end;
        self.value = value;
        self.open = true;
    }

    fn finish(&mut self, line: &mut String) {
        self.flush(line);
    }
}

/// fixedStep sections, one value per line.
#[derive(Default)]
struct WigSink {
    chr: String,
    span: i64,
    /// Where the next bin has to start for the open section to hold it.
    next_start: i64,
    open: bool,
}

impl BinSink for WigSink {
    fn bin(&mut self, line: &mut String, chr: &str, start: i64, end: i64, value: f32) {
        let span = end - start;
        if !self.open || span != self.span || start != self.next_start || self.chr != chr {
            // WIG coordinates are 1-based.
            let _ = writeln!(
                line,
                "fixedStep chrom={chr} start={} step={span} span={span}",
                start + 1
            );
            self.chr.clear();
            self.chr.push_str(chr);
            self.span = span;
            self.open = true;
        }
        super::text::push_float(line, value);
        line.push('\n');
        self.next_start = start + span;
    }

    /// Nothing is ever pending: a fixedStep section is closed by the next one
    /// or by the end of the file.
    fn finish(&mut self, _line: &mut String) {}
}

/// Write one formatted line, naming the output file if it fails.
fn write_line(
    writer: &mut std::io::BufWriter<std::fs::File>,
    line: &str,
    path: &std::path::Path,
) -> Result<()> {
    writer
        .write_all(line.as_bytes())
        .map_err(|e| Error::io(path.to_string_lossy(), e))
}

/// Flush explicitly rather than leaving it to `Drop`, which cannot report a
/// failure — and a truncated export that reported success is the worst outcome
/// here.
fn flush(writer: &mut std::io::BufWriter<std::fs::File>, path: &std::path::Path) -> Result<()> {
    writer
        .flush()
        .map_err(|e| Error::io(path.to_string_lossy(), e))
}