gwseq_io/genomic/

locs.rs

//! Loci: how a request names the regions it wants, and what a reader turns
//! that into.
//!
//! [`Locs::parse`] takes `starts`, `ends`, `centers` and `span` as optional
//! vectors and works out at runtime which combination was meant. It exists for
//! the Python layer, which takes all five keywords. Rust callers get the four
//! legal combinations as four constructors instead, so an illegal one cannot
//! be spelled.

use crate::error::{Error, Result};
use crate::genomic::ChrMap;

#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Strand {
    #[default]
    Forward,
    Reverse,
}

impl Strand {
    #[inline]
    pub fn is_reverse(self) -> bool {
        matches!(self, Strand::Reverse)
    }
}

impl std::str::FromStr for Strand {
    type Err = Error;

    /// `"+"`, `"-"`, and the `"."` a BED writes for an unstranded feature and
    /// the empty string, both of which count as forward.
    fn from_str(s: &str) -> Result<Self> {
        match s {
            "+" | "." | "" => Ok(Strand::Forward),
            "-" => Ok(Strand::Reverse),
            other => Err(Error::invalid(format!("Strand {other} invalid (+ or -)"))),
        }
    }
}

/// A request's regions, before they are resolved against a file.
#[derive(Debug, Clone, Default, PartialEq)]
pub struct Locs {
    pub chr_ids: Vec<String>,
    pub starts: Vec<i64>,
    pub ends: Vec<i64>,
    /// One per locus, or empty for all-forward.
    pub strands: Vec<Strand>,
}

impl Locs {
    /// The runtime-dispatched form, for the Python layer.
    ///
    /// Accepts one of: `starts` + `ends` with no span; or exactly one of
    /// `starts` / `ends` / `centers` together with `span`. `span` of `None`
    /// means no span was given, and so does `Some(-1)` — the Python layer has
    /// no `Option` to pass and uses -1 as the sentinel. Anything below that is
    /// refused: a caller reaching this from Rust with `Some(-5)` means
    /// something, and silently reading it as "no span" is not it.
    ///
    /// No loci at all is a request for nothing rather than a malformed one —
    /// every reader handles an empty batch — and that is checked first, because
    /// the tests below read an empty `starts`/`ends` as "not given".
    pub fn parse(
        chr_ids: &[String],
        starts: &[i64],
        ends: &[i64],
        centers: &[i64],
        span: Option<i64>,
    ) -> Result<Self> {
        if chr_ids.is_empty() && starts.is_empty() && ends.is_empty() && centers.is_empty() {
            return Ok(Self::default());
        }
        if span.is_some_and(|span| span < -1) {
            return Err(Error::invalid(format!(
                "span {} is negative (-1, or no span at all, means none was given)",
                span.expect("checked just above")
            )));
        }

        let (parsed_starts, parsed_ends) = match span {
            Some(span) if span >= 0 => {
                let given = [!starts.is_empty(), !ends.is_empty(), !centers.is_empty()];
                if given.iter().filter(|g| **g).count() != 1 {
                    return Err(Error::invalid(
                        "Exactly one of starts, ends or centers must be specified when using span.",
                    ));
                }
                if !starts.is_empty() {
                    (starts.to_vec(), starts.iter().map(|s| s + span).collect())
                } else if !ends.is_empty() {
                    (ends.iter().map(|e| e - span).collect(), ends.to_vec())
                } else {
                    // Asymmetric on an odd span, and the asymmetry is expected:
                    // [c - span/2, c + (span + 1)/2). Rust's `/` truncates
                    // toward zero, so a negative centre rounds the same way.
                    (
                        centers.iter().map(|c| c - span / 2).collect(),
                        centers.iter().map(|c| c + (span + 1) / 2).collect(),
                    )
                }
            }
            _ => {
                if starts.is_empty() || ends.is_empty() {
                    return Err(Error::invalid(
                        "Either starts+ends or exactly one of starts/ends/centers \
                         together with span must be specified.",
                    ));
                }
                (starts.to_vec(), ends.to_vec())
            }
        };

        if chr_ids.len() != parsed_starts.len() || chr_ids.len() != parsed_ends.len() {
            return Err(Error::invalid(format!(
                "Length mismatch between chr_ids ({}) and starts/ends/centers ({}/{})",
                chr_ids.len(),
                parsed_starts.len(),
                parsed_ends.len()
            )));
        }

        Ok(Self {
            chr_ids: chr_ids.to_vec(),
            starts: parsed_starts,
            ends: parsed_ends,
            strands: Vec::new(),
        })
    }

