diff --git a/Cargo.lock b/Cargo.lock index 28b615a68..b1c5b91f9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -849,6 +849,7 @@ dependencies = [ "criterion", "fgumi-dna", "fgumi-tag", + "memchr", "noodles", "proptest", "rstest", diff --git a/crates/fgumi-raw-bam/Cargo.toml b/crates/fgumi-raw-bam/Cargo.toml index 78492c5e1..d935aff3b 100644 --- a/crates/fgumi-raw-bam/Cargo.toml +++ b/crates/fgumi-raw-bam/Cargo.toml @@ -13,6 +13,7 @@ anyhow = { workspace = true, optional = true } bytemuck = { workspace = true } fgumi-dna = { workspace = true } fgumi-tag = { workspace = true } +memchr = { workspace = true } wide = { workspace = true } [features] diff --git a/crates/fgumi-raw-bam/src/cigar.rs b/crates/fgumi-raw-bam/src/cigar.rs index d80413ff4..199dfbaea 100644 --- a/crates/fgumi-raw-bam/src/cigar.rs +++ b/crates/fgumi-raw-bam/src/cigar.rs @@ -295,7 +295,7 @@ pub fn unclipped_5prime_raw(bam: &[u8], pos: i32, is_reverse: bool) -> i32 { /// For reverse strand mate: `unclipped_end` #[inline] #[must_use] -pub fn mate_unclipped_5prime(mate_pos: i32, mate_reverse: bool, mc_cigar: &str) -> i32 { +pub fn mate_unclipped_5prime(mate_pos: i32, mate_reverse: bool, mc_cigar: &[u8]) -> i32 { if mate_reverse { unclipped_other_end(mate_pos, mc_cigar) } else { @@ -311,7 +311,7 @@ pub fn mate_unclipped_5prime(mate_pos: i32, mate_reverse: bool, mc_cigar: &str) pub fn mate_unclipped_5prime_1based( mate_pos_0based: i32, mate_reverse: bool, - mc_cigar: &str, + mc_cigar: &[u8], ) -> i32 { let mate_pos_1based = mate_pos_0based.saturating_add(1); if mate_reverse { @@ -335,8 +335,8 @@ pub fn mate_unclipped_5prime_1based( /// `samtools sort --template-coordinate`. #[inline] #[must_use] -pub(crate) fn unclipped_other_start(mate_pos: i32, mc_cigar: &str) -> i32 { - mate_pos.saturating_sub(parse_leading_clips(mc_cigar.as_bytes())) +pub(crate) fn unclipped_other_start(mate_pos: i32, mc_cigar: &[u8]) -> i32 { + mate_pos.saturating_sub(parse_leading_clips(mc_cigar)) } /// Calculate mate's unclipped end from MC tag CIGAR string. @@ -349,8 +349,8 @@ pub(crate) fn unclipped_other_start(mate_pos: i32, mc_cigar: &str) -> i32 { /// the subtraction then walks the saturated value back. #[inline] #[must_use] -pub(crate) fn unclipped_other_end(mate_pos: i32, mc_cigar: &str) -> i32 { - let (ref_len, trailing_clips) = parse_ref_len_and_trailing_clips(mc_cigar.as_bytes()); +pub(crate) fn unclipped_other_end(mate_pos: i32, mc_cigar: &[u8]) -> i32 { + let (ref_len, trailing_clips) = parse_ref_len_and_trailing_clips(mc_cigar); mate_pos.saturating_sub(1).saturating_add(ref_len).saturating_add(trailing_clips) } @@ -1439,13 +1439,13 @@ mod tests { #[test] fn test_mate_unclipped_5prime_1based_forward() { // MC=5S10M: forward 5' = pos+1 - 5 = 96 - assert_eq!(mate_unclipped_5prime_1based(100, false, "5S10M"), 96); + assert_eq!(mate_unclipped_5prime_1based(100, false, b"5S10M"), 96); } #[test] fn test_mate_unclipped_5prime_1based_reverse() { // MC=10M5S: reverse 5' = pos+1 + 10 + 5 - 1 = 115 - assert_eq!(mate_unclipped_5prime_1based(100, true, "10M5S"), 115); + assert_eq!(mate_unclipped_5prime_1based(100, true, b"10M5S"), 115); } // ======================================================================== @@ -1506,52 +1506,52 @@ mod tests { /// -- for these inputs the old code panicked or wrapped, so it cannot serve /// as an oracle. #[rstest] - // Non-ASCII but valid UTF-8: `str::from_utf8` admits it, so it reaches the - // parser. Byte 4 falls inside 'e-acute', which panicked when the old code + // Non-ASCII bytes reach the parser directly now that `MC` is extracted as + // bytes. Byte 4 falls inside 'e-acute', which panicked when the old code // sliced the &str there. - #[case::non_ascii_utf8("10M\u{e9}", (10, 0))] - #[case::non_ascii_leading("\u{e9}10M", (10, 0))] + #[case::non_ascii_utf8("10M\u{e9}".as_bytes(), (10, 0))] + #[case::non_ascii_leading("\u{e9}10M".as_bytes(), (10, 0))] // Unknown ASCII operators are ignored, as before. - #[case::unknown_operator("10Mfoo5S", (10, 5))] - #[case::placeholder_star("*", (0, 0))] + #[case::unknown_operator(b"10Mfoo5S", (10, 5))] + #[case::placeholder_star(b"*", (0, 0))] // A single run wider than i32: `str::parse::` returns Err, and // `unwrap_or(0)` makes it 0. Saturation would wrongly give i32::MAX. - #[case::run_overflows_i32("99999999999M", (0, 0))] + #[case::run_overflows_i32(b"99999999999M", (0, 0))] // Totals that overflow across operators must clamp, not wrap to negative. - #[case::sum_overflows_i32("2000000000M2000000000M", (i32::MAX, 0))] - #[case::trailing_sum_overflows("1M2000000000S2000000000S", (1, i32::MAX))] + #[case::sum_overflows_i32(b"2000000000M2000000000M", (i32::MAX, 0))] + #[case::trailing_sum_overflows(b"1M2000000000S2000000000S", (1, i32::MAX))] // A non-ref, non-clip operator between the last ref op and the trailing // clips must not reset the run (matches htsjdk's getUnclippedEnd). - #[case::insertion_before_trailing_clip("10M5I3S", (10, 3))] - #[case::padding_before_trailing_clip("10M5P3S", (10, 3))] + #[case::insertion_before_trailing_clip(b"10M5I3S", (10, 3))] + #[case::padding_before_trailing_clip(b"10M5P3S", (10, 3))] // An unknown operator AFTER a trailing clip must not clear the run -- // only a reference-consuming operator resets it. - #[case::garbage_after_trailing_clip("10M3Sfoo", (10, 3))] + #[case::garbage_after_trailing_clip(b"10M3Sfoo", (10, 3))] // A run too wide for u64 itself: the accumulator must saturate, not wrap // back into i32 range. 18446744073709551617 == u64::MAX + 2. - #[case::run_overflows_u64("18446744073709551617M", (0, 0))] - #[case::leading_zeros("0007M", (7, 0))] - #[case::empty("", (0, 0))] - #[case::digits_only("10", (0, 0))] + #[case::run_overflows_u64(b"18446744073709551617M", (0, 0))] + #[case::leading_zeros(b"0007M", (7, 0))] + #[case::empty(b"", (0, 0))] + #[case::digits_only(b"10", (0, 0))] fn test_parse_ref_len_and_trailing_clips_malformed( - #[case] cigar: &str, + #[case] cigar: &[u8], #[case] expected: (i32, i32), ) { - assert_eq!(parse_ref_len_and_trailing_clips(cigar.as_bytes()), expected); + assert_eq!(parse_ref_len_and_trailing_clips(cigar), expected); } /// Same contract for the leading-clip parser. #[rstest] - #[case::non_ascii_utf8("\u{e9}10M", 0)] - #[case::run_overflows_i32("99999999999S", 0)] - #[case::run_overflows_u64("18446744073709551617S", 0)] - #[case::sum_overflows_i32("2000000000S2000000000S", i32::MAX)] - #[case::leading_zeros("0007S", 7)] - #[case::empty("", 0)] - #[case::digits_only("10", 0)] - #[case::stops_at_first_non_clip("5S10M5S", 5)] - fn test_parse_leading_clips_malformed(#[case] cigar: &str, #[case] expected: i32) { - assert_eq!(parse_leading_clips(cigar.as_bytes()), expected); + #[case::non_ascii_utf8("\u{e9}10M".as_bytes(), 0)] + #[case::run_overflows_i32(b"99999999999S", 0)] + #[case::run_overflows_u64(b"18446744073709551617S", 0)] + #[case::sum_overflows_i32(b"2000000000S2000000000S", i32::MAX)] + #[case::leading_zeros(b"0007S", 7)] + #[case::empty(b"", 0)] + #[case::digits_only(b"10", 0)] + #[case::stops_at_first_non_clip(b"5S10M5S", 5)] + fn test_parse_leading_clips_malformed(#[case] cigar: &[u8], #[case] expected: i32) { + assert_eq!(parse_leading_clips(cigar), expected); } /// One token of a generated MC value: a run of digits, a CIGAR operator, a @@ -1617,24 +1617,30 @@ mod tests { proptest::prop_assert!(trailing_clips >= 0); } - /// The public entry point must not panic either, on any MC that - /// survives `str::from_utf8` -- the filter these values pass through. + /// The public entry point must not panic either, on *any* byte string. + /// + /// This property used to be gated on `str::from_utf8` succeeding, + /// because the extractor discarded a non-UTF-8 `MC` before it could + /// reach here. These functions take CIGAR bytes now, so every generated + /// value is a real input and the gate would only hide the cases most + /// likely to break the parser. #[test] fn test_mate_unclipped_5prime_never_panics( cigar in cigar_ish_bytes(), mate_pos in proptest::prelude::any::(), mate_reverse in proptest::prelude::any::(), ) { - if let Ok(mc) = std::str::from_utf8(&cigar) { - // Only asserting it returns: the panic and the debug-build - // overflow are what this is guarding against. - let _ = mate_unclipped_5prime(mate_pos, mate_reverse, mc); - let _ = mate_unclipped_5prime_1based(mate_pos, mate_reverse, mc); - } + // Only asserting it returns: the panic and the debug-build + // overflow are what this is guarding against. + let _ = mate_unclipped_5prime(mate_pos, mate_reverse, &cigar); + let _ = mate_unclipped_5prime_1based(mate_pos, mate_reverse, &cigar); } } - /// The public entry points must not panic on a non-ASCII MC either. + /// The public entry points must not panic on a non-ASCII MC either. Since + /// these take CIGAR bytes, such a value now reaches the parser directly + /// rather than being discarded by a `str::from_utf8` gate first; the parser + /// reads the valid prefix and stops at the byte it cannot interpret. #[rstest] #[case::forward(false, 96)] #[case::reverse(true, 110)] @@ -1643,7 +1649,10 @@ mod tests { #[case] expected: i32, ) { // "5S10M\u{e9}": leading clips 5, ref_len 10, no trailing clips. - assert_eq!(mate_unclipped_5prime_1based(100, mate_reverse, "5S10M\u{e9}"), expected); + assert_eq!( + mate_unclipped_5prime_1based(100, mate_reverse, "5S10M\u{e9}".as_bytes()), + expected + ); } #[test] @@ -2430,25 +2439,25 @@ mod tests { #[test] fn test_unclipped_other_start_no_clips() { - assert_eq!(unclipped_other_start(100, "10M"), 100); + assert_eq!(unclipped_other_start(100, b"10M"), 100); } #[test] fn test_unclipped_other_end_no_trailing_clips() { // 10M: end = 100 + 10 + 0 - 1 = 109 - assert_eq!(unclipped_other_end(100, "10M"), 109); + assert_eq!(unclipped_other_end(100, b"10M"), 109); } #[test] fn test_unclipped_other_start_complex() { // 3H5S10M: leading = 3+5=8, start = 100 - 8 = 92 - assert_eq!(unclipped_other_start(100, "3H5S10M"), 92); + assert_eq!(unclipped_other_start(100, b"3H5S10M"), 92); } #[test] fn test_unclipped_other_end_complex() { // 10M5S3H: ref=10, trailing=5+3=8, end = 100 + 10 + 8 - 1 = 117 - assert_eq!(unclipped_other_end(100, "10M5S3H"), 117); + assert_eq!(unclipped_other_end(100, b"10M5S3H"), 117); } /// The inclusive-end `-1` must be applied before the saturating additions. @@ -2456,14 +2465,14 @@ mod tests { /// reporting `i32::MAX - 1` for a coordinate whose exact value is `i32::MAX`. #[rstest::rstest] // 1M at MAX: MAX + 1 - 1 == MAX exactly, so no saturation is warranted. - #[case::exactly_max(i32::MAX, "1M", i32::MAX)] + #[case::exactly_max(i32::MAX, b"1M", i32::MAX)] // Genuinely past the top: saturation is correct here. - #[case::past_max(i32::MAX, "10M", i32::MAX)] + #[case::past_max(i32::MAX, b"10M", i32::MAX)] // A malformed MC whose spans saturate must still clamp at the top, not wrap. - #[case::saturated_spans(100, "2000000000M2000000000M", i32::MAX)] + #[case::saturated_spans(100, b"2000000000M2000000000M", i32::MAX)] fn test_unclipped_other_end_saturates_at_the_coordinate_ceiling( #[case] mate_pos: i32, - #[case] mc_cigar: &str, + #[case] mc_cigar: &[u8], #[case] expected: i32, ) { assert_eq!(unclipped_other_end(mate_pos, mc_cigar), expected); @@ -2480,14 +2489,14 @@ mod tests { /// well-meaning `.max(0)`. #[rstest::rstest] // 1-based pos 3 with 10 leading clips: 3 - 10 = -7, before the contig start. - #[case::clips_past_contig_start(3, "10S5M", -7)] - #[case::clips_exactly_to_zero(10, "10S5M", 0)] + #[case::clips_past_contig_start(3, b"10S5M", -7)] + #[case::clips_exactly_to_zero(10, b"10S5M", 0)] // Even a malformed, saturated clip stays a (very negative) number rather // than wrapping: `saturating_sub` bottoms out at `i32::MIN`. - #[case::saturated_clip_does_not_wrap(100, "2000000000S2000000000S", 100 - i32::MAX)] + #[case::saturated_clip_does_not_wrap(100, b"2000000000S2000000000S", 100 - i32::MAX)] fn test_unclipped_other_start_is_not_clamped_to_the_contig_start( #[case] mate_pos: i32, - #[case] mc_cigar: &str, + #[case] mc_cigar: &[u8], #[case] expected: i32, ) { assert_eq!(unclipped_other_start(mate_pos, mc_cigar), expected); diff --git a/crates/fgumi-raw-bam/src/fields.rs b/crates/fgumi-raw-bam/src/fields.rs index 3bf16edcf..e9e0a34da 100644 --- a/crates/fgumi-raw-bam/src/fields.rs +++ b/crates/fgumi-raw-bam/src/fields.rs @@ -279,6 +279,26 @@ pub(crate) const TAG_FIXED_SIZES: [u8; 256] = { table }; +/// Offset of the first NUL byte in `bytes`, or `None` when there is no +/// terminator. +/// +/// Every `Z`- and `H`-typed aux value is NUL-terminated, so this search runs +/// once per such tag and is the inner loop of aux-tag extraction. It is worth a +/// vectorized implementation: on a 1000 Genomes WGS record about 62% of the ~120 +/// aux bytes sit inside `Z` values (`PG`, `MD`, `RG`, `MC`, and `XA` when +/// present), and walking them with a scalar `iter().position()` measured **12.3% +/// of the sort's serial Phase 1 thread** -- the thread that sets 60% of a +/// spill-heavy sort's wall clock. +/// +/// Returning `None` rather than the slice length is load-bearing: callers use it +/// to stop scanning a truncated record instead of treating the remainder as a +/// value. +#[inline] +#[must_use] +pub fn nul_offset(bytes: &[u8]) -> Option { + memchr::memchr(0, bytes) +} + /// Calculate the size of a tag value based on its type. // Inlining always is justified here: `tag_value_size` is the inner dispatch of // the hot `find_tag_position` loop and is tiny (a table lookup + match). @@ -292,7 +312,7 @@ pub fn tag_value_size(val_type: u8, data: &[u8]) -> Option { return Some(fixed as usize); } match val_type { - b'Z' | b'H' => Some(data.iter().position(|&b| b == 0)? + 1), + b'Z' | b'H' => Some(nul_offset(data)? + 1), b'B' => { if data.len() < 5 { return None; @@ -1368,4 +1388,37 @@ mod tests { assert_eq!(v.mate_pos(), 456); assert_eq!(v.template_length(), 300); } + + // ── nul_offset ─────────────────────────────────────────────────────────── + + #[test] + fn test_nul_offset_finds_the_first_terminator() { + assert_eq!(nul_offset(b"RG\0trailing"), Some(2)); + assert_eq!(nul_offset(b"\0"), Some(0)); + // Two terminators: the first one ends the value, the second belongs to + // whatever tag follows it. + assert_eq!(nul_offset(b"ab\0cd\0"), Some(2)); + } + + #[test] + fn test_nul_offset_reports_absence_rather_than_a_length() { + // A truncated aux value has no terminator at all, and the callers rely + // on `None` to stop scanning instead of reading past the record. + assert_eq!(nul_offset(b"unterminated"), None); + assert_eq!(nul_offset(b""), None); + } + + #[test] + fn test_nul_offset_is_exact_past_one_simd_block() { + // The vectorized search steps in blocks, so a terminator that lands just + // after a block boundary is the case a hand-rolled loop and `memchr` + // could disagree on. `XA:Z` values on the measured sample run to ~140 + // bytes, so this is the common length, not an edge case. + for len in [15_usize, 16, 17, 31, 32, 33, 63, 64, 65, 127, 128, 129] { + let mut value = vec![b'M'; len]; + value.push(0); + value.extend_from_slice(b"more"); + assert_eq!(nul_offset(&value), Some(len), "terminator at offset {len}"); + } + } } diff --git a/crates/fgumi-raw-bam/src/lib.rs b/crates/fgumi-raw-bam/src/lib.rs index 236d88b76..a996cb949 100644 --- a/crates/fgumi-raw-bam/src/lib.rs +++ b/crates/fgumi-raw-bam/src/lib.rs @@ -43,6 +43,7 @@ pub use fields::{ mapq, mate_pos, mate_ref_id, + nul_offset, pos, qual_offset, read_name, diff --git a/crates/fgumi-raw-bam/src/tags.rs b/crates/fgumi-raw-bam/src/tags.rs index 3fa6f1a1d..c1e1f7dee 100644 --- a/crates/fgumi-raw-bam/src/tags.rs +++ b/crates/fgumi-raw-bam/src/tags.rs @@ -2,7 +2,7 @@ use fgumi_tag::{AsTagBytes, SamTag}; use crate::fields::{ RawRecordMut, RawRecordView, TAG_FIXED_SIZES, aux_data_offset_from_record, aux_data_slice, - tag_value_size, + nul_offset, tag_value_size, }; /// Find a tag's position and type byte in auxiliary data. @@ -43,7 +43,7 @@ pub fn find_string_tag(aux_data: &[u8], tag: impl AsTagBytes) -> Option<&[u8]> { return None; } let start = p + 3; - let end = aux_data[start..].iter().position(|&b| b == 0)?; + let end = nul_offset(&aux_data[start..])?; Some(&aux_data[start..start + end]) } @@ -68,7 +68,7 @@ pub fn find_string_tag_position(aux_data: &[u8], tag: impl AsTagBytes) -> Option return None; } let start = p + 3; - let len = aux_data[start..].iter().position(|&b| b == 0)?; + let len = nul_offset(&aux_data[start..])?; Some((u32::try_from(start).ok()?, u16::try_from(len).ok()?)) } @@ -206,7 +206,7 @@ pub(crate) fn find_mi_tag(aux_data: &[u8]) -> Option<(u64, bool)> { if val_type == b'Z' { // String type - parse "12345" or "12345/A" or "12345/B" let start = pos + 3; - let end = aux_data[start..].iter().position(|&b| b == 0)?; + let end = nul_offset(&aux_data[start..])?; parse_mi_bytes(&aux_data[start..start + end]) } else { // Integer types: delegate to shared extractor, reject negative values @@ -257,15 +257,23 @@ fn parse_mi_bytes(s: &[u8]) -> Option<(u64, bool)> { /// Find the MC (mate CIGAR) tag in auxiliary data. /// -/// Returns the CIGAR string, or None if not found. +/// Returns the CIGAR bytes, or `None` if not found. +/// +/// No UTF-8 validation: `MC` is a CIGAR that every consumer parses as bytes, +/// and this accessor is the byte-level sibling of the batch extractors +/// ([`extract_aux_string_tags`], [`extract_template_aux_tags`]). Those two +/// hand the raw value over, so gating here would make the two paths disagree +/// on a non-UTF-8 value -- and `src/lib/grouper.rs` cross-checks them against +/// each other, with duplicate `MC` entries documented as the *only* case where +/// they may differ. #[must_use] -pub(crate) fn find_mc_tag(aux_data: &[u8]) -> Option<&str> { - find_string_tag(aux_data, SamTag::MC).and_then(|v| std::str::from_utf8(v).ok()) +pub(crate) fn find_mc_tag(aux_data: &[u8]) -> Option<&[u8]> { + find_string_tag(aux_data, SamTag::MC) } /// Find MC tag in a complete BAM record. #[must_use] -pub fn find_mc_tag_in_record(bam: &[u8]) -> Option<&str> { +pub fn find_mc_tag_in_record(bam: &[u8]) -> Option<&[u8]> { find_mc_tag(aux_data_slice(bam)) } @@ -274,7 +282,7 @@ pub fn find_mc_tag_in_record(bam: &[u8]) -> Option<&str> { pub struct AuxStringTags<'a> { pub rg: Option<&'a [u8]>, pub cell: Option<&'a [u8]>, - pub mc: Option<&'a str>, + pub mc: Option<&'a [u8]>, /// Aux-relative `(value_offset, value_len)` of the UMI tag's value bytes /// (excluding NUL terminator), when the caller supplied a `umi_tag` to look /// for. `None` when no `umi_tag` was supplied, or when the tag is absent @@ -304,7 +312,7 @@ pub fn extract_aux_string_tags( if val_type == b'Z' { let start = p + 3; - if let Some(end) = aux_data[start..].iter().position(|&b| b == 0) { + if let Some(end) = nul_offset(&aux_data[start..]) { let value = &aux_data[start..start + end]; if t == SamTag::RG { result.rg = Some(value); @@ -313,7 +321,7 @@ pub fn extract_aux_string_tags( result.cell = Some(value); found |= 2; } else if t == SamTag::MC { - result.mc = std::str::from_utf8(value).ok(); + result.mc = Some(value); found |= 4; } else if matches!(umi_tag, Some(ut) if t == ut) { // Width checks are defensive only (BAM records < 4 GiB, @@ -352,7 +360,7 @@ pub struct TemplateAuxTags<'a> { /// Cell barcode tag value. pub cell: Option<&'a [u8]>, /// MC (mate CIGAR) tag value. - pub mc: Option<&'a str>, + pub mc: Option<&'a [u8]>, } /// Extract MI, RG, cell barcode, and MC tags in a single pass over aux data. @@ -375,7 +383,7 @@ pub fn extract_template_aux_tags(bam: &[u8], cell_tag: Option) -> Templa if t == SamTag::MI { if val_type == b'Z' { let start = p + 3; - if let Some(end) = aux_data[start..].iter().position(|&b| b == 0) { + if let Some(end) = nul_offset(&aux_data[start..]) { result.mi = parse_mi_bytes(&aux_data[start..start + end]).unwrap_or((0, true)); p = start + end + 1; } else { @@ -407,7 +415,7 @@ pub fn extract_template_aux_tags(bam: &[u8], cell_tag: Option) -> Templa if val_type == b'Z' { let start = p + 3; - if let Some(end) = aux_data[start..].iter().position(|&b| b == 0) { + if let Some(end) = nul_offset(&aux_data[start..]) { let value = &aux_data[start..start + end]; if t == SamTag::RG { result.rg = Some(value); @@ -418,7 +426,7 @@ pub fn extract_template_aux_tags(bam: &[u8], cell_tag: Option) -> Templa found |= 4; } if t == SamTag::MC { - result.mc = std::str::from_utf8(value).ok(); + result.mc = Some(value); found |= 8; } if found & target_bits == target_bits { @@ -893,7 +901,7 @@ fn find_string_tag_range(record: &[u8], aux_offset: usize, tag: [u8; 2]) -> Opti return None; } let start = aux_offset + p + 3; - let nul_off = record[start..].iter().position(|&b| b == 0)?; + let nul_off = nul_offset(&record[start..])?; let end = start + nul_off; if end > start { Some((start, end)) } else { None } } @@ -1172,12 +1180,12 @@ impl<'a> RawTagsView<'a> { find_mi_tag(self.0) } - /// Returns the MC (mate CIGAR) tag as a `&str`. + /// Returns the MC (mate CIGAR) tag as raw CIGAR bytes. /// - /// Returns `None` if the tag is absent or is not valid UTF-8. + /// Returns `None` if the tag is absent or is not `Z`-typed. #[inline] #[must_use] - pub fn find_mc(&self) -> Option<&'a str> { + pub fn find_mc(&self) -> Option<&'a [u8]> { find_mc_tag(self.0) } @@ -1209,11 +1217,11 @@ impl<'a> RawTagsView<'a> { ]))) } b'Z' => { - let end = aux[start..].iter().position(|&b| b == 0)?; + let end = nul_offset(&aux[start..])?; Some(TagValue::String(&aux[start..start + end])) } b'H' => { - let end = aux[start..].iter().position(|&b| b == 0)?; + let end = nul_offset(&aux[start..])?; Some(TagValue::Hex(&aux[start..start + end])) } b'B' => parse_array_tag_at(aux, start).map(TagValue::Array), @@ -1499,7 +1507,7 @@ impl<'a> RawTagsMut<'a> { return false; } let start = p + 3; - let Some(nul_off) = self.0[start..].iter().position(|&b| b == 0) else { + let Some(nul_off) = nul_offset(&self.0[start..]) else { return false; }; if nul_off != value.len() { @@ -2673,7 +2681,7 @@ mod tests { fn test_find_mc_tag_present() { // MC:Z:10M5S\0 let aux = b"MCZ10M5S\x00"; - assert_eq!(find_mc_tag(aux), Some("10M5S")); + assert_eq!(find_mc_tag(aux), Some(b"10M5S".as_slice())); } #[test] @@ -2692,7 +2700,7 @@ mod tests { aux.extend_from_slice(b"NMC"); // tag NM, type C (unsigned byte) aux.push(5); // value aux.extend_from_slice(b"MCZ15M\x00"); // MC:Z:15M - assert_eq!(find_mc_tag(&aux), Some("15M")); + assert_eq!(find_mc_tag(&aux), Some(b"15M".as_slice())); } #[test] @@ -2731,7 +2739,7 @@ mod tests { let result = extract_aux_string_tags(&aux, SamTag::CB, None); assert_eq!(result.rg, Some(b"sample1".as_ref())); assert_eq!(result.cell, Some(b"cell42".as_ref())); - assert_eq!(result.mc, Some("10M5S")); + assert_eq!(result.mc, Some(b"10M5S".as_slice())); } #[test] @@ -2746,7 +2754,7 @@ mod tests { let result = extract_aux_string_tags(&aux, SamTag::CB, None); assert_eq!(result.rg, Some(b"rg1".as_ref())); assert_eq!(result.cell, Some(b"bc1".as_ref())); - assert_eq!(result.mc, Some("5M")); + assert_eq!(result.mc, Some(b"5M".as_slice())); } #[test] @@ -2770,7 +2778,7 @@ mod tests { let result = extract_aux_string_tags(&aux, SamTag::CB, None); assert_eq!(result.rg, Some(b"lib1".as_ref())); assert!(result.cell.is_none()); - assert_eq!(result.mc, Some("20M")); + assert_eq!(result.mc, Some(b"20M".as_slice())); } #[test] @@ -2781,15 +2789,32 @@ mod tests { assert!(result.mc.is_none()); } - #[test] - fn test_extract_aux_string_tags_mc_invalid_utf8() { - // MC tag with invalid UTF-8 bytes → mc should be None + /// A non-UTF-8 `MC` is handed over as bytes rather than discarded. + /// + /// This is a deliberate behavior change. Extraction used to run + /// `str::from_utf8(..).ok()` on the value, so a tag with any non-UTF-8 byte + /// became `None` and the mate position was used unadjusted. `MC` is a CIGAR + /// that every consumer immediately re-borrows as bytes, and validating it + /// per record cost 2.1% of `fgumi sort`'s serial Phase 1 thread, so the gate + /// is gone. + /// + /// The observable difference is narrow: the clip parsers read a leading run + /// of ASCII digits and stop at the first byte they cannot interpret, so a + /// value that is garbage from the start still yields zero clips -- the same + /// answer discarding it gave. Only a value with a *valid CIGAR prefix* + /// followed by invalid bytes now changes, and it changes toward parsing the + /// prefix, which is what samtools does. + #[test] + fn test_extract_aux_string_tags_keeps_a_non_utf8_mc_as_bytes() { let mut aux = Vec::new(); aux.extend_from_slice(b"MCZ"); - aux.extend_from_slice(&[0xFF, 0xFE, 0xFD]); // invalid UTF-8 - aux.push(0); // null terminator + aux.extend_from_slice(&[0xFF, 0xFE, 0xFD]); // not valid UTF-8 + aux.push(0); // NUL terminator let result = extract_aux_string_tags(&aux, SamTag::CB, None); - assert!(result.mc.is_none()); // from_utf8 fails + assert_eq!(result.mc, Some([0xFF, 0xFE, 0xFD].as_slice())); + // Garbage from the first byte still parses to no clips, so this value + // places the mate exactly where discarding the tag used to. + assert_eq!(crate::cigar::mate_unclipped_5prime(100, false, result.mc.expect("mc")), 100); } #[test] @@ -3525,7 +3550,7 @@ mod tests { assert_eq!(result.mi, (42, true)); assert_eq!(result.rg, Some(b"sample1".as_ref())); assert_eq!(result.cell, Some(b"cell99".as_ref())); - assert_eq!(result.mc, Some("10M5S")); + assert_eq!(result.mc, Some(b"10M5S".as_slice())); } #[test] @@ -3540,7 +3565,7 @@ mod tests { assert_eq!(result.mi, (7, false)); assert_eq!(result.rg, Some(b"lib1".as_ref())); assert!(result.cell.is_none()); - assert_eq!(result.mc, Some("20M")); + assert_eq!(result.mc, Some(b"20M".as_slice())); } #[test] @@ -3594,7 +3619,7 @@ mod tests { let result = extract_template_aux_tags(&rec, None); assert_eq!(result.mi, (100, true)); assert_eq!(result.rg, Some(b"lib2".as_ref())); - assert_eq!(result.mc, Some("5M")); + assert_eq!(result.mc, Some(b"5M".as_slice())); } // ======================================================================== @@ -3833,7 +3858,32 @@ mod tests { let s = RawRecordView::new(&rec).tags().extract_string_batch(SamTag::BC, None); assert_eq!(s.rg, Some(b"mygrp".as_slice())); assert_eq!(s.cell, Some(b"ACGT".as_slice())); - assert_eq!(s.mc, Some("50M")); + assert_eq!(s.mc, Some(b"50M".as_slice())); + } + + #[test] + fn test_mc_is_returned_as_bytes_without_a_utf8_gate() { + // `MC` is only ever consumed as CIGAR bytes, so extraction hands the raw + // value over rather than validating it as UTF-8 first. A value with a + // valid CIGAR prefix and a non-UTF-8 byte after it therefore reaches the + // parser, which reads the prefix and stops -- where the earlier + // `str::from_utf8(..).ok()` gate discarded the whole tag and silently + // left the mate position unadjusted. + let aux = b"MCZ10S40M\xff\0RGZrg1\0"; + let rec = make_bam_bytes(0, 0, 0, b"r", &[], 0, -1, -1, aux); + let tags = extract_template_aux_tags(&rec, None); + assert_eq!(tags.mc, Some(b"10S40M\xff".as_slice())); + assert_eq!(tags.rg, Some(b"rg1".as_slice()), "the scan continues past MC"); + } + + #[test] + fn test_mc_bytes_round_trip_through_the_clip_parser() { + // The point of dropping the gate is that the bytes still parse, so pin + // the value the sort key actually uses rather than only the field. + let aux = b"MCZ10S40M\0"; + let rec = make_bam_bytes(0, 0, 0, b"r", &[], 0, -1, -1, aux); + let mc = extract_template_aux_tags(&rec, None).mc.expect("MC present"); + assert_eq!(crate::cigar::mate_unclipped_5prime(100, false, mc), 90); } #[test] diff --git a/crates/fgumi-sam/src/lib.rs b/crates/fgumi-sam/src/lib.rs index 44f7e7ceb..595305616 100644 --- a/crates/fgumi-sam/src/lib.rs +++ b/crates/fgumi-sam/src/lib.rs @@ -29,7 +29,8 @@ //! - [`record_utils::read_pos_at_ref_pos`] - Map reference position to read position //! - [`record_utils::is_fr_pair_from_tags`] - Check if read is part of FR pair using tags //! - [`record_utils::mate_unclipped_start`] / [`record_utils::mate_unclipped_end`] - Get mate boundaries from MC tag -//! - [`record_utils::parse_cigar_string`] - Parse CIGAR string to operations +//! - [`record_utils::parse_cigar_string`] / [`record_utils::parse_cigar_bytes`] - +//! Parse a CIGAR to operations, from `&str` or from raw bytes pub mod alignment_tags; pub mod builder; @@ -116,9 +117,9 @@ pub use fgumi_tag::SamTag; pub use record_utils::{ PairOrientation, alignment_end, cigar_reference_length, get_pair_orientation, is_fr_pair, is_fr_pair_from_tags, leading_clipping, leading_soft_clipping, mate_unclipped_end, - mate_unclipped_start, parse_cigar_string, read_pos_at_ref_pos, reference_length, - trailing_clipping, trailing_soft_clipping, unclipped_end, unclipped_five_prime_position, - unclipped_start, + mate_unclipped_start, parse_cigar_bytes, parse_cigar_string, read_pos_at_ref_pos, + reference_length, trailing_clipping, trailing_soft_clipping, unclipped_end, + unclipped_five_prime_position, unclipped_start, }; pub use template_coordinate::{TC_TAG, TemplateCoordinateInfo}; diff --git a/crates/fgumi-sam/src/record_utils.rs b/crates/fgumi-sam/src/record_utils.rs index 15a0cabb7..48ed303fe 100644 --- a/crates/fgumi-sam/src/record_utils.rs +++ b/crates/fgumi-sam/src/record_utils.rs @@ -167,13 +167,17 @@ fn parse_mate_cigar(read: &RecordBuf) -> Option<(Vec<(Kind, usize)>, usize)> { let mc_tag = Tag::from([b'M', b'C']); let mc_value = read.data().get(&mc_tag)?; - let cigar_str = match mc_value { - Value::String(s) => std::str::from_utf8(s.as_ref()).ok()?, + let cigar_bytes = match mc_value { + // Parse the raw bytes directly rather than gating on UTF-8. A valid + // CIGAR prefix followed by non-UTF-8 trailing bytes is still usable, and + // this keeps the typed path in parity with the raw-byte sibling + // (`mate_unclipped_*_raw`), which parses the same bytes. + Value::String(s) => s.as_ref(), _ => return None, }; let mate_start = usize::from(read.mate_alignment_start()?); - let ops = parse_cigar_string(cigar_str); + let ops = parse_cigar_bytes(cigar_bytes); if ops.is_empty() { return None; } @@ -242,26 +246,38 @@ pub fn alignment_end(read: &RecordBuf) -> Option { /// This is useful for parsing the MC (mate CIGAR) tag value. #[must_use] pub fn parse_cigar_string(cigar_str: &str) -> Vec<(Kind, usize)> { + parse_cigar_bytes(cigar_str.as_bytes()) +} + +/// Parses a CIGAR from raw bytes, as [`parse_cigar_string`] does from `&str`. +/// +/// `MC` is extracted as bytes (`find_mc_tag_in_record`) so that the sort's +/// per-record ingest path does not pay a UTF-8 validation it cannot use, and +/// this is the entry point that consumes it. Operators are ASCII by +/// construction, so a non-ASCII byte falls into the same "unknown operator" +/// arm a non-ASCII `char` did: the two functions agree on every input. +#[must_use] +pub fn parse_cigar_bytes(cigar: &[u8]) -> Vec<(Kind, usize)> { let mut ops = Vec::new(); let mut num_str = String::new(); - for ch in cigar_str.chars() { - if ch.is_ascii_digit() { - num_str.push(ch); + for &b in cigar { + if b.is_ascii_digit() { + num_str.push(char::from(b)); } else { let len: usize = num_str.parse().unwrap_or(0); num_str.clear(); - let kind = match ch { - 'M' => Kind::Match, - 'I' => Kind::Insertion, - 'D' => Kind::Deletion, - 'N' => Kind::Skip, - 'S' => Kind::SoftClip, - 'H' => Kind::HardClip, - 'P' => Kind::Pad, - '=' => Kind::SequenceMatch, - 'X' => Kind::SequenceMismatch, + let kind = match b { + b'M' => Kind::Match, + b'I' => Kind::Insertion, + b'D' => Kind::Deletion, + b'N' => Kind::Skip, + b'S' => Kind::SoftClip, + b'H' => Kind::HardClip, + b'P' => Kind::Pad, + b'=' => Kind::SequenceMatch, + b'X' => Kind::SequenceMismatch, _ => continue, }; @@ -622,7 +638,7 @@ pub fn mate_unclipped_start_raw(bam: &[u8]) -> Option { if mate_start_0based < 0 { return None; } - let ops = parse_cigar_string(mc_cigar); + let ops = parse_cigar_bytes(mc_cigar); if ops.is_empty() { return None; } @@ -655,7 +671,7 @@ pub fn mate_unclipped_end_raw(bam: &[u8]) -> Option { reason = "mate_start_0based is verified non-negative by the guard above" )] let mate_start_1based = (mate_start_0based + 1) as usize; - let ops = parse_cigar_string(mc_cigar); + let ops = parse_cigar_bytes(mc_cigar); if ops.is_empty() { return None; } @@ -769,6 +785,7 @@ mod tests { use super::*; use crate::builder::{RecordBuilder, RecordPairBuilder}; use noodles::sam::alignment::record::Flags; + use rstest::rstest; // Flag constants for test readability const FLAG_PAIRED: u16 = 0x1; @@ -1190,6 +1207,51 @@ mod tests { assert_eq!(mate_unclipped_end(&record), None); } + /// A valid CIGAR prefix followed by non-UTF-8 trailing bytes must yield the + /// same mate coordinates on the typed `RecordBuf` path as on the raw-byte + /// sibling (which parses the bytes directly). The byte parser stops at the + /// prefix; the typed path used to reject the whole value on a UTF-8 gate, + /// disagreeing with the raw path on identical input. + #[rstest] + #[case::leading_soft_clip(b"10S50M", 190, 249)] + #[case::trailing_soft_clip(b"50M10S", 200, 259)] + fn test_mate_unclipped_valid_prefix_non_utf8_mc_matches_raw( + #[case] mc_prefix: &[u8], + #[case] expected_start: isize, + #[case] expected_end: usize, + ) { + let mc_tag = Tag::from([b'M', b'C']); + // The BAM encoder rejects non-UTF-8 strings, so encode with a valid + // placeholder trailing byte ('x'), then overwrite it with 0xFF in both + // the raw buffer and the typed record. 'x' and 0xFF are both non-CIGAR + // bytes the byte parser skips, so the coordinate is identical either + // way — the point is a genuinely non-UTF-8 MC value on both paths. + let mut placeholder = mc_prefix.to_vec(); + placeholder.push(b'x'); + + let mut record = create_mc_test_read("valid_prefix_non_utf8", 200, "50M"); + record.data_mut().insert(mc_tag, Value::String(placeholder.clone().into())); + let mut raw = to_raw(&record); + + // Inject the non-UTF-8 byte into the raw record's MC value. + let pos = raw + .windows(placeholder.len()) + .position(|w| w == placeholder.as_slice()) + .expect("placeholder MC value present in raw record"); + raw[pos + placeholder.len() - 1] = 0xFF; + + // And into the typed record. + let mut non_utf8 = mc_prefix.to_vec(); + non_utf8.push(0xFF); + record.data_mut().insert(mc_tag, Value::String(non_utf8.into())); + + assert_eq!(mate_unclipped_start(&record), Some(expected_start)); + assert_eq!(mate_unclipped_end(&record), Some(expected_end)); + // Parity with the raw-byte sibling on the identical non-UTF-8 input. + assert_eq!(mate_unclipped_start_raw(&raw), mate_unclipped_start(&record)); + assert_eq!(mate_unclipped_end_raw(&raw), mate_unclipped_end(&record)); + } + #[test] fn test_mate_unclipped_empty_cigar_mc_tag() { let record = create_mc_test_read("empty_cigar", 200, ""); diff --git a/crates/fgumi-sort/Cargo.toml b/crates/fgumi-sort/Cargo.toml index 54a505302..018280508 100644 --- a/crates/fgumi-sort/Cargo.toml +++ b/crates/fgumi-sort/Cargo.toml @@ -60,3 +60,7 @@ rstest = { workspace = true } proptest = { workspace = true } criterion = { workspace = true } tempfile = { workspace = true } + +[[bench]] +name = "template_key" +harness = false diff --git a/crates/fgumi-sort/benches/template_key.rs b/crates/fgumi-sort/benches/template_key.rs new file mode 100644 index 000000000..817ca90b1 --- /dev/null +++ b/crates/fgumi-sort/benches/template_key.rs @@ -0,0 +1,134 @@ +//! Per-record cost of the template-coordinate ingest path. +//! +//! Phase 1 of a spill-heavy sort is bound by its serial main thread: on +//! `1kg-wgs-HG00096` at 16 threads that thread is busy 219.5s of a 240.7s +//! phase, and `perf` puts 39.5% of it in `extract_template_key_inline` with a +//! further 12.3% in the slice-iterator `next()` the aux-tag scan walks bytes +//! with. This bench sizes that function directly, so a change to it is measured +//! before a seven-minute whole-genome run is spent on it. +//! +//! The records are built to the aux layout the measured sample actually carries +//! (`PG:Z AS:i XS:i MD:Z NM:i RG:Z MQ:i MC:Z`, 118 bytes) rather than a minimal +//! one: the scan's cost is per aux byte, so a record with two tags would report +//! a number that has nothing to do with the workload. `with_xa` adds the long +//! `XA:Z` alt-hit tag that a fraction of records carry, since it roughly doubles +//! the aux data and is where the scan's tail lives. + +use bstr::BString; +use criterion::{Criterion, criterion_group, criterion_main}; +use fgumi_raw_bam::SamTag; +use noodles::sam::Header; +use noodles::sam::header::record::value::Map; +use noodles::sam::header::record::value::map::ReadGroup; +use noodles::sam::header::record::value::map::read_group::tag as rg_tag; +use std::hint::black_box; + +const READ_NAME: &[u8] = b"A00132:53:HFHJKDSXX:1:1646:26467:33332\0"; +const RG_ID: &str = "HG00096_CGGACAAC-TCCGGATT_HFHJKDSXX_L001"; +const SEQ_LEN: usize = 151; + +/// Append a `Z`-typed aux tag (two tag bytes, `Z`, value, NUL). +fn push_z(aux: &mut Vec, tag: SamTag, value: &[u8]) { + aux.extend_from_slice(&*tag); + aux.push(b'Z'); + aux.extend_from_slice(value); + aux.push(0); +} + +/// Append an `i`-typed (32-bit signed) aux tag. +fn push_i(aux: &mut Vec, tag: SamTag, value: i32) { + aux.extend_from_slice(&*tag); + aux.push(b'i'); + aux.extend_from_slice(&value.to_le_bytes()); +} + +/// One BAM record body (from `ref_id`, i.e. without the `block_size` prefix), +/// carrying the aux layout of the measured 1000 Genomes WGS sample. +fn record(pos: i32, mate_pos: i32, with_xa: bool) -> Vec { + let mut aux = Vec::with_capacity(256); + push_z(&mut aux, SamTag::PG, b"MarkDuplicates"); + push_i(&mut aux, SamTag::AS, 64); + push_i(&mut aux, SamTag::XS, 61); + push_z(&mut aux, SamTag::MD, b"0N0N0N0N0N2A61"); + push_i(&mut aux, SamTag::NM, 6); + push_z(&mut aux, SamTag::RG, RG_ID.as_bytes()); + push_i(&mut aux, SamTag::MQ, 0); + push_z(&mut aux, SamTag::MC, b"81S69M"); + if with_xa { + push_z( + &mut aux, + SamTag::new(b'X', b'A'), + b"chr3,+198173832,34M116S,0;chr12,-108091,117S33M,0;chr1,-180805,118S32M,0;\ + chr1,-10052,118S32M,0;chr3,-10519,118S32M,0;", + ); + } + + let n_cigar_op: u16 = 2; + let mut rec = Vec::with_capacity(32 + READ_NAME.len() + 8 + SEQ_LEN * 2 + aux.len()); + rec.extend_from_slice(&0i32.to_le_bytes()); // ref_id + rec.extend_from_slice(&pos.to_le_bytes()); // pos + rec.push(u8::try_from(READ_NAME.len()).expect("read name fits a u8")); + rec.push(60); // mapq + rec.extend_from_slice(&0u16.to_le_bytes()); // bin + rec.extend_from_slice(&n_cigar_op.to_le_bytes()); + // PAIRED | PROPER_PAIR | REVERSE | FIRST_SEGMENT. The mate is forward, so + // the mate lane resolves through `unclipped_other_start` (leading clips). + rec.extend_from_slice(&0x0053u16.to_le_bytes()); + rec.extend_from_slice(&i32::try_from(SEQ_LEN).expect("seq len fits").to_le_bytes()); + rec.extend_from_slice(&0i32.to_le_bytes()); // next_ref_id + rec.extend_from_slice(&mate_pos.to_le_bytes()); + rec.extend_from_slice(&0i32.to_le_bytes()); // tlen + rec.extend_from_slice(READ_NAME); + // CIGAR 81S69M: (len << 4) | op, S = 4, M = 0. + rec.extend_from_slice(&((81u32 << 4) | 4).to_le_bytes()); + rec.extend_from_slice(&(69u32 << 4).to_le_bytes()); + rec.resize(rec.len() + SEQ_LEN.div_ceil(2), 0x11); // packed seq + rec.resize(rec.len() + SEQ_LEN, 30); // qual + rec.extend_from_slice(&aux); + rec +} + +/// A header declaring the sample's single read group, so the RG lookup resolves +/// rather than missing (a miss and a hit take different paths through the map). +fn header() -> Header { + let read_group = Map::::builder() + .insert(rg_tag::LIBRARY, String::from("lib1")) + .build() + .expect("a read group with only LB is valid"); + Header::builder().add_read_group(BString::from(RG_ID), read_group).build() +} + +fn bench_template_key(c: &mut Criterion) { + let lookup = fgumi_sort::LibraryLookup::from_header(&header()); + let hasher = fgumi_sort::cb_hasher(); + + for (label, with_xa) in [("aux-122b", false), ("aux-with-xa", true)] { + // A batch rather than one record: the scan is memory-bound over the aux + // bytes, and re-reading one cache-hot record would flatter it. + let records: Vec> = + (0..1024).map(|i| record(10_000 + i, 133_000_000 + i, with_xa)).collect(); + let aux_bytes: usize = + records.iter().map(|r| fgumi_raw_bam::aux_data_slice(r).len()).sum::() + / records.len(); + + let mut group = c.benchmark_group("extract_template_key_inline"); + group.throughput(criterion::Throughput::Elements(records.len() as u64)); + group.bench_function(format!("{label}-mean-aux-{aux_bytes}"), |b| { + b.iter(|| { + for rec in &records { + let key = fgumi_sort::extract_template_key_inline( + black_box(rec.as_slice()), + &lookup, + None, + &hasher, + ); + black_box(&key); + } + }); + }); + group.finish(); + } +} + +criterion_group!(benches, bench_template_key); +criterion_main!(benches); diff --git a/crates/fgumi-sort/src/external.rs b/crates/fgumi-sort/src/external.rs index 18ad969e3..645852f43 100644 --- a/crates/fgumi-sort/src/external.rs +++ b/crates/fgumi-sort/src/external.rs @@ -245,7 +245,13 @@ impl SortPhaseTimer { /// `max_temp_files` is reported so a run that consolidated says which limit /// it consolidated against. #[allow(clippy::cast_precision_loss)] - fn log_summary(&self, sort_threads: usize, merge_threads: usize, max_temp_files: usize) { + fn log_summary( + &self, + sort_threads: usize, + merge_threads: usize, + max_temp_files: usize, + phase1: Phase1FloorInputs, + ) { let overall = self.overall_start.map_or(0.0, |s| s.elapsed().as_secs_f64()); // Guard against division by zero when sort completes in negligible time. let overall_nonzero = if overall > 0.0 { overall } else { f64::EPSILON }; @@ -289,14 +295,167 @@ impl SortPhaseTimer { } stat!(" Total wall clock: {overall:.1}s"); stat!(" Threads: {}", format_thread_counts(sort_threads, merge_threads)); + self.log_phase1_floor(phase1); stat!("========================="); } + + /// Which of three limits Phase 1's *ingest* is against, and what is + /// recoverable without doing less work. + /// + /// Scoped to the read span rather than the whole phase, deliberately. The + /// in-memory sort is parallel (rayon) and the spill write-out overlaps the + /// next read, so folding them in would put parallel work on the same side of + /// the comparison as one thread's serial CPU and report the difference as + /// "coordination" -- naming a limit that is not there. The read span is the + /// part where one thread reads every record while the pool feeds it, which is + /// exactly the shape [`crate::merge_headroom`] models. + /// + /// Externally sampled, this phase is 60% of a whole-genome sort's wall clock + /// with its main thread 91% busy while 16 cores average 5.3. If that holds + /// in-process, the binding limit is the ingest thread and no amount of + /// additional worker capacity moves it. + fn log_phase1_floor(&self, phase1: Phase1FloorInputs) { + if self.read_secs <= 0.0 { + return; + } + let ingest = phase1.ingest; + let floors = crate::merge_headroom::MergeFloors { + loop_secs: self.read_secs, + consumer_secs: (self.read_secs - ingest.park_secs).max(0.0), + worker_busy_secs: phase1.input_busy_secs, + threads: phase1.threads.max(1), + }; + stat!(" Phase 1 ingest floor: {} is the limit", floors.binding().label()); + stat!( + " ingest serial {:.1}s | worker capacity {:.1}s ({} threads) | read span {:.1}s", + floors.consumer_secs, + floors.worker_floor_secs(), + floors.threads, + self.read_secs + ); + stat!( + " recoverable without doing less work: {:.1}s ({:.0}% of the read span)", + floors.recoverable_secs(), + 100.0 * floors.recoverable_share() + ); + if ingest.parks > 0 { + // Mean park separates a supply problem from a handoff problem: many + // short parks and few long ones need opposite fixes, and the totals + // alone cannot tell them apart. + stat!( + " ingest parked {:.1}s over {} parks ({:.0} us each): {} starved, {} head-of-line", + ingest.park_secs, + ingest.parks, + ingest.mean_park_micros().unwrap_or(0.0), + ingest.parks_starved, + ingest.parks_head_of_line + ); + } else { + stat!(" ingest never parked: the pool always had the next block ready"); + } + if ingest.spill_waits > 0 { + stat!( + " waited {:.1}s over {} spill handoffs (outside every phase bucket)", + ingest.spill_wait_secs, + ingest.spill_waits + ); + } + self.log_ingest_partition(phase1); + } + + /// What the ingest thread's serial CPU is made of, per segment. + /// + /// Correction happens at the sampled scale and scaling second: subtracting + /// one clock pair from the *scaled* total instead would understate the + /// correction by the scale factor, which on a 1-in-1021 sample is three + /// orders of magnitude. + #[allow(clippy::cast_precision_loss, reason = "record counts stay below 2^52")] + fn log_ingest_partition(&self, phase1: Phase1FloorInputs) { + if phase1.samples == 0 || phase1.records == 0 { + return; + } + let scale = phase1.records as f64 / phase1.samples as f64; + let segments = + phase1.sample.corrected(phase1.samples, phase1.clock_overhead_nanos).scaled(scale); + let partition = crate::phase1_stats::IngestPartition { + segments, + read_secs: self.read_secs, + park_secs: phase1.ingest.park_secs, + }; + let per_record = |secs: f64| 1e9 * secs / phase1.records as f64; + stat!( + " Ingest segments ({} samples of {} records, scaled {scale:.0}x, clock {}ns/pair)", + phase1.samples, + phase1.records, + phase1.clock_overhead_nanos + ); + stat!( + " fetch next record: {:.1}s ({:.0} ns/rec) [includes {:.1}s parked]", + segments.fetch, + per_record(segments.fetch), + partition.park_secs + ); + stat!( + " extract key: {:.1}s ({:.0} ns/rec)", + segments.key, + per_record(segments.key) + ); + stat!( + " verify lanes: {:.1}s ({:.0} ns/rec)", + segments.verify, + per_record(segments.verify) + ); + stat!( + " push to arena: {:.1}s ({:.0} ns/rec)", + segments.push, + per_record(segments.push) + ); + stat!( + " progress tick: {:.1}s ({:.0} ns/rec)", + segments.tick, + per_record(segments.tick) + ); + stat!( + " probe + mem check: {:.1}s ({:.0} ns/rec)", + segments.probe, + per_record(segments.probe) + ); + // Signed: a negative residual means the segments over-attribute, which is + // the failure mode this partition exists to make visible. + stat!( + " unattributed: {:+.1}s ({:+.0}% of the read span)", + partition.residual_secs(), + 100.0 * partition.residual_share() + ); + } } // ============================================================================ // Library Lookup for Template-Coordinate Sort // ============================================================================ +/// What the Phase 1 floor line needs that [`SortPhaseTimer`] cannot see: the +/// pool's worker time and the ingest thread's waits. +#[derive(Debug, Clone, Copy, Default)] +pub(crate) struct Phase1FloorInputs { + /// Worker seconds reading and decompressing input blocks. + pub(crate) input_busy_secs: f64, + /// Active Phase 1 worker threads the busy total is spread over. + pub(crate) threads: usize, + /// What the ingest thread waited for. + pub(crate) ingest: crate::phase1_stats::Phase1IngestReport, + /// Raw (unscaled, uncorrected) sub-phase sample from the ingest loop, when + /// that loop is instrumented. Left at zero by the orders that are not. + pub(crate) sample: crate::phase1_stats::IngestSample, + /// Records timed for `sample`. + pub(crate) samples: u64, + /// Records the loop processed, the numerator of the sampling scale. + pub(crate) records: u64, + /// Measured cost of one `Instant::now()`/`elapsed()` pair, subtracted once + /// per segment per sample. + pub(crate) clock_overhead_nanos: u64, +} + /// Deterministic hasher for cell barcode hashing in template-coordinate sort. /// /// Uses arbitrary fixed seeds so that hash values are reproducible across runs. @@ -504,13 +663,32 @@ fn process_umask() -> u32 { u32::from(previous) } +/// RG-id to library-ordinal map, hashed with `ahash` rather than std's `SipHash`. +/// +/// `ordinal_from_rg` runs once per record on the sort's serial Phase 1 thread -- +/// the thread that sets 60% of a spill-heavy sort's wall clock -- and hashing a +/// ~40-byte RG id with `SipHash` measured **6.1% of that thread** on +/// `1kg-wgs-HG00096` (`core::hash::sip::Hasher::write`, third-largest entry in +/// its profile). The keys come from the header of a BAM the caller chose to +/// sort, so there is no adversarial-input exposure that would argue for +/// `SipHash`'s collision resistance here. +/// +/// Unlike [`cb_hasher`] and [`LibraryLookup::hasher`], this one is left at +/// `ahash`'s default randomly-seeded state rather than `with_seeds`. Those two +/// feed hash values *into the sort key*, so they must be identical across +/// processes or the same input would sort differently run to run. This map's +/// hasher is never observed: it is probed only by `get`, and +/// `distinct_header_ordinals` collects its values into a set, so neither the +/// hash values nor the iteration order reaches an output. +type RgOrdinalMap = HashMap, u32, ahash::RandomState>; + /// Maps read group ID -> library ordinal for O(1) comparison. /// /// Pre-computes ordinals by sorting library names alphabetically. /// Empty/unknown library sorts first (ordinal 0). pub struct LibraryLookup { /// RG ID -> library ordinal - rg_to_ordinal: HashMap, u32>, + rg_to_ordinal: RgOrdinalMap, /// Deterministic hasher for read name hashing, constructed once for reuse. hasher: ahash::RandomState, } @@ -542,7 +720,7 @@ impl LibraryLookup { } // Build RG ID -> ordinal mapping - let rg_to_ordinal: HashMap, u32> = header + let rg_to_ordinal: RgOrdinalMap = header .read_groups() .iter() .map(|(id, rg)| { @@ -1242,10 +1420,42 @@ struct MergeConsumerDiag { reassembled: u64, } +/// Per-merge tallies of the record-fetch fast/slow split. +/// +/// [`RECORD_BORROWED`] and [`RECORD_REASSEMBLED`] used to be incremented per +/// record, directly from the merge consumer. That is the one serial thread that +/// touches every record of the merge, and it is the merge's binding floor on the +/// measured cell (`consumer serial 112.9s` against `worker capacity 89.4s` and a +/// 157.1s loop). A relaxed `fetch_add` there is not free: on aarch64 it is an +/// outline-atomics call into `__aarch64_ldadd8_relax`, the same helper that +/// measured 28% of this thread's cycles while the progress counter used it +/// per-record (see [`crate::progress_batch`]). +/// +/// Counting into locals and publishing once per merge keeps every number the +/// reports read -- they consume the statics as a before/after delta around the +/// loop -- and removes two atomics per record from the critical path. Publish +/// *before* that delta is read, or the merge reports zero. +#[derive(Default)] +struct RecordFetchCounts { + /// Records handed over borrowed from the current decompressed block. + borrowed: u64, + /// Records reassembled into scratch because they straddled a block boundary. + reassembled: u64, +} + +impl RecordFetchCounts { + /// Fold this merge's tallies into the process-wide totals. + fn publish(&self) { + RECORD_BORROWED.fetch_add(self.borrowed, std::sync::atomic::Ordering::Relaxed); + RECORD_REASSEMBLED.fetch_add(self.reassembled, std::sync::atomic::Ordering::Relaxed); + } +} + #[inline] fn winner_record_bytes<'a, K: RawSortKey + Default + 'static>( source: &'a ChunkSource, consumer: Option<&'a MainThreadChunkConsumer>, + counts: &mut RecordFetchCounts, ) -> Result<&'a [u8]> { match source { ChunkSource::PoolDisk { source_id, scratch } => { @@ -1262,12 +1472,14 @@ fn winner_record_bytes<'a, K: RawSortKey + Default + 'static>( // that touches every record, so a clock read here would cost more // than the step. What was unknown is how often the copy path fires // at all -- a frequency answers that, and the per-record cost is a - // memcpy of a ~100-byte record either way. + // memcpy of a ~100-byte record either way. The tally is a local + // (see [`RecordFetchCounts`]); an atomic per record on this thread + // is itself measurable. if let Some(borrowed) = consumer.current_record_bytes(*source_id) { - RECORD_BORROWED.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + counts.borrowed += 1; Ok(borrowed) } else { - RECORD_REASSEMBLED.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + counts.reassembled += 1; Ok(scratch.as_slice()) } } @@ -2707,7 +2919,16 @@ impl RawExternalSorter { pool: &std::sync::Arc, ) -> Result<()> { if let Some(prev) = pending.take() { + // Timed because it lands in no phase bucket: this sits between + // `end_read_span` and `time_sort`, so without its own counter the + // ingest thread's wall clock and its own CPU cannot be reconciled and + // the difference shows up only as an unexplained residual against + // total wall clock. + let waited_at = Instant::now(); prev.handle.wait()?; + pool.phase1_ingest_stats().record_spill_wait( + u64::try_from(waited_at.elapsed().as_nanos()).unwrap_or(u64::MAX), + ); timer.record_spill_growth(&prev.chunk_path, prev.size_before); // A chunk that extended an existing run is already represented in // `chunk_files`; only a chunk that started a run adds a merge source. @@ -3296,6 +3517,13 @@ impl RawExternalSorter { let rayon_pool = self.build_sort_rayon_pool()?; let progress = ProgressTracker::new("Read records").with_interval(1_000_000); + // Count records locally and forward in batches: `log_if_needed` does a + // relaxed `fetch_add` per call, which on aarch64 is an outline-atomics + // call into `__aarch64_ldadd8_relax`. On the ingest thread -- the serial + // thread that sets 60% of a spill-heavy sort's wall clock -- that helper + // measured 3.4% of the profile. The merge loops already batch for the + // same reason. + let mut progress_batch = crate::progress_batch::BatchedProgress::new(); debug!("Phase 1: Reading and sorting chunks (inline buffer, keyed output)..."); let mut probe = SpillProbe::new("phase1"); @@ -3305,7 +3533,7 @@ impl RawExternalSorter { // call, which is fine — `push_coordinate` copies the bytes into the buffer). while let Some(record) = record_source.next_record_borrowed()? { stats.total_records += 1; - progress.log_if_needed(1); + progress_batch.tick(&progress); // Push directly to buffer - key extracted inline from raw bytes buffer.push_coordinate(record)?; @@ -3366,6 +3594,7 @@ impl RawExternalSorter { } timer.end_read_span(); + progress_batch.flush(&progress); progress.log_final(); if let Some(err) = record_source.take_error() { return Err(anyhow::Error::from(err)); @@ -3383,6 +3612,15 @@ impl RawExternalSorter { probe.phase1_end(buffer.memory_usage() as u64); // Ingest is done: everything from here is Phase 2 (merge/write). + // Snapshot Phase 1's floor inputs before the pool is handed to the merge: + // the counters describe the phase that has just ended, and reading them + // later would both borrow a moved value and risk folding in Phase 2 work. + let phase1_floor = Phase1FloorInputs { + input_busy_secs: pool.phase1_input_busy_secs(), + threads: self.phase1_threads(), + ingest: pool.phase1_ingest_stats().snapshot(), + ..Phase1FloorInputs::default() + }; self.enter_output_phase(&pool); if chunk_files.is_empty() { @@ -3447,7 +3685,12 @@ impl RawExternalSorter { if let Ok(pool) = Arc::try_unwrap(pool) { pool.shutdown(); } - timer.log_summary(self.phase1_threads(), self.phase2_threads(), self.max_temp_files); + timer.log_summary( + self.phase1_threads(), + self.phase2_threads(), + self.max_temp_files, + phase1_floor, + ); debug!("Sort complete: {} records processed", stats.total_records); Ok(stats) @@ -3495,9 +3738,17 @@ impl RawExternalSorter { debug!("Phase 1: Reading and sorting chunks (inline buffer, keyed output)..."); let mut probe = SpillProbe::new("phase1"); - for record in record_source.by_ref() { + // Borrow each record's bytes straight out of the decompressed block, as the + // non-indexed sibling does. `RecordSource`'s `Iterator` impl yields an owned + // `RawRecord`, which is a heap allocation plus a full-record memcpy per + // record -- freed again as soon as `push_coordinate` has copied the bytes + // into the arena -- on the serial thread that sets 60% of a spill-heavy + // sort's wall clock. Nothing here needs the record to outlive the push: the + // BAI is built by the writer from BGZF offsets during the merge, not from + // this loop. + while let Some(record) = record_source.next_record_borrowed()? { stats.total_records += 1; - buffer.push_coordinate(record.as_ref())?; + buffer.push_coordinate(record)?; if probe.should_sample_read(stats.total_records) { probe.log_mid_read(probe_stats(&buffer), Some(pool.phase1_queue_depths())); @@ -3572,6 +3823,15 @@ impl RawExternalSorter { let output_header = self.create_output_header(header); // Ingest is done: everything from here is Phase 2 (merge/write). + // Snapshot Phase 1's floor inputs before the pool is handed to the merge: + // the counters describe the phase that has just ended, and reading them + // later would both borrow a moved value and risk folding in Phase 2 work. + let phase1_floor = Phase1FloorInputs { + input_busy_secs: pool.phase1_input_busy_secs(), + threads: self.phase1_threads(), + ingest: pool.phase1_ingest_stats().snapshot(), + ..Phase1FloorInputs::default() + }; self.enter_output_phase(&pool); if chunk_files.is_empty() { @@ -3646,7 +3906,12 @@ impl RawExternalSorter { if let Ok(pool) = Arc::try_unwrap(pool) { pool.shutdown(); } - timer.log_summary(self.phase1_threads(), self.phase2_threads(), self.max_temp_files); + timer.log_summary( + self.phase1_threads(), + self.phase2_threads(), + self.max_temp_files, + phase1_floor, + ); debug!("Sort complete: {} records processed", stats.total_records); Ok(stats) @@ -3718,12 +3983,19 @@ impl RawExternalSorter { let rayon_pool = self.build_sort_rayon_pool()?; let progress = ProgressTracker::new("Read records").with_interval(1_000_000); + // Count records locally and forward in batches: `log_if_needed` does a + // relaxed `fetch_add` per call, which on aarch64 is an outline-atomics + // call into `__aarch64_ldadd8_relax`. On the ingest thread -- the serial + // thread that sets 60% of a spill-heavy sort's wall clock -- that helper + // measured 3.4% of the profile. The merge loops already batch for the + // same reason. + let mut progress_batch = crate::progress_batch::BatchedProgress::new(); debug!("Phase 1: Reading and sorting chunks (keyed output)..."); let mut probe = SpillProbe::new("phase1"); for record in record_source.by_ref() { stats.total_records += 1; - progress.log_if_needed(1); + progress_batch.tick(&progress); // Extract key from raw bytes. Stamp the ingest position within this // chunk so the key is totally ordered: read name + flags alone is not @@ -3800,6 +4072,7 @@ impl RawExternalSorter { } timer.end_read_span(); + progress_batch.flush(&progress); progress.log_final(); if let Some(err) = record_source.take_error() { return Err(anyhow::Error::from(err)); @@ -3817,6 +4090,15 @@ impl RawExternalSorter { probe.phase1_end(memory_used as u64); // Ingest is done: everything from here is Phase 2 (merge/write). + // Snapshot Phase 1's floor inputs before the pool is handed to the merge: + // the counters describe the phase that has just ended, and reading them + // later would both borrow a moved value and risk folding in Phase 2 work. + let phase1_floor = Phase1FloorInputs { + input_busy_secs: pool.phase1_input_busy_secs(), + threads: self.phase1_threads(), + ingest: pool.phase1_ingest_stats().snapshot(), + ..Phase1FloorInputs::default() + }; self.enter_output_phase(&pool); if chunk_files.is_empty() { @@ -3921,7 +4203,12 @@ impl RawExternalSorter { if let Ok(pool) = Arc::try_unwrap(pool) { pool.shutdown(); } - timer.log_summary(self.phase1_threads(), self.phase2_threads(), self.max_temp_files); + timer.log_summary( + self.phase1_threads(), + self.phase2_threads(), + self.max_temp_files, + phase1_floor, + ); debug!("Sort complete: {} records processed", stats.total_records); Ok(stats) @@ -4080,6 +4367,13 @@ impl RawExternalSorter { let rayon_pool = self.build_sort_rayon_pool()?; let progress = ProgressTracker::new("Read records").with_interval(1_000_000); + // Count records locally and forward in batches: `log_if_needed` does a + // relaxed `fetch_add` per call, which on aarch64 is an outline-atomics + // call into `__aarch64_ldadd8_relax`. On the ingest thread -- the serial + // thread that sets 60% of a spill-heavy sort's wall clock -- that helper + // measured 3.4% of the profile. The merge loops already batch for the + // same reason. + let mut progress_batch = crate::progress_batch::BatchedProgress::new(); debug!("Phase 1: Reading and sorting chunks (inline buffer)..."); let mut probe = SpillProbe::new("phase1"); @@ -4089,7 +4383,7 @@ impl RawExternalSorter { // exceed the memory limit, so no spill check is needed here. if let Some(record) = first_record { stats.total_records += 1; - progress.log_if_needed(1); + progress_batch.tick(&progress); let bam_bytes = record.as_ref(); let full = extract_template_key_inline(bam_bytes, lib_lookup, self.cell_tag, cb_hasher); @@ -4103,32 +4397,84 @@ impl RawExternalSorter { buffer.push(bam_bytes, K::from_full(&full))?; } + // Sub-phase timing for the ingest thread's serial CPU, which the floor + // line identifies as this phase's binding limit. Sampled 1-in-N and + // clock-corrected, per `crate::phase1_stats::IngestSample`; the segments + // are checked against the measured read span with a signed residual, so a + // partition that over-attributes says so instead of looking tidy. + let ingest_sample_interval = crate::phase1_stats::INGEST_SAMPLE_INTERVAL; + let clock_overhead_nanos = crate::merge_headroom::measure_clock_overhead_nanos(); + let mut ingest_raw = crate::phase1_stats::IngestSample::default(); + let mut ingest_samples: u64 = 0; + let mut sample_countdown: u64 = 0; + // Borrow each record's bytes in place (see the coordinate ingest loop); // the key is extracted and the bytes copied into the buffer before the // borrow ends, so no owned `RawRecord` is needed here. - while let Some(bam_bytes) = record_source.next_record_borrowed()? { + loop { + // Decide sampling before the fetch, so every segment below is timed + // on the same records or on none. Timing a subset would bias the + // partition toward whichever step happened to be measured, and the + // partition's whole value is that its segments sum to the span. + let sample_this = sample_countdown == 0; + if sample_this { + sample_countdown = ingest_sample_interval - 1; + ingest_samples += 1; + } else { + sample_countdown -= 1; + } + + let t = sample_this.then(Instant::now); + let Some(bam_bytes) = record_source.next_record_borrowed()? else { break }; + if let Some(t0) = t { + ingest_raw.fetch += t0.elapsed().as_secs_f64(); + } + stats.total_records += 1; - progress.log_if_needed(1); + let t = sample_this.then(Instant::now); + progress_batch.tick(&progress); + if let Some(t0) = t { + ingest_raw.tick += t0.elapsed().as_secs_f64(); + } // Extract the full template key, verify the lanes the chosen variant // dropped are constant relative to the first record, then push the // narrowed key. + let t = sample_this.then(Instant::now); let full = extract_template_key_inline(bam_bytes, lib_lookup, self.cell_tag, cb_hasher); - if let Some(violation) = verify_dropped_lanes(&first, &full, variant) { + if let Some(t0) = t { + ingest_raw.key += t0.elapsed().as_secs_f64(); + } + let t = sample_this.then(Instant::now); + let violation = verify_dropped_lanes(&first, &full, variant); + if let Some(t0) = t { + ingest_raw.verify += t0.elapsed().as_secs_f64(); + } + if let Some(violation) = violation { let name = String::from_utf8_lossy( fgumi_raw_bam::RawRecordView::new(bam_bytes).read_name(), ) .into_owned(); return Err(dropped_lane_error(&name, violation)); } + let t = sample_this.then(Instant::now); buffer.push(bam_bytes, K::from_full(&full))?; + if let Some(t0) = t { + ingest_raw.push += t0.elapsed().as_secs_f64(); + } - if probe.should_sample_read(stats.total_records) { + let t = sample_this.then(Instant::now); + let should_probe = probe.should_sample_read(stats.total_records); + let over_limit = buffer.memory_usage() >= self.memory_limit; + if let Some(t0) = t { + ingest_raw.probe += t0.elapsed().as_secs_f64(); + } + if should_probe { probe.log_mid_read(probe_stats(&buffer), Some(pool.phase1_queue_depths())); } // Check memory usage - if buffer.memory_usage() >= self.memory_limit { + if over_limit { timer.end_read_span(); let bstats = probe_stats(&buffer); let depths = Some(pool.phase1_queue_depths()); @@ -4174,6 +4520,7 @@ impl RawExternalSorter { } timer.end_read_span(); + progress_batch.flush(&progress); progress.log_final(); if let Some(err) = record_source.take_error() { return Err(anyhow::Error::from(err)); @@ -4191,6 +4538,18 @@ impl RawExternalSorter { probe.phase1_end(buffer.memory_usage() as u64); // Ingest is done: everything from here is Phase 2 (merge/write). + // Snapshot Phase 1's floor inputs before the pool is handed to the merge: + // the counters describe the phase that has just ended, and reading them + // later would both borrow a moved value and risk folding in Phase 2 work. + let phase1_floor = Phase1FloorInputs { + input_busy_secs: pool.phase1_input_busy_secs(), + threads: self.phase1_threads(), + ingest: pool.phase1_ingest_stats().snapshot(), + sample: ingest_raw, + samples: ingest_samples, + records: stats.total_records, + clock_overhead_nanos, + }; self.enter_output_phase(&pool); if chunk_files.is_empty() { @@ -4251,7 +4610,12 @@ impl RawExternalSorter { if let Ok(pool) = Arc::try_unwrap(pool) { pool.shutdown(); } - timer.log_summary(self.phase1_threads(), self.phase2_threads(), self.max_temp_files); + timer.log_summary( + self.phase1_threads(), + self.phase2_threads(), + self.max_temp_files, + phase1_floor, + ); debug!("Sort complete: {} records processed", stats.total_records); Ok(stats) @@ -5459,6 +5823,7 @@ impl RawExternalSorter { // Snapshot the process-wide presentation counters so the sub-phase log // reports only this merge's delta, not totals accumulated by any prior // (sequential or concurrent) sort sharing the process. + let mut fetch_counts = RecordFetchCounts::default(); let borrowed_before = RECORD_BORROWED.load(std::sync::atomic::Ordering::Relaxed); let reassembled_before = RECORD_REASSEMBLED.load(std::sync::atomic::Ordering::Relaxed); let loop_start = Instant::now(); @@ -5497,7 +5862,8 @@ impl RawExternalSorter { } let t = sample_this.then(Instant::now); - let record_bytes = winner_record_bytes(&sources[src_idx], guard.consumer_ref())?; + let record_bytes = + winner_record_bytes(&sources[src_idx], guard.consumer_ref(), &mut fetch_counts)?; if let Some(t0) = t { merge_present_secs += t0.elapsed().as_secs_f64(); } @@ -5534,6 +5900,9 @@ impl RawExternalSorter { } } + // Before the delta below is read, or this merge reports zero. + fetch_counts.publish(); + let loop_total = loop_start.elapsed().as_secs_f64(); let borrowed_this_merge = RECORD_BORROWED.load(std::sync::atomic::Ordering::Relaxed) - borrowed_before; @@ -5697,6 +6066,7 @@ impl RawExternalSorter { let loop_start = Instant::now(); let mut records_merged: u64 = 0; let mut published_src: Option = None; + let mut fetch_counts = RecordFetchCounts::default(); while tree.winner_is_active() { let winner = tree.winner(); let src_idx = source_map[winner]; @@ -5712,7 +6082,8 @@ impl RawExternalSorter { published_src = Some(src_idx); pool.set_phase2_next_source(tree.runner_up().map(|w| source_map[w])); } - let record_bytes = winner_record_bytes(&sources[src_idx], guard.consumer_ref())?; + let record_bytes = + winner_record_bytes(&sources[src_idx], guard.consumer_ref(), &mut fetch_counts)?; writer.write_raw_record(record_bytes)?; records_merged += 1; merge_progress_batch.tick(&merge_progress); @@ -5734,6 +6105,7 @@ impl RawExternalSorter { // cannot change it, and the consumer's report is harvested before // `finish_output` releases the merge sources and with them the // consumer. Both describe the loop that has just ended. + fetch_counts.publish(); let loop_total = loop_start.elapsed().as_secs_f64(); let active_workers = pool.active_workers(); let stalls = { @@ -6090,6 +6462,7 @@ fn create_raw_bam_reader_pool_integrated>( pool.decompressed_input_done_flag(), pool.input_read_error_flag(), pool.decompress_error_flag(), + pool.phase1_ingest_stats(), ); // Deliberately not phrased as a header failure. The header was parsed @@ -7447,6 +7820,76 @@ mod tests { aux } + /// Build `MC:Z:` aux tag bytes. + fn mc_aux(value: &[u8]) -> Vec { + let mut aux = Vec::new(); + aux.extend_from_slice(b"MCZ"); + aux.extend_from_slice(value); + aux.push(0); // null terminator + aux + } + + /// Overwrite the mate position of a record built by `build_mapped_bam`, + /// which otherwise sets it equal to the record's own position. + fn with_mate_pos(mut bam: Vec, mate_pos: i32) -> Vec { + bam[24..28].copy_from_slice(&mate_pos.to_le_bytes()); + bam + } + + /// The mate lane resolves through `MC`, and it lands where a record whose + /// mate is already at the unclipped position lands. + /// + /// The template-coordinate key is output-identity-critical against + /// `samtools sort`, and the mate lane is the one part of it that comes from + /// parsing a tag rather than from a fixed field offset. Asserting the two + /// keys are equal pins the whole lane -- packing included -- rather than + /// just the parser, which `cigar.rs` already covers. + #[test] + fn test_extract_template_key_mate_lane_comes_from_mc() { + let header = Header::builder().build(); + let lib_lookup = LibraryLookup::from_header(&header); + + // Mate at 100 with 10 leading soft clips: unclipped 5' is 90. + let with_mc = with_mate_pos(build_mapped_bam(0, 50, b"read1", &mc_aux(b"10S40M")), 100); + // The same record whose mate is already reported at 90, and no MC to parse. + let without_mc = with_mate_pos(build_mapped_bam(0, 50, b"read1", &[]), 90); + + let keyed = extract_template_key_inline(&with_mc, &lib_lookup, None, &test_cb_hasher()); + let expected = + extract_template_key_inline(&without_mc, &lib_lookup, None, &test_cb_hasher()); + assert_eq!(keyed, expected, "MC-derived mate lane must equal the unclipped position"); + } + + /// A non-UTF-8 `MC` reaches the parser and its valid prefix still sets the + /// mate lane. + /// + /// Extraction hands `MC` over as raw bytes rather than validating it as + /// UTF-8 first, so this value is parsed where it was previously discarded + /// (leaving the mate lane at the raw mate position). This pins the change at + /// the key level, which is the level the sort order is defined at. + #[test] + fn test_extract_template_key_mate_lane_parses_a_non_utf8_mc_prefix() { + let header = Header::builder().build(); + let lib_lookup = LibraryLookup::from_header(&header); + + let with_bad_mc = + with_mate_pos(build_mapped_bam(0, 50, b"read1", &mc_aux(b"10S40M\xff")), 100); + let unclipped = with_mate_pos(build_mapped_bam(0, 50, b"read1", &[]), 90); + let raw_mate = with_mate_pos(build_mapped_bam(0, 50, b"read1", &[]), 100); + + let keyed = extract_template_key_inline(&with_bad_mc, &lib_lookup, None, &test_cb_hasher()); + assert_eq!( + keyed, + extract_template_key_inline(&unclipped, &lib_lookup, None, &test_cb_hasher()), + "the valid CIGAR prefix must still be applied" + ); + assert_ne!( + keyed, + extract_template_key_inline(&raw_mate, &lib_lookup, None, &test_cb_hasher()), + "discarding the tag would leave the mate lane at the raw mate position" + ); + } + #[test] fn test_extract_template_key_cb_present_has_nonzero_hash() { let header = Header::builder().build(); @@ -9183,7 +9626,7 @@ mod tests { // log_summary must not panic (output goes to log sink). `consolidate_count` // is 1 here, so the consolidation branch is exercised too. - timer.log_summary(4, 4, 64); + timer.log_summary(4, 4, 64, Phase1FloorInputs::default()); } // ======================================================================== diff --git a/crates/fgumi-sort/src/lib.rs b/crates/fgumi-sort/src/lib.rs index b50ea8a42..1add3947c 100644 --- a/crates/fgumi-sort/src/lib.rs +++ b/crates/fgumi-sort/src/lib.rs @@ -94,6 +94,7 @@ pub(crate) mod merge_headroom; pub(crate) mod merge_phases; pub(crate) mod merge_stalls; pub(crate) mod merge_trace; +pub(crate) mod phase1_stats; pub(crate) mod pipeline; pub(crate) mod pooled_bam_writer; pub(crate) mod pooled_chunk_writer; diff --git a/crates/fgumi-sort/src/phase1_stats.rs b/crates/fgumi-sort/src/phase1_stats.rs new file mode 100644 index 000000000..c35ec6839 --- /dev/null +++ b/crates/fgumi-sort/src/phase1_stats.rs @@ -0,0 +1,345 @@ +//! What Phase 1's serial ingest thread waits for. +//! +//! Phase 2 has a floor line ([`crate::merge_headroom`]) because its three +//! limits -- serial consumer, worker capacity, coordination -- imply unrelated +//! fixes and are routinely confused. Phase 1 has the same three limits and had +//! none of the accounting: its report was four wall-clock spans (read, sort, +//! spill write, merge) with no way to tell a thread that is busy from one that +//! is waiting. +//! +//! That gap matters more than it did for the merge. External sampling of a +//! 16-thread whole-genome sort put Phase 1 at **60% of total wall clock with its +//! main thread 91% busy** while all 16 cores averaged 5.3 -- so the phase is +//! bound by one thread's serial CPU, and every worker-side change is pushing on +//! a wall that is not there. In-process numbers should say that without needing +//! a `/proc` sampler attached from outside. +//! +//! # The two waits +//! +//! The ingest thread blocks in exactly two places, and neither was measured: +//! +//! 1. **Waiting for a decompressed block.** [`crate::read_ahead::PooledInputStream`] +//! parks when the next serial it needs has not arrived. Because blocks are +//! consumed in serial order through a reorder buffer, this can fire while +//! other blocks are ready -- head-of-line blocking -- which is a different +//! problem from an empty queue and is counted separately here. +//! 2. **Waiting for the previous spill to finish.** `drain_pending_spill` waits +//! on the prior chunk's write handle between the read span ending and the +//! in-memory sort starting, so that time lands in *no* phase bucket at all +//! and shows up only as an unexplained residual against total wall clock. +//! +//! Timing is exact rather than sampled: both waits are milliseconds against a +//! ~30 ns clock read, so the clock is 5 orders of magnitude below the quantity +//! and costs nothing to read (the same argument [`crate::merge_trace`] makes for +//! the merge's block pull, where an exact timer cost 0.15%). + +use std::sync::atomic::{AtomicU64, Ordering}; + +/// Why the ingest thread parked waiting for a block. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum ParkCause { + /// Nothing decompressed was available at all: the workers are behind. + Starved, + /// Blocks were buffered, but not the serial the consumer needs next. The + /// pipeline has work in flight and the consumer still cannot use it. + HeadOfLine, +} + +/// Counters for the ingest thread's waits, shared with the pool's input stream. +#[derive(Debug, Default)] +pub(crate) struct Phase1IngestStats { + /// Times the ingest thread parked waiting for its next block. + parks: AtomicU64, + /// Total nanoseconds parked waiting for a block. + park_nanos: AtomicU64, + /// Parks where the reorder buffer held nothing at all. + parks_starved: AtomicU64, + /// Parks where the reorder buffer held blocks, but not the next serial. + parks_head_of_line: AtomicU64, + /// Times the ingest thread waited on the previous spill to complete. + spill_waits: AtomicU64, + /// Total nanoseconds waiting on a previous spill. + spill_wait_nanos: AtomicU64, +} + +impl Phase1IngestStats { + /// Record one park and what the reorder buffer looked like when it happened. + pub(crate) fn record_park(&self, elapsed_nanos: u64, cause: ParkCause) { + self.parks.fetch_add(1, Ordering::Relaxed); + self.park_nanos.fetch_add(elapsed_nanos, Ordering::Relaxed); + match cause { + ParkCause::Starved => self.parks_starved.fetch_add(1, Ordering::Relaxed), + ParkCause::HeadOfLine => self.parks_head_of_line.fetch_add(1, Ordering::Relaxed), + }; + } + + /// Record one wait on the previous chunk's spill write. + pub(crate) fn record_spill_wait(&self, elapsed_nanos: u64) { + self.spill_waits.fetch_add(1, Ordering::Relaxed); + self.spill_wait_nanos.fetch_add(elapsed_nanos, Ordering::Relaxed); + } + + /// A consistent view of the counters, for reporting. + pub(crate) fn snapshot(&self) -> Phase1IngestReport { + Phase1IngestReport { + parks: self.parks.load(Ordering::Relaxed), + park_secs: secs(self.park_nanos.load(Ordering::Relaxed)), + parks_starved: self.parks_starved.load(Ordering::Relaxed), + parks_head_of_line: self.parks_head_of_line.load(Ordering::Relaxed), + spill_waits: self.spill_waits.load(Ordering::Relaxed), + spill_wait_secs: secs(self.spill_wait_nanos.load(Ordering::Relaxed)), + } + } +} + +#[allow(clippy::cast_precision_loss, reason = "nanosecond totals stay far below 2^52")] +fn secs(nanos: u64) -> f64 { + nanos as f64 / 1_000_000_000.0 +} + +/// What the ingest thread waited for, as reported. +#[derive(Debug, Clone, Copy, Default)] +pub(crate) struct Phase1IngestReport { + pub(crate) parks: u64, + pub(crate) park_secs: f64, + pub(crate) parks_starved: u64, + pub(crate) parks_head_of_line: u64, + pub(crate) spill_waits: u64, + pub(crate) spill_wait_secs: f64, +} + +impl Phase1IngestReport { + /// Mean park in microseconds, or `None` when it never parked. + /// + /// The mean is the discriminant between "parked rarely and long" (a supply + /// problem) and "parked constantly and briefly" (a handoff problem), which + /// the merge campaign found to be the difference between a fixable stall and + /// an unfixable one. + pub(crate) fn mean_park_micros(&self) -> Option { + #[allow(clippy::cast_precision_loss, reason = "park counts stay far below 2^52")] + (self.parks > 0).then(|| self.park_secs * 1_000_000.0 / self.parks as f64) + } +} + +/// One record in this many is timed for the sub-phase partition. +/// +/// Prime, so the sampled set cannot align with any periodic structure in the +/// input (read groups, tile boundaries, alternating mate records) and bias the +/// partition toward whichever records happen to be cheap. +pub(crate) const INGEST_SAMPLE_INTERVAL: u64 = 1021; + +/// Where the ingest thread's serial CPU goes, per record. +/// +/// The floor line says this thread *is* the limit -- on a 16-thread whole-genome +/// sort it is 137.2s of a 145.7s read span, against a worker-capacity floor of +/// 22.4s -- so the only question left is what the 137.2s is made of. Nothing else +/// in the phase can answer it: worker counters describe the pool, and wall-clock +/// spans describe the phase, and neither looks inside the loop. +/// +/// Sampled rather than timed on every record, for the same reason +/// [`crate::merge_headroom::ConsumerSample`] is: the loop runs at ~175 ns/record +/// and an `Instant::now()` pair costs 15-35 ns on aarch64, so timing five +/// segments on every record would cost more than several of the segments it +/// measures. One record in [`INGEST_SAMPLE_INTERVAL`] is timed and scaled, and +/// the scale is reported next to the result. +/// +/// Each field is **exactly one** timed region in the loop, which is what makes +/// [`Self::corrected`] valid: it subtracts one clock pair per field per sample, +/// so a field spanning two bracketed regions would be under-corrected by a +/// whole pair. That is why progress counting and spill probing are separate +/// fields rather than one "bookkeeping" bucket -- they sit at opposite ends of +/// the loop body and cannot share a bracket. +#[derive(Debug, Clone, Copy, Default)] +pub(crate) struct IngestSample { + /// Pulling the next record's bytes from the pool's decompressed stream. + /// **Includes park time**, so it is not pure CPU -- `park_secs` measures that + /// part exactly and separately. + pub(crate) fetch: f64, + /// Extracting the sort key from the record bytes. + pub(crate) key: f64, + /// Verifying that the lanes the chosen key variant drops are constant. + pub(crate) verify: f64, + /// Copying the record into the arena and appending its ref. + pub(crate) push: f64, + /// Counting the record toward the progress log. + pub(crate) tick: f64, + /// The spill probe's sample check and the memory-limit test that follows it. + pub(crate) probe: f64, +} + +impl IngestSample { + /// Every segment multiplied by the sampling scale. + #[must_use] + pub(crate) fn scaled(self, scale: f64) -> Self { + Self { + fetch: self.fetch * scale, + key: self.key * scale, + verify: self.verify * scale, + push: self.push * scale, + tick: self.tick * scale, + probe: self.probe * scale, + } + } + + /// Every segment with its own measurement overhead removed. + /// + /// One `Instant::now()`/`elapsed()` pair per segment per sampled record, and + /// that pair's cost lands inside the interval it times. Clamped at zero: a + /// segment cheaper than the clock measuring it cannot be resolved this way, + /// and zero says so where a negative would read as a bug. + #[must_use] + pub(crate) fn corrected(self, samples: u64, overhead_nanos: u64) -> Self { + if samples == 0 || overhead_nanos == 0 { + return self; + } + #[expect(clippy::cast_precision_loss, reason = "sample counts stay below 2^52")] + let per_segment = (samples * overhead_nanos) as f64 / 1e9; + let fix = |v: f64| (v - per_segment).max(0.0); + Self { + fetch: fix(self.fetch), + key: fix(self.key), + verify: fix(self.verify), + push: fix(self.push), + tick: fix(self.tick), + probe: fix(self.probe), + } + } + + /// The six segments summed. + #[must_use] + pub(crate) fn total(self) -> f64 { + self.fetch + self.key + self.verify + self.push + self.tick + self.probe + } +} + +/// A scaled ingest sample checked against the read span it should partition. +#[derive(Debug, Clone, Copy)] +pub(crate) struct IngestPartition { + /// Scaled, corrected per-segment seconds. + pub(crate) segments: IngestSample, + /// Measured read span, exact. + pub(crate) read_secs: f64, + /// Park time inside `segments.fetch`, measured exactly and separately. + pub(crate) park_secs: f64, +} + +impl IngestPartition { + /// Read-span time the segments do not account for. + /// + /// **Signed on purpose.** A negative residual means the sample + /// over-attributes -- clock overhead left inside the timed regions, or a + /// sampling bias -- and the merge's first partition did exactly that, + /// summing to 321.5s of a 189.3s loop. Only the sign made it visible; a + /// clamped residual would have reported a tidy zero and the partition would + /// have been believed. + #[must_use] + pub(crate) fn residual_secs(self) -> f64 { + self.read_secs - self.segments.total() + } + + /// Residual as a share of the read span, for judging whether the partition + /// is trustworthy at all. + #[must_use] + pub(crate) fn residual_share(self) -> f64 { + if self.read_secs > 0.0 { self.residual_secs() / self.read_secs } else { 0.0 } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_the_residual_is_signed_so_over_attribution_is_visible() { + // The merge's first partition summed to 321.5s of a 189.3s loop. A + // clamped residual would have shown 0.0 and the numbers would have been + // believed; the sign is what exposed the clock overhead inside them. + let over = IngestPartition { + segments: IngestSample { fetch: 200.0, ..IngestSample::default() }, + read_secs: 100.0, + park_secs: 0.0, + }; + assert!(over.residual_secs() < 0.0, "got {}", over.residual_secs()); + assert!((over.residual_share() + 1.0).abs() < 1e-9, "got {}", over.residual_share()); + } + + #[test] + fn test_clock_correction_subtracts_one_pair_per_segment_per_sample() { + // Ten samples, 20 ns per pair, five segments: each segment carries + // 10 x 20 ns = 200 ns of clock, and each is corrected independently. + let raw = IngestSample { + fetch: 1e-6, + key: 1e-6, + verify: 1e-7, + push: 1e-6, + tick: 1e-6, + probe: 1e-6, + }; + let fixed = raw.corrected(10, 20); + assert!((fixed.fetch - 0.8e-6).abs() < 1e-12, "got {}", fixed.fetch); + // A segment cheaper than the clock that measured it clamps to zero rather + // than going negative, which would read as a bug rather than as + // "unresolvable by this method". + assert!((fixed.verify - 0.0).abs() < 1e-12, "got {}", fixed.verify); + } + + #[test] + fn test_scaling_happens_before_correction_is_meaningful() { + // Scale multiplies the sampled segments up to the whole loop; correction + // works on the sampled scale. Applying them in the wrong order would + // subtract one pair's cost from the *scaled* total rather than from each + // sample, understating the correction by the scale factor. + let raw = IngestSample { key: 2e-6, ..IngestSample::default() }; + let corrected_then_scaled = raw.corrected(10, 20).scaled(1000.0); + let scaled_then_corrected = raw.scaled(1000.0).corrected(10, 20); + assert!(corrected_then_scaled.key < scaled_then_corrected.key); + } + + #[test] + fn test_park_causes_are_counted_separately() { + let stats = Phase1IngestStats::default(); + stats.record_park(1_000, ParkCause::Starved); + stats.record_park(3_000, ParkCause::HeadOfLine); + stats.record_park(6_000, ParkCause::HeadOfLine); + + let report = stats.snapshot(); + assert_eq!(report.parks, 3); + assert_eq!(report.parks_starved, 1); + assert_eq!(report.parks_head_of_line, 2); + // Starvation and head-of-line blocking have different fixes, so a report + // that only totalled them would not distinguish the two. + assert_eq!(report.parks_starved + report.parks_head_of_line, report.parks); + assert!((report.park_secs - 10e-6).abs() < 1e-12, "got {}", report.park_secs); + } + + #[test] + fn test_the_spill_handoff_wait_is_counted_on_its_own() { + let stats = Phase1IngestStats::default(); + stats.record_park(2_000_000_000, ParkCause::Starved); + stats.record_spill_wait(3_000_000_000); + + let report = stats.snapshot(); + // The spill wait sits between the read span ending and the sort starting, + // so it is in no phase bucket. It is reported separately from the park + // rather than summed into it: one is the pool failing to keep the ingest + // thread fed, the other is the previous chunk's write not being done, and + // a single "blocked" total would hide which. + assert!((report.spill_wait_secs - 3.0).abs() < 1e-9, "got {}", report.spill_wait_secs); + assert!((report.park_secs - 2.0).abs() < 1e-9, "got {}", report.park_secs); + assert_eq!(report.spill_waits, 1); + } + + #[test] + fn test_mean_park_is_absent_rather_than_zero_when_it_never_parked() { + let idle = Phase1IngestStats::default().snapshot(); + // A reported 0 us mean would read as "parked, instantly", which is the + // opposite of what never parking means. + assert_eq!(idle.mean_park_micros(), None); + + let stats = Phase1IngestStats::default(); + stats.record_park(4_000_000, ParkCause::Starved); + stats.record_park(6_000_000, ParkCause::Starved); + let mean = stats.snapshot().mean_park_micros().expect("parked twice"); + assert!((mean - 5_000.0).abs() < 1e-6, "got {mean}"); + } +} diff --git a/crates/fgumi-sort/src/read_ahead.rs b/crates/fgumi-sort/src/read_ahead.rs index 2d391ee11..f3cfa69b3 100644 --- a/crates/fgumi-sort/src/read_ahead.rs +++ b/crates/fgumi-sort/src/read_ahead.rs @@ -291,6 +291,8 @@ pub struct PooledInputStream { current_buf: Vec, /// Read position within `current_buf`. current_pos: usize, + /// Where this stream records the ingest thread's waits. + stats: std::sync::Arc, /// Reusable scratch buffer for records (or their length prefixes) that /// straddle a decompressed-block boundary and therefore cannot be borrowed /// directly out of `current_buf`. See [`PooledInputStream::next_record_borrowed`]. @@ -305,12 +307,14 @@ impl PooledInputStream { decompressed_input_done: std::sync::Arc, input_read_error: std::sync::Arc, decompression_error: std::sync::Arc, + stats: std::sync::Arc, ) -> Self { Self { decompressed_input, decompressed_input_done, input_read_error, decompression_error, + stats, reorder: fgumi_bam_io::ReorderBuffer::new(), current_buf: Vec::new(), current_pos: 0, @@ -383,8 +387,27 @@ impl PooledInputStream { return None; } - // Park until a worker pushes a block and calls unpark() + // Park until a worker pushes a block and calls unpark(). + // + // Timed exactly rather than sampled: a park is microseconds to + // milliseconds against a ~30 ns clock read, so the clock is orders of + // magnitude below the quantity it measures. The cause is captured + // before parking because it is not recoverable afterwards -- and the + // two causes have different fixes. An empty reorder buffer means the + // workers are behind; a non-empty one means blocks are ready and the + // serial order will not let the consumer have them, which no amount + // of extra decompression capacity would help. + let cause = if self.reorder.buffer_len() == 0 { + crate::phase1_stats::ParkCause::Starved + } else { + crate::phase1_stats::ParkCause::HeadOfLine + }; + let parked_at = std::time::Instant::now(); std::thread::park(); + self.stats.record_park( + u64::try_from(parked_at.elapsed().as_nanos()).unwrap_or(u64::MAX), + cause, + ); // After waking, check for errors before looping back to drain. // A worker may have set an error flag instead of pushing a block. @@ -806,6 +829,7 @@ mod tests { Arc::new(AtomicBool::new(true)), // decompressed_input_done Arc::new(AtomicBool::new(false)), // input_read_error Arc::new(AtomicBool::new(false)), // decompression_error + Arc::default(), // ingest wait counters ) } @@ -913,6 +937,7 @@ mod tests { Arc::new(AtomicBool::new(true)), Arc::new(AtomicBool::new(false)), Arc::new(AtomicBool::new(false)), + Arc::default(), ); let err = pooled.next_record_borrowed().expect_err("truncated body should error"); assert_eq!(err.kind(), std::io::ErrorKind::UnexpectedEof); diff --git a/crates/fgumi-sort/src/worker_pool.rs b/crates/fgumi-sort/src/worker_pool.rs index 97e994f94..cb78b70ed 100644 --- a/crates/fgumi-sort/src/worker_pool.rs +++ b/crates/fgumi-sort/src/worker_pool.rs @@ -1432,6 +1432,10 @@ pub(crate) struct SharedPipelineState { /// `do_shutdown` checks join results and sets this flag so the main thread /// does not park forever waiting for work that will never arrive. pub(crate) worker_panicked: Arc, + /// What Phase 1's serial ingest thread waited for. Owned here because the + /// input stream records into it and the sorter reports from it, and neither + /// owns the other. + pub(crate) phase1_ingest: Arc, /// Next serial for input block reading (atomic increment for ordering). input_read_serial: AtomicU64, /// Raw input blocks: `ReadInputBlocks` → `DecompressInput`. @@ -1627,6 +1631,7 @@ impl SharedPipelineState { decompression_error: Arc::new(AtomicBool::new(false)), chunk_read_error: Arc::new(AtomicBool::new(false)), worker_panicked: Arc::new(AtomicBool::new(false)), + phase1_ingest: Arc::new(crate::phase1_stats::Phase1IngestStats::default()), input_read_serial: AtomicU64::new(0), raw_input_blocks: Arc::new(ArrayQueue::new(data_queue_cap)), decompressed_input: Arc::new(ArrayQueue::new(data_queue_cap)), @@ -3521,6 +3526,28 @@ impl SortWorkerPool { self.shared.merge_phases.snapshot() } + /// Worker seconds spent feeding Phase 1's ingest thread: reading raw input + /// blocks off disk and decompressing them. + /// + /// Deliberately excludes `Compress`, which serves spill *and* output and so + /// cannot be attributed to the ingest span; the spill half is reported on its + /// own through [`crate::merge_phases`]. + pub(crate) fn phase1_input_busy_secs(&self) -> f64 { + let ns = self.pipeline_stats.step_ns[SortStep::ReadInputBlocks as usize] + .load(Ordering::Relaxed) + + self.pipeline_stats.step_ns[SortStep::DecompressInput as usize] + .load(Ordering::Relaxed); + #[allow(clippy::cast_precision_loss, reason = "nanosecond totals stay far below 2^52")] + { + ns as f64 / 1_000_000_000.0 + } + } + + /// Counters for what Phase 1's ingest thread waited on. + pub(crate) fn phase1_ingest_stats(&self) -> Arc { + Arc::clone(&self.shared.phase1_ingest) + } + /// The pool's shared state, for the merge consumer's own instrumentation. /// /// The consumer runs on the main thread and is not a pool worker, but it is diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index 75e7e3c45..2c1ef7084 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -1308,7 +1308,7 @@ mod tests { assert_eq!(supp.template_length(), -351); assert_ne!(supp.flags() & rflags::MATE_REVERSE, 0, "mate is reverse"); assert_eq!(supp.flags() & rflags::MATE_UNMAPPED, 0, "mate is mapped"); - assert_eq!(supp.tags().find_mc(), Some("50M")); + assert_eq!(supp.tags().find_mc(), Some(b"50M".as_slice())); assert_eq!(supp.tags().find_int(SamTag::MQ), Some(40)); } diff --git a/src/lib/grouper.rs b/src/lib/grouper.rs index bc988517c..dbbe05e57 100644 --- a/src/lib/grouper.rs +++ b/src/lib/grouper.rs @@ -419,7 +419,7 @@ impl RecordPositionGrouper { /// record's aux data: the Decode step already resolved the mate position from /// the `MC` tag, and falls back to a single-ended key when it could not (see /// `compute_group_key_from_raw`). Re-deriving that here cost a second full - /// aux-TLV walk plus a UTF-8 validation, per record, on this serial step. + /// aux-TLV walk per record, on this serial step. /// /// This makes validation agree with the value grouping actually uses. It also /// makes fgumi marginally more permissive for one malformed shape: a record diff --git a/src/lib/sam/mod.rs b/src/lib/sam/mod.rs index f46b63935..f1b70a864 100644 --- a/src/lib/sam/mod.rs +++ b/src/lib/sam/mod.rs @@ -35,7 +35,7 @@ pub use fgumi_sam::builder::{ pub use fgumi_sam::record_utils::{ PairOrientation, alignment_end, cigar_reference_length, get_pair_orientation, is_fr_pair, is_fr_pair_from_tags, leading_clipping, leading_soft_clipping, mate_unclipped_end, - mate_unclipped_start, parse_cigar_string, read_pos_at_ref_pos, reference_length, - trailing_clipping, trailing_soft_clipping, unclipped_end, unclipped_five_prime_position, - unclipped_start, + mate_unclipped_start, parse_cigar_bytes, parse_cigar_string, read_pos_at_ref_pos, + reference_length, trailing_clipping, trailing_soft_clipping, unclipped_end, + unclipped_five_prime_position, unclipped_start, };