pandora_lib_promethion/src/uniqueness/mod.rs

//! Genome-wide uniqueness index based on suffix arrays.
//!
//! This module computes, stores, and queries minimal unique sequence lengths
//! across a genome. Multi-contig references are indexed by concatenating the
//! selected contigs with reserved separator bytes, building one global suffix
//! array, and constraining all reported intervals to a single contig. The
//! coordinates exposed by [`GenomeIndex`] remain contig-local.
//!
//! For each genomic position `pos`, the index stores the shortest sequence
//! starting at `pos` that occurs exactly once in the genome:
//!
//! ```text
//! genome[pos .. pos + k)
//! ```
//!
//! where `k` is the minimal unique length at `pos`.
//!
//! Internally, uniqueness is computed from a suffix array and LCP
//! (Longest Common Prefix) array, allowing construction in `O(n log n)`
//! time and constant-time uniqueness queries.
//!
//! # Query types
//!
//! The index supports two complementary views of uniqueness:
//!
//! ## Unique intervals containing a position
//!
//! Given a genomic position, the index can return a unique interval that
//! contains that position.
//!
//! ```text
//! start <= pos < end
//! ```
//!
//! Several intervals may satisfy this condition; [`AnchorBias`] controls
//! how a candidate interval is selected.
//!
//! ## Independent left/right uniqueness
//!
//! The index can also report how much sequence is required on either side
//! of a position to obtain uniqueness independently.
//!
//! For a position `pos`:
//!
//! ```text
//! before: genome[pos - k .. pos)
//! after : genome[pos .. pos + k)
//! ```
//!
//! [`UniqueContext`] returns the smallest unique sequence ending at `pos`
//! (`before`) and the smallest unique sequence starting at `pos` (`after`).
//!
//! To support efficient left-side queries, the index stores uniqueness
//! information for both the forward genome and its reverse complement-free
//! reversal.
//!
//! # Caching
//!
//! Construction can be expensive for large genomes. The module therefore
//! supports disk-backed caching of precomputed uniqueness vectors.
//!
//! Cached data are keyed by genome metadata and loaded transparently when
//! available.
//!
//! # Main types
//!
//! - [`UniquenessIndex`] — uniqueness index for a single sequence.
//! - [`GenomeIndex`] — multi-contig genome-wide index.
//! - [`UniqueInterval`] — unique interval containing a queried position.
//! - [`UniqueContext`] — independent left/right uniqueness lengths.
//! - [`GenomicPos`] — genomic coordinate used by [`GenomeIndex`].

use log::{info, warn};
use rayon::prelude::*;
use std::io::{self};
use std::path::Path;
use std::time::Instant;

pub type GenomeCoord = u32;
pub type UniqueLen = u32;

pub const NOT_UNIQUE: UniqueLen = u32::MAX;

mod cache;
mod genome_index;
mod suffix;

pub use cache::CacheKey;
pub use genome_index::{ContigSpan, GenomeIndex, GenomicPos};

use cache::{load_from_cache, save_to_cache, MappedUniqueLenCache};
use suffix::build_suffix_array_and_lcp;

