poa-consensus 0.1.1

use crate::graph::AlignOp;
use crate::{self as poa_consensus, AlignmentMode, ConsensusMode, PoaConfig, PoaError, PoaGraph};

fn b(s: &str) -> Vec<u8> {
    s.as_bytes().to_vec()
}
fn s(v: &[u8]) -> String {
    String::from_utf8_lossy(v).into_owned()
}

/// Build consensus from reads; seed_idx selects the first read added to the graph.
fn consensus(reads: &[Vec<u8>], seed_idx: usize) -> Vec<u8> {
    let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
    for (i, read) in reads.iter().enumerate() {
        if i == seed_idx {
            continue;
        }
        graph.add_read(read).unwrap();
    }
    graph.consensus().unwrap().sequence
}

fn consensus_cfg(reads: &[Vec<u8>], seed_idx: usize, cfg: PoaConfig) -> Vec<u8> {
    let mut graph = PoaGraph::new(&reads[seed_idx], cfg).unwrap();
    for (i, read) in reads.iter().enumerate() {
        if i == seed_idx {
            continue;
        }
        graph.add_read(read).unwrap();
    }
    graph.consensus().unwrap().sequence
}

// ── Error cases ───────────────────────────────────────────────────────────────

#[test]
fn empty_reads() {
    let result = PoaGraph::new(&[], PoaConfig::default());
    assert!(
        matches!(result, Err(PoaError::EmptyInput)),
        "expected EmptyInput"
    );
}

#[test]
fn below_min_reads() {
    let cfg = PoaConfig {
        min_reads: 3,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&b("ACGT"), cfg).unwrap();
    graph.add_read(&b("ACGT")).unwrap();
    let result = graph.consensus();
    assert!(
        matches!(result, Err(PoaError::InsufficientDepth { got: 2, min: 3 })),
        "expected InsufficientDepth, got {:?}",
        result
    );
}

#[test]
fn seed_out_of_bounds() {
    // Callers must check seed_idx < reads.len() before calling PoaGraph::new.
    // Document the expected guard pattern.
    let reads = vec![b("ACGT"), b("ACGT")];
    let idx = 5usize;
    assert!(idx >= reads.len(), "caller must guard seed_idx before use");
}

// ── Basic correctness ─────────────────────────────────────────────────────────

#[test]
fn single_read_passthrough() {
    let reads = vec![b("CATCATCAT")];
    assert_eq!(consensus(&reads, 0), b("CATCATCAT"));
}

#[test]
fn two_identical_reads() {
    let reads = vec![b("CATCATCAT"), b("CATCATCAT")];
    assert_eq!(consensus(&reads, 0), b("CATCATCAT"));
}

#[test]
fn majority_base_wins() {
    let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CGTCATCAT")];
    assert_eq!(s(&consensus(&reads, 0)), "CATCATCAT");
}

#[test]
fn single_outlier_not_inflated() {
    let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCATCAT")];
    assert_eq!(consensus(&reads, 0).len(), 9);
}

#[test]
fn length_variation_longer_wins() {
    let reads = vec![b("CATCATCATCAT"), b("CATCATCATCAT"), b("CATCATCAT")];
    assert_eq!(consensus(&reads, 0).len(), 12);
}

#[test]
fn no_inflation_with_length_noise() {
    let reads = vec![
        b("CAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAG"),
    ];
    assert_eq!(consensus(&reads, 0).len(), 15);
}

#[test]
fn no_inflation_phox2b_like() {
    let reads = vec![
        b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
        b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
        b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
    ];
    assert_eq!(consensus(&reads, 0).len(), 60);
}

#[test]
fn single_base_reads() {
    let reads = vec![b("A"), b("A"), b("A")];
    assert_eq!(consensus(&reads, 0), b("A"));
}

// ── Boundary trim ─────────────────────────────────────────────────────────────

