mini_myers 0.1.2

use crate::constant::*;
use crate::tqueries::*;
use wide;
use wide::{i32x8, CmpEq, CmpGt};

/// Match information including position and cost for a query.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct MatchInfo {
    /// Index of the query (0-based)
    pub query_idx: usize,
    /// Minimum edit distance found (0 to k), including overhang penalties when used
    pub cost: f32,
    /// Position in target where best match ends (0-based)
    pub pos: i32,
}

/// Holds the reusable state vectors for the Myers search.
pub struct MyersSearchState {
    pv: Vec<i32x8>,
    mv: Vec<i32x8>,
    score: Vec<i32x8>,
    best_adjusted: Vec<i32x8>,
    /// The output vector.
    pub result: Vec<f32>,
}

impl MyersSearchState {
    /// Creates a new, empty state holder.
    pub fn new() -> Self {
        Self {
            pv: Vec::new(),
            mv: Vec::new(),
            score: Vec::new(),
            best_adjusted: Vec::new(),
            result: Vec::new(),
        }
    }

    /// Resizes and re-initializes state for `mini_search`.
    pub fn reset_for_mini_search(
        &mut self,
        nq: usize, // number of queries
        m: usize,  // query length
        initial_adjusted: i32,
    ) {
        let vectors_in_block = nq.div_ceil(SIMD_LANES);

        let all_ones = i32x8::splat(!0);
        let zero_v = i32x8::splat(0);
        let m_v = i32x8::splat(m as i32);
        let initial_adj_v = i32x8::splat(initial_adjusted);
        self.pv.resize(vectors_in_block, all_ones);
        self.pv.fill(all_ones);
        self.mv.resize(vectors_in_block, zero_v);
        self.mv.fill(zero_v);
        self.score.resize(vectors_in_block, m_v);
        self.score.fill(m_v);
        self.best_adjusted.resize(vectors_in_block, initial_adj_v);
        self.best_adjusted.fill(initial_adj_v);
        self.result.resize(nq, -1.0f32);
        self.result.fill(-1.0f32);
    }

    /// Resizes and re-initializes state for `mini_search_with_positions`.
    pub fn reset_for_position_search(
        &mut self,
        nq: usize, // number of queries
        m: usize,  // query length
    ) {
        let vectors_in_block = nq.div_ceil(SIMD_LANES);
        let all_ones = i32x8::splat(!0);
        let zero_v = i32x8::splat(0);
        let m_v = i32x8::splat(m as i32);
        self.pv.resize(vectors_in_block, all_ones);
        self.pv.fill(all_ones);
        self.mv.resize(vectors_in_block, zero_v);
        self.mv.fill(zero_v);
        self.score.resize(vectors_in_block, m_v);
        self.score.fill(m_v);
    }
}

impl Default for MyersSearchState {
    fn default() -> Self {
        Self::new()
    }
}

/// Searches for all queries in the target sequence using the Myers algorithm.
/// # Arguments
///
/// * `state` - A `MyersSearchState` to hold and reuse allocations.
/// * `transposed` - A `TQueries` structure containing the preprocessed queries
/// * `target` - The target DNA sequence to search in (as a byte slice)
/// * `k` - Maximum edit distance threshold.
/// * `alpha` - Optional overhang penalty.
///
/// # Returns
///
/// A slice `&[f32]` of the results, valid until `state` is next modified.
/// - The minimum edit distance found (0 to k)
/// - `-1.0` if no match within the threshold k was found
#[inline(always)]
pub fn mini_search<'a>(
    state: &'a mut MyersSearchState,
    transposed: &TQueries,
    target: &[u8],
    k: u8,
    alpha: Option<f32>,
) -> &'a [f32] {
    let nq = transposed.n_queries;
    if nq == 0 {
        return &[];
    }
    let m = transposed.query_length;
    assert!(m > 0 && m <= 32, "query length must be 1..=32");

    let alpha = alpha.unwrap_or(1.0).clamp(0.0, 1.0);
    const SCALE_SHIFT: u32 = 8;
    const SCALE: i32 = 1 << SCALE_SHIFT; // 256
    let alpha_scaled = ((alpha * (SCALE as f32)).round() as i32).clamp(0, SCALE);
    let initial_adjusted = (m as i32) * alpha_scaled;
    state.reset_for_mini_search(nq, m, initial_adjusted);
    search_simd_core(
        state,
        transposed,
        target,
        k,
        alpha_scaled,
        SCALE,
        SCALE_SHIFT,
    );
    &state.result[..nq]
}