/// Pure computation: genome bytes → min_unique_len vector.
///
/// `contig_end_at[p]` is the exclusive end coordinate of the contig containing
/// `p`. Separator positions must use `contig_end_at[p] = p`, which makes every
/// separator position non-queryable.
fn build_min_unique_len(genome: &[u8], contig_end_at: &[GenomeCoord]) -> Vec<UniqueLen> {
    let n = genome.len();
    assert_eq!(n, contig_end_at.len());
    assert!(
        n <= u32::MAX as usize,
        "uniqueness index only supports sequences up to u32::MAX bases"
    );

    info!(
        "[uniqueness] computing minimal unique lengths ({:.1} Mbp)",
        n as f64 / 1e6
    );
    let total_started = Instant::now();

    let (sa, lcp) = build_suffix_array_and_lcp(genome);

    info!("[uniqueness] building suffix rank array");
    let started = Instant::now();
    let mut rank = vec![0u32; n];
    for (i, &s) in sa.iter().enumerate() {
        rank[s as usize] = i as u32;
    }
    drop(sa);
    info!(
        "[uniqueness] suffix rank array built in {:.1}s",
        started.elapsed().as_secs_f64()
    );

    info!("[uniqueness] computing boundary-aware unique length vector");
    let started = Instant::now();
    let min_unique_len = (0..n)
        .into_par_iter()
        .map(|p| {
            let r = rank[p] as usize;
            let left = if r > 0 { lcp[r] } else { 0 };
            let right = if r + 1 < n { lcp[r + 1] } else { 0 };
            let Some(k) = left.max(right).checked_add(1) else {
                return NOT_UNIQUE;
            };
            let Some(end) = (p as GenomeCoord).checked_add(k) else {
                return NOT_UNIQUE;
            };
            if end <= contig_end_at[p] {
                k
            } else {
                NOT_UNIQUE
            }
        })
        .collect();
    drop(rank);
    drop(lcp);
    info!(
        "[uniqueness] unique length vector computed in {:.1}s (total {:.1}s)",
        started.elapsed().as_secs_f64(),
        total_started.elapsed().as_secs_f64()
    );
    min_unique_len
}

enum UniqueLenStore {
    Owned(Vec<UniqueLen>),
    Mapped(MappedUniqueLenCache),
}

impl UniqueLenStore {
    fn get(&self, index: usize) -> UniqueLen {
        match self {
            Self::Owned(values) => values[index],
            Self::Mapped(values) => values.get(index),
        }
    }

    #[cfg(test)]
    fn len(&self) -> usize {
        match self {
            Self::Owned(values) => values.len(),
            Self::Mapped(values) => values.len(),
        }
    }

    #[cfg(test)]
    fn to_vec(&self) -> Vec<UniqueLen> {
        (0..self.len()).map(|idx| self.get(idx)).collect()
    }
}

#[derive(Debug, Clone, Copy)]
pub struct UniqueContext {
    /// Bases required before `pos`.
    pub before: Option<UniqueLen>,

    /// Bases required from `pos` onward.
    pub after: Option<UniqueLen>,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct UniqueInterval {
    pub start: GenomeCoord,
    pub end: GenomeCoord,
}

impl UniqueInterval {
    pub fn len(&self) -> UniqueLen {
        self.end - self.start
    }
}

#[derive(Debug, Clone, Copy)]
pub enum AnchorBias {
    Left,
    Right,
    Auto,
}

pub struct UniquenessIndex {
    min_unique_len: UniqueLenStore,
    rev_min_unique_len: UniqueLenStore,
    genome_len: GenomeCoord,
}

impl UniquenessIndex {
    /// Build from genome bytes, using a disk cache keyed on `cache_key`.
    /// Cache lives in `cache_dir` (created if absent).
    /// If the cache file exists and is valid, loading is O(n) read — no recomputation.
    pub fn build_with_cache(
        genome: &[u8],
        contig_end_at: &[GenomeCoord],
        rev_contig_end_at: &[GenomeCoord],
        key: &CacheKey,
        cache_dir: &Path,
    ) -> io::Result<Self> {
        assert_eq!(genome.len(), contig_end_at.len());
        assert_eq!(genome.len(), rev_contig_end_at.len());
        std::fs::create_dir_all(cache_dir)?;

        let fwd_path = cache_dir.join(key.file_name());
        let rev_path = cache_dir.join(format!("rev.{}", key.file_name()));

        let min_unique_len = match load_from_cache(&fwd_path, genome.len()) {
            Ok(mul) => {
                info!(
                    "[cache] forward uniqueness cache hit: {}",
                    fwd_path.display()
                );
                mul
            }
            Err(e) => {
                if fwd_path.exists() {
                    warn!("[cache] invalid forward ({e}), rebuilding...");
                } else {
                    info!(
                        "[cache] forward uniqueness cache miss: {}",
                        fwd_path.display()
                    );
                }

                let mul = build_min_unique_len(genome, contig_end_at);
                info!(
                    "[cache] saving forward uniqueness cache: {}",
                    fwd_path.display()
                );
                save_to_cache(&fwd_path, &mul)?;
                mul
            }
        };

        let rev_min_unique_len = match load_from_cache(&rev_path, genome.len()) {
            Ok(mul) => {
                info!(
                    "[cache] reverse uniqueness cache hit: {}",
                    rev_path.display()
                );
                mul
            }
            Err(e) => {
                if rev_path.exists() {
                    warn!("[cache] invalid reverse ({e}), rebuilding...");
                } else {
                    info!(
                        "[cache] reverse uniqueness cache miss: {}",
                        rev_path.display()
                    );
                }

                info!("[uniqueness] reversing genome for reverse uniqueness pass");
                let started = Instant::now();
                let mut rev = genome.to_vec();
                rev.reverse();
                info!(
                    "[uniqueness] reversed genome prepared in {:.1}s",
                    started.elapsed().as_secs_f64()
                );

                let mul = build_min_unique_len(&rev, rev_contig_end_at);
                info!(
                    "[cache] saving reverse uniqueness cache: {}",
                    rev_path.display()
                );
                save_to_cache(&rev_path, &mul)?;
                mul
            }
        };

        Ok(Self::from_raw(min_unique_len, rev_min_unique_len))
    }

