rust_lapper/

lib.rs

//! This module provides a simple data structure for fast interval searches.
//! ## Features
//! - Extremely fast overlap queries on both ordinary genomic datasets and
//!   datasets with long intervals that engulf many shorter intervals.
//! - Extremely fast in order queries through the cursor-based `seek()` method.
//! - Extremely fast intersection counts based on the
//!   [BITS](https://arxiv.org/pdf/1208.3407.pdf) algorithm
//! - NEON acceleration on AArch64, runtime-detected AVX2 on x86-64, and an exact
//!   scalar fallback everywhere else.
//! - Parallel friendly. Queries are on an immutable structure, even for `seek()`.
//! - Consumer / Adapter paradigm. Iterators are returned and serve as the main
//!   API for interacting with the Lapper.
//!
//! ## Details:
//!
//! ```text
//!          0  1  2  3  4  5  6  7  8  9  10 11
//! [0, 10)  X  X  X  X  X  X  X  X  X  X
//! [2, 5)         X  X  X
//! [3, 8)            X  X  X  X  X
//! [3, 8)            X  X  X  X  X
//! [3, 8)            X  X  X  X  X
//! [3, 8)            X  X  X  X  X
//! [5, 9)                  X  X  X  X
//! [8, 11)                          X  X  X
//!
//! Query:  [8, 11)
//! Answer: [0, 10), [5, 9), [8, 11)
//! ```
//!
//! Most interaction with this crate will be through the [`Lapper`] struct. The
//! main methods are [`Lapper::find`], [`Lapper::seek`], and [`Lapper::count`].
//! `find()` handles independent queries, `seek()` reuses a caller-owned cursor
//! when query starts arrive in order, and `count()` is used when only the number
//! of overlaps is needed.
//!
//! Ranges are half-open: `[start, stop)`. Two ranges overlap when
//! `interval.start < query.stop` and `interval.stop > query.start`, so adjacent
//! ranges such as `[0, 10)` and `[10, 20)` do not overlap. This matches the
//! usual zero-based genomic coordinate system. Signed and unsigned primitive
//! coordinates are supported.
//!
//! Lapper does not use an interval tree. It keeps intervals sorted by start and
//! builds a small index over fixed blocks of 32 intervals. A prefix maximum
//! finds the first block that could overlap; each block's minimum and maximum
//! end positions then prove whether the block is a miss or a dense prefix, and
//! a next-greater link skips runs of blocks that cannot overlap. Mixed blocks
//! produce an exact 32-bit overlap mask with NEON, AVX2, or the scalar fallback.
//! Mask bits are drained from low to high, so results remain borrowed and in
//! ascending start order.
//!
//! The same block algorithm handles ordinary data and the old worst case where
//! one long interval engulfs many shorter intervals. There is no workload mode
//! to configure. `merge_overlaps()` remains useful when callers want merged
//! coverage, while `count()` remains the independent BITS implementation and is
//! fast regardless of interval shape.
//!
//! # Examples
//!
//! ```rust
//!    use rust_lapper::{Interval, Lapper};
//!    use std::cmp;
//!    type Iv = Interval<usize, u32>;
//!
//!    // create some fake data
//!    let data: Vec<Iv> = (0..20).step_by(5).map(|x| Iv{start: x, stop: x + 2, val: 0}).collect();
//!    println!("{:#?}", data);
//!
//!    // make lapper structure
//!    let laps = Lapper::new(data);
//!
//!    assert_eq!(laps.find(6, 11).next(), Some(&Iv{start: 5, stop: 7, val: 0}));
//!
//!    // Demonstration of seek function. By passing in the &mut cursor, seek can have thread local
//!    // cursors going
//!    let mut sim: usize = 0;
//!    let mut cursor = 0;
//!    // Calculate the overlap between the query and the found intervals, sum total overlap
//!    for i in (0..10).step_by(3) {
//!        sim += laps
//!            .seek(i, i + 2, &mut cursor)
//!            .map(|iv| cmp::min(i + 2, iv.stop) - cmp::max(i, iv.start))
//!            .sum::<usize>();
//!    }
//!    assert_eq!(sim, 4);
//! ```
use num_traits::{
    identities::{one, zero},
    PrimInt,
};
use std::cmp::Ordering::{self};
use std::collections::VecDeque;

mod simd;

use simd::{detect_backend, overlap_mask, MaskBackend, BLOCK_SIZE as INDEX_BLOCK_SIZE};

#[cfg(feature = "with_serde")]
use serde::{ser::SerializeStruct, Deserialize, Deserializer, Serialize, Serializer};

/// Represent a range from [start, stop)
/// Inclusive start, exclusive of stop
#[cfg_attr(feature = "with_serde", derive(Serialize, Deserialize))]
#[derive(Eq, Debug, Clone)]
pub struct Interval<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    pub start: I,
    pub stop: I,
    pub val: T,
}

/// Primary interval collection and query index.
///
/// The public interval vector is the canonical storage and can be read or used
/// to mutate payload values. Coordinate or structural changes must use
/// [`Lapper::insert`] or [`Lapper::merge_overlaps`] so the private query index
/// is rebuilt.
#[derive(Debug, Clone)]
pub struct Lapper<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    /// Intervals in ascending start order.
    ///
    /// Directly changing coordinates or vector length leaves the private query
    /// index stale. Payload-only changes are safe.
    pub intervals: Vec<Interval<I, T>>,
    /// Sorted list of start positions,
    starts: Vec<I>,
    /// Sorted list of end positions,
    stops: Vec<I>,
    /// End positions in the same order as `intervals`, for block-mask queries.
    stops_by_start: Vec<I>,
    /// Index of the next block with a strictly larger maximum end, or the
    /// number of blocks when no such block exists.
    block_index: Vec<usize>,
    /// Maximum end position in each fixed-size block.
    block_max_ends: Vec<I>,
    /// Minimum end position in each fixed-size block.
    block_min_ends: Vec<I>,
    /// Inclusive prefix maximum of `block_max_ends`, used to find the first
    /// candidate block for a query.
    block_prefix_max_ends: Vec<I>,
    /// The length of the longest interval
    max_len: I,
    /// Whether a valid interval length exceeds the positive range of `I`.
    max_len_overflowed: bool,
    /// The calculated number of positions covered by the intervals
    cov: Option<I>,
    /// Whether or not overlaps have been merged
    pub overlaps_merged: bool,
}

#[cfg(feature = "with_serde")]
impl<I, T> Serialize for Lapper<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync + Serialize,
    T: Eq + Clone + Send + Sync + Serialize,
{
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer,
    {
        let mut state = serializer.serialize_struct("Lapper", 6)?;
        state.serialize_field("intervals", &self.intervals)?;
        state.serialize_field("starts", &self.starts)?;
        state.serialize_field("stops", &self.stops)?;
        state.serialize_field("max_len", &self.max_len)?;
        state.serialize_field("cov", &self.cov)?;
        state.serialize_field("overlaps_merged", &self.overlaps_merged)?;
        state.end()
    }
}

