gwseq_io/bbi/

section.rs

//! The wig section buffer and its encoding heuristic.
//!
//! A bigWig stores values in sections, and each section is written in the
//! narrowest of three encodings that holds it: **fixedStep** while every value
//! shares one span and one step, four bytes a value; **variableStep** once the
//! starts turn irregular, eight; **bedGraph** once the spans differ too,
//! twelve. A section opens fixedStep and only ever widens.
//!
//! So every value that breaks the shape of the section it is being added to
//! poses the same question, and it is the only interesting question in this
//! file: **widen, or close the section here and open a narrower one?**
//!
//! - Widening costs the extra item width on everything already buffered *and*
//!   on everything the section goes on to hold, because it can never narrow
//!   again.
//! - Closing costs a section header, an R-tree leaf and a cold zlib stream —
//!   about 64 bytes — and leaves the next value to start a fresh section that
//!   may hit exactly the same problem.
//!
//! # What the answer is, and how it is known
//!
//! `examples/section_policy.rs` writes the same data through every rule in
//! [`SectionPolicy`] and reports the size of the file that comes out. That is
//! the only way to answer this: every rule produces a *valid* bigWig, so no
//! correctness test can distinguish them, and the difference between the best
//! and the worst is a factor of five.
//!
//! The measurement says two things.
//!
//! **The two obvious rules both fail catastrophically, in opposite
//! directions.** Always splitting is 356% worse on data that is irregular
//! throughout, because it writes one section per value. Always widening is 86%
//! worse on data that is regular with occasional breaks, because one odd value
//! makes the next thousand three times as wide. No fixed rule can do well on
//! both, and no arithmetic at the moment of the break can tell them apart:
//! the two cases look identical there.
//!
//! **What tells them apart is what happened last time.** A split that was
//! right is followed by a section that goes on to hold real work; a split that
//! was wrong is followed by one that breaks again immediately. So
//! [`SectionPolicy::Adaptive`] splits by default and stops as soon as two
//! sections in a row come out too short to be worth writing — and starts
//! again the moment one does not. Three details make it work, and each was
//! added because the measurement showed the file getting bigger without it:
//!
//! 1. **The tail is costed over the room left in the section**, not over what
//!    the current call happens to have left. A caller adding values one at a
//!    time has nothing left in the call, and the section still goes on filling
//!    in the wider encoding.
//! 2. **The streak only overrules the arithmetic for a section that is itself
//!    short.** A section holding a thousand uniform values is always worth
//!    closing, whatever the burst before it did.
//! 3. **A long run that cannot extend the open section flushes it first.**
//!    Nothing else in the writer is in a position to notice: after the first
//!    value widens a section, no *further* widening is needed, so the run
//!    pours in at bedGraph width and no decision is ever taken.
//!
//! Over six shapes of synthetic input that is **17% smaller** than the rule
//! this library shipped with, and best or within 0.05% of best on every one.
//! Over four real tracks it is best or tied on all four, and 0.8% smaller in
//! total. The two constants it adds are at the optimum and on a plateau: see
//! `--sweep`.
//!
//! # Reading a change to this file
//!
//! `section_counts` on the writer is the observable, and it is where to look
//! first: it reports how many sections went out in each encoding. A change
//! that makes files bigger shows up there before it shows up in a byte count.
//! `examples/section_policy.rs` is the check, and it takes seconds.

use crate::bbi::block::WigEncoding;
use crate::bbi::header::to_bbi_u32;
use crate::error::Result;

/// Hard ceiling on a section: its item count is a `u16` (Supp. Table 13).
pub const MAX_ITEMS_PER_SECTION: usize = 65535;

/// What closing a section early costs on disk: a fresh 24-byte section header,
/// a fresh 32-byte R-tree leaf, and the framing of a fresh zlib stream, since
/// Supp. Table 13 has every section compressed separately and so starting from
/// a cold window.
pub const SECTION_SPLIT_COST: i64 = 64;

/// How much of a widening survives deflate, as a percentage. The columns
/// varStep and bedGraph add are near-constant-stride coordinates, of which
/// deflate keeps very little, so widening a section costs far less on disk than
/// the difference in item size suggests.
pub const DOWNGRADE_COMPRESSION_PERCENT: i64 = 25;

/// Fewest items a section must hold to be worth writing on its own. Below
/// this, the 64 bytes of header, R-tree leaf and zlib framing outweigh what
/// the narrower encoding saves.
pub const MIN_SPLIT_ITEMS: i64 = 32;

/// How many too-short sections in a row [`SectionPolicy::Adaptive`] tolerates
/// before it stops splitting. Two: one runt can be a coincidence.
pub const RUNT_PATIENCE: u32 = 2;

