twitcher 0.6.9

use std::{collections::HashSet, sync::Arc};

use anyhow::{Context as _, bail};
use rust_htslib::bam::{self, FetchDefinition, Read as _};
use tracing::{error, warn};

use crate::{
    common::{bam_filters::BamFilters, coords::GenomeRegion},
    counter,
    vcf::pipeline::{clusterizer::phasing::GtAndPhase, record::InputRecord},
};

/// Confidence thresholds for local read-based phasing. These gate when an unphased
/// heterozygous record inside a cluster is allowed to be assigned to a haplotype from
/// read evidence. Only relevant when a phasing BAM (`--phasing-bam`) is supplied; a
/// record that clears neither gate stays unphased and continues to block its cluster.
#[derive(clap::Args, Debug, Clone)]
pub struct PhasingSettings {
    /// Minimum number of reads that must support the winning haplotype assignment.
    /// Assignments backed by fewer reads are rejected, guarding against decisions made on
    /// too little coverage; the record then stays unphased.
    #[arg(long = "phasing-min-reads", default_value_t = 3)]
    pub min_reads: u32,

    /// Minimum fraction of informative reads that must agree on the assignment [0.5–1.0].
    /// "Informative" means the read covers the position and carries a recognisable allele.
    /// Higher values reject conflicting or noisy evidence; lower values phase more
    /// aggressively. Records below this threshold remain unphased.
    #[arg(long = "phasing-confidence", default_value_t = 0.8)]
    pub confidence: f64,

    /// Minimum read mapping quality (MAPQ) to use a read as phasing evidence. Reads below this
    /// are skipped entirely, so mis-mapped reads cannot contribute spurious allele calls.
    #[arg(long = "phasing-min-mapq", default_value_t = 20)]
    pub min_mapq: u8,

    /// Minimum base quality (Phred) at the variant position for a read's allele call to count.
    /// A read whose base(s) at the locus fall below this are treated as uninformative there,
    /// so likely sequencing errors do not become votes. Deletions (no read base) are unaffected.
    #[arg(long = "phasing-min-base-qual", default_value_t = 20)]
    pub min_base_qual: u8,

    #[command(flatten)]
    pub bam_filters: BamFilters,
}

impl Default for PhasingSettings {
    fn default() -> Self {
        Self {
            min_reads: 3,
            confidence: 0.8,
            min_mapq: 20,
            min_base_qual: 20,
            bam_filters: BamFilters::default(),
        }
    }
}

/// Phasing orientation chosen for an unphased record relative to its GT allele slots `[g0, g1]`.
/// `Direct` keeps the slot order (`g0` on H0, `g1` on H1 → written `g0|g1`); `Flipped` swaps it
/// (`g1` on H0, `g0` on H1 → written `g1|g0`).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Orientation {
    Direct,
    Flipped,
}

/// Outcome of a per-record phasing decision, carrying the reason a record stayed unresolved
/// so it can be attributed to a granular counter (see [`record_outcome`]).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Decision {
    /// Assigned to a haplotype with the chosen orientation.
    Resolved(Orientation),
    /// No informative reads cover the record (missing evidence).
    NoEvidence,
    /// The winning orientation is supported by fewer than `min_reads` reads.
    InsufficientReads,
    /// Informative reads agree below the `confidence` threshold (conflict / noise).
    LowConfidence,
}

/// Increment the matching `phasing.records.*` counter for `decision` and map it to the
/// orientation to apply (`None` when the record stays unresolved).
fn record_outcome(decision: Decision, decisions: &mut Vec<Decision>) -> Option<Orientation> {
    decisions.push(decision);
    match decision {
        Decision::Resolved(o) => {
            counter!("phasing.records.resolved").inc(1);
            Some(o)
        }
        Decision::NoEvidence => {
            counter!("phasing.records.unresolved.no_evidence").inc(1);
            None
        }
        Decision::InsufficientReads => {
            counter!("phasing.records.unresolved.insufficient_reads").inc(1);
            None
        }
        Decision::LowConfidence => {
            counter!("phasing.records.unresolved.low_confidence").inc(1);
            None
        }
    }
}

/// Per-node metadata, precomputed once and indexed in lockstep with `node_indices`.
struct NodeInfo {
    pos: i64,
    /// GT allele indices of the two het slots (both present, distinct).
    slots: [u32; 2],
    /// Full allele list `[REF, ALT1, ...]` owned so the read loop need not re-parse.
    alleles: Vec<Vec<u8>>,
}

impl NodeInfo {
    /// Map the allele a read carries to `(slot, is_ref)`: slot `false` = slot 0, `true` = slot 1;
    /// `is_ref` = the carried allele is REF (index 0). `None` = the read carries neither slot
    /// allele (e.g. REF at a `1/2` site) or does not cover.
    fn slot_label(&self, read: &bam::Record, min_base_qual: u8) -> Option<(bool, bool)> {
        let refs: Vec<&[u8]> = self.alleles.iter().map(Vec::as_slice).collect();
        let idx = read_allele_index_at_pos(read, self.pos, &refs, min_base_qual)?;
        if idx == self.slots[0] {
            Some((false, idx == 0))
        } else if idx == self.slots[1] {
            Some((true, idx == 0))
        } else {
            None
        }
    }
}

#[derive(Debug)]
#[allow(dead_code)]
pub enum BamPhaseResolverError {
    MultiplePhaseSets(HashSet<i32>),
    BamReadError(rust_htslib::errors::Error),
    Other(anyhow::Error),
}

/// Attempt to resolve `UnphasedHet` records in `records` by read co-occurrence.
///
/// Reads covering the cluster `region` are fetched from `reader`. For each read the
/// supported allele (ALT or REF) at each record position is determined, and a
/// same/different evidence matrix is built. Records are assigned to haplotypes when
/// `settings.min_reads` reads support the dominant assignment AND at least
/// `settings.confidence` of informative reads agree.
///
/// Records that cannot be resolved remain `UnphasedHet` and will still block the cluster.
/// Records that are already `Phased` are used as anchors and are not changed.
pub fn resolve_phasing(
    records: &mut [Arc<InputRecord>],
    region: &GenomeRegion,
    reader: &mut bam::IndexedReader,
    settings: &PhasingSettings,
) -> Result<Vec<Decision>, BamPhaseResolverError> {
    let unphased_indices: Vec<usize> = records
        .iter()
        .enumerate()
        .filter_map(|(i, rec)| {
            let gtp = rec.gt();
            // Records with a missing allele (e.g. `1/.`) carry no resolvable second slot.
            (gtp.is_unphased_het() && !gtp.has_missing()).then_some(i)
        })
        .collect();

    if unphased_indices.is_empty() {
        return Ok(vec![]);
    }

    // Collect anchors: phased het records with their index in `records` and GtAndPhase.
    // Hom records are excluded — they carry no directional (H0 vs H1) evidence.
    let anchors: Vec<(usize, GtAndPhase)> = records
        .iter()
        .enumerate()
        .filter_map(|(i, rec)| {
            let gtp = rec.gt();
            if gtp.phase.is_some() && !gtp.is_hom() {
                Some((i, gtp))
            } else {
                None
            }
        })
        .collect();

    // If anchors span multiple distinct PSs, skip phasing — too ambiguous.
    let distinct_ps: HashSet<i32> = anchors.iter().filter_map(|(_, ps)| ps.phase).collect();
    if distinct_ps.len() > 1 {
        return Err(BamPhaseResolverError::MultiplePhaseSets(distinct_ps));
    }

    // Build the ordered list of nodes participating in evidence (anchors first, then unresolved).
    // Each node_indices[k] is an index into `records`.
    let node_indices: Vec<usize> = anchors
        .iter()
        .map(|(i, _)| *i)
        .chain(unphased_indices.iter().copied())
        .collect();
    let anchor_count = anchors.len();

    // Precompute per-node metadata once (avoids re-parsing genotypes / re-allocating `alleles()`
    // for every read). All nodes are het with both slots present, so the lookups below cannot
    // fail; they are checked anyway to keep a malformed cluster from taking the run down.
    let node_info: Vec<NodeInfo> = node_indices
        .iter()
        .map(|&ri| {
            let rec = records
                .get(ri)
                .with_context(|| format!("Node index {ri} is out of the cluster's records"))?;
            let gtp = rec.gt();
            let [Some(slot_0), Some(slot_1)] = gtp.alleles else {
                bail!("Node at record index {ri} is not a het genotype with both slots present");
            };
            Ok(NodeInfo {
                pos: rec.pos(),
                slots: [slot_0, slot_1],
                alleles: rec.alleles().iter().map(|a| a.to_vec()).collect(),
            })
        })
        .collect::<anyhow::Result<_>>()
        .map_err(BamPhaseResolverError::Other)?;

    let evidence = accumulate_evidence(reader, region, &node_info, settings)?;

    // Determine PS for resolved records.
    let (assigned_phaseset, has_anchors) =
        if let Some(ps) = anchors.first().and_then(|(_, gtp)| gtp.phase) {
            (ps, true)
        } else {
            // Without an anchor the cluster gets its own phase set, keyed on the first record.
            let ps = records
                .first()
                .map_or(0, |r| i32::try_from(r.pos() + 1).unwrap_or(i32::MAX));
            (ps, false)
        };

    let (assignment, decisions) = assign_orientations(
        &evidence,
        anchor_count,
        unphased_indices.len(),
        has_anchors,
        settings,
    );

    apply_assignments(
        records,
        &unphased_indices,
        &assignment,
        &node_info,
        anchor_count,
        assigned_phaseset,
    );

    Ok(decisions)
}

