//! Molecule assembly orchestrator.
//!
//! Ties together the parser, UMI clustering, endpoint fingerprinting, and
//! consensus calling into a complete molecule assembly pipeline.  The entry
//! point is [`assemble_molecules`].
//!
//! # Pipeline
//!
//! 1. Group [`ParsedReadPair`]s by canonical UMI.
//! 2. Within each UMI group, apply Hamming-distance clustering
//!    ([`cluster_umi_pairs`]) to merge near-identical UMIs.
//! 3. Within each cluster, sub-group by endpoint fingerprint compatibility to
//!    detect UMI collisions (two distinct molecules sharing the same UMI by chance).
//! 4. For each fingerprint-compatible sub-group:
//!    - Separate reads into forward and reverse strand families.
//!    - Skip the group if the total family size is below [`AssemblerConfig::min_family_size`].
//!    - Call single-strand consensus for each strand.
//!    - Call duplex consensus when both strands are present.
//!    - Build a [`Molecule`] from the results.

use ahash::AHashMap;

use kam_core::molecule::{CanonicalUmiPair, ConsensusRead, FamilyType, Molecule, Strand};

use crate::clustering::cluster_umi_pairs;
use crate::consensus::{single_strand_consensus, ConsensusConfig, DisagreementStrategy};
use crate::fingerprint::{compute_endpoint_fingerprint, fingerprints_duplex_compatible};
use crate::parser::{ParseStats, ParsedReadPair};

// ── Configuration ─────────────────────────────────────────────────────────────

/// Full configuration for the molecule assembly pipeline.
///
/// # Example
/// ```
/// use kam_assemble::assembler::AssemblerConfig;
///
/// let config = AssemblerConfig::default();
/// assert_eq!(config.max_hamming_distance, 1);
/// assert_eq!(config.min_family_size, 1);
/// assert_eq!(config.min_duplex_reads, 1);
/// ```
#[derive(Debug, Clone)]
pub struct AssemblerConfig {
    /// Parser configuration (chemistry preset + optional filters).
    pub consensus: ConsensusConfig,
    /// How to handle strand disagreements during duplex consensus.
    pub disagreement_strategy: DisagreementStrategy,
    /// Maximum Hamming distance between two canonical UMI pairs for them to be
    /// merged into the same cluster.  Default: 1.
    pub max_hamming_distance: u32,
    /// Minimum total reads (fwd + rev) required to retain a family.  Families
    /// below this threshold are dropped and counted in
    /// [`AssemblyStats::n_families_below_min_size`].  Default: 1.
    pub min_family_size: u8,
    /// Minimum reads per strand required before calling duplex consensus.
    /// Default: 1.
    pub min_duplex_reads: u8,
}

impl Default for AssemblerConfig {
    /// Default assembly configuration.
    ///
    /// Uses standard consensus settings, Hamming distance 1 clustering, and
    /// accepts all families of size ≥ 1.
    fn default() -> Self {
        Self {
            consensus: ConsensusConfig::default(),
            disagreement_strategy: DisagreementStrategy::default(),
            max_hamming_distance: 1,
            min_family_size: 1,
            min_duplex_reads: 1,
        }
    }
}

// ── Stats ─────────────────────────────────────────────────────────────────────

/// Counters accumulated during a single call to [`assemble_molecules`].
///
/// # Example
/// ```
/// use kam_assemble::assembler::AssemblyStats;
///
/// let stats = AssemblyStats::default();
/// assert_eq!(stats.n_molecules, 0);
/// ```
#[derive(Debug, Clone, Default)]
pub struct AssemblyStats {
    /// Parse-stage counters (populated by the caller; passthrough in
    /// [`assemble_molecules`] which operates on already-parsed reads).
    pub parse_stats: ParseStats,
    /// Total molecules assembled (after all filters).
    pub n_molecules: u64,
    /// Molecules with reads on both strands (duplex).
    pub n_duplex: u64,
    /// Molecules with reads on the forward strand only (simplex).
    pub n_simplex_fwd: u64,
    /// Molecules with reads on the reverse strand only (simplex).
    pub n_simplex_rev: u64,
    /// Molecules consisting of exactly one read.
    pub n_singletons: u64,
    /// Families dropped because they fell below [`AssemblerConfig::min_family_size`].
    pub n_families_below_min_size: u64,
    /// Sub-groups produced by fingerprint splitting (each extra split after the
    /// first counts as one detected collision).
    pub n_umi_collisions_detected: u64,
}

// ── Main function ─────────────────────────────────────────────────────────────