    /// `starts` and `ends` given directly.
    pub fn spans(chr_ids: &[String], starts: &[i64], ends: &[i64]) -> Result<Self> {
        Self::parse(chr_ids, starts, ends, &[], None)
    }

    /// `ends[i] = starts[i] + span`.
    pub fn from_starts(chr_ids: &[String], starts: &[i64], span: i64) -> Result<Self> {
        Self::parse(chr_ids, starts, &[], &[], Some(span))
    }

    /// `starts[i] = ends[i] - span`.
    pub fn from_ends(chr_ids: &[String], ends: &[i64], span: i64) -> Result<Self> {
        Self::parse(chr_ids, &[], ends, &[], Some(span))
    }

    /// `[c - span/2, c + (span + 1)/2)`.
    pub fn centered(chr_ids: &[String], centers: &[i64], span: i64) -> Result<Self> {
        Self::parse(chr_ids, &[], &[], centers, Some(span))
    }

    /// Chromosome names and nothing else.
    ///
    /// What the whole-file calls take: they resolve the regions themselves from
    /// the file's chromosome sizes, so there are no coordinates to give.
    /// [`Locs::parse`] refuses this shape — with no `starts` and no `span` there
    /// is nothing for it to build — which is right for a request that names
    /// windows and wrong for one that names chromosomes.
    pub fn chromosomes(chr_ids: Vec<String>) -> Self {
        Self {
            chr_ids,
            starts: Vec::new(),
            ends: Vec::new(),
            strands: Vec::new(),
        }
    }

    /// Every named chromosome end to end, for the `read_all_*` and `iter_all_*`
    /// paths. An empty `chr_ids` means all of them, in the file's order.
    pub fn whole_chromosomes(map: &ChrMap, chr_ids: &[String]) -> Result<Self> {
        let selected = map.select(chr_ids)?;
        Ok(Self {
            chr_ids: selected.iter().map(|e| e.id.clone()).collect(),
            starts: vec![0; selected.len()],
            ends: selected.iter().map(|e| e.size).collect(),
            strands: Vec::new(),
        })
    }

    /// Empty leaves every locus forward, which is what
    /// a request that named no strand gets; anything else must be one entry per
    /// locus.
    pub fn with_strands(mut self, strands: &[String]) -> Result<Self> {
        if strands.is_empty() {
            self.strands = Vec::new();
            return Ok(self);
        }
        if strands.len() != self.chr_ids.len() {
            return Err(Error::invalid(format!(
                "Length mismatch between strands ({}) and the {} loci of the request",
                strands.len(),
                self.chr_ids.len()
            )));
        }
        self.strands = strands
            .iter()
            .enumerate()
            .map(|(i, s)| {
                s.parse::<Strand>().map_err(|_| {
                    Error::invalid(format!("Strand {s} at index {i} invalid (+ or -)"))
                })
            })
            .collect::<Result<_>>()?;
        Ok(self)
    }

    /// Strand of locus `i`, forward when none were given.
    #[inline]
    pub fn strand(&self, i: usize) -> Strand {
        self.strands.get(i).copied().unwrap_or_default()
    }

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

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

/// A locus resolved against a file and against the request's binning: it knows
/// its chromosome's index, the bin grid it snaps to, and the slice of the
/// output row it owns.
#[derive(Debug, Clone)]
pub struct IndexedLoc {
    /// The chromosome's index **as the file numbers it**, which is what data
    /// records and R-tree items carry.
    pub chr_index: usize,
    pub start: i64,
    pub end: i64,
    /// Start and end snapped to the bin grid: the start always down, the end
    /// down as well unless `full_bin`, which rounds it up so every bin the
    /// locus touches is covered.
    pub binned_start: i64,
    pub binned_end: i64,
    /// Width of *this* locus's output bins, which is its binned span over the
    /// request's `bin_count` — equal to the request's `bin_size` only when the
    /// two happen to divide.
    pub bin_size: f64,
    pub reverse: bool,
    /// Half-open slice of the flat output buffer this locus fills, in the order
    /// the loci were **given** rather than the order they are read in. This is
    /// what carries the requested order through the sort into file order.
    pub output_start: usize,
    pub output_end: usize,
}

impl IndexedLoc {
    /// Bases the window reads, i.e. the span its bins share out.
    #[inline]
    pub fn span(&self) -> i64 {
        self.binned_end - self.binned_start
    }

