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

@@ -1,33 +1,48 @@
 // src/uniqueness/genome_index.rs
 
-use log::{info, warn};
+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::{FaiEntry, read_fai, read_single_contig_fasta, split_fasta},
-    uniqueness::{UniqueContext, cache::load_from_cache},
+    io::fasta::{read_fai, read_single_contig_fasta, split_fasta, FaiEntry},
+    uniqueness::UniqueContext,
 }; // ton module existant
 
+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: GenomeCoord,
+    pub end: GenomeCoord,
 }
 
 pub struct GenomeIndex {
-    contigs: HashMap<String, UniquenessIndex>,
+    index: UniquenessIndex,
+    contigs: Vec<ContigSpan>,
+    by_name: HashMap<String, usize>,
 }
 
 impl GenomeIndex {
     /// Build or load from cache.
     ///
     /// Uses your existing `split_fasta` to stream contig by contig into
-    /// a temp dir, builds/loads a UniquenessIndex per contig, then drops
-    /// the split files.
+    /// a temp dir, concatenates selected contigs with non-FASTA separator
+    /// bytes, builds/loads one global [`UniquenessIndex`], then drops the
+    /// split files.
     ///
     /// `filter`: if Some, only process these contig names.
     pub fn build(
@@ -43,103 +58,469 @@ impl GenomeIndex {
         let fai = read_fai(&fai_path)
             .with_context(|| format!("cannot read FAI {}", fai_path.display()))?;
 
-        // ── Premier passage : quels contigs ont besoin d'être (re)buildés ? ──
-        let mut need_build: Vec<&FaiEntry> = Vec::new();
-        let mut cached: HashMap<String, UniquenessIndex> = HashMap::new();
-
-        for entry in &fai {
-            if let Some(allowed) = filter {
-                if !allowed.contains(&entry.name.as_str()) {
-                    continue;
-                }
-            }
-
-            let key = CacheKey::from_fai(&entry.name, entry.len, fasta_path)?;
-            let fwd_path = cache_dir.join(key.file_name());
-            let rev_path = cache_dir.join(format!("rev.{}", key.file_name()));
-
-            if fwd_path.exists() && rev_path.exists() {
-                match (
-                    load_from_cache(&fwd_path, entry.len),
-                    load_from_cache(&rev_path, entry.len),
-                ) {
-                    (Ok(min_unique_len), Ok(rev_min_unique_len)) => {
-                        info!("[genome_index] cache hit  : {}", entry.name);
-                        cached.insert(
-                            entry.name.clone(),
-                            UniquenessIndex::from_raw(min_unique_len, rev_min_unique_len),
-                        );
-                    }
-                    (fwd, rev) => {
-                        warn!(
-                            "[genome_index] cache miss : {} (fwd={:?}, rev={:?})",
-                            entry.name,
-                            fwd.err(),
-                            rev.err()
-                        );
-                        need_build.push(entry);
-                    }
-                }
-            } else {
-                warn!("[genome_index] cache miss : {}", entry.name);
-                need_build.push(entry);
-            }
-        }
-
-        if need_build.is_empty() {
-            info!("[genome_index] all contigs loaded from cache, FASTA not read.");
-            return Ok(Self { contigs: cached });
-        }
+        let selected: Vec<&FaiEntry> = fai
+            .iter()
+            .filter(|entry| {
+                filter
+                    .map(|allowed| allowed.contains(&entry.name.as_str()))
+                    .unwrap_or(true)
+            })
+            .collect();
+        anyhow::ensure!(!selected.is_empty(), "no FASTA contigs selected");
 
-        let need_names: Vec<&str> = need_build.iter().map(|e| e.name.as_str()).collect();
+        let selected_names: Vec<&str> = selected.iter().map(|e| e.name.as_str()).collect();
         info!(
-            "[genome_index] building {} contig(s): {:?}",
-            need_names.len(),
-            need_names
+            "[genome_index] building global uniqueness over {} contig(s): {:?}",
+            selected_names.len(),
+            selected_names
         );
 
         let split_dir = cache_dir.join("_split_tmp");
         let contig_fastas = split_fasta(fasta_path, &split_dir)?;
+        let contig_fastas_by_name: HashMap<&str, &Path> = contig_fastas
+            .iter()
+            .map(|cf| (cf.name.as_str(), cf.fasta_path.as_path()))
+            .collect();
 
-        for cf in &contig_fastas {
-            if !need_names.contains(&cf.name.as_str()) {
-                continue;
+        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);
             }
 
-            let entry = need_build.iter().find(|e| e.name == cf.name).unwrap();
-            let seq = read_single_contig_fasta(&cf.fasta_path)?;
+            let fasta_path = contig_fastas_by_name
+                .get(entry.name.as_str())
+                .with_context(|| format!("split FASTA missing selected contig {}", entry.name))?;
+            let seq = read_single_contig_fasta(fasta_path)?;
+            anyhow::ensure!(
+                seq.len() == entry.len,
+                "FASTA/FAI length mismatch for {}: sequence has {}, FAI has {}",
+                entry.name,
+                seq.len(),
+                entry.len
+            );
+            anyhow::ensure!(
+                !seq.contains(&CONTIG_SEPARATOR),
+                "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] building SA : {} ({:.1} Mbp)",
-                cf.name,
+                "[genome_index] adding contig: {} ({:.1} Mbp)",
+                entry.name,
                 seq.len() as f64 / 1e6
             );
 
-            let key = CacheKey::from_fai(&cf.name, entry.len, fasta_path)?;
-            let idx = UniquenessIndex::build_with_cache(&seq, &key, cache_dir)
-                .with_context(|| format!("build_with_cache failed for {}", cf.name))?;
-            cached.insert(cf.name.clone(), idx);
+            let start = genome.len() as GenomeCoord;
+            genome.extend_from_slice(&seq);
+            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(),
+                start,
+                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);
+        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 { contigs: cached })
+        Ok(Self {
+            index,
+            contigs,
+            by_name,
+        })
     }
 
     pub fn contig_names(&self) -> impl Iterator<Item = &str> {
-        self.contigs.keys().map(String::as_str)
+        self.contigs.iter().map(|span| span.name.as_str())
     }
 