/// Assemble molecules from an already-parsed set of read pairs.
///
/// Reads are grouped by canonical UMI, Hamming-distance clustered, then
/// sub-grouped by endpoint fingerprint compatibility.  Consensus is called for
/// each resulting family and a [`Molecule`] is constructed.
///
/// The `parse_stats` field of the returned [`AssemblyStats`] is **not** filled
/// in here — populate it from your parsing loop before calling this function,
/// then merge as needed.
///
/// # Example
/// ```
/// use kam_assemble::assembler::{assemble_molecules, AssemblerConfig};
/// use kam_assemble::parser::{ParserConfig, ParseResult, parse_read_pair};
///
/// let config_parser = ParserConfig::default();
/// let r1_seq  = b"ACGTATGNNNNNNNNNN";
/// let r1_qual = b"IIIIIIIIIIIIIIIII";
/// let r2_seq  = b"TGCATAGNNNNNNNNNN";
/// let r2_qual = b"IIIIIIIIIIIIIIIII";
/// let pair = match parse_read_pair(r1_seq, r1_qual, r2_seq, r2_qual, &config_parser) {
///     ParseResult::Ok(p) => *p,
///     ParseResult::Dropped { .. } => panic!("unexpected drop"),
/// };
///
/// let (molecules, stats) = assemble_molecules(vec![pair], &AssemblerConfig::default());
/// assert_eq!(molecules.len(), 1);
/// assert_eq!(stats.n_singletons, 1);
/// ```
pub fn assemble_molecules(
    read_pairs: Vec<ParsedReadPair>,
    config: &AssemblerConfig,
) -> (Vec<Molecule>, AssemblyStats) {
    let mut stats = AssemblyStats::default();

    if read_pairs.is_empty() {
        return (Vec::new(), stats);
    }

    // ── Step 1: Collect all unique canonical UMIs with read counts ────────────
    // Hash-partition grouping: use a HashMap for O(1) UMI lookup instead of
    // the O(n) linear scan that made the original O(n²) over all read pairs.
    // This reduces step-1 complexity from O(n × u) to O(n), allowing processing
    // of 10–20M read pairs within a practical time budget.
    //
    // `unique_umis[i]` = (canonical_umi, read_count)
    // `read_to_umi[i]` = index into unique_umis for read_pairs[i]
    let mut umi_to_index: AHashMap<CanonicalUmiPair, usize> =
        AHashMap::with_capacity(read_pairs.len());
    let mut unique_umis: Vec<(CanonicalUmiPair, u32)> = Vec::new();
    let mut read_to_umi: Vec<usize> = vec![0; read_pairs.len()];

    for (idx, pair) in read_pairs.iter().enumerate() {
        let i = if let Some(&i) = umi_to_index.get(&pair.canonical_umi) {
            unique_umis[i].1 += 1;
            i
        } else {
            let i = unique_umis.len();
            umi_to_index.insert(pair.canonical_umi.clone(), i);
            unique_umis.push((pair.canonical_umi.clone(), 1));
            i
        };
        read_to_umi[idx] = i;
    }
    drop(umi_to_index); // free the HashMap before the O(u²) clustering step

    // ── Step 2: Hash-partition Hamming clustering ─────────────────────────────
    // Partition unique UMIs into 64 buckets by the first 3 bases of umi_a,
    // then cluster within each bucket.  This reduces the O(u²) clustering step
    // to O(u²/64) — a ~64× speedup at high depth — at the cost of missing
    // Hamming-1 merges where the mismatch falls in umi_a[0..3] (~3/10 of all
    // within-distance-1 pairs).  At practical depths (≤2M reads), the miss
    // rate is negligible: the affected pairs remain as separate clusters and
    // each still contributes independently to variant evidence.
    // When every UMI is unique (high molecule count, low grouping ratio),
    // Hamming clustering won't merge anything useful and costs O(u²/64).
    // Skip clustering when unique_umis > 200K.
    let effective_distance = if unique_umis.len() > 200_000 {
        0
    } else {
        config.max_hamming_distance
    };
    let clusters = partition_and_cluster(&unique_umis, effective_distance);

    // ── Step 3: Assign reads to clusters via inverted index ───────────────────
    // Avoids the O(n × c) scan of all reads for each cluster. Instead, build
    // a umi-to-cluster mapping and then assign each read in a single O(n) pass.
    let mut umi_to_cluster = vec![0usize; unique_umis.len()];
    for (cluster_idx, cluster_umi_indices) in clusters.iter().enumerate() {
        for &umi_idx in cluster_umi_indices {
            umi_to_cluster[umi_idx] = cluster_idx;
        }
    }

    let mut cluster_reads: Vec<Vec<usize>> = vec![Vec::new(); clusters.len()];
    for (read_idx, &umi_idx) in read_to_umi.iter().enumerate() {
        cluster_reads[umi_to_cluster[umi_idx]].push(read_idx);
    }

    // ── Step 4: Build molecules from each cluster's reads ────────────────────
    let mut molecules: Vec<Molecule> = Vec::new();

    for cluster_read_indices in &cluster_reads {
        if cluster_read_indices.is_empty() {
            continue;
        }

        // Sub-group by fingerprint compatibility to detect UMI collisions.
        let fingerprint_groups = split_by_fingerprint(cluster_read_indices, &read_pairs);

        // Count extra splits as collisions.
        if fingerprint_groups.len() > 1 {
            stats.n_umi_collisions_detected += (fingerprint_groups.len() - 1) as u64;
        }

        // Build a molecule for each fingerprint sub-group.
        for group_indices in fingerprint_groups {
            if let Some(mol) = build_molecule(&group_indices, &read_pairs, config, &mut stats) {
                molecules.push(mol);
            }
        }
    }

    // Sort molecules by id for deterministic output order (avoids HashMap
    // iteration non-determinism in output).
    molecules.sort_by_key(|m| m.id);

    stats.n_molecules = molecules.len() as u64;

    (molecules, stats)
}

// ── Internal helpers ──────────────────────────────────────────────────────────

/// Partition unique UMIs by the first 3 bases of `umi_a` (64 buckets), then
/// run Hamming-distance clustering within each bucket.
///
/// Two UMI pairs with a Hamming-distance-1 mismatch in `umi_a[0..3]` will
/// land in different buckets and will not be merged.  This is an acceptable
/// trade-off: ~3/10 of Hamming-1 pairs are missed, but the affected pairs
/// remain as separate clusters that each contribute independently to variant
/// evidence.  The speedup is ~64× for the O(u²) clustering step, enabling
/// processing of 10–20M read pairs within a practical time budget.
///
/// Returns `Vec<Vec<usize>>` where each inner vec is a cluster of indices into
/// the original `pairs` slice (same format as [`cluster_umi_pairs`]).
fn partition_and_cluster(pairs: &[(CanonicalUmiPair, u32)], max_distance: u32) -> Vec<Vec<usize>> {
    // Bucket global unique-UMI indices by their umi_a prefix (first 3 bases).
    let mut buckets: AHashMap<[u8; 3], Vec<usize>> = AHashMap::new();
    for (i, (umi, _count)) in pairs.iter().enumerate() {
        let key = [umi.umi_a[0], umi.umi_a[1], umi.umi_a[2]];
        buckets.entry(key).or_default().push(i);
    }

    let mut all_clusters: Vec<Vec<usize>> = Vec::new();

    for (_key, global_indices) in buckets {
        // Build a local pairs slice for this bucket.
        let local_pairs: Vec<(CanonicalUmiPair, u32)> =
            global_indices.iter().map(|&i| pairs[i].clone()).collect();

        // Cluster within this bucket — returned indices are into local_pairs.
        let local_clusters = cluster_umi_pairs(&local_pairs, max_distance);

        // Remap local indices to global unique-UMI indices.
        for local_cluster in local_clusters {
            let global_cluster: Vec<usize> = local_cluster
                .iter()
                .map(|&local_idx| global_indices[local_idx])
                .collect();
            all_clusters.push(global_cluster);
        }
    }

    all_clusters
}