/// A symmetric `n × n` tally over the nodes, stored flat in row-major order: a lookup is one
/// bounds-checked offset instead of a bounds-checked row lookup plus a pointer chase into a
/// separate per-row allocation.
struct CoOccurrence {
    counts: Vec<i32>,
    n: usize,
}

impl CoOccurrence {
    fn new(n: usize) -> Self {
        Self {
            counts: vec![0; n.saturating_mul(n)],
            n,
        }
    }

    /// The number of nodes this matrix covers.
    const fn len(&self) -> usize {
        self.n
    }

    /// The tally for the node pair `(i, j)`. An out-of-range pair reads as `0`, i.e. as "these
    /// two nodes were never seen together", which can only ever withhold a phasing decision.
    fn get(&self, i: usize, j: usize) -> i32 {
        self.counts.get(i * self.n + j).copied().unwrap_or(0)
    }

    /// Count one read for the node pair `(i, j)`, in both directions.
    fn tally(&mut self, i: usize, j: usize) {
        for ix in [i * self.n + j, j * self.n + i] {
            if let Some(count) = self.counts.get_mut(ix) {
                *count = count.saturating_add(1);
            }
        }
    }
}

/// Symmetric co-occurrence matrices over the nodes, indexed in lockstep with `node_indices`.
struct Evidence {
    /// `same.get(i, j)`: read carries the same slot at both i and j, with at least one being a
    /// non-REF allele (ALT/ALT, or ALT1/ALT1 at a 1/2 node). REF/REF agreement is excluded — it
    /// only says "neither mutation is present" and carries no haplotype linkage, so counting it as
    /// Direct evidence would let false-REF calls (e.g. unmatched indels in repeats) collapse
    /// records onto the anchor haplotype.
    same: CoOccurrence,
    /// `diff.get(i, j)`: read carries different slots at i and j (always involves a non-REF allele).
    diff: CoOccurrence,
}

/// Scan the reads covering `region` and tally per-read slot agreement between every pair of nodes.
fn accumulate_evidence(
    reader: &mut bam::IndexedReader,
    region: &GenomeRegion,
    node_info: &[NodeInfo],
    settings: &PhasingSettings,
) -> Result<Evidence, BamPhaseResolverError> {
    let n = node_info.len();
    let mut same = CoOccurrence::new(n);
    let mut diff = CoOccurrence::new(n);

    let fd = match FetchDefinition::try_from(region) {
        Ok(fd) => fd,
        Err(e) => {
            warn!("local_phasing: cannot build fetch definition: {e}");
            return Err(BamPhaseResolverError::Other(e));
        }
    };
    if let Err(e) = reader.fetch(fd) {
        warn!("local_phasing: fetch failed: {e}");
        return Err(BamPhaseResolverError::BamReadError(e));
    }

    let mut bam_rec = bam::Record::new();
    while let Some(result) = reader.read(&mut bam_rec) {
        if result.is_err() {
            warn!("local_phasing: BAM read error; stopping scan for this cluster");
            break;
        }

        if !settings.bam_filters.apply(&bam_rec) || bam_rec.mapq() < settings.min_mapq {
            continue;
        }

        // Per-node label for this read: (slot, is_ref); None = uninformative.
        let support: Vec<Option<(bool, bool)>> = node_info
            .iter()
            .map(|info| info.slot_label(&bam_rec, settings.min_base_qual))
            .collect();

        for (i, label_i) in support.iter().enumerate() {
            let Some((si, ri)) = *label_i else { continue };
            for (j, label_j) in support.iter().enumerate().skip(i + 1) {
                let Some((sj, rj)) = *label_j else { continue };
                if si == sj {
                    // REF/REF agreement carries no haplotype linkage → not directional evidence.
                    if ri && rj {
                        continue;
                    }
                    same.tally(i, j);
                } else {
                    diff.tally(i, j);
                }
            }
        }
    }

    Ok(Evidence { same, diff })
}

/// Turn the co-occurrence evidence into a per-unresolved-record orientation, plus the decision
/// (including the reason for staying unresolved) that produced it.
///
/// Nodes `0..anchor_count` are the phased anchors, `anchor_count..` the unresolved records in
/// the order of `unphased_indices`.
fn assign_orientations(
    Evidence { same, diff }: &Evidence,
    anchor_count: usize,
    unresolved_count: usize,
    has_anchors: bool,
    settings: &PhasingSettings,
) -> (Vec<Option<Orientation>>, Vec<Decision>) {
    let n = same.len();
    let mut assignment: Vec<Option<Orientation>> = Vec::with_capacity(unresolved_count);

    // Pick the orientation whose vote total clears the min-reads and confidence gates.
    // `v_direct` favours keeping the node's GT slot order, `v_flipped` favours swapping it.
    let decide = |v_direct: i32, v_flipped: i32| -> Decision {
        let total = (v_direct + v_flipped).max(0).unsigned_abs();
        if total == 0 {
            return Decision::NoEvidence;
        }
        let (winner, winner_count) = if v_direct >= v_flipped {
            (Orientation::Direct, v_direct.max(0).unsigned_abs())
        } else {
            (Orientation::Flipped, v_flipped.max(0).unsigned_abs())
        };
        if winner_count < settings.min_reads {
            return Decision::InsufficientReads;
        }
        let fraction = f64::from(winner_count) / f64::from(total);
        if fraction < settings.confidence {
            return Decision::LowConfidence;
        }
        Decision::Resolved(winner)
    };

    let mut decisions = Vec::new();

    for k in 0..unresolved_count {
        let node_k = anchor_count + k;

        if has_anchors {
            // A phased anchor's slot 0 is H0 and slot 1 is H1 by definition, and node labels use
            // the same slot convention. So a read agreeing on labels (same) across anchor and node
            // supports the Direct orientation; a disagreement (diff) supports Flipped — regardless
            // of the anchor's genotype.
            let mut v_direct = 0i32;
            let mut v_flipped = 0i32;
            for a_node in 0..anchor_count {
                v_direct += same.get(a_node, node_k);
                v_flipped += diff.get(a_node, node_k);
            }
            assignment.push(record_outcome(decide(v_direct, v_flipped), &mut decisions));
        } else {
            // No anchors: greedy partition.
            if k == 0 {
                // Seed record: assign Direct only when co-read evidence meets min_reads.
                // Use the maximum pairwise co-coverage with any other node as the read count.
                let max_pair_cov: u32 = (0..n)
                    .filter(|&j| j != node_k)
                    .map(|j| {
                        (same.get(node_k, j) + diff.get(node_k, j))
                            .max(0)
                            .unsigned_abs()
                    })
                    .max()
                    .unwrap_or(0);
                let decision = if max_pair_cov == 0 {
                    Decision::NoEvidence
                } else if max_pair_cov < settings.min_reads {
                    Decision::InsufficientReads
                } else {
                    Decision::Resolved(Orientation::Direct)
                };
                assignment.push(record_outcome(decision, &mut decisions));
                continue;
            }
            let mut v_direct = 0i32;
            let mut v_flipped = 0i32;
            for (j, item) in assignment.iter().enumerate() {
                let node_j = anchor_count + j;
                match item {
                    Some(Orientation::Direct) => {
                        v_direct += same.get(node_j, node_k);
                        v_flipped += diff.get(node_j, node_k);
                    }
                    Some(Orientation::Flipped) => {
                        v_flipped += same.get(node_j, node_k);
                        v_direct += diff.get(node_j, node_k);
                    }
                    None => {}
                }
            }
            assignment.push(record_outcome(decide(v_direct, v_flipped), &mut decisions));
        }
    }

    (assignment, decisions)
}