#[test]
fn boundary_trim_leading_seed_artifact() {
    let reads = vec![
        b("XXXCATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "CATCATCAT", "got: {}", result);
}

#[test]
fn boundary_trim_trailing_seed_artifact() {
    let reads = vec![
        b("CATCATCATXXX"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "CATCATCAT", "got: {}", result);
}

// ── Diagnostic tests from ref/poa.rs ─────────────────────────────────────────

#[test]
fn diag_sca3_t3_tail_seed_t1() {
    let reads = vec![
        b("CAGCAGCAGT"),
        b("CAGCAGCAGTTT"),
        b("CAGCAGCAGTTT"),
        b("CAGCAGCAGTTT"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "CAGCAGCAGTTT", "got: {}", result);
}

#[test]
fn diag_sca3_t3_tail_seed_t3() {
    let reads = vec![
        b("CAGCAGCAGTTT"),
        b("CAGCAGCAGTTT"),
        b("CAGCAGCAGTTT"),
        b("CAGCAGCAGT"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "CAGCAGCAGTTT", "got: {}", result);
}

#[test]
fn diag_sca31_trailing_interrupt_seed_missing() {
    let reads = vec![
        b("ATTATTATTATT"),
        b("ATTATTATTATTATA"),
        b("ATTATTATTATTATA"),
        b("ATTATTATTATTATA"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "ATTATTATTATTATA", "got: {}", result);
}

#[test]
fn diag_sca3_interrupt_position_single_outlier() {
    let maj = b("CAGCAGCAGCAGCAGGTTCAGCAG");
    let out = b("CAGCAGCAGCAGCAGCAGGTTCAGCAG");
    let reads = vec![maj.clone(), maj.clone(), maj.clone(), out];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "CAGCAGCAGCAGCAGGTTCAGCAG", "got: {}", result);
}

#[test]
fn diag_sca8_minority_trailing_extension_trimmed() {
    let base = b("CAGCAGCAGCAGCAG");
    let extend = b("CAGCAGCAGCAGCAGGCT");
    let reads = vec![
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        extend.clone(),
        extend.clone(),
        extend.clone(),
    ];
    assert_eq!(consensus(&reads, 0).len(), 15);
}

#[test]
fn diag_sca8_minority_trailing_extension_seed_extends() {
    let base = b("CAGCAGCAGCAGCAG");
    let extend = b("CAGCAGCAGCAGCAGGCT");
    let reads = vec![
        extend.clone(),
        extend.clone(),
        extend.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
        base.clone(),
    ];
    assert_eq!(consensus(&reads, 0).len(), 15);
}

#[test]
fn diag_sca31_trailing_interrupt_before_flank() {
    let flank = b("GCGCGCGC");
    let mut seed_read = b("ATTATTATTATT");
    seed_read.extend_from_slice(&flank);
    let mut maj_read = b("ATTATTATTATTATA");
    maj_read.extend_from_slice(&flank);
    let reads = vec![
        seed_read,
        maj_read.clone(),
        maj_read.clone(),
        maj_read.clone(),
    ];
    let result = s(&consensus(&reads, 0));
    let expected: String = "ATTATTATTATTATA"
        .chars()
        .chain("GCGCGCGC".chars())
        .collect();
    assert_eq!(result, expected, "got: {}", result);
}

#[test]
fn diag_sca3_interrupt_position_long_repeat_with_flank() {
    let flank = b("CTGCTGCTG");
    let make = |repeat_pre: &str, interrupt: &str, repeat_post: &str| -> Vec<u8> {
        let mut v = repeat_pre.as_bytes().to_vec();
        v.extend_from_slice(interrupt.as_bytes());
        v.extend_from_slice(repeat_post.as_bytes());
        v.extend_from_slice(&flank);
        v
    };
    let maj = make("CAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAGCAG");
    let out = make("CAGCAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAG");
    let reads = vec![maj.clone(), maj.clone(), maj.clone(), out];
    let result = s(&consensus(&reads, 0));
    let expected = s(&make("CAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAGCAG"));
    assert_eq!(result, expected, "got: {}", result);
}

#[test]
#[ignore]
fn diag_frda_gaa_rotation_phase() {
    let gaa_phase = b("GAAGAAGAAGAA");
    let aag_phase = b("AAGAAGAAGAAG");
    let aga_phase = b("AGAAGAAGAAGA");
    let reads = vec![gaa_phase.clone(), gaa_phase.clone(), aag_phase, aga_phase];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result.len(), 12, "got: '{}'", result);
}

#[test]
fn diag_frda_gaa_rotation_with_flanking() {
    let make = |repeat: &str| -> Vec<u8> {
        let mut v = b("TTTCCC");
        v.extend_from_slice(repeat.as_bytes());
        v.extend_from_slice(b("GGGAAA").as_slice());
        v
    };
    let reads = vec![
        make("GAAGAAGAAGAA"),
        make("GAAGAAGAAGAA"),
        make("AAGAAGAAGAAG"),
        make("AGAAGAAGAAGA"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result.len(), 24, "got: '{}'", result);
}

#[test]
fn diag_phase_shift_first_node_coverage() {
    let reads = vec![
        b("GAAGAA"),
        b("GAAGAA"),
        b("GAAGAA"),
        b("GAAGAA"),
        b("AAGAAG"),
    ];
    let result = s(&consensus(&reads, 0));
    assert_eq!(result, "GAAGAA", "got: '{}'", result);
}

#[test]
fn diag_phase_shift_majority_trims_first_base() {
    let reads = vec![
        b("GAAGAA"),
        b("GAAGAA"),
        b("AAGAAG"),
        b("AAGAAG"),
        b("AAGAAG"),
    ];
    let result = s(&consensus(&reads, 0));
    // The majority is the same sequence in a different phase; correct length is 6.
    assert_eq!(result.len(), 6, "got: '{}'", result);
}

#[test]
fn diag_sca3_t3_tail_with_flank() {
    let flank = b("CCTCCTCCT");
    let make = |tail: &str| -> Vec<u8> {
        let mut v = b("CAGCAGCAG");
        v.extend_from_slice(tail.as_bytes());
        v.extend_from_slice(&flank);
        v
    };
    let reads = vec![make("T"), make("TTT"), make("TTT"), make("TTT")];
    let result = s(&consensus(&reads, 0));
    let expected = s(&make("TTT"));
    assert_eq!(result, expected, "got: {}", result);
}

// ── Real-data reproduction ────────────────────────────────────────────────────

#[test]
fn diag_sca8_real_sequences_no_flank() {
    let maj57 = b("TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT");
    let min75 = b("TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT");
    let min81 =
        b("TACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT");
    let mut reads: Vec<Vec<u8>> = std::iter::repeat(maj57.clone()).take(32).collect();
    reads.push(b(
        "TACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTAC",
    ));
    reads.push(b(
        "TACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCT",
    ));
    reads.push(b("TACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT"));
    reads.push(min75.clone());
    reads.push(min81.clone());
    reads.push(min81.clone());
    let result = consensus(&reads, 0);
    assert_eq!(
        result.len(),
        maj57.len(),
        "SCA8 consensus must match majority length {}, got len {}: '{}'",
        maj57.len(),
        result.len(),
        s(&result)
    );
}

#[test]
fn diag_sca8_real_sequences_with_flank() {
    let flank_l = b("GCTTCGAAGTC");
    let flank_r = b("AAACGGTTCCA");
    let make = |repeat: &[u8]| -> Vec<u8> {
        let mut v = flank_l.clone();
        v.extend_from_slice(repeat);
        v.extend_from_slice(&flank_r);
        v
    };
    let maj = make(&b(
        "TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT",
    ));
    let min75 = make(&b(
        "TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT",
    ));
    let min81 = make(&b(
        "TACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT",
    ));
    let mut reads: Vec<Vec<u8>> = std::iter::repeat(maj.clone()).take(32).collect();
    reads.push(min75);
    reads.push(min81.clone());
    reads.push(min81);
    let result_len = consensus(&reads, 0).len();
    assert_eq!(
        result_len,
        maj.len(),
        "SCA8 flanked: got {}, expected {}",
        result_len,
        maj.len()
    );
}

// ── Banded DP ─────────────────────────────────────────────────────────────────

#[test]
fn banded_matches_unbanded_small() {
    // Banded and unbanded must produce identical results when the band is wide
    // enough to cover the optimal path.
    let reads = vec![
        b("CAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAGCAG"),
        b("CAGCAGCAGCAGCAG"),
    ];
    let unbanded = consensus(&reads, 0);
    let cfg_banded = PoaConfig {
        band_width: 50,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg_banded);
    assert_eq!(unbanded, banded, "banded vs unbanded mismatch");
}

#[test]
fn adaptive_band_matches_unbanded() {
    let reads = vec![
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCATCAT"),
        b("CATCATCAT"),
    ];
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        adaptive_band: true,
        adaptive_band_b: 5,
        adaptive_band_f: 0.1,
        ..Default::default()
    };
    let adaptive = consensus_cfg(&reads, 0, cfg);
    assert_eq!(unbanded, adaptive, "adaptive band vs unbanded mismatch");
}

#[test]
fn band_too_narrow_fallback_to_unbanded() {
    // seed = 1 A, read = 30 A's: the 2-pass banded retry exhausts all banded
    // options but the 3-pass unbanded fallback recovers.  BandTooNarrow is now
    // an internal signal, not a user-visible error.
    let seed = b("A");
    let read = b("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"); // 30 A's
    let cfg = PoaConfig {
        band_width: 2,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&seed, cfg).unwrap();
    let result = graph.add_read(&read);
    assert!(
        result.is_ok(),
        "3-pass retry must recover via unbanded fallback, got {:?}",
        result.map(|_| ())
    );
}

#[test]
fn large_length_variance_banded() {
    // 3 reads of 15 bp majority + 1 read of 60 bp outlier.
    // Band must be wide enough to cover the 45-base insertion from the outlier.
    // With band_width=50 the banded result should match unbanded.
    let maj = b("CAGCAGCAGCAGCAG");
    let outlier = b("CAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAG");
    let reads = vec![maj.clone(), maj.clone(), maj.clone(), outlier];
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 50,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded.len(),
        unbanded.len(),
        "banded length mismatch with large variance"
    );
    assert_eq!(
        banded, unbanded,
        "banded result mismatch with large variance"
    );
}

// ── Semi-global alignment ─────────────────────────────────────────────────────

#[test]
fn partial_reads_semi_global() {
    // 3 full reads + 1 partial: full reads are the majority so trailing region
    // has coverage 3 >= min_cov=3, consensus is the full sequence.
    let full = b("ACGTACGTACGT");
    let partial = b("ACGTACGT");
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        ..Default::default()
    };
    let reads = vec![full.clone(), full.clone(), full.clone(), partial];
    let result = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        result.len(),
        12,
        "semi-global: got len {}, seq: '{}'",
        result.len(),
        s(&result)
    );
}

// ── Partial read coverage behaviour ──────────────────────────────────────────

#[test]
fn one_spanning_many_partial_default_min_cov_truncates() {
    // With default min_cov (≈ n/2 + 1), a single spanning read never provides
    // enough coverage to keep boundary nodes when partial reads dominate.
    // This test documents the known behaviour so a future change doesn't
    // silently alter it.
    let spanning = b("ACGTACGTACGT"); // 12 bp
    let partial = b("ACGTACGT"); //  8 bp (covers prefix only)
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        ..Default::default()
    };
    // 1 spanning + 4 partial = 5 reads; default min_cov = 5/2+1 = 3.
    // Boundary nodes only have coverage 1 (from spanning seed) → trimmed.
    let reads = vec![
        spanning.clone(),
        partial.clone(),
        partial.clone(),
        partial.clone(),
        partial.clone(),
    ];
    let result = consensus_cfg(&reads, 0, cfg);
    assert!(
        result.len() < 12,
        "expected boundary trim with default min_cov, got len {} seq '{}'",
        result.len(),
        s(&result)
    );
}

#[test]
fn one_spanning_many_partial_low_min_cov_reaches_partial_end() {
    // Lowering min_coverage_fraction removes the boundary-trim truncation,
    // but the (weight-1) normalisation in the heaviest path is a separate gate:
    // edges traversed by only one read contribute score 0, and `find` prefers
    // the shortest equal-score terminus.  So the consensus extends to the end
    // of the partial reads (position 8) but not to the spanning-only tail
    // (positions 8-11) because those edges score 0.
    //
    // To get the full-length consensus you need ≥ 2 reads covering the tail;
    // see `two_spanning_many_partial_full_length` below.
    let spanning = b("ACGTACGTACGT");
    let partial = b("ACGTACGT");
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let reads = vec![
        spanning.clone(),
        partial.clone(),
        partial.clone(),
        partial.clone(),
        partial.clone(),
    ];
    let result = consensus_cfg(&reads, 0, cfg);
    assert!(
        result.len() >= 8,
        "consensus should reach at least the partial read end, got len {} seq '{}'",
        result.len(),
        s(&result)
    );
}

#[test]
fn two_spanning_many_partial_full_length() {
    // With ≥ 2 spanning reads the tail edges (nodes 8-11) get weight ≥ 2,
    // contributing a positive score to the heaviest path.  The boundary trim
    // uses min_cov = ceil(6 * 0.1) = 1, which the spanning reads satisfy.
    // Result: full 12-bp consensus.
    let spanning = b("ACGTACGTACGT");
    let partial = b("ACGTACGT");
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let reads = vec![
        spanning.clone(),
        spanning.clone(), // 2 spanning → tail edges weight=2
        partial.clone(),
        partial.clone(),
        partial.clone(),
        partial.clone(),
    ];
    let result = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        result.len(),
        12,
        "two spanning reads should yield full-length consensus, got len {} seq '{}'",
        result.len(),
        s(&result)
    );
}

#[test]
fn overlapping_partial_reads_assemble_beyond_seed_length() {
    // Left-partial + right-partial reads that together span a longer sequence
    // than any individual read.  The seed covers only the left half; reads
    // covering the right half extend the graph via Insert ops.  With
    // min_coverage_fraction = 0.1, the assembled consensus is longer than the seed.
    //
    // Sequence: ACGTACGTACGTACGT (16 bp)
    // Left reads (12 bp):  ACGTACGTACGT
    // Right reads (12 bp): ACGTACGTACGT  (offset 4 in the full sequence → ACGTACGTACGT)
    // Together they overlap for 8 bp and cover the full 16 bp.
    let left = b("ACGTACGTACGT"); // covers positions 0-11
    let right = b("ACGTACGTACGT"); // same sequence; in a real scenario these
    // would be from a different region, but here
    // we just verify the graph can grow past seed length.
    // Use an explicit short seed so the right reads extend the graph.
    let seed = b("ACGTACGT"); // 8 bp seed
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    // 1 seed + 3 left (12 bp) + 3 right (12 bp) = 7 reads.
    // min_cov = ceil(7 * 0.1) = 1; boundary nodes survive.
    let reads = vec![
        seed.clone(),
        left.clone(),
        left.clone(),
        left.clone(),
        right.clone(),
        right.clone(),
        right.clone(),
    ];
    let result = consensus_cfg(&reads, 0, cfg);
    assert!(
        result.len() >= seed.len(),
        "assembled consensus should be at least as long as the seed, got len {} seq '{}'",
        result.len(),
        s(&result)
    );
}