#[cfg(feature = "with_serde")]
#[derive(Deserialize)]
#[allow(dead_code)]
struct SerializedLapper<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    intervals: Vec<Interval<I, T>>,
    starts: Vec<I>,
    stops: Vec<I>,
    max_len: I,
    cov: Option<I>,
    overlaps_merged: bool,
}

#[cfg(feature = "with_serde")]
impl<'de, I, T> Deserialize<'de> for Lapper<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync + Deserialize<'de> + 'static,
    T: Eq + Clone + Send + Sync + Deserialize<'de>,
{
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: Deserializer<'de>,
    {
        let serialized = SerializedLapper::<I, T>::deserialize(deserializer)?;
        let mut lapper = Self::new(serialized.intervals);
        lapper.cov = serialized.cov;
        lapper.overlaps_merged = serialized.overlaps_merged;
        Ok(lapper)
    }
}

impl<I, T> Interval<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    /// Compute the intsect between two intervals
    #[inline]
    pub fn intersect(&self, other: &Interval<I, T>) -> I {
        std::cmp::min(self.stop, other.stop)
            .checked_sub(std::cmp::max(&self.start, &other.start))
            .unwrap_or_else(zero::<I>)
    }

    /// Check if two intervals overlap
    #[inline]
    pub fn overlap(&self, start: I, stop: I) -> bool {
        self.start < stop && self.stop > start
    }
}

impl<I, T> Ord for Interval<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    #[inline]
    fn cmp(&self, other: &Interval<I, T>) -> Ordering {
        match self.start.cmp(&other.start) {
            Ordering::Less => Ordering::Less,
            Ordering::Greater => Ordering::Greater,
            Ordering::Equal => self.stop.cmp(&other.stop),
        }
    }
}

impl<I, T> PartialOrd for Interval<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    #[inline]
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl<I, T> PartialEq for Interval<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync,
    T: Eq + Clone + Send + Sync,
{
    #[inline]
    fn eq(&self, other: &Interval<I, T>) -> bool {
        self.start == other.start && self.stop == other.stop
    }
}

