rustalign_fmindex/

ebwt.rs

//! Extended Burrows-Wheeler Transform (Ebwt) - Main FM-index structure
//!
//! This is the core data structure for RustAlign's sequence alignment.
//! It implements a memory-mapped FM-index for efficient pattern matching.

use crate::{EbwtParams, FmIndex};
use memmap2::Mmap;
use rustalign_common::{AlignError, AlignResult, Nuc, NucPair};
use std::fs::File;
use std::path::{Path, PathBuf};

/// Extended Burrows-Wheeler Transform (FM-index)
///
/// This is the main data structure used by RustAlign for ultrafast
/// sequence alignment. It is memory-mapped from .rai index files.
pub struct Ebwt {
    /// Index parameters
    params: EbwtParams,

    /// Path to the index base (without extension)
    path: PathBuf,

    /// Mmapped forward index (.1.rai file)
    ebwt_fw: Option<Mmap>,

    /// Mmapped suffix array samples (.2.rai file)
    offs_mmap: Option<Mmap>,

    /// Mmapped reverse index (not currently used)
    #[allow(dead_code)]
    ebwt_rc: Option<Mmap>,

    /// Whether this is a "large" index
    large: bool,

    /// Number of reference sequences (chromosomes)
    n_pat: u32,

    /// Number of fragments
    n_frag: u32,

    /// Length of each reference sequence
    plen: Vec<u32>,

    /// Reference sequence names (chromosome names)
    refnames: Vec<String>,

    /// rstarts array: [frag_start, text_id, frag_offset] for each fragment
    /// Stored as flat array of n_frag * 3 elements
    rstarts: Vec<u32>,

    /// Offset to ftab in the mmapped data
    ftab_off: usize,

    /// Offset to eftab in the mmapped data
    #[allow(dead_code)]
    eftab_off: usize,

    /// Offset to ebwt data in the mmapped data
    ebwt_off: usize,

    /// Offset to rstarts in the mmapped data
    #[allow(dead_code)]
    rstarts_off: usize,

    /// Position of '$' in BWT
    z_off: u32,

    /// First character array (C array) - start positions for each character
    fchr: [u64; 5],
}

impl Ebwt {
    /// Load an existing index from disk
    ///
    /// # Arguments
    /// * `path_base` - Path to the index (without .1.rai, .2.rai, etc. extensions)
    ///
    /// # Errors
    /// Returns an error if the index files don't exist or are corrupted.
    pub fn load<P: AsRef<Path>>(path_base: P) -> AlignResult<Self> {
        let path_base = path_base.as_ref();
        Self::load_with_options(path_base, false)
    }

    /// Load a "large" index variant
    pub fn load_large<P: AsRef<Path>>(path_base: P) -> AlignResult<Self> {
        let path_base = path_base.as_ref();
        Self::load_with_options(path_base, true)
    }

    /// Load index with size option
    fn load_with_options(path_base: &Path, large: bool) -> AlignResult<Self> {
        // Construct file extension based on size
        let ext = if large { ".l" } else { "" };
        let base_name = path_base
            .file_name()
            .ok_or_else(|| AlignError::InvalidFormat("Invalid index path".into()))?
            .to_str()
            .unwrap();

        // Load the first .1.rai file
        let file_1 = path_base.with_file_name(format!("{}.1.rai{}", base_name, ext));
        let file_2 = path_base.with_file_name(format!("{}.2.rai{}", base_name, ext));

        if !file_1.exists() || !file_2.exists() {
            return Err(AlignError::IndexCorrupted {
                path: path_base.to_path_buf(),
            });
        }

        // Load and parse parameters from file header
        let file = File::open(&file_1)?;
        let mmap = unsafe { Mmap::map(&file)? };

        // Load the suffix array samples from .2.rai file
        let file2 = File::open(&file_2)?;
        let offs_mmap = unsafe { Mmap::map(&file2)? };

        // Parse header and get offsets
        let (
            params,
            n_pat,
            n_frag,
            plen,
            rstarts,
            ftab_off,
            eftab_off,
            ebwt_off,
            rstarts_off,
            z_off,
            fchr,
            refnames,
        ) = Self::parse_header_and_offsets(&mmap)?;

        Ok(Self {
            params,
            path: path_base.to_path_buf(),
            ebwt_fw: Some(mmap),
            offs_mmap: Some(offs_mmap),
            ebwt_rc: None,
            large,
            n_pat,
            n_frag,
            plen,
            refnames,
            rstarts,
            ftab_off,
            eftab_off,
            ebwt_off,
            rstarts_off,
            z_off,
            fchr,
        })
    }