/// Attach the resolved phased genotype and the cluster's phase set to the records that
/// local phasing could orient.
fn apply_assignments(
    records: &mut [Arc<InputRecord>],
    unphased_indices: &[usize],
    assignment: &[Option<Orientation>],
    node_info: &[NodeInfo],
    anchor_count: usize,
    assigned_phaseset: i32,
) {
    // Node k occupies position anchor_count + k in node_info, in lockstep with the k-th
    // unresolved record and the k-th assignment.
    let nodes = node_info.iter().skip(anchor_count);
    for ((&rec_idx, orientation), info) in unphased_indices.iter().zip(assignment).zip(nodes) {
        let Some(orientation) = orientation else {
            continue;
        };
        let [g0, g1] = info.slots;

        // Direct: g0 on H0, g1 on H1 → "g0|g1". Flipped: swap → "g1|g0".
        let (new_a0, new_a1): (u32, u32) = match orientation {
            Orientation::Direct => (g0, g1),
            Orientation::Flipped => (g1, g0),
        };

        let Some(record) = records.get_mut(rec_idx) else {
            error!("Cannot apply local phasing: record index {rec_idx} is out of range");
            continue;
        };

        // The genotype is attached to the record rather than written into it: the record is
        // bound to the input header, which may not declare GT or PS at all. The writer emits
        // both onto the output records it builds.
        *record = Arc::new(record.with_local_gt(GtAndPhase {
            alleles: [Some(new_a0), Some(new_a1)],
            phase: Some(assigned_phaseset),
        }));
    }
}

/// Determine whether `read` carries the ALT or REF allele of a VCF record at `vcf_pos`.
///
/// Returns `Some(true)` = ALT, `Some(false)` = REF, `None` = read doesn't cover or ambiguous.
/// Supports SNVs, MNVs (equal-length multi-base substitutions), simple deletions, and
/// simple insertions. Complex alleles with unequal non-trivial lengths return `None`.
fn read_allele_at_pos(
    read: &bam::Record,
    vcf_pos: i64,
    ref_allele: &[u8],
    alt_allele: &[u8],
    min_base_qual: u8,
) -> Option<bool> {
    use rust_htslib::bam::record::Cigar;
    // Anchor-trim before the position guard: reads starting at the trimmed position must not
    // be dropped by comparing against the pre-trim (anchor) coordinate.
    let (vcf_pos, ref_allele, alt_allele) =
        match (ref_allele.split_first(), alt_allele.split_first()) {
            (Some((r_anchor, r_rest)), Some((a_anchor, a_rest))) if r_anchor == a_anchor => {
                (vcf_pos + 1, r_rest, a_rest)
            }
            _ => (vcf_pos, ref_allele, alt_allele),
        };

    if vcf_pos < read.pos() {
        return None;
    }

    let cigar = read.cigar();
    let seq = read.seq();
    let qual = read.qual();
    let mut ref_pos: i64 = read.pos();
    let mut query_pos: usize = 0;

    for op in &cigar {
        match op {
            Cigar::Match(raw_len) | Cigar::Equal(raw_len) | Cigar::Diff(raw_len) => {
                let op_len = i64::from(*raw_len);
                if ref_pos + op_len > vcf_pos {
                    return allele_within_match(
                        &MatchOp {
                            seq: &seq,
                            qual,
                            query_pos,
                            offset: usize::try_from(vcf_pos - ref_pos).ok()?,
                            op_len: *raw_len as usize,
                        },
                        ref_allele,
                        alt_allele,
                        min_base_qual,
                    );
                }
                ref_pos += op_len;
                query_pos += *raw_len as usize;
            }
            Cigar::Ins(ins_len) => {
                // Insertion at this ref_pos: REF is empty, ALT is non-empty after trim.
                if ref_allele.is_empty() && !alt_allele.is_empty() && ref_pos == vcf_pos {
                    // Only this ALT if the inserted bases match in length and sequence.
                    // Otherwise the read carries a *different* insertion → neither this ALT
                    // nor REF, so report a miss (lets a tri-allelic caller try the other ALT).
                    let il = *ins_len as usize;
                    if il == alt_allele.len()
                        && query_pos + il <= seq.len()
                        && alt_allele
                            .iter()
                            .enumerate()
                            .all(|(i, a)| seq[query_pos + i].eq_ignore_ascii_case(a))
                    {
                        if (0..il)
                            .any(|i| qual.get(query_pos + i).is_some_and(|&q| q < min_base_qual))
                        {
                            return None; // low-quality inserted base → unreliable ALT call
                        }
                        return Some(true);
                    }
                    return None;
                }
                query_pos += *ins_len as usize;
            }
            Cigar::Del(del_len) => {
                let dl = i64::from(*del_len);
                if !ref_allele.is_empty() && alt_allele.is_empty() {
                    // Looking for a deletion. Match only an exact-length deletion starting at
                    // our position. A spanning deletion of a different length is a *different*
                    // event → miss (not this ALT, not REF).
                    if ref_pos == vcf_pos
                        && usize::try_from(dl).is_ok_and(|dl| dl == ref_allele.len())
                    {
                        return Some(true);
                    } else if ref_pos <= vcf_pos && vcf_pos < ref_pos + dl {
                        return None;
                    }
                } else if ref_pos <= vcf_pos && vcf_pos < ref_pos + dl {
                    // Read has a deletion at this position but we're not looking for one.
                    return None;
                }
                ref_pos += dl;
            }
            Cigar::SoftClip(sc_len) => {
                query_pos += *sc_len as usize;
            }
            Cigar::RefSkip(rs_len) => {
                ref_pos += i64::from(*rs_len);
            }
            Cigar::HardClip(_) | Cigar::Pad(_) => {}
        }

        if ref_pos > vcf_pos {
            break;
        }
    }

    None
}

/// The aligned block a VCF position falls into, as located by [`read_allele_at_pos`].
struct MatchOp<'a> {
    seq: &'a bam::record::Seq<'a>,
    qual: &'a [u8],
    /// Query offset of the start of the CIGAR op.
    query_pos: usize,
    /// Offset of the VCF position within the CIGAR op.
    offset: usize,
    /// Length of the CIGAR op.
    op_len: usize,
}

/// Determine the allele of a read whose CIGAR aligns (M/=/X) across the VCF position.
///
/// Returns `Some(true)` = ALT, `Some(false)` = REF, `None` = ambiguous or unreliable.
/// The alleles are anchor-trimmed, so an empty side means an indel.
fn allele_within_match(
    op: &MatchOp<'_>,
    ref_allele: &[u8],
    alt_allele: &[u8],
    min_base_qual: u8,
) -> Option<bool> {
    let &MatchOp {
        seq,
        qual,
        query_pos,
        offset,
        op_len,
    } = op;
    match (ref_allele.len(), alt_allele.len()) {
        (_, 0) | (0, _) => {
            // After anchor-trim, one side is empty: deletion or insertion VCF.
            // Being inside a Match op means the event is absent → REF allele.
            // A low-quality base here is an unreliable REF call → uninformative.
            // `.get` guards against a record with no quality string (no panic).
            if qual
                .get(query_pos + offset)
                .is_some_and(|&q| q < min_base_qual)
            {
                return None;
            }
            Some(false)
        }
        (rlen, alen) if rlen == alen => {
            // SNV (rlen=1) or MNV (rlen>1).
            let avail = op_len.saturating_sub(offset);
            if rlen > avail {
                return None; // block extends beyond this CIGAR op
            }
            let qstart = query_pos + offset;
            if (0..rlen).any(|i| qual.get(qstart + i).is_some_and(|&q| q < min_base_qual)) {
                return None; // low-quality base → unreliable call
            }
            let mut matches_alt = true;
            let mut matches_ref = true;
            for (i, (r, a)) in ref_allele.iter().zip(alt_allele).enumerate() {
                let base_uc = seq[qstart + i].to_ascii_uppercase();
                if base_uc != a.to_ascii_uppercase() {
                    matches_alt = false;
                }
                if base_uc != r.to_ascii_uppercase() {
                    matches_ref = false;
                }
            }
            match (matches_alt, matches_ref) {
                (true, _) => Some(true),
                (false, true) => Some(false),
                _ => None,
            }
        }
        _ => None, // complex allele (unequal non-trivial lengths)
    }
}