impl<I, T> Lapper<I, T>
where
    I: PrimInt + Ord + Clone + Send + Sync + 'static,
    T: Eq + Clone + Send + Sync,
{
    /// Create a new instance of Lapper by passing in a vector of Intervals. This vector will
    /// immediately be sorted by start order.
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data = (0..20).step_by(5)
    ///                   .map(|x| Interval{start: x, stop: x + 10, val: true})
    ///                   .collect::<Vec<Interval<usize, bool>>>();
    /// let lapper = Lapper::new(data);
    /// ```
    pub fn new(mut intervals: Vec<Interval<I, T>>) -> Self {
        #[cfg(feature = "sort_unstable")]
        intervals.sort_unstable();

        #[cfg(not(feature = "sort_unstable"))]
        intervals.sort();

        let mut lapper = Lapper {
            intervals,
            starts: Vec::new(),
            stops: Vec::new(),
            stops_by_start: Vec::new(),
            block_index: Vec::new(),
            block_max_ends: Vec::new(),
            block_min_ends: Vec::new(),
            block_prefix_max_ends: Vec::new(),
            max_len: zero::<I>(),
            max_len_overflowed: false,
            cov: None,
            overlaps_merged: false,
        };
        lapper.rebuild_derived();
        lapper
    }

    fn rebuild_derived(&mut self) {
        let (starts, stops_by_start): (Vec<_>, Vec<_>) = self
            .intervals
            .iter()
            .map(|interval| (interval.start, interval.stop))
            .unzip();
        self.starts = starts;
        self.stops = stops_by_start.clone();
        self.stops_by_start = stops_by_start;

        self.max_len = zero::<I>();
        self.max_len_overflowed = false;
        for interval in &self.intervals {
            match interval.stop.checked_sub(&interval.start) {
                Some(length) => self.max_len = std::cmp::max(self.max_len, length),
                None if interval.stop >= interval.start => self.max_len_overflowed = true,
                None => {}
            }
        }

        #[cfg(feature = "sort_unstable")]
        self.stops.sort_unstable();

        #[cfg(not(feature = "sort_unstable"))]
        self.stops.sort();

        self.block_max_ends.clear();
        self.block_min_ends.clear();
        let interval_count = self.intervals.len();
        let block_count =
            interval_count / INDEX_BLOCK_SIZE + usize::from(interval_count % INDEX_BLOCK_SIZE != 0);
        self.block_max_ends.reserve(block_count);
        self.block_min_ends.reserve(block_count);
        for block in self.intervals.chunks(INDEX_BLOCK_SIZE) {
            let mut max_end = block[0].stop;
            let mut min_end = max_end;
            for interval in &block[1..] {
                max_end = std::cmp::max(max_end, interval.stop);
                min_end = std::cmp::min(min_end, interval.stop);
            }
            self.block_max_ends.push(max_end);
            self.block_min_ends.push(min_end);
        }

        self.block_index.clear();
        self.block_index.resize(block_count, block_count);
        let mut stack = Vec::<usize>::new();
        for block in (0..block_count).rev() {
            while stack.last().map_or(false, |&next| {
                self.block_max_ends[next] <= self.block_max_ends[block]
            }) {
                stack.pop();
            }
            self.block_index[block] = stack.last().copied().unwrap_or(block_count);
            stack.push(block);
        }

        self.block_prefix_max_ends.clone_from(&self.block_max_ends);
        for block in 1..block_count {
            self.block_prefix_max_ends[block] = std::cmp::max(
                self.block_prefix_max_ends[block - 1],
                self.block_prefix_max_ends[block],
            );
        }
    }

    /// Insert a new interval after the Lapper has been created. This is very
    /// inefficient and should be avoided if possible.
    ///
    /// SIDE EFFECTS: This clears cov() and overlaps_merged
    /// meaning that those will have to be recomputed after a insert
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data : Vec<Interval<usize, usize>>= vec!{
    ///     Interval{start:0,  stop:5,  val:1},
    ///     Interval{start:6,  stop:10, val:2},
    /// };
    /// let mut lapper = Lapper::new(data);
    /// lapper.insert(Interval{start:0, stop:20, val:5});
    /// assert_eq!(lapper.len(), 3);
    /// assert_eq!(lapper.find(1,3).collect::<Vec<&Interval<usize,usize>>>(),
    ///     vec![
    ///         &Interval{start:0, stop:5, val:1},
    ///         &Interval{start:0, stop:20, val:5},
    ///     ]
    /// );
    ///
    /// ```
    pub fn insert(&mut self, elem: Interval<I, T>) {
        let intervals_insert_index = Self::bsearch_seq_ref(&elem, &self.intervals);
        self.intervals.insert(intervals_insert_index, elem);
        self.rebuild_derived();
        self.cov = None;
        self.overlaps_merged = false;
    }

    /// Get the number over intervals in Lapper
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data = (0..20).step_by(5)
    ///                   .map(|x| Interval{start: x, stop: x + 10, val: true})
    ///                   .collect::<Vec<Interval<usize, bool>>>();
    /// let lapper = Lapper::new(data);
    /// assert_eq!(lapper.len(), 4);
    /// ```
    #[inline]
    pub fn len(&self) -> usize {
        self.intervals.len()
    }

    /// Check if lapper is empty
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data: Vec<Interval<usize, bool>> = vec![];
    /// let lapper = Lapper::new(data);
    /// assert_eq!(lapper.is_empty(), true);
    /// ```
    #[inline]
    pub fn is_empty(&self) -> bool {
        self.intervals.is_empty()
    }

    /// Get the number of positions covered by the intervals in Lapper. This provides immutable
    /// access if it has already been set, or on the fly calculation.
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data = (0..20).step_by(5)
    ///                   .map(|x| Interval{start: x, stop: x + 10, val: true})
    ///                   .collect::<Vec<Interval<usize, bool>>>();
    /// let lapper = Lapper::new(data);
    /// assert_eq!(lapper.cov(), 25);
    #[inline]
    pub fn cov(&self) -> I {
        match self.cov {
            None => self.calculate_coverage(),
            Some(cov) => cov,
        }
    }

    /// Get the number of positions covered by the intervals in Lapper and store it. If you are
    /// going to be using the coverage, you should set it to avoid calculating it over and over.
    pub fn set_cov(&mut self) -> I {
        let cov = self.calculate_coverage();
        self.cov = Some(cov);
        cov
    }

    /// Calculate the actual coverage behind the scenes.
    fn calculate_coverage(&self) -> I {
        let mut moving_interval = Interval {
            start: zero::<I>(),
            stop: zero::<I>(),
            val: zero::<I>(),
        };
        let mut cov = zero::<I>();

        for interval in self.intervals.iter() {
            // If it overlaps, embrace, extend, extinguish
            if moving_interval.overlap(interval.start, interval.stop) {
                moving_interval.start = std::cmp::min(moving_interval.start, interval.start);
                moving_interval.stop = std::cmp::max(moving_interval.stop, interval.stop);
            } else {
                // add the set and move on
                cov = cov + (moving_interval.stop - moving_interval.start);
                moving_interval.start = interval.start;
                moving_interval.stop = interval.stop;
            }
        }
        // add in the last bit
        cov = cov + (moving_interval.stop - moving_interval.start);
        cov
    }

    /// Return an iterator over the intervals in Lapper
    #[inline]
    pub fn iter(&self) -> IterLapper<'_, I, T> {
        IterLapper {
            inner: self,
            pos: 0,
        }
    }

    /// Merge any intervals that overlap with each other within the Lapper. This is an easy way to
    /// speed up queries.
    pub fn merge_overlaps(&mut self) {
        let mut stack: VecDeque<&mut Interval<I, T>> = VecDeque::new();
        let mut ivs = self.intervals.iter_mut();
        if let Some(first) = ivs.next() {
            stack.push_back(first);
            for interval in ivs {
                let top = stack.pop_back().unwrap();
                if top.stop < interval.start {
                    stack.push_back(top);
                    stack.push_back(interval);
                } else if top.stop < interval.stop {
                    top.stop = interval.stop;
                    //stack.pop_back();
                    stack.push_back(top);
                } else {
                    // they were equal
                    stack.push_back(top);
                }
            }
            self.overlaps_merged = true;
            self.intervals = stack
                .into_iter()
                .map(|x| Interval {
                    start: x.start,
                    stop: x.stop,
                    val: x.val.clone(),
                })
                .collect();
        }
        self.rebuild_derived();
    }

    /// Determine the first index that we should start checking for overlaps for via a binary
    /// search.
    /// Assumes that the maximum interval length in `intervals` has been subtracted from
    /// `start`, otherwise the result is undefined
    #[inline]
    pub fn lower_bound(start: I, intervals: &[Interval<I, T>]) -> usize {
        let mut size = intervals.len();
        let mut low = 0;

        while size > 0 {
            let half = size / 2;
            let other_half = size - half;
            let probe = low + half;
            let other_low = low + other_half;
            let v = &intervals[probe];
            size = half;
            low = if v.start < start { other_low } else { low }
        }
        low
    }

    #[inline]
    pub fn bsearch_seq<K>(key: K, elems: &[K]) -> usize
    where
        K: PartialEq + PartialOrd,
    {
        Self::bsearch_seq_ref(&key, elems)
    }

    #[inline]
    pub fn bsearch_seq_ref<K>(key: &K, elems: &[K]) -> usize
    where
        K: PartialEq + PartialOrd,
    {
        if elems.is_empty() || elems[0] >= *key {
            return 0;
        } else if elems[elems.len() - 1] < *key {
            return elems.len();
        }

        let mut cursor = 0;
        let mut length = elems.len();
        while length > 1 {
            let half = length >> 1;
            length -= half;
            cursor += (usize::from(elems[cursor + half - 1] < *key)) * half;
        }
        cursor
    }

    /// Return the number of positions in the union and intersection of two Lappers.
    ///
    /// The union counts each position covered by either Lapper once. The intersection counts each
    /// position covered by both Lappers once, regardless of how many intervals cover it.
    /// ``` rust
    /// use rust_lapper::{Lapper, Interval};
    /// type Iv = Interval<u32, u32>;
    /// let data1: Vec<Iv> = vec![
    ///     Iv{start: 70, stop: 120, val: 0}, // a long interval
    ///     Iv{start: 10, stop: 15, val: 0}, // exact overlap
    ///     Iv{start: 12, stop: 15, val: 0}, // inner overlap
    ///     Iv{start: 14, stop: 16, val: 0}, // overlap end
    ///     Iv{start: 68, stop: 71, val: 0}, // overlap start
    /// ];
    /// let data2: Vec<Iv> = vec![
    ///
    ///     Iv{start: 10, stop: 15, val: 0},
    ///     Iv{start: 40, stop: 45, val: 0},
    ///     Iv{start: 50, stop: 55, val: 0},
    ///     Iv{start: 60, stop: 65, val: 0},
    ///     Iv{start: 70, stop: 75, val: 0},
    /// ];
    ///
    /// let (mut lapper1, mut lapper2) = (Lapper::new(data1), Lapper::new(data2)) ;
    /// // Should be the same either way it's calculated
    /// let (union, intersect) = lapper1.union_and_intersect(&lapper2);
    /// assert_eq!(intersect, 10);
    /// assert_eq!(union, 73);
    /// let (union, intersect) = lapper2.union_and_intersect(&lapper1);
    /// assert_eq!(intersect, 10);
    /// assert_eq!(union, 73);
    /// lapper1.merge_overlaps();
    /// lapper1.set_cov();
    /// lapper2.merge_overlaps();
    /// lapper2.set_cov();
    ///
    /// // Should be the same either way it's calculated
    /// let (union, intersect) = lapper1.union_and_intersect(&lapper2);
    /// assert_eq!(intersect, 10);
    /// assert_eq!(union, 73);
    /// let (union, intersect) = lapper2.union_and_intersect(&lapper1);
    /// assert_eq!(intersect, 10);
    /// assert_eq!(union, 73);
    /// ```
    #[inline]
    pub fn union_and_intersect(&self, other: &Self) -> (I, I) {
        let mut cursor: usize = 0;

        if !self.overlaps_merged || !other.overlaps_merged {
            let mut intersections: Vec<Interval<I, bool>> = vec![];
            for self_iv in self.iter() {
                for other_iv in other.seek(self_iv.start, self_iv.stop, &mut cursor) {
                    let start = std::cmp::max(self_iv.start, other_iv.start);
                    let stop = std::cmp::min(self_iv.stop, other_iv.stop);
                    intersections.push(Interval {
                        start,
                        stop,
                        val: true,
                    });
                }
            }
            let mut temp_lapper = Lapper::new(intersections);
            temp_lapper.merge_overlaps();
            temp_lapper.set_cov();
            let union = self.cov() + other.cov() - temp_lapper.cov();
            (union, temp_lapper.cov())
        } else {
            let mut intersect = zero::<I>();
            for c1_iv in self.iter() {
                for c2_iv in other.seek(c1_iv.start, c1_iv.stop, &mut cursor) {
                    let local_intersect = c1_iv.intersect(c2_iv);
                    intersect = intersect + local_intersect;
                }
            }
            let union = self.cov() + other.cov() - intersect;
            (union, intersect)
        }
    }

    /// Find the intersect of two lapper objects.
    /// Intersect: The number of positions where both lappers intersect. Note that a position only
    /// counts one time, multiple Intervals covering the same position don't add up
    #[inline]
    pub fn intersect(&self, other: &Self) -> I {
        self.union_and_intersect(other).1
    }

    /// Find the union of two lapper objects.
    #[inline]
    pub fn union(&self, other: &Self) -> I {
        self.union_and_intersect(other).0
    }

    /// Return the contiguous intervals of coverage, `val` represents the number of intervals
    /// covering the returned interval.
    ///
    /// # Examples
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let data = (0..20).step_by(5)
    ///                   .map(|x| Interval{start: x, stop: x + 10, val: true})
    ///                   .collect::<Vec<Interval<usize, bool>>>();
    /// let lapper = Lapper::new(data);
    /// assert_eq!(lapper.depth().collect::<Vec<Interval<usize, usize>>>(), vec![
    ///             Interval { start: 0, stop: 5, val: 1 },
    ///             Interval { start: 5, stop: 20, val: 2 },
    ///             Interval { start: 20, stop: 25, val: 1 }]);
    /// ```
    #[inline]
    pub fn depth(&self) -> IterDepth<'_, I, T> {
        let mut merged_lapper = Lapper::new(
            self.intervals
                .iter()
                .map(|i| Interval {
                    start: i.start,
                    stop: i.stop,
                    val: true,
                })
                .collect::<Vec<Interval<I, bool>>>(),
        );
        merged_lapper.merge_overlaps();
        let merged_len = merged_lapper.intervals.len();
        IterDepth {
            inner: self,
            merged: merged_lapper,
            curr_merged_pos: zero::<I>(),
            initialized: false,
            curr_pos: 0,
            cursor: 0,
            end: merged_len,
        }
    }

    /// Count all intervals that overlap the half-open query `[start, stop)`.
    /// This performs two binary searches in order to
    /// find all the excluded elements, and then deduces the intersection from there. See
    /// [BITS](https://arxiv.org/pdf/1208.3407.pdf) for more details.
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let lapper = Lapper::new((0..100).step_by(5)
    ///                                 .map(|x| Interval{start: x, stop: x+2 , val: true})
    ///                                 .collect::<Vec<Interval<usize, bool>>>());
    /// assert_eq!(lapper.count(5, 11), 2);
    /// ```
    #[inline]
    pub fn count(&self, start: I, stop: I) -> usize {
        let len = self.intervals.len();
        let first = self
            .stops
            .partition_point(|interval_stop| *interval_stop <= start);
        let last = Self::bsearch_seq(stop, &self.starts);
        let num_cant_after = len - last;
        len - first - num_cant_after
    }

    /// Find all intervals that overlap the half-open query `[start, stop)`.
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let lapper = Lapper::new((0..100).step_by(5)
    ///                                 .map(|x| Interval{start: x, stop: x+2 , val: true})
    ///                                 .collect::<Vec<Interval<usize, bool>>>());
    /// assert_eq!(lapper.find(5, 11).count(), 2);
    /// ```
    #[inline]
    pub fn find(&self, start: I, stop: I) -> IterFind<'_, I, T> {
        let off = self
            .block_prefix_max_ends
            .partition_point(|max_end| *max_end <= start)
            * INDEX_BLOCK_SIZE;
        IterFind {
            inner: self,
            next_block_start: off,
            mask_block_start: 0,
            mask: 0,
            dense_next: 0,
            dense_end: 0,
            backend: detect_backend(),
            start,
            stop,
        }
    }

    /// Find all intervals that overlap the half-open query `[start, stop)`.
    ///
    /// Use this method when query starts arrive in nondecreasing order. A caller-owned cursor
    /// narrows the first candidate block, after which `seek()` uses the same block traversal as
    /// [`Lapper::find`]. Keeping the cursor outside `Lapper` allows immutable queries and preserves
    /// `Sync` when `T` and `I` are `Sync`.
    /// ```
    /// use rust_lapper::{Lapper, Interval};
    /// let lapper = Lapper::new((0..100).step_by(5)
    ///                                 .map(|x| Interval{start: x, stop: x+2 , val: true})
    ///                                 .collect::<Vec<Interval<usize, bool>>>());
    /// let mut cursor = 0;
    /// for i in lapper.iter() {
    ///    assert_eq!(lapper.seek(i.start, i.stop, &mut cursor).count(), 1);
    /// }
    /// ```
    #[inline]
    pub fn seek<'a>(&'a self, start: I, stop: I, cursor: &mut usize) -> IterFind<'a, I, T> {
        if self.max_len_overflowed {
            *cursor = self
                .block_prefix_max_ends
                .partition_point(|max_end| *max_end <= start)
                * INDEX_BLOCK_SIZE;
        } else {
            let earliest_start = start
                .checked_sub(&self.max_len)
                .unwrap_or_else(I::min_value);
            if *cursor == 0
                || (*cursor < self.intervals.len() && self.intervals[*cursor].start > start)
            {
                *cursor = Self::lower_bound(earliest_start, &self.intervals);
            }

            while *cursor + 1 < self.intervals.len()
                && self.intervals[*cursor + 1].start < earliest_start
            {
                *cursor += 1;
            }
        }

        IterFind {
            inner: self,
            next_block_start: (*cursor / INDEX_BLOCK_SIZE) * INDEX_BLOCK_SIZE,
            mask_block_start: 0,
            mask: 0,
            dense_next: 0,
            dense_end: 0,
            backend: detect_backend(),
            start,
            stop,
        }
    }
}