/// How a section decides between widening its encoding and closing where it
/// stands.
///
/// Every variant produces a valid bigWig; they differ only in how big it is.
/// [`SectionPolicy::Adaptive`] is the default and the only one a caller should
/// normally set — the rest exist so that "is this rule any good?" is a
/// question with a measured answer rather than an opinion. See
/// `examples/section_policy.rs`, which writes the same data through each and
/// reports the bytes; the answer, over six shapes of input, is that `Adaptive`
/// writes files **17% smaller** than `Cost` and is the smallest or within half
/// a percent of it on every one.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SectionPolicy {
    /// Weigh the bytes a widening costs against the bytes a split costs,
    /// counting the widening over what is buffered plus what the current call
    /// still has to give.
    ///
    /// The rule this library shipped with, and measurably not a good one: a
    /// caller adding values one at a time has nothing left in the call, so the
    /// second term vanishes and the widening is costed over the prefix alone
    /// — while the section goes on filling in the wider encoding and paying
    /// for every value of it. Kept as the baseline the others are measured
    /// against.
    Cost,
    /// Never widen: close the section and start a new one in the narrower
    /// encoding. Most sections, smallest items.
    Split,
    /// Never split for encoding: widen and keep going until `items_per_slot`.
    /// Fewest sections, widest items.
    Widen,
    /// Weigh only what is already buffered, ignoring what the call still has
    /// to hand over. The obvious rule, and what the cost model would be
    /// without its second term.
    PrefixOnly,
    /// Weigh the widening against the **room left in the section**, rather
    /// than against what the current call happens to have left.
    ///
    /// [`SectionPolicy::Cost`] costs the tail at `min(batch_remaining, room)`,
    /// and a caller adding values one at a time has `batch_remaining == 0` —
    /// so the tail term vanishes and the widening is costed over the prefix
    /// alone. But the section does not stop there: it goes on filling to
    /// `items_per_slot` in the wider encoding, paying the extra width for
    /// every one of them. This costs what will actually be paid.
    Room,
    /// [`SectionPolicy::Room`], and never split a section too short to be
    /// worth writing.
    ///
    /// The failure `Room` has on its own is that data which is irregular
    /// *throughout* breaks the section immediately, every time: split, reopen,
    /// break at item two, split again — sections of one item, and a file
    /// several times larger than it should be. Widening a section that has
    /// barely started costs almost nothing, so below the floor it always
    /// widens.
    Runt,
    /// [`SectionPolicy::Room`], and split unless the data has just shown that
    /// splitting does not work.
    ///
    /// The two failure modes pull opposite ways and no arithmetic at the
    /// break can tell them apart, because both look identical at that instant:
    ///
    /// - Data that is **regular with rare breaks** wants a split. Widening
    ///   makes the whole rest of the section wider to absorb one odd value.
    /// - Data that is **irregular throughout** wants a widening. Splitting
    ///   emits one section per value, an R-tree leaf and a zlib stream each.
    ///
    /// What does tell them apart is what happened last time. A split that was
    /// right is followed by a section that goes on to hold many values; a split
    /// that was wrong is followed by one that breaks again immediately. So this
    /// splits by default and stops as soon as two sections in a row have come
    /// out too short to be worth writing — and starts again the moment one
    /// does not.
    ///
    /// The cost is bounded: at most `runt_patience` short sections are written
    /// before it adapts, once per stretch of irregular data.
    ///
    /// Two things beyond the streak make this the default rather than
    /// [`SectionPolicy::Runt`]:
    ///
    /// - The streak only overrules the arithmetic for a section that is
    ///   *itself* short. A section holding a thousand uniform values is always
    ///   worth closing, whatever the burst before it did.
    /// - A long run that cannot extend the open section flushes it first — see
    ///   `WigSection::should_flush_for_run`. Without that, one odd value
    ///   makes the next thousand cost three times as much, and nothing else in
    ///   the writer is in a position to notice.
    #[default]
    Adaptive,
}

/// The two constants the [`SectionPolicy::Cost`] model is built from, exposed
/// so they can be swept rather than argued about.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct CostModel {
    /// Bytes a split costs. 64 by default: a fresh section header, a fresh
    /// R-tree leaf and the framing of a fresh zlib stream.
    pub split_cost: i64,
    /// Percentage of a widening that survives deflate. 25 by default: the
    /// columns a wider encoding adds are near-constant-stride coordinates, of
    /// which deflate keeps very little.
    pub compression_percent: i64,
    /// Fewest items a section must hold to count as worth writing. Below this
    /// [`SectionPolicy::Runt`] will not close a section, and
    /// [`SectionPolicy::Adaptive`] counts one against its patience.
    pub min_split_items: i64,
    /// How many too-short sections in a row [`SectionPolicy::Adaptive`] writes
    /// before it concludes that splitting is not working here.
    pub runt_patience: u32,
}

impl Default for CostModel {
    fn default() -> Self {
        Self {
            split_cost: SECTION_SPLIT_COST,
            compression_percent: DOWNGRADE_COMPRESSION_PERCENT,
            min_split_items: MIN_SPLIT_ITEMS,
            runt_patience: RUNT_PATIENCE,
        }
    }
}

/// The 24-byte section header (Supp. Table 13).
pub const WIG_HEADER_SIZE: usize = 24;

/// Values buffered for the section being built, and what shape they have taken.
///
/// `starts` and `ends` are only materialised once the encoding needs them:
/// while the section is still fixedStep every start is `first_start + i * step`
/// and every end that start plus `span`, so a contiguous run of values goes in
/// as one extend rather than as three pushes per item.
#[derive(Debug, Default)]
pub struct WigSection {
    pub chr_ix: u32,
    pub encoding: Option<WigEncoding>,
    pub first_start: i64,
    /// Start delta shared by every adjacent pair, `None` while fewer than two
    /// items are buffered.
    pub uniform_step: Option<i64>,
    /// Width shared by every item, meaningless once the encoding has reached
    /// bedGraph.
    pub uniform_span: i64,
    pub last_start: i64,
    pub last_end: i64,
    pub starts: Vec<i64>,
    pub ends: Vec<i64>,
    pub values: Vec<f32>,
    items_per_slot: usize,
    policy: SectionPolicy,
    cost: CostModel,
    /// Sections closed in a row holding fewer than `min_split_items`, which is
    /// what [`SectionPolicy::Adaptive`] reads the shape of the data from.
    /// Counted at [`WigSection::clear`], the one place a section ends.
    runt_streak: u32,
}