    /// Parse EbwtParams and calculate offsets from mmapped data
    ///
    /// The RustAlign index file format (.rai):
    /// - endian hint: int32 (should be 1 for little-endian)
    /// - len: uint32 (length of reference sequence)
    /// - lineRate: int32
    /// - linesPerSide: int32 (not used, always 2)
    /// - offRate: int32
    /// - ftabChars: int32
    /// - flags: int32 (colorspace and entire-reverse flags)
    /// - nPat: uint32 (number of reference sequences)
    /// - plen: uint32[nPat] (lengths of each sequence)
    /// - nFrag: uint32 (number of fragments)
    /// - rstarts: uint32[nFrag * 3] (fragment starts, text ids, offsets)
    /// - ebwt: uint8[ebwtTotLen] (BWT data with side structure)
    /// - zOff: uint32 (position of '$' in BWT)
    /// - fchr: uint32[5] (first character ranges: A, C, G, T, N)
    /// - ftab: uint32[4^ftabChars + 1] (frequency table)
    /// - eftab: uint32[ftabChars * 2] (extended frequency table)
    /// - refnames: null-terminated, newline-separated chromosome names
    #[allow(clippy::type_complexity)]
    fn parse_header_and_offsets(
        data: &[u8],
    ) -> AlignResult<(
        EbwtParams,
        u32,
        u32,
        Vec<u32>,
        Vec<u32>,
        usize,
        usize,
        usize,
        usize,
        u32,
        [u64; 5],
        Vec<String>,
    )> {
        use byteorder::{LittleEndian, ReadBytesExt};
        use std::io::Cursor;

        if data.len() < 28 {
            return Err(AlignError::InvalidFormat(
                "Index file too small for header".into(),
            ));
        }

        let mut cursor = Cursor::new(data);

        // Read endian hint (should be 1 for little-endian)
        let endian_hint = cursor
            .read_i32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read endian hint: {}", e)))?;

        if endian_hint != 1 {
            return Err(AlignError::InvalidFormat(format!(
                "Invalid endian hint: {} (expected 1)",
                endian_hint
            )));
        }

        // Read header fields
        let len = cursor
            .read_u32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read len: {}", e)))?
            as u64;

        let line_rate = cursor
            .read_i32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read lineRate: {}", e)))?;

        let _lines_per_side = cursor.read_i32::<LittleEndian>().map_err(|e| {
            AlignError::InvalidFormat(format!("Failed to read linesPerSide: {}", e))
        })?;

        let off_rate = cursor
            .read_i32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read offRate: {}", e)))?;

        let ftab_chars = cursor
            .read_i32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read ftabChars: {}", e)))?;

        let flags = cursor
            .read_i32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read flags: {}", e)))?;

        // Extract flags (from EBWT_FLAGS enum)
        // Flags are stored as negative values
        let color = flags < 0 && ((-flags) & 2) != 0;
        let entire_reverse = flags < 0 && ((-flags) & 4) != 0;

        // Create params using the parsed values
        let params =
            EbwtParams::with_options(len, line_rate, off_rate, ftab_chars, color, entire_reverse);

        if !params.is_valid() {
            return Err(AlignError::InvalidFormat(format!(
                "Invalid index parameters: len={}, lineRate={}, offRate={}, ftabChars={}",
                len, line_rate, off_rate, ftab_chars
            )));
        }

        // Read nPat (number of reference sequences)
        let n_pat = cursor
            .read_u32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read nPat: {}", e)))?;

        // Read plen array (n_pat * 4 bytes)
        let plen_off = cursor.position() as usize;
        let plen_size = n_pat as usize * 4;
        if plen_off + plen_size > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at plen".into(),
            ));
        }

        let mut plen = Vec::with_capacity(n_pat as usize);
        for i in 0..n_pat as usize {
            plen.push(Self::read_u32_at(data, plen_off + i * 4)?);
        }
        cursor.set_position(cursor.position() + plen_size as u64);

        // Read nFrag (number of fragments)
        let n_frag = cursor
            .read_u32::<LittleEndian>()
            .map_err(|e| AlignError::InvalidFormat(format!("Failed to read nFrag: {}", e)))?;

        // Read rstarts array (n_frag * 3 * 4 bytes)
        let rstarts_off = cursor.position() as usize;
        let rstarts_size = n_frag as usize * 3 * 4;
        if rstarts_off + rstarts_size > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at rstarts".into(),
            ));
        }

        let mut rstarts = Vec::with_capacity(n_frag as usize * 3);
        for i in 0..(n_frag as usize * 3) {
            rstarts.push(Self::read_u32_at(data, rstarts_off + i * 4)?);
        }
        cursor.set_position(cursor.position() + rstarts_size as u64);

        // Skip ebwt data (the BWT itself - this is the largest section)
        let ebwt_tot_len = params.ebwt_tot_len as usize;
        if cursor.position() as usize + ebwt_tot_len > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at ebwt".into(),
            ));
        }

        // Record ebwt offset before skipping
        let ebwt_off = cursor.position() as usize;
        cursor.set_position(cursor.position() + ebwt_tot_len as u64);

        // Read zOff (single uint32)
        if cursor.position() as usize + 4 > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at zOff".into(),
            ));
        }
        let z_off = Self::read_u32_at(data, cursor.position() as usize)?;
        cursor.set_position(cursor.position() + 4);

        // Now we're at fchr (first character array - 5 uint32 values)
        let fchr_off = cursor.position() as usize;
        if fchr_off + 20 > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at fchr".into(),
            ));
        }

        // Read fchr (first character ranges)
        // fchr[0] = start of A (always 0)
        // fchr[1] = start of C (count of A)
        // fchr[2] = start of G (count of A + C)
        // fchr[3] = start of T (count of A + C + G)
        // fchr[4] = end of T (bwt length)
        let fchr = [
            Self::read_u32_at(data, fchr_off)? as u64,      // A start
            Self::read_u32_at(data, fchr_off + 4)? as u64,  // C start
            Self::read_u32_at(data, fchr_off + 8)? as u64,  // G start
            Self::read_u32_at(data, fchr_off + 12)? as u64, // T start
            Self::read_u32_at(data, fchr_off + 16)? as u64, // N/end
        ];

        // Move past fchr to get to ftab
        cursor.set_position(fchr_off as u64 + 20);

        // ftab starts here
        let ftab_off = cursor.position() as usize;
        let ftab_len = params.ftab_len as usize;
        let ftab_size = ftab_len * 4;

        if ftab_off + ftab_size > data.len() {
            return Err(AlignError::InvalidFormat(
                "Index file truncated at ftab".into(),
            ));
        }

        // eftab follows ftab
        let eftab_off = ftab_off + ftab_size;
        let eftab_size = params.eftab_len as usize * 4;

        // Read chromosome names (refnames) - they come after eftab
        // The names are stored as the full FASTA header, but we only need
        // the first word (the actual chromosome name like "1", "Mt", "Pt")
        let refnames_off = eftab_off + eftab_size;
        let mut refnames = Vec::new();

        if refnames_off < data.len() {
            let mut current_name = String::new();
            let mut pos = refnames_off;

            while pos < data.len() {
                let c = data[pos];
                pos += 1;

                if c == 0 {
                    // Null terminator - end of names
                    if !current_name.is_empty() {
                        // Only take the first word (before any space)
                        let short_name = current_name.split_whitespace().next().unwrap_or("");
                        refnames.push(short_name.to_string());
                    }
                    break;
                } else if c == b'\n' {
                    // Newline - start new name
                    if !current_name.is_empty() {
                        // Only take the first word (before any space)
                        let short_name = current_name.split_whitespace().next().unwrap_or("");
                        refnames.push(short_name.to_string());
                        current_name.clear();
                    }
                } else {
                    current_name.push(c as char);
                }
            }

            // Don't forget the last name if not empty
            if !current_name.is_empty() {
                let short_name = current_name.split_whitespace().next().unwrap_or("");
                refnames.push(short_name.to_string());
            }
        }

        Ok((
            params,
            n_pat,
            n_frag,
            plen,
            rstarts,
            ftab_off,
            eftab_off,
            ebwt_off,
            rstarts_off,
            z_off,
            fchr,
            refnames,
        ))
    }

    /// Read a little-endian u32 at a given offset
    fn read_u32_at(data: &[u8], offset: usize) -> AlignResult<u32> {
        if offset + 4 > data.len() {
            return Err(AlignError::InvalidFormat("Offset out of bounds".into()));
        }
        let bytes: [u8; 4] = data[offset..offset + 4].try_into().unwrap();
        Ok(u32::from_le_bytes(bytes))
    }

    /// Read a u32 from ftab at given index
    fn ftab(&self, idx: usize) -> u64 {
        let offset = self.ftab_off + idx * 4;
        u32::from_le_bytes(
            self.ebwt_fw.as_ref().unwrap()[offset..offset + 4]
                .try_into()
                .unwrap(),
        ) as u64
    }

    /// Get the fchr (first character array) entry for a nucleotide
    fn fchr(&self, c: u8) -> u64 {
        match c {
            0 => self.fchr[0], // A
            1 => self.fchr[1], // C
            2 => self.fchr[2], // G
            3 => self.fchr[3], // T
            _ => self.fchr[4], // N/others
        }
    }

    /// Internal rank operation - count occurrences of character c up to position pos
    ///
    /// Uses RustAlign's side-based counting structure for efficiency.
    /// Each side has the format: [bwt_data (side_bwt_sz bytes)] [4 x u32 counts for A,C,G,T]
    /// The counts are cumulative - they represent all occurrences up to that side boundary.
    fn rank_internal(&self, c: usize, pos: u64) -> AlignResult<usize> {
        if pos == 0 {
            return Ok(0);
        }

        let data = self
            .ebwt_fw
            .as_ref()
            .ok_or_else(|| AlignError::InvalidFormat("Index not loaded".into()))?;

        // BWT is bit-packed: 2 bits per nucleotide, 4 nucleotides per byte
        let side_bwt_len = self.params.side_bwt_len as u64;
        let side_bwt_sz = self.params.side_bwt_sz as usize;
        let side_sz = self.params.side_sz as usize;

        // Determine which side the position falls in
        let side_num = pos / side_bwt_len;
        let char_off = (pos % side_bwt_len) as usize; // character offset within side

        // Read the cumulative count from the side boundary
        // Counts are stored at offset side_bwt_sz within each side
        let mut count = 0usize;

        if side_num > 0 {
            // Read cumulative counts from the start of the current side
            // These counts represent all occurrences in sides 0..side_num
            let side_start = self.ebwt_off + side_num as usize * side_sz;
            let count_offset = side_start + side_bwt_sz + c * 4; // 4 bytes per count

            if count_offset + 4 <= data.len() {
                count = u32::from_le_bytes(data[count_offset..count_offset + 4].try_into().unwrap())
                    as usize;
            }
        }

        // Count occurrences within the current side up to char_off
        // This is the countUpTo part
        if char_off > 0 {
            let side_start = self.ebwt_off + side_num as usize * side_sz;

            // Count nucleotides one by one from start of side to char_off
            for i in 0..char_off {
                let byte_off = side_start + (i >> 2); // i / 4
                let bit_off = (i & 3) * 2; // (i % 4) * 2

                if byte_off < data.len() {
                    let byte = data[byte_off];
                    let nuc = ((byte >> bit_off) & 0x03) as usize;
                    if nuc == c {
                        count += 1;
                    }
                }
            }
        }

        Ok(count)
    }

    /// Get the index parameters
    pub fn params(&self) -> &EbwtParams {
        &self.params
    }

    /// Get the path to this index
    pub fn path(&self) -> &Path {
        &self.path
    }

    /// Check if this is a large index
    pub fn is_large(&self) -> bool {
        self.large
    }

    /// Get the chromosome names
    pub fn refnames(&self) -> &[String] {
        &self.refnames
    }

    /// Get the number of chromosomes
    pub fn n_pat(&self) -> u32 {
        self.n_pat
    }

    /// Get the number of fragments
    pub fn n_frag(&self) -> u32 {
        self.n_frag
    }

    /// Get the rstarts array
    pub fn rstarts(&self) -> &[u32] {
        &self.rstarts
    }

    /// Get chromosome lengths
    pub fn plen(&self) -> &[u32] {
        &self.plen
    }

    /// Get the zOff value (position of $ in BWT)
    pub fn z_off(&self) -> u32 {
        self.z_off
    }

    /// Read a suffix array sample from the .2.rai file
    ///
    /// The offs array stores sampled suffix array values at positions
    /// that are multiples of (1 << offRate).
    #[allow(clippy::collapsible_if)]
    fn offs(&self, idx: usize) -> u32 {
        // Skip the 4-byte endianness hint at the start of .2.rai file
        let offset = 4 + idx * 4;
        if let Some(ref mmap) = self.offs_mmap {
            if offset + 4 <= mmap.len() {
                return u32::from_le_bytes(mmap[offset..offset + 4].try_into().unwrap());
            }
        }
        0
    }

    /// Convert a BWT row to a genome offset using suffix array samples and LF-mapping.
    ///
    /// This implements RustAlign's getOffset() algorithm:
    /// 1. If row is zOff, return 0
    /// 2. If row is at a sampled position, return the sampled offset
    /// 3. Otherwise, walk back using LF-mapping until hitting zOff or a sampled row
    pub fn get_offset(&self, row: u64) -> AlignResult<u64> {
        // Check for zOff (the $ character position)
        if row == self.z_off as u64 {
            return Ok(0);
        }

        // Check if row is at a sampled position (multiple of 2^offRate)
        // off_mask has all bits set except the lower offRate bits
        // (row & off_mask) == row means lower offRate bits are all 0
        if (row & self.params.off_mask) == row {
            let sample_idx = (row >> self.params.off_rate) as usize;
            let offset = self.offs(sample_idx);
            return Ok(offset as u64);
        }

        // Walk back using LF-mapping until hitting zOff or a sampled row
        let mut current_row = row;
        let mut jumps = 0u64;

        // No theoretical maximum - we just walk until we hit a known position
        loop {
            // Perform LF-mapping step
            let c = self.bwt_char(current_row)?;
            let fchr_c = self.fchr(c as u8);
            let rank_c = self.rank_internal(c, current_row)?;
            current_row = fchr_c + rank_c as u64;
            jumps += 1;

            // Check if we hit zOff
            if current_row == self.z_off as u64 {
                return Ok(jumps);
            }

            // Check if we hit a sampled row
            if (current_row & self.params.off_mask) == current_row {
                let sample_idx = (current_row >> self.params.off_rate) as usize;
                let offset = self.offs(sample_idx);
                return Ok(jumps + offset as u64);
            }

            // Safety limit to prevent infinite loops (should never hit this in practice)
            if jumps > self.params.len {
                return Err(AlignError::InvalidFormat(format!(
                    "getOffset failed to resolve row {} after {} steps (exceeded genome length)",
                    row, jumps
                )));
            }
        }
    }

    /// Get the character at a given BWT position
    fn bwt_char(&self, pos: u64) -> AlignResult<usize> {
        let data = self
            .ebwt_fw
            .as_ref()
            .ok_or_else(|| AlignError::InvalidFormat("Index not loaded".into()))?;

        let side_bwt_len = self.params.side_bwt_len as u64;
        let side_sz = self.params.side_sz as usize;

        let side_num = pos / side_bwt_len;
        let char_off = (pos % side_bwt_len) as usize;

        let side_start = self.ebwt_off + side_num as usize * side_sz;
        let byte_off = side_start + (char_off >> 2);
        let bit_off = (char_off & 3) * 2;

        if byte_off < data.len() {
            let byte = data[byte_off];
            Ok(((byte >> bit_off) & 0x03) as usize)
        } else {
            Err(AlignError::InvalidFormat(format!(
                "BWT position {} out of bounds",
                pos
            )))
        }
    }

    /// Convert a genome offset to chromosome index and position.
    ///
    /// This implements RustAlign's joinedToTextOff() algorithm using binary search
    /// on the rstarts array.
    ///
    /// Returns (chromosome_index, position_in_chromosome, chromosome_length)
    pub fn joined_to_text_off(
        &self,
        qlen: u64,
        off: u64,
        fw: bool,
    ) -> AlignResult<(u32, u64, u64)> {
        if self.rstarts.is_empty() {
            return Err(AlignError::InvalidFormat("rstarts array not loaded".into()));
        }

        let n_frag = self.n_frag as usize;

        // Binary search for the fragment containing this offset
        let mut top = 0usize;
        let mut bot = n_frag;

        while top < bot {
            let mid = top + ((bot - top) >> 1);

            // rstarts format: [frag_start, text_id, frag_offset] for each fragment
            let lower = self.rstarts[mid * 3] as u64;
            let upper = if mid + 1 >= n_frag {
                self.params.len
            } else {
                self.rstarts[(mid + 1) * 3] as u64
            };

            if lower <= off {
                if upper > off {
                    // Found the right fragment
                    let tidx = self.rstarts[mid * 3 + 1]; // chromosome index
                    let frag_off_base = self.rstarts[mid * 3 + 2]; // fragment offset within text

                    // Calculate position within the fragment
                    let frag_len = upper - lower;
                    let mut frag_off = off - lower;

                    // Handle reverse index orientation
                    if !fw {
                        if frag_off + qlen > frag_len {
                            // Would straddle fragment boundary
                            return Err(AlignError::InvalidFormat(
                                "Alignment straddles fragment boundary".into(),
                            ));
                        }
                        frag_off = frag_len - frag_off - 1;
                        frag_off = frag_off.saturating_sub(qlen - 1);
                    }

                    // Add fragment offset within the chromosome
                    let text_off = frag_off + frag_off_base as u64;

                    // Get chromosome length from plen
                    let tlen = if (tidx as usize) < self.plen.len() {
                        self.plen[tidx as usize] as u64
                    } else {
                        0
                    };

                    return Ok((tidx, text_off, tlen));
                } else {
                    top = mid + 1;
                }
            } else {
                bot = mid;
            }
        }

        Err(AlignError::InvalidFormat(format!(
            "Failed to find fragment for offset {}",
            off
        )))
    }

    /// Resolve a BWT row to chromosome name and position.
    ///
    /// This combines getOffset() and joinedToTextOff() to convert a BWT row
    /// directly to a chromosome name and position.
    ///
    /// Returns (chromosome_name, position, chromosome_length)
    pub fn resolve_position(
        &self,
        row: u64,
        qlen: u64,
        fw: bool,
    ) -> AlignResult<(String, u64, u64)> {
        // Step 1: Convert BWT row to genome offset
        let genome_off = self.get_offset(row)?;

        // Step 2: Convert genome offset to chromosome + position
        let (chr_idx, pos, chr_len) = self.joined_to_text_off(qlen, genome_off, fw)?;

        // Step 3: Get chromosome name
        let chr_name = if (chr_idx as usize) < self.refnames.len() {
            self.refnames[chr_idx as usize].clone()
        } else {
            format!("chr{}", chr_idx)
        };

        // SAM format uses 1-based positions
        Ok((chr_name, pos + 1, chr_len))
    }
}