/// Determine which allele index `read` carries at `vcf_pos`.
///
/// Returns `Some(0)` for REF, `Some(i)` for `alleles[i]` (an ALT), or `None` when the read does
/// not cover the position or the allele is ambiguous. `alleles` is `[REF, ALT1, ALT2, ...]`.
///
/// Each ALT is tested independently against REF via [`read_allele_at_pos`] (which applies its own
/// anchor-trim), so multi-allelic sites — including two distinct indel ALTs — are discriminated.
fn read_allele_index_at_pos(
    read: &bam::Record,
    vcf_pos: i64,
    alleles: &[&[u8]],
    min_base_qual: u8,
) -> Option<u32> {
    if alleles.len() < 2 {
        return None;
    }
    let reference = *alleles.first()?;
    let mut ref_seen = false;
    for (i, alt) in alleles.iter().enumerate().skip(1) {
        match read_allele_at_pos(read, vcf_pos, reference, alt, min_base_qual) {
            #[allow(clippy::cast_possible_truncation)]
            Some(true) => return Some(i as u32), // carries this ALT
            Some(false) => ref_seen = true, // genuine REF match
            None => {}
        }
    }
    ref_seen.then_some(0)
}

#[cfg(test)]
mod tests {
    use std::sync::Arc;

    use rust_htslib::bcf::record::GenotypeAllele;
    use rust_htslib::{
        bam::{self, record::CigarString},
        bcf::{self, Header, Writer},
    };

    use super::super::phasing::Haplotype;
    use super::*;
    use crate::common::{contig::ContigName, coords::GenomePosition};

    fn hap_of(gtp: &GtAndPhase) -> Option<Haplotype> {
        match gtp.alleles {
            [Some(a), _] if a > 0 => Some(Haplotype::H0),
            [Some(0), Some(b)] if b > 0 => Some(Haplotype::H1),
            _ => None,
        }
    }

    // ── helpers ───────────────────────────────────────────────────────────────