    /// Construct directly from a precomputed `min_unique_len` vector.
    /// Used when loading from disk cache — skips SA/LCP computation entirely.
    pub fn from_raw(min_unique_len: Vec<UniqueLen>, rev_min_unique_len: Vec<UniqueLen>) -> Self {
        assert_eq!(min_unique_len.len(), rev_min_unique_len.len());
        assert!(
            min_unique_len.len() <= u32::MAX as usize,
            "uniqueness index only supports sequences up to u32::MAX bases"
        );

        let genome_len = min_unique_len.len() as u32;

        Self {
            min_unique_len: UniqueLenStore::Owned(min_unique_len),
            rev_min_unique_len: UniqueLenStore::Owned(rev_min_unique_len),
            genome_len,
        }
    }

    pub(crate) fn from_mapped(
        min_unique_len: MappedUniqueLenCache,
        rev_min_unique_len: MappedUniqueLenCache,
    ) -> Self {
        assert_eq!(min_unique_len.len(), rev_min_unique_len.len());
        assert!(
            min_unique_len.len() <= u32::MAX as usize,
            "uniqueness index only supports sequences up to u32::MAX bases"
        );

        let genome_len = min_unique_len.len() as u32;
        Self {
            min_unique_len: UniqueLenStore::Mapped(min_unique_len),
            rev_min_unique_len: UniqueLenStore::Mapped(rev_min_unique_len),
            genome_len,
        }
    }

    /// Build without cache (e.g. for tests or one-shot use).
    pub fn build(genome: &[u8]) -> Self {
        let n = genome.len() as GenomeCoord;
        let contig_end_at = vec![n; genome.len()];
        let min_unique_len = build_min_unique_len(genome, &contig_end_at);

        let mut rev = genome.to_vec();
        rev.reverse();
        let rev_contig_end_at = vec![n; genome.len()];
        let rev_min_unique_len = build_min_unique_len(&rev, &rev_contig_end_at);

        Self::from_raw(min_unique_len, rev_min_unique_len)
    }