/// Split a set of read-pair indices into sub-groups where all members have
/// mutually compatible endpoint fingerprints.
///
/// Uses a greedy approach: each read is added to the first existing group whose
/// representative fingerprint is compatible with the read's fingerprint; if none
/// is found, a new group is started.
fn split_by_fingerprint(indices: &[usize], read_pairs: &[ParsedReadPair]) -> Vec<Vec<usize>> {
    // (representative_fingerprint, members)
    let mut groups: Vec<(u64, Vec<usize>)> = Vec::new();

    for &idx in indices {
        let rp = &read_pairs[idx];
        let fp = compute_endpoint_fingerprint(&rp.template_r1, &rp.template_r2);

        let found = groups
            .iter_mut()
            .find(|(rep_fp, _)| fingerprints_duplex_compatible(*rep_fp, fp));

        match found {
            Some((_, members)) => members.push(idx),
            None => groups.push((fp, vec![idx])),
        }
    }

    groups.into_iter().map(|(_, members)| members).collect()
}

/// Build a [`Molecule`] from a fingerprint-compatible group of read pairs,
/// updating `stats` in place.
///
/// Returns `None` if the family is below the minimum size threshold.
fn build_molecule(
    indices: &[usize],
    read_pairs: &[ParsedReadPair],
    config: &AssemblerConfig,
    stats: &mut AssemblyStats,
) -> Option<Molecule> {
    // Separate into forward and reverse strand reads.
    let mut fwd_reads: Vec<&ParsedReadPair> = Vec::new();
    let mut rev_reads: Vec<&ParsedReadPair> = Vec::new();

    for &idx in indices {
        let rp = &read_pairs[idx];
        match rp.strand {
            Strand::Forward => fwd_reads.push(rp),
            Strand::Reverse => rev_reads.push(rp),
        }
    }

    // Saturating cast: families larger than 255 reads are capped at 255.
    // Silent truncation via `as u8` would wrap, giving an incorrect count.
    let n_fwd = (fwd_reads.len().min(u8::MAX as usize)) as u8;
    let n_rev = (rev_reads.len().min(u8::MAX as usize)) as u8;
    let total = ((fwd_reads.len() + rev_reads.len()).min(u8::MAX as usize)) as u8;

    // ── min_family_size filter ─────────────────────────────────────────────
    if total < config.min_family_size {
        stats.n_families_below_min_size += 1;
        return None;
    }

    // ── Derive canonical UMI from the first read in the group ─────────────
    let canonical_umi = &read_pairs[indices[0]].canonical_umi;
    let umi_fwd = canonical_umi.umi_a.clone();
    let umi_rev = canonical_umi.umi_b.clone();

    // ── Single-strand consensus (forward) ─────────────────────────────────
    // Call SSC whenever there are forward reads, regardless of min_duplex_reads.
    // min_duplex_reads only gates whether both strands qualify for duplex consensus.
    let consensus_fwd = if !fwd_reads.is_empty() {
        call_ssc(&fwd_reads, n_fwd, n_rev, config, true)
    } else {
        None
    };

    // ── Single-strand consensus (reverse) ─────────────────────────────────
    let consensus_rev = if !rev_reads.is_empty() {
        call_ssc(&rev_reads, n_fwd, n_rev, config, false)
    } else {
        None
    };

    // ── Duplex consensus ───────────────────────────────────────────────────
    let duplex_consensus = if n_fwd >= config.min_duplex_reads && n_rev >= config.min_duplex_reads {
        match (&consensus_fwd, &consensus_rev) {
            (Some(fwd_cr), Some(rev_cr)) => {
                // Re-derive ConsensusBase vecs by re-calling SSC (we need
                // ConsensusBase, not ConsensusRead, for duplex_consensus).
                let fwd_bases = call_ssc_bases(&fwd_reads, config, true);
                let rev_bases = call_ssc_bases(&rev_reads, config, false);
                match (fwd_bases, rev_bases) {
                    (Some(fb), Some(rb)) => {
                        crate::consensus::duplex_consensus(&fb, &rb, config.disagreement_strategy)
                            .map(|duplex_bases| {
                                let sequence: Vec<u8> =
                                    duplex_bases.iter().map(|d| d.base).collect();
                                let per_base_error_prob: Vec<f32> =
                                    duplex_bases.iter().map(|d| d.error_prob).collect();
                                // Both fwd and rev depths per position.
                                let per_base_strand_support: Vec<(u8, u8)> = duplex_bases
                                    .iter()
                                    .map(|d| (d.fwd_depth, d.rev_depth))
                                    .collect();
                                // Suppress unused variable warnings — fwd_cr and
                                // rev_cr are used only to confirm both SSCs succeeded.
                                let _ = (fwd_cr, rev_cr);
                                ConsensusRead {
                                    sequence,
                                    per_base_error_prob,
                                    per_base_strand_support,
                                    family_size: (n_fwd, n_rev),
                                }
                            })
                    }
                    _ => None,
                }
            }
            _ => None,
        }
    } else {
        None
    };

    // ── Classify family type ───────────────────────────────────────────────
    let family_type = FamilyType::from_family_size((n_fwd, n_rev));

    match family_type {
        FamilyType::Duplex => stats.n_duplex += 1,
        FamilyType::SimplexFwd => stats.n_simplex_fwd += 1,
        FamilyType::SimplexRev => stats.n_simplex_rev += 1,
        FamilyType::Singleton => stats.n_singletons += 1,
    }

    // ── Stable molecule id from canonical UMI pair hash ───────────────────
    let id = hash_umi_pair(canonical_umi);

    Some(Molecule {
        id,
        umi_fwd,
        umi_rev,
        consensus_fwd,
        consensus_rev,
        duplex_consensus,
        evidence: None,
    })
}