impl WigSection {
    /// `items_per_slot`, and the rule the section decides widening with.
    /// See [`SectionPolicy`].
    pub fn with_policy(items_per_slot: usize, policy: SectionPolicy, cost: CostModel) -> Self {
        Self {
            items_per_slot,
            policy,
            cost,
            ..Default::default()
        }
    }

    pub fn len(&self) -> usize {
        self.values.len()
    }

    pub fn is_empty(&self) -> bool {
        self.values.is_empty()
    }

    pub fn is_full(&self) -> bool {
        self.values.len() >= self.items_per_slot
    }

    /// The encoding of a section that has items in it. A fresh section reports
    /// fixedStep, which is what it opens as.
    pub fn encoding(&self) -> WigEncoding {
        self.encoding.unwrap_or(WigEncoding::FixedStep)
    }

    pub fn clear(&mut self) {
        // The one place a section ends, so the one place to notice how it
        // went. A run of short ones is `Adaptive`'s signal that splitting is
        // not working on this data.
        //
        // A *long* section only clears that signal if it stayed fixedStep.
        // One that had to widen and then filled up is evidence the other way
        // — the data here is irregular and widening was the right call — so
        // clearing the streak on it would make the next section pay the same
        // couple of runt sections to learn the same thing, once per slot, for
        // as long as the irregularity lasts.
        if !self.is_empty() {
            if (self.len() as i64) < self.cost.min_split_items {
                self.runt_streak = self.runt_streak.saturating_add(1);
            } else if self.encoding() == WigEncoding::FixedStep {
                self.runt_streak = 0;
            }
        }
        // Not `*self = Self::new(..)`: the policy, the cost model and that
        // streak outlive a flush, as `items_per_slot` does.
        self.chr_ix = 0;
        self.encoding = None;
        self.first_start = -1;
        self.uniform_step = None;
        self.uniform_span = -1;
        self.last_start = -1;
        self.last_end = -1;
        // Cleared rather than replaced, so the capacity survives the flush and
        // a steady run of sections allocates nothing at all.
        self.starts.clear();
        self.ends.clear();
        self.values.clear();
    }

    /// Whether a section that has to widen from `old` to `new` should keep
    /// accumulating rather than close where it stands.
    ///
    /// Closing costs a section header, an R-tree leaf and a fresh zlib stream;
    /// widening costs the difference in item size over everything buffered, and
    /// over whatever is still to come in the call that broke the uniformity.
    /// Both are counted in bytes on disk, which is why the widening is
    /// discounted.
    fn should_widen(&self, old: WigEncoding, new: WigEncoding, batch_remaining: usize) -> bool {
        match self.policy {
            SectionPolicy::Split => return false,
            SectionPolicy::Widen => return true,
            SectionPolicy::Cost
            | SectionPolicy::PrefixOnly
            | SectionPolicy::Room
            | SectionPolicy::Runt
            | SectionPolicy::Adaptive => {}
        }
        let count = self.len() as i64;
        // A section too short to be worth writing is never worth closing: the
        // one it opens will meet the same value and close just as short.
        if self.policy == SectionPolicy::Runt && count < self.cost.min_split_items {
            return true;
        }
        // Splitting has already produced this many runts in a row, so it is not
        // going to work on this value either — *if* this section is itself
        // short. A section that has already collected real work is always
        // worth closing whatever the recent history: the streak is there to
        // stop a cascade of one-item sections, not to condemn a thousand
        // uniform values to bedGraph because the burst before them did not
        // split well.
        if self.policy == SectionPolicy::Adaptive
            && count < self.cost.min_split_items
            && self.runt_streak >= self.cost.runt_patience
        {
            return true;
        }
        let (old_size, new_size) = (old.item_size() as i64, new.item_size() as i64);
        let prefix_cost = (new_size - old_size) * count;
        // Against fixedStep, not against `old`: what the tail costs is the
        // whole widening from the narrowest form, since that is what those
        // items would otherwise have been written as in a section of their own.
        let room = self.items_per_slot as i64 - count;
        let widened_items = match self.policy {
            SectionPolicy::PrefixOnly => 0,
            SectionPolicy::Cost => (batch_remaining as i64).min(room),
            _ => room,
        };
        let tail_cost = (new_size - WigEncoding::FixedStep.item_size() as i64) * widened_items;
        (prefix_cost + tail_cost) * self.cost.compression_percent / 100 < self.cost.split_cost
    }

    /// Materialise the coordinate columns a wider encoding needs.
    fn widen(&mut self, target: WigEncoding) {
        let count = self.len();
        if self.encoding() == WigEncoding::FixedStep && target != WigEncoding::FixedStep {
            // With one item buffered there is no step yet, and the span stands
            // in for it — which is right, since the one start is `first_start`
            // either way.
            let step = if count >= 2 {
                self.uniform_step.unwrap_or(self.uniform_span)
            } else {
                self.uniform_span
            };
            self.starts.clear();
            self.starts
                .extend((0..count).map(|i| self.first_start + i as i64 * step));
            self.encoding = Some(WigEncoding::VarStep);
        }
        if self.encoding() == WigEncoding::VarStep && target == WigEncoding::BedGraph {
            self.ends.clear();
            let span = self.uniform_span;
            self.ends.extend(self.starts.iter().map(|s| s + span));
            self.encoding = Some(WigEncoding::BedGraph);
        }
    }