    /// Returns the number of nucleotides before and after `pos` needed to form
    /// a unique sequence containing `pos`.
    ///
    /// The returned flank defines the half-open interval:
    ///
    /// `[pos - front, pos + back)`
    ///
    /// where `front` is the number of nucleotides before `pos`, and `back` is the
    /// number of nucleotides from `pos` onward, including the nucleotide at `pos`.
    ///
    /// Returns `None` if no unique sequence containing `pos` can be found.
    ///
    /// # Panics
    ///
    /// Panics if `pos >= self.genome_len()`.
    pub fn minimal_unique_interval_containing(
        &self,
        pos: GenomeCoord,
        bias: AnchorBias,
    ) -> Option<UniqueInterval> {
        assert!(pos < self.genome_len);
        let pos_usize = pos as usize;

        match bias {
            AnchorBias::Left => {
                let k = self.min_unique_len.get(pos_usize);
                if k == NOT_UNIQUE {
                    return None;
                }
                Some(UniqueInterval {
                    start: pos,
                    end: pos + k,
                })
            }

            AnchorBias::Right => {
                for a in (0..=pos).rev() {
                    let k = self.min_unique_len.get(a as usize);
                    if k == NOT_UNIQUE {
                        continue;
                    }
                    if a + k > pos {
                        return Some(UniqueInterval {
                            start: a,
                            end: a + k,
                        });
                    }
                    // a + k <= pos and a is decreasing: gap can only grow, bail
                    if pos - a >= self.genome_len {
                        break;
                    }
                }
                None
            }

            AnchorBias::Auto => {
                let window = 500;
                let mut best: Option<UniqueInterval> = None;
                for a in pos.saturating_sub(window)..=pos {
                    let k = self.min_unique_len.get(a as usize);
                    if k == NOT_UNIQUE {
                        continue;
                    }
                    if a + k <= pos {
                        continue;
                    }
                    let c = UniqueInterval {
                        start: a,
                        end: a + k,
                    };
                    best = Some(match best {
                        None => c,
                        Some(b) => {
                            if c.len() < b.len() {
                                c
                            } else {
                                b
                            }
                        }
                    });
                }
                best
            }
        }
    }

    pub fn genome_len(&self) -> GenomeCoord {
        self.genome_len
    }

    /// Returns how much sequence is needed before and after `pos`
    /// to obtain uniqueness independently on each side.
    ///
    /// `before = Some(k)` means:
    ///
    /// `genome[pos - k .. pos)` is the shortest unique sequence ending at `pos`.
    ///
    /// `after = Some(k)` means:
    ///
    /// `genome[pos .. pos + k)` is the shortest unique sequence starting at `pos`.
    ///
    /// Returns `None` for a side if no unique sequence exists on that side.
    ///
    /// # Panics
    ///
    /// Panics if `pos > self.genome_len()`.
    pub fn unique_context_at(&self, pos: GenomeCoord) -> UniqueContext {
        assert!(pos <= self.genome_len);

        let after = if pos == self.genome_len {
            None
        } else {
            match self.min_unique_len.get(pos as usize) {
                NOT_UNIQUE => None,
                k => Some(k),
            }
        };

        let before = if pos == 0 {
            None
        } else {
            let rev_pos = self.genome_len - pos;

            match self.rev_min_unique_len.get(rev_pos as usize) {
                NOT_UNIQUE => None,
                k => Some(k),
            }
        };

        UniqueContext { before, after }
    }
}

// tests/uniqueness.rs

#[cfg(test)]
mod tests {
    use std::{fs::File};
    use std::io::{BufWriter, Write};

    use rust_htslib::bam::{self, Read};

    use crate::{
        helpers::get_genome_sizes,
        uniqueness::genome_index::{GenomeIndex, GenomicPos},
    };

    use super::*;

    fn toy_genome() -> &'static [u8] {
        // "banana" — well-known SA example, easy to verify by hand
        b"ACGTACGTNNACGT"
    }

    fn single_contig_end_at(genome: &[u8]) -> Vec<GenomeCoord> {
        vec![genome.len() as GenomeCoord; genome.len()]
    }

    #[test]
    fn test_cache_roundtrip() {
        let genome = toy_genome();
        let dir = std::path::Path::new("/home/t_steimle/tmp/");

        let key = CacheKey::from_genome("chr_test", genome);

        let contig_end_at = single_contig_end_at(genome);
        let idx1 =
            UniquenessIndex::build_with_cache(genome, &contig_end_at, &contig_end_at, &key, dir)
                .unwrap();
        // Second call must hit cache
        let idx2 =
            UniquenessIndex::build_with_cache(genome, &contig_end_at, &contig_end_at, &key, dir)
                .unwrap();

        assert_eq!(idx1.min_unique_len.to_vec(), idx2.min_unique_len.to_vec());
    }