/// Searches for all queries in the target sequence, returning detailed match information.
///
///
/// # Arguments
///
/// * `state` - A `MyersSearchState` to reuse allocations.
/// * `transposed` - A `TQueries` structure containing the preprocessed queries
/// * `target` - The target DNA sequence to search in (as a byte slice)
/// * `k` - Maximum edit distance threshold.
/// * `alpha` - Optional overhang penalty (see [`mini_search`]).
/// * `results` - Vector that will be populated with all matches found.
#[inline(always)]
pub fn mini_search_with_positions(
    state: &mut MyersSearchState,
    transposed: &TQueries,
    target: &[u8],
    k: u8,
    alpha: Option<f32>,
    results: &mut Vec<MatchInfo>,
) {
    let nq = transposed.n_queries;
    if nq == 0 {
        return;
    }
    let m = transposed.query_length;
    assert!(m > 0 && m <= 32, "query length must be 1..=32");

    let alpha = alpha.unwrap_or(1.0).clamp(0.0, 1.0);
    const SCALE_SHIFT: u32 = 8;
    const SCALE: i32 = 1 << SCALE_SHIFT;
    let alpha_scaled = ((alpha * (SCALE as f32)).round() as i32).clamp(0, SCALE);
    state.reset_for_position_search(nq, m);
    search_simd_with_positions_core(
        state,
        transposed,
        target,
        k,
        alpha_scaled,
        SCALE,
        SCALE_SHIFT,
        results,
    )
}

#[derive(Copy, Clone)]
struct OverhangColumnCtx {
    all_ones: i32x8,
    zero_v: i32x8,
    one_v: i32x8,
    mask_vec: i32x8,
    scale_shift: u32,
}

#[inline(always)]
fn advance_column(
    pv: &mut i32x8,
    mv: &mut i32x8,
    score: &mut i32x8,
    eq: i32x8,
    adjust_vec: i32x8,
    ctx: OverhangColumnCtx,
) -> i32x8 {
    let all_ones = ctx.all_ones;
    let zero_v = ctx.zero_v;
    let one_v = ctx.one_v;
    let mask_vec = ctx.mask_vec;
    let pv_val = *pv;
    let mv_val = *mv;
    let xv = eq | mv_val;
    let xh = (((eq & pv_val) + pv_val) ^ pv_val) | eq;
    let ph = mv_val | (all_ones ^ (xh | pv_val));
    let mh = pv_val & xh;
    let ph_bit_mask = (ph & mask_vec).simd_eq(zero_v);
    let ph_bit = (all_ones ^ ph_bit_mask) & one_v;
    let mh_bit_mask = (mh & mask_vec).simd_eq(zero_v);
    let mh_bit = (all_ones ^ mh_bit_mask) & one_v;
    let ph_shift = ph << 1;
    let new_pv = (mh << 1) | (all_ones ^ (xv | ph_shift));
    let new_mv = ph_shift & xv;
    *pv = new_pv;
    *mv = new_mv;
    let new_score = *score + (ph_bit - mh_bit);
    *score = new_score;
    let new_score_scaled = new_score << ctx.scale_shift;
    new_score_scaled - adjust_vec
}