/// Find Iterator
#[derive(Debug)]
pub struct IterFind<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    inner: &'a Lapper<I, T>,
    next_block_start: usize,
    mask_block_start: usize,
    mask: u32,
    dense_next: usize,
    dense_end: usize,
    backend: MaskBackend,
    start: I,
    stop: I,
}

impl<'a, I, T> IterFind<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync + 'static,
{
    #[inline(always)]
    fn next_blockwise(&mut self) -> Option<&'a Interval<I, T>> {
        loop {
            // Return the next pending match from the dense branch.
            if self.dense_next < self.dense_end {
                let index = self.dense_next;
                self.dense_next += 1;
                return Some(&self.inner.intervals[index]);
            }

            // Return the next pending match from a mixed block's mask.
            if self.mask != 0 {
                let lane = self.mask.trailing_zeros() as usize;
                self.mask &= self.mask - 1;
                return Some(&self.inner.intervals[self.mask_block_start + lane]);
            }

            let block_start = self.next_block_start;
            // No blocks remain, so there are no more matches.
            if block_start >= self.inner.starts.len() {
                return None;
            }

            // Private arrays are rebuilt together by constructors, mutations,
            // and deserialization. Public interval edits cannot extend this bound.
            // Starts are sorted, so no interval in this or a later block can match.
            if unsafe { *self.inner.starts.get_unchecked(block_start) } >= self.stop {
                return None;
            }

            let block = block_start / INDEX_BLOCK_SIZE;
            // No ends in this block reach past the query start. Skip to the next
            // block whose maximum end might match.
            if unsafe { *self.inner.block_max_ends.get_unchecked(block) } <= self.start {
                let next_block = unsafe { *self.inner.block_index.get_unchecked(block) };
                debug_assert!(
                    next_block <= self.inner.block_index.len(),
                    "Lapper block index is corrupt"
                );
                self.next_block_start = next_block * INDEX_BLOCK_SIZE;
                continue;
            }

            let block_end = (block_start + INDEX_BLOCK_SIZE).min(self.inner.starts.len());
            // All ends match. Take the dense branch for the prefix whose starts
            // are before the query stop.
            if unsafe { *self.inner.block_min_ends.get_unchecked(block) } > self.start {
                let starts = unsafe { self.inner.starts.get_unchecked(block_start..block_end) };
                // If the last start matches, the whole block is dense. Otherwise,
                // find the matching prefix.
                let active_len = if unsafe { *starts.get_unchecked(starts.len() - 1) } < self.stop {
                    starts.len()
                } else {
                    starts.partition_point(|lane_start| *lane_start < self.stop)
                };
                self.dense_next = block_start;
                self.dense_end = block_start + active_len;
                // Continue after a full block. A partial prefix means every later
                // start is outside the query.
                self.next_block_start = if active_len == starts.len() {
                    block_end
                } else {
                    self.inner.starts.len()
                };
                continue;
            }

            // Some ends match. Build the exact overlap mask for this mixed block.
            self.mask_block_start = block_start;
            self.next_block_start = block_end;
            self.mask = overlap_mask(
                self.backend,
                unsafe { self.inner.starts.get_unchecked(block_start..block_end) },
                unsafe {
                    self.inner
                        .stops_by_start
                        .get_unchecked(block_start..block_end)
                },
                self.start,
                self.stop,
            );
        }
    }
}