    #[test]
    fn test_cache_invalidated_on_length_mismatch() {
        let genome1 = b"ACGTACGT".as_slice();
        let genome2 = b"ACGTACGTNN".as_slice();
        let dir = std::path::Path::new("/home/t_steimle/tmp/");
        let key = CacheKey::from_genome("chr_test", genome1);

        let contig_end_at = single_contig_end_at(genome1);
        UniquenessIndex::build_with_cache(genome1, &contig_end_at, &contig_end_at, &key, dir)
            .unwrap();
        // genome2 has different len → load_from_cache returns Err → rebuild
        let key2 = CacheKey::from_genome("chr_test", genome2);
        // different hash → different file, no collision
        assert_ne!(key.file_name(), key2.file_name());
    }

    #[test]
    fn old_cache_format_is_rejected() {
        let path = std::env::temp_dir().join(format!(
            "pandora_old_uniqueness_cache_{}.uniqidx",
            std::process::id()
        ));
        let mut bytes = Vec::new();
        bytes.extend_from_slice(b"UNIQIDX1");
        bytes.extend_from_slice(&1u32.to_le_bytes());
        bytes.extend_from_slice(&1u64.to_le_bytes());
        bytes.extend_from_slice(&1u64.to_le_bytes());
        std::fs::write(&path, bytes).unwrap();

        assert!(cache::load_from_cache(&path, 1).is_err());
        let _ = std::fs::remove_file(path);
    }

    // ── helpers ───────────────────────────────────────────────────────────────

    /// Verify that a given interval is actually unique in the genome:
    /// the subsequence appears exactly once.
    fn assert_unique(genome: &[u8], iv: &UniqueInterval) {
        let needle = &genome[iv.start as usize..iv.end as usize];
        let count = genome
            .windows(needle.len())
            .filter(|w| *w == needle)
            .count();
        assert_eq!(
            count,
            1,
            "interval [{}, {}) = {:?} appears {count}x in genome — not unique",
            iv.start,
            iv.end,
            std::str::from_utf8(needle).unwrap_or("<binary>"),
        );
    }

    /// Verify that no strictly shorter sub-interval anchored at the same start
    /// is also unique (minimality check for Left bias).
    fn assert_minimal_left(genome: &[u8], iv: &UniqueInterval) {
        if iv.len() <= 1 {
            return;
        }
        let shorter = UniqueInterval {
            start: iv.start,
            end: iv.end - 1,
        };
        let needle = &genome[shorter.start as usize..shorter.end as usize];
        let count = genome
            .windows(needle.len())
            .filter(|w| *w == needle)
            .count();
        assert!(
            count > 1,
            "interval [{}, {}) = {:?} is shorter and already unique — minimality violated",
            shorter.start,
            shorter.end,
            std::str::from_utf8(needle).unwrap_or("<binary>"),
        );
    }

    /// Same for Right bias: no strictly shorter interval ending at the same end.
    fn assert_minimal_right(genome: &[u8], iv: &UniqueInterval) {
        if iv.len() <= 1 {
            return;
        }
        let shorter = UniqueInterval {
            start: iv.start + 1,
            end: iv.end,
        };
        let needle = &genome[shorter.start as usize..shorter.end as usize];
        let count = genome
            .windows(needle.len())
            .filter(|w| *w == needle)
            .count();
        assert!(
            count > 1,
            "interval [{}, {}) = {:?} is shorter and already unique — minimality violated",
            shorter.start,
            shorter.end,
            std::str::from_utf8(needle).unwrap_or("<binary>"),
        );
    }