    /// Output bins the window fills, which is the request's `bin_count`.
    #[inline]
    pub fn bin_count(&self) -> i64 {
        (self.output_end - self.output_start) as i64
    }

    /// Output bin the base at `pos` falls in.
    ///
    /// Exact integer arithmetic rather than `pos / bin_size`, whose `f64`
    /// quotient is not representable for most pairs. A base on an interior edge
    /// then landed on either side of it depending on the rounding, and
    /// [`bin_after`](Self::bin_after) disagreed with this floor about where the
    /// edge was — so a value one bin wide was counted into two.
    #[inline]
    pub fn bin_at(&self, pos: i64) -> i64 {
        ((pos - self.binned_start) * self.bin_count()) / self.span()
    }

    /// One past the last output bin the half-open range ending at `pos` reaches.
    #[inline]
    pub fn bin_after(&self, pos: i64) -> i64 {
        let span = self.span();
        ((pos - self.binned_start) * self.bin_count() + span - 1) / span
    }

    /// Bases of `[from, to)` that fall in output bin `b`, 0 when none do.
    ///
    /// Fractional: a bin is `span / bin_count` bases wide and drops below one
    /// base whenever more bins are asked for than the window has. Computed in
    /// units of one `bin_count`th of a base, where every bound is whole, so the
    /// result is exact rather than accumulated from a rounded bin width.
    ///
    /// The products are bounded by `span * bin_count` — 2.5e14 for a human
    /// chromosome in a million bins, against the 9.2e18 an `i64` holds.
    #[inline]
    pub fn bin_coverage(&self, b: i64, from: i64, to: i64) -> f64 {
        let n = self.bin_count();
        let s = self.span();
        let lo = ((from - self.binned_start) * n).max(b * s);
        let hi = ((to - self.binned_start) * n).min((b + 1) * s);
        if hi > lo {
            (hi - lo) as f64 / n as f64
        } else {
            0.0
        }
    }

    /// How much of output bin `b` the range `[from, to)` covers, from 0 to 1.
    ///
    /// [`bin_coverage`](Self::bin_coverage) in bins rather than bases, so that a
    /// pileup reports the depth over the bin however wide it is — an entry
    /// spanning the whole bin is a depth of 1, one covering a third of it a
    /// third.
    #[inline]
    pub fn bin_fraction(&self, b: i64, from: i64, to: i64) -> f64 {
        let n = self.bin_count();
        let s = self.span();
        let lo = ((from - self.binned_start) * n).max(b * s);
        let hi = ((to - self.binned_start) * n).min((b + 1) * s);
        if hi > lo {
            (hi - lo) as f64 / s as f64
        } else {
            0.0
        }
    }