-    pub fn contig_len(&self, contig: &str) -> Option<usize> {
-        self.contigs.get(contig).map(|idx| idx.genome_len())
+    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> {
-        Some(self.contigs.get(&pos.contig)?.unique_context_at(pos.pos))
+        let span = self.contig_span(&pos.contig)?;
+        let global_pos = self.global_boundary_pos(pos)?;
+        let context = self.index.unique_context_at(global_pos);
+
+        Some(UniqueContext {
+            before: context.before.filter(|&k| k <= global_pos - span.start),
+            after: context.after.filter(|&k| k <= span.end - global_pos),
+        })
     }
 
     pub fn query(&self, pos: &GenomicPos, bias: AnchorBias) -> Option<UniqueInterval> {
+        let span = self.contig_span(&pos.contig)?;
+        let global_pos = self.global_pos(pos)?;
+        let global = self.query_global_in_span(global_pos, span, bias)?;
+        Some(UniqueInterval {
+            start: global.start - span.start,
+            end: global.end - span.start,
+        })
+    }
+
+    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)
+        } else {
+            None
+        }
+    }
+
+    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,
+    ) -> Option<(String, UniqueInterval)> {
+        let span = self.contig_span_for_global_interval(interval)?;
+        Some((
+            span.name.clone(),
+            UniqueInterval {
+                start: interval.start - span.start,
+                end: interval.end - span.start,
+            },
+        ))
+    }
+
+    fn contig_span(&self, contig: &str) -> Option<&ContigSpan> {
+        self.by_name.get(contig).map(|&idx| &self.contigs[idx])
+    }
+
+    fn contig_span_for_global_interval(&self, interval: &UniqueInterval) -> Option<&ContigSpan> {
+        if interval.start >= interval.end {
+            return None;
+        }
+
         self.contigs
-            .get(&pos.contig)?
-            .minimal_unique_interval_containing(pos.pos, bias)
+            .iter()
+            .find(|span| span.start <= interval.start && interval.end <= span.end)
+    }
+
+    fn query_global_in_span(
+        &self,
+        global_pos: GenomeCoord,
+        span: &ContigSpan,
+        bias: AnchorBias,
+    ) -> Option<UniqueInterval> {
+        debug_assert!(span.start <= global_pos && global_pos < span.end);
+
+        match bias {
+            AnchorBias::Left => {
+                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 as usize];
+                    if let Some(iv) = self.interval_in_span(a, k, span) {
+                        if iv.end > global_pos {
+                            return Some(iv);
+                        }
+                    }
+                }
+                None
+            }
+            AnchorBias::Auto => {
+                let window = 500;
+                let mut best: Option<UniqueInterval> = None;
+                let start = global_pos.saturating_sub(window).max(span.start);
+
+                for a in start..=global_pos {
+                    let k = self.index.min_unique_len[a as usize];
+                    let Some(candidate) = self.interval_in_span(a, k, span) else {
+                        continue;
+                    };
+                    if candidate.end <= global_pos {
+                        continue;
+                    }
+
+                    best = Some(match best {
+                        None => candidate,
+                        Some(current) => {
+                            if candidate.len() < current.len() {
+                                candidate
+                            } else {
+                                current
+                            }
+                        }
+                    });
+                }
+
+                best
+            }
+        }
+    }
+
+    fn interval_in_span(
+        &self,
+        start: GenomeCoord,
+        len: UniqueLen,
+        span: &ContigSpan,
+    ) -> Option<UniqueInterval> {
+        if len == NOT_UNIQUE {
+            return None;
+        }
+
+        let end = start.checked_add(len)?;
+        if span.start <= start && start < span.end && end <= span.end {
+            Some(UniqueInterval { start, end })
+        } else {
+            None
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use std::{collections::HashMap, fs, path::Path};
+
+    use super::*;
+
+    fn write_fasta_and_fai(dir: &Path, records: &[(&str, &[u8])]) -> std::path::PathBuf {
+        let fasta = dir.join("test.fa");
+        let mut fasta_body = String::new();
+        let mut fai_body = String::new();
+
+        for (name, seq) in records {
+            fasta_body.push('>');
+            fasta_body.push_str(name);
+            fasta_body.push('\n');
+            fasta_body.push_str(std::str::from_utf8(seq).unwrap());
+            fasta_body.push('\n');
+
+            fai_body.push_str(name);
+            fai_body.push('\t');
+            fai_body.push_str(&seq.len().to_string());
+            fai_body.push_str("\t0\t0\t0\n");
+        }
+
+        fs::write(&fasta, fasta_body).unwrap();
+        fs::write(dir.join("test.fa.fai"), fai_body).unwrap();
+        fasta
+    }
+
+    #[test]
+    fn repeated_sequence_on_another_contig_is_not_unique() {
+        let dir =
+            std::env::temp_dir().join(format!("pandora_uniqueness_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"ACGT"), ("chr2", b"ACGT")]);
+        let index = GenomeIndex::build(&fasta, dir.join("cache"), None).unwrap();
+
+        for bias in [AnchorBias::Left, AnchorBias::Right, AnchorBias::Auto] {
+            assert_eq!(
+                index.query(
+                    &GenomicPos {
+                        contig: "chr1".to_owned(),
+                        pos: 0,
+                    },
+                    bias,
+                ),
+                None,
+                "{bias:?} should not use a separator-spanning interval to make chr1:0 unique"
+            );
+        }
+
+        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![
+            ContigSpan {
+                name: "chr1".to_owned(),
+                start: 0,
+                end: 3,
+            },
+            ContigSpan {
+                name: "chr2".to_owned(),
+                start: 4,
+                end: 6,
+            },
+        ];
+        let by_name = HashMap::from([("chr1".to_owned(), 0), ("chr2".to_owned(), 1)]);
+        let index = GenomeIndex {
+            index: UniquenessIndex::build(b"AAA\0TT"),
+            contigs,
+            by_name,
+        };
+
+        assert_eq!(
+            index.global_pos(&GenomicPos {
+                contig: "chr2".to_owned(),
+                pos: 1,
+            }),
+            Some(5)
+        );
+        assert_eq!(
+            index.global_interval_to_local(&UniqueInterval { start: 4, end: 6 }),
+            Some(("chr2".to_owned(), UniqueInterval { start: 0, end: 2 }))
+        );
+        assert_eq!(
+            index.global_interval_to_local(&UniqueInterval { start: 2, end: 5 }),
+            None
+        );
     }
 }

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,58 +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;
-    }
-
-    let mut mul = vec![usize::MAX; n];
-    for p in 0..n {
-        let r = rank[p];
-        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;
-        mul[p] = if p + k <= n { k } else { usize::MAX };
+        rank[s as usize] = i as u32;
     }
-    mul
+    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
 }
 
 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
     }
 }
@@ -131,45 +184,89 @@ 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 = match load_from_cache(&fwd_path, genome.len()) {
-            Ok(mul) => mul,
+            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);
+                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) => mul,
+            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);
+                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
             }
@@ -180,10 +277,14 @@ impl UniquenessIndex {
 
     /// 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,
@@ -194,11 +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 = build_min_unique_len(genome);
+        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_min_unique_len = build_min_unique_len(&rev);
+        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)
     }
@@ -220,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 {
@@ -239,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 {
@@ -261,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 {
@@ -285,7 +393,7 @@ impl UniquenessIndex {
         }
     }
 
-    pub fn genome_len(&self) -> usize {
+    pub fn genome_len(&self) -> GenomeCoord {
         self.genome_len
     }
 
@@ -304,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 {
@@ -318,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),
             }
         };
@@ -341,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();
@@ -348,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);
     }
@@ -362,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)
@@ -399,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)
@@ -422,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)
@@ -438,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)
@@ -565,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) {
@@ -610,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);
@@ -646,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,
@@ -666,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);
                 }
             }
@@ -683,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,
@@ -702,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),
@@ -723,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
 }