impl FmIndex for Ebwt {
    /// Exact match search using backward search on FM-index
    ///
    /// This implements the classic FM-index backward search algorithm:
    /// 1. Initialize range using ftab for the last ftabChars nucleotides
    /// 2. Extend left through remaining characters using LF-mapping
    /// 3. Return all positions in the final range
    fn exact_match(&self, pattern: &[Nuc]) -> AlignResult<Vec<usize>> {
        if pattern.is_empty() {
            return Ok(Vec::new());
        }

        // Filter out invalid nucleotides (N's)
        let valid_pattern: Vec<Nuc> = pattern.iter().filter(|&&n| n.is_valid()).copied().collect();
        if valid_pattern.is_empty() {
            return Ok(Vec::new());
        }

        let pat_len = valid_pattern.len();
        let ftab_len = self.params.ftab_chars as usize;

        // Step 1: Initialize range using ftab
        // Use the last min(ftab_len, pat_len) characters
        let init_len = ftab_len.min(pat_len);
        let ftab_start = pat_len - init_len;

        // Compute ftab index from the last init_len characters
        let mut ftab_idx: usize = 0;
        for i in 0..init_len {
            let c = valid_pattern[ftab_start + i] as u8;
            ftab_idx = (ftab_idx << 2) | (c as usize);
        }

        // Get initial range from ftab
        let top = self.ftab(ftab_idx);
        let bot = self.ftab(ftab_idx + 1);

        if bot <= top {
            // No matches
            return Ok(Vec::new());
        }

        // Step 2: Extend left through remaining characters using LF-mapping
        let mut top = top;
        let mut bot = bot;

        // Process characters from ftab_start-1 down to 0
        for i in (0..ftab_start).rev() {
            let c = valid_pattern[i] as u8;

            // LF-mapping: top' = fchr[c] + rank(c, top)
            //                       bot' = fchr[c] + rank(c, bot)
            let fchr_c = self.fchr(c);
            let rank_top = self.rank_internal(c as usize, top)?;
            let rank_bot = self.rank_internal(c as usize, bot)?;

            top = fchr_c + rank_top as u64;
            bot = fchr_c + rank_bot as u64;

            if bot <= top {
                // No matches
                return Ok(Vec::new());
            }
        }

        // Step 3: Convert range to positions
        // For now, return all positions in the range
        // In a full implementation, we'd resolve these via the suffix array
        Ok((top..bot).map(|i| i as usize).collect())
    }