    /// Row of this locus in the request, i.e. where its per-locus result
    /// belongs.
    ///
    /// `output_start` counts output *bins* rather than loci, so it is only the
    /// row itself while there is one bin per locus — which is what a
    /// quantification of a bigWig asks for, and no longer what one of a bigBed
    /// does.
    #[inline]
    pub fn row(&self, bin_count: usize) -> usize {
        self.output_start.checked_div(bin_count).unwrap_or(0)
    }
}

/// Windows of a bbi request, ordered by chromosome index then position.
#[derive(Debug, Clone, Default)]
pub struct IndexedLocs {
    /// Sorted by (chromosome index, binned start, binned end) — the order a bbi
    /// file is best read in. Each carries the slice of the output buffer it
    /// fills, which keeps the requested order.
    pub locs: Vec<IndexedLoc>,
    /// Output bins per locus, inferred from the widest locus when it was not
    /// requested.
    pub bin_count: usize,
    /// Size of the output buffer the windows fill.
    pub output_len: usize,
}

/// A contiguous run of loci one worker owns.
///
/// Half-open over [`IndexedLocs::locs`].
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct LocBatch {
    pub start: usize,
    pub end: usize,
}

impl LocBatch {
    #[inline]
    pub fn is_empty(&self) -> bool {
        self.start >= self.end
    }
    #[inline]
    pub fn len(&self) -> usize {
        self.end.saturating_sub(self.start)
    }
}

impl IndexedLocs {
    /// Resolve, bin and sort a request's loci against a file.
    ///
    /// Coordinates are left as they were asked for, so a window may reach
    /// outside its chromosome. Filling the bins no data covers is the reader's.
    pub fn build(
        map: &ChrMap,
        locs: &Locs,
        bin_size: f64,
        bin_count: Option<usize>,
        full_bin: bool,
    ) -> Result<Self> {
        // One place decides what a valid bin size is, and this is the door
        // every read goes through. `BinPlan::new` refuses a fractional one,
        // which is what lets the snapping below be integer arithmetic.
        let bin = crate::genomic::BinPlan::new(bin_size, bin_count, full_bin)?.whole_bin_size();
        if locs.is_empty() {
            return Ok(Self {
                locs: Vec::new(),
                bin_count: bin_count.unwrap_or(0),
                output_len: 0,
            });
        }

        let mut built = Vec::with_capacity(locs.len());
        let mut max_binned_span = 0i64;
        for i in 0..locs.len() {
            let entry = map.resolve(&locs.chr_ids[i])?;
            let (start, end) = (locs.starts[i], locs.ends[i]);
            if start > end {
                return Err(Error::invalid(format!(
                    "Locus {}:{start}-{end} at index {i} ends before it starts",
                    locs.chr_ids[i]
                )));
            }
            // Snapped to the genome-wide grid of `bin` bases, in integers.
            // Through `f64` this rounded: a coordinate past 2^53 is not
            // exactly representable, and `(start / b).floor() * b` then lands
            // a base or two off — which for a human chromosome is fine and for
            // a scaffolded plant genome is not.
            //
            // The window moves to the grid, it does not stay where it was
            // asked: `chr1:1010-2010` at `bin_size=100` reads `chr1:1000-2000`.
            // That is deliberate — it is what makes bin `n` mean the same
            // bases for every locus in a request, which is what a profile
            // needs — and it is documented under `bin_size` in the README.
            let binned_start = start.div_euclid(bin) * bin;
            let binned_end = if full_bin {
                end.saturating_add(bin - 1).div_euclid(bin) * bin
            } else {
                end.div_euclid(bin) * bin
            };
            max_binned_span = max_binned_span.max(binned_end - binned_start);
            built.push(IndexedLoc {
                chr_index: entry.index,
                start,
                end,
                binned_start,
                binned_end,
                bin_size: 0.0,
                reverse: locs.strand(i).is_reverse(),
                output_start: 0,
                output_end: 0,
            });
        }

        let bin_count = match bin_count {
            Some(n) => n,
            None => (max_binned_span / bin) as usize,
        };
        if bin_count < 1 {
            return Err(Error::invalid(format!(
                "No locus is a full bin wide (widest binned locus: {max_binned_span} bp, \
                 bin size: {bin_size})"
            )));
        }

        for (i, loc) in built.iter_mut().enumerate() {
            loc.bin_size = (loc.binned_end - loc.binned_start) as f64 / bin_count as f64;
            loc.output_start = i * bin_count;
            loc.output_end = loc.output_start + bin_count;
        }
        let output_len = built.len() * bin_count;

        // Stable, so two loci with identical bounds keep the order they were
        // asked in — which decides nothing about the data but keeps a request
        // reproducible.
        built.sort_by(|a, b| {
            (a.chr_index, a.binned_start, a.binned_end).cmp(&(
                b.chr_index,
                b.binned_start,
                b.binned_end,
            ))
        });

        Ok(Self {
            locs: built,
            bin_count,
            output_len,
        })
    }

    /// Base pairs the windows read: the sum of their binned spans.
    ///
    /// The work the request represents, not the region it spans — overlapping
    /// windows count their shared bases once per window, as the reader reads
    /// them once per window.
    pub fn coverage(&self) -> u64 {
        self.locs
            .iter()
            .map(|l| (l.binned_end - l.binned_start).max(0) as u64)
            .sum()
    }

