use of sais suffix array on whole genome

Cargo.lock

@@ -2378,6 +2378,28 @@ dependencies = [
  "libc",
 ]
 
+[[package]]
+name = "libsais"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a739d1c7c4b375c59e25802b577f6d94b0101b77ce478cbcc8a36327957f583c"
+dependencies = [
+ "bytemuck",
+ "either",
+ "libsais-sys",
+ "num-traits",
+]
+
+[[package]]
+name = "libsais-sys"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4ec495a2918bacd868428951b45b9b3fb9464840e2ac559462c4370d6fa2c404"
+dependencies = [
+ "cc",
+ "openmp-sys",
+]
+
 [[package]]
 name = "libsqlite3-sys"
 version = "0.36.0"
@@ -2762,6 +2784,15 @@ version = "0.3.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c08d65885ee38876c4f86fa503fb49d7b507c2b62552df7c70b2fce627e06381"
 
+[[package]]
+name = "openmp-sys"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "caef37114187ab867226e338413fb6cb23ccd785a63d4ee63e7335c88f1079e0"
+dependencies = [
+ "cc",
+]
+
 [[package]]
 name = "openssl-probe"
 version = "0.2.1"
@@ -2947,6 +2978,7 @@ dependencies = [
  "indicatif",
  "itertools 0.14.0",
  "lazy_static",
+ "libsais",
  "log",
  "noodles-bgzf",
  "noodles-core",

Cargo.toml

@@ -61,8 +61,8 @@ chainfile = "0.3.0"
 omics = "0.2.0"
 filetime = "0.2.27"
 pandora_lib_smb = { git = "https://git.t0m4.fr/Thomas/pandora_lib_smb" }
+libsais = { version = "0.2.0", features = ["openmp"] }
 
 [profile.dev]
 opt-level = 0
 debug = false
-

src/uniqueness/cache.rs

@@ -3,9 +3,11 @@ use std::io::{self, BufReader, BufWriter, Read, Write};
 use std::path::Path;
 use std::time::SystemTime;
 
+use super::UniqueLen;
+
 /// Magic + version guard
 const MAGIC: &[u8; 8] = b"UNIQIDX1";
-const VERSION: u32 = 1;
+const VERSION: u32 = 2;
 
 /// Identifies a cache entry unambiguously.
 /// Derive from FASTA path + contig name + genome length + content hash.
@@ -61,10 +63,10 @@ impl CacheKey {
 /// Binary format (little-endian):
 ///
 /// [0..8]   magic   = b"UNIQIDX1"
-/// [8..12]  version = 1u32
+/// [8..12]  version = 2u32
 /// [12..20] n       = genome_len as u64
-/// [20..]   min_unique_len: n × u64  (usize::MAX encoded as u64::MAX)
-pub fn save_to_cache(path: &Path, min_unique_len: &[usize]) -> io::Result<()> {
+/// [20..]   min_unique_len: n × u32  (`NOT_UNIQUE` means no unique interval)
+pub fn save_to_cache(path: &Path, min_unique_len: &[UniqueLen]) -> io::Result<()> {
     let tmp = path.with_extension("uniqidx.tmp");
     {
         let f = File::create(&tmp)?;
@@ -74,13 +76,10 @@ pub fn save_to_cache(path: &Path, min_unique_len: &[usize]) -> io::Result<()> {
         w.write_all(&VERSION.to_le_bytes())?;
         w.write_all(&(min_unique_len.len() as u64).to_le_bytes())?;
 
-        // Bulk-write as u64 array — avoids per-element call overhead
+        // Bulk-write as u32 array — avoids per-element call overhead.
         let buf: Vec<u8> = min_unique_len
             .iter()
-            .flat_map(|&v| {
-                let encoded: u64 = if v == usize::MAX { u64::MAX } else { v as u64 };
-                encoded.to_le_bytes()
-            })
+            .flat_map(|&v| v.to_le_bytes())
             .collect();
         w.write_all(&buf)?;
         w.flush()?;
@@ -90,7 +89,7 @@ pub fn save_to_cache(path: &Path, min_unique_len: &[usize]) -> io::Result<()> {
     Ok(())
 }
 
-pub fn load_from_cache(path: &Path, expected_len: usize) -> io::Result<Vec<usize>> {
+pub fn load_from_cache(path: &Path, expected_len: usize) -> io::Result<Vec<UniqueLen>> {
     let f = File::open(path)?;
     let mut r = BufReader::with_capacity(8 * 1024 * 1024, f);
 
@@ -123,20 +122,13 @@ pub fn load_from_cache(path: &Path, expected_len: usize) -> io::Result<Vec<usize
     }
 
     // Payload: read all at once then cast
-    let byte_len = n * 8;
+    let byte_len = n * 4;
     let mut buf = vec![0u8; byte_len];
     r.read_exact(&mut buf)?;
 
-    let result: Vec<usize> = buf
-        .chunks_exact(8)
-        .map(|c| {
-            let v = u64::from_le_bytes(c.try_into().unwrap());
-            if v == u64::MAX {
-                usize::MAX
-            } else {
-                v as usize
-            }
-        })
+    let result: Vec<UniqueLen> = buf
+        .chunks_exact(4)
+        .map(|c| u32::from_le_bytes(c.try_into().unwrap()))
         .collect();
 
     Ok(result)

src/uniqueness/genome_index.rs

@@ -3,10 +3,13 @@
 use log::info;
 use std::collections::HashMap;
 use std::path::Path;
+use std::time::Instant;
 
 use anyhow::Context;
 
-use super::{AnchorBias, CacheKey, UniqueInterval, UniquenessIndex};
+use super::{
+    AnchorBias, CacheKey, GenomeCoord, UniqueInterval, UniqueLen, UniquenessIndex, NOT_UNIQUE,
+};
 use crate::{
     io::fasta::{read_fai, read_single_contig_fasta, split_fasta, FaiEntry},
     uniqueness::UniqueContext,
@@ -17,14 +20,14 @@ const CONTIG_SEPARATOR: u8 = 0;
 #[derive(Debug, Clone, PartialEq, Eq)]
 pub struct GenomicPos {
     pub contig: String,
-    pub pos: usize, // 0-based
+    pub pos: GenomeCoord, // 0-based
 }
 
 #[derive(Debug, Clone, PartialEq, Eq)]
 pub struct ContigSpan {
     pub name: String,
-    pub start: usize,
-    pub end: usize,
+    pub start: GenomeCoord,
+    pub end: GenomeCoord,
 }
 
 pub struct GenomeIndex {
@@ -80,11 +83,18 @@ impl GenomeIndex {
             .collect();
 
         let mut genome = Vec::new();
+        let mut contig_end_at = Vec::new();
         let mut contigs = Vec::with_capacity(selected.len());
         let mut by_name = HashMap::with_capacity(selected.len());
+        let concat_started = Instant::now();
 
         for entry in selected {
             if !genome.is_empty() {
+                anyhow::ensure!(
+                    genome.len() < u32::MAX as usize,
+                    "concatenated genome exceeds u32::MAX bases"
+                );
+                contig_end_at.push(genome.len() as GenomeCoord);
                 genome.push(CONTIG_SEPARATOR);
             }
 
@@ -104,6 +114,10 @@ impl GenomeIndex {
                 "contig {} contains reserved separator byte",
                 entry.name
             );
+            anyhow::ensure!(
+                genome.len() + seq.len() <= u32::MAX as usize,
+                "concatenated genome exceeds u32::MAX bases"
+            );
 
             info!(
                 "[genome_index] adding contig: {} ({:.1} Mbp)",
@@ -111,9 +125,10 @@ impl GenomeIndex {
                 seq.len() as f64 / 1e6
             );
 
-            let start = genome.len();
+            let start = genome.len() as GenomeCoord;
             genome.extend_from_slice(&seq);
-            let end = genome.len();
+            let end = genome.len() as GenomeCoord;
+            contig_end_at.extend(std::iter::repeat(end).take(seq.len()));
             by_name.insert(entry.name.clone(), contigs.len());
             contigs.push(ContigSpan {
                 name: entry.name.clone(),
@@ -121,10 +136,44 @@ impl GenomeIndex {
                 end,
             });
         }
+        debug_assert_eq!(genome.len(), contig_end_at.len());
+        assert!(
+            genome.len() <= GenomeCoord::MAX as usize,
+            "concatenated genome length exceeds u32::MAX"
+        );
+        info!(
+            "[genome_index] concatenated genome: {} contig(s), {:.1} Mbp including separators ({:.1}s)",
+            contigs.len(),
+            genome.len() as f64 / 1e6,
+            concat_started.elapsed().as_secs_f64()
+        );
+
+        let genome_len = genome.len() as GenomeCoord;
+        info!("[genome_index] building reverse contig-boundary map");
+        let started = Instant::now();
+        let mut rev_contig_end_at: Vec<GenomeCoord> = (0..genome_len).collect();
+        for span in &contigs {
+            let rev_start = genome_len - span.end;
+            let rev_end = genome_len - span.start;
+            for p in rev_start..rev_end {
+                rev_contig_end_at[p as usize] = rev_end;
+            }
+        }
+        info!(
+            "[genome_index] reverse contig-boundary map built in {:.1}s",
+            started.elapsed().as_secs_f64()
+        );
 
         let key = CacheKey::from_genome("genome", &genome);
-        let index = UniquenessIndex::build_with_cache(&genome, &key, cache_dir)
-            .context("build_with_cache failed for concatenated genome")?;
+        info!("[genome_index] building/loading global uniqueness index");
+        let index = UniquenessIndex::build_with_cache(
+            &genome,
+            &contig_end_at,
+            &rev_contig_end_at,
+            &key,
+            cache_dir,
+        )
+        .context("build_with_cache failed for concatenated genome")?;
 
         let _ = std::fs::remove_dir_all(&split_dir);
         Ok(Self {
@@ -138,19 +187,19 @@ impl GenomeIndex {
         self.contigs.iter().map(|span| span.name.as_str())
     }
 
-    pub fn contig_len(&self, contig: &str) -> Option<usize> {
+    pub fn contig_len(&self, contig: &str) -> Option<GenomeCoord> {
         let span = self.contig_span(contig)?;
         Some(span.end - span.start)
     }
 
     pub fn unique_context_at(&self, pos: &GenomicPos) -> Option<UniqueContext> {
         let span = self.contig_span(&pos.contig)?;
-        let global_pos = self.global_pos(pos)?;
+        let global_pos = self.global_boundary_pos(pos)?;
         let context = self.index.unique_context_at(global_pos);
 
         Some(UniqueContext {
-            before: context.before.filter(|&k| global_pos >= span.start + k),
-            after: context.after.filter(|&k| global_pos + k <= span.end),
+            before: context.before.filter(|&k| k <= global_pos - span.start),
+            after: context.after.filter(|&k| k <= span.end - global_pos),
         })
     }
 
@@ -164,7 +213,7 @@ impl GenomeIndex {
         })
     }
 
-    pub fn global_pos(&self, pos: &GenomicPos) -> Option<usize> {
+    pub fn global_pos(&self, pos: &GenomicPos) -> Option<GenomeCoord> {
         let span = self.contig_span(&pos.contig)?;
         if pos.pos < span.end - span.start {
             Some(span.start + pos.pos)
@@ -173,6 +222,15 @@ impl GenomeIndex {
         }
     }
 
+    fn global_boundary_pos(&self, pos: &GenomicPos) -> Option<GenomeCoord> {
+        let span = self.contig_span(&pos.contig)?;
+        if pos.pos <= span.end - span.start {
+            Some(span.start + pos.pos)
+        } else {
+            None
+        }
+    }
+
     pub fn global_interval_to_local(
         &self,
         interval: &UniqueInterval,
@@ -203,7 +261,7 @@ impl GenomeIndex {
 
     fn query_global_in_span(
         &self,
-        global_pos: usize,
+        global_pos: GenomeCoord,
         span: &ContigSpan,
         bias: AnchorBias,
     ) -> Option<UniqueInterval> {
@@ -211,12 +269,12 @@ impl GenomeIndex {
 
         match bias {
             AnchorBias::Left => {
-                let k = self.index.min_unique_len[global_pos];
+                let k = self.index.min_unique_len[global_pos as usize];
                 self.interval_in_span(global_pos, k, span)
             }
             AnchorBias::Right => {
                 for a in (span.start..=global_pos).rev() {
-                    let k = self.index.min_unique_len[a];
+                    let k = self.index.min_unique_len[a as usize];
                     if let Some(iv) = self.interval_in_span(a, k, span) {
                         if iv.end > global_pos {
                             return Some(iv);
@@ -231,7 +289,7 @@ impl GenomeIndex {
                 let start = global_pos.saturating_sub(window).max(span.start);
 
                 for a in start..=global_pos {
-                    let k = self.index.min_unique_len[a];
+                    let k = self.index.min_unique_len[a as usize];
                     let Some(candidate) = self.interval_in_span(a, k, span) else {
                         continue;
                     };
@@ -258,11 +316,11 @@ impl GenomeIndex {
 
     fn interval_in_span(
         &self,
-        start: usize,
-        len: usize,
+        start: GenomeCoord,
+        len: UniqueLen,
         span: &ContigSpan,
     ) -> Option<UniqueInterval> {
-        if len == usize::MAX {
+        if len == NOT_UNIQUE {
             return None;
         }
 
@@ -331,6 +389,103 @@ mod tests {
         let _ = fs::remove_dir_all(&dir);
     }
 
+    #[test]
+    fn sequence_unique_across_all_contigs_is_reported_unique() {
+        let dir = std::env::temp_dir().join(format!(
+            "pandora_uniqueness_unique_test_{}",
+            std::process::id()
+        ));
+        let _ = fs::remove_dir_all(&dir);
+        fs::create_dir_all(&dir).unwrap();
+
+        let fasta = write_fasta_and_fai(&dir, &[("chr1", b"AAAAAC"), ("chr2", b"GGGGTT")]);
+        let index = GenomeIndex::build(&fasta, dir.join("cache"), None).unwrap();
+
+        assert_eq!(
+            index.query(
+                &GenomicPos {
+                    contig: "chr1".to_owned(),
+                    pos: 5,
+                },
+                AnchorBias::Left,
+            ),
+            Some(UniqueInterval { start: 5, end: 6 })
+        );
+
+        let _ = fs::remove_dir_all(&dir);
+    }
+
+    #[test]
+    fn unique_intervals_and_mappings_never_include_separators() {
+        let contigs = vec![
+            ContigSpan {
+                name: "chr1".to_owned(),
+                start: 0,
+                end: 4,
+            },
+            ContigSpan {
+                name: "chr2".to_owned(),
+                start: 5,
+                end: 9,
+            },
+        ];
+        let by_name = HashMap::from([("chr1".to_owned(), 0), ("chr2".to_owned(), 1)]);
+        let index = GenomeIndex {
+            index: UniquenessIndex::build(b"AAAA\0AAAC"),
+            contigs,
+            by_name,
+        };
+
+        assert_eq!(
+            index.global_interval_to_local(&UniqueInterval { start: 4, end: 5 }),
+            None
+        );
+        assert_eq!(
+            index.global_interval_to_local(&UniqueInterval { start: 3, end: 6 }),
+            None
+        );
+    }
+
+    #[test]
+    fn unique_context_respects_contig_boundaries() {
+        let dir = std::env::temp_dir().join(format!(
+            "pandora_uniqueness_context_test_{}",
+            std::process::id()
+        ));
+        let _ = fs::remove_dir_all(&dir);
+        fs::create_dir_all(&dir).unwrap();
+
+        let fasta = write_fasta_and_fai(&dir, &[("chr1", b"AAAAAC"), ("chr2", b"GGGGTT")]);
+        let index = GenomeIndex::build(&fasta, dir.join("cache"), None).unwrap();
+
+        let start = index
+            .unique_context_at(&GenomicPos {
+                contig: "chr1".to_owned(),
+                pos: 0,
+            })
+            .unwrap();
+        assert_eq!(start.before, None);
+
+        let c_base = index
+            .unique_context_at(&GenomicPos {
+                contig: "chr1".to_owned(),
+                pos: 5,
+            })
+            .unwrap();
+        assert_eq!(c_base.after, Some(1));
+
+        let contig_end = index
+            .unique_context_at(&GenomicPos {
+                contig: "chr1".to_owned(),
+                pos: 6,
+            })
+            .unwrap();
+        assert_eq!(contig_end.before, Some(1));
+        assert_eq!(contig_end.after, None);
+
+        let _ = fs::remove_dir_all(&dir);
+    }
+
     #[test]
     fn maps_between_local_and_global_coordinates() {
         let contigs = vec![

src/uniqueness/mod.rs

@@ -1,7 +1,10 @@
-//! Genome uniqueness index based on suffix arrays.
+//! Genome-wide uniqueness index based on suffix arrays.
 //!
 //! This module computes, stores, and queries minimal unique sequence lengths
-//! across a genome.
+//! 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:
@@ -67,64 +70,108 @@
 //! - [`UniqueContext`] — independent left/right uniqueness lengths.
 //! - [`GenomicPos`] — genomic coordinate used by [`GenomeIndex`].
 
-use log::warn;
+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};
 use suffix::{build_lcp_kasai, build_suffix_array};
 
 /// Pure computation: genome bytes → min_unique_len vector.
-/// Separated so GenomeIndex can call it without going through build_with_cache.
-fn build_min_unique_len(genome: &[u8]) -> Vec<usize> {
+///
+/// `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 = build_suffix_array(genome);
     let lcp = build_lcp_kasai(genome, &sa);
 
-    let mut rank = vec![0usize; n];
+    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] = i;
+        rank[s as usize] = i as u32;
     }
-
-    (0..n)
+    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];
+            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 k = left.max(right) + 1;
-            if p + k <= n {
+            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 {
-                usize::MAX
+                NOT_UNIQUE
             }
         })
-        .collect()
+        .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
 }
 
 pub struct UniqueContext {
     /// Bases required before `pos`.
-    pub before: Option<usize>,
+    pub before: Option<UniqueLen>,
 
     /// Bases required from `pos` onward.
-    pub after: Option<usize>,
+    pub after: Option<UniqueLen>,
 }
 
 #[derive(Debug, Clone, PartialEq, Eq)]
 pub struct UniqueInterval {
-    pub start: usize,
-    pub end: usize,
+    pub start: GenomeCoord,
+    pub end: GenomeCoord,
 }
 
 impl UniqueInterval {
-    pub fn len(&self) -> usize {
+    pub fn len(&self) -> UniqueLen {
         self.end - self.start
     }
 }
@@ -137,67 +184,107 @@ pub enum AnchorBias {
 }
 
 pub struct UniquenessIndex {
-    min_unique_len: Vec<usize>,
-    rev_min_unique_len: Vec<usize>,
-    genome_len: usize,
+    min_unique_len: Vec<UniqueLen>,
+    rev_min_unique_len: Vec<UniqueLen>,
+    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], key: &CacheKey, cache_dir: &Path) -> io::Result<Self> {
+    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, rev_min_unique_len) = rayon::join(
-            || -> io::Result<Vec<usize>> {
-                match load_from_cache(&fwd_path, genome.len()) {
-                    Ok(mul) => Ok(mul),
-                    Err(e) => {
-                        if fwd_path.exists() {
-                            warn!("[cache] invalid forward ({e}), rebuilding...");
-                        }
-
-                        let mul = build_min_unique_len(genome);
-                        save_to_cache(&fwd_path, &mul)?;
-                        Ok(mul)
-                    }
+        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()
+                    );
                 }
-            },
-            || -> io::Result<Vec<usize>> {
-                match load_from_cache(&rev_path, genome.len()) {
-                    Ok(mul) => Ok(mul),
-                    Err(e) => {
-                        if rev_path.exists() {
-                            warn!("[cache] invalid reverse ({e}), rebuilding...");
-                        }
 
-                        let mut rev = genome.to_vec();
-                        rev.reverse();
+                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 mul = build_min_unique_len(&rev);
-                        save_to_cache(&rev_path, &mul)?;
-                        Ok(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()
+                    );
                 }
-            },
-        );
 
-        let min_unique_len = min_unique_len?;
-        let rev_min_unique_len = rev_min_unique_len?;
+                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<usize>, rev_min_unique_len: Vec<usize>) -> Self {
+    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();
+        let genome_len = min_unique_len.len() as u32;
 
         Self {
             min_unique_len,
@@ -208,14 +295,14 @@ impl UniquenessIndex {
 
     /// Build without cache (e.g. for tests or one-shot use).
     pub fn build(genome: &[u8]) -> Self {
-        let (min_unique_len, rev_min_unique_len) = rayon::join(
-            || build_min_unique_len(genome),
-            || {
-                let mut rev = genome.to_vec();
-                rev.reverse();
-                build_min_unique_len(&rev)
-            },
-        );
+        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)
     }
@@ -237,15 +324,16 @@ impl UniquenessIndex {
     /// Panics if `pos >= self.genome_len()`.
     pub fn minimal_unique_interval_containing(
         &self,
-        pos: usize,
+        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[pos];
-                if k == usize::MAX {
+                let k = self.min_unique_len[pos_usize];
+                if k == NOT_UNIQUE {
                     return None;
                 }
                 Some(UniqueInterval {
@@ -256,8 +344,8 @@ impl UniquenessIndex {
 
             AnchorBias::Right => {
                 for a in (0..=pos).rev() {
-                    let k = self.min_unique_len[a];
-                    if k == usize::MAX {
+                    let k = self.min_unique_len[a as usize];
+                    if k == NOT_UNIQUE {
                         continue;
                     }
                     if a + k > pos {
@@ -278,8 +366,11 @@ impl UniquenessIndex {
                 let window = 500;
                 let mut best: Option<UniqueInterval> = None;
                 for a in pos.saturating_sub(window)..=pos {
-                    let k = self.min_unique_len[a];
-                    if k == usize::MAX || a + k <= pos {
+                    let k = self.min_unique_len[a as usize];
+                    if k == NOT_UNIQUE {
+                        continue;
+                    }
+                    if a + k <= pos {
                         continue;
                     }
                     let c = UniqueInterval {
@@ -302,7 +393,7 @@ impl UniquenessIndex {
         }
     }
 
-    pub fn genome_len(&self) -> usize {
+    pub fn genome_len(&self) -> GenomeCoord {
         self.genome_len
     }
 
@@ -321,13 +412,17 @@ impl UniquenessIndex {
     ///
     /// # Panics
     ///
-    /// Panics if `pos >= self.genome_len()`.
-    pub fn unique_context_at(&self, pos: usize) -> UniqueContext {
-        assert!(pos < self.genome_len);
+    /// Panics if `pos > self.genome_len()`.
+    pub fn unique_context_at(&self, pos: GenomeCoord) -> UniqueContext {
+        assert!(pos <= self.genome_len);
 
-        let after = match self.min_unique_len[pos] {
-            usize::MAX => None,
-            k => Some(k),
+        let after = if pos == self.genome_len {
+            None
+        } else {
+            match self.min_unique_len[pos as usize] {
+                NOT_UNIQUE => None,
+                k => Some(k),
+            }
         };
 
         let before = if pos == 0 {
@@ -335,8 +430,8 @@ impl UniquenessIndex {
         } else {
             let rev_pos = self.genome_len - pos;
 
-            match self.rev_min_unique_len[rev_pos] {
-                usize::MAX => None,
+            match self.rev_min_unique_len[rev_pos as usize] {
+                NOT_UNIQUE => None,
                 k => Some(k),
             }
         };
@@ -358,6 +453,10 @@ mod tests {
         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();
@@ -365,9 +464,14 @@ mod tests {
 
         let key = CacheKey::from_genome("chr_test", genome);
 
-        let idx1 = UniquenessIndex::build_with_cache(genome, &key, dir).unwrap();
+        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, &key, dir).unwrap();
+        let idx2 =
+            UniquenessIndex::build_with_cache(genome, &contig_end_at, &contig_end_at, &key, dir)
+                .unwrap();
 
         assert_eq!(idx1.min_unique_len, idx2.min_unique_len);
     }
@@ -379,19 +483,38 @@ mod tests {
         let dir = std::path::Path::new("/home/t_steimle/tmp/");
         let key = CacheKey::from_genome("chr_test", genome1);
 
-        UniquenessIndex::build_with_cache(genome1, &key, dir).unwrap();
+        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..iv.end];
+        let needle = &genome[iv.start as usize..iv.end as usize];
         let count = genome
             .windows(needle.len())
             .filter(|w| *w == needle)
@@ -416,7 +539,7 @@ mod tests {
             start: iv.start,
             end: iv.end - 1,
         };
-        let needle = &genome[shorter.start..shorter.end];
+        let needle = &genome[shorter.start as usize..shorter.end as usize];
         let count = genome
             .windows(needle.len())
             .filter(|w| *w == needle)
@@ -439,7 +562,7 @@ mod tests {
             start: iv.start + 1,
             end: iv.end,
         };
-        let needle = &genome[shorter.start..shorter.end];
+        let needle = &genome[shorter.start as usize..shorter.end as usize];
         let count = genome
             .windows(needle.len())
             .filter(|w| *w == needle)
@@ -455,16 +578,16 @@ mod tests {
 
     /// 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: usize) {
+    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() {
+            if a + 1 > genome.len() as GenomeCoord {
                 break;
             }
-            let max_end = (a + iv.len()).min(genome.len());
+            let max_end = (a + iv.len()).min(genome.len() as GenomeCoord);
             for end in (a + 1)..max_end {
-                let needle = &genome[a..end];
+                let needle = &genome[a as usize..end as usize];
                 let count = genome
                     .windows(needle.len())
                     .filter(|w| *w == needle)
@@ -582,7 +705,7 @@ mod tests {
         let idx = UniquenessIndex::build(&g);
 
         for pos in [5, 10, 15, 20] {
-            if pos >= g.len() {
+            if pos >= g.len() as GenomeCoord {
                 continue;
             }
             if let Some(iv) = idx.minimal_unique_interval_containing(pos, AnchorBias::Right) {
@@ -627,7 +750,7 @@ mod tests {
         let g = genome_biased();
         let idx = UniquenessIndex::build(&g);
 
-        for pos in 0..g.len() {
+        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);
@@ -663,7 +786,7 @@ mod tests {
         let g = genome_repeats();
         let idx = UniquenessIndex::build(&g);
 
-        for pos in 0..g.len() {
+        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,
@@ -683,12 +806,12 @@ mod tests {
         let g = genome_simple();
         let idx = UniquenessIndex::build(&g);
 
-        for pos in 0..g.len() {
+        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(), "end must be in bounds");
+                    assert!(iv.end <= g.len() as GenomeCoord, "end must be in bounds");
                     assert_unique(&g, &iv);
                 }
             }
@@ -700,7 +823,7 @@ mod tests {
         let g = genome_repeats();
         let idx = UniquenessIndex::build(&g);
 
-        for pos in 0..g.len() {
+        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,
@@ -719,10 +842,15 @@ mod tests {
         let dir = Path::new("/home/t_steimle/tmp/");
         let key = CacheKey::from_genome("test_contig", &g);
 
-        let idx_built = UniquenessIndex::build_with_cache(&g, &key, dir).unwrap();
-        let idx_cached = UniquenessIndex::build_with_cache(&g, &key, dir).unwrap();
+        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() {
+        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),
@@ -740,7 +868,7 @@ mod tests {
         let index = GenomeIndex::build(
             "/home/t_steimle/ref/hs1/chm13v2.0.fa",
             "/home/t_steimle/ref/hs1/hs1_uniqueness",
-            Some(&["chr1"]), // None = tout indexer
+            None, // None = tout indexer
         )?;
 
         let pos = GenomicPos {

src/uniqueness/suffix.rs

@@ -1,75 +1,78 @@
-use rayon::prelude::*;
+use super::{GenomeCoord, UniqueLen};
+use libsais::{SuffixArrayConstruction, ThreadCount};
+use log::info;
+use std::time::Instant;
 
-pub fn build_suffix_array(genome: &[u8]) -> Vec<usize> {
+pub fn build_suffix_array(genome: &[u8]) -> Vec<GenomeCoord> {
     let n = genome.len();
+    assert!(
+        n <= GenomeCoord::MAX as usize,
+        "suffix array only supports sequences up to u32::MAX bases"
+    );
     if n == 0 {
         return Vec::new();
     }
 
-    let mut sa: Vec<usize> = (0..n).collect();
-    let mut rank: Vec<i64> = genome.iter().map(|&b| b as i64).collect();
-    let mut tmp = vec![0i64; n];
-    let mut gap = 1usize;
+    let thread_count = rayon::current_num_threads().min(u16::MAX as usize) as u16;
+    info!(
+        "[uniqueness] building suffix array with libsais ({:.1} Mbp, {thread_count} thread(s))",
+        n as f64 / 1e6
+    );
+    let started = Instant::now();
 
-    while gap < n {
-        // ── Sort — main bottleneck, worth parallelising ───────────────────
-        sa.par_sort_unstable_by(|&a, &b| {
-            let second = |i: usize| if i + gap < n { rank[i + gap] } else { -1 };
-            (rank[a], second(a)).cmp(&(rank[b], second(b)))
-        });
+    let sa64: Vec<i64> = SuffixArrayConstruction::for_text(genome)
+        .in_owned_buffer64()
+        .multi_threaded(ThreadCount::fixed(thread_count))
+        .run()
+        .expect("libsais failed to construct suffix array")
+        .into_vec();
+    info!(
+        "[uniqueness] suffix array built in {:.1}s",
+        started.elapsed().as_secs_f64()
+    );
 
-        // ── Compute key-change flags in parallel ──────────────────────────
-        let differs: Vec<bool> = (1..n)
-            .into_par_iter()
-            .map(|i| {
-                let second = |i: usize| if i + gap < n { rank[i + gap] } else { -1i64 };
-                let prev = sa[i - 1];
-                let cur = sa[i];
-                (rank[prev], second(prev)) != (rank[cur], second(cur))
-            })
-            .collect();
-
-        // ── Prefix sum — sequential, O(n), cache-friendly ────────────────
-        tmp[sa[0]] = 0;
-        for i in 1..n {
-            tmp[sa[i]] = tmp[sa[i - 1]] + if differs[i - 1] { 1 } else { 0 };
-        }
-
-        rank.copy_from_slice(&tmp);
-
-        if rank[sa[n - 1]] == (n - 1) as i64 {
-            break;
-        }
-        gap *= 2;
-    }
-
-    sa
+    sa64.into_iter()
+        .map(|pos| {
+            assert!(pos >= 0, "libsais produced negative suffix position");
+            assert!(
+                pos <= GenomeCoord::MAX as i64,
+                "libsais produced suffix position outside u32 range"
+            );
+            pos as GenomeCoord
+        })
+        .collect()
 }
 
-pub fn build_lcp_kasai(genome: &[u8], sa: &[usize]) -> Vec<usize> {
+pub fn build_lcp_kasai(genome: &[u8], sa: &[GenomeCoord]) -> Vec<UniqueLen> {
     let n = genome.len();
-    let mut rank = vec![0usize; n];
+    info!("[uniqueness] building LCP array");
+    let started = Instant::now();
+    let mut rank = vec![0; n];
     for (i, &s) in sa.iter().enumerate() {
-        rank[s] = i;
+        rank[s as usize] = i as GenomeCoord;
     }
 
-    let mut lcp = vec![0usize; n];
+    let mut lcp = vec![0; n];
     let mut h = 0usize;
 
     for p in 0..n {
-        let r = rank[p];
+        let r = rank[p] as usize;
         if r == 0 {
             h = 0;
             continue;
         }
-        let prev = sa[r - 1];
+        let prev = sa[r - 1] as usize;
         while p + h < n && prev + h < n && genome[p + h] == genome[prev + h] {
             h += 1;
         }
-        lcp[r] = h;
+        lcp[r] = h as UniqueLen;
         if h > 0 {
             h -= 1;
         }
     }
+    info!(
+        "[uniqueness] LCP array built in {:.1}s",
+        started.elapsed().as_secs_f64()
+    );
     lcp
 }