    fn make_vcf_header() -> Header {
        let mut h = Header::new();
        h.push_record(b"##contig=<ID=chr1,length=1000000>");
        h.push_record(b"##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">");
        h.push_record(b"##FORMAT=<ID=PS,Number=1,Type=Integer,Description=\"Phase set\">");
        h.push_sample(b"S1");
        h
    }

    fn make_unphased_het(writer: &Writer, pos: i64, ref_a: &[u8], alt_a: &[u8]) -> InputRecord {
        let mut rec = writer.empty_record();
        let rid = rec.header().name2rid(b"chr1").unwrap();
        rec.set_rid(Some(rid));
        rec.set_pos(pos);
        rec.set_alleles(&[ref_a, alt_a]).unwrap();
        rec.push_genotypes(&[GenotypeAllele::Unphased(0), GenotypeAllele::Unphased(1)])
            .unwrap();
        InputRecord::new(rec)
    }

    fn make_phased_het(
        writer: &Writer,
        pos: i64,
        ref_a: &[u8],
        alt_a: &[u8],
        hap: Haplotype,
        ps: i32,
    ) -> InputRecord {
        let (a0, a1) = if hap == Haplotype::H1 {
            (GenotypeAllele::Unphased(0), GenotypeAllele::Phased(1))
        } else {
            (GenotypeAllele::Unphased(1), GenotypeAllele::Phased(0))
        };
        let mut rec = writer.empty_record();
        let rid = rec.header().name2rid(b"chr1").unwrap();
        rec.set_rid(Some(rid));
        rec.set_pos(pos);
        rec.set_alleles(&[ref_a, alt_a]).unwrap();
        rec.push_genotypes(&[a0, a1]).unwrap();
        rec.push_format_integer(b"PS", &[ps]).unwrap();
        InputRecord::new(rec)
    }

    /// Build a phased het record for a tri-allelic (or larger) site.
    /// `alleles` = [REF, ALT1, ALT2, ...]; `gt` = (`pos0_idx`, `pos1_idx`) as GT allele indices.
    fn make_phased_het_multi(
        writer: &Writer,
        pos: i64,
        alleles: &[&[u8]],
        gt: (u32, u32),
        ps: i32,
    ) -> InputRecord {
        let mut rec = writer.empty_record();
        let rid = rec.header().name2rid(b"chr1").unwrap();
        rec.set_rid(Some(rid));
        rec.set_pos(pos);
        rec.set_alleles(alleles).unwrap();
        rec.push_genotypes(&[
            GenotypeAllele::Unphased(gt.0 as i32),
            GenotypeAllele::Phased(gt.1 as i32),
        ])
        .unwrap();
        rec.push_format_integer(b"PS", &[ps]).unwrap();
        InputRecord::new(rec)
    }

    /// Build a minimal BAM record for CIGAR-based unit tests (no index needed). MAPQ 60 and
    /// uniform Phred-60 base qualities, so it passes the default MAPQ/base-quality gates.
    fn make_bam_record(pos: i64, seq: &[u8], cigar: &CigarString) -> bam::Record {
        make_bam_record_q(pos, seq, cigar, &vec![60u8; seq.len()])
    }

    /// Like [`make_bam_record`] but with explicit per-base qualities (Phred). MAPQ is 60.
    fn make_bam_record_q(pos: i64, seq: &[u8], cigar: &CigarString, quals: &[u8]) -> bam::Record {
        let mut rec = bam::Record::new();
        rec.set_pos(pos);
        rec.set_tid(0);
        rec.unset_unmapped();
        rec.set(b"read1", Some(cigar), seq, quals);
        rec.set_mapq(60);
        rec.cache_cigar();
        rec
    }

    fn make_bam_header() -> bam::Header {
        let mut header = bam::Header::new();
        let mut sq = bam::header::HeaderRecord::new(b"SQ");
        sq.push_tag(b"SN", "chr1");
        sq.push_tag(b"LN", 1_000_000i32);
        header.push_record(&sq);
        header
    }

    fn write_bam_and_index(mut records: Vec<bam::Record>, path: &std::path::Path) {
        let header = make_bam_header();
        let mut writer = bam::Writer::from_path(path, &header, bam::Format::Bam).unwrap();
        for rec in &mut records {
            rec.set_tid(0);
            writer.write(rec).unwrap();
        }
        drop(writer);
        bam::index::build(path, None, bam::index::Type::Bai, 1).unwrap();
    }

    fn make_indexed_reader(path: &std::path::Path) -> bam::IndexedReader {
        bam::IndexedReader::from_path(path).unwrap()
    }

    fn make_region(pos: usize, len: usize) -> GenomeRegion {
        GenomeRegion::new_bounded(GenomePosition::new_0(ContigName::new(b"chr1"), pos), len)
    }

    fn lenient_settings() -> PhasingSettings {
        PhasingSettings {
            min_reads: 1,
            confidence: 0.6,
            ..PhasingSettings::default()
        }
    }

    // ── read_allele_at_pos unit tests ─────────────────────────────────────────

    #[test]
    fn test_allele_support_snv_alt() {
        use rust_htslib::bam::record::Cigar;
        let rec = make_bam_record(100, b"TTTTT", &CigarString(vec![Cigar::Match(5)]));
        assert_eq!(read_allele_at_pos(&rec, 102, b"A", b"T", 0), Some(true));
    }

    #[test]
    fn test_allele_support_snv_ref() {
        use rust_htslib::bam::record::Cigar;
        let rec = make_bam_record(100, b"AAAAA", &CigarString(vec![Cigar::Match(5)]));
        assert_eq!(read_allele_at_pos(&rec, 102, b"A", b"T", 0), Some(false));
    }

    #[test]
    fn test_allele_support_snv_no_cover() {
        use rust_htslib::bam::record::Cigar;
        let rec = make_bam_record(100, b"AAAAA", &CigarString(vec![Cigar::Match(5)]));
        assert_eq!(read_allele_at_pos(&rec, 200, b"A", b"T", 0), None);
    }

    #[test]
    fn test_allele_support_mnv_alt() {
        use rust_htslib::bam::record::Cigar;
        // MNV: REF=AT (positions 102-103), ALT=GC
        // seq: AA GC AA (GC at offset 2-3)
        let rec = make_bam_record(100, b"AAGCAA", &CigarString(vec![Cigar::Match(6)]));
        assert_eq!(read_allele_at_pos(&rec, 102, b"AT", b"GC", 0), Some(true));
    }

    #[test]
    fn test_allele_support_mnv_ref() {
        use rust_htslib::bam::record::Cigar;
        let rec = make_bam_record(100, b"AAATAA", &CigarString(vec![Cigar::Match(6)]));
        assert_eq!(read_allele_at_pos(&rec, 102, b"AT", b"GC", 0), Some(false));
    }

    #[test]
    fn test_allele_support_deletion_alt() {
        use rust_htslib::bam::record::Cigar;
        // VCF: REF=AT, ALT=A → deletion of T at pos 101 (after anchor-trim: pos 101, REF=T, ALT="")
        // Read: 1M 1D 1M at pos=100
        let rec = make_bam_record(
            100,
            b"AA",
            &CigarString(vec![Cigar::Match(1), Cigar::Del(1), Cigar::Match(1)]),
        );
        assert_eq!(read_allele_at_pos(&rec, 100, b"AT", b"A", 0), Some(true));
    }

    #[test]
    fn test_allele_support_insertion_alt() {
        use rust_htslib::bam::record::Cigar;
        // VCF: REF=A, ALT=AT → insertion of T after pos=100
        // After anchor-trim: pos=101, REF="", ALT=T
        // Read: 1M 1I 1M at pos=100
        let rec = make_bam_record(
            100,
            b"ATB",
            &CigarString(vec![Cigar::Match(1), Cigar::Ins(1), Cigar::Match(1)]),
        );
        assert_eq!(read_allele_at_pos(&rec, 100, b"A", b"AT", 0), Some(true));
    }

    #[test]
    fn test_allele_support_deletion_ref() {
        use rust_htslib::bam::record::Cigar;
        // VCF: REF=AT, ALT=A → deletion at pos 101 (after trim).
        // Read has 2M — no deletion present → REF allele.
        let rec = make_bam_record(100, b"AA", &CigarString(vec![Cigar::Match(2)]));
        assert_eq!(read_allele_at_pos(&rec, 100, b"AT", b"A", 0), Some(false));
    }

    #[test]
    fn test_allele_support_insertion_ref() {
        use rust_htslib::bam::record::Cigar;
        // VCF: REF=A, ALT=AT → insertion after pos 100 (after trim: pos 101, ref="", alt="T").
        // Read has 2M — no insertion present → REF allele.
        let rec = make_bam_record(100, b"AA", &CigarString(vec![Cigar::Match(2)]));
        assert_eq!(read_allele_at_pos(&rec, 100, b"A", b"AT", 0), Some(false));
    }

    #[test]
    fn test_allele_support_deletion_read_starts_at_trimmed_pos() {
        use rust_htslib::bam::record::Cigar;
        // VCF: REF=AT, ALT=A at pos=100. After anchor-trim: vcf_pos=101, REF=T, ALT="".
        // Read starts at pos=101 with a Del(1) — deletion present → ALT.
        // Bug: pre-trim guard `100 < 101` drops this read before trimming.
        let rec = make_bam_record(
            101,
            b"A",
            &CigarString(vec![Cigar::Del(1), Cigar::Match(1)]),
        );
        assert_eq!(read_allele_at_pos(&rec, 100, b"AT", b"A", 0), Some(true));
    }

    #[test]
    fn test_allele_support_ref_read_starts_at_trimmed_pos() {
        use rust_htslib::bam::record::Cigar;
        // Same scenario, but read has no deletion → REF allele.
        let rec = make_bam_record(101, b"TA", &CigarString(vec![Cigar::Match(2)]));
        assert_eq!(read_allele_at_pos(&rec, 100, b"AT", b"A", 0), Some(false));
    }

    // ── resolve_phasing unit tests ────────────────────────────────────────────

    #[test]
    fn test_resolve_no_unphased_het() {
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(
            &mut input,
            &region,
            &mut reader,
            &PhasingSettings::default(),
        )
        .unwrap();
        assert_eq!(input.len(), 1);
        // Anchor record: already phased, should be unchanged.
        let gtp = input[0].gt();
        assert!(gtp.phase.is_some(), "anchor should remain phased");
        assert_eq!(hap_of(&gtp), Some(Haplotype::H1));
    }

    #[test]
    fn test_resolve_two_unphased_hets_no_coverage() {
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        let before = counter!("phasing.records.unresolved.no_evidence").get();
        resolve_phasing(
            &mut input,
            &region,
            &mut reader,
            &PhasingSettings::default(),
        )
        .unwrap();
        assert!(input[0].gt().is_unphased_het());
        assert!(input[1].gt().is_unphased_het());
        // Both records lack covering reads → both counted as no_evidence.
        assert_eq!(
            counter!("phasing.records.unresolved.no_evidence").get() - before,
            2
        );
    }

    #[test]
    fn test_resolve_two_unphased_hets_same_hap() {
        use rust_htslib::bam::record::Cigar;
        // r1 at pos=100 REF=A ALT=T, r2 at pos=110 REF=G ALT=C.
        // Two reads covering both positions: T at 100 AND C at 110 → same hap.
        let mut seq = vec![b'T'];
        seq.extend(b"AAAAAAAAA"); // positions 101-109
        seq.push(b'C'); // position 110
        seq.extend(b"AAAAAAAAA"); // positions 111-119
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let read1 = make_bam_record(100, &seq, &cigar);
        let read2 = make_bam_record(100, &seq, &cigar);

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![read1, read2], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        let before = counter!("phasing.records.resolved").get();
        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp0 = input[0].gt();
        let gtp1 = input[1].gt();
        assert!(gtp0.phase.is_some(), "r1 should be resolved");
        assert!(gtp1.phase.is_some(), "r2 should be resolved");
        assert_eq!(
            hap_of(&gtp0),
            hap_of(&gtp1),
            "both records should be on the same haplotype"
        );
        // Seed + dependent record both resolved.
        assert_eq!(counter!("phasing.records.resolved").get() - before, 2);
    }

    #[test]
    fn test_resolve_two_unphased_hets_opposite_haps() {
        use rust_htslib::bam::record::Cigar;
        // r1 at 100 REF=A ALT=T, r2 at 110 REF=G ALT=C.
        // hap0-reads: T at 100, G (ref) at 110.
        // hap1-reads: A (ref) at 100, C at 110.
        let mut s0 = vec![b'T'];
        s0.extend(b"AAAAAAAAA");
        s0.push(b'G'); // ref at 110
        s0.extend(b"AAAAAAAAA");

        let mut s1 = vec![b'A']; // ref at 100
        s1.extend(b"AAAAAAAAA");
        s1.push(b'C'); // alt at 110
        s1.extend(b"AAAAAAAAA");

        let cigar = CigarString(vec![Cigar::Match(20)]);
        let reads = vec![
            make_bam_record(100, &s0, &cigar),
            make_bam_record(100, &s0, &cigar),
            make_bam_record(100, &s1, &cigar),
            make_bam_record(100, &s1, &cigar),
        ];

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp0 = input[0].gt();
        let gtp1 = input[1].gt();
        assert!(
            gtp0.phase.is_some() && gtp1.phase.is_some(),
            "both should be resolved"
        );
        assert_ne!(
            hap_of(&gtp0),
            hap_of(&gtp1),
            "records should be on opposite haplotypes"
        );
    }

    #[test]
    fn test_resolve_below_min_reads_threshold() {
        use rust_htslib::bam::record::Cigar;
        // Only 1 read; min_reads=3 → the second record should remain UnphasedHet.
        let mut seq = vec![b'T'];
        seq.extend(b"AAAAAAAAA");
        seq.push(b'C');
        seq.extend(b"AAAAAAAAA");
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let read1 = make_bam_record(100, &seq, &cigar);

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![read1], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        let settings = PhasingSettings {
            min_reads: 3,
            confidence: 0.0,
            ..PhasingSettings::default()
        };
        let before = counter!("phasing.records.unresolved.insufficient_reads").get();
        resolve_phasing(&mut input, &region, &mut reader, &settings).unwrap();
        // Seed has 1 covering read < min_reads=3 → insufficient_reads.
        assert_eq!(
            counter!("phasing.records.unresolved.insufficient_reads").get() - before,
            1
        );
        assert!(
            input[1].gt().is_unphased_het(),
            "second record should remain unphased with only 1 read < min_reads=3"
        );
        assert!(
            input[0].gt().is_unphased_het(),
            "seed (first) record should also remain unphased: 1 read < min_reads=3"
        );
    }

    #[test]
    fn test_resolve_below_confidence_threshold() {
        use rust_htslib::bam::record::Cigar;
        // 2 same-hap reads + 2 diff-hap reads → 50% confidence < 0.8 → second stays unphased.
        let mut s_same = vec![b'T'];
        s_same.extend(b"AAAAAAAAA");
        s_same.push(b'C');
        s_same.extend(b"AAAAAAAAA");

        let mut s_diff = vec![b'T'];
        s_diff.extend(b"AAAAAAAAA");
        s_diff.push(b'G'); // ref at 110
        s_diff.extend(b"AAAAAAAAA");

        let cigar = CigarString(vec![Cigar::Match(20)]);
        let reads = vec![
            make_bam_record(100, &s_same, &cigar),
            make_bam_record(100, &s_same, &cigar),
            make_bam_record(100, &s_diff, &cigar),
            make_bam_record(100, &s_diff, &cigar),
        ];

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        let settings = PhasingSettings {
            min_reads: 1,
            confidence: 0.8,
            ..PhasingSettings::default()
        };
        let before = counter!("phasing.records.unresolved.low_confidence").get();
        resolve_phasing(&mut input, &region, &mut reader, &settings).unwrap();
        // 50% agreement < 0.8 confidence → low_confidence.
        assert_eq!(
            counter!("phasing.records.unresolved.low_confidence").get() - before,
            1
        );
        assert!(
            input[1].gt().is_unphased_het(),
            "second record should remain unphased below confidence threshold (50% < 0.8)"
        );
    }

    #[test]
    fn test_resolve_unphased_with_phased_anchor() {
        use rust_htslib::bam::record::Cigar;
        // r1 (pos=100): UnphasedHet, REF=A ALT=T.
        // r2 (pos=110): Phased H1, phaseset=100 (anchor).
        // Reads carry both ALTs → UnphasedHet should resolve to H1.
        let mut seq = vec![b'T']; // ALT at pos 100
        seq.extend(b"AAAAAAAAA");
        seq.push(b'C'); // ALT at pos 110
        seq.extend(b"AAAAAAAAA");
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let read1 = make_bam_record(100, &seq, &cigar);
        let read2 = make_bam_record(100, &seq, &cigar);

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het(&w, 100, b"A", b"T");
        let r2 = make_phased_het(&w, 110, b"G", b"C", Haplotype::H1, 100);
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![read1, read2], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp0 = input[0].gt();
        assert!(gtp0.phase.is_some(), "UnphasedHet should be resolved");
        assert_eq!(
            hap_of(&gtp0),
            Some(Haplotype::H1),
            "should resolve to same hap as anchor"
        );
        // Anchor (r2) unchanged:
        let gtp1 = input[1].gt();
        assert!(gtp1.phase.is_some(), "anchor should remain phased");
        assert_eq!(hap_of(&gtp1), Some(Haplotype::H1));
    }

    #[test]
    fn test_resolve_refref_only_does_not_phase() {
        use rust_htslib::bam::record::Cigar;
        // Anchor 0|1@100 (H1), unphased 0/1@110. Every read is REF at BOTH positions
        // (no read carries either ALT). Pure REF/REF co-occurrence carries no haplotype
        // linkage and must NOT resolve the unphased record. Regression for a real case
        // (chr1:212925716-212925734) where false-REF calls from unmatched indels in an STR
        // produced REF/REF-only evidence and wrongly phased the record onto the anchor.
        let mut seq = vec![b'A']; // pos100 REF
        seq.extend(b"AAAAAAAAA"); // 101..109
        seq.push(b'G'); // pos110 REF
        seq.extend(b"AAAAAAAAA"); // 111..119
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let reads: Vec<bam::Record> = (0..5).map(|_| make_bam_record(100, &seq, &cigar)).collect();

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let unphased = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unphased)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        let before = counter!("phasing.records.unresolved.no_evidence").get();
        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        assert!(
            input[1].gt().is_unphased_het(),
            "REF/REF-only evidence must not phase the record"
        );
        assert_eq!(
            counter!("phasing.records.unresolved.no_evidence").get() - before,
            1
        );
    }