    fn push(&mut self, start: i64, end: i64, value: f32) {
        self.values.push(value);
        if self.encoding() != WigEncoding::FixedStep {
            self.starts.push(start);
        }
        if self.encoding() == WigEncoding::BedGraph {
            self.ends.push(end);
        }
        self.last_start = start;
        self.last_end = end;
    }

    /// Open an empty section on `chr_ix` at `start`, with `span`.
    fn open(&mut self, chr_ix: u32, start: i64, span: i64) {
        self.chr_ix = chr_ix;
        self.first_start = start;
        self.uniform_span = span;
        self.uniform_step = None;
        self.encoding = Some(WigEncoding::FixedStep);
    }

    /// Offer one item to the open section.
    ///
    /// `batch_remaining` is how many items the caller's current call still has
    /// to hand over after this one: the heuristic costs a widening over what is
    /// buffered *and* over what is still coming, which is what tells a one-off
    /// break of uniformity from the start of a long run of a different shape,
    /// and why a run handed over in one call encodes better than the same
    /// values one at a time.
    ///
    /// Returns [`Accept::Flush`] when the section has to be written before this
    /// item can go anywhere — the item is *not* buffered, and the caller offers
    /// it again once the section is empty.
    pub fn offer(
        &mut self,
        chr_ix: u32,
        start: i64,
        span: i64,
        value: f32,
        batch_remaining: usize,
    ) -> Accept {
        if !self.is_empty() && chr_ix != self.chr_ix {
            return Accept::Flush;
        }
        if self.is_empty() {
            self.open(chr_ix, start, span);
            self.push(start, start + span, value);
            return Accept::Buffered;
        }

        let count = self.len();
        let mut needed = WigEncoding::FixedStep;
        if span != self.uniform_span {
            needed = WigEncoding::BedGraph;
        } else if count >= 2 && Some(start - self.last_start) != self.uniform_step {
            needed = WigEncoding::VarStep;
        }
        // With a single item buffered the step is not defined yet, so any
        // second item of the same span still leaves the section fixedStep.
        if needed.item_size() > self.encoding().item_size() {
            if !self.should_widen(self.encoding(), needed, batch_remaining) {
                return Accept::Flush;
            }
            self.widen(needed);
        }
        if count == 1 && self.encoding() == WigEncoding::FixedStep {
            self.uniform_step = Some(start - self.first_start);
        }
        self.push(start, start + span, value);
        Accept::Buffered
    }

    /// Whether a long contiguous run should be given a section of its own
    /// rather than poured into the one that is open.
    ///
    /// The case this catches has nothing to do with the run and everything to
    /// do with what is already there. A run that **cannot extend** the open
    /// section — because the section has widened, or holds a different span,
    /// or stopped somewhere else — is handed over one value at a time, and
    /// every one of those values forces the section wider. So a burst of
    /// irregular values, or a single value of an odd width, makes the *next*
    /// thousand uniform values cost twelve bytes each instead of four. Nothing
    /// else in the writer can notice, because after the first of them no
    /// further widening is needed: bedGraph accepts anything.
    ///
    /// The arithmetic is [`Self::should_widen`]'s read the other way round.
    /// The run saves the width difference on every value it holds and costs
    /// one split.
    pub fn should_flush_for_run(&self, chr_ix: u32, start: i64, span: i64, run_len: usize) -> bool {
        if !matches!(self.policy, SectionPolicy::Adaptive) {
            return false;
        }
        // Nothing to flush, or the run extends what is open and costs nothing.
        if self.is_empty() || self.extends_run(chr_ix, start, span) {
            return false;
        }
        // A change of chromosome flushes anyway, further up.
        if self.chr_ix != chr_ix {
            return false;
        }
        // Poured in one at a time, the run would end up at bedGraph width: a
        // second span always forces it, whatever the section holds now.
        let poured = WigEncoding::BedGraph.item_size() as i64;
        let fresh = WigEncoding::FixedStep.item_size() as i64;
        let room = (self.items_per_slot - self.len()) as i64;
        let saved = (poured - fresh) * (run_len as i64).min(room);
        saved * self.cost.compression_percent / 100 > self.cost.split_cost
    }

    /// Whether a contiguous fixedStep run of `span` starting at `start` extends
    /// the open section rather than opening one.
    ///
    /// The fast path's precondition, ported from `write_values`: the section
    /// has to be fixedStep, on this chromosome, of this span, either holding a
    /// single item (whose step is not fixed yet) or already stepping by `span`,
    /// and ending exactly where the run begins.
    pub fn extends_run(&self, chr_ix: u32, start: i64, span: i64) -> bool {
        self.encoding() == WigEncoding::FixedStep
            && !self.is_empty()
            && self.chr_ix == chr_ix
            && self.uniform_span == span
            && (self.len() == 1 || self.uniform_step == Some(span))
            && self.last_end == start
    }