/// Core SIMD Myers implementation re-using allocs
#[inline(always)]
pub(crate) fn search_simd_core(
    state: &mut MyersSearchState,
    transposed: &TQueries,
    target: &[u8],
    k: u8,
    alpha_scaled: i32,
    scale: i32,
    scale_shift: u32,
) {
    let nq = transposed.n_queries;
    let m = transposed.query_length;
    let vectors_in_block = nq.div_ceil(SIMD_LANES);

    let extra_penalty_scaled = scale - alpha_scaled;
    let k_scaled = (k as i32) << scale_shift;
    let inv_scale = 1.0f32 / (scale as f32);

    let all_ones = i32x8::splat(!0);
    let zero_v = i32x8::splat(0);
    let one_v = i32x8::splat(1);
    let zero_eq_vec = i32x8::splat(0);

    let high_bit: u32 = 1u32 << (m - 1);
    let mask_vec = i32x8::splat(high_bit as i32);
    let k_scaled_vec = i32x8::splat(k_scaled);
    let ctx = OverhangColumnCtx {
        all_ones,
        zero_v,
        one_v,
        mask_vec,
        scale_shift,
    };

    for (idx, &tb) in target.iter().enumerate() {
        let encoded = crate::iupac::get_encoded(tb);

        let overhang = m.saturating_sub(idx + 1) as i32;
        let adjust_vec = i32x8::splat(overhang * extra_penalty_scaled);
        let peq_slice = unsafe { transposed.peqs.get_unchecked(encoded as usize) };

        unsafe {
            for block_idx in 0..vectors_in_block {
                let eq = *peq_slice.get_unchecked(block_idx);
                let adjusted = advance_column(
                    state.pv.get_unchecked_mut(block_idx),
                    state.mv.get_unchecked_mut(block_idx),
                    state.score.get_unchecked_mut(block_idx),
                    eq,
                    adjust_vec,
                    ctx,
                );
                let best_slot = state.best_adjusted.get_unchecked_mut(block_idx);
                *best_slot = (*best_slot).min(adjusted);
            }
        }
    }

    if alpha_scaled < scale {
        for trailing in 1..=m {
            let adjust_vec = i32x8::splat((trailing as i32) * extra_penalty_scaled);
            unsafe {
                for block_idx in 0..vectors_in_block {
                    let adjusted = advance_column(
                        state.pv.get_unchecked_mut(block_idx),
                        state.mv.get_unchecked_mut(block_idx),
                        state.score.get_unchecked_mut(block_idx),
                        zero_eq_vec,
                        adjust_vec,
                        ctx,
                    );
                    let best_slot = state.best_adjusted.get_unchecked_mut(block_idx);
                    *best_slot = (*best_slot).min(adjusted);
                }
            }
        }
    }

    let neg_mask = i32x8::splat(i32::MAX);

    for (block_idx, &min_adj) in state.best_adjusted.iter().enumerate() {
        let mask = all_ones ^ min_adj.simd_gt(k_scaled_vec);
        let selected = mask.blend(min_adj, neg_mask);
        let base = block_idx * SIMD_LANES;
        let end = (base + SIMD_LANES).min(nq);
        let selected_arr = selected.to_array();

        for (lane_idx, &val) in selected_arr[..(end - base)].iter().enumerate() {
            if val != i32::MAX {
                unsafe {
                    *state.result.get_unchecked_mut(base + lane_idx) = (val as f32) * inv_scale;
                }
            }
        }
    }
}