    /// The narrowest output bin the request ends up with — what a zoom level
    /// has to be finer than.
    ///
    /// Not the `bin_size` argument: that is the grid the loci are snapped to,
    /// and stays at its default when the caller asked for a `bin_count`
    /// instead. The narrowest locus governs, since one narrower than the widest
    /// gets proportionally narrower bins, and reading it off a coarser level
    /// would give it the edges of the summaries rather than an average of them.
    ///
    /// `fallback` only stands in for a request that produced no locus at all.
    pub fn effective_bin_size(&self, fallback: f64) -> f64 {
        let narrowest = self
            .locs
            .iter()
            .map(|l| l.bin_size)
            .filter(|s| *s > 0.0)
            .fold(f64::INFINITY, f64::min);
        if narrowest.is_finite() {
            narrowest
        } else {
            fallback
        }
    }

    /// Split into `parallel` batches of roughly equal *coverage*, so workers
    /// reading a batch each have a comparable amount of work.
    ///
    /// Coverage-balanced, not count-balanced: one 10 Mb locus and a thousand
    /// 100 bp ones are not the same amount of work. The split is behaviour — it
    /// decides which blocks each worker touches, and so both how much of the
    /// cache is shared and what order the `f32` accumulators are summed in.
    /// Changing it changes the last bits of every mean.
    ///
    /// A window is never split, so fewer batches than asked for come back when
    /// the windows cannot be spread that thin, and never more. Batches stay
    /// contiguous, which keeps each worker reading the file forward.
    pub fn batches(&self, parallel: usize) -> (Vec<LocBatch>, u64) {
        let total = self.coverage();
        if self.locs.is_empty() {
            return (Vec::new(), total);
        }
        let parallel = parallel.max(1) as u64;
        // Rounded up, so filling every batch to it leaves no batch in excess.
        let per_batch = total.div_ceil(parallel).max(1);

        let mut batches = vec![LocBatch { start: 0, end: 0 }];
        let mut coverage = 0u64;
        for i in 0..self.locs.len().saturating_sub(1) {
            let loc = &self.locs[i];
            coverage += (loc.binned_end - loc.binned_start).max(0) as u64;
            if coverage < per_batch {
                continue;
            }
            // Windows the grid left no wider than a bin read nothing, and a run
            // of them at the end would otherwise open a batch beyond `parallel`
            // carrying no work. Once the count is reached everything left joins
            // the last batch, which is where rounding `per_batch` up already
            // expects the remainder to go.
            if batches.len() as u64 >= parallel {
                continue;
            }
            batches.last_mut().unwrap().end = i + 1;
            batches.push(LocBatch {
                start: i + 1,
                end: i + 1,
            });
            coverage = 0;
        }
        batches.last_mut().unwrap().end = self.locs.len();
        (batches, total)
    }

    /// Mirror the output slice of every reverse-strand locus, so its bins run
    /// from its end to its start.
    ///
    /// A whole row is turned around here rather than inside the extraction
    /// loop, where the mirror would be a branch per bin of a walk that has
    /// none. The file is read forward whatever the strand.
    pub fn reverse_output_rows(&self, output: &mut [f32]) {
        for loc in &self.locs {
            if loc.reverse {
                output[loc.output_start..loc.output_end].reverse();
            }
        }
    }
}

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

    fn ids(n: usize) -> Vec<String> {
        (0..n).map(|_| "chr1".to_string()).collect()
    }

    #[test]
    fn starts_and_ends_pass_through() {
        let l = Locs::spans(&ids(2), &[10, 20], &[15, 25]).unwrap();
        assert_eq!(l.starts, [10, 20]);
        assert_eq!(l.ends, [15, 25]);
    }

    #[test]
    fn span_extends_from_starts_and_back_from_ends() {
        assert_eq!(
            Locs::from_starts(&ids(2), &[10, 20], 5).unwrap().ends,
            [15, 25]
        );
        assert_eq!(
            Locs::from_ends(&ids(2), &[15, 25], 5).unwrap().starts,
            [10, 20]
        );
    }

    #[test]
    fn centered_windows_are_asymmetric_on_an_odd_span() {
        let l = Locs::centered(&ids(1), &[100], 5).unwrap();
        assert_eq!((l.starts[0], l.ends[0]), (98, 103));
        let l = Locs::centered(&ids(1), &[100], 4).unwrap();
        assert_eq!((l.starts[0], l.ends[0]), (98, 102));
    }

    #[test]
    fn no_loci_at_all_is_a_request_for_nothing() {
        let l = Locs::parse(&[], &[], &[], &[], None).unwrap();
        assert!(l.is_empty());
    }