    /// Take a contiguous fixedStep run whole: no starts materialised, one
    /// extend. The fast path the API documentation points callers at.
    ///
    /// The caller must have checked [`Self::extends_run`], or the section must
    /// be empty. `start` is the first item's start and the run is `values.len()`
    /// items of `span` each.
    pub fn extend_run(&mut self, chr_ix: u32, start: i64, span: i64, values: &[f32]) {
        debug_assert!(!values.is_empty());
        if self.is_empty() {
            self.open(chr_ix, start, span);
        }
        // A run laid down this way is contiguous, so its step is its span; a
        // section holding one item has none yet, and this is what gives it one.
        self.uniform_step = Some(span);
        self.values.extend_from_slice(values);
        self.last_start = start + span * (values.len() as i64 - 1);
        self.last_end = start + span * values.len() as i64;
    }

    /// Encode to the section's chosen form, ready to compress.
    pub fn encode(&self) -> Result<Vec<u8>> {
        let count = self.len();
        let encoding = self.encoding();
        let mut out = Vec::with_capacity(WIG_HEADER_SIZE + count * encoding.item_size());

        out.extend_from_slice(&to_bbi_u32(self.chr_ix as i64, "chromId")?.to_le_bytes());
        out.extend_from_slice(&to_bbi_u32(self.first_start, "chromStart")?.to_le_bytes());
        // The end of the last item, which is not the end of the last step: a
        // section of 10 bp items spaced 100 bp apart stops 90 bp short of it.
        out.extend_from_slice(&to_bbi_u32(self.last_end, "chromEnd")?.to_le_bytes());
        let step = if encoding == WigEncoding::FixedStep {
            if count >= 2 {
                self.uniform_step.unwrap_or(self.uniform_span)
            } else {
                self.uniform_span
            }
        } else {
            0
        };
        out.extend_from_slice(&to_bbi_u32(step, "itemStep")?.to_le_bytes());
        let span = if encoding == WigEncoding::BedGraph {
            0
        } else {
            self.uniform_span
        };
        out.extend_from_slice(&to_bbi_u32(span, "itemSpan")?.to_le_bytes());
        out.push(encoding as u8);
        out.push(0); // reserved
        if count > MAX_ITEMS_PER_SECTION {
            return Err(crate::error::Error::invalid(format!(
                "wig section of {count} items exceeds the {MAX_ITEMS_PER_SECTION} its item \
                 count is stored on"
            )));
        }
        out.extend_from_slice(&(count as u16).to_le_bytes());
        debug_assert_eq!(out.len(), WIG_HEADER_SIZE);

        match encoding {
            WigEncoding::FixedStep => {
                for value in &self.values {
                    out.extend_from_slice(&value.to_le_bytes());
                }
            }
            WigEncoding::VarStep => {
                for i in 0..count {
                    out.extend_from_slice(&to_bbi_u32(self.starts[i], "chromStart")?.to_le_bytes());
                    out.extend_from_slice(&self.values[i].to_le_bytes());
                }
            }
            WigEncoding::BedGraph => {
                for i in 0..count {
                    out.extend_from_slice(&to_bbi_u32(self.starts[i], "chromStart")?.to_le_bytes());
                    out.extend_from_slice(&to_bbi_u32(self.ends[i], "chromEnd")?.to_le_bytes());
                    out.extend_from_slice(&self.values[i].to_le_bytes());
                }
            }
        }
        Ok(out)
    }
}