    fn len(&self) -> usize {
        self.params.len as usize
    }

    fn is_empty(&self) -> bool {
        self.params.len == 0
    }

    fn bwt(&self, pos: usize) -> AlignResult<NucPair> {
        if pos >= self.params.bwt_len as usize {
            return Err(AlignError::InvalidArgument(
                "BWT position out of range".into(),
            ));
        }

        // In the full implementation, this would:
        // 1. Calculate SideLocus from position
        // 2. Extract the nucleotide pair from mmapped data
        Ok(NucPair::default())
    }

    fn rank(&self, c: Nuc, _pos: usize) -> AlignResult<usize> {
        if !c.is_valid() {
            return Err(AlignError::InvalidNucleotide(c as u8));
        }

        // In the full implementation, this would:
        // 1. Use the side-based rank structure
        // 2. Count occurrences of c up to position
        Ok(0)
    }

    fn select(&self, c: Nuc, _k: usize) -> AlignResult<usize> {
        if !c.is_valid() {
            return Err(AlignError::InvalidNucleotide(c as u8));
        }

        // In the full implementation, this would:
        // 1. Use the side-based select structure
        // 2. Find position of k-th occurrence of c
        Ok(0)
    }
}

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

    #[test]
    fn test_ebwt_default() {
        let params = EbwtParams::default();
        assert_eq!(params.len, 0);
    }

    #[test]
    fn test_ebwt_load_nonexistent() {
        let result = Ebwt::load("/nonexistent/path/lambda_virus");
        assert!(result.is_err());
    }

    // More tests would require actual index files or mock data
}