#[inline(always)]
fn search_simd_with_positions_core(
    state: &mut MyersSearchState,
    transposed: &TQueries,
    target: &[u8],
    k: u8,
    alpha_scaled: i32,
    scale: i32,
    scale_shift: u32,
    results: &mut Vec<MatchInfo>,
) {
    let all_ones = i32x8::splat(!0);
    let zero_v = i32x8::splat(0);
    let one_v = i32x8::splat(1);

    let nq = transposed.n_queries;
    let m = transposed.query_length;
    let vectors_in_block = nq.div_ceil(SIMD_LANES);

    let extra_penalty_scaled = scale - alpha_scaled;
    let k_scaled = (k as i32) << scale_shift;
    let k_scaled_vec = i32x8::splat(k_scaled);
    let inv_scale = 1.0f32 / (scale as f32);

    let high_bit: u32 = 1u32 << (m - 1);
    let mask_vec = i32x8::splat(high_bit as i32);
    let zero_eq_vec = i32x8::splat(0);
    let ctx = OverhangColumnCtx {
        all_ones,
        zero_v,
        one_v,
        mask_vec,
        scale_shift,
    };

    for (idx, &tb) in target.iter().enumerate() {
        let encoded = crate::iupac::get_encoded(tb);
        let overhang = m.saturating_sub(idx + 1) as i32;
        let adjust_vec = i32x8::splat(overhang * extra_penalty_scaled);
        let peq_slice = unsafe { transposed.peqs.get_unchecked(encoded as usize) };
        let pos = idx as i32;

        unsafe {
            for block_idx in 0..vectors_in_block {
                let eq = *peq_slice.get_unchecked(block_idx);
                let adjusted = advance_column(
                    state.pv.get_unchecked_mut(block_idx),
                    state.mv.get_unchecked_mut(block_idx),
                    state.score.get_unchecked_mut(block_idx),
                    eq,
                    adjust_vec,
                    ctx,
                );
                let match_mask = all_ones ^ adjusted.simd_gt(k_scaled_vec);
                let match_bits = match_mask.to_array();
                if match_bits.iter().any(|&b| b != 0) {
                    let adjusted_arr = adjusted.to_array();
                    let base = block_idx * SIMD_LANES;
                    let end = (base + SIMD_LANES).min(nq);

                    for (i, &score_scaled) in adjusted_arr[..(end - base)].iter().enumerate() {
                        if score_scaled <= k_scaled {
                            results.push(MatchInfo {
                                query_idx: base + i,
                                cost: score_scaled as f32 * inv_scale,
                                pos,
                            });
                        }
                    }
                }
            }
        }
    }

    if alpha_scaled < scale {
        for trailing in 1..=m {
            let adjust_vec = i32x8::splat((trailing as i32) * extra_penalty_scaled);
            let pos = (target.len() + trailing - 1) as i32;

            unsafe {
                for block_idx in 0..vectors_in_block {
                    let adjusted = advance_column(
                        state.pv.get_unchecked_mut(block_idx),
                        state.mv.get_unchecked_mut(block_idx),
                        state.score.get_unchecked_mut(block_idx),
                        zero_eq_vec,
                        adjust_vec,
                        ctx,
                    );
                    let match_masks = all_ones ^ adjusted.simd_gt(k_scaled_vec);
                    let match_bits = match_masks.to_array();
                    if match_bits.iter().any(|&b| b != 0) {
                        let adjusted_arr = adjusted.to_array();
                        let base = block_idx * SIMD_LANES;
                        let end = (base + SIMD_LANES).min(nq);

                        for (i, &score_scaled) in adjusted_arr[..(end - base)].iter().enumerate() {
                            if score_scaled <= k_scaled {
                                results.push(MatchInfo {
                                    query_idx: base + i,
                                    cost: score_scaled as f32 * inv_scale,
                                    pos,
                                });
                            }
                        }
                    }
                }
            }
        }
    }
}
#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_iupac_query_matches_standard_target() {
        let queries = vec![b"n".to_vec()];
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        assert_eq!(
            mini_search(&mut state, &transposed, b"A", 0, None),
            vec![0.0]
        );
        assert_eq!(
            mini_search(&mut state, &transposed, b"a", 0, None),
            vec![0.0]
        );
    }

    #[test]
    fn test_iupac_target_matches_standard_query() {
        let queries = vec![b"A".to_vec()];
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        assert_eq!(
            mini_search(&mut state, &transposed, b"N", 0, None),
            vec![0.0]
        );
        assert_eq!(
            mini_search(&mut state, &transposed, b"n", 0, None),
            vec![0.0]
        );
    }

    #[test]
    fn test_iupac_mismatch_requires_edit() {
        let queries = vec![b"R".to_vec()];
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        assert_eq!(
            mini_search(&mut state, &transposed, b"C", 0, None),
            vec![-1.0]
        );
        assert_eq!(
            mini_search(&mut state, &transposed, b"C", 1, None),
            vec![1.0]
        );
        assert_eq!(
            mini_search(&mut state, &transposed, b"c", 1, None),
            vec![1.0]
        );
    }

    #[test]
    #[should_panic(expected = "invalid IUPAC character")]
    fn test_invalid_query_panics() {
        let queries = vec![b"AZ".to_vec()];
        let _ = TQueries::new(&queries);
    }

    // #[test]
    // #[should_panic(expected = "invalid IUPAC character")]
    // fn test_invalid_target_panics() {
    //     let queries = vec![b"A".to_vec()];
    //     let transposed = TQueries::new(&queries);
    //     let mut state = MyersSearchState::new();
    //     // This should panic
    //     let _ = mini_search(&mut state, &transposed, b"Z", 1, None);
    // }

    #[test]
    fn test_mini_search() {
        let q = b"ATG".to_vec();
        let queries = vec![q];
        let t = b"CCCCCCCCCATGCCCCC";
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![0.0]); // Lowest edits, 0 for q at idx 0
    }

    #[test]
    fn test_double_search() {
        let q1 = b"ATG".to_vec();
        let q2 = b"TTG".to_vec();
        let queries = vec![q1, q2];
        let t = b"CCCTTGCCCCCCATGCCCCC";
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![0.0, 0.0]);
    }

    #[test]
    fn test_edit() {
        let q1 = b"ATCAGA";
        let queries = vec![q1.to_vec()];
        let t = b"ATCTGA"; // 1 edit
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![1.0]);
        let t = b"GTCTGA"; // 2 edits
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![2.0]);
        let t = b"GTTGA"; // 3 edits (1 del)
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![3.0]);
        // Match should not be recovered when k == 1
        let result = mini_search(&mut state, &transposed, t, 1, None);
        assert_eq!(result, vec![-1.0]);
    }

    #[test]
    fn test_lowest_edits_returned() {
        let q1 = b"GGGCATCGATGAC";
        let queries = vec![q1.to_vec()];
        let t = b"CCCCCCCGGGCATCGATGACCCCCCCCCCCCCCCGGGCTTCGATGAC";
        //                                                     ^^^^^^^^^^^^^ has one mutation
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![0.0]);
    }

    #[test]
    #[should_panic(expected = "All queries must have the same length")]
    fn test_error_unequal_lengths() {
        let q1 = b"GGGCATCGATGAC";
        let q2 = b"AAA";
        let t = b"CCC";
        let queries = vec![q1.to_vec(), q2.to_vec()];
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 4, None);
        assert_eq!(result, vec![0.0]);
    }

    #[test]
    fn read_example() {
        let queries = vec![b"ATG".to_vec(), b"TTG".to_vec()];
        let transposed = TQueries::new(&queries);
        let target = b"CCCTCGCCCCCCATGCCCCC";
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, target, 4, None);
        println!("Result: {:?}", result);
        assert_eq!(result, vec![0.0, 1.0]);
    }

    #[test]
    fn test_mini_search_with_positions() {
        let q = b"ATG".to_vec();
        let queries = vec![q];
        let t = b"CCCCCCCCCATGCCCCC";
        //             012345678901234567
        //                      ATG at pos 9-11 (ends at 11)
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let mut result = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 4, None, &mut result);

        // Should find at least one match
        assert!(!result.is_empty());
        // Find the exact match (cost 0)
        let exact_matches: Vec<_> = result.iter().filter(|m| m.cost == 0.0).collect();
        assert!(!exact_matches.is_empty());
        assert_eq!(exact_matches[0].query_idx, 0);
        assert_eq!(exact_matches[0].pos, 11); // Position where match ends
    }

    #[test]
    fn test_positions_multiple_occurrences() {
        let queries = vec![b"ATG".to_vec()];
        let t = b"ATGCCCATGCCCATG";
        //        012345678901234
        // ATG at positions 0-2 (ends at 2), 6-8 (ends at 8), 12-14 (ends at 14)
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let mut result = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 0, None, &mut result);

        // Should find all 3 exact matches
        let exact_matches: Vec<_> = result
            .iter()
            .filter(|m| (m.cost - 0.0).abs() < f32::EPSILON && m.query_idx == 0)
            .collect();
        assert_eq!(exact_matches.len(), 3);
        assert_eq!(exact_matches[0].pos, 2);
        assert_eq!(exact_matches[1].pos, 8);
        assert_eq!(exact_matches[2].pos, 14);
    }

    #[test]
    fn test_positions_match_costs() {
        let q1 = b"ATG".to_vec();
        let q2 = b"TTG".to_vec();
        let queries = vec![q1, q2];
        let t = b"CCCTTGCCCCCCATGCCCCC";
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();

        let result_basic = mini_search(&mut state, &transposed, t, 4, None).to_vec();
        let mut result_pos = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 4, None, &mut result_pos);

        // mini_search returns minimum cost per query
        // mini_search_with_positions returns all matches
        // Check that the minimum cost in positions matches mini_search result
        for (query_idx, _) in queries.iter().enumerate() {
            let min_cost_in_positions = result_pos
                .iter()
                .filter(|m| m.query_idx == query_idx)
                .map(|m| m.cost)
                .fold(f32::INFINITY, f32::min);
            let min_cost_in_positions = if min_cost_in_positions.is_finite() {
                min_cost_in_positions
            } else {
                -1.0
            };
            assert!((result_basic[query_idx] - min_cost_in_positions).abs() < 1e-6);
        }
    }

    #[test]
    fn test_positions_no_match() {
        let q = b"AAAAA".to_vec();
        let queries = vec![q];
        let t = b"CCCCCCCCCCCCC";
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let mut result = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 1, None, &mut result);

        // Should have no matches for query 0 with cost <= 1
        let matches_q0: Vec<_> = result.iter().filter(|m| m.query_idx == 0).collect();
        assert_eq!(matches_q0.len(), 0);
    }

    #[test]
    fn test_positions_multiple_queries() {
        let queries = vec![
            b"AAAA".to_vec(),
            b"TTTT".to_vec(),
            b"GGGG".to_vec(),
            b"CCCC".to_vec(),
        ];
        let t = b"AAGGGGTTTTCCCC";
        let transposed = TQueries::new(&queries);
        let mut state = MyersSearchState::new();
        let mut result = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 2, None, &mut result);

        // All queries should have at least one match with k=2
        for query_idx in 0..4 {
            let query_matches: Vec<_> =
                result.iter().filter(|m| m.query_idx == query_idx).collect();
            assert!(!query_matches.is_empty());
            // All matches should be within threshold
            for m in &query_matches {
                assert!(m.cost >= 0.0 && m.cost <= 2.0);
                assert!(m.pos >= 0);
            }
        }
    }
    #[test]
    fn test_overhang_half_penalty() {
        let queries = vec![b"AC".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"A";
        let mut state = MyersSearchState::new();
        let result = mini_search(&mut state, &transposed, t, 1, Some(0.5));
        assert_eq!(result.len(), 1);
        assert!(
            (result[0] - 0.5).abs() < 1e-6,
            "expected 0.5, got {}",
            result[0]
        );
    }

    #[test]
    fn test_overhang_matches_standard_when_alpha_one() {
        let queries = vec![b"ATG".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"CCCCCCCCCATGCCCCC";
        let mut state = MyersSearchState::new();
        let standard = mini_search(&mut state, &transposed, t, 4, None).to_vec(); // clone
        let overhang = mini_search(&mut state, &transposed, t, 4, Some(1.0));
        assert_eq!(standard.len(), overhang.len());
        for (i, &val) in standard.iter().enumerate() {
            assert!((overhang[i] - val).abs() < 1e-6);
        }
    }

    #[test]
    fn test_positions_overhang() {
        let queries = vec![b"AC".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"A";
        let mut state = MyersSearchState::new();
        let mut results = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 1, Some(0.5), &mut results);
        assert!(!results.is_empty());
        let min_cost_entry = results
            .iter()
            .min_by(|a, b| a.cost.partial_cmp(&b.cost).unwrap())
            .unwrap();
        assert_eq!(min_cost_entry.query_idx, 0);
        assert!((min_cost_entry.cost - 0.5).abs() < 1e-6);
    }

    #[test]
    fn test_longer_overhang_left() {
        let queries = vec![b"AAACCC".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"CCCGGGGGGGG";
        let mut state = MyersSearchState::new();
        let mut results = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 4, Some(0.5), &mut results);
        // Get minimum cost match
        let min_cost = results
            .iter()
            .map(|m| m.cost)
            .fold(f32::INFINITY, |a, b| a.min(b));
        assert!(
            (min_cost - 1.5).abs() < 1e-6,
            "expected 1.5, got {}",
            min_cost
        );
    }

    #[test]
    fn test_longer_overhang_right() {
        let queries = vec![b"AAACCC".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"GGGGGAAA";
        let mut state = MyersSearchState::new();
        let mut results = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 4, Some(0.5), &mut results);
        // Get minimum cost match
        let min_cost = results
            .iter()
            .map(|m| m.cost)
            .fold(f32::INFINITY, |a, b| a.min(b));
        assert!(
            (min_cost - 1.5).abs() < 1e-6,
            "expected 1.5, got {}",
            min_cost
        );
    }

    #[test]
    fn test_overhang_with_positions() {
        let queries = vec![b"AAACCC".to_vec()];
        let transposed = TQueries::new(&queries);
        let t = b"GGGGGAAA";
        let mut state = MyersSearchState::new();
        let mut results = Vec::new();
        mini_search_with_positions(&mut state, &transposed, t, 4, Some(0.5), &mut results);
        // Get minimum cost match
        let _min_cost = results
            .iter()
            .map(|m| m.cost)
            .fold(f32::INFINITY, |a, b| a.min(b));
        println!("results: {:?}", results);
        assert!(!results.is_empty());
    }
}