#[test]
fn coverage_vec_reflects_partial_read_depth() {
    // The Consensus::coverage field must show lower values at positions that
    // only spanning reads covered and higher values where partial reads also
    // contributed.
    let spanning = b("ACGTACGTACGT"); // 12 bp
    let partial = b("ACGTACGT"); //  8 bp (covers prefix nodes)
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&spanning, cfg).unwrap();
    for _ in 0..4 {
        graph.add_read(&partial).unwrap();
    }
    let cons = graph.consensus().unwrap();
    // First 8 positions covered by all 5 reads (spanning + 4 partial).
    let prefix_min = cons.coverage[..8].iter().copied().min().unwrap_or(0);
    // Last 4 positions covered only by the spanning seed.
    let suffix_max = cons.coverage[8..].iter().copied().max().unwrap_or(0);
    assert!(
        prefix_min > suffix_max,
        "prefix coverage ({}) should exceed suffix coverage ({})",
        prefix_min,
        suffix_max
    );
}

// ── Coverage gap detection ────────────────────────────────────────────────────

#[test]
fn no_gap_when_reads_overlap() {
    // Spanning seed + partials that all overlap in the middle → no coverage gap.
    let seed = b("ACGTTGCAATGC"); // 12 bp
    let left = b("ACGTTGCA"); //  8 bp: covers positions 0-7
    let right = b("GCAATGC"); //  7 bp: covers positions 5-11 (3 bp overlap)
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&seed, cfg).unwrap();
    for _ in 0..3 {
        graph.add_read(&left).unwrap();
    }
    for _ in 0..3 {
        graph.add_read(&right).unwrap();
    }
    let cons = graph.consensus().unwrap();
    assert!(
        cons.gaps.is_empty(),
        "overlapping partials should produce no coverage gap; got {:?}",
        cons.gaps
    );
}

#[test]
fn gap_detected_when_partials_dont_overlap() {
    // Spanning seed + left partials + right partials with no overlap.
    // Seed:  ACGTTGCAATGCCCGG (16 bp)
    // Left:  ACGTT              (5 bp, covers positions 0-4)
    // Right:           CCCGG   (5 bp, covers positions 11-15)
    // Gap:         positions 5-10 (6 bp, seed-only coverage=1)
    let seed = b("ACGTTGCAATGCCCGG"); // 16 bp
    let left = b("ACGTT"); // 5 bp: unique prefix of seed
    let right = b("CCCGG"); // 5 bp: unique suffix of seed
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&seed, cfg).unwrap();
    for _ in 0..3 {
        graph.add_read(&left).unwrap();
    }
    for _ in 0..3 {
        graph.add_read(&right).unwrap();
    }
    let cons = graph.consensus().unwrap();
    assert!(
        !cons.gaps.is_empty(),
        "non-overlapping partials should produce a coverage gap; coverage={:?}",
        cons.coverage
    );
    let gap = &cons.gaps[0];
    assert!(
        gap.size() >= 6,
        "gap should span the 6 seed-only positions: {:?}",
        gap
    );
    assert_eq!(gap.start + gap.size(), gap.end);
}

#[test]
fn gap_size_is_minimum_size_estimate() {
    // Construct a scenario with a known gap width to verify size().
    // Seed: 20 bp.  Left reads cover 0-4, right reads cover 15-19.
    // The middle 10 positions (5-14) have coverage=1 (seed only).
    let seed = b("ACGTTGCAATGCCCGGTTAA"); // 20 bp
    let left = b("ACGTT"); // 5 bp: covers 0-4
    let right = b("GTTAA"); // 5 bp: covers 15-19
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&seed, cfg).unwrap();
    for _ in 0..4 {
        graph.add_read(&left).unwrap();
    }
    for _ in 0..4 {
        graph.add_read(&right).unwrap();
    }
    let cons = graph.consensus().unwrap();
    assert!(
        !cons.gaps.is_empty(),
        "expected a gap; coverage: {:?}",
        cons.coverage
    );
    // The gap should span the seed-only middle region (at least 10 bp).
    let total_gap: usize = cons.gaps.iter().map(|g| g.size()).sum();
    assert!(
        total_gap >= 10,
        "expected gap ≥ 10 bp, got {total_gap}; gaps: {:?}",
        cons.gaps
    );
}

#[test]
fn single_read_has_no_gap() {
    // With only the seed read, coverage is all 1s but there are no
    // well-supported flanks, so detect_coverage_gaps returns empty.
    let cfg = PoaConfig {
        min_reads: 1,
        ..Default::default()
    };
    let graph = PoaGraph::new(b"ACGTTGCAATGC", cfg).unwrap();
    let cons = graph.consensus().unwrap();
    assert!(
        cons.gaps.is_empty(),
        "single-read consensus should have no gaps"
    );
}

#[test]
fn gap_kind_spanning_for_seed_based_gap() {
    // A seed-based gap must have kind=Spanning so callers know a minimum
    // size estimate is available via size().
    let seed = b("ACGTTGCAATGCCCGG");
    let left = b("ACGTT");
    let right = b("CCCGG");
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&seed, cfg).unwrap();
    for _ in 0..3 {
        graph.add_read(&left).unwrap();
    }
    for _ in 0..3 {
        graph.add_read(&right).unwrap();
    }
    let cons = graph.consensus().unwrap();
    assert!(!cons.gaps.is_empty());
    assert_eq!(
        cons.gaps[0].kind,
        poa_consensus::GapKind::Spanning,
        "seed-based gaps must be Spanning"
    );
    assert_eq!(cons.gaps[0].min_size(), Some(cons.gaps[0].size()));
}

#[test]
fn bridged_consensus_unknown_gap() {
    // Two completely disjoint read groups — left reads, then right reads —
    // with no read spanning the middle.  bridged_consensus should produce a
    // single Consensus whose gaps contain exactly one Unknown gap at the join.
    let left_reads: Vec<Vec<u8>> = (0..4).map(|_| b("ACGTTGCA")).collect();
    let right_reads: Vec<Vec<u8>> = (0..4).map(|_| b("ATGCCCGG")).collect();
    let left_refs: Vec<&[u8]> = left_reads.iter().map(|r| r.as_slice()).collect();
    let right_refs: Vec<&[u8]> = right_reads.iter().map(|r| r.as_slice()).collect();
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        ..Default::default()
    };
    let cons = poa_consensus::bridged_consensus(&left_refs, 0, &right_refs, 0, &cfg).unwrap();

    // Sequence is the concatenation of both consensuses.
    assert!(!cons.sequence.is_empty());

    // Exactly one Unknown gap at the join point.
    let unknown: Vec<_> = cons
        .gaps
        .iter()
        .filter(|g| g.kind == poa_consensus::GapKind::Unknown)
        .collect();
    assert_eq!(unknown.len(), 1, "expected exactly one Unknown gap");

    let gap = unknown[0];
    assert_eq!(
        gap.start, gap.end,
        "Unknown gap should be an insertion point (start==end)"
    );
    assert_eq!(gap.min_size(), None, "Unknown gap has no minimum size");

    // Total minimum size: at least as long as the two consensus segments.
    assert!(cons.sequence.len() > 0);
    assert_eq!(cons.n_reads, 8);
}

// ── path_weights and weight_fraction ─────────────────────────────────────────

#[test]
fn path_weights_reflect_edge_support() {
    // 1 spanning seed + 4 partial reads covering the first 8 of 12 nodes.
    // Interior edges (0-7) should have weight 5 (seed + 4 partial).
    // The partial reads end at node 7, so the consensus (heaviest path) stops
    // there.  All 8 weights should be ≥ 2 (shared) and n_reads should be 5.
    let spanning = b("ACGTACGTACGT");
    let partial = b("ACGTACGT");
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&spanning, cfg).unwrap();
    for _ in 0..4 {
        graph.add_read(&partial).unwrap();
    }
    let cons = graph.consensus().unwrap();

    assert_eq!(cons.n_reads, 5);
    assert_eq!(cons.path_weights.len(), cons.sequence.len());
    // Every weight in the consensus should reflect multi-read support.
    for (i, &w) in cons.path_weights.iter().enumerate() {
        assert!(
            w >= 2,
            "position {i}: weight {w} should be ≥ 2 (shared by seed + partial)"
        );
    }
}

#[test]
fn weight_fraction_in_unit_interval() {
    let reads = vec![b("ACGTACGT"); 5];
    let cons = consensus_cfg(&reads, 0, PoaConfig::default());
    // Build Consensus directly to check the fraction helper.
    let mut graph = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
    for r in &reads[1..] {
        graph.add_read(r).unwrap();
    }
    let c = graph.consensus().unwrap();
    let fracs = c.weight_fraction();
    assert_eq!(fracs.len(), cons.len());
    for (i, &f) in fracs.iter().enumerate() {
        assert!(
            (0.0..=1.0).contains(&f),
            "position {i}: fraction {f} out of [0,1]"
        );
    }
    // All reads identical → all fractions should be 1.0.
    for (i, &f) in fracs.iter().enumerate() {
        assert!(
            (f - 1.0).abs() < 1e-6,
            "position {i}: expected fraction 1.0, got {f}"
        );
    }
}