    /// Assert that no interval of the same length covering pos is shorter
    /// (Auto minimality: no other anchor yields a shorter unique window).
    fn assert_minimal_auto(genome: &[u8], iv: &UniqueInterval, pos: GenomeCoord) {
        // Try all anchors in [pos - len + 1 .. pos] and verify none gives shorter
        let window = iv.len().saturating_sub(1);
        for a in pos.saturating_sub(window)..=pos {
            if a + 1 > genome.len() as GenomeCoord {
                break;
            }
            let max_end = (a + iv.len()).min(genome.len() as GenomeCoord);
            for end in (a + 1)..max_end {
                let needle = &genome[a as usize..end as usize];
                let count = genome
                    .windows(needle.len())
                    .filter(|w| *w == needle)
                    .count();
                if count == 1 {
                    panic!(
                        "found shorter unique interval [{a}, {end}) len={} \
                         covering pos={pos}, but Auto returned [{}, {}) len={}",
                        end - a,
                        iv.start,
                        iv.end,
                        iv.len()
                    );
                }
            }
        }
    }

    // ── genome fixtures ───────────────────────────────────────────────────────

    /// Simple genome with obvious unique regions.
    ///
    /// Layout (positions):
    ///  0123456789...
    ///  AAAA TTTT AAAA GCTAGCTA NNNN GCTAGCTA
    ///                 ^unique^       ^repeat^
    ///
    /// "GCTAGCTA" appears twice → not unique
    /// "TTTT"     appears once  → unique but short
    /// The N-run should return None for all biases.
    fn genome_simple() -> Vec<u8> {
        //0         1         2         3
        //0123456789012345678901234567890123456789
        b"AAAATTTTAAAAGCTAGCTANNNNNNNNNNNGCTAGCTA".to_vec()
    }

    /// Genome where Left/Right/Auto produce verifiably different intervals.
    ///
    ///  ACGT ACGT TGCA TGCA AAAA
    ///  0    4    8    12   16
    ///
    /// "ACGTACGT" covers [0..8) — only one occurrence at pos 0
    /// "TGCATGCA" covers [8..16) — only one occurrence
    /// Repeated motifs force longer unique strings.
    fn genome_repeats() -> Vec<u8> {
        b"ACGTACGTTGCATGCAAAAACCCCCGGGGG".to_vec()
    }

    /// Genome constructed so that Left, Right and Auto give three distinct
    /// intervals for a single query position.
    ///
    ///  ATCGATCG GGGG ATCGATCG TTTT ATCG CCCC
    ///  0        8    12       20   24   28
    ///
    /// pos=6 (inside first ATCGATCG repeat region):
    ///   Left  → must extend right to break the repeat
    ///   Right → extends left into earlier unique context
    ///   Auto  → picks whichever anchor gives shortest total span
    fn genome_biased() -> Vec<u8> {
        //0         1         2         3
        //012345678901234567890123456789012345
        b"ATCGATCGGGGGGGGATCGATCGTTTTTATCGCCCCCCCC".to_vec()
    }

    // ── Left bias tests ───────────────────────────────────────────────────────

    #[test]
    fn left_interval_starts_at_pos() {
        let g = genome_repeats();
        let idx = UniquenessIndex::build(&g);

        for pos in [0, 4, 8, 12] {
            if let Some(iv) = idx.minimal_unique_interval_containing(pos, AnchorBias::Left) {
                assert_eq!(iv.start, pos, "Left: interval must start at pos={pos}");
                assert!(iv.end > pos, "Left: interval must extend past pos");
                assert_unique(&g, &iv);
                assert_minimal_left(&g, &iv);
            }
        }
    }

    #[test]
    fn left_interval_is_unique() {
        let g = genome_simple();
        let idx = UniquenessIndex::build(&g);

        // pos=4 is in the TTTT region — unique and short
        let iv = idx
            .minimal_unique_interval_containing(4, AnchorBias::Left)
            .expect("TTTT region should have a unique interval");
        assert_eq!(iv.start, 4);
        assert_unique(&g, &iv);
        assert_minimal_left(&g, &iv);
    }