/// What the buffer decided to do with an item.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Accept {
    /// Buffered.
    Buffered,
    /// The section must be flushed first; the item is *not* buffered.
    Flush,
}

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

    /// The default policy, which is what a writer gets.
    fn section() -> WigSection {
        WigSection::with_policy(1024, SectionPolicy::default(), CostModel::default())
    }

    /// A section pinned to the cost rule, for the tests that are about *that*
    /// rule's arithmetic rather than about what the writer does.
    fn cost_section() -> WigSection {
        WigSection::with_policy(1024, SectionPolicy::Cost, CostModel::default())
    }

    /// A section that always widens, for the tests about what an encoding
    /// *writes* rather than about when it is chosen. Pinning the policy keeps
    /// them from breaking every time the heuristic is measured again.
    fn widening_section() -> WigSection {
        WigSection::with_policy(1024, SectionPolicy::Widen, CostModel::default())
    }

    /// Offer an item, flushing and retrying once if the section says to — which
    /// is what the writer does, and what makes the encoding the observable.
    fn offer(s: &mut WigSection, start: i64, span: i64, value: f32, remaining: usize) -> bool {
        match s.offer(0, start, span, value, remaining) {
            Accept::Buffered => false,
            Accept::Flush => {
                s.clear();
                assert_eq!(s.offer(0, start, span, value, remaining), Accept::Buffered);
                true
            }
        }
    }

    #[test]
    fn a_regular_run_stays_fixed_step() {
        let mut s = section();
        for i in 0..10 {
            assert!(!offer(&mut s, i * 10, 10, i as f32, 0));
        }
        assert_eq!(s.encoding(), WigEncoding::FixedStep);
        assert_eq!(s.len(), 10);
        // No coordinate column was ever materialised.
        assert!(s.starts.is_empty() && s.ends.is_empty());
        assert_eq!(s.uniform_step, Some(10));
    }

    #[test]
    fn a_gap_widens_to_var_step_when_the_widening_is_cheap() {
        let mut s = cost_section();
        for i in 0..4 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        // One item out of step, nothing else coming: 4 items * 4 extra bytes
        // is 16, a quarter of which is 4, well under the 64 a split costs.
        assert!(!offer(&mut s, 100, 10, 1.0, 0));
        assert_eq!(s.encoding(), WigEncoding::VarStep);
        assert_eq!(s.starts, [0, 10, 20, 30, 100]);
    }

    #[test]
    fn a_different_span_widens_all_the_way_to_bedgraph() {
        let mut s = cost_section();
        for i in 0..3 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        assert!(!offer(&mut s, 30, 25, 1.0, 0));
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
        assert_eq!(s.starts, [0, 10, 20, 30]);
        assert_eq!(s.ends, [10, 20, 30, 55]);
    }

    #[test]
    fn a_long_tail_of_the_new_shape_splits_the_section_instead() {
        // The whole point of `batch_remaining`: the same irregularity that is
        // absorbed above is a split when a long run of it is still coming.
        let mut s = cost_section();
        for i in 0..4 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        // Widening to bedGraph costs 8 more bytes on each of 4 buffered items
        // (32) plus 8 on each of the 500 still to come (4000); a quarter of
        // 4032 is far past 64.
        assert!(offer(&mut s, 30, 25, 1.0, 500));
        // Flushed, and the item opened a fresh fixedStep section of its own.
        assert_eq!(s.encoding(), WigEncoding::FixedStep);
        assert_eq!(s.len(), 1);
        assert_eq!(s.uniform_span, 25);
    }

    #[test]
    fn the_tail_is_costed_against_fixed_step_not_against_the_current_encoding() {
        // A varStep section widening to bedGraph. Costed against varStep the
        // tail would be 4 bytes an item; against fixedStep it is 8, and only
        // the second answer splits here. This is the term that is easiest to
        // get wrong and the one nothing else would catch.
        let mut s = cost_section();
        offer(&mut s, 0, 10, 1.0, 0);
        offer(&mut s, 10, 10, 1.0, 0);
        offer(&mut s, 40, 10, 1.0, 0); // irregular start -> varStep
        assert_eq!(s.encoding(), WigEncoding::VarStep);

        // 3 buffered * 4 extra = 12; tail of 30 items * 8 = 240; (12+240)/4 =
        // 63, which is under 64 — so it widens, by one byte of margin.
        let mut narrow = cost_section();
        for (start, span) in [(0, 10), (10, 10), (40, 10)] {
            offer(&mut narrow, start, span, 1.0, 0);
        }
        assert!(!offer(&mut narrow, 50, 25, 1.0, 30));
        assert_eq!(narrow.encoding(), WigEncoding::BedGraph);

        // One more item in the tail and it does not: (12 + 248) / 4 = 65.
        let mut wide = cost_section();
        for (start, span) in [(0, 10), (10, 10), (40, 10)] {
            offer(&mut wide, start, span, 1.0, 0);
        }
        assert!(offer(&mut wide, 50, 25, 1.0, 31));
    }

    // ---- the default policy ------------------------------------------
    //
    // Three behaviours, and each is the fix for a case the measured comparison
    // in `examples/section_policy.rs` showed the old cost model losing.

    /// Data that breaks every other value must not produce one section per
    /// value. After `runt_patience` short sections the policy stops splitting.
    #[test]
    fn a_cascade_of_short_sections_stops_itself() {
        let mut s = section();
        let mut flushes = 0;
        let mut at = 0i64;
        for i in 0..200 {
            let span = if i % 2 == 0 { 10 } else { 11 };
            if offer(&mut s, at, span, 1.0, 0) {
                flushes += 1;
            }
            at += span;
        }
        // Two runts to learn it, and then nothing: the section widened and
        // took the rest.
        assert!(
            flushes <= RUNT_PATIENCE as usize + 1,
            "{flushes} sections for 200 values"
        );
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
        assert!(s.len() > 190, "{} values buffered", s.len());
    }

    /// ...but the streak must not condemn a section that has real work in it.
    /// A burst of irregular values followed by a long uniform run has to split
    /// at the end of the run, not widen it.
    #[test]
    fn a_long_section_splits_however_bad_the_recent_history() {
        let mut s = section();
        // A burst, which teaches the policy that splitting is not working.
        let mut at = 0i64;
        for i in 0..12 {
            let span = 7 + (i % 5);
            offer(&mut s, at, span, 1.0, 0);
            at += span + 1;
        }
        assert!(s.encoding() != WigEncoding::FixedStep);
        s.clear();

        // Now a long uniform run, and then one value that breaks it.
        let mut at = 10_000i64;
        for _ in 0..500 {
            offer(&mut s, at, 10, 1.0, 0);
            at += 10;
        }
        assert_eq!(s.encoding(), WigEncoding::FixedStep);
        assert_eq!(s.len(), 500);
        // 500 buffered values are worth far more than one split, so this
        // closes the section rather than making all 500 twelve bytes wide.
        assert!(offer(&mut s, at, 25, 1.0, 0), "a 500-value section widened");
        assert_eq!(s.len(), 1);
    }

    /// A long run that cannot extend the open section gets one of its own.
    ///
    /// Nothing else in the writer can notice this: once a section is bedGraph
    /// no further widening is needed, so the run pours in at twelve bytes a
    /// value and no decision is ever taken.
    #[test]
    fn a_long_run_flushes_a_section_it_cannot_extend() {
        let mut s = section();
        // A widened section holding a handful of values.
        let mut at = 0i64;
        for i in 0..6 {
            let span = 7 + (i % 3);
            offer(&mut s, at, span, 1.0, 0);
            at += span + 2;
        }
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
        // A thousand uniform values are worth a split many times over.
        assert!(s.should_flush_for_run(0, at, 10, 1000));
        // A handful are not.
        assert!(!s.should_flush_for_run(0, at, 10, 4));
        // Nor is a run on another chromosome, which flushes for its own reason.
        assert!(!s.should_flush_for_run(1, at, 10, 1000));

        // And a run that simply extends what is open costs nothing to keep.
        let mut fixed = section();
        for i in 0..50 {
            offer(&mut fixed, i * 10, 10, 1.0, 0);
        }
        assert!(!fixed.should_flush_for_run(0, 500, 10, 1000));
    }

    /// The floor is what "short" means, and both users of it agree.
    #[test]
    fn the_runt_floor_is_read_from_the_cost_model() {
        let cost = CostModel {
            min_split_items: 4,
            ..Default::default()
        };
        let mut s = WigSection::with_policy(1024, SectionPolicy::Adaptive, cost);
        // Three values, then flush: short, so the streak rises.
        for i in 0..3 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        s.clear();
        for i in 0..3 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        s.clear();
        // Two runts in a row at a floor of 4: the policy has stopped splitting.
        let mut at = 0i64;
        offer(&mut s, at, 10, 1.0, 0);
        at += 10;
        assert!(
            !offer(&mut s, at, 25, 1.0, 0),
            "still splitting after 2 runts"
        );
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
    }

    #[test]
    fn the_tail_is_capped_by_the_room_left_in_the_section() {
        // A batch of a million does not cost a million: only what still fits.
        let mut s = WigSection::with_policy(8, SectionPolicy::Cost, CostModel::default());
        for i in 0..4 {
            offer(&mut s, i * 10, 10, 1.0, 0);
        }
        // Room for 4 more, so the tail is 4 * 8 = 32, plus 4 * 8 = 32
        // buffered; a quarter of 64 is 16, under the split cost.
        assert!(!offer(&mut s, 30, 25, 1.0, 1_000_000));
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
    }

    #[test]
    fn a_change_of_chromosome_always_flushes() {
        let mut s = section();
        s.offer(0, 0, 10, 1.0, 0);
        assert_eq!(s.offer(1, 0, 10, 1.0, 0), Accept::Flush);
    }

    #[test]
    fn a_second_item_of_the_same_span_stays_fixed_step_whatever_its_start() {
        // With one item buffered there is no step to break.
        let mut s = section();
        offer(&mut s, 0, 10, 1.0, 0);
        assert!(!offer(&mut s, 900, 10, 1.0, 1_000_000));
        assert_eq!(s.encoding(), WigEncoding::FixedStep);
        assert_eq!(s.uniform_step, Some(900));
    }

    #[test]
    fn widening_a_single_item_section_uses_its_span_as_the_step() {
        // There is no step yet, and the one start is `first_start` either way.
        let mut s = widening_section();
        offer(&mut s, 40, 10, 1.0, 0);
        offer(&mut s, 60, 25, 1.0, 0); // different span -> bedGraph
        assert_eq!(s.encoding(), WigEncoding::BedGraph);
        assert_eq!(s.starts, [40, 60]);
        assert_eq!(s.ends, [50, 85]);
    }

    #[test]
    fn a_run_extends_the_open_section_only_when_it_lines_up() {
        let mut s = section();
        s.extend_run(0, 0, 10, &[1.0, 2.0, 3.0]);
        assert_eq!(s.len(), 3);
        assert_eq!(s.last_end, 30);
        assert!(s.extends_run(0, 30, 10));
        assert!(!s.extends_run(0, 40, 10), "a gap does not extend");
        assert!(!s.extends_run(0, 30, 20), "a different span does not");
        assert!(!s.extends_run(1, 30, 10), "another chromosome does not");
        s.extend_run(0, 30, 10, &[4.0]);
        assert_eq!(s.len(), 4);
        assert_eq!(s.encoding(), WigEncoding::FixedStep);
        assert!(s.starts.is_empty(), "the fast path materialises no starts");
    }

    #[test]
    fn a_single_item_section_is_extended_by_a_run_that_meets_it() {
        // `len() == 1` has no step yet, so a run of the same span extends it.
        let mut s = section();
        offer(&mut s, 0, 10, 1.0, 0);
        assert!(s.extends_run(0, 10, 10));
        s.extend_run(0, 10, 10, &[2.0, 3.0]);
        assert_eq!(s.uniform_step, Some(10));
        assert_eq!(s.len(), 3);
    }

    // ---- encoding --------------------------------------------------------

    fn header_of(bytes: &[u8]) -> (u32, u32, u32, u32, u32, u8, u16) {
        let u32_at = |o: usize| u32::from_le_bytes(bytes[o..o + 4].try_into().unwrap());
        (
            u32_at(0),
            u32_at(4),
            u32_at(8),
            u32_at(12),
            u32_at(16),
            bytes[20],
            u16::from_le_bytes(bytes[22..24].try_into().unwrap()),
        )
    }

    #[test]
    fn a_fixed_step_section_encodes_to_a_header_and_a_column_of_values() {
        let mut s = section();
        for i in 0..3 {
            offer(&mut s, 100 + i * 10, 10, i as f32, 0);
        }
        let bytes = s.encode().unwrap();
        assert_eq!(bytes.len(), WIG_HEADER_SIZE + 3 * 4);
        let (chr, start, end, step, span, kind, count) = header_of(&bytes);
        assert_eq!((chr, start, end), (0, 100, 130));
        assert_eq!((step, span), (10, 10));
        assert_eq!(kind, WigEncoding::FixedStep as u8);
        assert_eq!(count, 3);
    }

    #[test]
    fn a_sparse_fixed_step_section_ends_at_its_last_item_not_its_last_step() {
        // 10 bp items spaced 100 apart: chromEnd is 210, not 300.
        let mut s = section();
        for i in 0..3 {
            offer(&mut s, i * 100, 10, 1.0, 0);
        }
        let (_, start, end, step, span, _, _) = header_of(&s.encode().unwrap());
        assert_eq!((start, end), (0, 210));
        assert_eq!((step, span), (100, 10));
    }

    #[test]
    fn a_var_step_section_writes_starts_and_a_span_but_no_step() {
        let mut s = widening_section();
        offer(&mut s, 0, 10, 1.0, 0);
        offer(&mut s, 10, 10, 2.0, 0);
        offer(&mut s, 40, 10, 3.0, 0);
        let bytes = s.encode().unwrap();
        assert_eq!(bytes.len(), WIG_HEADER_SIZE + 3 * 8);
        let (_, _, end, step, span, kind, count) = header_of(&bytes);
        assert_eq!(kind, WigEncoding::VarStep as u8);
        assert_eq!((step, span, end, count), (0, 10, 50, 3));
        assert_eq!(u32::from_le_bytes(bytes[24..28].try_into().unwrap()), 0);
        assert_eq!(u32::from_le_bytes(bytes[32..36].try_into().unwrap()), 10);
        assert_eq!(u32::from_le_bytes(bytes[40..44].try_into().unwrap()), 40);
    }

    #[test]
    fn a_bedgraph_section_writes_both_coordinates_and_neither_step_nor_span() {
        let mut s = widening_section();
        offer(&mut s, 0, 10, 1.0, 0);
        offer(&mut s, 10, 25, 2.0, 0);
        let bytes = s.encode().unwrap();
        assert_eq!(bytes.len(), WIG_HEADER_SIZE + 2 * 12);
        let (_, _, end, step, span, kind, count) = header_of(&bytes);
        assert_eq!(kind, WigEncoding::BedGraph as u8);
        assert_eq!((step, span, end, count), (0, 0, 35, 2));
    }

    #[test]
    fn a_single_item_section_writes_its_span_as_its_step() {
        // There is no step, and a reader multiplies the step by the index —
        // which for the one item is zero either way.
        let mut s = section();
        offer(&mut s, 60, 15, 1.0, 0);
        let (_, start, end, step, span, _, count) = header_of(&s.encode().unwrap());
        assert_eq!((start, end, step, span, count), (60, 75, 15, 15, 1));
    }

    #[test]
    fn an_encoded_section_reads_back_through_the_reader() {
        // The round trip that matters: the writer's own reader has to decode
        // what it wrote, item for item.
        for shape in 0..3 {
            let mut s = section();
            match shape {
                0 => {
                    for i in 0..5 {
                        offer(&mut s, i * 10, 10, i as f32 * 1.5, 0);
                    }
                }
                1 => {
                    offer(&mut s, 0, 10, 1.0, 0);
                    offer(&mut s, 10, 10, 2.0, 0);
                    offer(&mut s, 45, 10, 3.0, 0);
                }
                _ => {
                    offer(&mut s, 0, 10, 1.0, 0);
                    offer(&mut s, 10, 25, 2.0, 0);
                    offer(&mut s, 60, 5, 3.0, 0);
                }
            }
            let want: Vec<(i64, i64, f32)> = (0..s.len())
                .map(|i| {
                    let start = if s.encoding() == WigEncoding::FixedStep {
                        s.first_start + i as i64 * s.uniform_step.unwrap_or(s.uniform_span)
                    } else {
                        s.starts[i]
                    };
                    let end = if s.encoding() == WigEncoding::BedGraph {
                        s.ends[i]
                    } else {
                        start + s.uniform_span
                    };
                    (start, end, s.values[i])
                })
                .collect();

            let bytes = bytes::Bytes::from(s.encode().unwrap());
            let header = crate::bbi::block::read_wig_header(&bytes, "test.bigwig").unwrap();
            assert_eq!(header.item_count as usize, s.len(), "shape {shape}");
            let got: Vec<(i64, i64, f32)> = (0..header.item_count as usize)
                .map(|i| {
                    let item = crate::bbi::block::read_wig_item(&bytes, &header, i, "test.bigwig")
                        .unwrap();
                    (item.start, item.end, item.value)
                })
                .collect();
            assert_eq!(got, want, "shape {shape}");
        }
    }

    #[test]
    fn a_coordinate_past_32_bits_is_refused_rather_than_truncated() {
        let mut s = section();
        offer(&mut s, 5_000_000_000, 10, 1.0, 0);
        let err = s.encode().unwrap_err().to_string();
        assert!(err.contains("chromStart 5000000000"), "{err}");
    }
}