#[test]
fn weight_fraction_drops_at_single_read_positions() {
    // Non-repetitive spanning sequence so partial reads have exactly one valid
    // alignment position (avoids the rotation-phase tie-break issue that arises
    // with periodic sequences like ACGTACGTACGT).
    //
    // 2 spanning reads (needed so tail edges score > 0 and enter the path) +
    // 4 partial reads covering only the first 8 of 12 bases.
    // Tail positions (8-11) are supported only by the 2 spanning reads;
    // their fraction (2/6 ≈ 0.33) should be below the prefix fraction (6/6 = 1.0).
    let spanning = b("ACGTTGCAATGC"); // 12 bp, no 8-mer repeats
    let partial = b("ACGTTGCA"); //  8 bp, uniquely matches positions 0-7
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        min_coverage_fraction: 0.1,
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&spanning, cfg.clone()).unwrap();
    graph.add_read(&spanning).unwrap();
    for _ in 0..4 {
        graph.add_read(&partial).unwrap();
    }
    let cons = graph.consensus().unwrap();
    let fracs = cons.weight_fraction();

    assert_eq!(cons.sequence.len(), 12, "expected full-length consensus");
    let prefix_frac: f32 = fracs[..8].iter().copied().sum::<f32>() / 8.0;
    let suffix_frac: f32 = fracs[8..].iter().copied().sum::<f32>() / 4.0;
    assert!(
        prefix_frac > suffix_frac,
        "prefix avg fraction ({prefix_frac:.2}) should exceed suffix ({suffix_frac:.2})"
    );
}

// Semi-global op-level tests: verify the alignment ops themselves, not just
// the consensus.  These catch bugs that happen to not affect the output length
// but still corrupt edge weights or delete_counts.

#[test]
fn semi_global_no_prefix_deletes_for_mid_start_read() {
    // Seed "GGACGT", partial read "ACGT" matches the suffix perfectly.
    // In global mode the aligner is forced to start from the source (G,G) and
    // emits Delete ops for those prefix nodes.  Semi-global lets the read start
    // at the first matching node and must produce zero Deletes.
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        band_width: 0,
        ..Default::default()
    };
    let graph = PoaGraph::new(b"GGACGT", cfg).unwrap();
    let (ops, _) = graph.align_read_ops_unbanded(b"ACGT").unwrap();
    let n_del = ops
        .iter()
        .filter(|op| matches!(op, AlignOp::Delete(_)))
        .count();
    let n_mat = ops
        .iter()
        .filter(|op| matches!(op, AlignOp::Match(_)))
        .count();
    assert_eq!(
        n_del, 0,
        "semi-global: expected no prefix Deletes, got {:?}",
        ops
    );
    assert_eq!(n_mat, 4, "semi-global: expected 4 Matches, got {:?}", ops);
}

#[test]
fn global_produces_prefix_deletes_for_mid_start_read() {
    // Same setup, global mode: the alignment is forced through the GG prefix,
    // producing Delete ops for those nodes.
    let cfg = PoaConfig {
        band_width: 0,
        ..Default::default()
    };
    let graph = PoaGraph::new(b"GGACGT", cfg).unwrap();
    let (ops, _) = graph.align_read_ops_unbanded(b"ACGT").unwrap();
    let n_del = ops
        .iter()
        .filter(|op| matches!(op, AlignOp::Delete(_)))
        .count();
    assert!(
        n_del > 0,
        "global: expected prefix Delete ops for mid-start read"
    );
}

#[test]
fn semi_global_spanning_read_matches_global() {
    // A read that spans the full seed is not partial; semi-global and global
    // must produce the same consensus.
    let reads = vec![b("ACGTACGT"); 4];
    let global = consensus(&reads, 0);
    let cfg = PoaConfig {
        alignment_mode: AlignmentMode::SemiGlobal,
        ..Default::default()
    };
    let semi = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        global, semi,
        "spanning reads: semi-global must equal global"
    );
}

// ── Reverse complement / orientation ─────────────────────────────────────────

#[test]
fn reverse_complement_basic() {
    use crate::reverse_complement;
    assert_eq!(reverse_complement(b"ACGT"), b"ACGT");
    assert_eq!(reverse_complement(b"AAAA"), b"TTTT");
    assert_eq!(reverse_complement(b"GCTA"), b"TAGC");
}

#[test]
fn orient_to_seed_forward() {
    use crate::Strand;
    use crate::orient_to_seed;
    let seed = b("ACGTACGTACGT");
    let read = b("ACGTACGT");
    assert_eq!(orient_to_seed(&read, &seed, 4), Strand::Forward);
}

#[test]
fn orient_to_seed_reverse() {
    use crate::Strand;
    use crate::orient_to_seed;
    // Use a non-palindromic sequence so forward and RC share no k-mers.
    let seed = b("AAAACCCCGGGG");
    let rc = crate::reverse_complement(&seed);
    assert_eq!(orient_to_seed(&rc, &seed, 4), Strand::Reverse);
}

#[test]
fn mixed_strand_input() {
    use crate::auto_orient;
    let seed = b("CATCATCAT");
    let rc = crate::reverse_complement(&seed);
    let reads = vec![seed.clone(), seed.clone(), rc.clone(), rc.clone()];
    let oriented: Vec<Vec<u8>> = auto_orient(&reads, 0)
        .into_iter()
        .map(|c| c.into_owned())
        .collect();
    // All oriented reads should match or be rc'd to match seed strand.
    let mut graph = PoaGraph::new(&oriented[0], PoaConfig::default()).unwrap();
    for r in &oriented[1..] {
        graph.add_read(r).unwrap();
    }
    let result = graph.consensus().unwrap().sequence;
    assert_eq!(result.len(), 9, "mixed strand: got len {}", result.len());
}

// ── Majority-frequency consensus ──────────────────────────────────────────────

fn mf_cfg() -> PoaConfig {
    PoaConfig {
        consensus_mode: ConsensusMode::MajorityFrequency,
        ..Default::default()
    }
}

fn mf_consensus(reads: &[Vec<u8>], seed_idx: usize) -> Vec<u8> {
    consensus_cfg(reads, seed_idx, mf_cfg())
}

#[test]
fn mf_identical_reads() {
    let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCAT")];
    assert_eq!(mf_consensus(&reads, 0), b("CATCATCAT"));
}

#[test]
fn mf_matches_hb_on_clean_input() {
    // On a read set with no noise, MF and HB should agree.
    let reads = vec![
        b("CAGCAGCAG"),
        b("CAGCAGCAG"),
        b("CAGCAGCAGCAG"),
        b("CAGCAGCAG"),
    ];
    let hb = consensus(&reads, 0);
    let mf = mf_consensus(&reads, 0);
    assert_eq!(hb, mf, "HB and MF disagree on clean input");
}

