symcurve/curve/

iters.rs

//! This module provides data structures and Iterator implementations for the calculation of DNA
//! curvature.
//!
//! It includes the necessary data structures for representing the DNA data and the traits and
//! implementations for iterating over this data. The iterators provided allow for efficient and
//! convenient traversal and manipulation of the DNA data for the purpose of curvature calculation.
use crate::curve::matrix;
use std::collections::VecDeque;
use std::f64::consts::{PI, TAU};
use std::iter::{FusedIterator, Iterator};

/// The step a triplet contributes to the traced path.
///
/// The roll displaces along the accumulated twist and the tilt along the perpendicular:
///
/// ```text
/// dx = roll * sin(T) + tilt * sin(T - pi/2)
/// dy = roll * cos(T) + tilt * cos(T - pi/2)
/// ```
///
/// The tilt term is `T - pi/2`, matching both the reference implementation and the
/// published equations. `T + pi/2` is the opposite perpendicular, which is the same as
/// negating the tilt, so the sign here is not free to choose. It has no effect while the
/// supplied tilt matrix is uniformly zero, which is why it is pinned by a test on this
/// function rather than by one on the iterator.
fn step(roll: f64, tilt: f64, twist_sum: f64) -> TripletData {
    let perpendicular = twist_sum - PI / 2.0;
    TripletData {
        dx: (roll * twist_sum.sin()) + (tilt * perpendicular.sin()),
        dy: (roll * twist_sum.cos()) + (tilt * perpendicular.cos()),
    }
}

/// How many items a sliding window will still yield.
///
/// A window of `window` items over `remaining` inputs yields `remaining - window + 1`, and
/// `buffered` items are already held. Reported so that collecting into a `Vec` can
/// allocate once rather than growing repeatedly, which matters at genome scale.
fn window_size_hint(
    inner: (usize, Option<usize>),
    buffered: usize,
    window: usize,
) -> (usize, Option<usize>) {
    let available = |n: usize| (n + buffered + 1).saturating_sub(window);
    (available(inner.0), inner.1.map(available))
}

/// The step a triplet of nucleotides contributes to the traced path.
///
/// The twist, roll and tilt looked up for the triplet are combined into this step as it is
/// produced; only the step itself is carried forward, since nothing downstream reads the
/// individual parameters.
///
/// # Fields
///
/// * `dx`: The delta x value, calculated from the roll, tilt and accumulated twist.
/// * `dy`: The delta y value, calculated from the roll, tilt and accumulated twist.
#[derive(Clone, Copy, Debug)]
struct TripletData {
    dx: f64,
    dy: f64,
}

/// An iterator-wrapping struct that yields TripletData from an inner `u8` iterator.
///
/// `TripletWindowsIter` wraps around another iterator that yields `u8` (representing nucleotides),
/// and looks up the roll/tilt/twist values for each triplet of nucleotides in the inner iterator,
/// then populates a `TripletData` struct with these values.  While iterating, it also keeps track
/// of the sum of the twist values for the current window of triplets.
///
/// # Type Parameters
///
/// * `I`: The type of the inner iterator. Must be an iterator over `u8`.
///
/// # Fields
///
/// * `base_buffer`: A buffer that stores the current triplet of nucleotides.
/// * `inner`: The inner iterator that yields `u8`.
/// * `twist_sum`: The sum of the twist values for the current triplet.
/// * `roll_type`: The current roll type.
struct TripletWindowsIter<I> {
    base_buffer: VecDeque<u8>,
    inner: I,
    twist_sum: f64,
    roll_type: matrix::RollType,
    /// Set once the inner iterator has returned None, so it is not polled again.
    /// Iterator only guarantees anything about repeated calls after exhaustion for
    /// iterators that are Fused, and the inner one need not be.
    inner_done: bool,
}