    #[test]
    fn left_returns_none_for_n_run() {
        let g = genome_simple();
        let idx = UniquenessIndex::build(&g);
        // N-run starts at pos 20 in genome_simple
        // Ns produce repeated k-mers → no unique interval anchored there
        let result = idx.minimal_unique_interval_containing(22, AnchorBias::Left);
        // We don't assert None strictly (an N-run touching unique flanks might
        // resolve), but if Some, it must be genuinely unique.
        if let Some(iv) = result {
            assert_unique(&g, &iv);
        }
    }

    // ── Right bias tests ──────────────────────────────────────────────────────

    #[test]
    fn right_interval_covers_pos() {
        let g = genome_repeats();
        let idx = UniquenessIndex::build(&g);

        for pos in [5, 10, 15, 20] {
            if pos >= g.len() as GenomeCoord {
                continue;
            }
            if let Some(iv) = idx.minimal_unique_interval_containing(pos, AnchorBias::Right) {
                assert!(iv.start <= pos, "Right: interval must start ≤ pos={pos}");
                assert!(iv.end > pos, "Right: interval must end   > pos={pos}");
                assert_unique(&g, &iv);
                assert_minimal_right(&g, &iv);
            }
        }
    }

    #[test]
    fn right_interval_ends_as_early_as_possible() {
        let g = genome_biased();
        let idx = UniquenessIndex::build(&g);
        let pos = 6;

        let iv_right = idx
            .minimal_unique_interval_containing(pos, AnchorBias::Right)
            .expect("should find a unique interval with Right bias");
        let iv_left = idx
            .minimal_unique_interval_containing(pos, AnchorBias::Left)
            .expect("should find a unique interval with Left bias");

        // Right-biased interval uses a left anchor → end is generally earlier
        // than Left-biased interval's end (which is anchored at pos and must
        // extend further right to become unique).
        // This is not guaranteed in all genomes, but holds for genome_biased.
        assert!(
            iv_right.end <= iv_left.end,
            "Right end={} should be ≤ Left end={} for this genome",
            iv_right.end,
            iv_left.end,
        );
        assert_unique(&g, &iv_right);
    }

    // ── Auto bias tests ───────────────────────────────────────────────────────

    #[test]
    fn auto_is_shortest_among_all_biases() {
        let g = genome_biased();
        let idx = UniquenessIndex::build(&g);

        for pos in 0..g.len() as GenomeCoord {
            let auto = idx.minimal_unique_interval_containing(pos, AnchorBias::Auto);
            let left = idx.minimal_unique_interval_containing(pos, AnchorBias::Left);
            let right = idx.minimal_unique_interval_containing(pos, AnchorBias::Right);

            if let Some(ref a) = auto {
                assert_unique(&g, a);
                assert_minimal_auto(&g, a, pos);

                // Auto must be ≤ Left in length
                if let Some(ref l) = left {
                    assert!(
                        a.len() <= l.len(),
                        "pos={pos}: Auto len={} > Left len={} — Auto must be minimal",
                        a.len(),
                        l.len()
                    );
                }
                // Auto must be ≤ Right in length
                if let Some(ref r) = right {
                    assert!(
                        a.len() <= r.len(),
                        "pos={pos}: Auto len={} > Right len={} — Auto must be minimal",
                        a.len(),
                        r.len()
                    );
                }
            }
        }
    }

    #[test]
    fn auto_covers_pos() {
        let g = genome_repeats();
        let idx = UniquenessIndex::build(&g);

        for pos in 0..g.len() as GenomeCoord {
            if let Some(iv) = idx.minimal_unique_interval_containing(pos, AnchorBias::Auto) {
                assert!(
                    iv.start <= pos && iv.end > pos,
                    "Auto interval [{}, {}) must contain pos={pos}",
                    iv.start,
                    iv.end
                );
                assert_unique(&g, &iv);
            }
        }
    }

    // ── Cross-bias consistency ────────────────────────────────────────────────