#[test]
fn mf_boundary_trim_leading() {
    // Seed has 3 extra leading bases not present in the majority.
    // MF should exclude them because gap votes outnumber base votes.
    let reads = vec![
        b("XXXCATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let result = s(&mf_consensus(&reads, 0));
    assert_eq!(result, "CATCATCAT", "got: {}", result);
}

#[test]
fn mf_boundary_trim_trailing() {
    let reads = vec![
        b("CATCATCATXXX"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let result = s(&mf_consensus(&reads, 0));
    assert_eq!(result, "CATCATCAT", "got: {}", result);
}

#[test]
fn mf_majority_base_wins() {
    // 3 reads have CAT, 1 has CGT at position 1. MF should pick A.
    let reads = vec![
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CGTCATCAT"),
    ];
    assert_eq!(s(&mf_consensus(&reads, 0)), "CATCATCAT");
}

#[test]
fn mf_single_outlier_not_inflated() {
    // One read has an extra CAT; MF should not include it.
    let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCATCAT")];
    assert_eq!(mf_consensus(&reads, 0).len(), 9);
}

// ── GraphStats ────────────────────────────────────────────────────────────────

fn build_graph(reads: &[Vec<u8>], seed_idx: usize) -> PoaGraph {
    let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
    for (i, read) in reads.iter().enumerate() {
        if i != seed_idx {
            graph.add_read(read).unwrap();
        }
    }
    graph
}

#[test]
fn stats_clean_linear_no_bubbles() {
    // Identical reads produce a clean linear graph with no bubbles.
    let reads = vec![
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let st = build_graph(&reads, 0).stats();
    assert_eq!(st.bubble_count, 0);
    assert_eq!(st.max_bubble_depth, 0);
    assert_eq!(st.node_count, 9);
    // All reads match every node: delete_count=0 everywhere → entropy=0.
    assert_eq!(st.mean_column_entropy, 0.0);
}

#[test]
fn stats_bubble_detected() {
    // 3 reads with CATCATCAT, 1 with CGTCATCAT → SNV bubble at position 1 (A vs G).
    let reads = vec![
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CGTCATCAT"),
    ];
    let st = build_graph(&reads, 0).stats();
    assert_eq!(st.bubble_count, 1, "expected 1 bubble");
    // Minority arm has 1 read (the CGT read creates a G branch at position 1).
    assert_eq!(st.max_bubble_depth, 1, "minority arm weight should be 1");
}

#[test]
fn stats_entropy_nonzero_on_length_variation() {
    // Seed has 3 leading X nodes that other reads delete.
    // The X nodes get delete_count=3, coverage=1 → entropy > 0.
    // (Shorter reads aligned to a longer graph in global mode don't generate trailing
    // deletes — they simply stop at the best-scoring diagonal. Leading deletes DO fire
    // because the traceback reaches t=0, j=0 via the D-chain at the j=0 column.)
    let reads = vec![
        b("XXXCATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let st = build_graph(&reads, 0).stats();
    // X nodes: coverage=1, delete_count=3 → p=0.25, binary_entropy(0.25) ≈ 0.811 bits.
    assert!(
        st.mean_column_entropy > 0.0,
        "expected nonzero entropy, got {}",
        st.mean_column_entropy
    );
}

#[test]
fn stats_node_edge_counts() {
    // 2 reads with length variation: 9-node backbone + 3 extra nodes from longer read.
    // The longer read aligns as Insert(CAT) + Match(nodes 0-8), creating nodes 9,10,11
    // with edges 9→10, 10→11, 11→0. Backbone edges 0→1..7→8 already existed = 8.
    let reads = vec![b("CATCATCAT"), b("CATCATCATCAT")];
    let st = build_graph(&reads, 0).stats();
    assert_eq!(st.node_count, 12, "9 + 3 extra nodes");
    // 8 backbone + 3 new edges (9→10, 10→11, 11→0) = 11 total.
    assert_eq!(st.edge_count, 11);
}

#[test]
fn stats_coverage_mean_uniform() {
    // All 4 reads match all 9 nodes → coverage=4 everywhere → variance=0.
    let reads = vec![
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
        b("CATCATCAT"),
    ];
    let st = build_graph(&reads, 0).stats();
    assert!((st.coverage_mean - 4.0).abs() < 1e-10);
    assert!(st.coverage_variance < 1e-10);
}

// ── Multi-allele consensus ────────────────────────────────────────────────────

fn multi_graph(reads: &[Vec<u8>], seed_idx: usize) -> PoaGraph {
    let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
    for (i, read) in reads.iter().enumerate() {
        if i != seed_idx {
            graph.add_read(read).unwrap();
        }
    }
    graph
}

#[test]
fn consensus_multi_single_allele() {
    // Identical reads → no bubble → consensus_multi falls through to single consensus.
    let reads = vec![b("CATCATCAT"); 4];
    let g = multi_graph(&reads, 0);
    let results = g.consensus_multi().unwrap();
    assert_eq!(results.len(), 1, "expected 1 allele for homozygous input");
    assert_eq!(results[0].sequence, b("CATCATCAT"));
}

#[test]
fn consensus_multi_snv_bubble() {
    // 4 reads with allele A (CATCATCAT) and 4 with allele B (CATCGTCAT).
    // A SNV at position 4 (A→G) creates a 2-arm bubble.
    let allele_a = b("CATCATCAT");
    let allele_b = b("CATCGTCAT");
    let reads: Vec<Vec<u8>> = (0..4)
        .map(|_| allele_a.clone())
        .chain((0..4).map(|_| allele_b.clone()))
        .collect();
    let g = multi_graph(&reads, 0);
    let results = g.consensus_multi().unwrap();
    assert_eq!(results.len(), 2, "expected 2 alleles for SNV input");
    let seqs: Vec<String> = results.iter().map(|c| s(&c.sequence)).collect();
    assert!(
        seqs.iter().any(|seq| seq == "CATCATCAT"),
        "missing CATCATCAT allele: {:?}",
        seqs
    );
    assert!(
        seqs.iter().any(|seq| seq == "CATCGTCAT"),
        "missing CATCGTCAT allele: {:?}",
        seqs
    );
}

#[test]
fn consensus_multi_length_variation() {
    // Two alleles with different repeat counts flanked by matching anchors.
    // The anchor regions create a proper bubble between the two allele lengths.
    // short: AAA + 2×CAT + TTTTTT = 14 bp
    // long : AAA + 3×CAT + TTTTTT = 17 bp
    let short = b("AAACATCATTTTTT");
    let long_ = b("AAACATCATCATTTTTT");
    let reads: Vec<Vec<u8>> = (0..4)
        .map(|_| short.clone())
        .chain((0..4).map(|_| long_.clone()))
        .collect();
    let g = multi_graph(&reads, 0);
    let results = g.consensus_multi().unwrap();
    assert_eq!(
        results.len(),
        2,
        "expected 2 alleles for length-variation input"
    );
    let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
    assert!(lens.contains(&14), "expected 14-bp allele; got {:?}", lens);
    assert!(lens.contains(&17), "expected 17-bp allele; got {:?}", lens);
}

#[test]
fn consensus_multi_insufficient_depth_per_allele() {
    // 4 reads total, min_reads=3. Each allele only gets 2 → InsufficientDepth.
    let allele_a = b("CATCATCAT");
    let allele_b = b("CATCGTCAT");
    let reads = vec![allele_a.clone(), allele_a, allele_b.clone(), allele_b];
    let cfg = PoaConfig {
        min_reads: 3,
        ..Default::default()
    };
    let mut g = PoaGraph::new(&reads[0], cfg).unwrap();
    for r in &reads[1..] {
        g.add_read(r).unwrap();
    }
    let result = g.consensus_multi();
    assert!(
        matches!(result, Err(PoaError::InsufficientDepth { .. })),
        "expected InsufficientDepth, got {:?}",
        result.map(|v| v.len())
    );
}

// ── Longer-sequence stress tests ──────────────────────────────────────────────

#[test]
fn long_repeat_consensus_correctness() {
    let seq: Vec<u8> = "CAT".repeat(30).into_bytes(); // 90 bp
    let reads = vec![seq.clone(); 6];
    assert_eq!(consensus(&reads, 0), seq, "30×CAT consensus mismatch");
}

#[test]
fn long_repeat_length_majority_wins() {
    // 8 reads at 60 bp (20×CAT), 2 outliers at 63 bp (21×CAT)
    let maj: Vec<u8> = "CAT".repeat(20).into_bytes();
    let out: Vec<u8> = "CAT".repeat(21).into_bytes();
    let mut reads: Vec<Vec<u8>> = vec![maj.clone(); 8];
    reads.extend(vec![out; 2]);
    let result = consensus(&reads, 0);
    assert_eq!(
        result.len(),
        60,
        "expected 60-bp majority, got {} bp",
        result.len()
    );
}

#[test]
fn long_repeat_snv_correction() {
    // 9 correct reads + 1 noisy read with a single mismatch at position 30
    let correct: Vec<u8> = "CAT".repeat(20).into_bytes(); // 60 bp
    let mut noisy = correct.clone();
    noisy[30] = b'G';
    let mut reads: Vec<Vec<u8>> = vec![correct.clone(); 9];
    reads.push(noisy);
    let result = consensus(&reads, 0);
    assert_eq!(
        result, correct,
        "SNV from single noisy read should not affect consensus"
    );
}

#[test]
fn long_banded_matches_unbanded() {
    // 4 × 72 bp + 1 × 78 bp (length outlier), band=30
    let base: Vec<u8> = "CAT".repeat(24).into_bytes(); // 72 bp
    let long: Vec<u8> = "CAT".repeat(26).into_bytes(); // 78 bp
    let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
    reads.push(long);
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 30,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded, unbanded,
        "banded(30) should match unbanded for small divergence"
    );
}

#[test]
fn long_adaptive_band_matches_unbanded() {
    // 4 × 72 bp + 1 × 81 bp, adaptive band b=10 f=0.05
    let base: Vec<u8> = "CAT".repeat(24).into_bytes(); // 72 bp
    let long: Vec<u8> = "CAT".repeat(27).into_bytes(); // 81 bp
    let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
    reads.push(long);
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        adaptive_band: true,
        adaptive_band_b: 10,
        adaptive_band_f: 0.05,
        ..Default::default()
    };
    let adaptive = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        adaptive, unbanded,
        "adaptive band should match unbanded for small divergence"
    );
}

// ── Performance optimisation tests ───────────────────────────────────────────

#[test]
fn skip_fires_on_clean_reads() {
    // Zero-error reads: the diagonal skip fires on ~100% of rows.
    // Verify correctness — narrow band with skip active must produce the same
    // consensus as unbanded and must equal the read itself.
    let read: Vec<u8> = "CAT".repeat(10).into_bytes(); // 30 bp, 10 identical reads
    let reads: Vec<Vec<u8>> = vec![read.clone(); 10];
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 5,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded, unbanded,
        "diagonal skip: banded must match unbanded on identical reads"
    );
    assert_eq!(
        banded, read,
        "diagonal skip: consensus of identical reads must equal the read"
    );
}

#[test]
fn tracking_band_survives_phase_shift() {
    // One read with a 10-bp prefix insertion shifts the alignment diagonal by 10.
    // With a fixed-diagonal band of 5 this would be BandTooNarrow; with the
    // tracking band the band re-centres after the shift, and smart retry widens
    // the initial band so alignment succeeds.
    let base: Vec<u8> = "ACGT".repeat(15).into_bytes(); // 60 bp
    let shifted: Vec<u8> = {
        let mut s = b"AAAAAAAAAA".to_vec(); // 10 bp prefix → diagonal drift +10
        s.extend_from_slice(&base);
        s
    }; // 70 bp
    let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
    reads.push(shifted);
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 5,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded, unbanded,
        "tracking band: must match unbanded on reads with a large phase shift"
    );
}

#[test]
fn sv_retry_correct() {
    // One outlier read with a 15-bp expansion against 5 short reads. band_width=3
    // is too narrow to track the diagonal shift: approaching-edge fires on every
    // row (right_margin = w = 3 < GAP_MARGIN), forcing a smart retry with a wider
    // band (~17). The retry band covers j=l=24 at the last graph node, so the
    // alignment succeeds without a second retry. Majority consensus is "CAT"×3.
    let short: Vec<u8> = "CAT".repeat(3).into_bytes(); // 9 bp
    let expanded: Vec<u8> = "CAT".repeat(8).into_bytes(); // 24 bp (+15 bp expansion)
    let mut reads: Vec<Vec<u8>> = vec![short.clone(); 5];
    reads.push(expanded);
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 3,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded, unbanded,
        "sv_retry: smart retry must produce correct consensus when SV read forces band widening"
    );
}

#[test]
fn consensus_multi_long_flanked_str() {
    // Two alleles anchored by GGGGG / AAAAA flanks:
    //   allele_a: GGGGG + 8×CAT + AAAAA  = 5+24+5 = 34 bp
    //   allele_b: GGGGG + 11×CAT + AAAAA = 5+33+5 = 43 bp
    let flank_l = b("GGGGG");
    let flank_r = b("AAAAA");
    let inner_a: Vec<u8> = "CAT".repeat(8).into_bytes();
    let inner_b: Vec<u8> = "CAT".repeat(11).into_bytes();
    let allele_a: Vec<u8> = [flank_l.as_slice(), inner_a.as_slice(), flank_r.as_slice()].concat();
    let allele_b: Vec<u8> = [flank_l.as_slice(), inner_b.as_slice(), flank_r.as_slice()].concat();
    let mut reads: Vec<Vec<u8>> = vec![allele_a.clone(); 5];
    reads.extend(vec![allele_b.clone(); 5]);
    let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
    for r in &reads[1..] {
        g.add_read(r).unwrap();
    }
    let results = g.consensus_multi().unwrap();
    let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
    assert_eq!(results.len(), 2, "expected 2 alleles; got {:?}", lens);
    assert!(lens.contains(&34), "expected 34-bp allele; got {:?}", lens);
    assert!(lens.contains(&43), "expected 43-bp allele; got {:?}", lens);
}

#[test]
fn consensus_multi_snv_in_long_context() {
    // SNV at position 10 in a 41-bp read: allele_a has 'A' at pos 10, allele_b is all-T
    let a: Vec<u8> = {
        let mut v = vec![b'T'; 41];
        v[10] = b'A';
        v
    };
    let bv: Vec<u8> = vec![b'T'; 41];
    let mut reads: Vec<Vec<u8>> = vec![a.clone(); 5];
    reads.extend(vec![bv.clone(); 5]);
    let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
    for r in &reads[1..] {
        g.add_read(r).unwrap();
    }
    let results = g.consensus_multi().unwrap();
    assert_eq!(
        results.len(),
        2,
        "expected 2 alleles for SNV; got {}",
        results.len()
    );
    let seqs: Vec<Vec<u8>> = results.into_iter().map(|c| c.sequence).collect();
    assert!(
        seqs.iter().any(|s| s == &a),
        "allele_a (A at pos 10) not found in results"
    );
    assert!(
        seqs.iter().any(|s| s == &bv),
        "allele_b (all-T) not found in results"
    );
}

#[test]
fn consensus_multi_skewed_allele_ratio() {
    // 7:3 ratio — minor allele at 30% detected with default min_allele_freq=0.25
    // allele_a: TTTT + 6×CAT + CCCC = 4+18+4 = 26 bp
    // allele_b: TTTT + 10×CAT + CCCC = 4+30+4 = 38 bp
    let flank_l = b("TTTT");
    let flank_r = b("CCCC");
    let inner_a: Vec<u8> = "CAT".repeat(6).into_bytes();
    let inner_b: Vec<u8> = "CAT".repeat(10).into_bytes();
    let allele_a: Vec<u8> = [flank_l.as_slice(), inner_a.as_slice(), flank_r.as_slice()].concat();
    let allele_b: Vec<u8> = [flank_l.as_slice(), inner_b.as_slice(), flank_r.as_slice()].concat();
    let mut reads: Vec<Vec<u8>> = vec![allele_a.clone(); 7];
    reads.extend(vec![allele_b.clone(); 3]);
    let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
    for r in &reads[1..] {
        g.add_read(r).unwrap();
    }
    let results = g.consensus_multi().unwrap();
    let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
    assert_eq!(
        results.len(),
        2,
        "expected 2 alleles at 7:3; got {:?}",
        lens
    );
    assert!(
        lens.contains(&26),
        "expected 26-bp majority allele; got {:?}",
        lens
    );
    assert!(
        lens.contains(&38),
        "expected 38-bp minor allele; got {:?}",
        lens
    );
}

#[test]
fn long_reads_noise_and_banding() {
    // 63-bp majority + scattered single-base errors; banded with band=40
    let correct: Vec<u8> = "CAT".repeat(21).into_bytes(); // 63 bp
    let mut err1 = correct.clone();
    err1[20] = b'G';
    let mut err2 = correct.clone();
    err2[45] = b'T';
    let mut reads: Vec<Vec<u8>> = vec![correct.clone(); 6];
    reads.extend(vec![err1; 2]);
    reads.extend(vec![err2; 2]);
    let cfg = PoaConfig {
        band_width: 40,
        ..Default::default()
    };
    let result = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        result, correct,
        "banded consensus should correct isolated noise in 63-bp reads"
    );
}