/// Call single-strand consensus and convert the result to a [`ConsensusRead`].
///
/// Returns `None` when the read list is empty or SSC returns no result.
fn call_ssc(
    reads: &[&ParsedReadPair],
    n_fwd: u8,
    n_rev: u8,
    config: &AssemblerConfig,
    is_forward: bool,
) -> Option<ConsensusRead> {
    if reads.is_empty() {
        return None;
    }

    // Collect sequence and quality slices for the relevant template.
    // For a forward-strand read we use template_r1; for reverse, template_r2.
    // This matches the orientation the read came in on.
    let seqs: Vec<&[u8]> = reads
        .iter()
        .map(|rp| {
            if is_forward {
                rp.template_r1.as_slice()
            } else {
                rp.template_r2.as_slice()
            }
        })
        .collect();
    let quals: Vec<&[u8]> = reads
        .iter()
        .map(|rp| {
            if is_forward {
                rp.qual_r1.as_slice()
            } else {
                rp.qual_r2.as_slice()
            }
        })
        .collect();

    let bases = single_strand_consensus(&seqs, &quals, &config.consensus)?;

    let sequence: Vec<u8> = bases.iter().map(|b| b.base).collect();
    let per_base_error_prob: Vec<f32> = bases.iter().map(|b| b.error_prob).collect();
    let per_base_strand_support: Vec<(u8, u8)> = bases
        .iter()
        .map(|b| {
            if is_forward {
                (b.depth, 0)
            } else {
                (0, b.depth)
            }
        })
        .collect();

    Some(ConsensusRead {
        sequence,
        per_base_error_prob,
        per_base_strand_support,
        family_size: (n_fwd, n_rev),
    })
}

/// Call single-strand consensus and return the raw [`ConsensusBase`] vec
/// (used internally to feed [`crate::consensus::duplex_consensus`]).
///
/// `is_forward` selects the correct template: `template_r1` for forward reads,
/// `template_r2` for reverse reads.  This mirrors the selection in [`call_ssc`]
/// and ensures duplex consensus compares R1-derived bases against R2-derived
/// bases rather than R1 against R1.
fn call_ssc_bases(
    reads: &[&ParsedReadPair],
    config: &AssemblerConfig,
    is_forward: bool,
) -> Option<Vec<crate::consensus::ConsensusBase>> {
    if reads.is_empty() {
        return None;
    }

    let seqs: Vec<&[u8]> = reads
        .iter()
        .map(|rp| {
            if is_forward {
                rp.template_r1.as_slice()
            } else {
                rp.template_r2.as_slice()
            }
        })
        .collect();
    let quals: Vec<&[u8]> = reads
        .iter()
        .map(|rp| {
            if is_forward {
                rp.qual_r1.as_slice()
            } else {
                rp.qual_r2.as_slice()
            }
        })
        .collect();

    single_strand_consensus(&seqs, &quals, &config.consensus)
}

/// Compute a stable 64-bit molecule id from a canonical UMI pair.
///
// Uses FNV-1a: streams umi_a, a 0xFF separator (which cannot appear in a
/// valid UMI), then umi_b.  The separator prevents "ACGT" + "TTTTT" from
/// colliding with "ACGTT" + "TTTT" if UMI lengths ever differ.
/// Supports UMIs of any length and produces deterministic output across
/// Rust versions and platforms.
fn hash_umi_pair(pair: &CanonicalUmiPair) -> u64 {
    const OFFSET_BASIS: u64 = 14_695_981_039_346_656_037;
    const PRIME: u64 = 1_099_511_628_211;

    let mut hash = OFFSET_BASIS;
    for &byte in &pair.umi_a {
        hash ^= byte as u64;
        hash = hash.wrapping_mul(PRIME);
    }
    // Separator — 0xFF is not a valid ASCII base, so it cannot appear in a UMI.
    hash ^= 0xFF_u64;
    hash = hash.wrapping_mul(PRIME);
    for &byte in &pair.umi_b {
        hash ^= byte as u64;
        hash = hash.wrapping_mul(PRIME);
    }
    hash
}