impl<'a, I, T> Iterator for IterFind<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync + 'static,
{
    type Item = &'a Interval<I, T>;

    #[inline]
    fn next(&mut self) -> Option<Self::Item> {
        self.next_blockwise()
    }
}

/// Depth Iterator
#[derive(Debug)]
pub struct IterDepth<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    inner: &'a Lapper<I, T>,
    merged: Lapper<I, bool>, // A lapper that is the merged_lapper of inner
    curr_merged_pos: I,      // Current start position in current interval
    initialized: bool,
    curr_pos: usize, // In merged list of non-overlapping intervals
    cursor: usize,   // cursor for seek over inner lapper
    end: usize,      // len of merged
}

impl<'a, I, T> Iterator for IterDepth<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync + 'static,
{
    type Item = Interval<I, I>;

    #[inline]
    fn next(&mut self) -> Option<Self::Item> {
        let mut interval: &Interval<I, bool> = &self.merged.intervals[self.curr_pos];
        if !self.initialized {
            self.curr_merged_pos = interval.start;
            self.initialized = true;
        }
        if interval.stop == self.curr_merged_pos {
            if self.curr_pos + 1 != self.end {
                self.curr_pos += 1;
                interval = &self.merged.intervals[self.curr_pos];
                self.curr_merged_pos = interval.start;
            } else {
                return None;
            }
        }
        let start = self.curr_merged_pos;
        let depth_at_point = self
            .inner
            .seek(
                self.curr_merged_pos,
                self.curr_merged_pos + one::<I>(),
                &mut self.cursor,
            )
            .count();
        let mut new_depth_at_point = depth_at_point;
        while new_depth_at_point == depth_at_point && self.curr_merged_pos < interval.stop {
            self.curr_merged_pos = self.curr_merged_pos + one::<I>();
            if self.curr_merged_pos == interval.stop {
                break;
            }
            new_depth_at_point = self
                .inner
                .seek(
                    self.curr_merged_pos,
                    self.curr_merged_pos + one::<I>(),
                    &mut self.cursor,
                )
                .count();
        }
        Some(Interval {
            start,
            stop: self.curr_merged_pos,
            val: I::from(depth_at_point).unwrap(), // from usize should always work
        })
    }
}
/// Lapper Iterator
pub struct IterLapper<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    inner: &'a Lapper<I, T>,
    pos: usize,
}

impl<'a, I, T> Iterator for IterLapper<'a, I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    type Item = &'a Interval<I, T>;

    fn next(&mut self) -> Option<Self::Item> {
        if self.pos >= self.inner.intervals.len() {
            None
        } else {
            self.pos += 1;
            self.inner.intervals.get(self.pos - 1)
        }
    }
}

impl<I, T> IntoIterator for Lapper<I, T>
where
    T: Eq + Clone + Send + Sync,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    type Item = Interval<I, T>;
    type IntoIter = ::std::vec::IntoIter<Self::Item>;

    fn into_iter(self) -> Self::IntoIter {
        self.intervals.into_iter()
    }
}

impl<'a, I, T> IntoIterator for &'a Lapper<I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    type Item = &'a Interval<I, T>;
    type IntoIter = std::slice::Iter<'a, Interval<I, T>>;

    fn into_iter(self) -> std::slice::Iter<'a, Interval<I, T>> {
        self.intervals.iter()
    }
}