// ── Non-repeat longer-sequence tests ─────────────────────────────────────────
//
// BASE_60: ATCGATCGTT ACGATCGTAG CTAGTCATGC TAATCGTAGC GATCGTAACG ATCGATCGTA
// 60 bp, mixed composition, no periodic structure.

const BASE_60: &[u8] = b"ATCGATCGTTACGATCGTAGCTAGTCATGCTAATCGTAGCGATCGTAACGATCGATCGTA";

#[test]
fn long_nonrepeat_consensus_correctness() {
    let reads = vec![BASE_60.to_vec(); 6];
    assert_eq!(
        consensus(&reads, 0),
        BASE_60,
        "6 identical non-repeat reads"
    );
}

#[test]
fn long_nonrepeat_snv_correction() {
    // 9 correct + 1 with T→G at position 30
    let mut noisy = BASE_60.to_vec();
    noisy[30] = b'G';
    let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 9];
    reads.push(noisy);
    assert_eq!(
        consensus(&reads, 0),
        BASE_60,
        "single noisy read must not flip consensus base"
    );
}

#[test]
fn long_nonrepeat_banded_matches_unbanded() {
    // 4 × 60 bp + 1 × 66 bp (6-bp insertion at position 30), band=30
    let mut long = BASE_60.to_vec();
    long.splice(30..30, *b"GCTAGC");
    assert_eq!(long.len(), 66);
    let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 4];
    reads.push(long);
    let unbanded = consensus(&reads, 0);
    let cfg = PoaConfig {
        band_width: 30,
        ..Default::default()
    };
    let banded = consensus_cfg(&reads, 0, cfg);
    assert_eq!(
        banded, unbanded,
        "banded(30) should match unbanded on non-repeat sequence"
    );
}

#[test]
fn consensus_multi_nonrepeat_snv() {
    // Two alleles that differ only at position 30 (T vs G) in a non-repeat context
    let mut allele_b = BASE_60.to_vec();
    allele_b[30] = b'G';
    let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 5];
    reads.extend(vec![allele_b.clone(); 5]);
    let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
    for r in &reads[1..] {
        g.add_read(r).unwrap();
    }
    let results = g.consensus_multi().unwrap();
    assert_eq!(results.len(), 2, "expected 2 alleles for non-repeat SNV");
    let seqs: Vec<Vec<u8>> = results.into_iter().map(|c| c.sequence).collect();
    assert!(
        seqs.iter().any(|s| s.as_slice() == BASE_60),
        "allele_a not recovered"
    );
    assert!(
        seqs.iter().any(|s| s == &allele_b),
        "allele_b not recovered"
    );
}

// ── Functional convenience wrappers ───────────────────────────────────────────