// ── Tests ─────────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;
    use crate::parser::{parse_read_pair, ParseResult, ParserConfig};

    // ── Test 0: hash_umi_pair produces a stable, deterministic value ─────────
    // A fixed input must always produce the same output. If this fails,
    // the hash function has changed and molecule IDs will not be reproducible.

    #[test]
    fn hash_umi_pair_is_stable() {
        let pair_a = CanonicalUmiPair::new(b"ACGTA".to_vec(), b"TGCAT".to_vec());
        let pair_b = CanonicalUmiPair::new(b"AAAAA".to_vec(), b"TTTTT".to_vec());
        // Two distinct pairs must produce different hashes.
        assert_ne!(hash_umi_pair(&pair_a), hash_umi_pair(&pair_b));
        // The same pair must always produce the same hash.
        assert_eq!(hash_umi_pair(&pair_a), hash_umi_pair(&pair_a));
    }

    /// Build a [`ParsedReadPair`] with the given UMIs and template sequences.
    ///
    /// Template length must be at least 8 bases so fingerprinting is
    /// meaningful.  Both templates are set to the same value for simplicity.
    /// UMIs may be any length.
    fn make_pair(umi_r1: &[u8], umi_r2: &[u8], template: &[u8]) -> ParsedReadPair {
        let mut r1_seq = Vec::new();
        r1_seq.extend_from_slice(umi_r1);
        r1_seq.extend_from_slice(b"TG"); // skip
        r1_seq.extend_from_slice(template);

        let mut r2_seq = Vec::new();
        r2_seq.extend_from_slice(umi_r2);
        r2_seq.extend_from_slice(b"TG"); // skip
        r2_seq.extend_from_slice(template);

        let r1_qual = vec![b'I'; r1_seq.len()]; // Phred 40
        let r2_qual = vec![b'I'; r2_seq.len()];

        let config = ParserConfig::default();
        match parse_read_pair(&r1_seq, &r1_qual, &r2_seq, &r2_qual, &config) {
            ParseResult::Ok(p) => *p,
            ParseResult::Dropped { reason, detail } => {
                panic!("make_pair produced Dropped({reason:?}): {detail}");
            }
        }
    }

    // ── Test 1: Single read pair → singleton molecule ─────────────────────

    #[test]
    fn single_read_pair_is_singleton() {
        let pair = make_pair(b"ACGTA", b"TGCAT", b"NNNNNNNNNN");
        let (molecules, stats) = assemble_molecules(vec![pair], &AssemblerConfig::default());

        assert_eq!(molecules.len(), 1, "expected exactly one molecule");
        let mol = &molecules[0];
        assert_eq!(mol.umi_fwd, b"ACGTA");
        assert_eq!(mol.umi_rev, b"TGCAT");
        assert_eq!(stats.n_singletons, 1);
        assert_eq!(stats.n_duplex, 0);
    }

    // ── Test 2: Two read pairs with same UMI, same strand → simplex family ──

    #[test]
    fn two_same_strand_reads_form_simplex() {
        // Both reads have ACGTA as R1 UMI and TGCAT as R2 UMI.
        // "ACGTA" < "TGCAT" → both are Forward strand reads.
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAA");
        let p2 = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAA");

        let (molecules, stats) = assemble_molecules(vec![p1, p2], &AssemblerConfig::default());

        assert_eq!(molecules.len(), 1);
        let mol = &molecules[0];
        assert!(mol.consensus_fwd.is_some(), "expected forward consensus");
        assert!(mol.consensus_rev.is_none(), "no reverse reads");
        assert!(
            mol.duplex_consensus.is_none(),
            "cannot be duplex with one strand"
        );
        assert_eq!(stats.n_simplex_fwd, 1);
        assert_eq!(stats.n_duplex, 0);
    }

    // ── Test 3: Two read pairs with same UMI, opposite strands → duplex ─────

    #[test]
    fn opposite_strand_reads_form_duplex() {
        // Forward strand read: R1 UMI = canonical umi_a = ACGTA
        let fwd = make_pair(b"ACGTA", b"TGCAT", b"ACGTACGTACGTACGTACGT");
        // Reverse strand read: R1 UMI = canonical umi_b = TGCAT
        let rev = make_pair(b"TGCAT", b"ACGTA", b"ACGTACGTACGTACGTACGT");

        let (molecules, stats) = assemble_molecules(vec![fwd, rev], &AssemblerConfig::default());

        assert_eq!(molecules.len(), 1, "should collapse to one molecule");
        let mol = &molecules[0];
        assert!(mol.consensus_fwd.is_some(), "expected fwd consensus");
        assert!(mol.consensus_rev.is_some(), "expected rev consensus");
        assert!(mol.duplex_consensus.is_some(), "expected duplex consensus");
        assert_eq!(stats.n_duplex, 1);
        assert_eq!(stats.n_singletons, 0);
    }

    // ── Test 4: Same UMI, different fingerprints → two molecules ─────────────

    #[test]
    fn different_fingerprints_split_into_two_molecules() {
        // Both pairs share the same canonical UMI.
        // template_a is all A's; template_b is all T's — these fingerprints are
        // maximally different (all bits flip: A=00, T=11 → 64 bits XOR = all 1s).
        let p_a = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");
        let p_b = make_pair(b"ACGTA", b"TGCAT", b"TTTTTTTTTTTTTTTTTTTTTTTTT");

        let (molecules, stats) = assemble_molecules(vec![p_a, p_b], &AssemblerConfig::default());

        assert_eq!(
            molecules.len(),
            2,
            "fingerprint split should yield two molecules"
        );
        assert_eq!(stats.n_umi_collisions_detected, 1);
    }

    // ── Test 5: UMI Hamming distance 1 → merged into one molecule ────────────

    #[test]
    fn umi_hamming_distance_one_merged() {
        // Two pairs: UMIs differ by one base in umi_a.
        // "ACGTA" and "ACGTT" differ at position 4 (A→T) — Hamming dist = 1.
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");
        let p2 = make_pair(b"ACGTT", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");

        let config = AssemblerConfig {
            max_hamming_distance: 1,
            ..AssemblerConfig::default()
        };
        let (molecules, _stats) = assemble_molecules(vec![p1, p2], &config);

        // Both reads carry "AAAA..." templates → same fingerprint group → one molecule.
        assert_eq!(molecules.len(), 1, "Hamming-1 UMIs should be merged");
    }

    // ── Test 6: min_family_size filter ────────────────────────────────────────

    #[test]
    fn min_family_size_filters_singletons() {
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"NNNNNNNNNNNNNNNNNNNNNNNNNN");
        let p2 = make_pair(b"GGGGG", b"CCCCC", b"NNNNNNNNNNNNNNNNNNNNNNNNNN");

        // min_family_size = 2 → both singletons should be dropped.
        let config = AssemblerConfig {
            min_family_size: 2,
            ..AssemblerConfig::default()
        };
        let (molecules, stats) = assemble_molecules(vec![p1, p2], &config);

        assert_eq!(
            molecules.len(),
            0,
            "both families are below min_family_size=2"
        );
        assert_eq!(stats.n_families_below_min_size, 2);
    }

    // ── Test 7: AssemblyStats counts are correct ──────────────────────────────

    #[test]
    fn assembly_stats_counts_correct() {
        // Three distinct UMI pairs.  We construct:
        //  - One duplex molecule (fwd + rev)
        //  - One simplex-fwd molecule (two fwd reads, same UMI)
        //  - One singleton

        // Duplex pair
        let d_fwd = make_pair(b"ACGTA", b"TGCAT", b"ACGTACGTACGTACGTACGT");
        let d_rev = make_pair(b"TGCAT", b"ACGTA", b"ACGTACGTACGTACGTACGT");

        // Simplex-fwd pair (GGGGG < TTTTT, so both are Forward strand)
        let s1 = make_pair(b"GGGGG", b"TTTTT", b"CCCCCCCCCCCCCCCCCCCC");
        let s2 = make_pair(b"GGGGG", b"TTTTT", b"CCCCCCCCCCCCCCCCCCCC");

        // Singleton
        let single = make_pair(b"CCCCC", b"AAAAA", b"GGGGGGGGGGGGGGGGGGGG");
        // Note: "AAAAA" < "CCCCC" so canonical is (AAAAA, CCCCC);
        // R1=CCCCC → Reverse strand → singleton on rev.

        let (molecules, stats) = assemble_molecules(
            vec![d_fwd, d_rev, s1, s2, single],
            &AssemblerConfig::default(),
        );

        assert_eq!(molecules.len(), 3, "expected 3 molecules");
        assert_eq!(stats.n_molecules, 3);
        assert_eq!(stats.n_duplex, 1);
        assert_eq!(stats.n_simplex_fwd, 1);
        assert_eq!(stats.n_singletons, 1);
    }

    // ── Test 8: Genuine duplex pair with rotated fingerprints → one molecule ──
    // Before the rotation fix, forward and reverse reads from the same molecule
    // would fail fingerprint compatibility and be split into two separate
    // molecules, suppressing duplex calling.  The duplex-aware check must group
    // them into a single duplex molecule.

    #[test]
    fn duplex_pair_with_rotated_fingerprints_forms_one_molecule() {
        use kam_core::molecule::{CanonicalUmiPair, Strand};

        // For a genuine duplex pair the forward read has
        //   template_r1 = T,      template_r2 = RC(T)
        // and the reverse read has
        //   template_r1 = RC(T),  template_r2 = T.
        //
        // This gives fp_fwd and fp_rev = rotate_left(fp_fwd, 32), so direct
        // fingerprint comparison fails (~32 bits differ).
        //
        // Use T = AAAAAAAAGGGGGGGG and RC(T) = CCCCCCCCTTTTTTTT to get clearly
        // distinct templates whose fingerprint XOR has exactly 32 bits set.

        let t: Vec<u8> = b"AAAAAAAAGGGGGGGG".to_vec();
        let rc_t: Vec<u8> = b"CCCCCCCCTTTTTTTT".to_vec();

        // UMIs: ACGTA < TGCAT → canonical = (ACGTA, TGCAT).
        // Forward read: R1_UMI = umi_a = ACGTA → strand = Forward.
        // Reverse read: R1_UMI = umi_b = TGCAT → strand = Reverse.
        let canonical_umi = CanonicalUmiPair::new(b"ACGTA".to_vec(), b"TGCAT".to_vec());
        let qual = vec![b'I'; t.len()];

        let fwd_pair = crate::parser::ParsedReadPair {
            umi_r1: b"ACGTA".to_vec(),
            umi_r2: b"TGCAT".to_vec(),
            skip_r1: b"TG".to_vec(),
            skip_r2: b"TG".to_vec(),
            template_r1: t.clone(),
            template_r2: rc_t.clone(),
            qual_r1: qual.clone(),
            qual_r2: qual.clone(),
            umi_qual_r1: vec![b'I'; 5],
            umi_qual_r2: vec![b'I'; 5],
            canonical_umi: canonical_umi.clone(),
            strand: Strand::Forward,
        };

        let rev_pair = crate::parser::ParsedReadPair {
            umi_r1: b"TGCAT".to_vec(),
            umi_r2: b"ACGTA".to_vec(),
            skip_r1: b"TG".to_vec(),
            skip_r2: b"TG".to_vec(),
            template_r1: rc_t.clone(),
            template_r2: t.clone(),
            qual_r1: qual.clone(),
            qual_r2: qual.clone(),
            umi_qual_r1: vec![b'I'; 5],
            umi_qual_r2: vec![b'I'; 5],
            canonical_umi: canonical_umi.clone(),
            strand: Strand::Reverse,
        };

        let (molecules, stats) =
            assemble_molecules(vec![fwd_pair, rev_pair], &AssemblerConfig::default());

        assert_eq!(
            molecules.len(),
            1,
            "genuine duplex pair must be grouped into one molecule, not split"
        );
        assert_eq!(stats.n_duplex, 1, "should be classified as duplex");
        assert_eq!(
            stats.n_umi_collisions_detected, 0,
            "no UMI collision should be reported for a genuine duplex pair"
        );
    }

    // ── Test 9: Duplex consensus uses R1 for forward SSC, R2 for reverse SSC ──
    // This is a regression test for the bug where call_ssc_bases always read
    // template_r1, making duplex consensus compare R1 against R1 instead of
    // R1 against R2.  The test gives each read deliberately distinct R1 and R2
    // template sequences so the wrong selection would produce a different duplex
    // consensus sequence.
    //
    // Setup:
    //   fwd read: template_r1 = AAAAAAAAAAAAAAAA, template_r2 = CCCCCCCCCCCCCCCC
    //   rev read: template_r1 = CCCCCCCCCCCCCCCC, template_r2 = AAAAAAAAAAAAAAAA
    //
    // Correct: fwd SSC uses R1 → "AAAA...", rev SSC uses R2 → "AAAA...".
    //          Both strands agree → duplex consensus = "AAAA...".
    //
    // Bug: fwd SSC uses R1 → "AAAA...", rev SSC also uses R1 → "CCCC...".
    //      Both strands disagree → duplex consensus would be N-masked or absent.

    #[test]
    fn duplex_consensus_uses_r1_for_fwd_and_r2_for_rev() {
        use kam_core::molecule::{CanonicalUmiPair, Strand};

        let fwd_template: Vec<u8> = b"AAAAAAAAAAAAAAAA".to_vec();
        let rev_template: Vec<u8> = b"CCCCCCCCCCCCCCCC".to_vec();
        let qual = vec![b'I'; fwd_template.len()]; // Phred 40

        let canonical_umi = CanonicalUmiPair::new(b"ACGTA".to_vec(), b"TGCAT".to_vec());

        // Forward read: R1_UMI = umi_a = ACGTA → strand = Forward.
        // template_r1 = "AAAA..." (the molecule's sense strand template).
        // template_r2 = "CCCC..." (the antisense template, present on this read).
        let fwd_pair = crate::parser::ParsedReadPair {
            umi_r1: b"ACGTA".to_vec(),
            umi_r2: b"TGCAT".to_vec(),
            skip_r1: b"TG".to_vec(),
            skip_r2: b"TG".to_vec(),
            template_r1: fwd_template.clone(),
            template_r2: rev_template.clone(),
            qual_r1: qual.clone(),
            qual_r2: qual.clone(),
            umi_qual_r1: vec![b'I'; 5],
            umi_qual_r2: vec![b'I'; 5],
            canonical_umi: canonical_umi.clone(),
            strand: Strand::Forward,
        };

        // Reverse read: R1_UMI = umi_b = TGCAT → strand = Reverse.
        // template_r1 = "CCCC..." (antisense strand, as seen on this read's R1).
        // template_r2 = "AAAA..." (sense strand, as seen on this read's R2).
        let rev_pair = crate::parser::ParsedReadPair {
            umi_r1: b"TGCAT".to_vec(),
            umi_r2: b"ACGTA".to_vec(),
            skip_r1: b"TG".to_vec(),
            skip_r2: b"TG".to_vec(),
            template_r1: rev_template.clone(),
            template_r2: fwd_template.clone(),
            qual_r1: qual.clone(),
            qual_r2: qual.clone(),
            umi_qual_r1: vec![b'I'; 5],
            umi_qual_r2: vec![b'I'; 5],
            canonical_umi: canonical_umi.clone(),
            strand: Strand::Reverse,
        };

        let (molecules, stats) =
            assemble_molecules(vec![fwd_pair, rev_pair], &AssemblerConfig::default());

        assert_eq!(molecules.len(), 1, "should form one duplex molecule");
        assert_eq!(stats.n_duplex, 1);

        let mol = &molecules[0];
        let duplex = mol
            .duplex_consensus
            .as_ref()
            .expect("duplex consensus must be present");

        // With the correct template selection, fwd uses R1 = "AAAA..." and rev
        // uses R2 = "AAAA..." — both agree, so every base must be 'A'.
        assert!(
            duplex.sequence.iter().all(|&b| b == b'A'),
            "duplex consensus must be all-A (both strands agree on 'A'); \
             got {:?} — likely template_r1 was used for the reverse SSC instead of template_r2",
            std::str::from_utf8(&duplex.sequence).unwrap_or("<invalid utf8>"),
        );
    }

    // ── Test 10: Empty input returns empty vec ────────────────────────────────

    #[test]
    fn empty_input_returns_empty() {
        let (molecules, stats) = assemble_molecules(vec![], &AssemblerConfig::default());
        assert!(molecules.is_empty());
        assert_eq!(stats.n_molecules, 0);
        assert_eq!(stats.n_duplex, 0);
        assert_eq!(stats.n_singletons, 0);
    }

    // ── Test 11: 12 bp UMI — full assemble pipeline works ────────────────────
    // Exercises the arbitrary-UMI-length path introduced in CHEM-003.

    #[test]
    fn twelve_bp_umi_assembles_correctly() {
        use kam_core::chemistry::ReadStructure;

        let config_parser = ParserConfig {
            read_structure: ReadStructure {
                umi_length: 12,
                skip_length: 2,
            },
            ..ParserConfig::default()
        };

        // Build two reads from opposite strands of the same molecule.
        // R1 UMI forward: "ACGTACGTACGT" < "TGCATGCATGCA" → umi_a = first one → Forward.
        // R1 UMI reverse: "TGCATGCATGCA" → Reverse.
        let template = b"AAAAAAAAAAAAAAAA"; // 16 bp — long enough for fingerprinting

        let mut r1_fwd = b"ACGTACGTACGT".to_vec(); // 12 bp UMI
        r1_fwd.extend_from_slice(b"TG"); // skip
        r1_fwd.extend_from_slice(template);
        let q_fwd = vec![b'I'; r1_fwd.len()];

        let mut r2_fwd = b"TGCATGCATGCA".to_vec();
        r2_fwd.extend_from_slice(b"TG");
        r2_fwd.extend_from_slice(template);
        let q2_fwd = vec![b'I'; r2_fwd.len()];

        let fwd_pair = match parse_read_pair(&r1_fwd, &q_fwd, &r2_fwd, &q2_fwd, &config_parser) {
            ParseResult::Ok(p) => *p,
            ParseResult::Dropped { reason, detail } => {
                panic!("Dropped({reason:?}): {detail}");
            }
        };

        let mut r1_rev = b"TGCATGCATGCA".to_vec();
        r1_rev.extend_from_slice(b"TG");
        r1_rev.extend_from_slice(template);
        let q_rev = vec![b'I'; r1_rev.len()];

        let mut r2_rev = b"ACGTACGTACGT".to_vec();
        r2_rev.extend_from_slice(b"TG");
        r2_rev.extend_from_slice(template);
        let q2_rev = vec![b'I'; r2_rev.len()];

        let rev_pair = match parse_read_pair(&r1_rev, &q_rev, &r2_rev, &q2_rev, &config_parser) {
            ParseResult::Ok(p) => *p,
            ParseResult::Dropped { reason, detail } => {
                panic!("Dropped({reason:?}): {detail}");
            }
        };

        let (molecules, stats) =
            assemble_molecules(vec![fwd_pair, rev_pair], &AssemblerConfig::default());

        assert_eq!(
            molecules.len(),
            1,
            "both reads should form one duplex molecule"
        );
        assert_eq!(stats.n_duplex, 1);

        let mol = &molecules[0];
        assert_eq!(mol.umi_fwd.len(), 12, "umi_fwd should be 12 bp");
        assert_eq!(mol.umi_rev.len(), 12, "umi_rev should be 12 bp");
        assert_eq!(mol.umi_fwd, b"ACGTACGTACGT");
        assert_eq!(mol.umi_rev, b"TGCATGCATGCA");
    }

    // ── New edge-case tests ──────────────────────────────────────────────────

    // Two reads with same canonical UMI, same fingerprint: form one molecule.
    #[test]
    fn same_umi_same_fingerprint_one_molecule() {
        let template = b"ACGTACGTACGTACGTACGT";
        let p1 = make_pair(b"GGGGA", b"TTTTT", template);
        let p2 = make_pair(b"GGGGA", b"TTTTT", template);
        let (molecules, stats) = assemble_molecules(vec![p1, p2], &AssemblerConfig::default());
        assert_eq!(
            molecules.len(),
            1,
            "same UMI + same fingerprint = one molecule"
        );
        assert_eq!(stats.n_umi_collisions_detected, 0);
    }

    // min_family_size = 1: singletons are kept.
    #[test]
    fn min_family_size_one_keeps_singletons() {
        let pair = make_pair(b"ACGTA", b"TGCAT", b"NNNNNNNNNNNNNNNNNNNNNNNNNN");
        let config = AssemblerConfig {
            min_family_size: 1,
            ..AssemblerConfig::default()
        };
        let (molecules, stats) = assemble_molecules(vec![pair], &config);
        assert_eq!(
            molecules.len(),
            1,
            "singleton should be kept with min_family_size=1"
        );
        assert_eq!(stats.n_families_below_min_size, 0);
    }

    // min_family_size = 2: a pair with 2 reads passes, singleton is filtered.
    #[test]
    fn min_family_size_two_pair_passes_singleton_filtered() {
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"ACGTACGTACGTACGTACGT");
        let p2 = make_pair(b"ACGTA", b"TGCAT", b"ACGTACGTACGTACGTACGT");
        let single = make_pair(b"GGGGG", b"CCCCC", b"NNNNNNNNNNNNNNNNNNNNNNNNNN");

        let config = AssemblerConfig {
            min_family_size: 2,
            ..AssemblerConfig::default()
        };
        let (molecules, stats) = assemble_molecules(vec![p1, p2, single], &config);
        assert_eq!(molecules.len(), 1, "only the family of 2 should survive");
        assert_eq!(
            stats.n_families_below_min_size, 1,
            "singleton should be filtered"
        );
    }

    // Hamming distance = 2 at threshold 1: UMIs are NOT merged.
    #[test]
    fn umi_hamming_distance_two_not_merged_at_threshold_one() {
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");
        let p2 = make_pair(b"ACGCC", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");

        let config = AssemblerConfig {
            max_hamming_distance: 1,
            ..AssemblerConfig::default()
        };
        let (molecules, _stats) = assemble_molecules(vec![p1, p2], &config);
        assert!(
            molecules.len() >= 2,
            "Hamming distance 2 at threshold 1 should not be merged: got {} molecules",
            molecules.len()
        );
    }

    // All reads have unique UMIs: each forms its own singleton.
    #[test]
    fn all_unique_umis_each_singleton() {
        let umis: Vec<(&[u8], &[u8])> = vec![
            (b"ACGTA", b"TGCAT"),
            (b"GGGGG", b"CCCCC"),
            (b"AAAAA", b"TTTTT"),
            (b"CCCCG", b"GGGGA"),
        ];
        let pairs: Vec<ParsedReadPair> = umis
            .into_iter()
            .map(|(u1, u2)| make_pair(u1, u2, b"NNNNNNNNNNNNNNNNNNNNNNNNNN"))
            .collect();
        let n = pairs.len();
        let (molecules, stats) = assemble_molecules(pairs, &AssemblerConfig::default());
        assert_eq!(
            molecules.len(),
            n,
            "all unique UMIs should produce one molecule each"
        );
        assert_eq!(stats.n_singletons, n as u64);
    }

    // Stress test: 100 reads with the same UMI. Should not panic.
    #[test]
    fn hundred_reads_same_umi_no_panic() {
        let template = b"ACGTACGTACGTACGTACGTACGTACGT";
        let pairs: Vec<ParsedReadPair> = (0..100)
            .map(|_| make_pair(b"ACGTA", b"TGCAT", template))
            .collect();
        let (molecules, stats) = assemble_molecules(pairs, &AssemblerConfig::default());
        assert_eq!(
            molecules.len(),
            1,
            "100 reads with same UMI should form one molecule"
        );
        assert_eq!(stats.n_molecules, 1);
    }

    // Hamming distance = 0 clustering: no merging even for distance-1 pairs.
    #[test]
    fn hamming_zero_no_merging() {
        let p1 = make_pair(b"ACGTA", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");
        let p2 = make_pair(b"ACGTT", b"TGCAT", b"AAAAAAAAAAAAAAAAAAAAAAAAA");

        let config = AssemblerConfig {
            max_hamming_distance: 0,
            ..AssemblerConfig::default()
        };
        let (molecules, _stats) = assemble_molecules(vec![p1, p2], &config);
        assert_eq!(
            molecules.len(),
            2,
            "Hamming distance 0 should not merge distance-1 pairs"
        );
    }

    // Molecule id is deterministic: same input always produces the same id.
    #[test]
    fn molecule_id_deterministic() {
        let template = b"ACGTACGTACGTACGTACGT";
        let p1 = make_pair(b"ACGTA", b"TGCAT", template);
        let p2 = make_pair(b"ACGTA", b"TGCAT", template);
        let (mols_a, _) = assemble_molecules(vec![p1], &AssemblerConfig::default());
        let (mols_b, _) = assemble_molecules(vec![p2], &AssemblerConfig::default());
        assert_eq!(
            mols_a[0].id, mols_b[0].id,
            "same canonical UMI must produce the same molecule id"
        );
    }

    // Molecules are sorted by id for deterministic output.
    #[test]
    fn output_sorted_by_molecule_id() {
        let pairs = vec![
            make_pair(b"ZZZZZ", b"YYYYY", b"NNNNNNNNNNNNNNNNNNNNNNNNNN"),
            make_pair(b"AAAAA", b"BBBBB", b"NNNNNNNNNNNNNNNNNNNNNNNNNN"),
            make_pair(b"GGGGG", b"HHHHH", b"NNNNNNNNNNNNNNNNNNNNNNNNNN"),
        ];
        let (molecules, _) = assemble_molecules(pairs, &AssemblerConfig::default());
        for window in molecules.windows(2) {
            assert!(
                window[0].id <= window[1].id,
                "molecules must be sorted by id: {} > {}",
                window[0].id,
                window[1].id
            );
        }
    }

    // Three reads: two form one strand, third forms the other strand. Duplex.
    #[test]
    fn three_reads_two_fwd_one_rev_duplex() {
        let template = b"ACGTACGTACGTACGTACGT";
        let fwd1 = make_pair(b"ACGTA", b"TGCAT", template);
        let fwd2 = make_pair(b"ACGTA", b"TGCAT", template);
        let rev = make_pair(b"TGCAT", b"ACGTA", template);

        let (molecules, stats) =
            assemble_molecules(vec![fwd1, fwd2, rev], &AssemblerConfig::default());
        assert_eq!(molecules.len(), 1, "should form one duplex molecule");
        assert_eq!(stats.n_duplex, 1);
        let mol = &molecules[0];
        assert!(mol.consensus_fwd.is_some());
        assert!(mol.consensus_rev.is_some());
        let fwd_cr = mol.consensus_fwd.as_ref().expect("fwd consensus");
        assert_eq!(fwd_cr.family_size, (2, 1));
    }

    // min_duplex_reads = 2: a pair with 1 fwd and 1 rev should NOT produce
    // duplex consensus.
    #[test]
    fn min_duplex_reads_two_suppresses_duplex_with_one_per_strand() {
        let template = b"ACGTACGTACGTACGTACGT";
        let fwd = make_pair(b"ACGTA", b"TGCAT", template);
        let rev = make_pair(b"TGCAT", b"ACGTA", template);

        let config = AssemblerConfig {
            min_duplex_reads: 2,
            ..AssemblerConfig::default()
        };
        let (molecules, _stats) = assemble_molecules(vec![fwd, rev], &config);
        assert_eq!(molecules.len(), 1);
        let mol = &molecules[0];
        assert!(
            mol.duplex_consensus.is_none(),
            "min_duplex_reads=2 should suppress duplex with only 1 read per strand"
        );
        assert!(mol.consensus_fwd.is_some());
        assert!(mol.consensus_rev.is_some());
    }
}