    #[test]
    fn test_resolve_refref_reads_do_not_inflate() {
        use rust_htslib::bam::record::Cigar;
        // 3 ALT/ALT reads give genuine Direct linkage; 5 REF/REF reads carry no linkage and
        // must be ignored (not counted toward confidence). Node still resolves to H1.
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let mut alt = vec![b'T'];
        alt.extend(b"AAAAAAAAA");
        alt.push(b'C');
        alt.extend(b"AAAAAAAAA");
        let mut refr = vec![b'A'];
        refr.extend(b"AAAAAAAAA");
        refr.push(b'G');
        refr.extend(b"AAAAAAAAA");
        let mut reads: Vec<bam::Record> =
            (0..3).map(|_| make_bam_record(100, &alt, &cigar)).collect();
        reads.extend((0..5).map(|_| make_bam_record(100, &refr, &cigar)));

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let unphased = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unphased)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        // confidence=0.8 would be met either way; the point is the 5 REF/REF reads do not
        // participate, so resolution rests purely on the 3 ALT/ALT reads (3/3 = 1.0).
        let settings = PhasingSettings {
            min_reads: 3,
            confidence: 0.8,
            ..PhasingSettings::default()
        };
        resolve_phasing(&mut input, &region, &mut reader, &settings).unwrap();
        let gtp = input[1].gt();
        assert!(
            gtp.phase.is_some(),
            "genuine ALT/ALT linkage should still resolve"
        );
        assert_eq!(hap_of(&gtp), Some(Haplotype::H1));
    }

    #[test]
    fn test_resolve_low_mapq_reads_ignored() {
        use rust_htslib::bam::record::Cigar;
        // Anchor 0|1@100 (H1), unphased 0/1@110. Reads carry genuine ALT/ALT linkage that would
        // resolve the record, but all have MAPQ 5 (< default 20). Mis-mappable reads must not
        // count, so the record stays unphased.
        let mut seq = vec![b'T']; // ALT@100
        seq.extend(b"AAAAAAAAA");
        seq.push(b'C'); // ALT@110
        seq.extend(b"AAAAAAAAA");
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let reads: Vec<bam::Record> = (0..5)
            .map(|_| {
                let mut r = make_bam_record(100, &seq, &cigar);
                r.set_mapq(5);
                r
            })
            .collect();

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let unphased = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unphased)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        // Default settings → min_mapq=20. lenient_settings keeps min_mapq=20 via ..default too,
        // so use it for the lenient min_reads/confidence but the MAPQ gate still applies.
        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        assert!(
            input[1].gt().is_unphased_het(),
            "reads below the MAPQ floor must not phase the record"
        );
    }

    #[test]
    fn test_resolve_low_base_quality_at_variant_ignored() {
        use rust_htslib::bam::record::Cigar;
        // Anchor 0|1@100 (H1), unphased 0/1@110. Reads carry ALT/ALT linkage, but the base AT the
        // unphased variant (pos 110, query index 10) has Phred 5 (< default 20). That call is
        // uninformative, so no evidence reaches the node and it stays unphased. The anchor base
        // (index 0) stays high-quality.
        let mut seq = vec![b'T']; // ALT@100
        seq.extend(b"AAAAAAAAA");
        seq.push(b'C'); // ALT@110
        seq.extend(b"AAAAAAAAA");
        let cigar = CigarString(vec![Cigar::Match(20)]);
        let mut quals = vec![60u8; seq.len()];
        quals[10] = 5; // low quality at pos 110
        let reads: Vec<bam::Record> = (0..5)
            .map(|_| make_bam_record_q(100, &seq, &cigar, &quals))
            .collect();

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let unphased = make_unphased_het(&w, 110, b"G", b"C");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unphased)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        // min_base_qual default is 20; lenient_settings keeps it via ..default.
        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        assert!(
            input[1].gt().is_unphased_het(),
            "a low-quality base at the variant must not phase the record"
        );
    }

    #[test]
    fn test_resolve_multiple_phaseset_early_exit() {
        // Two phased records with different PSs → ambiguous → return unchanged immediately.
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let r2 = make_phased_het(&w, 110, b"G", b"C", Haplotype::H0, 200); // different PS
        let r3 = make_unphased_het(&w, 120, b"C", b"G");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2), Arc::new(r3)];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 70);

        let err = resolve_phasing(
            &mut input,
            &region,
            &mut reader,
            &PhasingSettings::default(),
        )
        .unwrap_err();
        assert!(
            matches!(err, BamPhaseResolverError::MultiplePhaseSets(_)),
            "anchors spanning two distinct phasesets should be rejected"
        );
        assert!(
            input[2].gt().is_unphased_het(),
            "multiple PSs should cause early exit leaving UnphasedHet unchanged"
        );
    }

    #[test]
    fn test_resolve_2_1_anchor_correct_direction() {
        use rust_htslib::bam::record::Cigar;
        // Anchor: 2|1 at pos=100, alleles=[A, T, C], GT=(2,1) → H0=C(ALT2), H1=T(ALT1).
        // read_allele_at_pos uses alleles[0]=A and alleles[1]=T, so ALT1=T reads are visible
        // (Some(true)=H1). With the old `> 0` formula anchor_is_h0=true and ALT1 reads would
        // be credited to H0 — inverted. With `== Some(1)` anchor_is_h0=false → H1. ✓
        //
        // Unresolved: 0/1 at pos=110, REF=G, ALT=A.
        // Reads: 3 reads carrying T at pos=100 (ALT1=H1 of anchor) and A at pos=110.
        // Expected: unresolved resolves to H1.

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het_multi(&w, 100, &[b"A", b"T", b"C"], (2, 1), 100);
        let unresolved = make_unphased_het(&w, 110, b"G", b"A");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unresolved)];

        // seq[0]=T → ref-pos 100 (ALT1 of 2|1 anchor), seq[10]=A → ref-pos 110 (ALT of unresolved).
        let mut seq = vec![b'T'];
        seq.extend(b"AAAAAAAAA"); // pos 101..109
        seq.push(b'A'); // pos 110 = ALT of unresolved (ALT=A, REF=G, so A≠G confirms alt)
        let cigar = CigarString(vec![Cigar::Match(seq.len() as u32)]);
        let reads: Vec<bam::Record> = (0..3).map(|_| make_bam_record(100, &seq, &cigar)).collect();

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[1].gt();
        assert!(gtp.phase.is_some(), "unresolved should be phased");
        assert_eq!(
            hap_of(&gtp),
            Some(Haplotype::H1),
            "should be H1: ALT1 of 2|1 anchor is on H1"
        );
    }

    #[test]
    fn test_resolve_1_2_anchor_h0_evidence_correct() {
        use rust_htslib::bam::record::Cigar;
        // Anchor: 1|2 at pos=100, alleles=[A, T, C], GT=(1,2) → H0=T(ALT1), H1=C(ALT2).
        // ALT1=T reads are visible (Some(true)=H0). ALT2=C reads return None (invisible).
        // Unresolved: 0/1 at pos=110, REF=G, ALT=A.
        // Reads: 3 reads carrying T at pos=100 (H0) and A at pos=110.
        // Expected: unresolved resolves to H0.

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het_multi(&w, 100, &[b"A", b"T", b"C"], (1, 2), 100);
        let unresolved = make_unphased_het(&w, 110, b"G", b"A");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unresolved)];

        // seq[0]=T → ref-pos 100 (ALT1 of 1|2 anchor = H0), seq[10]=A → ref-pos 110.
        let mut seq = vec![b'T'];
        seq.extend(b"AAAAAAAAA"); // pos 101..109
        seq.push(b'A'); // pos 110 = ALT of unresolved (REF=G, ALT=A)
        let cigar = CigarString(vec![Cigar::Match(seq.len() as u32)]);
        let reads: Vec<bam::Record> = (0..3).map(|_| make_bam_record(100, &seq, &cigar)).collect();

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[1].gt();
        assert!(gtp.phase.is_some(), "unresolved should be phased");
        assert_eq!(
            hap_of(&gtp),
            Some(Haplotype::H0),
            "should be H0: ALT1 of 1|2 anchor is on H0"
        );
    }

    #[test]
    fn test_resolve_missing_allele_not_modified() {
        // 1/. record: is_unphased_het()=true but has_missing()=true.
        // With the has_missing guard it must be excluded from unphased_indices and
        // left completely unchanged (no PS field, no modified GT).
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();

        let mut rec_missing = w.empty_record();
        let rid = rec_missing.header().name2rid(b"chr1").unwrap();
        rec_missing.set_rid(Some(rid));
        rec_missing.set_pos(100);
        rec_missing.set_alleles(&[b"A", b"T"]).unwrap();
        rec_missing
            .push_genotypes(&[GenotypeAllele::Unphased(1), GenotypeAllele::UnphasedMissing])
            .unwrap();

        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(InputRecord::new(rec_missing))];

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(vec![], &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 30);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[0].gt();
        // has_missing → excluded from resolution → no PS written → still not phased
        assert!(
            !gtp.is_phased(),
            "1/. record must not be phased by resolve_phasing"
        );
        assert!(
            gtp.has_missing(),
            "alleles must remain unchanged (still has missing)"
        );
    }

    #[test]
    fn test_resolve_1_2_unphased_ref_reads_excluded() {
        use rust_htslib::bam::record::Cigar;
        // Anchor: 0|1 at pos=100, REF=A, ALT=T → H0=REF, H1=ALT.
        // Unresolved: 1/2 at pos=110, alleles=[G, C, X].
        //   Both GT allele indices > 0 → REF reads at this node must be excluded (Some(false)→None).
        //
        // Reads: 3 reads carry T at pos=100 (anchor H1=ALT) AND C at pos=110 (ALT1 of 1/2).
        //   → These should accumulate evidence toward H1 for the 1/2 node.
        // Noise reads: 3 reads carry T at pos=100 AND G (REF) at pos=110.
        //   → With the fix, REF at the 1/2 node returns None → not counted.
        //   → Without the fix, they would inflate diff and potentially prevent resolution.
        //
        // Expected: 1/2 node resolves to H1 (ALT1=C on H1 → written as 2|1).

        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);

        let mut rec_12 = w.empty_record();
        let rid = rec_12.header().name2rid(b"chr1").unwrap();
        rec_12.set_rid(Some(rid));
        rec_12.set_pos(110);
        rec_12.set_alleles(&[b"G", b"C", b"X"]).unwrap();
        rec_12
            .push_genotypes(&[GenotypeAllele::Unphased(1), GenotypeAllele::Unphased(2)])
            .unwrap();
        let mut input: Vec<Arc<InputRecord>> =
            vec![Arc::new(anchor), Arc::new(InputRecord::new(rec_12))];

        // Signal reads: seq[0]=T (ALT=H1 of anchor), seq[10]=C (ALT1 of 1/2).
        let mut sig_seq = vec![b'T'];
        sig_seq.extend(b"AAAAAAAAA"); // pos 101..109
        sig_seq.push(b'C'); // pos 110 = ALT1 of 1/2
        let cigar = CigarString(vec![Cigar::Match(sig_seq.len() as u32)]);
        let mut reads: Vec<bam::Record> = (0..3)
            .map(|_| make_bam_record(100, &sig_seq, &cigar))
            .collect();

        // Noise reads: T at 100, G=REF at 110. REF at 1/2 node must be filtered (Some(false)→None).
        let mut noise_seq = vec![b'T'];
        noise_seq.extend(b"AAAAAAAAA"); // pos 101..109
        noise_seq.push(b'G'); // pos 110 = REF of 1/2
        let noise_cigar = CigarString(vec![Cigar::Match(noise_seq.len() as u32)]);
        reads.extend((0..3).map(|_| make_bam_record(100, &noise_seq, &noise_cigar)));

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[1].gt();
        assert!(gtp.phase.is_some(), "1/2 unresolved should be phased");
        // hap_of checks alleles[0]>0 which is biallelic-only. For 2|1 that would return H0.
        // Check directly: ALT1 (index 1 in the allele array) must be at phased position alleles[1].
        assert_eq!(
            gtp.alleles[1],
            Some(1),
            "ALT1(C) of 1/2 should be at phased GT position (alleles[1]=1), meaning it is on H1"
        );
    }

    /// Build an unphased multi-allelic het record (e.g. `1/2`). `alleles` = [REF, ALT1, ...].
    fn make_unphased_het_multi(
        writer: &Writer,
        pos: i64,
        alleles: &[&[u8]],
        gt: (u32, u32),
    ) -> InputRecord {
        let mut rec = writer.empty_record();
        let rid = rec.header().name2rid(b"chr1").unwrap();
        rec.set_rid(Some(rid));
        rec.set_pos(pos);
        rec.set_alleles(alleles).unwrap();
        rec.push_genotypes(&[
            GenotypeAllele::Unphased(gt.0 as i32),
            GenotypeAllele::Unphased(gt.1 as i32),
        ])
        .unwrap();
        InputRecord::new(rec)
    }

    // ── read_allele_index_at_pos unit tests ───────────────────────────────────

    #[test]
    fn test_index_snv_picks_each_allele() {
        use rust_htslib::bam::record::Cigar;
        let alleles: &[&[u8]] = &[b"A", b"T", b"C"];
        let cigar = CigarString(vec![Cigar::Match(5)]);
        // offset 2 → pos 102
        let alt1 = make_bam_record(100, b"AATAA", &cigar);
        let alt2 = make_bam_record(100, b"AACAA", &cigar);
        let refr = make_bam_record(100, b"AAAAA", &cigar);
        let other = make_bam_record(100, b"AAGAA", &cigar);
        assert_eq!(read_allele_index_at_pos(&alt1, 102, alleles, 0), Some(1));
        assert_eq!(read_allele_index_at_pos(&alt2, 102, alleles, 0), Some(2));
        assert_eq!(read_allele_index_at_pos(&refr, 102, alleles, 0), Some(0));
        assert_eq!(read_allele_index_at_pos(&other, 102, alleles, 0), None);
    }

    #[test]
    fn test_index_alt2_not_mistaken_for_alt1() {
        use rust_htslib::bam::record::Cigar;
        // A read carrying ALT2 must not be reported as ALT1.
        let alleles: &[&[u8]] = &[b"A", b"T", b"C"];
        let rec = make_bam_record(100, b"AACAA", &CigarString(vec![Cigar::Match(5)]));
        assert_eq!(read_allele_index_at_pos(&rec, 102, alleles, 0), Some(2));
    }

    #[test]
    fn test_index_insertion_discrimination() {
        use rust_htslib::bam::record::Cigar;
        // REF=A, ALT1=AT (ins T), ALT2=AG (ins G) at pos=100.
        let alleles: &[&[u8]] = &[b"A", b"AT", b"AG"];
        let cigar = CigarString(vec![Cigar::Match(1), Cigar::Ins(1), Cigar::Match(1)]);
        let ins_t = make_bam_record(100, b"ATB", &cigar);
        let ins_g = make_bam_record(100, b"AGB", &cigar);
        assert_eq!(read_allele_index_at_pos(&ins_t, 100, alleles, 0), Some(1));
        assert_eq!(read_allele_index_at_pos(&ins_g, 100, alleles, 0), Some(2));
    }

    #[test]
    fn test_index_deletion_length_discrimination() {
        use rust_htslib::bam::record::Cigar;
        // REF=ATG, ALT1=A → deletion of "TG" (len 2 @101) after single-base anchor-trim.
        let alleles: &[&[u8]] = &[b"ATG", b"A"];
        // Exact-length deletion (Del 2) → ALT1.
        let del2 = make_bam_record(
            100,
            b"AC",
            &CigarString(vec![Cigar::Match(1), Cigar::Del(2), Cigar::Match(1)]),
        );
        // Wrong-length deletion (Del 1) must NOT be credited to the len-2 ALT.
        let del1 = make_bam_record(
            100,
            b"AG",
            &CigarString(vec![Cigar::Match(1), Cigar::Del(1), Cigar::Match(1)]),
        );
        assert_eq!(read_allele_index_at_pos(&del2, 100, alleles, 0), Some(1));
        assert_eq!(
            read_allele_index_at_pos(&del1, 100, alleles, 0),
            None,
            "a different-length deletion is not this ALT"
        );
    }

    // ── native multi-allelic resolve_phasing tests ────────────────────────────

    #[test]
    fn test_resolve_1_2_anchor_alt2_evidence() {
        use rust_htslib::bam::record::Cigar;
        // Anchor: 1|2 at pos=100, alleles=[A,T,C] → H0=T(ALT1), H1=C(ALT2).
        // Only ALT2 (C) reads cover the anchor — invisible to the old biallelic algorithm.
        // Unresolved: 0/1 at pos=110, REF=G ALT=A. Reads carry C@100 (anchor H1) + A@110.
        // Expected: unresolved resolves to H1.
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het_multi(&w, 100, &[b"A", b"T", b"C"], (1, 2), 100);
        let unresolved = make_unphased_het(&w, 110, b"G", b"A");
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unresolved)];

        let mut seq = vec![b'C']; // ALT2 of anchor = H1
        seq.extend(b"AAAAAAAAA");
        seq.push(b'A'); // ALT of unresolved
        let cigar = CigarString(vec![Cigar::Match(seq.len() as u32)]);
        let reads: Vec<bam::Record> = (0..3)
            .map(|_| make_bam_record(100, &seq, &cigar.clone()))
            .collect();

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[1].gt();
        assert!(gtp.phase.is_some(), "unresolved should be phased");
        assert_eq!(
            hap_of(&gtp),
            Some(Haplotype::H1),
            "ALT2 of 1|2 anchor is on H1 → unresolved ALT should be H1"
        );
    }

    #[test]
    fn test_resolve_native_1_2_unresolved_by_alt2() {
        use rust_htslib::bam::record::Cigar;
        // Anchor: 0|1 at pos=100, REF=A ALT=T → H1=ALT1=T.
        // Unresolved: 1/2 at pos=110, alleles=[G,C,A] → ALT1=C, ALT2=A.
        // Reads carry T@100 (anchor H1) and ALT2=A@110.
        // Expected: ALT2(A) lands on H1 → GT written 1|2 → alleles[1]=Some(2).
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let anchor = make_phased_het(&w, 100, b"A", b"T", Haplotype::H1, 100);
        let unresolved = make_unphased_het_multi(&w, 110, &[b"G", b"C", b"A"], (1, 2));
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(anchor), Arc::new(unresolved)];

        let mut seq = vec![b'T']; // ALT1 of anchor = H1
        seq.extend(b"AAAAAAAAA");
        seq.push(b'A'); // ALT2 of 1/2
        let cigar = CigarString(vec![Cigar::Match(seq.len() as u32)]);
        let reads: Vec<bam::Record> = (0..3)
            .map(|_| make_bam_record(100, &seq, &cigar.clone()))
            .collect();

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        let gtp = input[1].gt();
        assert!(gtp.phase.is_some(), "1/2 unresolved should be phased");
        assert_eq!(
            gtp.alleles[1],
            Some(2),
            "ALT2(A) co-occurs with anchor H1 → ALT2 at phased position (alleles[1]=2)"
        );
    }

    #[test]
    fn test_resolve_no_anchor_pair_of_1_2() {
        use rust_htslib::bam::record::Cigar;
        // Two 1/2 records, no anchors. r1=[A,T,C]@100, r2=[G,C,A]@110.
        // Reads carry ALT1 of each (T@100, C@110) → co-occurrence links them.
        let h = make_vcf_header();
        let w = Writer::from_stdout(&h, true, bcf::Format::Vcf).unwrap();
        let r1 = make_unphased_het_multi(&w, 100, &[b"A", b"T", b"C"], (1, 2));
        let r2 = make_unphased_het_multi(&w, 110, &[b"G", b"C", b"A"], (1, 2));
        let mut input: Vec<Arc<InputRecord>> = vec![Arc::new(r1), Arc::new(r2)];

        let mut seq = vec![b'T']; // ALT1 of r1
        seq.extend(b"AAAAAAAAA");
        seq.push(b'C'); // ALT1 of r2
        let cigar = CigarString(vec![Cigar::Match(seq.len() as u32)]);
        let reads: Vec<bam::Record> = (0..3)
            .map(|_| make_bam_record(100, &seq, &cigar.clone()))
            .collect();

        let dir = tempfile::tempdir().unwrap();
        let bam_path = dir.path().join("test.bam");
        write_bam_and_index(reads, &bam_path);
        let mut reader = make_indexed_reader(&bam_path);
        let region = make_region(90, 50);

        resolve_phasing(&mut input, &region, &mut reader, &lenient_settings()).unwrap();
        assert!(input[0].gt().phase.is_some(), "r1 should be resolved");
        assert!(input[1].gt().phase.is_some(), "r2 should be resolved");
    }
}