    #[test]
    fn all_biases_return_unique_intervals() {
        let g = genome_simple();
        let idx = UniquenessIndex::build(&g);

        for pos in 0..g.len() as GenomeCoord {
            for bias in [AnchorBias::Left, AnchorBias::Right, AnchorBias::Auto] {
                if let Some(iv) = idx.minimal_unique_interval_containing(pos, bias) {
                    assert!(iv.start <= pos, "start must be ≤ pos");
                    assert!(iv.end > pos, "end must be > pos");
                    assert!(iv.end <= g.len() as GenomeCoord, "end must be in bounds");
                    assert_unique(&g, &iv);
                }
            }
        }
    }

    #[test]
    fn left_start_constraint_holds_for_all_positions() {
        let g = genome_repeats();
        let idx = UniquenessIndex::build(&g);

        for pos in 0..g.len() as GenomeCoord {
            if let Some(iv) = idx.minimal_unique_interval_containing(pos, AnchorBias::Left) {
                assert_eq!(
                    iv.start, pos,
                    "Left bias must anchor at pos={pos}, got start={}",
                    iv.start
                );
            }
        }
    }

    // ── Cache integration ─────────────────────────────────────────────────────

    #[test]
    fn cache_preserves_query_results() {
        let g = genome_biased();
        let dir = Path::new("/home/t_steimle/tmp/");
        let key = CacheKey::from_genome("test_contig", &g);

        let contig_end_at = single_contig_end_at(&g);
        let idx_built =
            UniquenessIndex::build_with_cache(&g, &contig_end_at, &contig_end_at, &key, dir)
                .unwrap();
        let idx_cached =
            UniquenessIndex::build_with_cache(&g, &contig_end_at, &contig_end_at, &key, dir)
                .unwrap();

        for pos in 0..g.len() as GenomeCoord {
            for bias in [AnchorBias::Left, AnchorBias::Right, AnchorBias::Auto] {
                assert_eq!(
                    idx_built.minimal_unique_interval_containing(pos, bias),
                    idx_cached.minimal_unique_interval_containing(pos, bias),
                    "cache mismatch at pos={pos} bias={bias:?}",
                );
            }
        }
    }
    use crate::helpers::test_init;

    fn mean(v: &[u32]) -> Option<f64> {
        (!v.is_empty()).then(|| v.iter().sum::<u32>() as f64 / v.len() as f64)
    }

    #[test]
    fn uniqueness_genome() -> anyhow::Result<()> {
        test_init();
        let index = GenomeIndex::build(
            "/home/t_steimle/ref/hs1/chm13v2.0.fa",
            "/home/t_steimle/ref/hs1/hs1_uniqueness",
            None, // None = tout indexer
        )?;

        let mean_read_len = 150u32;
        let interval = 1_000_000u32;

        let reader = bam::Reader::from_path(
            "/mnt/beegfs02/scratch/t_steimle/data/wgs/DUMCO/diag/DUMCO_diag_hs1.bam",
        )
        .unwrap();
        let header = bam::Header::from_template(reader.header());
        let genome_sizes = get_genome_sizes(&header).unwrap();

        let file = File::create("res_150.tsv")?;
        let mut writer = BufWriter::new(file);

        for (contig, size) in genome_sizes {
            let from = 0u32;
            let to = size as u32;

            let half_mean_read_len = mean_read_len / 2;

            let mut current_interval = from;
            loop {
                let mut n_above_half = 0;
                for step in 0..interval {
                    let pos = GenomicPos {
                        contig: contig.clone(),
                        pos: current_interval + step,
                    };

                    if let Some(UniqueContext {
                        before: Some(before),
                        after: Some(after),
                    }) = index.unique_context_at(&pos)
                    {
                        if before > half_mean_read_len || after > half_mean_read_len {
                            n_above_half += 1;
                        }
                        // left.push(before);
                        // right.push(after);
                    }
                }

                println!(
                    "{contig}:{current_interval}-{}\t{n_above_half}",
                    current_interval + interval
                );

                writeln!(
                    writer,
                    "{contig}:{current_interval}-{}\t{n_above_half}",
                    current_interval + interval
                )?;

                current_interval += interval;
                if current_interval >= to {
                    break;
                }
            }
        }

        writer.flush()?;

        Ok(())
    }
}