    #[test]
    fn ambiguous_and_missing_combinations_are_refused() {
        // Two of the three with a span.
        let e = Locs::parse(&ids(1), &[1], &[2], &[], Some(5))
            .unwrap_err()
            .to_string();
        assert!(e.contains("Exactly one of starts, ends or centers"), "{e}");
        // Centers with no span.
        let e = Locs::parse(&ids(1), &[], &[], &[5], None)
            .unwrap_err()
            .to_string();
        assert!(e.contains("Either starts+ends"), "{e}");
        // Starts without ends and without a span.
        assert!(Locs::parse(&ids(1), &[1], &[], &[], None).is_err());
    }

    #[test]
    fn length_mismatch_names_both_counts() {
        let e = Locs::spans(&ids(3), &[1, 2], &[2, 3])
            .unwrap_err()
            .to_string();
        assert!(e.contains("chr_ids (3)") && e.contains("(2/2)"), "{e}");
    }

    #[test]
    fn strands_default_forward_and_accept_dot_and_empty() {
        let l = Locs::spans(&ids(3), &[0, 0, 0], &[1, 1, 1]).unwrap();
        assert_eq!(l.strand(0), Strand::Forward);
        let l = l
            .with_strands(&["+".into(), "-".into(), ".".into()])
            .unwrap();
        assert_eq!(l.strand(1), Strand::Reverse);
        assert_eq!(l.strand(2), Strand::Forward);
        let l = Locs::spans(&ids(1), &[0], &[1])
            .unwrap()
            .with_strands(&["".into()])
            .unwrap();
        assert_eq!(l.strand(0), Strand::Forward);
    }

    #[test]
    fn a_bad_strand_names_its_index() {
        let e = Locs::spans(&ids(2), &[0, 0], &[1, 1])
            .unwrap()
            .with_strands(&["+".into(), "?".into()])
            .unwrap_err()
            .to_string();
        assert!(e.contains("Strand ? at index 1"), "{e}");
    }

    #[test]
    fn strand_count_must_match_loci() {
        let e = Locs::spans(&ids(2), &[0, 0], &[1, 1])
            .unwrap()
            .with_strands(&["+".into()])
            .unwrap_err()
            .to_string();
        assert!(e.contains("strands (1)") && e.contains("2 loci"), "{e}");
    }

    #[test]
    fn chromosomes_names_without_coordinates() {
        let l = Locs::chromosomes(vec!["chr1".into(), "chr2".into()]);
        assert_eq!(l.chr_ids, ["chr1", "chr2"]);
        assert!(l.starts.is_empty() && l.ends.is_empty());
        // The shape `parse` refuses, which is why this constructor exists.
        assert!(Locs::parse(&l.chr_ids, &[], &[], &[], None).is_err());
    }

    #[test]
    fn whole_chromosomes_span_each_one() {
        let map = ChrMap::from_entries([("chr1".into(), 100), ("chr2".into(), 50)]);
        let l = Locs::whole_chromosomes(&map, &[]).unwrap();
        assert_eq!(l.chr_ids, ["chr1", "chr2"]);
        assert_eq!(l.starts, [0, 0]);
        assert_eq!(l.ends, [100, 50]);
        let l = Locs::whole_chromosomes(&map, &["2".into()]).unwrap();
        assert_eq!(l.chr_ids, ["chr2"]);
        assert_eq!(l.ends, [50]);
    }

    #[test]
    fn reverse_output_rows_mirrors_only_reverse_loci() {
        let locs = IndexedLocs {
            locs: vec![
                IndexedLoc {
                    chr_index: 0,
                    start: 0,
                    end: 3,
                    binned_start: 0,
                    binned_end: 3,
                    bin_size: 1.0,
                    reverse: false,
                    output_start: 0,
                    output_end: 3,
                },
                IndexedLoc {
                    chr_index: 0,
                    start: 0,
                    end: 3,
                    binned_start: 0,
                    binned_end: 3,
                    bin_size: 1.0,
                    reverse: true,
                    output_start: 3,
                    output_end: 6,
                },
            ],
            bin_count: 3,
            output_len: 6,
        };
        let mut out = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        locs.reverse_output_rows(&mut out);
        assert_eq!(out, [1.0, 2.0, 3.0, 6.0, 5.0, 4.0]);
    }
}