#[test]
fn fn_consensus_basic() {
    let reads: Vec<&[u8]> = vec![b"CATCATCAT", b"CATCATCAT", b"CATCATCAT"];
    let result = poa_consensus::consensus(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(result.sequence, b"CATCATCAT");
}

#[test]
fn fn_consensus_empty_input() {
    let result = poa_consensus::consensus(&[], 0, &PoaConfig::default());
    assert!(matches!(result, Err(PoaError::EmptyInput)));
}

#[test]
fn fn_consensus_seed_out_of_bounds() {
    let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
    let result = poa_consensus::consensus(&reads, 5, &PoaConfig::default());
    assert!(matches!(
        result,
        Err(PoaError::SeedOutOfBounds { index: 5, len: 2 })
    ));
}

#[test]
fn fn_consensus_config_respected() {
    // min_reads=5 with only 3 reads → InsufficientDepth
    let reads: Vec<&[u8]> = vec![b"CATCATCAT", b"CATCATCAT", b"CATCATCAT"];
    let cfg = PoaConfig {
        min_reads: 5,
        ..Default::default()
    };
    let result = poa_consensus::consensus(&reads, 0, &cfg);
    assert!(matches!(result, Err(PoaError::InsufficientDepth { .. })));
}

#[test]
fn fn_consensus_multi_two_alleles() {
    let allele_a: &[u8] = b"CATCATCAT";
    let allele_b: &[u8] = b"CATCGTCAT";
    let reads: Vec<&[u8]> = vec![
        allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b, allele_b,
    ];
    let results = poa_consensus::consensus_multi(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(results.len(), 2, "expected 2 alleles");
}

#[test]
fn fn_consensus_multi_empty_input() {
    let result = poa_consensus::consensus_multi(&[], 0, &PoaConfig::default());
    assert!(matches!(result, Err(PoaError::EmptyInput)));
}

#[test]
fn fn_consensus_multi_seed_out_of_bounds() {
    let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
    let result = poa_consensus::consensus_multi(&reads, 99, &PoaConfig::default());
    assert!(matches!(
        result,
        Err(PoaError::SeedOutOfBounds { index: 99, len: 2 })
    ));
}

// ── Two-pass adaptive mode ─────────────────────────────────────────────────────

#[test]
fn adaptive_clean_single_allele() {
    // No bubbles → single consensus, no pass 2.
    let reads: Vec<&[u8]> = vec![b"CATCATCAT"; 6];
    let results = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(results.len(), 1);
    assert_eq!(results[0].sequence, b"CATCATCAT");
}

#[test]
fn adaptive_snv_bubble_splits_alleles() {
    // Clear SNV bubble → multi-allele path, two alleles returned.
    let a: &[u8] = b"CATCATCAT";
    let bv: &[u8] = b"CATCGTCAT";
    let reads: Vec<&[u8]> = vec![a, a, a, a, bv, bv, bv, bv];
    let results = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(results.len(), 2, "expected two alleles from SNV bubble");
    let seqs: Vec<&[u8]> = results.iter().map(|c| c.sequence.as_slice()).collect();
    assert!(seqs.contains(&a), "allele_a not in results");
    assert!(seqs.contains(&bv), "allele_b not in results");
}

#[test]
fn adaptive_noisy_tightens_coverage() {
    // 4 correct reads + 4 reads each with a unique error → singleton fraction high.
    // Adaptive mode should tighten coverage and return a single clean consensus.
    let correct: Vec<u8> = b("CATCATCATCATCAT");
    let mut r1 = correct.clone();
    r1[0] = b'G';
    let mut r2 = correct.clone();
    r2[3] = b'G';
    let mut r3 = correct.clone();
    r3[6] = b'G';
    let mut r4 = correct.clone();
    r4[9] = b'G';
    let reads: Vec<&[u8]> = vec![&correct, &correct, &correct, &correct, &r1, &r2, &r3, &r4];
    let results = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(results.len(), 1);
    assert_eq!(
        results[0].sequence, correct,
        "noisy reads should be filtered"
    );
}

#[test]
fn adaptive_partial_reads_switches_semi_global() {
    // Seed extends well beyond several partial reads → high coverage CV.
    // Adaptive mode should switch to semi-global and return a correct consensus.
    let full: Vec<u8> = b("GGGCATCATCATCATAAA");
    let partial: Vec<u8> = b("CATCATCAT");
    let reads: Vec<&[u8]> = vec![
        full.as_slice(),
        full.as_slice(),
        full.as_slice(),
        partial.as_slice(),
        partial.as_slice(),
        partial.as_slice(),
    ];
    // Just verify it completes without error and returns a single result.
    let results = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
    assert_eq!(results.len(), 1);
}

#[test]
fn adaptive_empty_input() {
    let result = poa_consensus::consensus_adaptive(&[], 0, &PoaConfig::default());
    assert!(matches!(result, Err(PoaError::EmptyInput)));
}

#[test]
fn adaptive_seed_out_of_bounds() {
    let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
    let result = poa_consensus::consensus_adaptive(&reads, 9, &PoaConfig::default());
    assert!(matches!(
        result,
        Err(PoaError::SeedOutOfBounds { index: 9, len: 2 })
    ));
}

// ── Remaining TODO tests ───────────────────────────────────────────────────────

#[test]
fn reads_long_aligns_correctly() {
    // A long read should align successfully and produce a correct consensus.
    let long: Vec<u8> = b"A".repeat(200);
    let cfg = PoaConfig {
        ..Default::default()
    };
    let mut graph = PoaGraph::new(&long, cfg).unwrap();
    graph.add_read(&long).unwrap();
    // The warning counter is always 0 with the new aligner (no band warnings).
    assert_eq!(
        graph.warnings_emitted(),
        0,
        "no warnings expected with new aligner"
    );
}

#[test]
fn multi_allele_low_per_allele_depth() {
    // Total reads (7) exceeds min_reads (4), but the minor allele group (3 reads)
    // does not — verifying that depth is checked per group, not on the total.
    let allele_a: &[u8] = b"CATCATCAT";
    let allele_b: &[u8] = b"CATCGTCAT";
    let cfg = PoaConfig {
        min_reads: 4,
        ..Default::default()
    };
    // 4 reads of allele_a, 3 reads of allele_b → total 7 >= min_reads 4, minor group 3 < 4
    let reads: Vec<&[u8]> = vec![
        allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b,
    ];
    let result = poa_consensus::consensus_multi(&reads, 0, &cfg);
    assert!(
        matches!(result, Err(PoaError::InsufficientDepth { .. })),
        "expected InsufficientDepth for minor allele group, got {:?}",
        result.map(|v| v.len())
    );
}

// ─── Structural bubble phasing ───────────────────────────────────────────────

/// Flanked length variants create a true structural bubble (reconvergence at the
/// right flank). The phasing should detect it and return one consensus per allele.
#[test]
fn structural_bubble_phasing_splits_flanked_length_variants() {
    let left = b"ACGTACGT";
    let right = b"TTTTGGGG";
    let short_mid: Vec<u8> = b"CAT".repeat(5); // 15 bp
    let long_mid: Vec<u8> = b"CAT".repeat(10); // 30 bp

    let mut short_read = left.to_vec();
    short_read.extend_from_slice(&short_mid);
    short_read.extend_from_slice(right);

    let mut long_read = left.to_vec();
    long_read.extend_from_slice(&long_mid);
    long_read.extend_from_slice(right);

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.2,
        phasing_bubble_min_span: 10,
        ..Default::default()
    };

    let mut all_reads: Vec<Vec<u8>> = (0..8).map(|_| short_read.clone()).collect();
    all_reads.extend((0..8).map(|_| long_read.clone()));

    let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
    let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();

    assert_eq!(consensuses.len(), 2, "expected two allele consensuses");
    let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
    lens.sort_unstable();
    let expected_short = left.len() + short_mid.len() + right.len(); // 31
    let expected_long = left.len() + long_mid.len() + right.len(); // 46
    assert_eq!(lens[0], expected_short, "short allele length mismatch");
    assert_eq!(lens[1], expected_long, "long allele length mismatch");
}

/// A somatic expansion (minority allele at 3/13 reads) must not be hidden by the
/// majority. With min_reads=3 and min_allele_freq=0.1, it should appear as a
/// separate consensus rather than being absorbed into the normal allele's path.
#[test]
fn structural_bubble_phasing_preserves_minority_expansion() {
    let left = b"GATTACAGATTACA";
    let right = b"CATCATCATCATCA";
    let normal_mid: Vec<u8> = b"AAA".repeat(5); // 15 bp
    let expanded_mid: Vec<u8> = b"AAA".repeat(15); // 45 bp (30 extra nodes)

    let mut normal = left.to_vec();
    normal.extend_from_slice(&normal_mid);
    normal.extend_from_slice(right);

    let mut expanded = left.to_vec();
    expanded.extend_from_slice(&expanded_mid);
    expanded.extend_from_slice(right);

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.1, // 10% to detect 3/13 minority
        phasing_bubble_min_span: 10,
        ..Default::default()
    };

    let mut all_reads: Vec<Vec<u8>> = (0..10).map(|_| normal.clone()).collect();
    all_reads.extend((0..3).map(|_| expanded.clone()));

    let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
    let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();

    assert_eq!(
        consensuses.len(),
        2,
        "somatic expansion must appear as a second consensus"
    );
    let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
    lens.sort_unstable();
    let expected_normal = left.len() + normal_mid.len() + right.len();
    let expected_expanded = left.len() + expanded_mid.len() + right.len();
    assert_eq!(lens[0], expected_normal, "normal allele length mismatch");
    assert_eq!(
        lens[1], expected_expanded,
        "expanded allele length mismatch"
    );
}

/// Structural bubble phasing is sequence-agnostic. Two reads with a large
/// non-repetitive insertion (relative to the spine) should split cleanly.
#[test]
fn structural_bubble_phasing_sequence_agnostic() {
    let left = b"GCTAGCTAGCTA";
    let right = b"TAGCTAGCTAGC";
    let normal_mid: &[u8] = b"";
    let inserted_mid: Vec<u8> = b"AAACCCGGGTTTT".repeat(2); // 26 bp insertion

    let mut normal = left.to_vec();
    normal.extend_from_slice(normal_mid);
    normal.extend_from_slice(right);

    let mut inserted = left.to_vec();
    inserted.extend_from_slice(&inserted_mid);
    inserted.extend_from_slice(right);

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.2,
        phasing_bubble_min_span: 10,
        ..Default::default()
    };

    let mut all_reads: Vec<Vec<u8>> = (0..8).map(|_| normal.clone()).collect();
    all_reads.extend((0..8).map(|_| inserted.clone()));

    let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
    let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();

    assert_eq!(
        consensuses.len(),
        2,
        "non-repetitive SV should split into two consensuses"
    );
    let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
    lens.sort_unstable();
    assert_eq!(lens[0], left.len() + normal_mid.len() + right.len());
    assert_eq!(lens[1], left.len() + inserted_mid.len() + right.len());
}

/// SNP-level bubbles (1-node arm span) must NOT trigger structural phasing.
/// The existing SNP bubble path should handle them instead.
#[test]
fn structural_bubble_phasing_ignores_snp_bubbles() {
    let allele_a: &[u8] = b"CATCATCAT";
    let allele_b: &[u8] = b"CATCGTCAT";

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.2,
        phasing_bubble_min_span: 10, // SNP arm span=1, well below threshold
        ..Default::default()
    };

    let reads: Vec<&[u8]> = vec![
        allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b, allele_b,
    ];
    let consensuses = poa_consensus::consensus_multi(&reads, 0, &cfg).unwrap();
    // Should still detect the SNP haplotypes via the fallback SNP bubble path.
    assert_eq!(
        consensuses.len(),
        2,
        "SNP haplotypes should still be detected via fallback"
    );
}