impl<'a, I, T> IntoIterator for &'a mut Lapper<I, T>
where
    T: Eq + Clone + Send + Sync + 'a,
    I: PrimInt + Ord + Clone + Send + Sync,
{
    type Item = &'a mut Interval<I, T>;
    type IntoIter = std::slice::IterMut<'a, Interval<I, T>>;

    fn into_iter(self) -> std::slice::IterMut<'a, Interval<I, T>> {
        self.intervals.iter_mut()
    }
}

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

    type Iv = Interval<usize, u32>;
    fn setup_nonoverlapping() -> Lapper<usize, u32> {
        let data: Vec<Iv> = (0..100)
            .step_by(20)
            .map(|x| Iv {
                start: x,
                stop: x + 10,
                val: 0,
            })
            .collect();
        Lapper::new(data)
    }

    fn setup_overlapping() -> Lapper<usize, u32> {
        let data: Vec<Iv> = (0..100)
            .step_by(10)
            .map(|x| Iv {
                start: x,
                stop: x + 15,
                val: 0,
            })
            .collect();
        Lapper::new(data)
    }

    fn setup_badlapper() -> Lapper<usize, u32> {
        let data: Vec<Iv> = vec![
            Iv{start: 70, stop: 120, val: 0}, // max_len = 50
            Iv{start: 10, stop: 15, val: 0},
            Iv{start: 10, stop: 15, val: 0}, // exact overlap
            Iv{start: 12, stop: 15, val: 0}, // inner overlap
            Iv{start: 14, stop: 16, val: 0}, // overlap end
            Iv{start: 40, stop: 45, val: 0},
            Iv{start: 50, stop: 55, val: 0},
            Iv{start: 60, stop: 65, val: 0},
            Iv{start: 68, stop: 71, val: 0}, // overlap start
            Iv{start: 70, stop: 75, val: 0},
        ];
        Lapper::new(data)
    }

    fn setup_single() -> Lapper<usize, u32> {
        let data: Vec<Iv> = vec![Iv {
            start: 10,
            stop: 35,
            val: 0,
        }];
        Lapper::new(data)
    }

    // Test that inserting data ends up with the same lapper (nonoverlapping)
    #[test]
    fn insert_equality_nonoverlapping() {
        let data: Vec<Iv> = (0..100)
            .step_by(20)
            .map(|x| Iv {
                start: x,
                stop: x + 10,
                val: 0,
            })
            .collect();
        let new_lapper = Lapper::new(data.clone());
        let mut insert_lapper = Lapper::new(vec![]);
        for elem in data {
            insert_lapper.insert(elem);
        }
        assert_eq!(new_lapper.starts, insert_lapper.starts);
        assert_eq!(new_lapper.stops, insert_lapper.stops);
        assert_eq!(new_lapper.intervals, insert_lapper.intervals);
        assert_eq!(new_lapper.max_len, insert_lapper.max_len);
    }

    // Test that inserting data ends up with the same lapper (overlapping)
    #[test]
    fn insert_equality_overlapping() {
        let data: Vec<Iv> = (0..100)
            .step_by(10)
            .map(|x| Iv {
                start: x,
                stop: x + 15,
                val: 0,
            })
            .collect();
        let new_lapper = Lapper::new(data.clone());
        let mut insert_lapper = Lapper::new(vec![]);
        for elem in data {
            insert_lapper.insert(elem);
        }
        assert_eq!(new_lapper.starts, insert_lapper.starts);
        assert_eq!(new_lapper.stops, insert_lapper.stops);
        assert_eq!(new_lapper.intervals, insert_lapper.intervals);
        assert_eq!(new_lapper.max_len, insert_lapper.max_len);
    }

    // Test that inserting data half with new and half with insert
    // ends up with the same lapper
    #[test]
    fn insert_equality_half_and_half() {
        let data: Vec<Iv> = (0..100)
            .step_by(1)
            .map(|x| Iv {
                start: x,
                stop: x + 15,
                val: 0,
            })
            .collect();
        let new_lapper = Lapper::new(data.clone());
        let (new_data, insert_data) = data.split_at(50);
        let mut insert_lapper = Lapper::new(new_data.to_vec());
        let mut insert_data = insert_data.to_vec();
        insert_data.reverse();
        for elem in insert_data {
            insert_lapper.insert(elem);
        }
        assert_eq!(new_lapper.starts, insert_lapper.starts);
        assert_eq!(new_lapper.stops, insert_lapper.stops);
        assert_eq!(new_lapper.intervals, insert_lapper.intervals);
        assert_eq!(new_lapper.max_len, insert_lapper.max_len);
    }

    // Test that inserting data ends up with the same lapper (badlapper)
    #[test]
    fn insert_equality_badlapper() {
        let data: Vec<Iv> = vec![
            Iv{start: 70, stop: 120, val: 0}, // max_len = 50
            Iv{start: 10, stop: 15, val: 0},
            Iv{start: 10, stop: 15, val: 0}, // exact overlap
            Iv{start: 12, stop: 15, val: 0}, // inner overlap
            Iv{start: 14, stop: 16, val: 0}, // overlap end
            Iv{start: 40, stop: 45, val: 0},
            Iv{start: 50, stop: 55, val: 0},
            Iv{start: 60, stop: 65, val: 0},
            Iv{start: 68, stop: 71, val: 0}, // overlap start
            Iv{start: 70, stop: 75, val: 0},
        ];
        let new_lapper = Lapper::new(data.clone());
        let mut insert_lapper = Lapper::new(vec![]);
        for elem in data {
            insert_lapper.insert(elem);
        }
        assert_eq!(new_lapper.starts, insert_lapper.starts);
        assert_eq!(new_lapper.stops, insert_lapper.stops);
        assert_eq!(new_lapper.intervals, insert_lapper.intervals);
        assert_eq!(new_lapper.max_len, insert_lapper.max_len);
    }

    // Test that inserting data ends up with the same lapper (single)
    #[test]
    fn insert_equality_single() {
        let data: Vec<Iv> = vec![Iv {
            start: 10,
            stop: 35,
            val: 0,
        }];
        let new_lapper = Lapper::new(data.clone());
        let mut insert_lapper = Lapper::new(vec![]);
        for elem in data {
            insert_lapper.insert(elem);
        }
        assert_eq!(new_lapper.starts, insert_lapper.starts);
        assert_eq!(new_lapper.stops, insert_lapper.stops);
        assert_eq!(new_lapper.intervals, insert_lapper.intervals);
        assert_eq!(new_lapper.max_len, insert_lapper.max_len);
    }

    // Test that a query stop that hits an interval start returns no interval
    #[test]
    fn test_query_stop_interval_start() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        assert_eq!(None, lapper.find(15, 20).next());
        assert_eq!(None, lapper.seek(15, 20, &mut cursor).next());
        assert_eq!(lapper.find(15, 20).count(), lapper.count(15, 20));
    }

    // Test that a query start that hits an interval end returns no interval
    #[test]
    fn test_query_start_interval_stop() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        assert_eq!(None, lapper.find(30, 35).next());
        assert_eq!(None, lapper.seek(30, 35, &mut cursor).next());
        assert_eq!(lapper.find(30, 35).count(), lapper.count(30, 35));
    }

    // Test that a query that overlaps the start of an interval returns that interval
    #[test]
    fn test_query_overlaps_interval_start() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        let expected = Iv {
            start: 20,
            stop: 30,
            val: 0,
        };
        assert_eq!(Some(&expected), lapper.find(15, 25).next());
        assert_eq!(Some(&expected), lapper.seek(15, 25, &mut cursor).next());
        assert_eq!(lapper.find(15, 25).count(), lapper.count(15, 25));
    }

    // Test that a query that overlaps the stop of an interval returns that interval
    #[test]
    fn test_query_overlaps_interval_stop() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        let expected = Iv {
            start: 20,
            stop: 30,
            val: 0,
        };
        assert_eq!(Some(&expected), lapper.find(25, 35).next());
        assert_eq!(Some(&expected), lapper.seek(25, 35, &mut cursor).next());
        assert_eq!(lapper.find(25, 35).count(), lapper.count(25, 35));
    }

    // Test that a query that is enveloped by interval returns interval
    #[test]
    fn test_interval_envelops_query() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        let expected = Iv {
            start: 20,
            stop: 30,
            val: 0,
        };
        assert_eq!(Some(&expected), lapper.find(22, 27).next());
        assert_eq!(Some(&expected), lapper.seek(22, 27, &mut cursor).next());
        assert_eq!(lapper.find(22, 27).count(), lapper.count(22, 27));
    }

    // Test that a query that envolops an interval returns that interval
    #[test]
    fn test_query_envolops_interval() {
        let lapper = setup_nonoverlapping();
        let mut cursor = 0;
        let expected = Iv {
            start: 20,
            stop: 30,
            val: 0,
        };
        assert_eq!(Some(&expected), lapper.find(15, 35).next());
        assert_eq!(Some(&expected), lapper.seek(15, 35, &mut cursor).next());
        assert_eq!(lapper.find(15, 35).count(), lapper.count(15, 35));
    }

    #[test]
    fn test_overlapping_intervals() {
        let lapper = setup_overlapping();
        let mut cursor = 0;
        let e1 = Iv {
            start: 0,
            stop: 15,
            val: 0,
        };
        let e2 = Iv {
            start: 10,
            stop: 25,
            val: 0,
        };
        assert_eq!(vec![&e1, &e2], lapper.find(8, 20).collect::<Vec<&Iv>>());
        assert_eq!(
            vec![&e1, &e2],
            lapper.seek(8, 20, &mut cursor).collect::<Vec<&Iv>>()
        );
        assert_eq!(lapper.count(8, 20), 2);
    }

    #[test]
    fn test_merge_overlaps() {
        let mut lapper = setup_badlapper();
        let expected: Vec<&Iv> = vec![
            &Iv{start: 10, stop: 16, val: 0},
            &Iv{start: 40, stop: 45, val: 0},
            &Iv{start: 50, stop: 55, val: 0},
            &Iv{start: 60, stop: 65, val: 0},
            &Iv{start: 68, stop: 120, val: 0}, // max_len = 50
        ];
        assert_eq!(lapper.intervals.len(), lapper.starts.len());
        lapper.merge_overlaps();
        assert_eq!(expected, lapper.iter().collect::<Vec<&Iv>>());
        assert_eq!(lapper.intervals.len(), lapper.starts.len())

    }

    // This test was added because this breakage was found in a library user's code, where after
    // calling merge_overlaps(), the find() call returned an empty iterator.
    #[test]
    fn test_merge_overlaps_find() {
        let data = vec![
                Iv{start: 2, stop: 3, val: 0},
                Iv{start: 3, stop: 4, val: 0},
                Iv{start: 4, stop: 6, val: 0},
                Iv{start: 6, stop: 7, val: 0},
                Iv{start: 7, stop: 8, val: 0},
        ];
        let mut lapper = Lapper::new(data);

        let found = lapper.find(7, 9).collect::<Vec<&Interval<_,_>>>();
        assert_eq!(found, vec![
            &Iv{start:7, stop: 8, val: 0},
        ]);

        // merge_overlaps should merge all intervals to one, which should be returned in the find call.
        lapper.merge_overlaps();

        let found = lapper.find(7, 9).collect::<Vec<&Interval<_,_>>>();
        assert_eq!(found, vec![
            &Iv{start:2, stop: 8, val: 0},
        ]);
    }

    #[test]
    fn test_lapper_cov() {
        let mut lapper = setup_badlapper();
        let before = lapper.cov();
        lapper.merge_overlaps();
        let after = lapper.cov();
        assert_eq!(before, after);

        let mut lapper = setup_nonoverlapping();
        lapper.set_cov();
        assert_eq!(lapper.cov(), 50);
    }

    #[test]
    fn test_interval_intersects() {
        let i1 = Iv{start: 70, stop: 120, val: 0}; // max_len = 50
        let i2 = Iv{start: 10, stop: 15, val: 0};
        let i3 = Iv{start: 10, stop: 15, val: 0}; // exact overlap
        let i4 = Iv{start: 12, stop: 15, val: 0}; // inner overlap
        let i5 = Iv{start: 14, stop: 16, val: 0}; // overlap end
        let i6 = Iv{start: 40, stop: 50, val: 0};
        let i7 = Iv{start: 50, stop: 55, val: 0};
        let i_8 = Iv{start: 60, stop: 65, val: 0};
        let i9 = Iv{start: 68, stop: 71, val: 0}; // overlap start
        let i10 = Iv{start: 70, stop: 75, val: 0};

        assert_eq!(i2.intersect(&i3), 5); // exact match
        assert_eq!(i2.intersect(&i4), 3); // inner intersect
        assert_eq!(i2.intersect(&i5), 1); // end intersect
        assert_eq!(i9.intersect(&i10), 1); // start intersect
        assert_eq!(i7.intersect(&i_8), 0); // no intersect
        assert_eq!(i6.intersect(&i7), 0); // no intersect stop = start
        assert_eq!(i1.intersect(&i10), 5); // inner intersect at start
    }

    #[test]
    fn test_union_and_intersect() {
        let data1: Vec<Iv> = vec![
            Iv{start: 70, stop: 120, val: 0}, // max_len = 50
            Iv{start: 10, stop: 15, val: 0}, // exact overlap
            Iv{start: 12, stop: 15, val: 0}, // inner overlap
            Iv{start: 14, stop: 16, val: 0}, // overlap end
            Iv{start: 68, stop: 71, val: 0}, // overlap start
        ];
        let data2: Vec<Iv> = vec![

            Iv{start: 10, stop: 15, val: 0},
            Iv{start: 40, stop: 45, val: 0},
            Iv{start: 50, stop: 55, val: 0},
            Iv{start: 60, stop: 65, val: 0},
            Iv{start: 70, stop: 75, val: 0},
        ];

        let (mut lapper1, mut lapper2) = (Lapper::new(data1), Lapper::new(data2)) ;
        // Should be the same either way it's calculated
        let (union, intersect) = lapper1.union_and_intersect(&lapper2);
        assert_eq!(intersect, 10);
        assert_eq!(union, 73);
        let (union, intersect) = lapper2.union_and_intersect(&lapper1);
        assert_eq!(intersect, 10);
        assert_eq!(union, 73);
        lapper1.merge_overlaps();
        lapper1.set_cov();
        lapper2.merge_overlaps();
        lapper2.set_cov();

        // Should be the same either way it's calculated
        let (union, intersect) = lapper1.union_and_intersect(&lapper2);
        assert_eq!(intersect, 10);
        assert_eq!(union, 73);
        let (union, intersect) = lapper2.union_and_intersect(&lapper1);
        assert_eq!(intersect, 10);
        assert_eq!(union, 73);
    }

    #[test]
    fn test_find_overlaps_in_large_intervals() {
        let data1: Vec<Iv> = vec![
            Iv{start: 0, stop: 8, val: 0},
            Iv{start: 1, stop: 10, val: 0},
            Iv{start: 2, stop: 5, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 4, stop: 7, val: 0},
            Iv{start: 5, stop: 8, val: 0},
            Iv{start: 8, stop: 8, val: 0},
            Iv{start: 9, stop: 11, val: 0},
            Iv{start: 10, stop: 13, val: 0},
            Iv{start: 100, stop: 200, val: 0},
            Iv{start: 110, stop: 120, val: 0},
            Iv{start: 110, stop: 124, val: 0},
            Iv{start: 111, stop: 160, val: 0},
            Iv{start: 150, stop: 200, val: 0},
        ];
        let lapper = Lapper::new(data1);
        let found = lapper.find(8, 11).collect::<Vec<&Iv>>();
        assert_eq!(found, vec![
            &Iv{start: 1, stop: 10, val: 0},
            &Iv{start: 9, stop: 11, val: 0},
            &Iv{start: 10, stop: 13, val: 0},
        ]);
        assert_eq!(lapper.count(8, 11), 3);
        let found = lapper.find(145, 151).collect::<Vec<&Iv>>();
        assert_eq!(found, vec![
            &Iv{start: 100, stop: 200, val: 0},
            &Iv{start: 111, stop: 160, val: 0},
            &Iv{start: 150, stop: 200, val: 0},
        ]);

        assert_eq!(lapper.count(145, 151), 3);
    }

    #[test]
    fn test_depth_sanity() {
        let data1: Vec<Iv> = vec![
            Iv{start: 0, stop: 10, val: 0},
            Iv{start: 5, stop: 10, val: 0}
        ];
        let lapper = Lapper::new(data1);
        let found = lapper.depth().collect::<Vec<Interval<usize, usize>>>();
        assert_eq!(found, vec![
                   Interval{start: 0, stop: 5, val: 1},
                   Interval{start: 5, stop: 10, val: 2}
        ]);
    }

    #[test]
    fn test_depth_hard() {
        let data1: Vec<Iv> = vec![
            Iv{start: 1, stop: 10, val: 0},
            Iv{start: 2, stop: 5, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 5, stop: 8, val: 0},
            Iv{start: 9, stop: 11, val: 0},
        ];
        let lapper = Lapper::new(data1);
        let found = lapper.depth().collect::<Vec<Interval<usize, usize>>>();
        assert_eq!(found, vec![
                   Interval{start: 1, stop: 2, val: 1},
                   Interval{start: 2, stop: 3, val: 2},
                   Interval{start: 3, stop: 8, val: 5},
                   Interval{start: 8, stop: 9, val: 1},
                   Interval{start: 9, stop: 10, val: 2},
                   Interval{start: 10, stop: 11, val: 1},
        ]);
    }
    #[test]
    fn test_depth_harder() {
        let data1: Vec<Iv> = vec![
            Iv{start: 1, stop: 10, val: 0},
            Iv{start: 2, stop: 5, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 3, stop: 8, val: 0},
            Iv{start: 5, stop: 8, val: 0},
            Iv{start: 9, stop: 11, val: 0},
            Iv{start: 15, stop: 20, val: 0},
        ];
        let lapper = Lapper::new(data1);
        let found = lapper.depth().collect::<Vec<Interval<usize, usize>>>();
        assert_eq!(found, vec![
                   Interval{start: 1, stop: 2, val: 1},
                   Interval{start: 2, stop: 3, val: 2},
                   Interval{start: 3, stop: 8, val: 5},
                   Interval{start: 8, stop: 9, val: 1},
                   Interval{start: 9, stop: 10, val: 2},
                   Interval{start: 10, stop: 11, val: 1},
                   Interval{start: 15, stop: 20, val: 1},
        ]);
    }
    // BUG TESTS - these are tests that came from real life

    // Test that it's not possible to induce index out of bounds by pushing the cursor past the end
    // of the lapper.
    #[test]
    fn test_seek_over_len() {
        let lapper = setup_nonoverlapping();
        let single = setup_single();
        let mut cursor: usize = 0;

        for interval in lapper.iter() {
            for o_interval in single.seek(interval.start, interval.stop, &mut cursor) {
                println!("{:#?}", o_interval);
            }
        }
    }

    // Test that if lower_bound puts us before the first match, we still return a match
    #[test]
    fn test_find_over_behind_first_match() {
        let lapper = setup_badlapper();
        let e1 = Iv {start: 50, stop: 55, val: 0};
        let found = lapper.find(50, 55).next();
        assert_eq!(found, Some(&e1));
        assert_eq!(lapper.find(50, 55).count(), lapper.count(50,55));
    }

    // When there is a very long interval that spans many little intervals, test that the little
    // Intervals still get returned properly.
    #[test]
    fn test_bad_skips() {
        let data = vec![
            Iv{start:25264912, stop: 25264986, val: 0},
            Iv{start:27273024, stop: 27273065	, val: 0},
            Iv{start:27440273, stop: 27440318	, val: 0},
            Iv{start:27488033, stop: 27488125	, val: 0},
            Iv{start:27938410, stop: 27938470	, val: 0},
            Iv{start:27959118, stop: 27959171	, val: 0},
            Iv{start:28866309, stop: 33141404	, val: 0},
        ];
        let lapper = Lapper::new(data);

        let found = lapper.find(28974798, 33141355).collect::<Vec<&Iv>>();
        assert_eq!(found, vec![
            &Iv{start:28866309, stop: 33141404	, val: 0},
        ]);
        assert_eq!(lapper.count(28974798, 33141355), 1);
    }

    #[cfg(feature = "with_serde")]
    #[derive(Debug, PartialEq, serde::Serialize, serde::Deserialize)]
    struct LegacyLapper {
        intervals: Vec<Iv>,
        starts: Vec<usize>,
        stops: Vec<usize>,
        max_len: usize,
        cov: Option<usize>,
        overlaps_merged: bool,
    }

    #[cfg(feature = "with_serde")]
    #[test]
    fn serde_keeps_the_v1_six_field_representation() {
        let data = vec![
            Iv{start:25264912, stop: 25264986, val: 0},
            Iv{start:27273024, stop: 27273065	, val: 0},
            Iv{start:27440273, stop: 27440318	, val: 0},
            Iv{start:27488033, stop: 27488125	, val: 0},
            Iv{start:27938410, stop: 27938470	, val: 0},
            Iv{start:27959118, stop: 27959171	, val: 0},
            Iv{start:28866309, stop: 33141404	, val: 0},
        ];
        let lapper = Lapper::new(data);

        let legacy = LegacyLapper {
            intervals: lapper.intervals.clone(),
            starts: lapper.starts.clone(),
            stops: lapper.stops.clone(),
            max_len: lapper.max_len,
            cov: lapper.cov,
            overlaps_merged: lapper.overlaps_merged,
        };
        let legacy_bytes = bincode::serialize(&legacy).unwrap();
        let deserialized: Lapper<usize, u32> = bincode::deserialize(&legacy_bytes).unwrap();
        let current_bytes = bincode::serialize(&deserialized).unwrap();
        assert_eq!(current_bytes, legacy_bytes);
        let legacy_again: LegacyLapper = bincode::deserialize(&current_bytes).unwrap();
        assert_eq!(legacy_again, legacy);

        let found = deserialized.find(28974798, 33141355).collect::<Vec<&Iv>>();
        assert_eq!(found, vec![
            &Iv{start:28866309, stop: 33141404	, val: 0},
        ]);
        assert_eq!(deserialized.count(28974798, 33141355), 1);
    }

}