/// Implementation of the `Iterator` trait for `TripletWindowsIter` struct.
///
/// This iterator yields `TripletData` items, which are calculated based on the next three bases
/// as a sliding window from the inner iterator, as well as the current roll type.
///
/// # Type Parameters
///
/// * `I`: The type of the inner iterator. Must be an iterator over `u8`.
///
/// # Returns
///
/// The `next` method returns `Some(TripletData)` if there are enough items left in the inner
/// iterator, or `None` if there are not.
impl<I> Iterator for TripletWindowsIter<I>
where
    I: Iterator<Item = u8>,
{
    type Item = TripletData;

    fn next(&mut self) -> Option<Self::Item> {
        // Fill the buffer with the next three items from the inner iterator.
        while !self.inner_done && self.base_buffer.len() < matrix::TRIPLET_SIZE {
            match self.inner.next() {
                Some(item) => self.base_buffer.push_back(item),
                None => self.inner_done = true,
            }
        }
        // When the buffer is full, calculate the twist, roll, and tilt values.
        if self.base_buffer.len() >= matrix::TRIPLET_SIZE {
            // Fixed-size, so no allocation: this runs once per base.
            let triplet: [u8; matrix::TRIPLET_SIZE] = [
                self.base_buffer[0],
                self.base_buffer[1],
                self.base_buffer[2],
            ];
            // Decode the ASCII bases once, then index each matrix with the result.
            let ixs = matrix::triplet_indices(&triplet).unwrap();
            let twist = matrix::lookup_by_index(&ixs, &matrix::TWIST);
            let roll = match self.roll_type {
                matrix::RollType::Simple => matrix::lookup_by_index(&ixs, &matrix::ROLL_SIMPLE),
                matrix::RollType::Active => matrix::lookup_by_index(&ixs, &matrix::ROLL_ACTIVE),
            };
            let tilt = matrix::lookup_by_index(&ixs, &matrix::TILT);
            self.twist_sum += twist;
            // Only sin and cos of this are ever used, so keeping it in [0, TAU) changes
            // nothing mathematically while stopping it from growing without bound. Left
            // to accumulate it reaches ~1.5e8 radians over a chromosome, where an ulp is
            // 3e-8 radians and the per-step rounding has compounded far past that.
            if !(0.0..TAU).contains(&self.twist_sum) {
                self.twist_sum = self.twist_sum.rem_euclid(TAU);
            }
            // Create a TripletData instance and return it.
            let window = step(roll, tilt, self.twist_sum);
            self.base_buffer.pop_front();
            Some(window)
        } else {
            None
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        window_size_hint(
            self.inner.size_hint(),
            self.base_buffer.len(),
            matrix::TRIPLET_SIZE,
        )
    }
}

impl<I> FusedIterator for TripletWindowsIter<I> where I: Iterator<Item = u8> {}

/// A trait for `u8` Iterators to yield `TripletData`.
///
/// `TripletWindowsIterator` is a trait for iterators over `u8` that provides a method for
/// transforming the iterator into a `TripletWindowsIter`. This allows for convenient conversion
/// of any iterator over `u8` into an iterator that yields triplets of nucleotides. This is
/// **layer 1** of the iterator stack.
///
/// # Type Parameters
///
/// * `Self`: The type implementing this trait. Must be an iterator over `u8`.
///
/// # Methods
///
/// * `triplet_windows_iter`: Takes a `RollType` and returns a `TripletWindowsIter` that yields
///   triplets of nucleotides from the original iterator.
trait TripletWindowsIterator: Iterator<Item = u8> + Sized {
    fn triplet_windows_iter(self, roll_type: matrix::RollType) -> TripletWindowsIter<Self> {
        TripletWindowsIter {
            base_buffer: VecDeque::new(),
            inner: self,
            twist_sum: 0.0,
            roll_type,
            inner_done: false,
        }
    }
}

impl<I: Iterator<Item = u8>> TripletWindowsIterator for I {}

/// A point on the path traced by the accumulated steps.
///
/// # Fields
///
/// * `x`: The x coordinate.
/// * `y`: The y coordinate.
#[derive(Clone, Copy, Debug)]
struct CoordsData {
    x: f64,
    y: f64,
}

impl CoordsData {
    /// Constructor for `CoordsData`.
    fn new(x: f64, y: f64) -> Self {
        CoordsData { x, y }
    }
}

/// An iterator-wrapping struct that yields `CoordsData` from another iterator.
///
/// `CoordsIter` wraps around another iterator that yields `TripletData`, and yields `CoordsData`
/// calculated from the `TripletData` and the previous coordinates and deltas. It also keeps track
/// of whether it has yielded the tail coordinates yet.
///
/// # Type Parameters
///
/// * `I`: The type of the inner iterator. Must be an iterator over `TripletData`.
///
/// # Fields
///
/// * `inner`: The inner iterator that yields `TripletData`.
/// * `head`: A boolean that indicates whether the first `CoordsData` has been yielded yet.
/// * `tail`: A boolean that indicates whether the end of the iterator has been reached,
///   at which point one more `CoordsData` is yielded with no associated `TripletData`.
/// * `prev_x_coord`: The x coordinate from the previous `CoordsData`.
/// * `prev_y_coord`: The y coordinate from the previous `CoordsData`.
/// * `prev_dx`: The delta x from the previous `TripletData`.
/// * `prev_dy`: The delta y from the previous `TripletData`.
struct CoordsIter<I: Iterator> {
    inner: I,
    head: bool,
    tail: bool,
    prev_x_coord: f64,
    prev_y_coord: f64,
    prev_dx: f64,
    prev_dy: f64,
}

impl<I> Iterator for CoordsIter<I>
where
    I: Iterator<Item = TripletData>,
{
    type Item = CoordsData;

    /// Advance the traced path by one step.
    ///
    /// The first point the path would yield is the origin before any step has been
    /// applied, which carries no information, so it is skipped. Once the inner iterator
    /// runs out one final point is emitted, applying the last step.
    fn next(&mut self) -> Option<Self::Item> {
        loop {
            match self.inner.next() {
                Some(triplet_data) => {
                    let point = self.step();
                    self.prev_dx = triplet_data.dx;
                    self.prev_dy = triplet_data.dy;
                    if self.head {
                        return Some(point);
                    }
                    // Discard the origin and go round again rather than recursing.
                    self.head = true;
                }
                None if !self.tail => {
                    self.tail = true;
                    return Some(self.step());
                }
                None => return None,
            }
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        // One point is dropped from the front and one added at the end, so the count
        // matches the inner iterator's, give or take what has already been consumed.
        let (lower, upper) = self.inner.size_hint();
        let pending = usize::from(!self.tail);
        (
            lower.saturating_add(pending).saturating_sub(1),
            upper.and_then(|u| u.checked_add(pending)),
        )
    }
}

impl<I> FusedIterator for CoordsIter<I> where I: Iterator<Item = TripletData> {}

impl<I> CoordsIter<I>
where
    I: Iterator<Item = TripletData>,
{
    /// Apply the previous step to the current position and return the point reached.
    fn step(&mut self) -> CoordsData {
        self.prev_x_coord += self.prev_dx;
        self.prev_y_coord += self.prev_dy;
        CoordsData::new(self.prev_x_coord, self.prev_y_coord)
    }
}

/// A trait for `TripletData` Iterators to yield `CoordsData`.
///
/// `CoordsIterator` is a trait for iterators over `TripletData` that provides a method for
/// transforming the iterator into a `CoordsIter`. This allows for convenient conversion
/// of any iterator over `TripletData` into an iterator that yields `CoordsData`. This
/// is **layer 2** of the iterator stack.
///
/// # Type Parameters
///
/// * `Self`: The type implementing this trait. Must be an iterator over `TripletData`.
///
/// # Methods
///
/// * `coords_iter`: Returns a `CoordsIter` that yields `CoordsData` calculated from the
///   `TripletData` yielded by the original iterator.
trait CoordsIterator: Iterator<Item = TripletData> + Sized {
    fn coords_iter(self) -> CoordsIter<Self> {
        CoordsIter {
            inner: self,
            head: false,
            tail: false,
            prev_x_coord: 0.0,
            prev_y_coord: 0.0,
            prev_dx: 0.0,
            prev_dy: 0.0,
        }
    }
}

impl<I: Iterator<Item = TripletData>> CoordsIterator for I {}

/// Represents the data for a rolling mean of the x and y coordinates.
///
/// # Fields
///
/// * `x_bar`: The weighted mean of the x coordinates.
/// * `y_bar`: The weighted mean of the y coordinates.
struct RollMeanData {
    x_bar: f64,
    y_bar: f64,
}

/// How many items may pass before the rolling sums are rebuilt from the buffer.
///
/// A running sum that is added to and subtracted from never sheds the rounding of the
/// values that have left it, so its error ratchets upward. Rebuilding costs one pass over
/// a window of about a hundred items, so amortised over this interval it is a fraction of
/// a percent of the work.
const ROLL_SUM_REBUILD_INTERVAL: usize = 1 << 16;

/// Represents the data for a rolling mean of the x and y coordinates.
///
/// The `RollMeanData` struct contains the weighted x and y means for a window of coordinates
/// that is 2 * `step_size` + 1 in length.
///
/// # Fields
///
/// * `inner`: The inner iterator that yields `CoordsData`.
/// * `buffer`: A buffer that stores the current window of coordinates.
/// * `step_size`: Half the size of the window minus one.  In other words,
///   2 * `step_size` + 1 is the size of the window.
/// * `x_roll_sum`: The sum of the x coordinates in the current window.
/// * `y_roll_sum`: The sum of the y coordinates in the current window.
/// * `since_rebuild`: Items processed since the rolling sums were last rebuilt.
struct RollMeanIter<I> {
    inner: I,
    buffer: VecDeque<CoordsData>,
    step_size: usize,
    x_roll_sum: f64,
    y_roll_sum: f64,
    since_rebuild: usize,
    /// Set once the inner iterator has returned None, so it is not polled again.
    inner_done: bool,
}

/// Implementation of the `Iterator` trait for `RollMeanIter`.
///
/// This iterator wraps another iterator of items of type `CoordsData` and computes
/// a rolling mean of the `x` and `y` values of the items.
impl<I> Iterator for RollMeanIter<I>
where
    I: Iterator<Item = CoordsData>,
{
    type Item = RollMeanData;

    /// Computes the next item of the rolling mean iterator.
    ///
    /// This method computes the rolling mean of the `x` and `y` values of the next
    /// `window_size` items from the inner iterator, where `window_size` is `step_size * 2 + 1`.
    ///
    /// The method returns `Some(RollMeanData)` if there are enough items in the inner iterator,
    /// and `None` otherwise.
    fn next(&mut self) -> Option<Self::Item> {
        // Fill the buffer with the next three items from the inner iterator.
        let window_size = self.step_size * 2 + 1;
        while !self.inner_done && self.buffer.len() < window_size {
            match self.inner.next() {
                Some(item) => {
                    self.x_roll_sum += item.x;
                    self.y_roll_sum += item.y;
                    self.buffer.push_back(item);
                }
                None => self.inner_done = true,
            }
        }
        if self.buffer.len() >= window_size {
            self.since_rebuild += 1;
            if self.since_rebuild >= ROLL_SUM_REBUILD_INTERVAL {
                self.x_roll_sum = self.buffer.iter().map(|item| item.x).sum();
                self.y_roll_sum = self.buffer.iter().map(|item| item.y).sum();
                self.since_rebuild = 0;
            }
            // get the fron/back items without removing them and adjust the roll sum
            let adj_x_roll_sum = self.x_roll_sum
                - (0.5 * self.buffer.front().unwrap().x)
                - (0.5 * self.buffer.back().unwrap().x);
            let adj_y_roll_sum = self.y_roll_sum
                - (0.5 * self.buffer.front().unwrap().y)
                - (0.5 * self.buffer.back().unwrap().y);
            let x_bar = adj_x_roll_sum / (window_size as f64 - 1_f64);
            let y_bar = adj_y_roll_sum / (window_size as f64 - 1_f64);
            let result = Some(RollMeanData { x_bar, y_bar });
            let item = self.buffer.pop_front().unwrap();
            self.x_roll_sum -= item.x;
            self.y_roll_sum -= item.y;
            result
        } else {
            None
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        window_size_hint(
            self.inner.size_hint(),
            self.buffer.len(),
            self.step_size * 2 + 1,
        )
    }
}

impl<I> FusedIterator for RollMeanIter<I> where I: Iterator<Item = CoordsData> {}

/// A trait for iterators that can compute a rolling mean of `CoordsData`.
///
/// This trait extends the `Iterator` trait, adding a `roll_mean_iter` method that
/// wraps the iterator in a `RollMeanIter`. The `RollMeanIter` computes a rolling mean
/// of the `x` and `y` values of the items from the original iterator.
trait RollMeanIterator: Iterator<Item = CoordsData> + Sized {
    /// Wraps the iterator in a `RollMeanIter`.
    ///
    /// This method takes ownership of the iterator and returns a `RollMeanIter` that
    /// computes a rolling mean of the `x` and `y` values of the items from the original iterator.
    ///
    /// # Parameters
    ///
    /// * `step_size`: half of the window size minus one. In other words, 2 * `step_size` + 1 is
    ///   the size of the window.
    ///
    /// # Returns
    ///
    /// A `RollMeanIter` that computes a rolling mean of the `x` and `y` values of the items.
    fn roll_mean_iter(self, step_size: usize) -> RollMeanIter<Self> {
        RollMeanIter {
            inner: self,
            buffer: VecDeque::new(),
            step_size,
            x_roll_sum: 0.0,
            y_roll_sum: 0.0,
            since_rebuild: 0,
            inner_done: false,
        }
    }
}

impl<I: Iterator<Item = CoordsData>> RollMeanIterator for I {}

/// An iterator that computes the Euclidean distance between pairs of items from an inner iterator.
///
/// `EucDistIter` wraps another iterator that yields `RollMeanData`. It computes the Euclidean
/// distance between each pair of items from the inner iterator.
///
/// # Fields
///
/// * `inner`: The inner iterator that yields `RollMeanData`.
///
/// * `buffer`: A buffer that stores 2 * `curve_step_size` + 1 items from the inner iterator.
///
/// * `curve_step_size`: The distance from the midpoint base in the window.  
struct EucDistIter<I> {
    inner: I,
    buffer: VecDeque<RollMeanData>,
    curve_step_size: usize,
    /// Set once the inner iterator has returned None, so it is not polled again.
    inner_done: bool,
}

impl<I> Iterator for EucDistIter<I>
where
    I: Iterator<Item = RollMeanData>,
{
    type Item = f64;

    /// Computes the next item of the Euclidean distance iterator.
    ///
    /// This method computes the Euclidean distance between each pair of consecutive items
    /// from the inner iterator. The Euclidean distance is computed as the square root of
    /// the sum of the squares of the differences of the `x_bar` and `y_bar` values of the items.
    ///
    /// The method returns `Some(f64)` if there are enough items in the inner iterator,
    /// and `None` otherwise.
    fn next(&mut self) -> Option<Self::Item> {
        // Fill the buffer with the next three items from the inner iterator.
        let window_size = self.curve_step_size * 2 + 1;
        while !self.inner_done && self.buffer.len() < window_size {
            match self.inner.next() {
                Some(item) => self.buffer.push_back(item),
                None => self.inner_done = true,
            }
        }
        if self.buffer.len() >= window_size {
            let left = self.buffer.front().unwrap();
            let right = self.buffer.back().unwrap();
            let curve =
                ((right.y_bar - left.y_bar).powi(2) + (right.x_bar - left.x_bar).powi(2)).sqrt();
            self.buffer.pop_front();
            Some(curve)
        } else {
            None
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        window_size_hint(
            self.inner.size_hint(),
            self.buffer.len(),
            self.curve_step_size * 2 + 1,
        )
    }
}

impl<I> FusedIterator for EucDistIter<I> where I: Iterator<Item = RollMeanData> {}

trait EucDistIterator: Iterator<Item = RollMeanData> + Sized {
    fn euc_dist_iter(self, curve_step_size: usize) -> EucDistIter<Self> {
        EucDistIter {
            inner: self,
            buffer: VecDeque::new(),
            curve_step_size,
            inner_done: false,
        }
    }
}

impl<I: Iterator<Item = RollMeanData>> EucDistIterator for I {}

/// The value assigned when a dyad's symmetry component comes out exactly zero.
///
/// A zero sum means every mirrored pair around the dyad was exactly equal, so the
/// reciprocal is undefined. The reference implementation substitutes 100 and carries on,
/// and callers downstream treat that as a saturated score rather than a real one.
pub const DEGENERATE_SYMMETRY: f64 = 100.0;

/// An iterator that computes symmetry of curvature around each candidate dyad.
///
/// This is the last stage: it consumes curvature values and yields one symmetry score per
/// dyad. A score is non-zero only where the curvature has a strict local minimum, which is
/// what a nucleosome dyad is expected to look like, and rises the more symmetric the
/// curvature is on either side of it.
///
/// # Fields
///
/// * `inner`: The inner iterator that yields curvature values.
/// * `buffer`: A buffer holding 2 * `win` + 1 curvature values, the dyad at its centre.
/// * `win`: The margin kept on each side of the dyad.
/// * `step`: The stride, both between dyads and between the mirrored pairs summed at each.
/// * `inner_done`: Set once the inner iterator has returned None, so it is not polled again.
pub struct SymCurveIter<I> {
    inner: I,
    buffer: VecDeque<f64>,
    win: usize,
    step: usize,
    inner_done: bool,
}

impl<I> Iterator for SymCurveIter<I>
where
    I: Iterator<Item = f64>,
{
    type Item = f64;

    /// Computes the symmetry score for the next dyad.
    ///
    /// Following the reference implementation, for a dyad \(d\):
    ///
    /// ```text
    /// sum    = SUM over m of |curv[d + m] - curv[d - m]|,  m = 0, step, 2*step, ... <= win/2
    /// slope  = (curv[d-1] - curv[d]) + (curv[d+1] - curv[d])
    /// weight = 1/slope   if curv[d] is a strict local minimum and slope >= 0.01
    ///          0         otherwise
    /// score  = weight / sum
    /// ```
    ///
    /// The mirrored sum runs out to `win / 2`, but a dyad is only considered once `win`
    /// values are available on each side. That wider margin is the reference's, and it
    /// means the first and last `win` curvature values yield no score even though only
    /// `win / 2` are read. Reproduced here so the output positions match.
    fn next(&mut self) -> Option<Self::Item> {
        let span = 2 * self.win + 1;
        while !self.inner_done && self.buffer.len() < span {
            match self.inner.next() {
                Some(value) => self.buffer.push_back(value),
                None => self.inner_done = true,
            }
        }
        if self.buffer.len() < span {
            return None;
        }

        let dyad = self.win;
        let half = self.win / 2;
        let mut sum = 0.0;
        let mut offset = 0;
        while offset <= half {
            sum += (self.buffer[dyad + offset] - self.buffer[dyad - offset]).abs();
            offset += self.step;
        }

        let current = self.buffer[dyad];
        let before = self.buffer[dyad - 1];
        let after = self.buffer[dyad + 1];
        let slope = (before - current) + (after - current);
        let weight = if current < before && current < after && slope >= 0.01 {
            1.0 / slope
        } else {
            0.0
        };

        // Compared against zero exactly, as the reference does: the substitution is for a
        // sum that is precisely zero, not one that is merely small.
        let score = if sum != 0.0 {
            weight / sum
        } else {
            DEGENERATE_SYMMETRY
        };

        // Advance by `step`: drop what the buffer holds and skip the rest at the source,
        // so a stride longer than the span still lands on the reference's next dyad.
        let held = self.buffer.len().min(self.step);
        self.buffer.drain(..held);
        for _ in held..self.step {
            if self.inner.next().is_none() {
                self.inner_done = true;
                break;
            }
        }
        Some(score)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        let (lower, upper) =
            window_size_hint(self.inner.size_hint(), self.buffer.len(), span_of(self.win));
        let strided = |n: usize| n.div_ceil(self.step);
        (strided(lower), upper.map(strided))
    }
}

impl<I> FusedIterator for SymCurveIter<I> where I: Iterator<Item = f64> {}

/// The number of curvature values a dyad needs in view: `win` on each side, plus itself.
fn span_of(win: usize) -> usize {
    2 * win + 1
}

/// A trait for curvature iterators to yield symmetry scores.
///
/// This is **layer 5** of the iterator stack, sitting on the curvature values that
/// [`CurveIter`] produces.
pub trait SymCurveIterator: Iterator<Item = f64> + Sized {
    /// Wraps the iterator in a [`SymCurveIter`].
    ///
    /// # Parameters
    ///
    /// * `win`: The margin kept on each side of a dyad. The reference uses 101.
    /// * `step`: The stride between dyads and between mirrored pairs. Values below 1 are
    ///   treated as 1, since a stride of zero would never advance.
    fn sym_curve_iter(self, win: usize, step: usize) -> SymCurveIter<Self> {
        SymCurveIter {
            inner: self,
            buffer: VecDeque::new(),
            win,
            step: step.max(1),
            inner_done: false,
        }
    }
}

impl<I: Iterator<Item = f64>> SymCurveIterator for I {}

/// An iterator that computes the curvature of a DNA sequence.
///
/// `CurveIter` wraps an iterator that yields `u8` and computes the curvature of the DNA sequence
/// represented by the nucleotides.
///
/// # Fields
///
/// * `inner`: The inner iterator that yields `u8`.
pub struct CurveIter<I: Iterator<Item = u8>> {
    inner: EucDistIter<RollMeanIter<CoordsIter<TripletWindowsIter<I>>>>,
    curve_scale: f64,
}

impl<I: Iterator<Item = u8>> Iterator for CurveIter<I> {
    type Item = f64;

    /// Computes the next item of the curvature iterator.
    fn next(&mut self) -> Option<Self::Item> {
        self.inner.next().map(|x| x * self.curve_scale)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        // Scaling is one-to-one, so the stack below reports the count unchanged.
        self.inner.size_hint()
    }
}

impl<I: Iterator<Item = u8>> FusedIterator for CurveIter<I> {}

/// Construct a `CurveIter` from an iterator that yields `u8`.
///
/// This function constructs a `CurveIter` from an iterator that yields `u8`. The `CurveIter`
/// computes the curvature of the DNA sequence represented by the nucleotides.
///
/// # Parameters
///
/// * `seq_iter`: An iterator that yields `u8`.
/// * `roll_type`: The type of roll (either simple or activated).
/// * `step_b`: Half of the window size minus one. In other words, 2 * `step_size` + 1 is
///   the size of the window.
/// * `step_c`: The distance from the midpoint base to the sides in the curve window.
impl<I: Iterator<Item = u8>> CurveIter<I> {
    /// Build a curvature iterator over a stream of bases.
    ///
    /// Bases must be A, C, G or T in either case; anything else will panic, so split a
    /// record with [`crate::fasta::split_seq_by_gaps`] first. For scoring whole records
    /// use [`crate::curve::scan`], which handles splitting, positions and threading.
    ///
    /// The first `step_b + step_c + 1` bases and the last of the same produce no value,
    /// because the windows need context on both sides.
    ///
    /// ```
    /// use symcurve::curve::iters::CurveIter;
    /// use symcurve::curve::matrix::RollType;
    ///
    /// let seq = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
    /// let curves: Vec<f64> =
    ///     CurveIter::new(seq.iter().copied(), RollType::Simple, 5, 15, 0.33335).collect();
    /// assert_eq!(curves.len(), seq.len() - 2 * (5 + 15 + 1));
    /// ```
    pub fn new(
        seq_iter: I,
        roll_type: matrix::RollType,
        step_b: usize,
        step_c: usize,
        curve_scale: f64,
    ) -> Self {
        Self {
            inner: seq_iter
                .triplet_windows_iter(roll_type)
                .coords_iter()
                .roll_mean_iter(step_b)
                .euc_dist_iter(step_c),
            curve_scale,
        }
    }
}

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

    /// Below is a table of some of the expected values for the triplet iterator over the DNA
    ///
    /// | pos|nuc|trip | ixs |  twist |  roll_s |   tilt |twist_sum| dx_simp | dy_simp |
    /// | --:| -:| --: | --: | -----: | ------: | -----: | ------: | ------: | ------: |
    /// |  0 | C | CCA | 330 | 0.5986 |  0.7000 | 0.0000 |  0.5986 |  0.3945 |  0.5783 |
    /// |  1 | C | CAA | 300 | 0.5986 |  6.2000 | 0.0000 |  1.1973 |  5.7725 |  2.2622 |
    /// |  2 | A | AAC | 003 | 0.5986 |  1.6000 | 0.0000 |  1.7959 |  1.5596 | -0.3572 |
    /// |  3 | A | ACA | 030 | 0.5986 |  5.8000 | 0.0000 |  2.3946 |  3.9408 | -4.2556 |
    /// |  4 | C | CAT | 301 | 0.5986 |  8.7000 | 0.0000 |  2.9932 |  1.2860 | -8.6044 |
    /// |  5 | A | ATT | 011 | 0.5986 |  0.0000 | 0.0000 |  3.5919 |  0.0000 |  0.0000 |
    /// |  6 | T | TTT | 111 | 0.5986 |  0.1000 | 0.0000 |  4.1905 | -0.0867 | -0.0498 |
    /// |  7 | T | TTT | 111 | 0.5986 |  0.1000 | 0.0000 |  4.7892 | -0.0997 |  0.0077 |
    /// |  8 | T | TTG | 112 | 0.5986 |  6.2000 | 0.0000 |  5.3878 | -4.8387 |  3.8765 |
    /// |  9 | T | TGA | 120 | 0.5986 | 10.0000 | 0.0000 |  5.9865 | -2.9238 |  9.5630 |
    /// | 10 | G | GAC | 203 | 0.5986 |  5.6000 | 0.0000 |  6.5851 |  1.6653 |  5.3467 |
    /// | 11 | A | ACT | 031 | 0.5986 |  2.0000 | 0.0000 |  7.1838 |  1.5674 |  1.2423 |
    /// | 12 | C | CTT | 311 | 0.5986 |  4.2000 | 0.0000 |  7.7824 |  4.1892 |  0.3003 |
    /// | 13 | T | TTT | 111 | 0.5986 |  0.1000 | 0.0000 |  8.3811 |  0.0864 | -0.0503 |
    /// | 14 | T | TTT | 111 | 0.5986 |  0.1000 | 0.0000 |  8.9797 |  0.0431 | -0.0903 |
    /// | 15 | T | TTT | 111 | 0.5986 |  0.1000 | 0.0000 |  9.5784 | -0.0153 | -0.0988 |
    /// | 16 | T | TTG | 112 | 0.5986 |  6.2000 | 0.0000 | 10.1770 | -4.2363 | -4.5270 |
    /// | 17 | T | TGG | 122 | 0.5986 |  0.7000 | 0.0000 | 10.7757 | -0.6831 | -0.1527 |
    /// | 18 | G | GGG | 222 | 0.5986 |  5.7000 | 0.0000 | 11.3743 | -5.2961 |  2.1075 |
    /// | 19 | G | GGA | 220 | 0.5986 |  6.2000 | 0.0000 | 11.9729 | -3.4670 |  5.1400 |
    /// | 20 | G | GAG | 202 | 0.5986 |  6.6000 | 0.0000 | 12.5716 |  0.0345 |  6.5999 |
    /// | 21 | A | AGG | 022 | 0.5986 |  4.7000 | 0.0000 | 13.1702 |  2.6688 |  3.8688 |
    /// | 22 | G | GGG | 222 | 0.5986 |  5.7000 | 0.0000 | 13.7689 |  5.3178 |  2.0520 |
    /// | 23 | G | GGC | 223 | 0.5986 |  8.2000 | 0.0000 | 14.3675 |  7.9834 | -1.8724 |
    /// | 24 | G | GCA | 230 | 0.5986 |  7.5000 | 0.0000 | 14.9662 |  5.0670 | -5.5295 |
    /// | 25 | C | CAC | 303 | 0.5986 |  6.8000 | 0.0000 | 15.5648 |  0.9700 | -6.7305 |
    /// | 26 | A | ACT | 031 | 0.5986 |  2.0000 | 0.0000 | 16.1635 | -0.8799 | -1.7961 |
    /// | 27 | C | CTA | 310 | 0.5986 |  7.8000 | 0.0000 | 16.7621 | -6.7820 | -3.8528 |
    /// | 28 | T | TAG | 102 | 0.5986 |  7.8000 | 0.0000 | 17.3608 | -7.7738 |  0.6390 |
    /// | 29 | A | AGC | 023 | 0.5986 |  6.3000 | 0.0000 | 17.9594 | -4.8961 |  3.9646 |
    /// | 30 | G | GCA | 230 | 0.5986 |  7.5000 | 0.0000 | 18.5581 | -2.1553 |  7.1836 |
    /// | 31 | C | CAC | 303 | 0.5986 |  6.8000 | 0.0000 | 19.1567 |  2.0560 |  6.4817 |
    /// | 32 | A | ACC | 033 | 0.5986 |  5.2000 | 0.0000 | 19.7554 |  4.0920 |  3.2087 |
    /// | 33 | C | CCT | 331 | 0.5986 |  4.7000 | 0.0000 | 20.3540 |  4.6897 |  0.3116 |
    /// | 34 | C | CTA | 310 | 0.5986 |  7.8000 | 0.0000 | 20.9527 |  6.7208 | -3.9587 |
    /// | 35 | T | TAT | 101 | 0.5986 |  9.7000 | 0.0000 | 21.5513 |  4.1302 | -8.7767 |
    /// | 36 | A | ATC | 013 | 0.5986 |  3.6000 | 0.0000 | 22.1500 | -0.5693 | -3.5547 |
    /// | 37 | T | TCT | 131 | 0.5986 |  6.5000 | 0.0000 | 22.7486 | -4.4660 | -4.7228 |
    /// | 38 | C | CTA | 310 | 0.5986 |  7.8000 | 0.0000 | 23.3472 | -7.6209 | -1.6618 |
    /// | 39 | T | TAC | 103 | 0.5986 |  6.4000 | 0.0000 | 23.9459 | -5.9340 |  2.3974 |
    /// | 40 | A | ACC | 033 | 0.5986 |  5.2000 | 0.0000 | 24.5445 | -2.8853 |  4.3261 |
    /// | 41 | C | CCC | 333 | 0.5986 |  5.7000 | 0.0000 | 25.1432 |  0.0596 |  5.6997 |
    /// | 42 | C | CCT | 331 | 0.5986 |  4.7000 | 0.0000 | 25.7418 |  2.6890 |  3.8548 |
    /// | 43 | C | CTG | 312 | 0.5986 |  9.6000 | 0.0000 | 26.3405 |  8.9743 |  3.4092 |
    /// | 44 | T | TGA | 120 | 0.5986 | 10.0000 | 0.0000 | 26.9391 |  9.7238 | -2.3342 |
    /// | 45 | G | GAA | 200 | 0.5986 |  5.1000 | 0.0000 | 27.5378 |  3.4259 | -3.7780 |
    /// | 46 | A | AAT | 001 | 0.5986 |  0.0000 | 0.0000 | 28.1364 |  0.0000 |  0.0000 |
    /// | 47 | A | ATC | 013 | 0.5986 |  3.6000 | 0.0000 | 28.7351 | -1.6006 | -3.2246 |
    /// | 48 | T |     |     |         |        |        |         |         |         |
    /// | 49 | C |     |     |         |        |        |         |         |         |
    #[test]
    fn test_triplet_iter_long() {
        let dna = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let windows: Vec<TripletData> = dna
            .iter()
            .cloned()
            .triplet_windows_iter(matrix::RollType::Simple)
            .collect();
        assert_eq!(windows.len(), dna.len() - 2);
        // check first two
        assert_relative_eq!(windows[0].dx, 0.3945, epsilon = 1e-4);
        assert_relative_eq!(windows[0].dy, 0.5783, epsilon = 1e-4);
        assert_relative_eq!(windows[1].dx, 5.7725, epsilon = 1e-4);
        assert_relative_eq!(windows[1].dy, 2.2622, epsilon = 1e-4);
        // check last two
        assert_relative_eq!(windows[46].dx, 0.0000, epsilon = 1e-4);
        assert_relative_eq!(windows[46].dy, 0.0000, epsilon = 1e-4);
        assert_relative_eq!(windows[47].dx, -1.6006, epsilon = 1e-4);
        assert_relative_eq!(windows[47].dy, -3.2246, epsilon = 1e-4);
    }

    #[test]
    fn test_triplet_iter_too_short() {
        let dna = b"AC";
        let windows: Vec<TripletData> = dna
            .iter()
            .cloned()
            .triplet_windows_iter(matrix::RollType::Simple)
            .collect();
        assert_eq!(windows.len(), 0);
    }

    /// Below is a table of some of the expected values for the coords iterator over the DNA
    ///
    /// | pos|nuc|trip | dx_simp | dy_simp |  x_coord |  y_coord |
    /// | --:| -:| --: | ------: | ------: | -------: | -------: |
    /// |  0 | C | CCA |  0.3945 |  0.5783 |          |          |
    /// |  1 | C | CAA |  5.7725 |  2.2622 |   0.3945 |   0.5783 |
    /// |  2 | A | AAC |  1.5596 | -0.3572 |   6.1670 |   2.8405 |
    /// |  3 | A | ACA |  3.9408 | -4.2556 |   7.7266 |   2.4833 |
    /// |  4 | C | CAT |  1.2860 | -8.6044 |  11.6674 |  -1.7723 |
    /// |  5 | A | ATT |  0.0000 |  0.0000 |  12.9534 | -10.3767 |
    /// |  6 | T | TTT | -0.0867 | -0.0498 |  12.9534 | -10.3767 |
    /// |  7 | T | TTT | -0.0997 |  0.0077 |  12.8667 | -10.4266 |
    /// |  8 | T | TTG | -4.8387 |  3.8765 |  12.7670 | -10.4189 |
    /// |  9 | T | TGA | -2.9238 |  9.5630 |   7.9283 |  -6.5424 |
    /// | 10 | G | GAC |  1.6653 |  5.3467 |   5.0045 |   3.0206 |
    /// | 11 | A | ACT |  1.5674 |  1.2423 |   6.6698 |   8.3673 |
    /// | 12 | C | CTT |  4.1892 |  0.3003 |   8.2372 |   9.6096 |
    /// | 13 | T | TTT |  0.0864 | -0.0503 |  12.4264 |   9.9099 |
    /// | 14 | T | TTT |  0.0431 | -0.0903 |  12.5128 |   9.8596 |
    /// | 15 | T | TTT | -0.0153 | -0.0988 |  12.5559 |   9.7693 |
    /// | 16 | T | TTG | -4.2363 | -4.5270 |  12.5406 |   9.6705 |
    /// | 17 | T | TGG | -0.6831 | -0.1527 |   8.3043 |   5.1435 |
    /// | 18 | G | GGG | -5.2961 |  2.1075 |   7.6212 |   4.9908 |
    /// | 19 | G | GGA | -3.4670 |  5.1400 |   2.3251 |   7.0983 |
    /// | 20 | G | GAG |  0.0345 |  6.5999 |  -1.1419 |  12.2383 |
    /// | 21 | A | AGG |  2.6688 |  3.8688 |  -1.1074 |  18.8382 |
    /// | 22 | G | GGG |  5.3178 |  2.0520 |   1.5614 |  22.7069 |
    /// | 23 | G | GGC |  7.9834 | -1.8724 |   6.8792 |  24.7590 |
    /// | 24 | G | GCA |  5.0670 | -5.5295 |  14.8626 |  22.8866 |
    /// | 25 | C | CAC |  0.9700 | -6.7305 |  19.9296 |  17.3571 |
    /// | 26 | A | ACT | -0.8799 | -1.7961 |  20.8995 |  10.6266 |
    /// | 27 | C | CTA | -6.7820 | -3.8528 |  20.0197 |   8.8305 |
    /// | 28 | T | TAG | -7.7738 |  0.6390 |  13.2377 |   4.9777 |
    /// | 29 | A | AGC | -4.8961 |  3.9646 |   5.4639 |   5.6167 |
    /// | 30 | G | GCA | -2.1553 |  7.1836 |   0.5678 |   9.5814 |
    /// | 31 | C | CAC |  2.0560 |  6.4817 |  -1.5875 |  16.7650 |
    /// | 32 | A | ACC |  4.0920 |  3.2087 |   0.4685 |  23.2467 |
    /// | 33 | C | CCT |  4.6897 |  0.3116 |   4.5605 |  26.4554 |
    /// | 34 | C | CTA |  6.7208 | -3.9587 |   9.2502 |  26.7669 |
    /// | 35 | T | TAT |  4.1302 | -8.7767 |  15.9709 |  22.8083 |
    /// | 36 | A | ATC | -0.5693 | -3.5547 |  20.1012 |  14.0315 |
    /// | 37 | T | TCT | -4.4660 | -4.7228 |  19.5319 |  10.4768 |
    /// | 38 | C | CTA | -7.6209 | -1.6618 |  15.0659 |   5.7540 |
    /// | 39 | T | TAC | -5.9340 |  2.3974 |   7.4450 |   4.0922 |
    /// | 40 | A | ACC | -2.8853 |  4.3261 |   1.5109 |   6.4896 |
    /// | 41 | C | CCC |  0.0596 |  5.6997 |  -1.3743 |  10.8157 |
    /// | 42 | C | CCT |  2.6890 |  3.8548 |  -1.3148 |  16.5154 |
    /// | 43 | C | CTG |  8.9743 |  3.4092 |   1.3742 |  20.3701 |
    /// | 44 | T | TGA |  9.7238 | -2.3342 |  10.3485 |  23.7794 |
    /// | 45 | G | GAA |  3.4259 | -3.7780 |  20.0722 |  21.4451 |
    /// | 46 | A | AAT |  0.0000 |  0.0000 |  23.4981 |  17.6671 |
    /// | 47 | A | ATC | -1.6006 | -3.2246 |  23.4981 |  17.6671 |
    /// | 48 | T |     |         |         |  21.8975 |  14.4425 |
    /// | 49 | C |     |         |         |          |          |
    #[test]
    fn test_coords_iter() {
        let dna = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let coords: Vec<CoordsData> = dna
            .iter()
            .cloned()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .collect();
        let coords_len = coords.len();
        assert_eq!(coords_len, dna.len() - 2);
        // check first two
        assert_relative_eq!(coords[0].x, 0.3945, epsilon = 1e-4);
        assert_relative_eq!(coords[0].y, 0.5783, epsilon = 1e-4);
        assert_relative_eq!(coords[1].x, 6.1670, epsilon = 1e-4);
        assert_relative_eq!(coords[1].y, 2.8405, epsilon = 1e-4);
        // check last two
        assert_relative_eq!(coords[coords_len - 2].x, 23.4981, epsilon = 1e-4);
        assert_relative_eq!(coords[coords_len - 2].y, 17.6671, epsilon = 1e-4);
        assert_relative_eq!(coords[coords_len - 1].x, 21.8975, epsilon = 1e-4);
        assert_relative_eq!(coords[coords_len - 1].y, 14.4425, epsilon = 1e-4);
    }

    /// Helper for test_rollmean_iter()
    fn get_some_coords() -> Vec<CoordsData> {
        let x_values = vec![
            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
        ];
        let y_values = vec![
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0,
        ];

        x_values
            .into_iter()
            .zip(y_values)
            .map(|(x, y)| CoordsData::new(x, y))
            .collect()
    }

    #[test]
    fn test_rollmean_iter() {
        let rolls: Vec<_> = get_some_coords().into_iter().roll_mean_iter(2).collect();
        assert_eq!(rolls.len(), 8);
        // x̄₃ = (½x₁ + x₂ + x₃ + x₄ + ½x₅)/4
        // x̄₃ = (0.5 + 2 + 3 + 4 + 2.5)/4 = 3
        assert_relative_eq!(rolls[0].x_bar, 3.0, epsilon = 1e-4);
        assert_relative_eq!(rolls[0].y_bar, 0.0, epsilon = 1e-4);
        // x̄₃ = (½x₂ + x₃ + x₄ + x₅ + ½x₆)/4
        // x̄₃ = (1 + 3 + 4 + 5 + 3)/4 = 16 / 4 = 4
        assert_relative_eq!(rolls[1].x_bar, 4.0, epsilon = 1e-4);
        assert_relative_eq!(rolls[2].x_bar, 5.0, epsilon = 1e-4);
        assert_relative_eq!(rolls[7].y_bar, 10.0, epsilon = 1e-4);
        let rolls: Vec<_> = get_some_coords().into_iter().roll_mean_iter(3).collect();
        // x̄₃ = (½x₁ + x₂ + x₃ + x₄ + x₅ + x₆+ ½x₇)/6
        // x̄₃ = (0.5 + 2 + 3 + 4 + 5 + 6 + 3.5)/6 = 24 / 6 = 4
        assert_relative_eq!(rolls[0].x_bar, 4.0, epsilon = 1e-4);
        assert_eq!(rolls.len(), 6);
    }

    /// | pos|nuc|trip |  x_coord |  y_coord |    x_bar |    y_bar |
    /// | --:| -:| --: | -------: | -------: | -------: | -------: |
    /// |  0 | C | CCA |          |          |          |          |
    /// |  1 | C | CAA |   0.3945 |   0.5783 |          |          |
    /// |  2 | A | AAC |   6.1670 |   2.8405 |          |          |
    /// |  3 | A | ACA |   7.7266 |   2.4833 |          |          |
    /// |  4 | C | CAT |  11.6674 |  -1.7723 |          |          |
    /// |  5 | A | ATT |  12.9534 | -10.3767 |          |          |
    /// |  6 | T | TTT |  12.9534 | -10.3767 |   9.3566 |  -3.7097 |
    /// |  7 | T | TTT |  12.8667 | -10.4266 |   9.7739 |  -2.9818 |
    /// |  8 | T | TTG |  12.7670 | -10.4189 |  10.1124 |  -2.2720 |
    /// |  9 | T | TGA |   7.9283 |  -6.5424 |  10.3897 |  -1.3191 |
    /// | 10 | G | GAC |   5.0045 |   3.0206 |  10.4121 |   0.2698 |
    /// | 11 | A | ACT |   6.6698 |   8.3673 |  10.3716 |   2.2795 |
    /// | 12 | C | CTT |   8.2372 |   9.6096 |  10.1228 |   4.0604 |
    /// | 13 | T | TTT |  12.4264 |   9.9099 |   9.6374 |   5.6094 |
    /// | 14 | T | TTT |  12.5128 |   9.8596 |   9.0999 |   7.0619 |
    /// | 15 | T | TTT |  12.5559 |   9.7693 |   8.5125 |   8.2048 |
    /// | 16 | T | TTG |  12.5406 |   9.6705 |   7.8163 |   9.1892 |
    /// | 17 | T | TGG |   8.3043 |   5.1435 |   7.0936 |  10.3676 |
    /// | 18 | G | GGG |   7.6212 |   4.9908 |   6.4825 |  11.7650 |
    /// | 19 | G | GGA |   2.3251 |   7.0983 |   6.3226 |  13.1588 |
    /// | 20 | G | GAG |  -1.1419 |  12.2383 |   6.8088 |  14.1895 |
    /// | 21 | A | AGG |  -1.1074 |  18.8382 |   7.5954 |  14.6167 |
    /// | 22 | G | GGG |   1.5614 |  22.7069 |   8.5991 |  14.8489 |
    /// | 23 | G | GGC |   6.8792 |  24.7590 |   9.4657 |  15.0326 |
    /// | 24 | G | GCA |  14.8626 |  22.8866 |   9.9035 |  14.9578 |
    /// | 25 | C | CAC |  19.9296 |  17.3571 |  10.1459 |  14.7509 |
    /// | 26 | A | ACT |  20.8995 |  10.6266 |  10.2074 |  14.5144 |
    /// | 27 | C | CTA |  20.0197 |   8.8305 |  10.1287 |  14.4377 |
    /// | 28 | T | TAG |  13.2377 |   4.9777 |   9.9582 |  14.5496 |
    /// | 29 | A | AGC |   5.4639 |   5.6167 |   9.5616 |  14.8284 |
    /// | 30 | G | GCA |   0.5678 |   9.5814 |   9.0830 |  15.2950 |
    /// | 31 | C | CAC |  -1.5875 |  16.7650 |   8.8452 |  15.7378 |
    /// | 32 | A | ACC |   0.4685 |  23.2467 |   8.7809 |  15.9903 |
    /// | 33 | C | CCT |   4.5605 |  26.4554 |   8.8479 |  16.1115 |
    /// | 34 | C | CTA |   9.2502 |  26.7669 |   9.0384 |  16.0740 |
    /// | 35 | T | TAT |  15.9709 |  22.8083 |   9.1846 |  15.8432 |
    /// | 36 | A | ATC |  20.1012 |  14.0315 |   9.2424 |  15.3912 |
    /// | 37 | T | TCT |  19.5319 |  10.4768 |   9.1639 |  14.7571 |
    /// | 38 | C | CTA |  15.0659 |   5.7540 |   8.9154 |  14.1163 |
    /// | 39 | T | TAC |   7.4450 |   4.0922 |   8.8110 |  13.6627 |
    /// | 40 | A | ACC |   1.5109 |   6.4896 |   9.0710 |  13.4451 |
    /// | 41 | C | CCC |  -1.3743 |  10.8157 |   9.4459 |  13.5588 |
    /// | 42 | C | CCT |  -1.3148 |  16.5154 |   9.8141 |  14.1000 |
    /// | 43 | C | CTG |   1.3742 |  20.3701 |  10.3540 |  14.8940 |
    /// | 44 | T | TGA |  10.3485 |  23.7794 |          |          |
    /// | 45 | G | GAA |  20.0722 |  21.4451 |          |          |
    /// | 46 | A | AAT |  23.4981 |  17.6671 |          |          |
    /// | 47 | A | ATC |  23.4981 |  17.6671 |          |          |
    /// | 48 | T |     |  21.8975 |  14.4425 |          |          |
    /// | 49 | C |     |          |          |          |          |
    #[test]
    fn test_rollmeans_from_seq() {
        let dna = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let step_size = 5;
        let means: Vec<RollMeanData> = dna
            .iter()
            .cloned()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .roll_mean_iter(step_size)
            .collect();
        let means_len = means.len();
        assert_eq!(means_len, dna.len() - 2 - 2 * step_size);
        // check first two
        assert_relative_eq!(means[0].x_bar, 9.3566, epsilon = 1e-4);
        assert_relative_eq!(means[0].y_bar, -3.7097, epsilon = 1e-4);
        assert_relative_eq!(means[1].x_bar, 9.7739, epsilon = 1e-4);
        assert_relative_eq!(means[1].y_bar, -2.9818, epsilon = 1e-4);
        // check last two
        assert_relative_eq!(means[means_len - 2].x_bar, 9.8141, epsilon = 1e-4);
        assert_relative_eq!(means[means_len - 2].y_bar, 14.1000, epsilon = 1e-4);
        assert_relative_eq!(means[means_len - 1].x_bar, 10.3540, epsilon = 1e-4);
        assert_relative_eq!(means[means_len - 1].y_bar, 14.8940, epsilon = 1e-4);
    }

    /// Helper for test_eucdist_iter()
    fn get_some_means() -> Vec<RollMeanData> {
        let x_values = vec![3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.0, 5.0, 17.0];
        let y_values = vec![0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0];

        x_values
            .into_iter()
            .zip(y_values)
            .map(|(x_bar, y_bar)| RollMeanData { x_bar, y_bar })
            .collect()
    }

    #[test]
    fn test_eucdist_iter() {
        let mean_rolls: Vec<_> = get_some_means();
        let vec_size = mean_rolls.len();
        let curve_step_size = 2;
        let euc_dists: Vec<_> = mean_rolls
            .into_iter()
            .euc_dist_iter(curve_step_size)
            .collect();
        // check curve_step_size number of items on both flanks are discarded
        assert_eq!(euc_dists.len(), vec_size - 2 * (curve_step_size));
        // √((7.0-3.0)² + (10.0-0.0)²) = √116 = 10.770329614269007
        assert_relative_eq!(euc_dists[0], 10.7703, epsilon = 1e-4);
        // √((8.0-4.0)² + (10.0-0.0)²) = √116 = 10.770329614269007
        assert_relative_eq!(euc_dists[1], 10.7703, epsilon = 1e-4);
        // √((8.0-5.0)² + (10.0-0.0)²) = √109 = 10.44031
        assert_relative_eq!(euc_dists[2], 10.44031, epsilon = 1e-4);
        // √((5.0-6.0)² + (10.0-0.0)²) = √101 = 10.04988
        assert_relative_eq!(euc_dists[3], 10.04988, epsilon = 1e-4);
        // √((17.0-7.0)² + (10.0-10.0)²) = √100 = 10.0
        assert_relative_eq!(euc_dists[4], 10.0, epsilon = 1e-4);
    }
    /// | pos|nuc|trip |    x_bar |    y_bar |   curve |
    /// | --:| -:| --: | -------: | -------: | ------: |
    /// |  0 | C | CCA |          |          |         |
    /// |  1 | C | CAA |          |          |         |
    /// |  2 | A | AAC |          |          |         |
    /// |  3 | A | ACA |          |          |         |
    /// |  4 | C | CAT |          |          |         |
    /// |  5 | A | ATT |          |          |         |
    /// |  6 | T | TTT |   9.3566 |  -3.7097 |         |
    /// |  7 | T | TTT |   9.7739 |  -2.9818 |         |
    /// |  8 | T | TTG |  10.1124 |  -2.2720 |         |
    /// |  9 | T | TGA |  10.3897 |  -1.3191 |         |
    /// | 10 | G | GAC |  10.4121 |   0.2698 |         |
    /// | 11 | A | ACT |  10.3716 |   2.2795 |         |
    /// | 12 | C | CTT |  10.1228 |   4.0604 |         |
    /// | 13 | T | TTT |   9.6374 |   5.6094 |         |
    /// | 14 | T | TTT |   9.0999 |   7.0619 |         |
    /// | 15 | T | TTT |   8.5125 |   8.2048 |         |
    /// | 16 | T | TTG |   7.8163 |   9.1892 |         |
    /// | 17 | T | TGG |   7.0936 |  10.3676 |         |
    /// | 18 | G | GGG |   6.4825 |  11.7650 |         |
    /// | 19 | G | GGA |   6.3226 |  13.1588 |         |
    /// | 20 | G | GAG |   6.8088 |  14.1895 |         |
    /// | 21 | A | AGG |   7.5954 |  14.6167 | 19.1012 |
    /// | 22 | G | GGG |   8.5991 |  14.8489 | 17.7494 |
    /// | 23 | G | GGC |   9.4657 |  15.0326 | 16.4319 |
    /// | 24 | G | GCA |   9.9035 |  14.9578 | 15.0647 |
    /// | 25 | C | CAC |  10.1459 |  14.7509 | 13.2434 |
    /// | 26 | A | ACT |  10.2074 |  14.5144 | 11.3172 |
    /// | 27 | C | CTA |  10.1287 |  14.4377 | 10.0444 |
    /// | 28 | T | TAG |   9.9582 |  14.5496 |  9.3122 |
    /// | 29 | A | AGC |   9.5616 |  14.8284 |         |
    /// | 30 | G | GCA |   9.0830 |  15.2950 |         |
    /// | 31 | C | CAC |   8.8452 |  15.7378 |         |
    /// | 32 | A | ACC |   8.7809 |  15.9903 |         |
    /// | 33 | C | CCT |   8.8479 |  16.1115 |         |
    /// | 34 | C | CTA |   9.0384 |  16.0740 |         |
    /// | 35 | T | TAT |   9.1846 |  15.8432 |         |
    /// | 36 | A | ATC |   9.2424 |  15.3912 |         |
    /// | 37 | T | TCT |   9.1639 |  14.7571 |         |
    /// | 38 | C | CTA |   8.9154 |  14.1163 |         |
    /// | 39 | T | TAC |   8.8110 |  13.6627 |         |
    /// | 40 | A | ACC |   9.0710 |  13.4451 |         |
    /// | 41 | C | CCC |   9.4459 |  13.5588 |         |
    /// | 42 | C | CCT |   9.8141 |  14.1000 |         |
    /// | 43 | C | CTG |  10.3540 |  14.8940 |         |
    /// | 44 | T | TGA |          |          |         |
    /// | 45 | G | GAA |          |          |         |
    /// | 46 | A | AAT |          |          |         |
    /// | 47 | A | ATC |          |          |         |
    /// | 48 | T |     |          |          |         |
    /// | 49 | C |     |          |          |         |
    #[test]
    fn test_eucdist_iter_from_seq() {
        let dna = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let step_size = 5;
        let curve_step = 15;
        let curves: Vec<_> = dna
            .iter()
            .cloned()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .roll_mean_iter(step_size)
            .euc_dist_iter(curve_step)
            .collect();
        let curves_len = curves.len();
        assert_eq!(
            curves_len,
            dna.len() - 2 - (2 * step_size) - (2 * curve_step)
        );
        // check all
        assert_relative_eq!(curves[0], 19.1012, epsilon = 1e-4);
        assert_relative_eq!(curves[1], 17.7494, epsilon = 1e-4);
        assert_relative_eq!(curves[2], 16.4319, epsilon = 1e-4);
        assert_relative_eq!(curves[3], 15.0647, epsilon = 1e-4);
        assert_relative_eq!(curves[4], 13.2434, epsilon = 1e-4);
        assert_relative_eq!(curves[5], 11.3172, epsilon = 1e-4);
        assert_relative_eq!(curves[6], 10.0444, epsilon = 1e-4);
        assert_relative_eq!(curves[7], 9.3122, epsilon = 1e-4);
    }

    #[test]
    fn test_curve_iter() {
        let seq = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let seq_len = seq.len();
        let curves: Vec<_> = CurveIter::new(
            seq.iter().cloned(),
            matrix::RollType::Simple,
            5,
            15,
            0.33335,
        )
        .collect();
        assert_eq!(curves.len(), seq_len - (21 * 2));
        assert_relative_eq!(curves[0], 6.3674, epsilon = 1e-4);
        assert_relative_eq!(curves[1], 5.9168, epsilon = 1e-4);
        assert_relative_eq!(curves[2], 5.4776, epsilon = 1e-4);
        assert_relative_eq!(curves[3], 5.0218, epsilon = 1e-4);
        assert_relative_eq!(curves[4], 4.4147, epsilon = 1e-4);
        assert_relative_eq!(curves[5], 3.7726, epsilon = 1e-4);
        assert_relative_eq!(curves[6], 3.3483, epsilon = 1e-4);
        assert_relative_eq!(curves[7], 3.1042, epsilon = 1e-4);
    }

    /// A direct transcription of the reference implementation's SYMCURV subroutine,
    /// kept deliberately unidiomatic so it reads against the Perl line by line and can
    /// serve as an oracle for the iterator.
    ///
    /// ```perl
    /// for (my $dyad = $win ; $dyad < scalar(@Curv) - $win ; $dyad += $step) {
    ///     my $weight = 0; my $sum = 0;
    ///     for (my $j = $dyad, my $k = $dyad ;
    ///          $j < $dyad + int($win/2) + 1, $k > $dyad - int($win/2) - 1 ;
    ///          $j += $step, $k -= $step) {
    ///         $sum += abs($Curv[$j] - $Curv[$k]);
    ///     }
    ///     if (($Curv[$dyad] < $Curv[$dyad-1]) and ($Curv[$dyad] < $Curv[$dyad+1])
    ///         and ((($Curv[$dyad-1]-$Curv[$dyad]) + ($Curv[$dyad+1]-$Curv[$dyad])) >= 0.01)) {
    ///         $weight = 1/(($Curv[$dyad-1]-$Curv[$dyad]) + ($Curv[$dyad+1]-$Curv[$dyad]))
    ///     } else { $weight = 0; }
    ///     if ($sum != 0) { $symcurv[$dyad] = (1/$sum) * $weight; }
    ///     else           { $symcurv[$dyad] = 100; }
    /// }
    /// ```
    ///
    /// Note the inner loop's comma operator: Perl evaluates only the last condition, so
    /// the `$j` bound is dead and `$k` alone terminates the loop. Transcribed as written.
    fn perl_symcurv(curv: &[f64], win: usize, step: usize) -> Vec<(usize, f64)> {
        let mut out = Vec::new();
        if curv.len() < 2 * win + 1 {
            return out;
        }
        let half = win / 2;
        let mut dyad = win;
        while dyad < curv.len() - win {
            let mut sum = 0.0;
            let (mut j, mut k) = (dyad, dyad);
            // `$k > $dyad - int($win/2) - 1`
            while k + half + 1 > dyad {
                sum += (curv[j] - curv[k]).abs();
                j += step;
                if k < step {
                    break;
                }
                k -= step;
            }
            let weight = if curv[dyad] < curv[dyad - 1]
                && curv[dyad] < curv[dyad + 1]
                && ((curv[dyad - 1] - curv[dyad]) + (curv[dyad + 1] - curv[dyad])) >= 0.01
            {
                1.0 / ((curv[dyad - 1] - curv[dyad]) + (curv[dyad + 1] - curv[dyad]))
            } else {
                0.0
            };
            let value = if sum != 0.0 {
                (1.0 / sum) * weight
            } else {
                100.0
            };
            out.push((dyad, value));
            dyad += step;
        }
        out
    }

    fn synthetic_curvature(n: usize, seed: u64) -> Vec<f64> {
        // Smooth-ish with genuine local minima, so the minimum test is actually exercised.
        let mut x = seed;
        (0..n)
            .map(|i| {
                x ^= x << 13;
                x ^= x >> 7;
                x ^= x << 17;
                let noise = (x >> 11) as f64 / (1u64 << 53) as f64;
                2.0 + (i as f64 / 7.0).sin() + 0.5 * (i as f64 / 3.0).cos() + 0.05 * noise
            })
            .collect()
    }

    #[test]
    fn test_sym_curve_matches_the_reference_transcription() {
        for (n, win, step) in [
            (600usize, 101usize, 1usize),
            (600, 101, 3),
            (400, 51, 1),
            (300, 20, 1),
            (300, 20, 7),
            (250, 101, 1), // too short: no dyads at all
        ] {
            let curv = synthetic_curvature(n, 0x2545F4914F6CDD1D ^ n as u64);
            let expected = perl_symcurv(&curv, win, step);
            let got: Vec<f64> = curv.iter().copied().sym_curve_iter(win, step).collect();
            assert_eq!(
                got.len(),
                expected.len(),
                "count differs for n={n} win={win} step={step}"
            );
            for (i, (&value, &(dyad, want))) in got.iter().zip(&expected).enumerate() {
                assert_relative_eq!(value, want, epsilon = 1e-12, max_relative = 1e-12);
                // The first score belongs to curvature index `win`, then every `step`.
                assert_eq!(dyad, win + i * step);
            }
        }
    }

    #[test]
    fn test_sym_curve_is_zero_away_from_local_minima() {
        // A strictly increasing curve has no local minimum, so every dyad scores zero.
        let curv: Vec<f64> = (0..500).map(|i| i as f64 * 0.01).collect();
        let got: Vec<f64> = curv.iter().copied().sym_curve_iter(101, 1).collect();
        assert!(!got.is_empty());
        assert!(
            got.iter().all(|&v| v == 0.0),
            "expected all zero, got {got:?}"
        );
    }

    #[test]
    fn test_sym_curve_saturates_when_perfectly_symmetric() {
        // Constant curvature makes every mirrored pair equal, so the sum is exactly zero
        // and the reference substitutes 100.
        let curv = vec![1.25f64; 500];
        let got: Vec<f64> = curv.iter().copied().sym_curve_iter(101, 1).collect();
        assert!(!got.is_empty());
        assert!(
            got.iter().all(|&v| v == DEGENERATE_SYMMETRY),
            "expected the saturated value"
        );
    }

    #[test]
    fn test_sym_curve_scores_a_symmetric_minimum() {
        // A V shape centred in the window: a genuine local minimum with perfectly
        // symmetric sides, so weight is positive and the score is finite and positive.
        let win = 20usize;
        let n = 2 * win + 1;
        let centre = win as f64;
        let curv: Vec<f64> = (0..n)
            .map(|i| 1.0 + (i as f64 - centre).abs() * 0.1)
            .collect();
        let got: Vec<f64> = curv.iter().copied().sym_curve_iter(win, 1).collect();
        assert_eq!(got.len(), 1);
        // Mirrored pairs are equal by construction, so the sum is zero and it saturates.
        assert_eq!(got[0], DEGENERATE_SYMMETRY);

        // Break the symmetry slightly: now the sum is non-zero and the score is finite.
        let mut skewed = curv.clone();
        skewed[win + 3] += 0.4;
        let got: Vec<f64> = skewed.iter().copied().sym_curve_iter(win, 1).collect();
        assert_eq!(got.len(), 1);
        assert!(got[0] > 0.0 && got[0].is_finite(), "got {}", got[0]);
        assert!(got[0] < DEGENERATE_SYMMETRY);
    }

    #[test]
    fn test_sym_curve_needs_win_on_both_sides() {
        for len in [0usize, 1, 100, 202, 203, 204] {
            let curv = synthetic_curvature(len, 7);
            let got: Vec<f64> = curv.iter().copied().sym_curve_iter(101, 1).collect();
            let expected = len.saturating_sub(2 * 101);
            assert_eq!(got.len(), expected, "len {len}");
        }
    }

    #[test]
    fn test_sym_curve_stacks_onto_the_curve_iterator() {
        // The whole pipeline, bases through to symmetry.
        let unit = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let mut seq = Vec::new();
        for _ in 0..20 {
            seq.extend_from_slice(unit);
        }
        let (step_b, step_c, win) = (5usize, 15usize, 101usize);
        let curves: Vec<f64> = CurveIter::new(
            seq.iter().copied(),
            matrix::RollType::Simple,
            step_b,
            step_c,
            0.33335,
        )
        .collect();
        let direct: Vec<f64> = curves.iter().copied().sym_curve_iter(win, 1).collect();
        let stacked: Vec<f64> = CurveIter::new(
            seq.iter().copied(),
            matrix::RollType::Simple,
            step_b,
            step_c,
            0.33335,
        )
        .sym_curve_iter(win, 1)
        .collect();
        assert_eq!(direct.len(), curves.len() - 2 * win);
        assert_eq!(direct, stacked);
        assert!(direct.iter().any(|&v| v > 0.0), "no dyad scored at all");
    }

    #[test]
    fn test_step_places_tilt_on_the_reference_side() {
        // The supplied tilt matrix is uniformly zero, so no test driving the iterator can
        // tell T - pi/2 from T + pi/2. Exercise the formula directly with a non-zero tilt.
        //
        // Expected values come from the reference implementation's expression,
        //   dx = roll*sin(T) + tilt*sin(T - pi/2)
        //   dy = roll*cos(T) + tilt*cos(T - pi/2)
        // written out here as its trigonometric identity so the test does not simply
        // restate the code: sin(T - pi/2) = -cos(T) and cos(T - pi/2) = sin(T).
        for &(roll, tilt, twist) in &[
            (5.0, 2.0, 0.0),
            (0.7, 1.5, 0.598647428),
            (3.05865, -0.25, 1.7),
            (0.0, 1.0, 3.0),
            (6.2, 0.0, 2.5),
        ] {
            let got = step(roll, tilt, twist);
            let expected_dx = roll * twist.sin() - tilt * twist.cos();
            let expected_dy = roll * twist.cos() + tilt * twist.sin();
            assert_relative_eq!(got.dx, expected_dx, epsilon = 1e-12);
            assert_relative_eq!(got.dy, expected_dy, epsilon = 1e-12);
        }
    }

    #[test]
    fn test_step_tilt_sign_is_not_the_opposite_perpendicular() {
        // T + pi/2 is the exact negation of T - pi/2, so a sign slip is invisible unless
        // something checks it. With a non-zero tilt the two differ by twice the tilt term.
        let (roll, tilt, twist) = (2.0, 1.0, 0.9);
        let got = step(roll, tilt, twist);
        let wrong_dx = roll * twist.sin() + tilt * (twist + PI / 2.0).sin();
        assert!(
            (got.dx - wrong_dx).abs() > 1.0,
            "dx {} matches the opposite perpendicular {}",
            got.dx,
            wrong_dx
        );
        // Tilt contributes nothing when it is zero, whichever convention is used.
        assert_relative_eq!(
            step(roll, 0.0, twist).dx,
            roll * twist.sin(),
            epsilon = 1e-12
        );
    }

    #[test]
    fn test_step_matches_the_pipeline() {
        // The iterator must actually be using this function.
        let seq = b"CCAACATTTT";
        let twist = matrix::matrix_lookup(b"CCA", &matrix::TWIST).unwrap();
        let roll = matrix::matrix_lookup(b"CCA", &matrix::ROLL_SIMPLE).unwrap();
        let tilt = matrix::matrix_lookup(b"CCA", &matrix::TILT).unwrap();
        let first = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple)
            .next()
            .unwrap();
        let expected = step(roll, tilt, twist);
        assert_relative_eq!(first.dx, expected.dx, epsilon = 1e-12);
        assert_relative_eq!(first.dy, expected.dy, epsilon = 1e-12);
    }

    #[test]
    fn test_roll_type_selects_the_matching_matrix() {
        // Guards the pairing that the ROLL_SIMPLE / ROLL_ACTIVE doc comments describe.
        // The two matrices disagree at CCA (0.7 simple, 3.05865 active), so this fails if
        // the constants are ever swapped to "fix" a mismatch with their docs.
        //
        // Checked through dx rather than a stored roll value: after one triplet the
        // accumulated twist is a single TWIST entry, so dx is roll * sin(twist), which
        // pins the routing and the step formula together.
        let cca = b"CCA";
        let twist = matrix::matrix_lookup(cca, &matrix::TWIST).unwrap();
        let dx_for = |roll_type: matrix::RollType| -> f64 {
            cca.iter()
                .copied()
                .triplet_windows_iter(roll_type)
                .next()
                .unwrap()
                .dx
        };
        assert_relative_eq!(
            dx_for(matrix::RollType::Simple),
            0.7 * twist.sin(),
            epsilon = 1e-12
        );
        assert_relative_eq!(
            dx_for(matrix::RollType::Active),
            3.05865 * twist.sin(),
            epsilon = 1e-12
        );
        assert_relative_eq!(
            matrix::matrix_lookup(cca, &matrix::ROLL_SIMPLE).unwrap(),
            0.7,
            epsilon = 1e-9
        );
        assert_relative_eq!(
            matrix::matrix_lookup(cca, &matrix::ROLL_ACTIVE).unwrap(),
            3.05865,
            epsilon = 1e-9
        );
    }

    #[test]
    fn test_long_runs_agree_with_locally_computed_scores() {
        // A score depends only on the bases within lead_in of it, so computing one at the
        // far end of a long sequence must match computing it from a short slice around
        // that position. Accumulated state is what breaks this: before twist was kept
        // bounded and the rolling sums rebuilt, the same comparison drifted to 1.6e-10
        // over this input, so the threshold here fails against that behaviour.
        let bases = *b"ACGT";
        let mut x: u64 = 0x5555AAAA33337777;
        let n = 300_000usize;
        let (step_b, step_c) = (5usize, 15usize);
        let lead = step_b + step_c + 1;
        let seq: Vec<u8> = (0..n)
            .map(|_| {
                x ^= x << 13;
                x ^= x >> 7;
                x ^= x << 17;
                bases[(x % 4) as usize]
            })
            .collect();

        let score = |s: &[u8]| -> Vec<f64> {
            CurveIter::new(
                s.iter().copied(),
                matrix::RollType::Simple,
                step_b,
                step_c,
                0.33335,
            )
            .collect()
        };

        let long = score(&seq);
        for &p in &[long.len() - 1, long.len() / 2, long.len() - 1000] {
            let local = score(&seq[p..p + 2 * lead + 1]);
            let rel = (long[p] - local[0]).abs() / long[p].abs().max(1e-12);
            assert!(
                rel < 1e-11,
                "score {p} drifted: long {} vs local {} (rel {rel:.3e})",
                long[p],
                local[0]
            );
        }
    }

    #[test]
    fn test_rolling_sums_stay_equal_to_a_fresh_sum() {
        // Run past the rebuild interval so the rebuild path is exercised, and check the
        // rolling means still match ones computed directly over each window.
        let bases = *b"ACGT";
        let mut x: u64 = 0x0F1E2D3C4B5A6978;
        let n = ROLL_SUM_REBUILD_INTERVAL + 5_000;
        let seq: Vec<u8> = (0..n)
            .map(|_| {
                x ^= x << 13;
                x ^= x >> 7;
                x ^= x << 17;
                bases[(x % 4) as usize]
            })
            .collect();

        let step_size = 5usize;
        let window = 2 * step_size + 1;
        let coords: Vec<CoordsData> = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .collect();
        let means: Vec<RollMeanData> = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .roll_mean_iter(step_size)
            .collect();

        assert!(means.len() > ROLL_SUM_REBUILD_INTERVAL);
        for &i in &[
            0usize,
            1,
            ROLL_SUM_REBUILD_INTERVAL - 1,
            ROLL_SUM_REBUILD_INTERVAL + 1,
            means.len() - 1,
        ] {
            // The trapezoidal mean: interior at full weight, the two ends at half.
            let slice = &coords[i..i + window];
            let x_sum: f64 = slice.iter().map(|c| c.x).sum::<f64>()
                - 0.5 * slice[0].x
                - 0.5 * slice[window - 1].x;
            let expected = x_sum / (window as f64 - 1.0);
            assert_relative_eq!(
                means[i].x_bar,
                expected,
                epsilon = 1e-9,
                max_relative = 1e-12
            );
        }
    }

    #[test]
    fn test_size_hint_is_accurate_at_every_layer() {
        // A size_hint that lies is worse than none, since callers preallocate from it.
        // Check each layer against the count it actually produces, before and partway
        // through iteration.
        let seq = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let (step_b, step_c) = (5usize, 15usize);

        let triplets = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple);
        let (lo, hi) = triplets.size_hint();
        let actual = triplets.count();
        assert_eq!(actual, seq.len() - 2);
        assert!(
            lo <= actual && hi == Some(actual),
            "triplets: {lo}..{hi:?} vs {actual}"
        );

        let coords = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter();
        let (lo, hi) = coords.size_hint();
        let actual = coords.count();
        assert!(
            lo <= actual && hi.is_some_and(|h| h >= actual),
            "coords: {lo}..{hi:?} vs {actual}"
        );

        let means = seq
            .iter()
            .copied()
            .triplet_windows_iter(matrix::RollType::Simple)
            .coords_iter()
            .roll_mean_iter(step_b);
        let (lo, hi) = means.size_hint();
        let actual = means.count();
        assert!(
            lo <= actual && hi.is_some_and(|h| h >= actual),
            "means: {lo}..{hi:?} vs {actual}"
        );

        let curves = CurveIter::new(
            seq.iter().copied(),
            matrix::RollType::Simple,
            step_b,
            step_c,
            0.33335,
        );
        let (lo, hi) = curves.size_hint();
        let collected: Vec<f64> = curves.collect();
        assert_eq!(collected.len(), seq.len() - 2 * (step_b + step_c + 1));
        assert!(
            lo <= collected.len() && hi.is_some_and(|h| h >= collected.len()),
            "curves: {lo}..{hi:?} vs {}",
            collected.len()
        );

        // Partway through, the hint must still bound what is left.
        let mut it = CurveIter::new(
            seq.iter().copied(),
            matrix::RollType::Simple,
            step_b,
            step_c,
            0.33335,
        );
        it.next();
        it.next();
        let (lo, hi) = it.size_hint();
        let left = it.count();
        assert!(
            lo <= left && hi.is_some_and(|h| h >= left),
            "partway: {lo}..{hi:?} vs {left}"
        );
    }

    #[test]
    fn test_iterators_keep_returning_none_after_exhaustion() {
        // FusedIterator is a promise, and the layers previously kept polling an inner
        // iterator that had already finished, which Iterator does not define for a
        // non-fused source.
        let seq = b"CCAACATTTTGACTTTTTGGGAGGG";
        let mut curves =
            CurveIter::new(seq.iter().copied(), matrix::RollType::Simple, 2, 2, 0.33335);
        while curves.next().is_some() {}
        for _ in 0..5 {
            assert!(curves.next().is_none(), "yielded again after finishing");
        }

        // Too short to produce anything at all: still None, repeatedly.
        let mut empty = CurveIter::new(
            b"ACG".iter().copied(),
            matrix::RollType::Simple,
            5,
            15,
            0.33335,
        );
        for _ in 0..5 {
            assert!(empty.next().is_none());
        }
    }

    #[test]
    fn test_curve_iter_is_case_insensitive() {
        // Soft-masked sequence must score identically to the same sequence unmasked.
        // Before non-ACGT handling was added, the lowercase run panicked in the
        // matrix lookup rather than producing a value at all.
        let upper = b"CCAACATTTTGACTTTTTGGGAGGGCACTAGCACCTATCTACCCTGAATC";
        let mixed = b"CCAACATTTTgacttttTGGGAGGGCACTagcacctatcTACCCTGAATC";
        assert_eq!(upper.len(), mixed.len());

        let curve = |seq: &[u8]| -> Vec<f64> {
            CurveIter::new(
                seq.iter().copied(),
                matrix::RollType::Simple,
                5,
                15,
                0.33335,
            )
            .collect()
        };

        let from_upper = curve(upper);
        let from_mixed = curve(mixed);
        assert_eq!(from_upper.len(), from_mixed.len());
        assert!(!from_upper.is_empty());
        for (u, m) in from_upper.iter().zip(&from_mixed) {
            assert_relative_eq!(u, m, epsilon = 1e-12);
        }
    }

    #[test]
    fn test_sym_curve_step_larger_than_the_span_matches_the_reference() {
        // A stride longer than the buffered span drains the buffer entirely, so the
        // remainder of the stride has to be skipped in the source as well, or the dyads
        // drift away from the reference's `win + i * step`.
        let curv = synthetic_curvature(80, 0x9E3779B97F4A7C15);
        for (win, step) in [(2usize, 9usize), (3, 7), (5, 30)] {
            let expected = perl_symcurv(&curv, win, step);
            let got: Vec<f64> = curv.iter().copied().sym_curve_iter(win, step).collect();
            assert_eq!(got.len(), expected.len(), "count for win={win} step={step}");
            for (i, (&value, &(dyad, want))) in got.iter().zip(&expected).enumerate() {
                assert_eq!(dyad, win + i * step);
                assert_relative_eq!(value, want, epsilon = 1e-12, max_relative = 1e-12);
            }
        }
    }
}