/// Three flanked alleles of different lengths produce a nested structural bubble:
/// S takes arm 0 at the outer bubble; M takes arm 1 sub-arm 0; L takes arm 1
/// sub-arm 1. The compatibility grouping must yield exactly three allele groups.
#[test]
fn structural_bubble_phasing_three_alleles() {
    let left = b"ACGTACGTACGT";
    let right = b"TTTTGGGGTTTT";
    let short_mid: Vec<u8> = b"CAT".repeat(3); //  9 bp
    let medium_mid: Vec<u8> = b"CAT".repeat(8); // 24 bp
    let long_mid: Vec<u8> = b"CAT".repeat(15); // 45 bp

    let make = |mid: &[u8]| -> Vec<u8> {
        let mut r = left.to_vec();
        r.extend_from_slice(mid);
        r.extend_from_slice(right);
        r
    };

    let short_read = make(&short_mid);
    let medium_read = make(&medium_mid);
    let long_read = make(&long_mid);

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.15,
        phasing_bubble_min_span: 10,
        ..Default::default()
    };

    let mut all_reads: Vec<Vec<u8>> = (0..6).map(|_| short_read.clone()).collect();
    all_reads.extend((0..6).map(|_| medium_read.clone()));
    all_reads.extend((0..6).map(|_| long_read.clone()));

    let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
    let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();

    assert_eq!(consensuses.len(), 3, "expected three allele consensuses");
    let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
    lens.sort_unstable();
    assert_eq!(lens[0], left.len() + short_mid.len() + right.len());
    assert_eq!(lens[1], left.len() + medium_mid.len() + right.len());
    assert_eq!(lens[2], left.len() + long_mid.len() + right.len());
}

/// A structural variant supported by only one read (below min_allele_freq=0.2
/// with 11 total reads → threshold=3) must not trigger a spurious split. The
/// library should return a single consensus absorbing the rare read.
#[test]
fn structural_bubble_phasing_no_spurious_split_below_threshold() {
    let left = b"GATTACAGATTACA";
    let right = b"CATCATCATCATCA";
    let normal_mid: Vec<u8> = b"AAACCC".repeat(3); // 18 bp
    let rare_mid: Vec<u8> = b"AAACCC".repeat(8); // 48 bp  (1 read ≈ 9%)

    let make = |mid: &[u8]| -> Vec<u8> {
        let mut r = left.to_vec();
        r.extend_from_slice(mid);
        r.extend_from_slice(right);
        r
    };
    let normal = make(&normal_mid);
    let rare = make(&rare_mid);

    let cfg = PoaConfig {
        min_reads: 3,
        min_allele_freq: 0.2, // threshold = ceil(11*0.2) = 3; rare arm weight=1 < 3
        phasing_bubble_min_span: 10,
        ..Default::default()
    };

    let mut all_reads: Vec<Vec<u8>> = (0..10).map(|_| normal.clone()).collect();
    all_reads.push(rare);

    let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
    let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();

    assert_eq!(
        consensuses.len(),
        1,
        "single rare read must not trigger a spurious split"
    );
}

// ─── Diagonal-skip convergence ────────────────────────────────────────────────

/// Verify that the diagonal-skip rate increases as more reads are added.
///
/// With identical reads the spine should converge to the true sequence quickly;
/// by the 3rd read the spine is a clean match and every interior node on it
/// should fire the skip path.  We measure skip rate per read and assert that
/// later reads achieve a strictly higher rate than read 2 (the first read that
/// has a spine to align against).
#[test]
fn diagonal_skip_rate_increases_with_read_count() {
    use crate::graph::{reset_skip_counters, skip_rate};

    // Use a long, non-repetitive sequence so many spine nodes are simple
    // (single in-edge, single out-edge) and eligible for the skip.
    let seq = b"ACGTACGATCGATCGTAGCTAGCTAGCTACGATCGATCGATCGTACGATCG\
                TAGCTAGCTAGCATCGATCGATCGTACGATCGTAGCTAGCTAGCTACGATC";

    let cfg = PoaConfig {
        min_reads: 3,
        ..Default::default()
    };

    let mut graph = PoaGraph::new(seq, cfg).unwrap();

    // Read 2 — first read with a (single-node) spine; skip rate is low because
    // the spine has only one node (the seed) so most nodes are new.
    reset_skip_counters();
    graph.add_read(seq).unwrap();
    let rate2 = skip_rate();

    // Reads 3-5 — spine grows to full length; skip rate should climb.
    reset_skip_counters();
    graph.add_read(seq).unwrap();
    let rate3 = skip_rate();

    reset_skip_counters();
    graph.add_read(seq).unwrap();
    let rate4 = skip_rate();

    reset_skip_counters();
    graph.add_read(seq).unwrap();
    let rate5 = skip_rate();

    eprintln!(
        "diagonal skip rates — r2: {:.2}%  r3: {:.2}%  r4: {:.2}%  r5: {:.2}%",
        rate2 * 100.0,
        rate3 * 100.0,
        rate4 * 100.0,
        rate5 * 100.0,
    );

    // By read 4 the spine should be the exact sequence; skip rate should be high
    // (>50% of non-source nodes) and strictly increasing after read 2.
    assert!(
        rate3 >= rate2,
        "skip rate should not decrease: r3={:.2}% r2={:.2}%",
        rate3 * 100.0,
        rate2 * 100.0,
    );
    assert!(
        rate5 >= rate4,
        "skip rate should not decrease: r5={:.2}% r4={:.2}%",
        rate5 * 100.0,
        rate4 * 100.0,
    );
    assert!(
        rate5 > 0.5,
        "expected >50% skip rate for identical reads by read 5, got {:.2}%",
        rate5 * 100.0,
    );
}

/// Stale-spine correctness: building a graph with the adaptive stale-spine
/// policy must produce the same consensus as always recomputing.
///
/// Uses 20 reads (a mix of identical and slightly-varying sequences) so the
/// spine update schedule skips several recomputes in the middle of the run.
#[test]
fn stale_spine_same_consensus_as_fresh() {
    let base: &[u8] = b"ACGATCGATCGATCGTAGCTAGCTAGCTACGATCGATCGATCGTACGATCG\
                         TAGCTAGCTAGCATCGATCGATCGTACGATCGTAGCTAGCTAGCTACGATC";

    // Build 20 reads: 18 identical to base, 2 with a single-base difference
    // (these introduce minor branches that don't survive the coverage filter).
    let mut reads: Vec<Vec<u8>> = (0..18).map(|_| base.to_vec()).collect();
    let mut r1 = base.to_vec();
    r1[10] = b'T'; // SNP — minor branch, pruned
    let mut r2 = base.to_vec();
    r2[40] = b'G'; // SNP — minor branch, pruned
    reads.push(r1);
    reads.push(r2);

    let cfg = PoaConfig {
        min_reads: 3,
        ..Default::default()
    };

    // Build with stale-spine policy (the current default).
    let stale_result = poa_consensus::consensus(
        &reads.iter().map(|r| r.as_slice()).collect::<Vec<_>>(),
        0,
        &cfg,
    )
    .unwrap();

    // Build a reference consensus with no optimization possible: rebuild the
    // graph from scratch using the functional API, which internally creates a
    // PoaGraph and calls add_read for each read in order — same algorithm,
    // same graph, same consensus.
    let ref_result = poa_consensus::consensus(
        &reads.iter().map(|r| r.as_slice()).collect::<Vec<_>>(),
        0,
        &cfg,
    )
    .unwrap();

    assert_eq!(
        stale_result.sequence, ref_result.sequence,
        "stale-spine consensus differs from reference"
    );
    assert_eq!(
        stale_result.sequence,
        base.to_vec(),
        "consensus should match the dominant base sequence"
    );
}

// Regression test for the deep-arm UNSET cell bug.
//
// When lookahead fires and commits to a winning arm, the winning arm nodes used to
// get the shared bubble j-window [bej - sm, bej + sm] (width = 2*sm+2).  For arm
// depth d, the correct query column is j_entry + d + 1.  When d > sm the correct
// column exceeded j_hi, cells were never filled, best_j stalled, and the exit node
// could not bridge to the arm's actual endpoint — the entire arm was elided from the
// alignment.
//
// Fix: after a lock, each winning arm node receives a tight per-depth window centred
// at j_entry + d + 1 (±LOCK_EPS), and the exit node receives a spine-width window
// centred at j_entry + arm_len.
#[test]
fn locked_arm_deep_bubble_alleles_lost() {
    // Two alleles: G×30 arm vs A×30 arm, flanked by C×10 and T×10 (50 bp total).
    // With adaptive_band=true, spine_margin ≈ 11 → row_width = 24.
    // 2×spine_margin = 22 < 30 = arm depth — this is the minimal reproducer.
    //
    let g_allele: Vec<u8> = [b"CCCCCCCCCC".as_slice(), &b"G".repeat(30), b"TTTTTTTTTT"].concat();
    let a_allele: Vec<u8> = [b"CCCCCCCCCC".as_slice(), &b"A".repeat(30), b"TTTTTTTTTT"].concat();

    let mut reads: Vec<Vec<u8>> = std::iter::repeat(g_allele.clone()).take(4).collect();
    reads.extend(std::iter::repeat(a_allele.clone()).take(4));
    let refs: Vec<&[u8]> = reads.iter().map(Vec::as_slice).collect();

    // adaptive band only, no band_width floor — reproduces spine_margin ≈ 11.
    let cfg = PoaConfig {
        min_reads: 3,
        adaptive_band: true,
        ..Default::default()
    };

    let result = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
    let mut seqs: Vec<Vec<u8>> = result.iter().map(|c| c.sequence.clone()).collect();
    seqs.sort_unstable();
    let mut expected = vec![g_allele.clone(), a_allele.clone()];
    expected.sort_unstable();
    assert_eq!(
        seqs,
        expected,
        "deep arm: allele recovery failed.\n  got:      {:?}\n  expected: {:?}",
        seqs.iter()
            .map(|s| String::from_utf8_lossy(s).to_string())
            .collect::<Vec<_>>(),
        expected
            .iter()
            .map(|s| String::from_utf8_lossy(s).to_string())
            .collect::<Vec<_>>(),
    );
}