From d53a2388f125a23dfdd4f6817de83ed0513e80d5 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 09:30:49 -0700 Subject: [PATCH 1/8] feat(micromzpaf): pack IonAnnot into a u32 and add neutral losses IonAnnot becomes a packed `u32` instead of a struct of fields, so equality is one word compare and `(IonAnnot, f32)` stays 8 bytes. A spectral library carries one annotation per fragment, so the type is replicated millions of times in a loaded arena. The new capacity is spent on the annotations an mzSpecLib reader needs: neutral losses, internal fragments, bare immonium ions, and the mass-error suffix. Losses are keyed by atomic composition rather than by spelling, because libraries write the same chemical loss different ways (`-CH3SOH` in NIST, `-CH4OS` in SpectraST) and fragment labels must be unique within a precursor. Losses outside the table are rejected, never coerced onto a nearby representable ion. Drops the rustyms dependency: `from_fragment` had no callers. `UnknownIonCounter` replaces the counter three readers each hand-rolled. Skyline's used `saturating_add`, so past 255 unknown ions it reissued `?255` and made every later peak carrying that label unreachable. `IonSeriesTerminality` is gone -- nothing called `terminality()`. `Series` takes its place in the re-export and resolves the standalone-series-enum TODO in `FragmentLabel`. --- Cargo.lock | 2 +- rust/micromzpaf/Cargo.toml | 6 +- rust/micromzpaf/src/lib.rs | 1200 +++++++++++++------ rust/micromzpaf/src/loss.rs | 449 +++++++ rust/timsquery/src/lib.rs | 3 +- rust/timsquery/src/serde/diann_io.rs | 11 +- rust/timsquery/src/serde/skyline_io.rs | 6 +- rust/timsquery/src/serde/spectronaut_io.rs | 6 +- rust/timsquery/src/traits/fragment_label.rs | 12 +- rust/timsseek/src/lib.rs | 1 - 10 files changed, 1325 insertions(+), 371 deletions(-) create mode 100644 rust/micromzpaf/src/loss.rs diff --git a/Cargo.lock b/Cargo.lock index 2eb4c5b7..dc2df80d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4019,8 +4019,8 @@ dependencies = [ name = "micromzpaf" version = "0.33.0" dependencies = [ - "rustyms", "serde", + "serde_json", "thiserror 2.0.18", ] diff --git a/rust/micromzpaf/Cargo.toml b/rust/micromzpaf/Cargo.toml index 83c47f4c..e168ab37 100644 --- a/rust/micromzpaf/Cargo.toml +++ b/rust/micromzpaf/Cargo.toml @@ -8,5 +8,7 @@ license.workspace = true serde = { workspace = true, features = ["derive"] } thiserror = { workspace = true } -# Workspace-inherited deps -rustyms = { workspace = true } +[dev-dependencies] +# Only to pin that `IonAnnot` serialises as its mzPAF string rather than the +# packed word. +serde_json = { workspace = true } diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 551e4b24..9301ff4a 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -1,19 +1,52 @@ //! Compact representation of fragment ion annotations for mass spectrometry. //! -//! This crate provides a minimal, memory-efficient encoding of fragment ion annotations -//! that pack ion series, charge state, ordinal position, and isotope offset into a -//! compact structure suitable for high-throughput proteomics workflows. +//! A spectral library carries one annotation per fragment, so this type is +//! replicated millions of times in a loaded arena and compared on the scoring +//! hot path. It is therefore a packed `u32` rather than a struct of fields: +//! equality is a single word compare, and `(IonAnnot, f32)` stays 8 bytes with +//! no padding. //! -//! # mzPAF Format Compliance +//! # Bit layout //! -//! This implementation supports a subset of the mzPAF annotation format: -//! - Common ion series (a, b, c, d, v, w, x, y, z, precursor) -//! - Charge states (^N notation) -//! - Positive isotope offsets (+Ni notation) +//! ```text +//! bit: 31 30 29 18 17 12 11 8 7 4 3 0 +//! ┌───────┬────────────────┬──────────┬────────┬────────┬───────┐ +//! │ spare │ payload │ loss │isotope │ charge │ kind │ +//! │ 2b │ 12b │ 6b │ 4b zz │ 4b zz │ 4b │ +//! └───────┴────────────────┴──────────┴────────┴────────┴───────┘ +//! ``` +//! +//! `payload` is reinterpreted per `kind` -- a tagged union inside the word: +//! +//! | kind | payload | +//! |---|---| +//! | backbone (a/b/c/d/v/w/x/y/z), unknown | ordinal, 8b (0..=255) | +//! | internal | start 6b │ end 6b (peptides to 63 residues) | +//! | immonium | residue index 5b | +//! | precursor | unused | +//! +//! `unknown` is discriminant 0 and `charge` is stored biased by one, so the +//! all-zero word is the valid annotation `?0` at charge 1. That matters +//! because `IonAnnot: Default` is forced by `tinyvec::Array` and a default can +//! reach any serde path. //! -//! NOTABLY boes not support: -//! - Negative isotope offsets (not yet implemented) -//! - Complex neutral losses or modifications +//! `charge` and `isotope` are zigzag-encoded so they stay signed in 4 bits. +//! Their ranges (±7) are far wider than anything observed: the HUPO-PSI corpus +//! tops out at charge 3 and isotope 3, with no negative charges at all. Because +//! the field truncates rather than wrapping loudly, every constructor +//! range-checks -- see [`IonAnnot::try_new`]. +//! +//! # mzPAF compliance +//! +//! Supported: the a/b/c/d/v/w/x/y/z series, precursor (`p`), unknown (`?`), +//! internal fragments (`m:`), bare immonium (`IA`), charge (`^N`), +//! positive isotopes (`+Ni`), neutral losses from [`NeutralLoss`], and the +//! mass-error suffix (`/-0.0003`, `/1.2ppm`). +//! +//! Not supported: negative isotope offsets, modified immonium +//! (`IC[Carbamidomethyl]` carries an arbitrary mod string), and losses outside +//! the [`NeutralLoss`] table. These are reported as errors, never coerced into +//! a nearby representable ion. //! //! # Examples //! @@ -29,44 +62,79 @@ //! assert_eq!(ion.get_charge(), 3); //! ``` -use rustyms::fragment::FragmentType; +pub mod loss; + +pub use loss::NeutralLoss; use serde::{ Deserialize, Serialize, }; use std::fmt::Display; use std::hash::Hash; -use std::str::FromStr; use thiserror::Error; -/// Compact representation of fragment annotations. -/// -/// This is a very compressed representation of a fragment -/// ion annotation. Essentially we are packing in 32 bytes -/// the ion series (b, y, ...), charge (+1 / -1 ...), -/// ordinal (12 in the ion series) and isotope. -/// -/// It is not meant to represent all possible ions but rather have -/// a very compact representation of the common ones. -/// -/// # Invariants -/// -/// - **charge**: Must be non-zero (±1 to ±127). Zero charge is invalid and rejected by constructors. -/// - **ordinal**: Limited to u8 range (1-255). Peptides with >255 residues cannot be represented. -/// - **isotope**: Isotope offset relative to monoisotopic peak (M+0), range -128 to +127. +// ── Bit layout ─────────────────────────────────────────────────────────────── + +const KIND_SHIFT: u32 = 0; +const KIND_BITS: u32 = 4; +const CHARGE_SHIFT: u32 = 4; +const CHARGE_BITS: u32 = 4; +const ISOTOPE_SHIFT: u32 = 8; +const ISOTOPE_BITS: u32 = 4; +const LOSS_SHIFT: u32 = 12; +const LOSS_BITS: u32 = 6; +const PAYLOAD_SHIFT: u32 = 18; +const PAYLOAD_BITS: u32 = 12; + +/// Widest charge the 4-bit zigzag field holds. Observed maximum is 3. +pub const CHARGE_MIN: i8 = -7; +pub const CHARGE_MAX: i8 = 7; +/// Widest isotope offset the 4-bit zigzag field holds. Observed maximum is 3. +pub const ISOTOPE_MIN: i8 = -7; +pub const ISOTOPE_MAX: i8 = 7; +/// Width of each internal-fragment endpoint inside `payload`. +const INTERNAL_POS_BITS: u32 = 6; +/// Width of the immonium residue index inside `payload`. +const IMMONIUM_BITS: u32 = 5; +/// Widest residue index an internal fragment endpoint holds. +pub(crate) const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; + +#[inline] +const fn mask(bits: u32) -> u32 { + (1u32 << bits) - 1 +} +/// Zigzag: map a small signed value onto an unsigned one without losing the +/// sign bit to the field width. +#[inline] +const fn zigzag(v: i8) -> u32 { + (((v as i32) << 1) ^ ((v as i32) >> 31)) as u32 +} +#[inline] +const fn unzigzag(u: u32) -> i8 { + (((u >> 1) as i32) ^ -((u & 1) as i32)) as i8 +} + +/// Charge is stored biased by one, so the zero field decodes to charge 1. /// -/// # Memory Layout +/// Charge 0 is rejected by every constructor, so it is not a value the field +/// needs to represent -- and spending the zero word on it would make +/// `IonAnnot::default()` render an annotation that cannot be parsed back. +/// `CHARGE_MIN..=CHARGE_MAX` minus one still zigzags inside 4 bits. +#[inline] +const fn zigzag_charge(charge: i8) -> u32 { + zigzag(charge - 1) +} +#[inline] +const fn unzigzag_charge(u: u32) -> i8 { + unzigzag(u) + 1 +} + +/// Compact representation of fragment annotations. /// -/// The struct fits in 32 bytes with optimal packing: -/// - `IonSeriesOrdinal`: 2 bytes (enum discriminant + u8 ordinal) -/// - `charge`: 1 byte (i8) -/// - `isotope`: 1 byte (i8) +/// A packed `u32`; see the crate docs for the bit layout. Ordering is by the +/// packed word, not field-by-field. #[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] -pub struct IonAnnot { - series_ordinal: IonSeriesOrdinal, - charge: i8, - isotope: i8, -} +pub struct IonAnnot(u32); impl Serialize for IonAnnot { fn serialize(&self, serializer: S) -> Result @@ -83,6 +151,9 @@ impl Serialize for IonAnnot { /// b12+i^3 -> b12 charge 3 isotope 1 /// b12+3i^3 -> b12 charge 3 isotope 2 /// b13 -> b13 (implicit charge 1 and isotope 0) +/// +/// The wire format is the mzPAF string, not the packed word, so the bit layout +/// can change without breaking existing files. impl<'de> Deserialize<'de> for IonAnnot { fn deserialize(deserializer: D) -> Result where @@ -94,66 +165,241 @@ impl<'de> Deserialize<'de> for IonAnnot { } impl IonAnnot { + /// Build a backbone / precursor / unknown annotation. + /// + /// Errors when `charge` is zero or either of `charge`/`isotope` falls + /// outside the packed field's range. That range check is load-bearing: the + /// bit field truncates silently, so an unchecked value would corrupt the + /// annotation rather than fail. pub fn try_new( ion_type: char, ordinal: Option, charge: i8, isotope: i8, ) -> Result { - Ok(Self { - series_ordinal: IonSeriesOrdinal::try_new(ion_type, ordinal)?, + Self::try_new_with_loss(ion_type, ordinal, charge, isotope, NeutralLoss::None) + } + + /// As [`Self::try_new`], carrying a neutral loss. + pub fn try_new_with_loss( + ion_type: char, + ordinal: Option, + charge: i8, + isotope: i8, + loss: NeutralLoss, + ) -> Result { + Self::pack( + IonSeriesOrdinal::from_series_char(ion_type, ordinal)?, + loss, charge, isotope, - }) + ) } - pub fn from_fragment( - frag: FragmentType, + /// Build an internal fragment spanning residues `start..=end`. + /// + /// Endpoints are capped at [`INTERNAL_POS_MAX`], narrower than a backbone + /// ordinal, because two of them share the 12-bit payload. Internal + /// fragments are bounded by peptide length, so 63 is well past tryptic. + pub fn try_new_internal( + start: u8, + end: u8, charge: i8, isotope: i8, + loss: NeutralLoss, ) -> Result { - Ok(Self { - series_ordinal: IonSeriesOrdinal::try_from(frag)?, + if start > INTERNAL_POS_MAX || end > INTERNAL_POS_MAX { + return Err(IonParsingError::OrdinalOutOfRange { + ordinal: start.max(end), + series: 'm', + }); + } + Self::pack( + IonSeriesOrdinal::internal { start, end }, + loss, charge, isotope, - }) + ) } - pub fn terminality(&self) -> IonSeriesTerminality { - self.series_ordinal.terminality() + /// Build a bare immonium ion for an uppercase residue code. + /// + /// Modified immonium (`IC[Carbamidomethyl]`) is not representable at any + /// field width and is rejected by the parser rather than silently losing + /// the modification. + pub fn try_new_immonium( + residue: char, + charge: i8, + isotope: i8, + loss: NeutralLoss, + ) -> Result { + if !residue.is_ascii_uppercase() { + return Err(IonParsingError::UnsupportedFragmentType { + fragment_type: residue, + }); + } + Self::pack( + IonSeriesOrdinal::immonium { residue }, + loss, + charge, + isotope, + ) } - pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { - let new_isotope = - self.isotope - .checked_add(offset_neutrons) - .ok_or_else(|| IonParsingError::Custom { - error: format!( - "Isotope offset overflow: {} + {} exceeds i8 range", - self.isotope, offset_neutrons - ), - })?; + fn pack( + series: IonSeriesOrdinal, + loss: NeutralLoss, + charge: i8, + isotope: i8, + ) -> Result { + if charge == 0 { + return Err(IonParsingError::ChargeCannotBeZero); + } + if !(CHARGE_MIN..=CHARGE_MAX).contains(&charge) { + return Err(IonParsingError::ChargeOutOfRange { charge }); + } + if !(ISOTOPE_MIN..=ISOTOPE_MAX).contains(&isotope) { + return Err(IonParsingError::IsotopeOutOfRange { isotope }); + } + let (kind, payload) = series.to_parts(); + debug_assert!(payload <= mask(PAYLOAD_BITS), "payload overflows its field"); + Ok(IonAnnot( + (kind << KIND_SHIFT) + | ((zigzag_charge(charge) & mask(CHARGE_BITS)) << CHARGE_SHIFT) + | ((zigzag(isotope) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT) + | ((loss as u32 & mask(LOSS_BITS)) << LOSS_SHIFT) + | ((payload & mask(PAYLOAD_BITS)) << PAYLOAD_SHIFT), + )) + } - Ok(Self { - series_ordinal: self.series_ordinal, - charge: self.charge, - isotope: new_isotope, - }) + #[inline] + fn payload(self) -> u32 { + (self.0 >> PAYLOAD_SHIFT) & mask(PAYLOAD_BITS) } + #[inline] pub fn get_charge(&self) -> i8 { - self.charge + unzigzag_charge((self.0 >> CHARGE_SHIFT) & mask(CHARGE_BITS)) + } + + #[inline] + pub fn get_isotope(&self) -> i8 { + unzigzag((self.0 >> ISOTOPE_SHIFT) & mask(ISOTOPE_BITS)) + } + + /// The neutral loss this ion carries, [`NeutralLoss::None`] if it carries + /// none. + #[inline] + pub fn loss(&self) -> NeutralLoss { + NeutralLoss::from_discriminant(((self.0 >> LOSS_SHIFT) & mask(LOSS_BITS)) as u8) + } + + /// Shift the isotope by `offset_neutrons`. + /// + /// Errors when the result leaves [`ISOTOPE_MIN`]..=[`ISOTOPE_MAX`]. That + /// bound is an order of magnitude past any observed isotope offset. + pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { + // Saturating rather than checked: the field is far narrower than `i8`, + // so the range check below is the one that matters and it reports the + // bound that actually applies. + let new_isotope = self.get_isotope().saturating_add(offset_neutrons); + if !(ISOTOPE_MIN..=ISOTOPE_MAX).contains(&new_isotope) { + return Err(IonParsingError::IsotopeOutOfRange { + isotope: new_isotope, + }); + } + Ok(IonAnnot( + (self.0 & !(mask(ISOTOPE_BITS) << ISOTOPE_SHIFT)) + | ((zigzag(new_isotope) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT), + )) } + /// The series ordinal, for the kinds that have one. + /// + /// `None` for precursor, unknown, internal and immonium: `?1` is a + /// uniqueness counter, not a position in a ladder. pub fn try_get_ordinal(&self) -> Option { - self.series_ordinal.try_get_ordinal() + use IonSeriesOrdinal as S; + match self.series_ordinal() { + S::backbone { ordinal, .. } => Some(ordinal), + S::unknown { .. } | S::precursor | S::internal { .. } | S::immonium { .. } => None, + } + } + + /// The backbone series this annotation belongs to, if any. + /// + /// `None` for precursor, unknown, internal and immonium ions, none of + /// which sit on a backbone ladder. + pub fn try_get_series(&self) -> Option { + match self.series_ordinal() { + IonSeriesOrdinal::backbone { series, .. } => Some(series), + _ => None, + } + } + + /// The logical series-and-payload view of this annotation. + pub fn series_ordinal(&self) -> IonSeriesOrdinal { + IonSeriesOrdinal::from_parts((self.0 >> KIND_SHIFT) & mask(KIND_BITS), self.payload()) } } -impl TryFrom<&str> for IonAnnot { - type Error = IonParsingError; +/// A mass-error suffix on an mzPAF annotation: observed minus theoretical. +/// +/// mzSpecLib peak lists carry *observed* m/z, so the theoretical value a +/// library reader wants is `observed - error`. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum MassError { + /// Absolute, in daltons. + Da(f64), + /// Relative, in parts per million. + Ppm(f64), +} - fn try_from(value: &str) -> Result { +impl MassError { + /// Recover the theoretical m/z from the observed one. + pub fn theoretical_from_observed(&self, observed: f64) -> f64 { + match self { + MassError::Da(d) => observed - d, + MassError::Ppm(p) => observed / (1.0 + p * 1e-6), + } + } +} + +/// Split the trailing `/[ppm]` off an annotation, if present. +/// +/// Public because a library reader needs the error even when the ion itself is +/// unrepresentable: the error is what recovers theoretical m/z, so a peak that +/// ends up with an unknown label still gets a correct mass. +pub fn split_mass_error(s: &str) -> Result<(&str, Option), IonParsingError> { + let Some((head, tail)) = s.rsplit_once('/') else { + return Ok((s, None)); + }; + let (num, is_ppm) = match tail.strip_suffix("ppm") { + Some(n) => (n, true), + None => (tail, false), + }; + let v: f64 = num.parse().map_err(|_| IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Unable to parse the mass-error suffix"), + })?; + Ok(( + head, + Some(if is_ppm { + MassError::Ppm(v) + } else { + MassError::Da(v) + }), + )) +} + +impl IonAnnot { + /// Parse the ion itself, with no mass-error suffix. + /// + /// A comma-separated list of alternatives is NOT handled here -- that is + /// ambiguity, and resolving it needs the caller's policy. Split on `,` and + /// parse each alternative. + fn parse_ion(value: &str) -> Result { + // charge: trailing ^N let (rest, charge) = match value.split_once('^') { Some((rest, charge)) => { let charge = charge @@ -167,20 +413,15 @@ impl TryFrom<&str> for IonAnnot { None => (value, 1), }; - // TODO: Implement 'negative isotopes' parsing ... right now I dont use them - // for serialization ... - // Note that this is not 100% compliant with mzPAF + // isotope: +Ni. Negative isotope offsets are not supported. let (rest, isotope) = match rest.split_once('+') { Some((rest, adducts)) => { - // Make sure the adduct is only +{number}?i let adducts = adducts .strip_suffix('i') .ok_or(IonParsingError::ParsingError { error: adducts.to_string(), context: Some("Unsupported adduct found"), })?; - // If its empty its an implicit 1 isotope - // Since we stripped the 'i' from '+i' let isotope = if adducts.is_empty() { 1 } else { @@ -195,33 +436,99 @@ impl TryFrom<&str> for IonAnnot { } None => (rest, 0), }; - let series_ord = IonSeriesOrdinal::from_str(rest)?; - if charge == 0 { - return Err(IonParsingError::ParsingError { + + // neutral loss: everything from the first '-'. Internal fragments use + // 'm:' and never contain '-' before the loss, so splitting + // at the first '-' is unambiguous. + let (core, loss) = match rest.split_once('-') { + Some((core, loss_expr)) => { + let loss = NeutralLoss::from_expression(loss_expr)?.ok_or_else(|| { + IonParsingError::UnsupportedNeutralLoss { + loss: loss_expr.to_string(), + } + })?; + (core, loss) + } + None => (rest, NeutralLoss::None), + }; + + // internal fragment: m: + if let Some(spans) = core.strip_prefix('m') + && let Some((a, b)) = spans.split_once(':') + { + let start = a.parse::().map_err(|_| IonParsingError::ParsingError { error: value.to_string(), - context: Some("Charge cannot be 0"), - }); + context: Some("Unable to parse internal-fragment start"), + })?; + let end = b.parse::().map_err(|_| IonParsingError::ParsingError { + error: value.to_string(), + context: Some("Unable to parse internal-fragment end"), + })?; + return Self::try_new_internal(start, end, charge, isotope, loss); } - Ok(Self { - series_ordinal: series_ord, - charge, - isotope, - }) + + // immonium: I, bare only. + if let Some(res) = core.strip_prefix('I') { + let mut ch = res.chars(); + return match (ch.next(), ch.next()) { + (Some(r), None) => Self::try_new_immonium(r, charge, isotope, loss), + // `IC[Carbamidomethyl]` and friends carry a mod string that no + // fixed-width field can hold. + _ => Err(IonParsingError::UnsupportedModifiedImmonium { + annotation: value.to_string(), + }), + }; + } + + // Backbone / precursor / unknown: a series char then an ordinal. + let mut chars = core.chars(); + let series = chars.next().ok_or(IonParsingError::ParsingError { + error: value.to_string(), + context: Some("Empty string"), + })?; + let rest = chars.as_str(); + let ordinal = if rest.is_empty() { + None + } else { + Some( + rest.parse::() + .map_err(|e| IonParsingError::ParsingError { + error: format!("{rest} -> {e:?}"), + context: Some("Unable to parse the ordinal number"), + })?, + ) + }; + Self::try_new_with_loss(series, ordinal, charge, isotope, loss) + } +} + +impl TryFrom<&str> for IonAnnot { + type Error = IonParsingError; + + /// Parses an annotation, discarding any mass-error suffix. Use + /// [`split_mass_error`] first to keep it. + fn try_from(value: &str) -> Result { + Self::parse_ion(split_mass_error(value)?.0) } } impl Display for IonAnnot { + /// Renders the canonical mzPAF spelling. Not byte-inverse to parsing: a + /// non-canonical loss spelling (`-CH3SOH`) renders canonically (`-CH4OS`). + /// The mass-error suffix is not part of the annotation and is not rendered. fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "{}", self.series_ordinal)?; + write!(f, "{}", self.series_ordinal())?; + write!(f, "{}", self.loss())?; - match self.isotope { + match self.get_isotope() { 0 => {} 1 => write!(f, "+i")?, i => write!(f, "+{}i", i)?, } - if self.charge != 1 { - write!(f, "^{}", self.charge)?; + let charge = self.get_charge(); + if charge != 1 { + write!(f, "^{}", charge)?; } Ok(()) @@ -230,322 +537,513 @@ impl Display for IonAnnot { #[derive(Debug, Error)] pub enum IonParsingError { - #[error("Ordinal {ordinal} out of range for series '{series:?}'")] - OrdinalOutOfRange { ordinal: i32, series: Option }, + #[error("Ordinal {ordinal} out of range for series '{series}'")] + OrdinalOutOfRange { ordinal: u8, series: char }, + #[error("Series '{series}' requires an ordinal")] + MissingOrdinal { series: char }, #[error("Unsupported fragment type: '{fragment_type}'")] UnsupportedFragmentType { fragment_type: char }, + #[error("Charge cannot be 0")] + ChargeCannotBeZero, + #[error("Charge {charge} outside the representable range")] + ChargeOutOfRange { charge: i8 }, + #[error("Isotope offset {isotope} outside the representable range")] + IsotopeOutOfRange { isotope: i8 }, + #[error("Neutral loss '{loss}' is not representable")] + UnsupportedNeutralLoss { loss: String }, + #[error("Modified immonium ions are not representable: '{annotation}'")] + UnsupportedModifiedImmonium { annotation: String }, + #[error("Ran out of distinct unknown-ion labels: the 8-bit ordinal is exhausted")] + UnknownIonsExhausted, #[error("Parsing error: {error}{}", .context.map(|c| format!(" ({})", c)).unwrap_or_default())] ParsingError { error: String, context: Option<&'static str>, }, - #[error("{error}")] - Custom { error: String }, } -/// Refers to what terminus of the original peptide retains the -/// charge after a fragmentation. -#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] -pub enum IonSeriesTerminality { - NTerm, - CTerm, - None, +/// Hands out `?1`, `?2`, ... for peaks whose annotation this crate cannot +/// represent. +/// +/// Fragment labels must be unique within a precursor -- lookup is by first +/// match, so a repeated label makes every later peak carrying it unreachable. +/// A monotonic counter makes that uniqueness structural, and returning an +/// error once the 8-bit ordinal is spent keeps the overflow from being +/// something each reader has to remember to check. +#[derive(Debug, Default, Clone, Copy)] +pub struct UnknownIonCounter(u8); + +impl UnknownIonCounter { + pub fn new() -> Self { + Self::default() + } + + /// The next unused unknown label at `charge`. + pub fn next(&mut self, charge: i8) -> Result { + let ordinal = self + .0 + .checked_add(1) + .ok_or(IonParsingError::UnknownIonsExhausted)?; + let annot = IonAnnot::try_new('?', Some(ordinal), charge, 0)?; + self.0 = ordinal; + Ok(annot) + } +} + +/// A backbone fragment ion series. +/// +/// The nine mzPAF backbone series differ only by their letter, so they are one +/// enum with one letter table rather than nine variants repeated across every +/// match in this module. +#[derive(Clone, Copy, Eq, PartialEq, Ord, PartialOrd, Hash, Debug)] +#[allow(non_camel_case_types)] +#[repr(u8)] +pub enum Series { + a = 1, + b, + c, + d, + v, + w, + x, + y, + z, +} + +impl Series { + /// Every series, in discriminant order. Parallel to [`Self::CHARS`]. + pub const ALL: [Self; 9] = [ + Self::a, + Self::b, + Self::c, + Self::d, + Self::v, + Self::w, + Self::x, + Self::y, + Self::z, + ]; + /// The mzPAF letters, in discriminant order. The single place the + /// letter↔discriminant pairing lives. + const CHARS: &'static [u8; 9] = b"abcdvwxyz"; + + /// The mzPAF letter for this series. + pub const fn as_char(self) -> char { + Self::CHARS[self as usize - 1] as char + } + + /// The series for an mzPAF letter, or `None` if it names no backbone series. + fn from_char(c: char) -> Option { + let idx = Self::CHARS.iter().position(|&b| b as char == c)?; + Some(Self::ALL[idx]) + } +} + +impl Display for Series { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.as_char()) + } } -#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy, Default)] +/// The logical series-and-payload view of an [`IonAnnot`]. +/// +/// This is a *view*: `IonAnnot` stores a packed word and reconstructs this on +/// demand. Constructing one directly does not allocate an annotation. +#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy)] #[allow(non_camel_case_types)] pub enum IonSeriesOrdinal { - a { - ordinal: u8, - }, - b { - ordinal: u8, - }, - c { - ordinal: u8, - }, - d { - ordinal: u8, - }, - v { - ordinal: u8, - }, - w { - ordinal: u8, - }, - x { - ordinal: u8, - }, - y { - ordinal: u8, - }, - z { + /// One of the nine backbone series, at position `ordinal` in the ladder. + backbone { + series: Series, ordinal: u8, }, + /// An unannotated peak. `ordinal` is a uniqueness counter, not a position. unknown { ordinal: u8, }, precursor, - - /// This variant should not be used directly ... its mainly added to satisfy trait constraints by TinyVec - #[default] - None, + /// An internal fragment spanning residues `start..=end`. + internal { + start: u8, + end: u8, + }, + /// A bare immonium ion for an uppercase residue code. + immonium { + residue: char, + }, } impl IonSeriesOrdinal { - pub fn try_new(series: char, ordinal: Option) -> Result { - let tmp = match (series, ordinal) { - ('a', Some(ordinal)) => Self::a { ordinal }, - ('b', Some(ordinal)) => Self::b { ordinal }, - ('c', Some(ordinal)) => Self::c { ordinal }, - ('d', Some(ordinal)) => Self::d { ordinal }, - ('v', Some(ordinal)) => Self::v { ordinal }, - ('w', Some(ordinal)) => Self::w { ordinal }, - ('x', Some(ordinal)) => Self::x { ordinal }, - ('y', Some(ordinal)) => Self::y { ordinal }, - ('z', Some(ordinal)) => Self::z { ordinal }, - ('?', Some(ordinal)) => Self::unknown { ordinal }, - ('p', None) => Self::precursor, - ('p', Some(ordinal)) => { - return Err(IonParsingError::OrdinalOutOfRange { - ordinal: ordinal as i32, - series: Some(series), - }); - } - _ => { - return Err(IonParsingError::UnsupportedFragmentType { - fragment_type: series, - }); + /// Split into the `kind` discriminant and its `payload`, the two fields + /// [`IonAnnot`] packs. + /// + /// This and [`Self::from_parts`] are the only place the numbering lives. + /// The nine backbone series share one arm, so their discriminants come + /// from [`Series`] itself and cannot drift out of step with their letters. + const fn to_parts(self) -> (u32, u32) { + match self { + Self::backbone { series, ordinal } => (series as u32, ordinal as u32), + // Discriminant 0, so the all-zero word -- `IonAnnot::default()` -- + // is the unknown ion `?0` rather than an undecodable value. + Self::unknown { ordinal } => (0, ordinal as u32), + Self::precursor => (10, 0), + Self::internal { start, end } => { + (11, (start as u32) | ((end as u32) << INTERNAL_POS_BITS)) } - }; - - Ok(tmp) + Self::immonium { residue } => (12, (residue as u8 - b'A') as u32), + } } - pub fn terminality(&self) -> IonSeriesTerminality { - match self { - IonSeriesOrdinal::a { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::b { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::c { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::d { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::v { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::w { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::x { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::y { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::z { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::unknown { ordinal: _ } => IonSeriesTerminality::None, - IonSeriesOrdinal::precursor => IonSeriesTerminality::None, - IonSeriesOrdinal::None => panic!("IonSeriesOrdinal::None should not be used directly"), + /// Inverse of [`Self::to_parts`]. Total by construction: `unknown` is the + /// catch-all discriminant, so a value this build does not recognise -- only + /// reachable from a corrupted word -- decodes as an unknown ion instead of + /// panicking on a path `Display` (and therefore `Serialize`) reaches. + const fn from_parts(kind: u32, payload: u32) -> Self { + let ordinal = payload as u8; + match kind { + 1..=9 => Self::backbone { + // `kind` is in range, so this is the inverse of `series as u32`. + series: Series::ALL[kind as usize - 1], + ordinal, + }, + 10 => Self::precursor, + 11 => Self::internal { + start: (payload & mask(INTERNAL_POS_BITS)) as u8, + end: ((payload >> INTERNAL_POS_BITS) & mask(INTERNAL_POS_BITS)) as u8, + }, + 12 => Self::immonium { + residue: (b'A' + (payload & mask(IMMONIUM_BITS)) as u8) as char, + }, + _ => Self::unknown { ordinal }, } } - pub fn try_get_ordinal(&self) -> Option { - match self { - IonSeriesOrdinal::a { ordinal } => Some(*ordinal), - IonSeriesOrdinal::b { ordinal } => Some(*ordinal), - IonSeriesOrdinal::c { ordinal } => Some(*ordinal), - IonSeriesOrdinal::d { ordinal } => Some(*ordinal), - IonSeriesOrdinal::v { ordinal } => Some(*ordinal), - IonSeriesOrdinal::w { ordinal } => Some(*ordinal), - IonSeriesOrdinal::x { ordinal } => Some(*ordinal), - IonSeriesOrdinal::y { ordinal } => Some(*ordinal), - IonSeriesOrdinal::z { ordinal } => Some(*ordinal), - IonSeriesOrdinal::unknown { .. } => None, - // ?1 does not mean its an ordinal, just a placeholder - IonSeriesOrdinal::precursor => None, - IonSeriesOrdinal::None => None, + /// Build the series view from an mzPAF series letter and its ordinal. + /// + /// `p` is the only letter that takes no ordinal; every other one requires + /// one. Internal fragments and immonium ions are spelled differently and + /// have their own constructors. + fn from_series_char(c: char, ordinal: Option) -> Result { + if c == 'p' { + return match ordinal { + None => Ok(Self::precursor), + Some(ordinal) => Err(IonParsingError::OrdinalOutOfRange { ordinal, series: c }), + }; } + let ordinal = ordinal.ok_or(IonParsingError::MissingOrdinal { series: c })?; + if c == '?' { + return Ok(Self::unknown { ordinal }); + } + Series::from_char(c) + .map(|series| Self::backbone { series, ordinal }) + .ok_or(IonParsingError::UnsupportedFragmentType { fragment_type: c }) } } impl Display for IonSeriesOrdinal { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { match self { - IonSeriesOrdinal::a { ordinal } => write!(f, "a{}", ordinal), - IonSeriesOrdinal::b { ordinal } => write!(f, "b{}", ordinal), - IonSeriesOrdinal::c { ordinal } => write!(f, "c{}", ordinal), - IonSeriesOrdinal::d { ordinal } => write!(f, "d{}", ordinal), - IonSeriesOrdinal::v { ordinal } => write!(f, "v{}", ordinal), - IonSeriesOrdinal::w { ordinal } => write!(f, "w{}", ordinal), - IonSeriesOrdinal::x { ordinal } => write!(f, "x{}", ordinal), - IonSeriesOrdinal::y { ordinal } => write!(f, "y{}", ordinal), - IonSeriesOrdinal::z { ordinal } => write!(f, "z{}", ordinal), - IonSeriesOrdinal::unknown { ordinal } => write!(f, "?{}", ordinal), - IonSeriesOrdinal::precursor => write!(f, "p"), - IonSeriesOrdinal::None => panic!("IonSeriesOrdinal::None should not be used directly"), + Self::backbone { series, ordinal } => write!(f, "{}{}", series.as_char(), ordinal), + Self::unknown { ordinal } => write!(f, "?{}", ordinal), + Self::precursor => write!(f, "p"), + Self::internal { start, end } => write!(f, "m{}:{}", start, end), + Self::immonium { residue } => write!(f, "I{}", residue), } } } +#[cfg(test)] +mod tests { + use super::*; -impl FromStr for IonSeriesOrdinal { - type Err = IonParsingError; - - fn from_str(s: &str) -> Result { - // "b12" split into "b" and "12" - let (series_chunk, ordinal_chunk) = s.split_at(1); - let series_id = series_chunk.chars().next(); - let series_ordinal = ordinal_chunk.parse::(); - - match (series_id, series_ordinal) { - (None, _) => Err(IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Empty string"), - }), - (Some(x), Ok(y)) => IonSeriesOrdinal::try_new(x, Some(y)), - (Some('p'), Err(_)) => Ok(IonSeriesOrdinal::precursor), - (Some(_), Err(err)) => Err(IonParsingError::ParsingError { - error: format!("{ordinal_chunk} -> {err:?}"), - context: Some("Unable to parse the ordinal number"), - }), - } + fn ion(s: &str) -> IonAnnot { + IonAnnot::try_from(s).unwrap_or_else(|e| panic!("{s:?} must parse: {e}")) } -} -impl TryFrom for IonSeriesOrdinal { - type Error = IonParsingError; + /// The whole point of the packed representation. + #[test] + fn packs_into_one_word() { + assert_eq!(size_of::(), 4); + // Paired with an intensity on the scoring hot path; padding here would + // grow the inline TinyVec storage in `ExpectedIntensities`. + assert_eq!(size_of::<(IonAnnot, f32)>(), 8); + } + + /// `IonAnnot: Default` is not optional -- `tinyvec::Array` requires + /// `Item: Default`, `TimsElutionGroup` stores labels in a `TinyVec`, and + /// timsquery's `KeyLike` propagates the bound. So a default can reach any + /// serde path, and the zero word has to mean something. + /// + /// It means `?0`: charge is stored as `zigzag(charge - 1)`, so the zero + /// field is charge 1 rather than the impossible charge 0. The default is + /// therefore a real annotation and `Serialize`/`Deserialize` are inverses + /// on it, instead of rendering a value no constructor accepts. + #[test] + fn the_default_annotation_is_a_real_annotation() { + let d = IonAnnot::default(); + assert_eq!(d.to_string(), "?0"); + assert_eq!(d.get_charge(), 1); + assert_eq!(serde_json::to_string(&d).unwrap(), "\"?0\""); + assert_eq!(IonAnnot::try_from("?0").expect("the default parses"), d); + } - fn try_from(value: FragmentType) -> Result { - fn try_convert_ordinal(ordinal: usize, series: char) -> Result { - ordinal - .try_into() - .map_err(|_| IonParsingError::OrdinalOutOfRange { - ordinal: ordinal as i32, - series: Some(series), - }) + #[test] + fn test_deserialize() { + let cases = [ + ("b12", 'b', Some(12u8), 1i8, 0i8), + ("b12^3", 'b', Some(12), 3, 0), + ("y12^3", 'y', Some(12), 3, 0), + ("b12+i^3", 'b', Some(12), 3, 1), + ("b12+3i^3", 'b', Some(12), 3, 3), + ("b13", 'b', Some(13), 1, 0), + ("p^2", 'p', None, 2, 0), + ("p", 'p', None, 1, 0), + ("?12^2", '?', Some(12), 2, 0), + ]; + + for (input, series, ordinal, charge, isotope) in cases { + let annot = ion(input); + let expected = IonAnnot::try_new(series, ordinal, charge, isotope).unwrap(); + assert_eq!(annot, expected, "{input}"); + assert_eq!(annot.get_charge(), charge, "{input}"); + assert_eq!(annot.get_isotope(), isotope, "{input}"); + // Round-trips byte-identically when no loss is involved. + assert_eq!(format!("{}", annot), input); } - let tmp = match value { - FragmentType::a(ordinal, _) => IonSeriesOrdinal::a { - ordinal: try_convert_ordinal(ordinal.series_number, 'a')?, - }, - FragmentType::b(ordinal, _) => IonSeriesOrdinal::b { - ordinal: try_convert_ordinal(ordinal.series_number, 'b')?, - }, - FragmentType::c(ordinal, _) => IonSeriesOrdinal::c { - ordinal: try_convert_ordinal(ordinal.series_number, 'c')?, - }, - FragmentType::d(ordinal, _, _, _, _) => IonSeriesOrdinal::d { - ordinal: try_convert_ordinal(ordinal.series_number, 'd')?, - }, - FragmentType::v(ordinal, _, _, _) => IonSeriesOrdinal::v { - ordinal: try_convert_ordinal(ordinal.series_number, 'v')?, - }, - FragmentType::w(ordinal, _, _, _, _) => IonSeriesOrdinal::w { - ordinal: try_convert_ordinal(ordinal.series_number, 'w')?, - }, - FragmentType::x(ordinal, _) => IonSeriesOrdinal::x { - ordinal: try_convert_ordinal(ordinal.series_number, 'x')?, - }, - FragmentType::y(ordinal, _) => IonSeriesOrdinal::y { - ordinal: try_convert_ordinal(ordinal.series_number, 'y')?, - }, - FragmentType::z(ordinal, _) => IonSeriesOrdinal::z { - ordinal: try_convert_ordinal(ordinal.series_number, 'z')?, - }, - FragmentType::Precursor => IonSeriesOrdinal::precursor, - _ => { - return Err(IonParsingError::Custom { - error: format!("Unsupported fragment type: {value:?}"), - }); + } + + /// Every field must survive the pack/unpack at its extremes. + #[test] + fn every_field_round_trips_at_its_extremes() { + for charge in CHARGE_MIN..=CHARGE_MAX { + if charge == 0 { + continue; } - }; - Ok(tmp) + for isotope in ISOTOPE_MIN..=ISOTOPE_MAX { + for ordinal in [1u8, 2, 127, 255] { + for loss in [ + NeutralLoss::None, + NeutralLoss::Water, + NeutralLoss::PhosphoricAcidWater, + ] { + let a = + IonAnnot::try_new_with_loss('y', Some(ordinal), charge, isotope, loss) + .expect("in range"); + assert_eq!(a.get_charge(), charge); + assert_eq!(a.get_isotope(), isotope); + assert_eq!(a.try_get_ordinal(), Some(ordinal)); + assert_eq!(a.loss(), loss); + } + } + } + } } -} -#[cfg(test)] -mod tests { - use super::*; + #[test] + fn out_of_range_is_rejected_not_truncated() { + assert!(matches!( + IonAnnot::try_new('y', Some(1), CHARGE_MAX + 1, 0), + Err(IonParsingError::ChargeOutOfRange { .. }) + )); + assert!(matches!( + IonAnnot::try_new('y', Some(1), 1, ISOTOPE_MAX + 1), + Err(IonParsingError::IsotopeOutOfRange { .. }) + )); + assert!(IonAnnot::try_new('y', Some(1), 0, 0).is_err()); + assert!(matches!( + IonAnnot::try_new_internal(INTERNAL_POS_MAX + 1, 1, 1, 0, NeutralLoss::None), + Err(IonParsingError::OrdinalOutOfRange { .. }) + )); + } #[test] - fn test_ion_series_ord_from_str() { - let ion: IonSeriesOrdinal = IonSeriesOrdinal::from_str("b12").unwrap(); - assert_eq!(ion, IonSeriesOrdinal::b { ordinal: 12 }); + fn isotope_offset_respects_the_field_bound() { + let a = ion("y5"); + assert_eq!(a.try_with_offset_neutrons(2).unwrap().get_isotope(), 2); + assert!(a.try_with_offset_neutrons(ISOTOPE_MAX + 1).is_err()); } #[test] - fn test_deserialize() { - let serde_pairs = vec![ - ( - "b12", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 1, - isotope: 0, - }, - ), - ( - "b12^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 0, - }, - ), - ( - "y12^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::y { ordinal: 12 }, - charge: 3, - isotope: 0, - }, - ), - ( - "b12+i^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 1, - }, - ), - ( - "b12+3i^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 3, - }, - ), - ( - "b13", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 13 }, - charge: 1, - isotope: 0, - }, - ), - ( - "p^2", - IonAnnot { - series_ordinal: IonSeriesOrdinal::precursor, - charge: 2, - isotope: 0, - }, - ), - ( - "p", - IonAnnot { - series_ordinal: IonSeriesOrdinal::precursor, - charge: 1, - isotope: 0, - }, - ), - ( - "?12^2", - IonAnnot { - series_ordinal: IonSeriesOrdinal::unknown { ordinal: 12 }, - charge: 2, - isotope: 0, - }, - ), - ]; + fn parses_neutral_losses() { + let a = ion("y5-H2O"); + assert_eq!(a.loss(), NeutralLoss::Water); + assert_eq!(a.try_get_ordinal(), Some(5)); + assert_eq!(format!("{}", a), "y5-H2O"); + + // Non-canonical spelling resolves to the same annotation, and renders + // canonically -- so this pair is equal, which is the property that + // upholds per-precursor label uniqueness. + assert_eq!(ion("y5-CH3SOH"), ion("y5-CH4OS")); + assert_eq!(format!("{}", ion("y5-CH3SOH")), "y5-CH4OS"); + + // Loss combines with charge and isotope. + let b = ion("y10-NH3+i^2"); + assert_eq!(b.loss(), NeutralLoss::Ammonia); + assert_eq!(b.get_isotope(), 1); + assert_eq!(b.get_charge(), 2); + } + + #[test] + fn parses_internal_fragments() { + let a = ion("m2:11"); + assert_eq!( + a.series_ordinal(), + IonSeriesOrdinal::internal { start: 2, end: 11 } + ); + // An internal fragment has no ladder position. + assert_eq!(a.try_get_ordinal(), None); + assert_eq!(format!("{}", a), "m2:11"); - for (input, expected) in serde_pairs { - let annot = IonAnnot::try_from(input).unwrap(); - assert_eq!(annot, expected); + let b = ion("m11:12-CO"); + assert_eq!(b.loss(), NeutralLoss::CarbonMonoxide); + assert_eq!(format!("{}", b), "m11:12-CO"); + } - // Re-serialize and check that its the same - let serialized = format!("{}", annot); - assert_eq!(serialized, input); + #[test] + fn parses_bare_immonium_and_rejects_modified() { + let a = ion("IA"); + assert_eq!( + a.series_ordinal(), + IonSeriesOrdinal::immonium { residue: 'A' } + ); + assert_eq!(format!("{}", a), "IA"); + + // Carries an arbitrary mod string -- no field width represents it, so + // it must fail rather than silently degrade to a bare immonium. + assert!(matches!( + IonAnnot::try_from("IC[Carbamidomethyl]"), + Err(IonParsingError::UnsupportedModifiedImmonium { .. }) + )); + } + + #[test] + fn parses_and_applies_the_mass_error_suffix() { + let (rest, err) = split_mass_error("y1/-0.0005").unwrap(); + assert_eq!(ion(rest), ion("y1")); + assert_eq!(err, Some(MassError::Da(-0.0005))); + // theoretical = observed - error; verified against a real SpectraST + // peak: y1 for C-terminal R, observed 175.1184, theoretical 175.11895. + let theo = err.unwrap().theoretical_from_observed(175.1184); + assert!((theo - 175.1189).abs() < 1e-9, "got {theo}"); + + let (_, err) = split_mass_error("y6/1.2ppm").unwrap(); + assert_eq!(err, Some(MassError::Ppm(1.2))); + let theo = err.unwrap().theoretical_from_observed(700.0); + assert!((theo - 699.99916).abs() < 1e-4, "got {theo}"); + + // Absent suffix is not an error. + assert_eq!(split_mass_error("y6").unwrap(), ("y6", None)); + // TryFrom discards it rather than failing. + assert_eq!(IonAnnot::try_from("y1/-0.0005").unwrap(), ion("y1")); + } + + /// One representative of every `IonSeriesOrdinal` case: all nine backbone + /// series (each with a distinct ordinal, so a transposition in + /// `to_parts`/`from_parts` cannot cancel out) plus the four others, at + /// their field boundaries. + fn all_series() -> Vec { + let mut out: Vec = Series::ALL + .iter() + .enumerate() + .map(|(i, &series)| IonSeriesOrdinal::backbone { + series, + ordinal: i as u8 + 1, + }) + .collect(); + out.extend([ + IonSeriesOrdinal::precursor, + IonSeriesOrdinal::unknown { ordinal: 10 }, + IonSeriesOrdinal::internal { start: 2, end: 11 }, + // Both endpoints at their 6-bit ceiling: the widest payload the + // 12 bits hold, and the case a narrowed field would silently clip. + IonSeriesOrdinal::internal { + start: INTERNAL_POS_MAX, + end: INTERNAL_POS_MAX, + }, + IonSeriesOrdinal::immonium { residue: 'A' }, + // Highest residue index the 5-bit immonium field must hold. + IonSeriesOrdinal::immonium { residue: 'Z' }, + ]); + out + } + + /// `to_parts` and `from_parts` are hand-written inverses. Without this, + /// swapping two arms (`v` encoding as `w`) mislabels a whole ion series + /// and every other test still passes. + #[test] + fn every_series_variant_round_trips_through_the_packed_word() { + for series in all_series() { + let (kind, payload) = series.to_parts(); + assert_eq!(IonSeriesOrdinal::from_parts(kind, payload), series); + assert!(kind <= mask(KIND_BITS), "{series:?} kind overflows"); + assert!( + payload <= mask(PAYLOAD_BITS), + "{series:?} payload overflows" + ); + } + + // Every backbone series must land on its own discriminant; the other + // four are singletons and are covered by the round trip above. + let mut kinds: Vec = Series::ALL + .iter() + .map(|&series| { + IonSeriesOrdinal::backbone { series, ordinal: 1 } + .to_parts() + .0 + }) + .collect(); + kinds.sort_unstable(); + kinds.dedup(); + assert_eq!(kinds.len(), Series::ALL.len(), "two series share a kind"); + } + + /// `from_parts` claims to be total, and three of the sixteen `kind` + /// values are unassigned. Renumbering the assigned ones is safe only + /// while the unassigned ones stay decodable, since a corrupted word + /// reaches this on a `Display` (and so `Serialize`) path. + #[test] + fn every_kind_bit_pattern_decodes_without_panicking() { + for kind in 0..=mask(KIND_BITS) { + for payload in [0, 1, mask(PAYLOAD_BITS)] { + // Also exercises `Display`, which is where a partial decode + // would have panicked. + let _ = IonSeriesOrdinal::from_parts(kind, payload).to_string(); + } + } + } + + /// `Display`, `from_series_char` and the parser must agree with the + /// packing for every case. + #[test] + fn every_series_variant_round_trips_through_its_mzpaf_spelling() { + for series in all_series() { + let annot = IonAnnot::pack(series, NeutralLoss::None, 1, 0).expect("valid"); + assert_eq!(annot.series_ordinal(), series); + let text = annot.to_string(); + assert_eq!(ion(&text).series_ordinal(), series, "{text}"); + } + } + + /// The letters are the mzPAF spelling of the discriminants, and + /// `from_char`/`as_char` are hand-written inverses of each other. + #[test] + fn series_letters_and_discriminants_agree() { + assert_eq!(Series::ALL.len(), Series::CHARS.len()); + for &series in &Series::ALL { + assert_eq!(Series::from_char(series.as_char()), Some(series)); + } + for c in ['p', '?', 'm', 'I', 'q', 'A'] { + assert_eq!(Series::from_char(c), None, "{c} is not a backbone series"); } } + + /// An unrepresentable loss must fail loudly. Parsing `y1-HCOOH` as plain + /// `y1` would put a loss peak's m/z on the `y1` label and collide with the + /// real `y1`. + #[test] + fn unrepresentable_loss_is_rejected_not_stripped() { + assert!(matches!( + IonAnnot::try_from("y1-HCOOH"), + Err(IonParsingError::UnsupportedNeutralLoss { .. }) + )); + } } diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs new file mode 100644 index 00000000..2fcea755 --- /dev/null +++ b/rust/micromzpaf/src/loss.rs @@ -0,0 +1,449 @@ +//! Neutral losses, keyed by atomic composition rather than by spelling. +//! +//! Libraries write the same chemical loss different ways. Two real examples +//! from the HUPO-PSI mzSpecLib corpus: +//! +//! | written | library | composition | +//! |---|---|---| +//! | `-CH3SOH` | NIST | C1H4O1S1 | +//! | `-CH4OS` | SpectraST | C1H4O1S1 | +//! | `-NH2-CO-CH2SH` | NIST | C2H5N1O1S1 | +//! | `-C2H5NOS` | SpectraST | C2H5N1O1S1 | +//! +//! Keying on the string would make `y5-CH4OS` and `y5-CH3SOH` distinct labels +//! for one ion. Since fragment labels must be unique within a precursor (see +//! `ExpectedIntensities::try_from_pairs` in timsseek), that is exactly the +//! wrong direction: it hides a genuine duplicate behind two spellings. +//! +//! So parsing goes `text -> composition -> discriminant`, and the composition +//! is a *parse-time* concept only. What gets stored on an `IonAnnot` is the +//! discriminant, so there is no per-annotation cost at runtime. + +use std::fmt::Display; + +use crate::IonParsingError; + +/// Slots in [`Composition`], in the order [`Composition::new`] takes them. +const C: usize = 0; +const H: usize = 1; +const N: usize = 2; +const O: usize = 3; +const S: usize = 4; +const P: usize = 5; + +/// Atom counts for the elements that appear in peptide neutral losses. +/// +/// Deliberately not a general chemical formula: these losses only ever draw +/// from C/H/N/O/S/P, and keeping it to six `u8`s makes equality a single +/// 6-byte compare during the parse-time table lookup. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub(crate) struct Composition([u8; 6]); + +impl Composition { + pub(crate) const fn new(c: u8, h: u8, n: u8, o: u8, s: u8, p: u8) -> Self { + Self([c, h, n, o, s, p]) + } + + /// Multiply every count, saturating. Used for the `2H2O` multiplier form. + fn scaled(self, k: u8) -> Self { + Self(self.0.map(|n| n.saturating_mul(k))) + } + + fn plus(self, other: Self) -> Self { + Self(std::array::from_fn(|i| { + self.0[i].saturating_add(other.0[i]) + })) + } + + /// Parse a bare formula like `H2O`, `CH4OS`, `C2H5NOS`. + /// + /// Only single-letter C/H/N/O/S/P are recognized; anything else is an + /// error rather than a silent skip, so an unsupported loss surfaces as + /// "not representable" instead of being mistaken for a smaller one. + fn parse_formula(s: &str) -> Result { + if s.is_empty() { + return Err(IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Empty neutral-loss formula"), + }); + } + let mut out = Composition::default(); + let b = s.as_bytes(); + let mut i = 0; + while i < b.len() { + let elem = b[i]; + i += 1; + let start = i; + while i < b.len() && b[i].is_ascii_digit() { + i += 1; + } + let count: u8 = if start == i { + 1 + } else { + s[start..i] + .parse() + .map_err(|_| IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Neutral-loss atom count out of range"), + })? + }; + let slot = match elem { + b'C' => C, + b'H' => H, + b'N' => N, + b'O' => O, + b'S' => S, + b'P' => P, + _ => { + return Err(IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Unsupported element in neutral loss"), + }); + } + }; + out.0[slot] = out.0[slot].saturating_add(count); + } + Ok(out) + } + + /// Parse a full loss expression: `-` separated terms, each optionally + /// prefixed by a repeat count. `2H2O`, `H2O-NH3`, `NH2-CO-CH2SH`. + /// + /// Because terms are summed, ordering and multiplier spelling collapse for + /// free: `H2O-NH3` == `NH3-H2O`, and `2H2O` == `H2O-H2O`. + pub(crate) fn parse_expression(s: &str) -> Result { + let mut total = Composition::default(); + for term in s.split('-') { + let term = term.trim(); + if term.is_empty() { + return Err(IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Empty term in neutral-loss expression"), + }); + } + // Leading digits are a repeat count for the whole term. + let digits = term + .find(|c: char| !c.is_ascii_digit()) + .unwrap_or(term.len()); + let (mult, formula) = if digits > 0 { + let m: u8 = term[..digits] + .parse() + .map_err(|_| IonParsingError::ParsingError { + error: s.to_string(), + context: Some("Neutral-loss repeat count out of range"), + })?; + (m, &term[digits..]) + } else { + (1, term) + }; + total = total.plus(Composition::parse_formula(formula)?.scaled(mult)); + } + Ok(total) + } +} + +/// The neutral losses this crate can represent, as a packed discriminant. +/// +/// Scoped deliberately: DIA-NN emits none, Spectronaut two (`-H2O`, `-NH3`), +/// NIST eight, plus the phospho losses that no non-phospho corpus can show. +/// SpectraST's wider combinatorics are NOT here -- an unlisted loss parses to a +/// composition that misses the table and is reported as unrepresentable, which +/// routes the peak to an unknown label rather than silently mislabelling it. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] +#[repr(u8)] +pub enum NeutralLoss { + #[default] + None = 0, + /// H2O + Water = 1, + /// NH3 + Ammonia = 2, + /// CO + CarbonMonoxide = 3, + /// CO2 + CarbonDioxide = 4, + /// 2 H2O + WaterX2 = 5, + /// 2 NH3 + AmmoniaX2 = 6, + /// H2O + NH3 + WaterAmmonia = 7, + /// CH4OS -- methanesulfenic acid, off oxidized Met. Also spelled `CH3SOH`. + Methanesulfenic = 8, + /// C2H5NOS -- also spelled `NH2-CO-CH2SH`. + Carbamidomethylthiol = 9, + /// H3PO4 -- phospho-Ser/Thr. + PhosphoricAcid = 10, + /// HPO3 -- phospho-Tyr, and phospho-Ser/Thr. + Metaphosphoric = 11, + /// H3PO4 + H2O + PhosphoricAcidWater = 12, +} + +/// `(composition, discriminant, canonical spelling)`. +/// +/// The canonical spelling is what `Display` emits, so a non-canonical input +/// (`-CH3SOH`) round-trips to the canonical form (`-CH4OS`). Round-trip tests +/// must therefore compare parsed values, not bytes. +const TABLE: &[(Composition, NeutralLoss, &str)] = &[ + ( + Composition::new(0, 2, 0, 1, 0, 0), + NeutralLoss::Water, + "H2O", + ), + ( + Composition::new(0, 3, 1, 0, 0, 0), + NeutralLoss::Ammonia, + "NH3", + ), + ( + Composition::new(1, 0, 0, 1, 0, 0), + NeutralLoss::CarbonMonoxide, + "CO", + ), + ( + Composition::new(1, 0, 0, 2, 0, 0), + NeutralLoss::CarbonDioxide, + "CO2", + ), + ( + Composition::new(0, 4, 0, 2, 0, 0), + NeutralLoss::WaterX2, + "2H2O", + ), + ( + Composition::new(0, 6, 2, 0, 0, 0), + NeutralLoss::AmmoniaX2, + "2NH3", + ), + ( + Composition::new(0, 5, 1, 1, 0, 0), + NeutralLoss::WaterAmmonia, + "H2O-NH3", + ), + ( + Composition::new(1, 4, 0, 1, 1, 0), + NeutralLoss::Methanesulfenic, + "CH4OS", + ), + ( + Composition::new(2, 5, 1, 1, 1, 0), + NeutralLoss::Carbamidomethylthiol, + "C2H5NOS", + ), + ( + Composition::new(0, 3, 0, 4, 0, 1), + NeutralLoss::PhosphoricAcid, + "H3PO4", + ), + ( + Composition::new(0, 1, 0, 3, 0, 1), + NeutralLoss::Metaphosphoric, + "HPO3", + ), + ( + Composition::new(0, 5, 0, 5, 0, 1), + NeutralLoss::PhosphoricAcidWater, + "H3PO4-H2O", + ), +]; + +impl NeutralLoss { + /// Inverse of the `#[repr(u8)]` discriminant, for unpacking out of a bit + /// field. Lives next to the enum so the two cannot drift apart. + /// + /// An unrecognized value maps to [`Self::None`]: the only way to produce + /// one is a reserved discriminant, which no constructor emits. + pub(crate) fn from_discriminant(d: u8) -> Self { + match d { + 1 => Self::Water, + 2 => Self::Ammonia, + 3 => Self::CarbonMonoxide, + 4 => Self::CarbonDioxide, + 5 => Self::WaterX2, + 6 => Self::AmmoniaX2, + 7 => Self::WaterAmmonia, + 8 => Self::Methanesulfenic, + 9 => Self::Carbamidomethylthiol, + 10 => Self::PhosphoricAcid, + 11 => Self::Metaphosphoric, + 12 => Self::PhosphoricAcidWater, + _ => Self::None, + } + } + + /// Resolve a loss expression (without the leading `-`) to a discriminant. + /// + /// `None` means "parsed as a valid composition, but not one we represent" -- + /// distinct from `Err`, which means the text was not a loss expression at + /// all. Callers route the former to an unknown label and the latter to a + /// parse failure. + pub(crate) fn from_expression(s: &str) -> Result, IonParsingError> { + let comp = Composition::parse_expression(s)?; + Ok(TABLE + .iter() + .find(|(c, _, _)| *c == comp) + .map(|(_, l, _)| *l)) + } + + /// Canonical spelling, without the leading `-`. Empty for [`Self::None`]. + pub(crate) fn canonical(self) -> &'static str { + if self == NeutralLoss::None { + return ""; + } + TABLE + .iter() + .find(|(_, l, _)| *l == self) + .map(|(_, _, s)| *s) + .unwrap_or("") + } +} + +impl Display for NeutralLoss { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + if *self == NeutralLoss::None { + return Ok(()); + } + write!(f, "-{}", self.canonical()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The count parser: implicit 1, multi-digit counts, and a count past the + /// `u8` slot. Single-digit cases are covered by every TABLE row. + #[test] + fn atom_counts_parse_and_reject_out_of_range() { + assert_eq!( + Composition::parse_expression("H2O").unwrap(), + Composition::new(0, 2, 0, 1, 0, 0), + "an element with no digits is one atom" + ); + assert_eq!( + Composition::parse_expression("C10H12").unwrap(), + Composition::new(10, 12, 0, 0, 0, 0), + "counts are multi-digit, not one digit per element" + ); + assert_eq!( + Composition::parse_expression("C255").unwrap(), + Composition::new(255, 0, 0, 0, 0, 0), + "the u8 slot is full at 255" + ); + assert!( + Composition::parse_expression("C256").is_err(), + "one past the slot must be an error, not a wrap to 0" + ); + } + + /// The two cross-library spelling collisions this module exists for. + #[test] + fn different_spellings_resolve_to_one_loss() { + // NIST vs SpectraST, methanesulfenic acid. + assert_eq!( + NeutralLoss::from_expression("CH3SOH").unwrap(), + Some(NeutralLoss::Methanesulfenic) + ); + assert_eq!( + NeutralLoss::from_expression("CH4OS").unwrap(), + Some(NeutralLoss::Methanesulfenic) + ); + // NIST structural notation vs SpectraST molecular notation. + assert_eq!( + NeutralLoss::from_expression("NH2-CO-CH2SH").unwrap(), + Some(NeutralLoss::Carbamidomethylthiol) + ); + assert_eq!( + NeutralLoss::from_expression("C2H5NOS").unwrap(), + Some(NeutralLoss::Carbamidomethylthiol) + ); + } + + /// Summing terms collapses ordering and multiplier spelling for free. + #[test] + fn ordering_and_multipliers_normalize() { + assert_eq!( + NeutralLoss::from_expression("H2O-NH3").unwrap(), + NeutralLoss::from_expression("NH3-H2O").unwrap() + ); + assert_eq!( + NeutralLoss::from_expression("2H2O").unwrap(), + NeutralLoss::from_expression("H2O-H2O").unwrap() + ); + assert_eq!( + NeutralLoss::from_expression("2H2O").unwrap(), + Some(NeutralLoss::WaterX2) + ); + } + + /// A well-formed composition outside the table is `Ok(None)` -- "valid but + /// not representable" -- while malformed text is `Err`. Callers need to + /// tell those apart to route one to an unknown label and the other to a + /// parse failure. + #[test] + fn unrepresentable_is_distinct_from_malformed() { + assert_eq!(NeutralLoss::from_expression("HCOOH").unwrap(), None); + assert!(NeutralLoss::from_expression("Xe2").is_err()); + assert!(NeutralLoss::from_expression("").is_err()); + assert!(NeutralLoss::from_expression("H2O-").is_err()); + } + + /// Every table entry must survive canonical -> composition -> discriminant, + /// and back out through the bit field. `from_discriminant` hand-mirrors the + /// `#[repr(u8)]` values, so adding a loss without updating it would + /// silently decode as [`NeutralLoss::None`] -- this is what catches that. + #[test] + fn table_round_trips_through_canonical_spelling() { + for (_comp, loss, canon) in TABLE { + assert_eq!( + NeutralLoss::from_expression(canon).unwrap(), + Some(*loss), + "canonical spelling {canon} must resolve to its own loss" + ); + assert_eq!(loss.canonical(), *canon); + assert_eq!( + NeutralLoss::from_discriminant(*loss as u8), + *loss, + "{canon} does not survive the discriminant round trip" + ); + } + assert_eq!(NeutralLoss::from_discriminant(0), NeutralLoss::None); + // Every non-None variant must be in TABLE, or it has no spelling and no + // composition and could never be produced by parsing. + assert_eq!( + TABLE.len(), + (1..=u8::MAX) + .filter(|d| NeutralLoss::from_discriminant(*d) != NeutralLoss::None) + .count(), + "a variant is decodable but missing from TABLE" + ); + } + + /// Compositions must be unique: two entries sharing one would make the + /// table lookup order-dependent. + #[test] + fn table_compositions_are_unique() { + for (i, (a, _, sa)) in TABLE.iter().enumerate() { + for (b, _, sb) in TABLE.iter().skip(i + 1) { + assert_ne!(a, b, "{sa} and {sb} share a composition"); + } + // `scaled` and `plus` saturate, so an absurd input like `200H2O` + // lands on 255 in some slot. That is only safe to report as + // unrepresentable while no table entry holds a saturated count -- + // otherwise the saturation would alias onto a real loss. + assert!( + a.0.iter().all(|&n| n < u8::MAX), + "{sa} holds a saturated atom count" + ); + } + } + + /// The two things the TABLE round trip cannot see: the `-` prefix Display + /// adds, and that `None` renders as nothing at all rather than "-". + #[test] + fn display_adds_the_prefix_and_none_renders_empty() { + assert_eq!(NeutralLoss::None.to_string(), ""); + assert_eq!(NeutralLoss::Water.to_string(), "-H2O"); + } +} diff --git a/rust/timsquery/src/lib.rs b/rust/timsquery/src/lib.rs index e0f2c8d9..84fd7bf8 100644 --- a/rust/timsquery/src/lib.rs +++ b/rust/timsquery/src/lib.rs @@ -58,7 +58,8 @@ pub mod ion { IonAnnot, IonParsingError, IonSeriesOrdinal, - IonSeriesTerminality, + Series, + UnknownIonCounter, }; } diff --git a/rust/timsquery/src/serde/diann_io.rs b/rust/timsquery/src/serde/diann_io.rs index 9d18d66b..2392a41b 100644 --- a/rust/timsquery/src/serde/diann_io.rs +++ b/rust/timsquery/src/serde/diann_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use crate::models::OwnedSourceId; use arrow::array::{ @@ -413,7 +414,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(rows.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::new(); for (i, row) in rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -433,8 +434,7 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions.next(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -743,7 +743,7 @@ fn parse_precursor_group_from_parquet( let mut fragment_mzs = Vec::with_capacity(indices.len()); buffers.fragment_labels.clear(); let mut rel_intensities = Vec::with_capacity(indices.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::new(); for (i, &idx) in indices.iter().enumerate() { let fragment_mz = columns.product_mzs[idx] as f64; @@ -765,8 +765,7 @@ fn parse_precursor_group_from_parquet( columns.fragment_loss_types[idx], i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions.next(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); rel_intensities.push((ion_annot, rel_intensity)); diff --git a/rust/timsquery/src/serde/skyline_io.rs b/rust/timsquery/src/serde/skyline_io.rs index f91c4203..a18da490 100644 --- a/rust/timsquery/src/serde/skyline_io.rs +++ b/rust/timsquery/src/serde/skyline_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use serde::{ Deserialize, @@ -309,7 +310,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(fragment_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(fragment_rows.len()); - let mut num_unknown_losses = 0u8; + let mut unknown_ions = UnknownIonCounter::new(); for (i, row) in fragment_rows.iter().enumerate() { let fragment_mz = row.product_mz; @@ -348,8 +349,7 @@ fn parse_precursor_group( falling back to unknown ion", frag_type, i ); - num_unknown_losses = num_unknown_losses.saturating_add(1); - IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)? + unknown_ions.next(frag_charge as i8)? } }; diff --git a/rust/timsquery/src/serde/spectronaut_io.rs b/rust/timsquery/src/serde/spectronaut_io.rs index c850ce94..2c1a6680 100644 --- a/rust/timsquery/src/serde/spectronaut_io.rs +++ b/rust/timsquery/src/serde/spectronaut_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use serde::Deserialize; use std::path::Path; @@ -283,7 +284,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(included_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(included_rows.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::new(); for (i, row) in included_rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -303,8 +304,7 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions.next(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); diff --git a/rust/timsquery/src/traits/fragment_label.rs b/rust/timsquery/src/traits/fragment_label.rs index e415a698..335e08be 100644 --- a/rust/timsquery/src/traits/fragment_label.rs +++ b/rust/timsquery/src/traits/fragment_label.rs @@ -1,5 +1,6 @@ use crate::IonAnnot; use crate::traits::KeyLike; +use micromzpaf::Series; use std::sync::Arc; /// Capability to apply a decoy m/z shift to a fragment label's m/z. @@ -20,11 +21,12 @@ pub trait DecoyShift { /// `FragmentLabel`. NOTE: decoy m/z shifting is gated on `DecoyShift`, not this /// trait -- a decoy variant computes its shift through `DecoyShift::decoy_shift_mz`, /// which every label implements. -// TODO(fragment-features): add series() accessor once a standalone series -// enum exists (today series and ordinal are fused in `IonSeriesOrdinal`), and -// wire this bound into the score that consumes it. +// TODO(fragment-features): wire this bound into the score that consumes it. pub trait FragmentLabel: KeyLike + DecoyShift { fn try_get_ordinal(&self) -> Option; + /// The backbone series, or `None` for ions that sit on no ladder + /// (precursor, unknown, internal, immonium). + fn try_get_series(&self) -> Option; fn get_charge(&self) -> i8; } @@ -45,6 +47,10 @@ impl FragmentLabel for IonAnnot { IonAnnot::try_get_ordinal(self) } + fn try_get_series(&self) -> Option { + IonAnnot::try_get_series(self) + } + fn get_charge(&self) -> i8 { IonAnnot::get_charge(self) } diff --git a/rust/timsseek/src/lib.rs b/rust/timsseek/src/lib.rs index 2ea374f6..88908d39 100644 --- a/rust/timsseek/src/lib.rs +++ b/rust/timsseek/src/lib.rs @@ -33,6 +33,5 @@ pub use scoring::{ pub use timsquery::ion::{ IonAnnot, IonParsingError, - IonSeriesTerminality, }; pub use traits::ScorerQueriable; From d92399c6d58d42be740146c2f7b60e368aa020f1 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 10:56:46 -0700 Subject: [PATCH 2/8] fix(micromzpaf): make the layout self-checking and the docs true Review follow-ups, all mechanical. The crate doc led with a false claim: the predecessor was already 4 bytes, so `size_of` pins nothing on its own. The argument is against the alternative -- bolting a loss byte onto the old struct would have made a 5-byte key that pads to a 12-byte tuple, growing timsseek's inline storage 104 -> 156. The test now asserts the added fields cost no space rather than restating a size that never changed. Isotope offsets become unsigned. `ISOTOPE_MIN = -7` let `Display` emit `y5+-2i`, which is not mzPAF and reparsed happily -- so a library could round-trip through serde into text no other reader accepts, contradicting the crate's own doc. Rejecting negatives makes that spelling unreachable and frees the sign bit, widening the useful range from 0..=7 to 0..=15. `const` assertions replace the prose contract: total width, abutting shifts, two internal endpoints per payload, the immonium alphabet, both zigzag bounds, and that the loss table has not outgrown `LOSS_BITS` (`pack` masks the discriminant, so a 64th loss would decode as a different one). `TABLE` becomes the index for `from_discriminant` and `canonical`, deleting a 13-arm hand-mirrored match and turning two linear scans into `TABLE.get`. The drift the old round-trip test guarded against is now one invariant. `IonAnnot` loses `Ord`. The packed word sorts by ordinal, then loss, then isotope, then charge, then series, which is not an order anyone means; `KeyLike` never required it, and the three downstream `PartialOrd` derives had no call sites. `UnknownIonCounter` loses `Copy`/`Clone`: a forked uniqueness counter reissues labels, which is the bug it exists to prevent. `ElutionGroupInput::try_fill_labels_annot` mints through that counter instead of casting an index, so running past the label space errors rather than wrapping into duplicates; it also no longer stamps isotope 1 onto every placeholder. `try_fill_labels_u8` rejects >255 fragments for the same reason. `ParsingError.context` was an `Option` that was never `None`, and two sites reported the wrong thing: a bare `I` as a *modified* immonium, and `p1` as an out-of-range ordinal for a series that takes none. Fixes the two `cargo doc` failures: public docs linked private items. --- rust/micromzpaf/src/lib.rs | 300 +++++++++++++----- rust/micromzpaf/src/loss.rs | 116 ++++--- rust/timsquery/src/serde/diann_io.rs | 55 ++-- rust/timsquery/src/serde/diann_speclib_io.rs | 4 +- .../src/serde/elution_group_inputs.rs | 63 +++- rust/timsquery/src/serde/skyline_io.rs | 6 +- rust/timsquery/src/serde/spectronaut_io.rs | 6 +- 7 files changed, 351 insertions(+), 199 deletions(-) diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 9301ff4a..55f0e7dd 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -1,21 +1,31 @@ //! Compact representation of fragment ion annotations for mass spectrometry. //! //! A spectral library carries one annotation per fragment, so this type is -//! replicated millions of times in a loaded arena and compared on the scoring -//! hot path. It is therefore a packed `u32` rather than a struct of fields: -//! equality is a single word compare, and `(IonAnnot, f32)` stays 8 bytes with -//! no padding. +//! replicated millions of times in a loaded arena. It is a packed `u32` so that +//! the annotations an mzSpecLib-shaped library needs -- neutral losses, +//! internal fragments, immonium ions -- fit *without* growing the type. +//! +//! That is the whole argument, and it is about the alternative rather than +//! about the previous representation. The predecessor was a 4-byte struct of +//! `(series+ordinal, charge, isotope)`; bolting a loss field onto it would have +//! cost three bytes, not one, because a 5-byte key paired with an `f32` pads to +//! a 12-byte tuple. That would have grown the inline `TinyVec` storage in +//! timsseek's `ExpectedIntensities` from 104 to 156 bytes. Packed, the tuple +//! stays 8 bytes and the loss rides along in bits nobody was using. //! //! # Bit layout //! //! ```text -//! bit: 31 30 29 18 17 12 11 8 7 4 3 0 +//! bit: 31 30 29 18 17 12 11 8 7 4 3 0 //! ┌───────┬────────────────┬──────────┬────────┬────────┬───────┐ //! │ spare │ payload │ loss │isotope │ charge │ kind │ -//! │ 2b │ 12b │ 6b │ 4b zz │ 4b zz │ 4b │ +//! │ 2b │ 12b │ 6b │ 4b │ 4b zz │ 4b │ //! └───────┴────────────────┴──────────┴────────┴────────┴───────┘ //! ``` //! +//! The widths and their shifts are checked by `const` assertions next to the +//! constants, so this diagram cannot drift away from the code. +//! //! `payload` is reinterpreted per `kind` -- a tagged union inside the word: //! //! | kind | payload | @@ -30,11 +40,12 @@ //! because `IonAnnot: Default` is forced by `tinyvec::Array` and a default can //! reach any serde path. //! -//! `charge` and `isotope` are zigzag-encoded so they stay signed in 4 bits. -//! Their ranges (±7) are far wider than anything observed: the HUPO-PSI corpus -//! tops out at charge 3 and isotope 3, with no negative charges at all. Because -//! the field truncates rather than wrapping loudly, every constructor -//! range-checks -- see [`IonAnnot::try_new`]. +//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned, +//! because a negative offset has no spelling this crate can emit (see +//! [`ISOTOPE_MIN`]). Both ranges are far wider than anything observed: the +//! HUPO-PSI corpus tops out at charge 3 and isotope 3, with no negative charges +//! at all. Because a bit field truncates rather than wrapping loudly, every +//! constructor range-checks -- see [`IonAnnot::try_new`]. //! //! # mzPAF compliance //! @@ -46,7 +57,8 @@ //! Not supported: negative isotope offsets, modified immonium //! (`IC[Carbamidomethyl]` carries an arbitrary mod string), and losses outside //! the [`NeutralLoss`] table. These are reported as errors, never coerced into -//! a nearby representable ion. +//! a nearby representable ion -- and unrepresentable in the field too, so +//! `Display` cannot emit a spelling the parser would reject. //! //! # Examples //! @@ -87,17 +99,69 @@ const PAYLOAD_SHIFT: u32 = 18; const PAYLOAD_BITS: u32 = 12; /// Widest charge the 4-bit zigzag field holds. Observed maximum is 3. +/// +/// The bias applied when packing makes `-7..=8` representable, but the range is +/// capped symmetrically: nothing has ever needed charge 8, and an asymmetric +/// public bound invites the reader to check the arithmetic rather than trust it. pub const CHARGE_MIN: i8 = -7; pub const CHARGE_MAX: i8 = 7; -/// Widest isotope offset the 4-bit zigzag field holds. Observed maximum is 3. -pub const ISOTOPE_MIN: i8 = -7; -pub const ISOTOPE_MAX: i8 = 7; +/// Widest isotope offset the 4-bit field holds. Observed maximum is 3. +/// +/// Isotope offsets are unsigned. mzPAF spells a negative offset `-Ni`, which +/// this crate does not parse, so allowing one in the field would let `Display` +/// emit `y5+-2i` -- text no other mzPAF reader accepts, through a type whose +/// wire format *is* that text. Rejecting it here keeps the invalid spelling +/// unreachable and spends the whole 4 bits on offsets that can round-trip. +pub const ISOTOPE_MIN: i8 = 0; +pub const ISOTOPE_MAX: i8 = mask(ISOTOPE_BITS) as i8; /// Width of each internal-fragment endpoint inside `payload`. const INTERNAL_POS_BITS: u32 = 6; /// Width of the immonium residue index inside `payload`. const IMMONIUM_BITS: u32 = 5; /// Widest residue index an internal fragment endpoint holds. -pub(crate) const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; +pub const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; + +// The layout is otherwise enforced by prose and a diagram. These make it +// self-checking, so a field width cannot be widened without the build failing. +const _: () = assert!( + KIND_BITS + CHARGE_BITS + ISOTOPE_BITS + LOSS_BITS + PAYLOAD_BITS <= 32, + "the fields overflow the word" +); +const _: () = assert!(KIND_SHIFT == 0); +const _: () = assert!(CHARGE_SHIFT == KIND_SHIFT + KIND_BITS, "fields must abut"); +const _: () = assert!( + ISOTOPE_SHIFT == CHARGE_SHIFT + CHARGE_BITS, + "fields must abut" +); +const _: () = assert!( + LOSS_SHIFT == ISOTOPE_SHIFT + ISOTOPE_BITS, + "fields must abut" +); +const _: () = assert!(PAYLOAD_SHIFT == LOSS_SHIFT + LOSS_BITS, "fields must abut"); +const _: () = assert!( + 2 * INTERNAL_POS_BITS <= PAYLOAD_BITS, + "two internal-fragment endpoints must fit in one payload" +); +const _: () = assert!( + mask(IMMONIUM_BITS) >= (b'Z' - b'A') as u32, + "the immonium field must hold every uppercase residue" +); +// The zigzag bounds are hand-derived; these are what make them checked. +const _: () = assert!( + zigzag_charge(CHARGE_MIN) <= mask(CHARGE_BITS) + && zigzag_charge(CHARGE_MAX) <= mask(CHARGE_BITS), + "the charge range does not survive zigzag inside CHARGE_BITS" +); +const _: () = assert!( + ISOTOPE_MIN >= 0 && ISOTOPE_MAX as u32 <= mask(ISOTOPE_BITS), + "the isotope range does not fit ISOTOPE_BITS" +); +// `pack` masks the loss discriminant, so a table grown past the field would +// truncate into a *different* loss rather than fail. +const _: () = assert!( + NeutralLoss::COUNT as u32 <= mask(LOSS_BITS), + "the loss table has outgrown LOSS_BITS" +); #[inline] const fn mask(bits: u32) -> u32 { @@ -131,9 +195,13 @@ const fn unzigzag_charge(u: u32) -> i8 { /// Compact representation of fragment annotations. /// -/// A packed `u32`; see the crate docs for the bit layout. Ordering is by the -/// packed word, not field-by-field. -#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] +/// A packed `u32`; see the crate docs for the bit layout. +/// +/// Deliberately not `Ord`. Ordering the packed word sorts by ordinal first, +/// then loss, then isotope, then charge, then series -- an order nobody means. +/// Nothing in the workspace sorts annotations, and `KeyLike` does not require +/// it, so the trait is not offered rather than offered and meaningless. +#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, Default)] pub struct IonAnnot(u32); impl Serialize for IonAnnot { @@ -266,7 +334,7 @@ impl IonAnnot { Ok(IonAnnot( (kind << KIND_SHIFT) | ((zigzag_charge(charge) & mask(CHARGE_BITS)) << CHARGE_SHIFT) - | ((zigzag(isotope) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT) + | (((isotope as u32) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT) | ((loss as u32 & mask(LOSS_BITS)) << LOSS_SHIFT) | ((payload & mask(PAYLOAD_BITS)) << PAYLOAD_SHIFT), )) @@ -284,7 +352,7 @@ impl IonAnnot { #[inline] pub fn get_isotope(&self) -> i8 { - unzigzag((self.0 >> ISOTOPE_SHIFT) & mask(ISOTOPE_BITS)) + ((self.0 >> ISOTOPE_SHIFT) & mask(ISOTOPE_BITS)) as i8 } /// The neutral loss this ion carries, [`NeutralLoss::None`] if it carries @@ -310,7 +378,7 @@ impl IonAnnot { } Ok(IonAnnot( (self.0 & !(mask(ISOTOPE_BITS) << ISOTOPE_SHIFT)) - | ((zigzag(new_isotope) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT), + | (((new_isotope as u32) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT), )) } @@ -378,10 +446,9 @@ pub fn split_mass_error(s: &str) -> Result<(&str, Option), IonParsing Some(n) => (n, true), None => (tail, false), }; - let v: f64 = num.parse().map_err(|_| IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Unable to parse the mass-error suffix"), - })?; + let v: f64 = num + .parse() + .map_err(|_| IonParsingError::parse(s, "Unable to parse the mass-error suffix"))?; Ok(( head, Some(if is_ppm { @@ -402,12 +469,9 @@ impl IonAnnot { // charge: trailing ^N let (rest, charge) = match value.split_once('^') { Some((rest, charge)) => { - let charge = charge - .parse::() - .map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse the charge number"), - })?; + let charge = charge.parse::().map_err(|_| { + IonParsingError::parse(value, "Unable to parse the charge number") + })?; (rest, charge) } None => (value, 1), @@ -416,21 +480,16 @@ impl IonAnnot { // isotope: +Ni. Negative isotope offsets are not supported. let (rest, isotope) = match rest.split_once('+') { Some((rest, adducts)) => { - let adducts = adducts - .strip_suffix('i') - .ok_or(IonParsingError::ParsingError { - error: adducts.to_string(), - context: Some("Unsupported adduct found"), - })?; + let adducts = adducts.strip_suffix('i').ok_or(IonParsingError::parse( + value, + "Only the isotope adduct '+Ni' is supported", + ))?; let isotope = if adducts.is_empty() { 1 } else { - adducts - .parse::() - .map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse the isotope number"), - })? + adducts.parse::().map_err(|_| { + IonParsingError::parse(value, "Unable to parse the isotope number") + })? }; (rest, isotope) } @@ -456,13 +515,11 @@ impl IonAnnot { if let Some(spans) = core.strip_prefix('m') && let Some((a, b)) = spans.split_once(':') { - let start = a.parse::().map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse internal-fragment start"), + let start = a.parse::().map_err(|_| { + IonParsingError::parse(value, "Unable to parse internal-fragment start") })?; - let end = b.parse::().map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse internal-fragment end"), + let end = b.parse::().map_err(|_| { + IonParsingError::parse(value, "Unable to parse internal-fragment end") })?; return Self::try_new_internal(start, end, charge, isotope, loss); } @@ -472,8 +529,13 @@ impl IonAnnot { let mut ch = res.chars(); return match (ch.next(), ch.next()) { (Some(r), None) => Self::try_new_immonium(r, charge, isotope, loss), - // `IC[Carbamidomethyl]` and friends carry a mod string that no - // fixed-width field can hold. + // A bare `I` names no residue at all; `IC[Carbamidomethyl]` and + // friends carry a mod string no fixed-width field can hold. + // Both are unrepresentable, and the message says which it is. + (None, _) => Err(IonParsingError::parse( + value, + "Immonium ion names no residue", + )), _ => Err(IonParsingError::UnsupportedModifiedImmonium { annotation: value.to_string(), }), @@ -482,22 +544,18 @@ impl IonAnnot { // Backbone / precursor / unknown: a series char then an ordinal. let mut chars = core.chars(); - let series = chars.next().ok_or(IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Empty string"), - })?; + let series = chars + .next() + .ok_or(IonParsingError::parse(value, "Empty string"))?; let rest = chars.as_str(); - let ordinal = if rest.is_empty() { - None - } else { - Some( - rest.parse::() - .map_err(|e| IonParsingError::ParsingError { - error: format!("{rest} -> {e:?}"), - context: Some("Unable to parse the ordinal number"), - })?, - ) - }; + let ordinal = + if rest.is_empty() { + None + } else { + Some(rest.parse::().map_err(|_| { + IonParsingError::parse(value, "Ordinal is not a number in 0..=255") + })?) + }; Self::try_new_with_loss(series, ordinal, charge, isotope, loss) } } @@ -535,33 +593,56 @@ impl Display for IonAnnot { } } +/// Why an annotation could not be parsed or represented. +/// +/// Never matched outside this crate today, so the variants exist for the +/// message a user sees when a library fails to load. Each one names the +/// offending value, because the readers that surface these discard the variant +/// and show only the rendered string. #[derive(Debug, Error)] pub enum IonParsingError { #[error("Ordinal {ordinal} out of range for series '{series}'")] OrdinalOutOfRange { ordinal: u8, series: char }, #[error("Series '{series}' requires an ordinal")] MissingOrdinal { series: char }, + #[error("Series '{series}' takes no ordinal, got {ordinal}")] + UnexpectedOrdinal { series: char, ordinal: u8 }, #[error("Unsupported fragment type: '{fragment_type}'")] UnsupportedFragmentType { fragment_type: char }, #[error("Charge cannot be 0")] ChargeCannotBeZero, - #[error("Charge {charge} outside the representable range")] + #[error("Charge {charge} outside the representable range {CHARGE_MIN}..={CHARGE_MAX}")] ChargeOutOfRange { charge: i8 }, - #[error("Isotope offset {isotope} outside the representable range")] + #[error( + "Isotope offset {isotope} outside the representable range {ISOTOPE_MIN}..={ISOTOPE_MAX}" + )] IsotopeOutOfRange { isotope: i8 }, #[error("Neutral loss '{loss}' is not representable")] UnsupportedNeutralLoss { loss: String }, - #[error("Modified immonium ions are not representable: '{annotation}'")] + #[error("Immonium ions must be a bare uppercase residue, got '{annotation}'")] UnsupportedModifiedImmonium { annotation: String }, #[error("Ran out of distinct unknown-ion labels: the 8-bit ordinal is exhausted")] UnknownIonsExhausted, - #[error("Parsing error: {error}{}", .context.map(|c| format!(" ({})", c)).unwrap_or_default())] + #[error("Could not parse '{annotation}': {context}")] ParsingError { - error: String, - context: Option<&'static str>, + /// The whole annotation, not the fragment of it that failed: a reader + /// reports this with no row index, so the full text is the only handle + /// the user gets on which row broke. + annotation: String, + context: &'static str, }, } +impl IonParsingError { + /// Shorthand for [`Self::ParsingError`], which is built at thirteen sites. + fn parse(annotation: &str, context: &'static str) -> Self { + Self::ParsingError { + annotation: annotation.to_string(), + context, + } + } +} + /// Hands out `?1`, `?2`, ... for peaks whose annotation this crate cannot /// represent. /// @@ -570,16 +651,16 @@ pub enum IonParsingError { /// A monotonic counter makes that uniqueness structural, and returning an /// error once the 8-bit ordinal is spent keeps the overflow from being /// something each reader has to remember to check. -#[derive(Debug, Default, Clone, Copy)] +/// +/// Deliberately neither `Copy` nor `Clone`: a duplicated counter forks, and +/// each fork reissues labels the other already handed out -- exactly the bug +/// this type exists to prevent. Pass it by `&mut`. +#[derive(Debug, Default)] pub struct UnknownIonCounter(u8); impl UnknownIonCounter { - pub fn new() -> Self { - Self::default() - } - /// The next unused unknown label at `charge`. - pub fn next(&mut self, charge: i8) -> Result { + pub fn next_unknown(&mut self, charge: i8) -> Result { let ordinal = self .0 .checked_add(1) @@ -592,9 +673,9 @@ impl UnknownIonCounter { /// A backbone fragment ion series. /// -/// The nine mzPAF backbone series differ only by their letter, so they are one -/// enum with one letter table rather than nine variants repeated across every -/// match in this module. +/// The nine mzPAF backbone series differ only by their letter, so the +/// letter-to-discriminant pairing lives in exactly one private table and every +/// match over them is a single arm. #[derive(Clone, Copy, Eq, PartialEq, Ord, PartialOrd, Hash, Debug)] #[allow(non_camel_case_types)] #[repr(u8)] @@ -611,7 +692,8 @@ pub enum Series { } impl Series { - /// Every series, in discriminant order. Parallel to [`Self::CHARS`]. + /// Every series, in discriminant order, and parallel to the private letter + /// table -- `series_letters_and_discriminants_agree` pins that pairing. pub const ALL: [Self; 9] = [ Self::a, Self::b, @@ -727,7 +809,7 @@ impl IonSeriesOrdinal { if c == 'p' { return match ordinal { None => Ok(Self::precursor), - Some(ordinal) => Err(IonParsingError::OrdinalOutOfRange { ordinal, series: c }), + Some(ordinal) => Err(IonParsingError::UnexpectedOrdinal { series: c, ordinal }), }; } let ordinal = ordinal.ok_or(IonParsingError::MissingOrdinal { series: c })?; @@ -759,13 +841,32 @@ mod tests { IonAnnot::try_from(s).unwrap_or_else(|e| panic!("{s:?} must parse: {e}")) } - /// The whole point of the packed representation. + /// The size claim the design rests on: a loss, an internal-fragment span + /// and an immonium residue all ride along without growing the tuple that + /// timsseek stores inline. + /// + /// The predecessor was also 4 bytes, so `size_of` alone pins nothing -- + /// what this asserts is that the *added* fields cost nothing, which is only + /// meaningful together with the round-trip tests proving they are really in + /// there. #[test] - fn packs_into_one_word() { + fn the_added_fields_cost_no_space() { assert_eq!(size_of::(), 4); // Paired with an intensity on the scoring hot path; padding here would - // grow the inline TinyVec storage in `ExpectedIntensities`. + // grow the inline TinyVec storage in `ExpectedIntensities` from 104 to + // 156 bytes at the current inline capacity of 13. assert_eq!(size_of::<(IonAnnot, f32)>(), 8); + + // All four of these are new capacity, and none of them widened the word. + let loaded = IonAnnot::try_new_internal(2, 11, 3, 2, NeutralLoss::PhosphoricAcidWater) + .expect("every field at once"); + assert_eq!(size_of_val(&loaded), 4); + assert_eq!(loaded.loss(), NeutralLoss::PhosphoricAcidWater); + assert_eq!( + loaded.series_ordinal(), + IonSeriesOrdinal::internal { start: 2, end: 11 } + ); + assert_eq!((loaded.get_charge(), loaded.get_isotope()), (3, 2)); } /// `IonAnnot: Default` is not optional -- `tinyvec::Array` requires @@ -855,6 +956,31 @@ mod tests { )); } + /// The crate docs say negative isotope offsets are unsupported. When the + /// field held them, `Display` emitted `y5+-2i` -- not mzPAF, and reparsed + /// happily, so a library could round-trip through serde into text no other + /// mzPAF reader accepts. The field is unsigned so that spelling is + /// unreachable rather than merely undocumented. + #[test] + fn negative_isotopes_are_unrepresentable_not_just_unparsed() { + assert!(matches!( + IonAnnot::try_new('y', Some(5), 1, -1), + Err(IonParsingError::IsotopeOutOfRange { isotope: -1 }) + )); + assert!(ion("y5").try_with_offset_neutrons(-1).is_err()); + // mzPAF's own spelling for a negative offset is still not parsed. + assert!(IonAnnot::try_from("y5-2i").is_err()); + // And the spelling that used to leak out is not accepted either. + assert!(IonAnnot::try_from("y5+-2i").is_err()); + // Every representable isotope renders as something that parses back. + for isotope in ISOTOPE_MIN..=ISOTOPE_MAX { + let a = IonAnnot::try_new('y', Some(5), 1, isotope).expect("in range"); + let text = a.to_string(); + assert!(!text.contains("+-"), "{text} is not mzPAF"); + assert_eq!(ion(&text), a, "{text}"); + } + } + #[test] fn isotope_offset_respects_the_field_bound() { let a = ion("y5"); diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs index 2fcea755..b5d3c2a5 100644 --- a/rust/micromzpaf/src/loss.rs +++ b/rust/micromzpaf/src/loss.rs @@ -62,10 +62,7 @@ impl Composition { /// "not representable" instead of being mistaken for a smaller one. fn parse_formula(s: &str) -> Result { if s.is_empty() { - return Err(IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Empty neutral-loss formula"), - }); + return Err(IonParsingError::parse(s, "Empty neutral-loss formula")); } let mut out = Composition::default(); let b = s.as_bytes(); @@ -80,12 +77,9 @@ impl Composition { let count: u8 = if start == i { 1 } else { - s[start..i] - .parse() - .map_err(|_| IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Neutral-loss atom count out of range"), - })? + s[start..i].parse().map_err(|_| { + IonParsingError::parse(s, "Neutral-loss atom count out of range") + })? }; let slot = match elem { b'C' => C, @@ -95,10 +89,10 @@ impl Composition { b'S' => S, b'P' => P, _ => { - return Err(IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Unsupported element in neutral loss"), - }); + return Err(IonParsingError::parse( + s, + "Unsupported element in neutral loss", + )); } }; out.0[slot] = out.0[slot].saturating_add(count); @@ -116,22 +110,19 @@ impl Composition { for term in s.split('-') { let term = term.trim(); if term.is_empty() { - return Err(IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Empty term in neutral-loss expression"), - }); + return Err(IonParsingError::parse( + s, + "Empty term in neutral-loss expression", + )); } // Leading digits are a repeat count for the whole term. let digits = term .find(|c: char| !c.is_ascii_digit()) .unwrap_or(term.len()); let (mult, formula) = if digits > 0 { - let m: u8 = term[..digits] - .parse() - .map_err(|_| IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Neutral-loss repeat count out of range"), - })?; + let m: u8 = term[..digits].parse().map_err(|_| { + IonParsingError::parse(s, "Neutral-loss repeat count out of range") + })?; (m, &term[digits..]) } else { (1, term) @@ -249,27 +240,28 @@ const TABLE: &[(Composition, NeutralLoss, &str)] = &[ ]; impl NeutralLoss { + /// How many losses the table names, so the width of the packed field can be + /// checked against it rather than assumed to be roomy. + pub(crate) const COUNT: u8 = TABLE.len() as u8; + /// Inverse of the `#[repr(u8)]` discriminant, for unpacking out of a bit /// field. Lives next to the enum so the two cannot drift apart. /// /// An unrecognized value maps to [`Self::None`]: the only way to produce /// one is a reserved discriminant, which no constructor emits. pub(crate) fn from_discriminant(d: u8) -> Self { - match d { - 1 => Self::Water, - 2 => Self::Ammonia, - 3 => Self::CarbonMonoxide, - 4 => Self::CarbonDioxide, - 5 => Self::WaterX2, - 6 => Self::AmmoniaX2, - 7 => Self::WaterAmmonia, - 8 => Self::Methanesulfenic, - 9 => Self::Carbamidomethylthiol, - 10 => Self::PhosphoricAcid, - 11 => Self::Metaphosphoric, - 12 => Self::PhosphoricAcidWater, - _ => Self::None, - } + Self::at(d).map_or(Self::None, |(_, loss, _)| *loss) + } + + /// The [`TABLE`] row a discriminant names, or `None` for [`Self::None`] and + /// for the reserved values above the table. + /// + /// `TABLE` is in discriminant order starting at 1, which + /// `table_is_indexed_by_discriminant` pins. Indexing it rather than + /// re-listing the variants is what keeps the discriminant, the composition + /// and the spelling from drifting apart. + fn at(d: u8) -> Option<&'static (Composition, NeutralLoss, &'static str)> { + TABLE.get((d as usize).checked_sub(1)?) } /// Resolve a loss expression (without the leading `-`) to a discriminant. @@ -288,14 +280,7 @@ impl NeutralLoss { /// Canonical spelling, without the leading `-`. Empty for [`Self::None`]. pub(crate) fn canonical(self) -> &'static str { - if self == NeutralLoss::None { - return ""; - } - TABLE - .iter() - .find(|(_, l, _)| *l == self) - .map(|(_, _, s)| *s) - .unwrap_or("") + Self::at(self as u8).map_or("", |(_, _, spelling)| *spelling) } } @@ -389,10 +374,31 @@ mod tests { assert!(NeutralLoss::from_expression("H2O-").is_err()); } + /// `from_discriminant` and `canonical` both index [`TABLE`] by discriminant, + /// so a row sitting at the wrong offset would decode as its neighbour. This + /// is the one invariant that keeps the enum, the composition and the + /// spelling in step. + #[test] + fn table_is_indexed_by_discriminant() { + for (i, (_comp, loss, canon)) in TABLE.iter().enumerate() { + assert_eq!( + *loss as usize, + i + 1, + "{canon} sits at offset {i} but its discriminant is {}", + *loss as u8 + ); + } + // Nothing outside the table decodes to a loss, in either direction. + assert_eq!(NeutralLoss::from_discriminant(0), NeutralLoss::None); + assert_eq!( + NeutralLoss::from_discriminant(TABLE.len() as u8 + 1), + NeutralLoss::None + ); + assert_eq!(NeutralLoss::None.canonical(), ""); + } + /// Every table entry must survive canonical -> composition -> discriminant, - /// and back out through the bit field. `from_discriminant` hand-mirrors the - /// `#[repr(u8)]` values, so adding a loss without updating it would - /// silently decode as [`NeutralLoss::None`] -- this is what catches that. + /// and back out through the bit field. #[test] fn table_round_trips_through_canonical_spelling() { for (_comp, loss, canon) in TABLE { @@ -408,16 +414,6 @@ mod tests { "{canon} does not survive the discriminant round trip" ); } - assert_eq!(NeutralLoss::from_discriminant(0), NeutralLoss::None); - // Every non-None variant must be in TABLE, or it has no spelling and no - // composition and could never be produced by parsing. - assert_eq!( - TABLE.len(), - (1..=u8::MAX) - .filter(|d| NeutralLoss::from_discriminant(*d) != NeutralLoss::None) - .count(), - "a variant is decodable but missing from TABLE" - ); } /// Compositions must be unique: two entries sharing one would make the diff --git a/rust/timsquery/src/serde/diann_io.rs b/rust/timsquery/src/serde/diann_io.rs index 2392a41b..2fc2757b 100644 --- a/rust/timsquery/src/serde/diann_io.rs +++ b/rust/timsquery/src/serde/diann_io.rs @@ -86,7 +86,7 @@ impl From for DiannReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct DiannPrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -414,7 +414,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(rows.len()); - let mut unknown_ions = UnknownIonCounter::new(); + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -434,7 +434,7 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - let ion_annot = unknown_ions.next(frag_charge as i8)?; + let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -743,7 +743,7 @@ fn parse_precursor_group_from_parquet( let mut fragment_mzs = Vec::with_capacity(indices.len()); buffers.fragment_labels.clear(); let mut rel_intensities = Vec::with_capacity(indices.len()); - let mut unknown_ions = UnknownIonCounter::new(); + let mut unknown_ions = UnknownIonCounter::default(); for (i, &idx) in indices.iter().enumerate() { let fragment_mz = columns.product_mzs[idx] as f64; @@ -765,7 +765,7 @@ fn parse_precursor_group_from_parquet( columns.fragment_loss_types[idx], i ); - let ion_annot = unknown_ions.next(frag_charge as i8)?; + let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); rel_intensities.push((ion_annot, rel_intensity)); @@ -965,29 +965,28 @@ AAAAAAALQAK\tAAAAAAALQAK\t478.7\t2\t11.0\t0.9\tP2\t0\t300.0\ty\t3\t1\tnoloss\t1. let mgrysgk = &elution_groups[1].0; let hgdtgrr = &elution_groups[0].0; - let mut mgrysgk_expected_labels = vec!["y6", "b6", "b3", "b5", "y4"] - .into_iter() - .map(|s| IonAnnot::try_from(s).unwrap()) - .collect::>(); - let mut actual_labels: Vec = mgrysgk - .iter_fragments() - .map(|(label, _mz)| *label) - .collect(); - mgrysgk_expected_labels.sort(); - actual_labels.sort(); - assert_eq!(actual_labels, mgrysgk_expected_labels); - - let mut hgdtgrr_expected_labels = vec!["y6", "b3", "y4", "y5"] - .into_iter() - .map(|s| IonAnnot::try_from(s).unwrap()) - .collect::>(); - let mut actual_labels: Vec = hgdtgrr - .iter_fragments() - .map(|(label, _mz)| *label) - .collect(); - hgdtgrr_expected_labels.sort(); - actual_labels.sort(); - assert_eq!(actual_labels, hgdtgrr_expected_labels); + // Labels as a set: the reader's fragment order is not part of the + // contract being tested here. Keyed on the mzPAF spelling because + // `IonAnnot` is deliberately not `Ord`. + fn label_set(eg: &Target) -> Vec { + let mut out: Vec = eg + .iter_fragments() + .map(|(label, _mz)| label.to_string()) + .collect(); + out.sort(); + out + } + fn expected_set(labels: &[&str]) -> Vec { + let mut out: Vec = labels.iter().map(|s| s.to_string()).collect(); + out.sort(); + out + } + + assert_eq!( + label_set(mgrysgk), + expected_set(&["y6", "b6", "b3", "b5", "y4"]) + ); + assert_eq!(label_set(hgdtgrr), expected_set(&["y6", "b3", "y4", "y5"])); } #[test] diff --git a/rust/timsquery/src/serde/diann_speclib_io.rs b/rust/timsquery/src/serde/diann_speclib_io.rs index 6c36efa3..1ab6c73e 100644 --- a/rust/timsquery/src/serde/diann_speclib_io.rs +++ b/rust/timsquery/src/serde/diann_speclib_io.rs @@ -764,7 +764,9 @@ fn map_entry( if f.typ() & 0x80 != 0 { stats.exclude_flagged += 1; } - // IonAnnot cannot represent neutral loss; drop lossy fragments. + // `IonAnnot` can represent neutral losses, but the map from DIA-NN's + // loss code byte to `NeutralLoss` is not written, so these are dropped. + // The `loss_dropped` assertion in the tests below measures the cost. #105 if f.loss() != 0 { stats.loss_dropped += 1; continue; diff --git a/rust/timsquery/src/serde/elution_group_inputs.rs b/rust/timsquery/src/serde/elution_group_inputs.rs index 97b6533c..0c05675d 100644 --- a/rust/timsquery/src/serde/elution_group_inputs.rs +++ b/rust/timsquery/src/serde/elution_group_inputs.rs @@ -1,4 +1,7 @@ -use crate::ion::IonAnnot; +use crate::ion::{ + IonAnnot, + UnknownIonCounter, +}; use crate::tinyvec::{ TinyVec, tiny_vec, @@ -10,10 +13,21 @@ use crate::{ #[derive(Debug)] pub enum ElutionGroupInputError { - MismatchedFragmentLabelsLength { expected: usize, found: usize }, + MismatchedFragmentLabelsLength { + expected: usize, + found: usize, + }, AlreadyHasFragmentLabels, - IonConversionError { inner: String }, + IonConversionError { + inner: String, + }, MissingFragmentLabels, + /// More fragments than the `u8` label space holds. Wrapping instead would + /// mint a duplicate label, and lookup is by first match, so every later + /// fragment sharing it would be unreachable. + TooManyFragmentsToLabel { + found: usize, + }, } /// User-friendly format for specifying elution groups in an input file @@ -44,6 +58,11 @@ impl ElutionGroupInput { if self.fragment_labels.is_some() { return Err(ElutionGroupInputError::AlreadyHasFragmentLabels); } + if num_fragments > u8::MAX as usize { + return Err(ElutionGroupInputError::TooManyFragmentsToLabel { + found: num_fragments, + }); + } let fragment_labels: Vec = (0..num_fragments).map(|i| i as u8).collect(); Ok(ElutionGroupInput { @@ -58,24 +77,34 @@ impl ElutionGroupInput { }) } + /// Fill in `?1`, `?2`, ... for an input that named no fragment labels. + /// + /// Minted through [`UnknownIonCounter`] rather than an index cast, so + /// running past the label space is an error rather than a wrap into + /// duplicate labels. pub fn try_fill_labels_annot( self, ) -> Result, ElutionGroupInputError> { - let tmp = self.try_fill_labels_u8()?; - let new_frags = tmp - .fragment_labels - .unwrap() - .into_iter() - .map(|lbl| IonAnnot::try_new('?', Some(lbl), 1, 1).unwrap()) - .collect(); + if self.fragment_labels.is_some() { + return Err(ElutionGroupInputError::AlreadyHasFragmentLabels); + } + let mut unknown_ions = UnknownIonCounter::default(); + let new_frags = self + .fragments + .iter() + .map(|_| unknown_ions.next_unknown(1)) + .collect::, _>>() + .map_err(|e| ElutionGroupInputError::IonConversionError { + inner: e.to_string(), + })?; Ok(ElutionGroupInput { - id: tmp.id, - mobility: tmp.mobility, - rt_seconds: tmp.rt_seconds, - precursor: tmp.precursor, - precursor_charge: tmp.precursor_charge, - precursor_isotopes: tmp.precursor_isotopes, - fragments: tmp.fragments, + id: self.id, + mobility: self.mobility, + rt_seconds: self.rt_seconds, + precursor: self.precursor, + precursor_charge: self.precursor_charge, + precursor_isotopes: self.precursor_isotopes, + fragments: self.fragments, fragment_labels: Some(new_frags), }) } diff --git a/rust/timsquery/src/serde/skyline_io.rs b/rust/timsquery/src/serde/skyline_io.rs index a18da490..84e966de 100644 --- a/rust/timsquery/src/serde/skyline_io.rs +++ b/rust/timsquery/src/serde/skyline_io.rs @@ -80,7 +80,7 @@ impl From for SkylineReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct SkylinePrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -310,7 +310,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(fragment_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(fragment_rows.len()); - let mut unknown_ions = UnknownIonCounter::new(); + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in fragment_rows.iter().enumerate() { let fragment_mz = row.product_mz; @@ -349,7 +349,7 @@ fn parse_precursor_group( falling back to unknown ion", frag_type, i ); - unknown_ions.next(frag_charge as i8)? + unknown_ions.next_unknown(frag_charge as i8)? } }; diff --git a/rust/timsquery/src/serde/spectronaut_io.rs b/rust/timsquery/src/serde/spectronaut_io.rs index 2c1a6680..f12a0ba3 100644 --- a/rust/timsquery/src/serde/spectronaut_io.rs +++ b/rust/timsquery/src/serde/spectronaut_io.rs @@ -80,7 +80,7 @@ impl From for SpectronautReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct SpectronautPrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -284,7 +284,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(included_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(included_rows.len()); - let mut unknown_ions = UnknownIonCounter::new(); + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in included_rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -304,7 +304,7 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - let ion_annot = unknown_ions.next(frag_charge as i8)?; + let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); From e412151e5a6747b0ec2ce175e0ee4a5e8032a3dc Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 11:06:02 -0700 Subject: [PATCH 3/8] refactor(micromzpaf): name the atoms in Composition, trim the crate doc PR review notes. `Composition` becomes a struct of named `u8`s instead of `[u8; 6]` plus six index constants, so `TABLE` reads as chemistry: `H2O` is `h: 2, o: 1` rather than `Composition::new(0, 2, 0, 1, 0, 0)`, where a transposition is invisible on review and would alias one loss onto another. The symbol-to-field mapping now lives once, in `count_mut`, which the parser indexes instead of keeping its own list of slots. Rows use a local `C!` macro: a plain struct literal has to spell all six counts, which is the positional noise the struct is meant to remove, and `..ZERO` is not allowed bare in a `const` item. Crate doc: drop the sentence announcing what the argument is instead of making it, and retitle "mzPAF compliance" -- this parses a subset, so the heading overstated it. The examples now show a neutral loss with its two spellings, an internal fragment, an immonium ion, the mass-error suffix, and two rejections, rather than two variations on `b12`. --- rust/micromzpaf/src/lib.rs | 63 +++++++----- rust/micromzpaf/src/loss.rs | 191 ++++++++++++++++++++---------------- 2 files changed, 140 insertions(+), 114 deletions(-) diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 55f0e7dd..69a63c52 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -5,13 +5,10 @@ //! the annotations an mzSpecLib-shaped library needs -- neutral losses, //! internal fragments, immonium ions -- fit *without* growing the type. //! -//! That is the whole argument, and it is about the alternative rather than -//! about the previous representation. The predecessor was a 4-byte struct of -//! `(series+ordinal, charge, isotope)`; bolting a loss field onto it would have -//! cost three bytes, not one, because a 5-byte key paired with an `f32` pads to -//! a 12-byte tuple. That would have grown the inline `TinyVec` storage in -//! timsseek's `ExpectedIntensities` from 104 to 156 bytes. Packed, the tuple -//! stays 8 bytes and the loss rides along in bits nobody was using. +//! Bolting a loss field onto a struct of fields would have cost three bytes, +//! not one: a 5-byte key paired with an `f32` pads to a 12-byte tuple, growing +//! the inline `TinyVec` storage in timsseek's `ExpectedIntensities` from 104 to +//! 156 bytes. Packed, the tuple stays 8 bytes. //! //! # Bit layout //! @@ -40,38 +37,50 @@ //! because `IonAnnot: Default` is forced by `tinyvec::Array` and a default can //! reach any serde path. //! -//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned, -//! because a negative offset has no spelling this crate can emit (see -//! [`ISOTOPE_MIN`]). Both ranges are far wider than anything observed: the -//! HUPO-PSI corpus tops out at charge 3 and isotope 3, with no negative charges -//! at all. Because a bit field truncates rather than wrapping loudly, every -//! constructor range-checks -- see [`IonAnnot::try_new`]. +//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned +//! (see [`ISOTOPE_MIN`]). A bit field truncates rather than wrapping loudly, so +//! every constructor range-checks -- see [`IonAnnot::try_new`]. //! -//! # mzPAF compliance +//! # The mzPAF subset //! -//! Supported: the a/b/c/d/v/w/x/y/z series, precursor (`p`), unknown (`?`), +//! Parses: the a/b/c/d/v/w/x/y/z series, precursor (`p`), unknown (`?`), //! internal fragments (`m:`), bare immonium (`IA`), charge (`^N`), //! positive isotopes (`+Ni`), neutral losses from [`NeutralLoss`], and the //! mass-error suffix (`/-0.0003`, `/1.2ppm`). //! -//! Not supported: negative isotope offsets, modified immonium -//! (`IC[Carbamidomethyl]` carries an arbitrary mod string), and losses outside -//! the [`NeutralLoss`] table. These are reported as errors, never coerced into -//! a nearby representable ion -- and unrepresentable in the field too, so -//! `Display` cannot emit a spelling the parser would reject. +//! Rejects, rather than coercing onto a nearby representable ion: negative +//! isotope offsets, modified immonium (`IC[Carbamidomethyl]` carries an +//! arbitrary mod string), and losses outside the [`NeutralLoss`] table. //! //! # Examples //! //! ``` -//! use micromzpaf::IonAnnot; +//! use micromzpaf::{IonAnnot, NeutralLoss, split_mass_error}; //! -//! // Parse a simple b-ion annotation -//! let ion: IonAnnot = "b12".try_into().unwrap(); -//! assert_eq!(format!("{}", ion), "b12"); -//! -//! // Parse with charge and isotope +//! // A backbone ion, with charge and isotope suffixes. //! let ion: IonAnnot = "b12+i^3".try_into().unwrap(); -//! assert_eq!(ion.get_charge(), 3); +//! assert_eq!(ion.try_get_ordinal(), Some(12)); +//! assert_eq!((ion.get_charge(), ion.get_isotope()), (3, 1)); +//! +//! // A neutral loss. Two spellings of one chemical loss are one annotation, +//! // and render canonically. +//! let ion: IonAnnot = "y5-CH3SOH".try_into().unwrap(); +//! assert_eq!(ion.loss(), NeutralLoss::Methanesulfenic); +//! assert_eq!(ion.to_string(), "y5-CH4OS"); +//! +//! // An internal fragment spanning residues 2..=11, and a bare immonium ion. +//! // Neither sits on a ladder, so neither has an ordinal. +//! assert_eq!(IonAnnot::try_from("m2:11").unwrap().try_get_ordinal(), None); +//! assert_eq!(IonAnnot::try_from("IA").unwrap().to_string(), "IA"); +//! +//! // mzSpecLib peaks carry observed m/z; the suffix recovers theoretical. +//! let (ion, err) = split_mass_error("y1/-0.0005").unwrap(); +//! assert_eq!(ion, "y1"); +//! assert!((err.unwrap().theoretical_from_observed(175.1184) - 175.1189).abs() < 1e-4); +//! +//! // Unrepresentable annotations fail rather than losing the detail. +//! assert!(IonAnnot::try_from("y1-HCOOH").is_err()); +//! assert!(IonAnnot::try_from("IC[Carbamidomethyl]").is_err()); //! ``` pub mod loss; diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs index b5d3c2a5..7009372b 100644 --- a/rust/micromzpaf/src/loss.rs +++ b/rust/micromzpaf/src/loss.rs @@ -23,36 +23,93 @@ use std::fmt::Display; use crate::IonParsingError; -/// Slots in [`Composition`], in the order [`Composition::new`] takes them. -const C: usize = 0; -const H: usize = 1; -const N: usize = 2; -const O: usize = 3; -const S: usize = 4; -const P: usize = 5; - /// Atom counts for the elements that appear in peptide neutral losses. /// /// Deliberately not a general chemical formula: these losses only ever draw -/// from C/H/N/O/S/P, and keeping it to six `u8`s makes equality a single -/// 6-byte compare during the parse-time table lookup. +/// from C/H/N/O/S/P, and six named `u8`s make equality a 6-byte compare during +/// the parse-time table lookup. +/// +/// Named fields rather than `[u8; 6]` so that [`TABLE`] reads as chemistry +/// (`H2O` is `h: 2, o: 1`) instead of six positional numbers, where a +/// transposition would be invisible on review and would silently alias one +/// loss onto another. #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] -pub(crate) struct Composition([u8; 6]); +pub(crate) struct Composition { + c: u8, + h: u8, + n: u8, + o: u8, + s: u8, + p: u8, +} + +/// A [`Composition`] naming only the elements it contains: `C!(h: 2, o: 1)`. +/// +/// A plain struct literal would have to spell all six counts, which is the +/// positional noise this struct exists to remove; functional update syntax +/// (`..ZERO`) is not permitted in a `const` item. +macro_rules! C { + ($($field:ident: $count:expr),+ $(,)?) => { + Composition { $($field: $count,)+ ..Composition::ZERO } + }; +} impl Composition { - pub(crate) const fn new(c: u8, h: u8, n: u8, o: u8, s: u8, p: u8) -> Self { - Self([c, h, n, o, s, p]) + /// All-zero, for [`TABLE`] rows to fill in only the elements they contain. + const ZERO: Self = Self { + c: 0, + h: 0, + n: 0, + o: 0, + s: 0, + p: 0, + }; + + /// The count for one element symbol, or `None` if this crate does not + /// represent that element. + /// + /// The single place the symbol-to-field mapping lives, so the parser cannot + /// disagree with the struct about which letter means which count. + fn count_mut(&mut self, symbol: u8) -> Option<&mut u8> { + Some(match symbol { + b'C' => &mut self.c, + b'H' => &mut self.h, + b'N' => &mut self.n, + b'O' => &mut self.o, + b'S' => &mut self.s, + b'P' => &mut self.p, + _ => return None, + }) } - /// Multiply every count, saturating. Used for the `2H2O` multiplier form. + /// The six counts, for tests that need to range over them. + #[cfg(test)] + fn counts(self) -> [u8; 6] { + [self.c, self.h, self.n, self.o, self.s, self.p] + } + + /// Combine element-wise. Saturating: the only inputs that reach the ceiling + /// are absurd (`200H2O`), and no [`TABLE`] row holds a saturated count, so a + /// saturated result cannot alias onto a real loss. Pinned by + /// `table_compositions_are_unique`. + fn zip(self, other: Self, f: impl Fn(u8, u8) -> u8) -> Self { + Self { + c: f(self.c, other.c), + h: f(self.h, other.h), + n: f(self.n, other.n), + o: f(self.o, other.o), + s: f(self.s, other.s), + p: f(self.p, other.p), + } + } + + /// Multiply every count. Used for the `2H2O` multiplier form. fn scaled(self, k: u8) -> Self { - Self(self.0.map(|n| n.saturating_mul(k))) + self.zip(Self::ZERO, |n, _| n.saturating_mul(k)) } fn plus(self, other: Self) -> Self { - Self(std::array::from_fn(|i| { - self.0[i].saturating_add(other.0[i]) - })) + self.zip(other, u8::saturating_add) } /// Parse a bare formula like `H2O`, `CH4OS`, `C2H5NOS`. @@ -68,7 +125,7 @@ impl Composition { let b = s.as_bytes(); let mut i = 0; while i < b.len() { - let elem = b[i]; + let symbol = b[i]; i += 1; let start = i; while i < b.len() && b[i].is_ascii_digit() { @@ -81,21 +138,10 @@ impl Composition { IonParsingError::parse(s, "Neutral-loss atom count out of range") })? }; - let slot = match elem { - b'C' => C, - b'H' => H, - b'N' => N, - b'O' => O, - b'S' => S, - b'P' => P, - _ => { - return Err(IonParsingError::parse( - s, - "Unsupported element in neutral loss", - )); - } - }; - out.0[slot] = out.0[slot].saturating_add(count); + let slot = out + .count_mut(symbol) + .ok_or_else(|| IonParsingError::parse(s, "Unsupported element in neutral loss"))?; + *slot = slot.saturating_add(count); } Ok(out) } @@ -171,69 +217,40 @@ pub enum NeutralLoss { PhosphoricAcidWater = 12, } -/// `(composition, discriminant, canonical spelling)`. +/// `(composition, discriminant, canonical spelling)`, indexed by discriminant. +/// +/// Rows must stay in discriminant order starting at 1 -- +/// `table_is_indexed_by_discriminant` pins that, and both `from_discriminant` +/// and `canonical` index straight into this. /// /// The canonical spelling is what `Display` emits, so a non-canonical input /// (`-CH3SOH`) round-trips to the canonical form (`-CH4OS`). Round-trip tests /// must therefore compare parsed values, not bytes. +/// +/// Compositions name their elements, so each row can be read against its own +/// spelling. const TABLE: &[(Composition, NeutralLoss, &str)] = &[ + (C!(h: 2, o: 1), NeutralLoss::Water, "H2O"), + (C!(h: 3, n: 1), NeutralLoss::Ammonia, "NH3"), + (C!(c: 1, o: 1), NeutralLoss::CarbonMonoxide, "CO"), + (C!(c: 1, o: 2), NeutralLoss::CarbonDioxide, "CO2"), + (C!(h: 4, o: 2), NeutralLoss::WaterX2, "2H2O"), + (C!(h: 6, n: 2), NeutralLoss::AmmoniaX2, "2NH3"), + (C!(h: 5, n: 1, o: 1), NeutralLoss::WaterAmmonia, "H2O-NH3"), ( - Composition::new(0, 2, 0, 1, 0, 0), - NeutralLoss::Water, - "H2O", - ), - ( - Composition::new(0, 3, 1, 0, 0, 0), - NeutralLoss::Ammonia, - "NH3", - ), - ( - Composition::new(1, 0, 0, 1, 0, 0), - NeutralLoss::CarbonMonoxide, - "CO", - ), - ( - Composition::new(1, 0, 0, 2, 0, 0), - NeutralLoss::CarbonDioxide, - "CO2", - ), - ( - Composition::new(0, 4, 0, 2, 0, 0), - NeutralLoss::WaterX2, - "2H2O", - ), - ( - Composition::new(0, 6, 2, 0, 0, 0), - NeutralLoss::AmmoniaX2, - "2NH3", - ), - ( - Composition::new(0, 5, 1, 1, 0, 0), - NeutralLoss::WaterAmmonia, - "H2O-NH3", - ), - ( - Composition::new(1, 4, 0, 1, 1, 0), + C!(c: 1, h: 4, o: 1, s: 1), NeutralLoss::Methanesulfenic, "CH4OS", ), ( - Composition::new(2, 5, 1, 1, 1, 0), + C!(c: 2, h: 5, n: 1, o: 1, s: 1), NeutralLoss::Carbamidomethylthiol, "C2H5NOS", ), + (C!(h: 3, o: 4, p: 1), NeutralLoss::PhosphoricAcid, "H3PO4"), + (C!(h: 1, o: 3, p: 1), NeutralLoss::Metaphosphoric, "HPO3"), ( - Composition::new(0, 3, 0, 4, 0, 1), - NeutralLoss::PhosphoricAcid, - "H3PO4", - ), - ( - Composition::new(0, 1, 0, 3, 0, 1), - NeutralLoss::Metaphosphoric, - "HPO3", - ), - ( - Composition::new(0, 5, 0, 5, 0, 1), + C!(h: 5, o: 5, p: 1), NeutralLoss::PhosphoricAcidWater, "H3PO4-H2O", ), @@ -303,17 +320,17 @@ mod tests { fn atom_counts_parse_and_reject_out_of_range() { assert_eq!( Composition::parse_expression("H2O").unwrap(), - Composition::new(0, 2, 0, 1, 0, 0), + C!(h: 2, o: 1), "an element with no digits is one atom" ); assert_eq!( Composition::parse_expression("C10H12").unwrap(), - Composition::new(10, 12, 0, 0, 0, 0), + C!(c: 10, h: 12), "counts are multi-digit, not one digit per element" ); assert_eq!( Composition::parse_expression("C255").unwrap(), - Composition::new(255, 0, 0, 0, 0, 0), + C!(c: 255), "the u8 slot is full at 255" ); assert!( @@ -429,7 +446,7 @@ mod tests { // unrepresentable while no table entry holds a saturated count -- // otherwise the saturation would alias onto a real loss. assert!( - a.0.iter().all(|&n| n < u8::MAX), + a.counts().iter().all(|&n| n < u8::MAX), "{sa} holds a saturated atom count" ); } From d7381b6025d9a41f890520856b12d56458ac5627 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 11:55:14 -0700 Subject: [PATCH 4/8] fix: name the row that failed, and make cargo doc pass Reader errors said "DIA-NN precursor parsing error" and nothing else. The `From` impl that flattened them had no way to see the row, so it is replaced by a `map_err` adaptor that stamps one on, and the three `*ReadingError`s now carry the precursor error instead of discarding it. A bad row aborts the whole load, so the row index is the only handle the user gets. `CHARGE_MAX` is 8, not 7. The bias means the field holds `-7..=8`; capping at 7 rejected a charge the encoding can represent. The `const` assertion is what keeps this honest, and it is why `CHARGE_MIN` cannot follow suit: `-8` zigzags to 17 and truncates to 1, decoding as charge 0. The `.speclib` fragment charge is a raw wire byte and was cast with `as i8`, turning anything above 127 into a plausible negative charge. It is range-checked now. The sibling `typ()` is masked `& 0x7F` for DIA-NN flag bits; if this byte carries flags too, this is where it surfaces (#105). The isotope-offset failure said the situation "should never happen" instead of what the limit is. It now names the fragment and the representable range. `TryFrom<&str>` still drops the mass-error suffix -- the suffix belongs to one observed peak, not to the ion, and storing it would make two `b12`s unequal -- but it warns once per process rather than doing it silently. `cargo doc --workspace` now passes and has a `task doc` target. It was failing in five crates: ten links to private items, and ten "broken" links that were never links at all (`[U:4]`, `["Carbamidomethyl@C"]`, `[0,1]`, `[Apply]`). Also drops `sort_vecs_by_first!`, which had no callers outside its own tests. --- Cargo.lock | 1 + Taskfile.yml | 6 ++ rust/apex_sim/src/plots.rs | 2 +- rust/micromzpaf/Cargo.toml | 2 + rust/micromzpaf/src/lib.rs | 42 ++++++++-- rust/rescore_dash/src/app.rs | 2 +- rust/rescore_dash/src/lib.rs | 3 +- rust/speclib_build_cli/src/config.rs | 6 +- rust/speclib_build_cli/src/entry.rs | 2 +- rust/speclib_build_cli/src/mods.rs | 4 +- rust/timsquery/src/serde/diann_io.rs | 62 +++++++++++---- rust/timsquery/src/serde/diann_speclib_io.rs | 41 +++++++--- rust/timsquery/src/serde/skyline_io.rs | 43 +++++++--- rust/timsquery/src/serde/spectronaut_io.rs | 44 +++++++--- rust/timsquery/src/utils/sorting.rs | 84 -------------------- rust/timsquery_viewer/src/calibration.rs | 2 +- rust/timsseek/src/data_sources/speclib.rs | 2 +- rust/timsseek/src/ml/cv.rs | 12 +-- rust/timsseek/src/ml/qvalues.rs | 24 +++--- rust/timsseek/src/utils/elution_group_ops.rs | 9 ++- 20 files changed, 221 insertions(+), 172 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index dc2df80d..865d7112 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4022,6 +4022,7 @@ dependencies = [ "serde", "serde_json", "thiserror 2.0.18", + "tracing", ] [[package]] diff --git a/Taskfile.yml b/Taskfile.yml index 676dc1d6..809d37a0 100644 --- a/Taskfile.yml +++ b/Taskfile.yml @@ -30,6 +30,12 @@ tasks: # dashboard.rs only exists under the opt-in feature; CI only `cargo check`s it. - cargo clippy -p timsseek_cli --features dashboard --all-targets -- -D warnings + # Broken and private intra-doc links only bite the reader, so nothing else + # catches them: CI does not build docs. + doc: + cmds: + - RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps + todos: cmds: - grep -nH -R TODO rust/tims*/src diff --git a/rust/apex_sim/src/plots.rs b/rust/apex_sim/src/plots.rs index 79015208..936f2645 100644 --- a/rust/apex_sim/src/plots.rs +++ b/rust/apex_sim/src/plots.rs @@ -278,7 +278,7 @@ impl Default for TraceToggles { } /// Bottom panel: the apex-finder intermediate traces, min-max normalized to -/// [0,1] so heterogeneous scales overlay legibly. Apex markers overlaid. +/// `[0,1]` so heterogeneous scales overlay legibly. Apex markers overlaid. pub fn traces(ui: &mut egui::Ui, score: &ScoreResult, tog: &TraceToggles, true_apex: f32) { let t = &score.traces; let mut series: Vec<(&str, &[f32])> = Vec::new(); diff --git a/rust/micromzpaf/Cargo.toml b/rust/micromzpaf/Cargo.toml index e168ab37..168e12a8 100644 --- a/rust/micromzpaf/Cargo.toml +++ b/rust/micromzpaf/Cargo.toml @@ -7,6 +7,8 @@ license.workspace = true [dependencies] serde = { workspace = true, features = ["derive"] } thiserror = { workspace = true } +# Only to warn once when the lossy `TryFrom<&str>` drops a mass-error suffix. +tracing = { workspace = true } [dev-dependencies] # Only to pin that `IonAnnot` serialises as its mzPAF string rather than the diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 69a63c52..858edf1a 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -107,13 +107,15 @@ const LOSS_BITS: u32 = 6; const PAYLOAD_SHIFT: u32 = 18; const PAYLOAD_BITS: u32 = 12; -/// Widest charge the 4-bit zigzag field holds. Observed maximum is 3. +/// Charge bounds. Observed maximum in the HUPO-PSI corpus is 3. /// -/// The bias applied when packing makes `-7..=8` representable, but the range is -/// capped symmetrically: nothing has ever needed charge 8, and an asymmetric -/// public bound invites the reader to check the arithmetic rather than trust it. +/// Asymmetric because the field stores `charge - 1` zigzagged: that shifts the +/// whole window up by one, so `8` fits and `-8` does not. The `const` assertion +/// below is what holds these to the field rather than to this comment -- +/// `CHARGE_MIN = -8` would zigzag to 17 and truncate to 1, silently decoding as +/// charge 0. pub const CHARGE_MIN: i8 = -7; -pub const CHARGE_MAX: i8 = 7; +pub const CHARGE_MAX: i8 = 8; /// Widest isotope offset the 4-bit field holds. Observed maximum is 3. /// /// Isotope offsets are unsigned. mzPAF spells a negative offset `-Ni`, which @@ -569,13 +571,37 @@ impl IonAnnot { } } +/// Guards the one-shot warning in [`IonAnnot::try_from`]. +static MASS_ERROR_DISCARDED: std::sync::Once = std::sync::Once::new(); + impl TryFrom<&str> for IonAnnot { type Error = IonParsingError; - /// Parses an annotation, discarding any mass-error suffix. Use - /// [`split_mass_error`] first to keep it. + /// Parses an annotation, **discarding any mass-error suffix**, and warns + /// once per process the first time it discards a real one. + /// + /// The suffix is a property of one observed peak, not of the ion: two peaks + /// annotated `b12` in different spectra carry different errors, so keeping + /// it on the annotation would make two `b12`s unequal and break the + /// per-precursor label uniqueness the whole crate keys on. + /// + /// A caller that needs the error wants it for the *m/z*, not the label: + /// call [`split_mass_error`] first, recover the theoretical m/z with + /// [`MassError::theoretical_from_observed`], and store that. Nothing is + /// lost that way -- only the residual, which belongs to the measurement. fn try_from(value: &str) -> Result { - Self::parse_ion(split_mass_error(value)?.0) + let (ion, mass_error) = split_mass_error(value)?; + if mass_error.is_some() { + MASS_ERROR_DISCARDED.call_once(|| { + tracing::warn!( + "annotation {value:?} carries a mass-error suffix, which \ + `IonAnnot` does not store; the m/z used is whatever the \ + caller supplied. Use `split_mass_error` to recover the \ + theoretical m/z. Warning once per process." + ); + }); + } + Self::parse_ion(ion) } } diff --git a/rust/rescore_dash/src/app.rs b/rust/rescore_dash/src/app.rs index 85468857..26e3a27d 100644 --- a/rust/rescore_dash/src/app.rs +++ b/rust/rescore_dash/src/app.rs @@ -294,7 +294,7 @@ fn available() -> bool { /// would unwind past the caller's warn-only `if let Err` -- after the results /// have already been written to disk. /// -/// The terminal is restored on every returning path. [`catch_panics`] covers +/// The terminal is restored on every returning path. `catch_panics` covers /// the event loop, but only in a debug build; see its doc. pub fn run(dash: Dashboard) -> std::io::Result<()> { if !available() { diff --git a/rust/rescore_dash/src/lib.rs b/rust/rescore_dash/src/lib.rs index e90fed8c..a747f7c4 100644 --- a/rust/rescore_dash/src/lib.rs +++ b/rust/rescore_dash/src/lib.rs @@ -2,7 +2,8 @@ //! target/decoy separation, FDR curve and decoy calibration. //! //! Two steps, deliberately separate. [`Dashboard::build`] materializes -//! everything on screen from a [`RescoreView`] -- see [`precompute`] for what is +//! everything on screen from a [`RescoreView`] -- see the `precompute` module +//! for what is //! exact and what is sampled -- and [`run`] opens the TUI over the result. //! Splitting them lets the caller drop the feature matrix, gigabytes at a //! realistic library size, before the TUI blocks for as long as the user leaves diff --git a/rust/speclib_build_cli/src/config.rs b/rust/speclib_build_cli/src/config.rs index b399c92c..7d301154 100644 --- a/rust/speclib_build_cli/src/config.rs +++ b/rust/speclib_build_cli/src/config.rs @@ -91,9 +91,11 @@ impl Default for DigestionConfig { #[derive(Debug, Deserialize)] #[serde(default, deny_unknown_fields)] pub struct ModificationsConfig { - /// Fixed modifications applied to every matching residue, e.g. ["Carbamidomethyl@C"]. + /// Fixed modifications applied to every matching residue, e.g. + /// `["Carbamidomethyl@C"]`. pub fixed: Vec, - /// Variable modifications considered during peptide generation, e.g. ["Oxidation@M"]. + /// Variable modifications considered during peptide generation, e.g. + /// `["Oxidation@M"]`. pub variable: Vec, /// Maximum number of variable modifications per peptide. pub max_variable: usize, diff --git a/rust/speclib_build_cli/src/entry.rs b/rust/speclib_build_cli/src/entry.rs index 68a7668e..3d253596 100644 --- a/rust/speclib_build_cli/src/entry.rs +++ b/rust/speclib_build_cli/src/entry.rs @@ -30,7 +30,7 @@ pub struct EntryFilters { /// Strip bracket-enclosed modifications from a sequence. /// -/// "PEPTC[U:4]IDEK" → "PEPTCIDEK" +/// `"PEPTC[U:4]IDEK"` -> `"PEPTCIDEK"` pub fn strip_mods(seq: &str) -> String { let mut out = String::with_capacity(seq.len()); let mut depth = 0usize; diff --git a/rust/speclib_build_cli/src/mods.rs b/rust/speclib_build_cli/src/mods.rs index a1c7d7cb..f7edfb69 100644 --- a/rust/speclib_build_cli/src/mods.rs +++ b/rust/speclib_build_cli/src/mods.rs @@ -1,6 +1,6 @@ /// Proforma-like modification parsing and application for speclib_build. /// -/// Supports notations like "C[U:4]", "M[U:35]", "S[U:21]", "M[+15.995]". +/// Supports notations like `C[U:4]`, `M[U:35]`, `S[U:21]`, `M[+15.995]`. /// Fixed mods are inserted after every matching residue in the sequence. /// Variable mods generate all combinations up to `max_mods` sites. @@ -10,7 +10,7 @@ pub struct Modification { /// Uppercase single-letter amino acid code. pub target_residue: char, - /// Bracket notation including brackets, e.g. "[U:4]" or "[+15.995]". + /// Bracket notation including brackets, e.g. `[U:4]` or `[+15.995]`. pub notation: String, } diff --git a/rust/timsquery/src/serde/diann_io.rs b/rust/timsquery/src/serde/diann_io.rs index 2fc2757b..b979348a 100644 --- a/rust/timsquery/src/serde/diann_io.rs +++ b/rust/timsquery/src/serde/diann_io.rs @@ -27,7 +27,7 @@ use tracing::{ pub enum DiannReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(DiannPrecursorParsingError), Parquet(String), Arrow(String), } @@ -37,8 +37,8 @@ impl std::fmt::Display for DiannReadingError { match self { DiannReadingError::Io => write!(f, "IO error"), DiannReadingError::Csv => write!(f, "CSV parsing error"), - DiannReadingError::PrecursorParsing => { - write!(f, "DIA-NN precursor parsing error") + DiannReadingError::PrecursorParsing(err) => { + write!(f, "DIA-NN precursor parsing error: {}", err) } DiannReadingError::Parquet(msg) => write!(f, "Parquet error: {}", msg), DiannReadingError::Arrow(msg) => write!(f, "Arrow error: {}", msg), @@ -50,7 +50,13 @@ impl std::error::Error for DiannReadingError {} #[derive(Debug)] pub enum DiannPrecursorParsingError { - IonParsingError, + /// Which fragment row failed, and why. The row index is the only handle the + /// user gets on a bad row: one failure aborts the whole library load, and + /// nothing upstream knows where it happened. + IonParsing { + row: usize, + source: IonParsingError, + }, /// A library that names its other precursors left this one blank. See /// [`Naming`]. UnnamedPrecursor, @@ -59,16 +65,34 @@ pub enum DiannPrecursorParsingError { Other, } -impl From for DiannPrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - DiannPrecursorParsingError::IonParsingError +impl std::fmt::Display for DiannPrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => { + write!(f, "fragment row {}: {}", row, source) + } + Self::UnnamedPrecursor => write!( + f, + "a library that names its other precursors left this one blank" + ), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::EmptyIonString => write!(f, "empty FragmentType"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl DiannPrecursorParsingError { + /// `map_err` adaptor that stamps the fragment row onto an ion-parsing + /// failure. Replaces a `From` impl, which had no way to see the row. + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for DiannReadingError { - fn from(_err: DiannPrecursorParsingError) -> Self { - DiannReadingError::PrecursorParsing + fn from(err: DiannPrecursorParsingError) -> Self { + DiannReadingError::PrecursorParsing(err) } } @@ -434,7 +458,9 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -457,7 +483,8 @@ fn parse_precursor_group( DiannPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); @@ -765,7 +792,9 @@ fn parse_precursor_group_from_parquet( columns.fragment_loss_types[idx], i ); - let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); rel_intensities.push((ion_annot, rel_intensity)); @@ -788,7 +817,8 @@ fn parse_precursor_group_from_parquet( DiannPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); @@ -867,7 +897,9 @@ AAAAAAALQAK\tAAAAAAALQAK\t478.7\t2\t11.0\t0.9\tP2\t0\t300.0\ty\t3\t1\tnoloss\t1. assert!( matches!( read_targets(blank.path()), - Err(DiannReadingError::PrecursorParsing) + Err(DiannReadingError::PrecursorParsing( + DiannPrecursorParsingError::UnnamedPrecursor + )) ), "a blank name in a naming library must not load" ); diff --git a/rust/timsquery/src/serde/diann_speclib_io.rs b/rust/timsquery/src/serde/diann_speclib_io.rs index 1ab6c73e..34a58f8a 100644 --- a/rust/timsquery/src/serde/diann_speclib_io.rs +++ b/rust/timsquery/src/serde/diann_speclib_io.rs @@ -8,12 +8,13 @@ //! (little-endian, no padding, no section table); it must be parsed strictly in //! order, so there is no random access across sections. //! -//! The reader is three layers: -//! 1. Decode ([`Cursor`], [`Fragment`], [`Peptide`]): typed, zero-copy views +//! The reader is three layers. The types named here are private to this module, +//! so they are spelled rather than linked: +//! 1. Decode (`Cursor`, `Fragment`, `Peptide`): typed, zero-copy views //! over byte ranges of the in-memory file. No domain knowledge. -//! 2. Emit ([`SpecLib`] + [`EntryIter`]): a pull parser that walks the -//! variable-length envelope and yields one [`EntryView`] per precursor. -//! 3. Map ([`map_entry`]): [`EntryView`] -> rows pushed directly into a +//! 2. Emit (`SpecLib` + `EntryIter`): a pull parser that walks the +//! variable-length envelope and yields one `EntryView` per precursor. +//! 3. Map (`map_entry`): `EntryView` -> rows pushed directly into a //! columnar `TargetColumns` (plus a parallel reference-intensity //! sidecar). This is where the ion-type table, dedup, and drop stats live. //! @@ -569,9 +570,11 @@ fn bounded_capacity(count: usize, remaining_bytes: usize, min_record_bytes: usiz #[derive(Debug, Default, Clone, Copy)] pub struct SpeclibDecodeStats { /// Fragments flagged `ExcludeFromAssay` (`type & 0x80`). Kept (see - /// [`map_entry`]); counted for reporting only. + /// `map_entry`); counted for reporting only. pub exclude_flagged: usize, - /// Fragments with a neutral loss (`loss != 0`) -- `IonAnnot` cannot hold loss. + /// Fragments with a neutral loss (`loss != 0`). `IonAnnot` can represent + /// these; what is missing is the DIA-NN loss-code mapping (see #105), so + /// they are dropped and counted to keep the cost visible. pub loss_dropped: usize, /// Fragments whose ion-type code or recovered series was unusable. pub unknown_ion_dropped: usize, @@ -809,11 +812,29 @@ fn map_entry( continue; } - let ion = match IonAnnot::try_new(type_char, Some(series as u8), f.charge() as i8, 0) { + // `charge` is a raw wire byte. `as i8` would reinterpret anything above + // 127 as negative and hand a plausible-looking charge to the + // constructor, so it is range-checked instead. Note the sibling `typ()` + // is masked `& 0x7F` because DIA-NN sets a flag bit there; if this byte + // carries flags too, this is where it will surface (see #105). + let charge = match i8::try_from(f.charge()) { + Ok(charge) => charge, + Err(_) => { + warn!( + "speclib entry {:?}: fragment charge byte {} is not a charge; dropping", + name, + f.charge() + ); + stats.unknown_ion_dropped += 1; + continue; + } + }; + + let ion = match IonAnnot::try_new(type_char, Some(series as u8), charge, 0) { Ok(ion) => ion, Err(e) => { warn!( - "speclib entry {:?}: failed to build IonAnnot ({:?}); dropping fragment", + "speclib entry {:?}: failed to build IonAnnot ({}); dropping fragment", name, e ); stats.unknown_ion_dropped += 1; @@ -874,7 +895,7 @@ type ParsedSpeclib = (TargetColumns, Vec, SpeclibDecodeStats, boo /// structural desync. /// /// Reads the file, parses the header, then parses+maps entries in parallel -/// ([`SpecLib::open`] + [`SpecLib::parse_parallel`]). +/// (`SpecLib::open` + `SpecLib::parse_parallel`). pub fn parse_speclib_reader(reader: R) -> Result { SpecLib::open(reader)?.parse_parallel() } diff --git a/rust/timsquery/src/serde/skyline_io.rs b/rust/timsquery/src/serde/skyline_io.rs index 84e966de..d4630390 100644 --- a/rust/timsquery/src/serde/skyline_io.rs +++ b/rust/timsquery/src/serde/skyline_io.rs @@ -21,7 +21,7 @@ use tracing::{ pub enum SkylineReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(SkylinePrecursorParsingError), } #[derive(Debug)] @@ -37,8 +37,8 @@ impl std::fmt::Display for SkylineReadingError { match self { SkylineReadingError::Io => write!(f, "IO error"), SkylineReadingError::Csv => write!(f, "CSV parsing error"), - SkylineReadingError::PrecursorParsing => { - write!(f, "Skyline precursor parsing error") + SkylineReadingError::PrecursorParsing(err) => { + write!(f, "Skyline precursor parsing error: {}", err) } } } @@ -48,21 +48,37 @@ impl std::error::Error for SkylineReadingError {} #[derive(Debug)] pub enum SkylinePrecursorParsingError { - IonParsingError, + /// Which fragment row failed, and why. One failure aborts the whole library + /// load, and nothing upstream knows where it happened. + IonParsing { + row: usize, + source: IonParsingError, + }, IonOverCapacity, Other, } -impl From for SkylinePrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - SkylinePrecursorParsingError::IonParsingError +impl std::fmt::Display for SkylinePrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => write!(f, "fragment row {}: {}", row, source), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl SkylinePrecursorParsingError { + /// `map_err` adaptor that stamps the fragment row onto an ion-parsing + /// failure. Replaces a `From` impl, which had no way to see the row. + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for SkylineReadingError { - fn from(_err: SkylinePrecursorParsingError) -> Self { - SkylineReadingError::PrecursorParsing + fn from(err: SkylinePrecursorParsingError) -> Self { + SkylineReadingError::PrecursorParsing(err) } } @@ -341,7 +357,8 @@ fn parse_precursor_group( ); SkylinePrecursorParsingError::IonOverCapacity })?; - IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)? + IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(SkylinePrecursorParsingError::ion(i))? } None => { warn!( @@ -349,7 +366,9 @@ fn parse_precursor_group( falling back to unknown ion", frag_type, i ); - unknown_ions.next_unknown(frag_charge as i8)? + unknown_ions + .next_unknown(frag_charge as i8) + .map_err(SkylinePrecursorParsingError::ion(i))? } }; diff --git a/rust/timsquery/src/serde/spectronaut_io.rs b/rust/timsquery/src/serde/spectronaut_io.rs index f12a0ba3..0b3ccf58 100644 --- a/rust/timsquery/src/serde/spectronaut_io.rs +++ b/rust/timsquery/src/serde/spectronaut_io.rs @@ -17,7 +17,7 @@ use tracing::{ pub enum SpectronautReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(SpectronautPrecursorParsingError), } /// Error type for library format detection (sniffing) @@ -36,8 +36,8 @@ impl std::fmt::Display for SpectronautReadingError { match self { SpectronautReadingError::Io => write!(f, "IO error"), SpectronautReadingError::Csv => write!(f, "CSV parsing error"), - SpectronautReadingError::PrecursorParsing => { - write!(f, "Spectronaut precursor parsing error") + SpectronautReadingError::PrecursorParsing(err) => { + write!(f, "Spectronaut precursor parsing error: {}", err) } } } @@ -47,22 +47,39 @@ impl std::error::Error for SpectronautReadingError {} #[derive(Debug)] pub enum SpectronautPrecursorParsingError { - IonParsingError, + /// Which fragment row failed, and why. One failure aborts the whole library + /// load, and nothing upstream knows where it happened. + IonParsing { + row: usize, + source: IonParsingError, + }, IonOverCapacity, EmptyIonString, Other, } -impl From for SpectronautPrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - SpectronautPrecursorParsingError::IonParsingError +impl std::fmt::Display for SpectronautPrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => write!(f, "fragment row {}: {}", row, source), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::EmptyIonString => write!(f, "empty FragmentType"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl SpectronautPrecursorParsingError { + /// `map_err` adaptor that stamps the fragment row onto an ion-parsing + /// failure. Replaces a `From` impl, which had no way to see the row. + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for SpectronautReadingError { - fn from(_err: SpectronautPrecursorParsingError) -> Self { - SpectronautReadingError::PrecursorParsing + fn from(err: SpectronautPrecursorParsingError) -> Self { + SpectronautReadingError::PrecursorParsing(err) } } @@ -304,7 +321,9 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - let ion_annot = unknown_ions.next_unknown(frag_charge as i8)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(SpectronautPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -327,7 +346,8 @@ fn parse_precursor_group( SpectronautPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(SpectronautPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); diff --git a/rust/timsquery/src/utils/sorting.rs b/rust/timsquery/src/utils/sorting.rs index 4bac6dea..d9ca8e04 100644 --- a/rust/timsquery/src/utils/sorting.rs +++ b/rust/timsquery/src/utils/sorting.rs @@ -1,48 +1,3 @@ -/// Macro that sorts an arbitrary number of vecs by a the values -/// first one. -/// -/// NOTE: This macro creates a new ordered vec for each one. -/// In theory its possible to have this happen in-place (see commit history) -/// but it seems ineficient for the most part when I benchmarked it. -/// -/// TODO: Make a variant that sort in parallel. -/// TODO: Add a parameter to specify ascending or descending. -/// -/// # Example -/// ``` -/// use timsquery::sort_vecs_by_first; -/// -/// let va = vec![9, 8, 7]; -/// let vb = vec![1, 2, 3]; -/// let vc = vec!['a', 'b', 'c']; -/// let out = sort_vecs_by_first!(&va, &vb, &vc); -/// -/// assert_eq!(out.0, vec![7, 8, 9]); -/// assert_eq!(out.1, vec![3, 2, 1]); -/// assert_eq!(out.2, vec!['c', 'b', 'a']); -/// ``` -#[macro_export] -macro_rules! sort_vecs_by_first { - ($first:expr $(,$rest:expr)*) => {{ - let first_vec = $first; - let len = first_vec.len(); - - // Create and sort indices - let mut indices: Vec<_> = (0..len).collect(); - indices.sort_unstable_by_key(|&i| &first_vec[i]); - - // Reorder first vector - let sorted_first: Vec<_> = indices.iter().map(|&i| first_vec[i]).collect(); - - // Reorder all other vectors - (sorted_first, $( { - let other_vec = $rest; - assert_eq!(other_vec.len(), len, "All vectors must have the same length"); - indices.iter().map(|&i| other_vec[i]).collect::>() - }, )*) - }}; -} - /// Returns the top n elements of a slice. /// /// The indices of the elements are also returned. @@ -94,45 +49,6 @@ pub fn top_n(slice: &[T], n: usize) -> (Vec, Vec mod tests { use super::*; - #[test] - fn test_sort_two_vecs() { - let v1 = vec![3, 1, 4, 1, 5]; - let v2 = vec!['a', 'b', 'c', 'd', 'e']; - - let (sorted_v1, sorted_v2) = sort_vecs_by_first!(v1, v2); - - assert_eq!(sorted_v1, vec![1, 1, 3, 4, 5]); - assert_eq!(sorted_v2, vec!['b', 'd', 'a', 'c', 'e']); - } - - #[test] - fn test_sort_three_vecs() { - let v1 = vec![3, 1, 4]; - let v2 = vec!['x', 'y', 'z']; - let v3 = vec![true, false, true]; - - let (sorted_v1, sorted_v2, sorted_v3) = sort_vecs_by_first!(v1, v2, v3); - - assert_eq!(sorted_v1, vec![1, 3, 4]); - assert_eq!(sorted_v2, vec!['y', 'x', 'z']); - assert_eq!(sorted_v3, vec![false, true, true]); - } - - #[test] - fn test_sort_four_vecs() { - let v1 = vec![3, 1, 4]; - let v2 = vec!['x', 'y', 'z']; - let v3 = vec![true, false, true]; - let v4 = vec![1.0, 2.0, 3.0]; - - let (sorted_v1, sorted_v2, sorted_v3, sorted_v4) = sort_vecs_by_first!(v1, v2, v3, v4); - - assert_eq!(sorted_v1, vec![1, 3, 4]); - assert_eq!(sorted_v2, vec!['y', 'x', 'z']); - assert_eq!(sorted_v3, vec![false, true, true]); - assert_eq!(sorted_v4, vec![2.0, 1.0, 3.0]); - } - #[test] fn test_top_n() { let v = vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; diff --git a/rust/timsquery_viewer/src/calibration.rs b/rust/timsquery_viewer/src/calibration.rs index 392e9c2e..729500c6 100644 --- a/rust/timsquery_viewer/src/calibration.rs +++ b/rust/timsquery_viewer/src/calibration.rs @@ -653,7 +653,7 @@ impl ViewerCalibrationState { /// /// `indexed_data` and `elution_groups` are needed to enable the Start /// button (we need both loaded). `tolerance` is written when the user - /// clicks [Apply]. + /// clicks `Apply`. pub fn render_panel( &mut self, ui: &mut egui::Ui, diff --git a/rust/timsseek/src/data_sources/speclib.rs b/rust/timsseek/src/data_sources/speclib.rs index 1905ea8d..45909066 100644 --- a/rust/timsseek/src/data_sources/speclib.rs +++ b/rust/timsseek/src/data_sources/speclib.rs @@ -176,7 +176,7 @@ fn strip_mods(s: &str) -> String { out } -/// Summary of a [`finalize_reference_library`] call, for load-time logging. +/// Summary of a `finalize_reference_library` call, for load-time logging. #[derive(Debug, Clone, Copy)] pub struct LoadReport { /// Physical stored rows (pre decoy expansion), i.e. `TargetColumns::n_rows`. diff --git a/rust/timsseek/src/ml/cv.rs b/rust/timsseek/src/ml/cv.rs index 8ea27621..cbd24006 100644 --- a/rust/timsseek/src/ml/cv.rs +++ b/rust/timsseek/src/ml/cv.rs @@ -612,7 +612,7 @@ pub(crate) fn fold_weights(data: &D, rows: &[usize]) -> Vec /// [`FoldModel`] adapter for `forust`'s [`GradientBooster`]. A newtype only /// because `GradientBooster` is a foreign type; it also carries the lane width -/// so [`FoldModel::importance`] can return a full-width, lane-indexed vector +/// so `FoldModel::importance` can return a full-width, lane-indexed vector /// (forust reports only the columns it split on). pub(crate) struct GbmFoldModel { booster: GradientBooster, @@ -667,7 +667,7 @@ impl FoldModel for GbmFoldModel { /// forust reports only the columns it actually split on. Every other lane /// is `NAN` -- "not reported" -- rather than `0.0`, because this model has no /// gain measurement for a column it never used. See the - /// [`FoldModel::importance`] contract. + /// `FoldModel::importance` contract. fn importance(&self) -> Vec { let raw = self .booster @@ -755,10 +755,10 @@ pub struct FeatureStat { /// Per-fold feature statistics. /// /// `feature_stats` is in the dataset's own column order -/// ([`FoldDataset::column_names`], i.e. the matrix's column order), one entry per +/// (`FoldDataset::column_names`, i.e. the matrix's column order), one entry per /// column. `feature_importance` is sorted by importance DESCENDING (top features /// first), and carries only the columns the model reported a finite value for -- -/// see [`FoldModel::importance`] -- so it is generally shorter than +/// see `FoldModel::importance` -- so it is generally shorter than /// `feature_stats` and in a different order. /// /// Descending because the two TSV sidecars `timsseek_cli` writes @@ -798,7 +798,7 @@ pub(crate) fn fold_feature_stats( for (fold, rows) in fold_rows.iter().enumerate() { // --- Importance, back to the (name, gain) sidecar shape --- // UNREPORTED (`NAN`) columns are dropped; every REPORTED column is - // emitted, `0.0` included. See the [`FoldModel::importance`] contract + // emitted, `0.0` included. See the `FoldModel::importance` contract // for why those are different things. The drop matters because the // sidecar and the dashboard's fold-averaged gain both treat a feature // as "reported by this fold" simply by being present, so a model's @@ -1472,7 +1472,7 @@ mod test { } } - /// THE [`FoldModel::importance`] sentinel contract, at the sidecar boundary: + /// THE `FoldModel::importance` sentinel contract, at the sidecar boundary: /// `NAN` means "this model reports nothing for this column" and is dropped; /// every FINITE value reaches the sidecar, `0.0` included. /// diff --git a/rust/timsseek/src/ml/qvalues.rs b/rust/timsseek/src/ml/qvalues.rs index 0fefd4a9..e16a08fe 100644 --- a/rust/timsseek/src/ml/qvalues.rs +++ b/rust/timsseek/src/ml/qvalues.rs @@ -230,7 +230,7 @@ fn finalize( /// A cross-fit (leak-free) model over the canonical rescore fold partition. /// -/// See [`crossfit`] for the partition contract. +/// See `crossfit` for the partition contract. struct CrossFit { /// Held-out score per row, row-aligned with the input matrix. scores: Vec, @@ -259,7 +259,7 @@ impl CrossFit { /// Cross-fit a [`FoldModel`] over the canonical rescore fold partition and /// return each row's HELD-OUT score. /// -/// The shared statement of leak-freedom for [`crossfit_lda`], `rescore_lda`, +/// The shared statement of leak-freedom for `crossfit_lda`, `rescore_lda`, /// and `rescore_hybrid`. /// /// Generic over the model because the partition is independent of the fitted @@ -380,7 +380,7 @@ where }) } -/// Cross-fit an [`LdaModel`] on [`LdaConfig::default`] -- see [`crossfit`] for +/// Cross-fit an [`LdaModel`] on [`LdaConfig::default`] -- see `crossfit` for /// the partition and the failure policy. /// /// `N_RESCORE_FOLDS` fits instead of 1 is affordable precisely because the LDA @@ -451,10 +451,10 @@ pub fn rescore(mut data: Vec) -> RescoreResult { /// only the discriminant score source changes. /// /// CROSS-FIT, not a single in-sample fit: every row's score comes from an LDA -/// fitted without that row, via [`crossfit_lda`] -- see [`crossfit`] for the +/// fitted without that row, via `crossfit_lda` -- see `crossfit` for the /// partition and why it is mandatory. /// -/// Returns PER-FOLD [`FoldStats`] (one per fold, like the GBM path): feature +/// Returns PER-FOLD `FoldStats` (one per fold, like the GBM path): feature /// means/NaN ratios over each fold's held-out rows, `|coef|` importance from /// that fold's model. /// @@ -502,7 +502,7 @@ pub fn rescore_lda(mut data: Vec) -> RescoreResult { /// The fold count is `N_RESCORE_FOLDS`, i.e. the SAME /// `get_fold` the GBM's [`CrossValidatedScorer`] derives its partition from /// below. That shared definition is what makes "the same fold assignment on both -/// sides" structural rather than a comment; see [`crossfit`]. +/// sides" structural rather than a comment; see `crossfit`. fn hybrid_linear_dataset(data: &[CompetedCandidate]) -> StreamingDataset<'_, CompetedCandidate> { StreamingDataset::new( data, @@ -559,7 +559,7 @@ fn hybrid_frame( /// `lda_score` as one extra column into the NONLINEAR lane, then train the GBM /// CV on `nonlinear + lda_score` instead of the full feature frame. /// -/// LEAK-FREEDOM: `lda_score` is cross-fit via [`crossfit`] -- see there for the +/// LEAK-FREEDOM: `lda_score` is cross-fit via `crossfit` -- see there for the /// partition, why a label-aware feature fed to a CV'd GBM in particular must be /// leak-free, and why the fold ASSIGNMENT has to match the one /// `CrossValidatedScorer` derives its own partition from. ASSIGNMENT, not @@ -645,7 +645,7 @@ fn rescore_mlp_with(mut data: Vec, config: MlpConfig) -> Resc /// [`rescore`] gives the GBM, so the two are directly comparable. The one /// [`crate::ml::RescoreModel::Mlp`] selects. /// -/// See [`rescore_mlp_with`] for the cross-fit and determinism contracts. +/// See `rescore_mlp_with` for the cross-fit and determinism contracts. /// /// Runtime and sensitivity comparisons are not constant across candidate counts. pub fn rescore_mlp(data: Vec) -> RescoreResult { @@ -817,7 +817,7 @@ fn build_all_matrix<'a>( out } -/// Competed candidates in the shape [`build_all_matrix`] consumes. +/// Competed candidates in the shape `build_all_matrix` consumes. fn competed_rows( data: &[CompetedCandidate], ) -> impl ExactSizeIterator { @@ -833,7 +833,7 @@ pub fn feature_frame(data: &[FinalResult]) -> (Vec>, Vec) { (all_feature_name_set(), build_all_matrix(rows)) } -/// LINEAR-lane feature names (LDA), in [`project_linear_row`]'s order. +/// LINEAR-lane feature names (LDA), in `project_linear_row`'s order. pub fn linear_feature_name_set() -> Vec> { let mut n = NameSink::new(); ::linear_feature_names(&mut n); @@ -842,7 +842,7 @@ pub fn linear_feature_name_set() -> Vec> { n.into_names() } -/// NONLINEAR-lane feature names, in [`project_nonlinear_row`]'s order. The +/// NONLINEAR-lane feature names, in `project_nonlinear_row`'s order. The /// `sequence_counts` names are unconditional -- a peptide with no parsed /// sequence contributes NaN values under them, not a shorter row. pub fn nonlinear_feature_name_set() -> Vec> { @@ -855,7 +855,7 @@ pub fn nonlinear_feature_name_set() -> Vec> { } /// The ALL-lane feature names (GBM) = linear ++ nonlinear, matching -/// [`build_all_matrix`]'s column order. +/// `build_all_matrix`'s column order. pub fn all_feature_name_set() -> Vec> { let mut v = linear_feature_name_set(); v.extend(nonlinear_feature_name_set()); diff --git a/rust/timsseek/src/utils/elution_group_ops.rs b/rust/timsseek/src/utils/elution_group_ops.rs index dbd8f3e2..44fea8e7 100644 --- a/rust/timsseek/src/utils/elution_group_ops.rs +++ b/rust/timsseek/src/utils/elution_group_ops.rs @@ -16,9 +16,12 @@ pub fn apply_isotope_offset_fragments_into( ) { dst.reset_from(src); for (k, v) in dst.iter_fragments_refs_mut() { - let new_ions = k.try_with_offset_neutrons(offset).expect( - "Isotope offset overflow - this should never happen with realistic isotope offsets", - ); + // The error names the representable range, so a library carrying an + // isotope offset too close to the ceiling says what the limit is rather + // than asserting the situation is impossible. + let new_ions = k.try_with_offset_neutrons(offset).unwrap_or_else(|e| { + panic!("fragment {k} cannot take a {offset:+} isotope offset: {e}") + }); let mz_offset = (C13_C12_MASS_DIFF / k.get_charge() as f64) * offset as f64; *v += mz_offset; *k = new_ions; From f908845a10266f938a1c3b018a0f4a9d94dbc50a Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 12:06:56 -0700 Subject: [PATCH 5/8] refactor(micromzpaf): give each inverse pair one home `Display for IonSeriesOrdinal` emitted `m2:11` and `IA`, but the code parsing those forms was inlined in `IonAnnot::parse_ion` -- an inverse pair split across a type boundary. Likewise `INTERNAL_POS_MAX` was checked in `IonAnnot`, 400 lines from the 6-bit packing that makes 63 the right bound; the `debug_assert` in `pack` was a debug-only backstop for an invariant two unrelated functions maintained by hand. `series` now owns all three: the spelling (`Display` next to `parse`), the payload layout (`to_parts` next to `from_parts`), and the bounds that layout imposes (`try_internal`, `try_immonium`). It exports `KIND_BITS`/`PAYLOAD_BITS` so the word can be assembled without knowing how they are filled, and asserts its own payload bounds. `to_parts` is total now -- every arm masks to its own width, so the immonium `- b'A'` cannot underflow for a variant built by hand. `IonAnnot` keeps two constructors instead of four: `new(series, loss, charge, isotope)`, and `try_new(char, ordinal, charge, isotope)` for the shape a file reader actually has. `try_new_with_loss`, `try_new_internal` and `try_new_immonium` are gone. `parse_ion` is now the three suffixes peeled in the one order that is unambiguous, then a delegation, and the module doc draws which delimiter belongs to which stage -- previously that ordering was load-bearing and only the `-` case said so. `lib.rs` drops 1210 -> 705 lines, with `error.rs` and `parse.rs` taking the rest. Tests moved with their subjects: the packing and spelling round trips are in `series`, where they can see the encoding they are pinning. The constructors are `try_internal`/`try_immonium` rather than `internal`/`immonium` because a constructor sharing a name with its variant is ambiguous to rustdoc, and reads ambiguously to a person too. --- rust/micromzpaf/src/error.rs | 60 +++ rust/micromzpaf/src/lib.rs | 665 ++++------------------------------ rust/micromzpaf/src/parse.rs | 144 ++++++++ rust/micromzpaf/src/series.rs | 458 +++++++++++++++++++++++ 4 files changed, 742 insertions(+), 585 deletions(-) create mode 100644 rust/micromzpaf/src/error.rs create mode 100644 rust/micromzpaf/src/parse.rs create mode 100644 rust/micromzpaf/src/series.rs diff --git a/rust/micromzpaf/src/error.rs b/rust/micromzpaf/src/error.rs new file mode 100644 index 00000000..052d0035 --- /dev/null +++ b/rust/micromzpaf/src/error.rs @@ -0,0 +1,60 @@ +//! Why an annotation could not be parsed or represented. + +use thiserror::Error; + +use crate::{ + CHARGE_MAX, + CHARGE_MIN, + ISOTOPE_MAX, + ISOTOPE_MIN, +}; + +/// Why an annotation could not be parsed or represented. +/// +/// Never matched outside this crate today, so the variants exist for the message +/// a user sees when a library fails to load. Each one names the offending value +/// and, where there is one, the bound it missed -- a reader surfaces this with +/// no other context. +#[derive(Debug, Error)] +pub enum IonParsingError { + #[error("Ordinal {ordinal} out of range for series '{series}'")] + OrdinalOutOfRange { ordinal: u8, series: char }, + #[error("Series '{series}' requires an ordinal")] + MissingOrdinal { series: char }, + #[error("Series '{series}' takes no ordinal, got {ordinal}")] + UnexpectedOrdinal { series: char, ordinal: u8 }, + #[error("Unsupported fragment type: '{fragment_type}'")] + UnsupportedFragmentType { fragment_type: char }, + #[error("Charge cannot be 0")] + ChargeCannotBeZero, + #[error("Charge {charge} outside the representable range {CHARGE_MIN}..={CHARGE_MAX}")] + ChargeOutOfRange { charge: i8 }, + #[error( + "Isotope offset {isotope} outside the representable range {ISOTOPE_MIN}..={ISOTOPE_MAX}" + )] + IsotopeOutOfRange { isotope: i8 }, + #[error("Neutral loss '{loss}' is not representable")] + UnsupportedNeutralLoss { loss: String }, + #[error("Immonium ions must be a bare uppercase residue, got '{annotation}'")] + UnsupportedModifiedImmonium { annotation: String }, + #[error("Ran out of distinct unknown-ion labels: the 8-bit ordinal is exhausted")] + UnknownIonsExhausted, + #[error("Could not parse '{annotation}': {context}")] + ParsingError { + /// The whole annotation, not the fragment of it that failed: a reader + /// reports this with no row index, so the full text is the only handle + /// the user gets on which row broke. + annotation: String, + context: &'static str, + }, +} + +impl IonParsingError { + /// Shorthand for [`Self::ParsingError`], which is built at a dozen sites. + pub(crate) fn parse(annotation: &str, context: &'static str) -> Self { + Self::ParsingError { + annotation: annotation.to_string(), + context, + } + } +} diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 858edf1a..67b74833 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -83,21 +83,40 @@ //! assert!(IonAnnot::try_from("IC[Carbamidomethyl]").is_err()); //! ``` +pub mod error; pub mod loss; +pub mod parse; +pub mod series; +pub use error::IonParsingError; pub use loss::NeutralLoss; +pub use parse::{ + MassError, + split_mass_error, +}; +pub use series::{ + INTERNAL_POS_MAX, + IonSeriesOrdinal, + Series, +}; + use serde::{ Deserialize, Serialize, }; +use series::{ + KIND_BITS, + PAYLOAD_BITS, +}; use std::fmt::Display; use std::hash::Hash; -use thiserror::Error; -// ── Bit layout ─────────────────────────────────────────────────────────────── +// ── The word ───────────────────────────────────────────────────────────────── +// +// `kind` and `payload` are `series`'s to define; this module only decides where +// they sit and what surrounds them. const KIND_SHIFT: u32 = 0; -const KIND_BITS: u32 = 4; const CHARGE_SHIFT: u32 = 4; const CHARGE_BITS: u32 = 4; const ISOTOPE_SHIFT: u32 = 8; @@ -105,7 +124,6 @@ const ISOTOPE_BITS: u32 = 4; const LOSS_SHIFT: u32 = 12; const LOSS_BITS: u32 = 6; const PAYLOAD_SHIFT: u32 = 18; -const PAYLOAD_BITS: u32 = 12; /// Charge bounds. Observed maximum in the HUPO-PSI corpus is 3. /// @@ -125,15 +143,10 @@ pub const CHARGE_MAX: i8 = 8; /// unreachable and spends the whole 4 bits on offsets that can round-trip. pub const ISOTOPE_MIN: i8 = 0; pub const ISOTOPE_MAX: i8 = mask(ISOTOPE_BITS) as i8; -/// Width of each internal-fragment endpoint inside `payload`. -const INTERNAL_POS_BITS: u32 = 6; -/// Width of the immonium residue index inside `payload`. -const IMMONIUM_BITS: u32 = 5; -/// Widest residue index an internal fragment endpoint holds. -pub const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; // The layout is otherwise enforced by prose and a diagram. These make it -// self-checking, so a field width cannot be widened without the build failing. +// self-checking, so a field cannot be widened or moved without the build +// failing. `series` asserts its own payload bounds. const _: () = assert!( KIND_BITS + CHARGE_BITS + ISOTOPE_BITS + LOSS_BITS + PAYLOAD_BITS <= 32, "the fields overflow the word" @@ -149,14 +162,6 @@ const _: () = assert!( "fields must abut" ); const _: () = assert!(PAYLOAD_SHIFT == LOSS_SHIFT + LOSS_BITS, "fields must abut"); -const _: () = assert!( - 2 * INTERNAL_POS_BITS <= PAYLOAD_BITS, - "two internal-fragment endpoints must fit in one payload" -); -const _: () = assert!( - mask(IMMONIUM_BITS) >= (b'Z' - b'A') as u32, - "the immonium field must hold every uppercase residue" -); // The zigzag bounds are hand-derived; these are what make them checked. const _: () = assert!( zigzag_charge(CHARGE_MIN) <= mask(CHARGE_BITS) @@ -175,9 +180,10 @@ const _: () = assert!( ); #[inline] -const fn mask(bits: u32) -> u32 { +pub(crate) const fn mask(bits: u32) -> u32 { (1u32 << bits) - 1 } + /// Zigzag: map a small signed value onto an unsigned one without losing the /// sign bit to the field width. #[inline] @@ -194,7 +200,6 @@ const fn unzigzag(u: u32) -> i8 { /// Charge 0 is rejected by every constructor, so it is not a value the field /// needs to represent -- and spending the zero word on it would make /// `IonAnnot::default()` render an annotation that cannot be parsed back. -/// `CHARGE_MIN..=CHARGE_MAX` minus one still zigzags inside 4 bits. #[inline] const fn zigzag_charge(charge: i8) -> u32 { zigzag(charge - 1) @@ -224,15 +229,10 @@ impl Serialize for IonAnnot { } } -/// Deserializes simple annotations for fragments. -/// -/// b12^3 -> b12 charge 3 (implicit 0 isotope) -/// b12+i^3 -> b12 charge 3 isotope 1 -/// b12+3i^3 -> b12 charge 3 isotope 2 -/// b13 -> b13 (implicit charge 1 and isotope 0) +/// Deserializes an annotation from its mzPAF string. /// -/// The wire format is the mzPAF string, not the packed word, so the bit layout -/// can change without breaking existing files. +/// The wire format is the string, not the packed word, so the bit layout can +/// change without breaking existing files. impl<'de> Deserialize<'de> for IonAnnot { fn deserialize(deserializer: D) -> Result where @@ -244,88 +244,14 @@ impl<'de> Deserialize<'de> for IonAnnot { } impl IonAnnot { - /// Build a backbone / precursor / unknown annotation. - /// - /// Errors when `charge` is zero or either of `charge`/`isotope` falls - /// outside the packed field's range. That range check is load-bearing: the - /// bit field truncates silently, so an unchecked value would corrupt the - /// annotation rather than fail. - pub fn try_new( - ion_type: char, - ordinal: Option, - charge: i8, - isotope: i8, - ) -> Result { - Self::try_new_with_loss(ion_type, ordinal, charge, isotope, NeutralLoss::None) - } - - /// As [`Self::try_new`], carrying a neutral loss. - pub fn try_new_with_loss( - ion_type: char, - ordinal: Option, - charge: i8, - isotope: i8, - loss: NeutralLoss, - ) -> Result { - Self::pack( - IonSeriesOrdinal::from_series_char(ion_type, ordinal)?, - loss, - charge, - isotope, - ) - } - - /// Build an internal fragment spanning residues `start..=end`. - /// - /// Endpoints are capped at [`INTERNAL_POS_MAX`], narrower than a backbone - /// ordinal, because two of them share the 12-bit payload. Internal - /// fragments are bounded by peptide length, so 63 is well past tryptic. - pub fn try_new_internal( - start: u8, - end: u8, - charge: i8, - isotope: i8, - loss: NeutralLoss, - ) -> Result { - if start > INTERNAL_POS_MAX || end > INTERNAL_POS_MAX { - return Err(IonParsingError::OrdinalOutOfRange { - ordinal: start.max(end), - series: 'm', - }); - } - Self::pack( - IonSeriesOrdinal::internal { start, end }, - loss, - charge, - isotope, - ) - } - - /// Build a bare immonium ion for an uppercase residue code. + /// Build an annotation from its parts. /// - /// Modified immonium (`IC[Carbamidomethyl]`) is not representable at any - /// field width and is rejected by the parser rather than silently losing - /// the modification. - pub fn try_new_immonium( - residue: char, - charge: i8, - isotope: i8, - loss: NeutralLoss, - ) -> Result { - if !residue.is_ascii_uppercase() { - return Err(IonParsingError::UnsupportedFragmentType { - fragment_type: residue, - }); - } - Self::pack( - IonSeriesOrdinal::immonium { residue }, - loss, - charge, - isotope, - ) - } - - fn pack( + /// The one constructor. `series` has already been validated against the + /// payload width by [`IonSeriesOrdinal`]'s own constructors, so all this + /// checks is `charge` and `isotope` -- and it must, because a bit field + /// truncates silently, so an unchecked value would corrupt the annotation + /// rather than fail. + pub fn new( series: IonSeriesOrdinal, loss: NeutralLoss, charge: i8, @@ -341,7 +267,6 @@ impl IonAnnot { return Err(IonParsingError::IsotopeOutOfRange { isotope }); } let (kind, payload) = series.to_parts(); - debug_assert!(payload <= mask(PAYLOAD_BITS), "payload overflows its field"); Ok(IonAnnot( (kind << KIND_SHIFT) | ((zigzag_charge(charge) & mask(CHARGE_BITS)) << CHARGE_SHIFT) @@ -351,6 +276,26 @@ impl IonAnnot { )) } + /// A backbone / precursor / unknown annotation from its mzPAF letter. + /// + /// The shape a file reader has: a series character straight out of a column. + /// Internal fragments and immonium ions are spelled differently and go + /// through [`IonSeriesOrdinal::try_internal`] / [`IonSeriesOrdinal::try_immonium`] + /// and [`Self::new`]. + pub fn try_new( + ion_type: char, + ordinal: Option, + charge: i8, + isotope: i8, + ) -> Result { + Self::new( + IonSeriesOrdinal::from_series_char(ion_type, ordinal)?, + NeutralLoss::None, + charge, + isotope, + ) + } + #[inline] fn payload(self) -> u32 { (self.0 >> PAYLOAD_SHIFT) & mask(PAYLOAD_BITS) @@ -375,8 +320,8 @@ impl IonAnnot { /// Shift the isotope by `offset_neutrons`. /// - /// Errors when the result leaves [`ISOTOPE_MIN`]..=[`ISOTOPE_MAX`]. That - /// bound is an order of magnitude past any observed isotope offset. + /// Errors when the result leaves [`ISOTOPE_MIN`]..=[`ISOTOPE_MAX`], naming + /// that range, since the caller cannot see the field width. pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { // Saturating rather than checked: the field is far narrower than `i8`, // so the range check below is the one that matters and it reports the @@ -416,161 +361,12 @@ impl IonAnnot { } } - /// The logical series-and-payload view of this annotation. + /// What kind of ion this is, and whatever that kind carries. pub fn series_ordinal(&self) -> IonSeriesOrdinal { IonSeriesOrdinal::from_parts((self.0 >> KIND_SHIFT) & mask(KIND_BITS), self.payload()) } } -/// A mass-error suffix on an mzPAF annotation: observed minus theoretical. -/// -/// mzSpecLib peak lists carry *observed* m/z, so the theoretical value a -/// library reader wants is `observed - error`. -#[derive(Debug, Clone, Copy, PartialEq)] -pub enum MassError { - /// Absolute, in daltons. - Da(f64), - /// Relative, in parts per million. - Ppm(f64), -} - -impl MassError { - /// Recover the theoretical m/z from the observed one. - pub fn theoretical_from_observed(&self, observed: f64) -> f64 { - match self { - MassError::Da(d) => observed - d, - MassError::Ppm(p) => observed / (1.0 + p * 1e-6), - } - } -} - -/// Split the trailing `/[ppm]` off an annotation, if present. -/// -/// Public because a library reader needs the error even when the ion itself is -/// unrepresentable: the error is what recovers theoretical m/z, so a peak that -/// ends up with an unknown label still gets a correct mass. -pub fn split_mass_error(s: &str) -> Result<(&str, Option), IonParsingError> { - let Some((head, tail)) = s.rsplit_once('/') else { - return Ok((s, None)); - }; - let (num, is_ppm) = match tail.strip_suffix("ppm") { - Some(n) => (n, true), - None => (tail, false), - }; - let v: f64 = num - .parse() - .map_err(|_| IonParsingError::parse(s, "Unable to parse the mass-error suffix"))?; - Ok(( - head, - Some(if is_ppm { - MassError::Ppm(v) - } else { - MassError::Da(v) - }), - )) -} - -impl IonAnnot { - /// Parse the ion itself, with no mass-error suffix. - /// - /// A comma-separated list of alternatives is NOT handled here -- that is - /// ambiguity, and resolving it needs the caller's policy. Split on `,` and - /// parse each alternative. - fn parse_ion(value: &str) -> Result { - // charge: trailing ^N - let (rest, charge) = match value.split_once('^') { - Some((rest, charge)) => { - let charge = charge.parse::().map_err(|_| { - IonParsingError::parse(value, "Unable to parse the charge number") - })?; - (rest, charge) - } - None => (value, 1), - }; - - // isotope: +Ni. Negative isotope offsets are not supported. - let (rest, isotope) = match rest.split_once('+') { - Some((rest, adducts)) => { - let adducts = adducts.strip_suffix('i').ok_or(IonParsingError::parse( - value, - "Only the isotope adduct '+Ni' is supported", - ))?; - let isotope = if adducts.is_empty() { - 1 - } else { - adducts.parse::().map_err(|_| { - IonParsingError::parse(value, "Unable to parse the isotope number") - })? - }; - (rest, isotope) - } - None => (rest, 0), - }; - - // neutral loss: everything from the first '-'. Internal fragments use - // 'm:' and never contain '-' before the loss, so splitting - // at the first '-' is unambiguous. - let (core, loss) = match rest.split_once('-') { - Some((core, loss_expr)) => { - let loss = NeutralLoss::from_expression(loss_expr)?.ok_or_else(|| { - IonParsingError::UnsupportedNeutralLoss { - loss: loss_expr.to_string(), - } - })?; - (core, loss) - } - None => (rest, NeutralLoss::None), - }; - - // internal fragment: m: - if let Some(spans) = core.strip_prefix('m') - && let Some((a, b)) = spans.split_once(':') - { - let start = a.parse::().map_err(|_| { - IonParsingError::parse(value, "Unable to parse internal-fragment start") - })?; - let end = b.parse::().map_err(|_| { - IonParsingError::parse(value, "Unable to parse internal-fragment end") - })?; - return Self::try_new_internal(start, end, charge, isotope, loss); - } - - // immonium: I, bare only. - if let Some(res) = core.strip_prefix('I') { - let mut ch = res.chars(); - return match (ch.next(), ch.next()) { - (Some(r), None) => Self::try_new_immonium(r, charge, isotope, loss), - // A bare `I` names no residue at all; `IC[Carbamidomethyl]` and - // friends carry a mod string no fixed-width field can hold. - // Both are unrepresentable, and the message says which it is. - (None, _) => Err(IonParsingError::parse( - value, - "Immonium ion names no residue", - )), - _ => Err(IonParsingError::UnsupportedModifiedImmonium { - annotation: value.to_string(), - }), - }; - } - - // Backbone / precursor / unknown: a series char then an ordinal. - let mut chars = core.chars(); - let series = chars - .next() - .ok_or(IonParsingError::parse(value, "Empty string"))?; - let rest = chars.as_str(); - let ordinal = - if rest.is_empty() { - None - } else { - Some(rest.parse::().map_err(|_| { - IonParsingError::parse(value, "Ordinal is not a number in 0..=255") - })?) - }; - Self::try_new_with_loss(series, ordinal, charge, isotope, loss) - } -} - /// Guards the one-shot warning in [`IonAnnot::try_from`]. static MASS_ERROR_DISCARDED: std::sync::Once = std::sync::Once::new(); @@ -601,14 +397,15 @@ impl TryFrom<&str> for IonAnnot { ); }); } - Self::parse_ion(ion) + let (series, suffixes) = parse::parse_ion(ion)?; + Self::new(series, suffixes.loss, suffixes.charge, suffixes.isotope) } } impl Display for IonAnnot { - /// Renders the canonical mzPAF spelling. Not byte-inverse to parsing: a - /// non-canonical loss spelling (`-CH3SOH`) renders canonically (`-CH4OS`). - /// The mass-error suffix is not part of the annotation and is not rendered. + /// Renders the canonical mzPAF spelling, the exact inverse of parsing except + /// that a non-canonical loss spelling (`-CH3SOH`) renders canonically + /// (`-CH4OS`). The mass-error suffix is not part of the annotation. fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { write!(f, "{}", self.series_ordinal())?; write!(f, "{}", self.loss())?; @@ -628,56 +425,6 @@ impl Display for IonAnnot { } } -/// Why an annotation could not be parsed or represented. -/// -/// Never matched outside this crate today, so the variants exist for the -/// message a user sees when a library fails to load. Each one names the -/// offending value, because the readers that surface these discard the variant -/// and show only the rendered string. -#[derive(Debug, Error)] -pub enum IonParsingError { - #[error("Ordinal {ordinal} out of range for series '{series}'")] - OrdinalOutOfRange { ordinal: u8, series: char }, - #[error("Series '{series}' requires an ordinal")] - MissingOrdinal { series: char }, - #[error("Series '{series}' takes no ordinal, got {ordinal}")] - UnexpectedOrdinal { series: char, ordinal: u8 }, - #[error("Unsupported fragment type: '{fragment_type}'")] - UnsupportedFragmentType { fragment_type: char }, - #[error("Charge cannot be 0")] - ChargeCannotBeZero, - #[error("Charge {charge} outside the representable range {CHARGE_MIN}..={CHARGE_MAX}")] - ChargeOutOfRange { charge: i8 }, - #[error( - "Isotope offset {isotope} outside the representable range {ISOTOPE_MIN}..={ISOTOPE_MAX}" - )] - IsotopeOutOfRange { isotope: i8 }, - #[error("Neutral loss '{loss}' is not representable")] - UnsupportedNeutralLoss { loss: String }, - #[error("Immonium ions must be a bare uppercase residue, got '{annotation}'")] - UnsupportedModifiedImmonium { annotation: String }, - #[error("Ran out of distinct unknown-ion labels: the 8-bit ordinal is exhausted")] - UnknownIonsExhausted, - #[error("Could not parse '{annotation}': {context}")] - ParsingError { - /// The whole annotation, not the fragment of it that failed: a reader - /// reports this with no row index, so the full text is the only handle - /// the user gets on which row broke. - annotation: String, - context: &'static str, - }, -} - -impl IonParsingError { - /// Shorthand for [`Self::ParsingError`], which is built at thirteen sites. - fn parse(annotation: &str, context: &'static str) -> Self { - Self::ParsingError { - annotation: annotation.to_string(), - context, - } - } -} - /// Hands out `?1`, `?2`, ... for peaks whose annotation this crate cannot /// represent. /// @@ -706,168 +453,6 @@ impl UnknownIonCounter { } } -/// A backbone fragment ion series. -/// -/// The nine mzPAF backbone series differ only by their letter, so the -/// letter-to-discriminant pairing lives in exactly one private table and every -/// match over them is a single arm. -#[derive(Clone, Copy, Eq, PartialEq, Ord, PartialOrd, Hash, Debug)] -#[allow(non_camel_case_types)] -#[repr(u8)] -pub enum Series { - a = 1, - b, - c, - d, - v, - w, - x, - y, - z, -} - -impl Series { - /// Every series, in discriminant order, and parallel to the private letter - /// table -- `series_letters_and_discriminants_agree` pins that pairing. - pub const ALL: [Self; 9] = [ - Self::a, - Self::b, - Self::c, - Self::d, - Self::v, - Self::w, - Self::x, - Self::y, - Self::z, - ]; - /// The mzPAF letters, in discriminant order. The single place the - /// letter↔discriminant pairing lives. - const CHARS: &'static [u8; 9] = b"abcdvwxyz"; - - /// The mzPAF letter for this series. - pub const fn as_char(self) -> char { - Self::CHARS[self as usize - 1] as char - } - - /// The series for an mzPAF letter, or `None` if it names no backbone series. - fn from_char(c: char) -> Option { - let idx = Self::CHARS.iter().position(|&b| b as char == c)?; - Some(Self::ALL[idx]) - } -} - -impl Display for Series { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "{}", self.as_char()) - } -} - -/// The logical series-and-payload view of an [`IonAnnot`]. -/// -/// This is a *view*: `IonAnnot` stores a packed word and reconstructs this on -/// demand. Constructing one directly does not allocate an annotation. -#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy)] -#[allow(non_camel_case_types)] -pub enum IonSeriesOrdinal { - /// One of the nine backbone series, at position `ordinal` in the ladder. - backbone { - series: Series, - ordinal: u8, - }, - /// An unannotated peak. `ordinal` is a uniqueness counter, not a position. - unknown { - ordinal: u8, - }, - precursor, - /// An internal fragment spanning residues `start..=end`. - internal { - start: u8, - end: u8, - }, - /// A bare immonium ion for an uppercase residue code. - immonium { - residue: char, - }, -} - -impl IonSeriesOrdinal { - /// Split into the `kind` discriminant and its `payload`, the two fields - /// [`IonAnnot`] packs. - /// - /// This and [`Self::from_parts`] are the only place the numbering lives. - /// The nine backbone series share one arm, so their discriminants come - /// from [`Series`] itself and cannot drift out of step with their letters. - const fn to_parts(self) -> (u32, u32) { - match self { - Self::backbone { series, ordinal } => (series as u32, ordinal as u32), - // Discriminant 0, so the all-zero word -- `IonAnnot::default()` -- - // is the unknown ion `?0` rather than an undecodable value. - Self::unknown { ordinal } => (0, ordinal as u32), - Self::precursor => (10, 0), - Self::internal { start, end } => { - (11, (start as u32) | ((end as u32) << INTERNAL_POS_BITS)) - } - Self::immonium { residue } => (12, (residue as u8 - b'A') as u32), - } - } - - /// Inverse of [`Self::to_parts`]. Total by construction: `unknown` is the - /// catch-all discriminant, so a value this build does not recognise -- only - /// reachable from a corrupted word -- decodes as an unknown ion instead of - /// panicking on a path `Display` (and therefore `Serialize`) reaches. - const fn from_parts(kind: u32, payload: u32) -> Self { - let ordinal = payload as u8; - match kind { - 1..=9 => Self::backbone { - // `kind` is in range, so this is the inverse of `series as u32`. - series: Series::ALL[kind as usize - 1], - ordinal, - }, - 10 => Self::precursor, - 11 => Self::internal { - start: (payload & mask(INTERNAL_POS_BITS)) as u8, - end: ((payload >> INTERNAL_POS_BITS) & mask(INTERNAL_POS_BITS)) as u8, - }, - 12 => Self::immonium { - residue: (b'A' + (payload & mask(IMMONIUM_BITS)) as u8) as char, - }, - _ => Self::unknown { ordinal }, - } - } - - /// Build the series view from an mzPAF series letter and its ordinal. - /// - /// `p` is the only letter that takes no ordinal; every other one requires - /// one. Internal fragments and immonium ions are spelled differently and - /// have their own constructors. - fn from_series_char(c: char, ordinal: Option) -> Result { - if c == 'p' { - return match ordinal { - None => Ok(Self::precursor), - Some(ordinal) => Err(IonParsingError::UnexpectedOrdinal { series: c, ordinal }), - }; - } - let ordinal = ordinal.ok_or(IonParsingError::MissingOrdinal { series: c })?; - if c == '?' { - return Ok(Self::unknown { ordinal }); - } - Series::from_char(c) - .map(|series| Self::backbone { series, ordinal }) - .ok_or(IonParsingError::UnsupportedFragmentType { fragment_type: c }) - } -} - -impl Display for IonSeriesOrdinal { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - match self { - Self::backbone { series, ordinal } => write!(f, "{}{}", series.as_char(), ordinal), - Self::unknown { ordinal } => write!(f, "?{}", ordinal), - Self::precursor => write!(f, "p"), - Self::internal { start, end } => write!(f, "m{}:{}", start, end), - Self::immonium { residue } => write!(f, "I{}", residue), - } - } -} #[cfg(test)] mod tests { use super::*; @@ -893,8 +478,13 @@ mod tests { assert_eq!(size_of::<(IonAnnot, f32)>(), 8); // All four of these are new capacity, and none of them widened the word. - let loaded = IonAnnot::try_new_internal(2, 11, 3, 2, NeutralLoss::PhosphoricAcidWater) - .expect("every field at once"); + let loaded = IonAnnot::new( + IonSeriesOrdinal::try_internal(2, 11).expect("in range"), + NeutralLoss::PhosphoricAcidWater, + 3, + 2, + ) + .expect("every field at once"); assert_eq!(size_of_val(&loaded), 4); assert_eq!(loaded.loss(), NeutralLoss::PhosphoricAcidWater); assert_eq!( @@ -961,9 +551,11 @@ mod tests { NeutralLoss::Water, NeutralLoss::PhosphoricAcidWater, ] { - let a = - IonAnnot::try_new_with_loss('y', Some(ordinal), charge, isotope, loss) - .expect("in range"); + let series = IonSeriesOrdinal::backbone { + series: Series::y, + ordinal, + }; + let a = IonAnnot::new(series, loss, charge, isotope).expect("in range"); assert_eq!(a.get_charge(), charge); assert_eq!(a.get_isotope(), isotope); assert_eq!(a.try_get_ordinal(), Some(ordinal)); @@ -985,8 +577,11 @@ mod tests { Err(IonParsingError::IsotopeOutOfRange { .. }) )); assert!(IonAnnot::try_new('y', Some(1), 0, 0).is_err()); + // The payload bounds are `series`'s to enforce; see + // `constructors_reject_what_the_payload_cannot_hold` there. This pins + // that an out-of-range span cannot reach a word through the parser. assert!(matches!( - IonAnnot::try_new_internal(INTERNAL_POS_MAX + 1, 1, 1, 0, NeutralLoss::None), + IonAnnot::try_from("m64:1"), Err(IonParsingError::OrdinalOutOfRange { .. }) )); } @@ -1097,106 +692,6 @@ mod tests { assert_eq!(IonAnnot::try_from("y1/-0.0005").unwrap(), ion("y1")); } - /// One representative of every `IonSeriesOrdinal` case: all nine backbone - /// series (each with a distinct ordinal, so a transposition in - /// `to_parts`/`from_parts` cannot cancel out) plus the four others, at - /// their field boundaries. - fn all_series() -> Vec { - let mut out: Vec = Series::ALL - .iter() - .enumerate() - .map(|(i, &series)| IonSeriesOrdinal::backbone { - series, - ordinal: i as u8 + 1, - }) - .collect(); - out.extend([ - IonSeriesOrdinal::precursor, - IonSeriesOrdinal::unknown { ordinal: 10 }, - IonSeriesOrdinal::internal { start: 2, end: 11 }, - // Both endpoints at their 6-bit ceiling: the widest payload the - // 12 bits hold, and the case a narrowed field would silently clip. - IonSeriesOrdinal::internal { - start: INTERNAL_POS_MAX, - end: INTERNAL_POS_MAX, - }, - IonSeriesOrdinal::immonium { residue: 'A' }, - // Highest residue index the 5-bit immonium field must hold. - IonSeriesOrdinal::immonium { residue: 'Z' }, - ]); - out - } - - /// `to_parts` and `from_parts` are hand-written inverses. Without this, - /// swapping two arms (`v` encoding as `w`) mislabels a whole ion series - /// and every other test still passes. - #[test] - fn every_series_variant_round_trips_through_the_packed_word() { - for series in all_series() { - let (kind, payload) = series.to_parts(); - assert_eq!(IonSeriesOrdinal::from_parts(kind, payload), series); - assert!(kind <= mask(KIND_BITS), "{series:?} kind overflows"); - assert!( - payload <= mask(PAYLOAD_BITS), - "{series:?} payload overflows" - ); - } - - // Every backbone series must land on its own discriminant; the other - // four are singletons and are covered by the round trip above. - let mut kinds: Vec = Series::ALL - .iter() - .map(|&series| { - IonSeriesOrdinal::backbone { series, ordinal: 1 } - .to_parts() - .0 - }) - .collect(); - kinds.sort_unstable(); - kinds.dedup(); - assert_eq!(kinds.len(), Series::ALL.len(), "two series share a kind"); - } - - /// `from_parts` claims to be total, and three of the sixteen `kind` - /// values are unassigned. Renumbering the assigned ones is safe only - /// while the unassigned ones stay decodable, since a corrupted word - /// reaches this on a `Display` (and so `Serialize`) path. - #[test] - fn every_kind_bit_pattern_decodes_without_panicking() { - for kind in 0..=mask(KIND_BITS) { - for payload in [0, 1, mask(PAYLOAD_BITS)] { - // Also exercises `Display`, which is where a partial decode - // would have panicked. - let _ = IonSeriesOrdinal::from_parts(kind, payload).to_string(); - } - } - } - - /// `Display`, `from_series_char` and the parser must agree with the - /// packing for every case. - #[test] - fn every_series_variant_round_trips_through_its_mzpaf_spelling() { - for series in all_series() { - let annot = IonAnnot::pack(series, NeutralLoss::None, 1, 0).expect("valid"); - assert_eq!(annot.series_ordinal(), series); - let text = annot.to_string(); - assert_eq!(ion(&text).series_ordinal(), series, "{text}"); - } - } - - /// The letters are the mzPAF spelling of the discriminants, and - /// `from_char`/`as_char` are hand-written inverses of each other. - #[test] - fn series_letters_and_discriminants_agree() { - assert_eq!(Series::ALL.len(), Series::CHARS.len()); - for &series in &Series::ALL { - assert_eq!(Series::from_char(series.as_char()), Some(series)); - } - for c in ['p', '?', 'm', 'I', 'q', 'A'] { - assert_eq!(Series::from_char(c), None, "{c} is not a backbone series"); - } - } - /// An unrepresentable loss must fail loudly. Parsing `y1-HCOOH` as plain /// `y1` would put a loss peak's m/z on the `y1` label and collide with the /// real `y1`. diff --git a/rust/micromzpaf/src/parse.rs b/rust/micromzpaf/src/parse.rs new file mode 100644 index 00000000..875078f2 --- /dev/null +++ b/rust/micromzpaf/src/parse.rs @@ -0,0 +1,144 @@ +//! Peeling an mzPAF annotation apart, and the mass-error suffix. +//! +//! The suffixes come off right to left, and that order is not arbitrary: each +//! delimiter can only appear to the right of the previous one, so peeling in +//! this sequence is unambiguous while any other order is not. +//! +//! ```text +//! y5 - H2O + 2i ^ 3 / -0.0003 +//! └core┘└loss┘└iso┘└ch┘└error┘ +//! 4 3 2 1 0 <- order peeled +//! ``` +//! +//! `0` is [`split_mass_error`], `1..=3` are `peel_suffixes`, and `4` is +//! `IonSeriesOrdinal::parse`, which owns the spelling of the ion itself. + +use crate::series::IonSeriesOrdinal; +use crate::{ + IonParsingError, + NeutralLoss, +}; + +/// A mass-error suffix on an mzPAF annotation: observed minus theoretical. +/// +/// mzSpecLib peak lists carry *observed* m/z, so the theoretical value a +/// library reader wants is `observed - error`. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum MassError { + /// Absolute, in daltons. + Da(f64), + /// Relative, in parts per million. + Ppm(f64), +} + +impl MassError { + /// Recover the theoretical m/z from the observed one. + pub fn theoretical_from_observed(&self, observed: f64) -> f64 { + match self { + MassError::Da(d) => observed - d, + MassError::Ppm(p) => observed / (1.0 + p * 1e-6), + } + } +} + +/// Split the trailing `/[ppm]` off an annotation, if present. +/// +/// Public because a library reader needs the error even when the ion itself is +/// unrepresentable: the error is what recovers theoretical m/z, so a peak that +/// ends up with an unknown label still gets a correct mass. +pub fn split_mass_error(s: &str) -> Result<(&str, Option), IonParsingError> { + let Some((head, tail)) = s.rsplit_once('/') else { + return Ok((s, None)); + }; + let (num, is_ppm) = match tail.strip_suffix("ppm") { + Some(n) => (n, true), + None => (tail, false), + }; + let v: f64 = num + .parse() + .map_err(|_| IonParsingError::parse(s, "mass-error suffix is not a number"))?; + Ok(( + head, + Some(if is_ppm { + MassError::Ppm(v) + } else { + MassError::Da(v) + }), + )) +} + +/// The three suffixes an annotation can carry after its ion kind. +pub(crate) struct Suffixes { + pub(crate) charge: i8, + pub(crate) isotope: i8, + pub(crate) loss: NeutralLoss, +} + +/// Peel `^N`, `+Ni` and `-` off `value`, returning the ion core and the +/// suffixes. `value` must already have had any mass-error suffix removed. +/// +/// Defaults when a suffix is absent: charge 1, isotope 0, no loss. +pub(crate) fn peel_suffixes(value: &str) -> Result<(&str, Suffixes), IonParsingError> { + // Charge: `^N`, rightmost of the three. + let (rest, charge) = match value.split_once('^') { + Some((rest, charge)) => ( + rest, + charge + .parse::() + .map_err(|_| IonParsingError::parse(value, "charge is not a number"))?, + ), + None => (value, 1), + }; + + // Isotope: `+Ni`, or bare `+i` for one. + let (rest, isotope) = match rest.split_once('+') { + Some((rest, adducts)) => { + let digits = adducts.strip_suffix('i').ok_or(IonParsingError::parse( + value, + "only the isotope adduct '+Ni' is supported", + ))?; + let isotope = if digits.is_empty() { + 1 + } else { + digits + .parse::() + .map_err(|_| IonParsingError::parse(value, "isotope count is not a number"))? + }; + (rest, isotope) + } + None => (rest, 0), + }; + + // Loss: everything from the first `-`. Internal fragments spell themselves + // `m:` and so never contain a `-` left of the loss. + let (core, loss) = match rest.split_once('-') { + Some((core, expr)) => { + let loss = NeutralLoss::from_expression(expr)?.ok_or_else(|| { + IonParsingError::UnsupportedNeutralLoss { + loss: expr.to_string(), + } + })?; + (core, loss) + } + None => (rest, NeutralLoss::None), + }; + + Ok(( + core, + Suffixes { + charge, + isotope, + loss, + }, + )) +} + +/// Parse a whole annotation, minus any mass-error suffix, into its parts. +/// +/// A comma-separated list of alternatives is NOT handled here -- that is +/// ambiguity, and resolving it needs the caller's policy. Split on `,` and +/// parse each alternative. +pub(crate) fn parse_ion(value: &str) -> Result<(IonSeriesOrdinal, Suffixes), IonParsingError> { + let (core, suffixes) = peel_suffixes(value)?; + Ok((IonSeriesOrdinal::parse(core)?, suffixes)) +} diff --git a/rust/micromzpaf/src/series.rs b/rust/micromzpaf/src/series.rs new file mode 100644 index 00000000..f3834ba8 --- /dev/null +++ b/rust/micromzpaf/src/series.rs @@ -0,0 +1,458 @@ +//! What kind of ion an annotation names, and how that fits in a packed payload. +//! +//! This module owns three things that have to agree, and only agree if they sit +//! together: +//! +//! 1. the mzPAF **spelling** of each kind (`Display` and `IonSeriesOrdinal::parse`, +//! which are exact inverses), +//! 2. the **payload layout** each kind packs into (`IonSeriesOrdinal::to_parts` +//! and `from_parts`, also inverses), +//! 3. the **bounds** that layout imposes, enforced by the constructors. +//! +//! Splitting any of those across a module boundary is what lets a `v` render as +//! a `w`, or an internal-fragment endpoint get validated against a width it no +//! longer has. [`IonAnnot`](crate::IonAnnot) owns the surrounding word -- +//! charge, isotope, loss, and where these two fields sit in it -- and nothing +//! else. + +use std::fmt::Display; + +use crate::{ + IonParsingError, + mask, +}; + +/// Width of the `kind` discriminant in the packed word. +pub(crate) const KIND_BITS: u32 = 4; +/// Width of the `payload` this module encodes into. +pub(crate) const PAYLOAD_BITS: u32 = 12; +/// Width of each internal-fragment endpoint inside `payload`. Two share it. +const INTERNAL_POS_BITS: u32 = 6; +/// Width of the immonium residue index inside `payload`. +const IMMONIUM_BITS: u32 = 5; + +/// Widest residue position an internal fragment endpoint holds. +/// +/// Narrower than a backbone ordinal because two endpoints share one payload. +/// Internal fragments are bounded by peptide length, so 63 is well past tryptic. +pub const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; + +const _: () = assert!( + 2 * INTERNAL_POS_BITS <= PAYLOAD_BITS, + "two internal-fragment endpoints must fit in one payload" +); +const _: () = assert!( + IMMONIUM_BITS <= PAYLOAD_BITS && mask(IMMONIUM_BITS) >= (b'Z' - b'A') as u32, + "the immonium field must fit the payload and hold every uppercase residue" +); + +/// `kind` discriminants. The numbering lives only here; [`IonSeriesOrdinal::to_parts`] +/// and `from_parts` are its only readers. +/// +/// Backbone ions take 1..=9 from [`Series`] itself, so a series can never drift +/// out of step with its letter. `unknown` is 0 so that the all-zero word decodes +/// to a real annotation. 13..=15 are unassigned and decode as `unknown`. +const KIND_UNKNOWN: u32 = 0; +const KIND_PRECURSOR: u32 = 10; +const KIND_INTERNAL: u32 = 11; +const KIND_IMMONIUM: u32 = 12; + +const _: () = assert!( + KIND_IMMONIUM <= mask(KIND_BITS), + "a kind discriminant overflows KIND_BITS" +); + +/// A backbone fragment ion series. +/// +/// The nine mzPAF backbone series differ only by their letter, so the +/// letter-to-discriminant pairing lives in exactly one private table and every +/// match over them is a single arm. +#[derive(Clone, Copy, Eq, PartialEq, Ord, PartialOrd, Hash, Debug)] +#[allow(non_camel_case_types)] +#[repr(u8)] +pub enum Series { + a = 1, + b, + c, + d, + v, + w, + x, + y, + z, +} + +impl Series { + /// Every series, in discriminant order, and parallel to the private letter + /// table -- `series_letters_and_discriminants_agree` pins that pairing. + pub const ALL: [Self; 9] = [ + Self::a, + Self::b, + Self::c, + Self::d, + Self::v, + Self::w, + Self::x, + Self::y, + Self::z, + ]; + /// The mzPAF letters, in discriminant order. The single place the + /// letter-to-discriminant pairing lives. + const CHARS: &'static [u8; 9] = b"abcdvwxyz"; + + /// The mzPAF letter for this series. + pub const fn as_char(self) -> char { + Self::CHARS[self as usize - 1] as char + } + + /// The series for an mzPAF letter, or `None` if it names no backbone series. + fn from_char(c: char) -> Option { + let idx = Self::CHARS.iter().position(|&b| b as char == c)?; + Some(Self::ALL[idx]) + } +} + +impl Display for Series { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.as_char()) + } +} + +/// What kind of ion an annotation names, with whatever that kind carries. +/// +/// A *view*: [`IonAnnot`](crate::IonAnnot) stores a packed word and reconstructs +/// this on demand, so building one allocates no annotation. Fields are public +/// because reading them is the point; the validating constructors +/// ([`Self::try_internal`], [`Self::try_immonium`]) are how you get one that is known to +/// fit the payload. +#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy)] +#[allow(non_camel_case_types)] +pub enum IonSeriesOrdinal { + /// One of the nine backbone series, at position `ordinal` in the ladder. + backbone { + series: Series, + ordinal: u8, + }, + /// An unannotated peak. `ordinal` is a uniqueness counter, not a position. + unknown { + ordinal: u8, + }, + precursor, + /// An internal fragment spanning residues `start..=end`. + internal { + start: u8, + end: u8, + }, + /// A bare immonium ion for an uppercase residue code. + immonium { + residue: char, + }, +} + +impl IonSeriesOrdinal { + /// An internal fragment spanning residues `start..=end`. + /// + /// Endpoints are bounded by [`INTERNAL_POS_MAX`], which is this module's + /// payload width and is checked here rather than by the caller. + pub fn try_internal(start: u8, end: u8) -> Result { + if start > INTERNAL_POS_MAX || end > INTERNAL_POS_MAX { + return Err(IonParsingError::OrdinalOutOfRange { + ordinal: start.max(end), + series: 'm', + }); + } + Ok(Self::internal { start, end }) + } + + /// A bare immonium ion for an uppercase residue code. + /// + /// Modified immonium (`IC[Carbamidomethyl]`) carries an arbitrary mod string + /// that no fixed-width field holds, so only a bare uppercase residue is + /// accepted. + pub fn try_immonium(residue: char) -> Result { + if !residue.is_ascii_uppercase() { + return Err(IonParsingError::UnsupportedFragmentType { + fragment_type: residue, + }); + } + Ok(Self::immonium { residue }) + } + + /// A backbone / precursor / unknown ion from its mzPAF letter. + /// + /// `p` is the only letter that takes no ordinal; every other one requires + /// one. Internal fragments and immonium ions are spelled differently and + /// have their own constructors. + pub fn from_series_char(c: char, ordinal: Option) -> Result { + if c == 'p' { + return match ordinal { + None => Ok(Self::precursor), + Some(ordinal) => Err(IonParsingError::UnexpectedOrdinal { series: c, ordinal }), + }; + } + let ordinal = ordinal.ok_or(IonParsingError::MissingOrdinal { series: c })?; + if c == '?' { + return Ok(Self::unknown { ordinal }); + } + Series::from_char(c) + .map(|series| Self::backbone { series, ordinal }) + .ok_or(IonParsingError::UnsupportedFragmentType { fragment_type: c }) + } + + /// Parse the ion-kind part of an annotation: everything left of the charge, + /// isotope and loss suffixes. + /// + /// The exact inverse of this type's [`Display`], which is why the two live + /// side by side -- `every_series_variant_round_trips_through_its_mzpaf_spelling` + /// pins that. + pub(crate) fn parse(core: &str) -> Result { + // Internal fragment: `m:`. Checked before the single-letter + // forms because `m` is not a backbone letter, so there is no ambiguity. + if let Some(spans) = core.strip_prefix('m') + && let Some((a, b)) = spans.split_once(':') + { + let start = a + .parse::() + .map_err(|_| IonParsingError::parse(core, "internal-fragment start is not a u8"))?; + let end = b + .parse::() + .map_err(|_| IonParsingError::parse(core, "internal-fragment end is not a u8"))?; + return Self::try_internal(start, end); + } + + // Immonium: `I`, bare only. + if let Some(res) = core.strip_prefix('I') { + let mut ch = res.chars(); + return match (ch.next(), ch.next()) { + (Some(r), None) => Self::try_immonium(r), + (None, _) => Err(IonParsingError::parse( + core, + "immonium ion names no residue", + )), + // Carries a mod string no fixed-width field can hold. + _ => Err(IonParsingError::UnsupportedModifiedImmonium { + annotation: core.to_string(), + }), + }; + } + + // Backbone / precursor / unknown: a series letter then an ordinal. + let mut chars = core.chars(); + let series = chars + .next() + .ok_or(IonParsingError::parse(core, "empty annotation"))?; + let rest = chars.as_str(); + let ordinal = + if rest.is_empty() { + None + } else { + Some(rest.parse::().map_err(|_| { + IonParsingError::parse(core, "ordinal is not a number in 0..=255") + })?) + }; + Self::from_series_char(series, ordinal) + } + + /// Split into the `kind` discriminant and its `payload`, the two fields + /// [`IonAnnot`](crate::IonAnnot) packs. + /// + /// Total: every arm masks to its own width, so a value built by hand rather + /// than through the constructors encodes as something decodable instead of + /// panicking on a path `Display` (and so `Serialize`) reaches. + pub(crate) const fn to_parts(self) -> (u32, u32) { + match self { + Self::backbone { series, ordinal } => (series as u32, ordinal as u32), + Self::unknown { ordinal } => (KIND_UNKNOWN, ordinal as u32), + Self::precursor => (KIND_PRECURSOR, 0), + Self::internal { start, end } => ( + KIND_INTERNAL, + ((start as u32) & mask(INTERNAL_POS_BITS)) + | (((end as u32) & mask(INTERNAL_POS_BITS)) << INTERNAL_POS_BITS), + ), + // `saturating_sub` rather than `-`: the constructor rejects a + // non-uppercase residue, but the variant is public and this must not + // underflow for one built directly. + Self::immonium { residue } => ( + KIND_IMMONIUM, + (residue as u32).saturating_sub(b'A' as u32) & mask(IMMONIUM_BITS), + ), + } + } + + /// Inverse of [`Self::to_parts`]. Total by construction: `unknown` is the + /// catch-all discriminant, so a value this build does not recognise -- only + /// reachable from a corrupted word -- decodes as an unknown ion. + pub(crate) const fn from_parts(kind: u32, payload: u32) -> Self { + let ordinal = payload as u8; + match kind { + 1..=9 => Self::backbone { + // `kind` is in range, so this is the inverse of `series as u32`. + series: Series::ALL[kind as usize - 1], + ordinal, + }, + KIND_PRECURSOR => Self::precursor, + KIND_INTERNAL => Self::internal { + start: (payload & mask(INTERNAL_POS_BITS)) as u8, + end: ((payload >> INTERNAL_POS_BITS) & mask(INTERNAL_POS_BITS)) as u8, + }, + KIND_IMMONIUM => Self::immonium { + residue: (b'A' + (payload & mask(IMMONIUM_BITS)) as u8) as char, + }, + _ => Self::unknown { ordinal }, + } + } +} + +impl Display for IonSeriesOrdinal { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::backbone { series, ordinal } => write!(f, "{}{}", series.as_char(), ordinal), + Self::unknown { ordinal } => write!(f, "?{}", ordinal), + Self::precursor => write!(f, "p"), + Self::internal { start, end } => write!(f, "m{}:{}", start, end), + Self::immonium { residue } => write!(f, "I{}", residue), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// One representative of every [`IonSeriesOrdinal`] case: all nine backbone + /// series (each with a distinct ordinal, so a transposition in + /// `to_parts`/`from_parts` cannot cancel out) plus the four others, at their + /// field boundaries. + fn all_series() -> Vec { + let mut out: Vec = Series::ALL + .iter() + .enumerate() + .map(|(i, &series)| IonSeriesOrdinal::backbone { + series, + ordinal: i as u8 + 1, + }) + .collect(); + out.extend([ + IonSeriesOrdinal::precursor, + IonSeriesOrdinal::unknown { ordinal: 10 }, + IonSeriesOrdinal::internal { start: 2, end: 11 }, + // Both endpoints at their 6-bit ceiling: the widest payload the + // 12 bits hold, and the case a narrowed field would silently clip. + IonSeriesOrdinal::internal { + start: INTERNAL_POS_MAX, + end: INTERNAL_POS_MAX, + }, + IonSeriesOrdinal::immonium { residue: 'A' }, + // Highest residue index the 5-bit immonium field must hold. + IonSeriesOrdinal::immonium { residue: 'Z' }, + ]); + out + } + + /// `to_parts` and `from_parts` are hand-written inverses. Without this, + /// swapping two arms (`v` encoding as `w`) mislabels a whole ion series and + /// every other test still passes. + #[test] + fn every_variant_round_trips_through_the_packed_parts() { + for series in all_series() { + let (kind, payload) = series.to_parts(); + assert_eq!(IonSeriesOrdinal::from_parts(kind, payload), series); + assert!(kind <= mask(KIND_BITS), "{series:?} kind overflows"); + assert!( + payload <= mask(PAYLOAD_BITS), + "{series:?} payload overflows" + ); + } + + // Every backbone series must land on its own discriminant; the other + // four are singletons and are covered by the round trip above. + let mut kinds: Vec = Series::ALL + .iter() + .map(|&series| { + IonSeriesOrdinal::backbone { series, ordinal: 1 } + .to_parts() + .0 + }) + .collect(); + kinds.sort_unstable(); + kinds.dedup(); + assert_eq!(kinds.len(), Series::ALL.len(), "two series share a kind"); + } + + /// `from_parts` claims to be total, and three of the sixteen `kind` values + /// are unassigned. Renumbering the assigned ones is safe only while the + /// unassigned ones stay decodable, since a corrupted word reaches this on a + /// `Display` (and so `Serialize`) path. + #[test] + fn every_kind_bit_pattern_decodes_without_panicking() { + for kind in 0..=mask(KIND_BITS) { + for payload in [0, 1, mask(PAYLOAD_BITS)] { + // Also exercises `Display`, which is where a partial decode + // would have panicked. + let _ = IonSeriesOrdinal::from_parts(kind, payload).to_string(); + } + } + } + + /// `Display` and [`IonSeriesOrdinal::parse`] are hand-written inverses, and + /// this module exists so they cannot drift apart. + #[test] + fn every_variant_round_trips_through_its_mzpaf_spelling() { + for series in all_series() { + let text = series.to_string(); + assert_eq!( + IonSeriesOrdinal::parse(&text).expect("own spelling parses"), + series, + "{text}" + ); + } + } + + /// The letters are the mzPAF spelling of the discriminants, and + /// `from_char`/`as_char` are hand-written inverses of each other. + #[test] + fn series_letters_and_discriminants_agree() { + assert_eq!(Series::ALL.len(), Series::CHARS.len()); + for &series in &Series::ALL { + assert_eq!(Series::from_char(series.as_char()), Some(series)); + } + for c in ['p', '?', 'm', 'I', 'q', 'A'] { + assert_eq!(Series::from_char(c), None, "{c} is not a backbone series"); + } + } + + /// The bounds this module owns, checked where they are enforced rather than + /// through a packed word. + #[test] + fn constructors_reject_what_the_payload_cannot_hold() { + assert!(IonSeriesOrdinal::try_internal(INTERNAL_POS_MAX, INTERNAL_POS_MAX).is_ok()); + assert!(matches!( + IonSeriesOrdinal::try_internal(INTERNAL_POS_MAX + 1, 1), + Err(IonParsingError::OrdinalOutOfRange { .. }) + )); + assert!(IonSeriesOrdinal::try_immonium('A').is_ok()); + assert!(IonSeriesOrdinal::try_immonium('Z').is_ok()); + // Lowercase would underflow the `- b'A'` in `to_parts`. + assert!(IonSeriesOrdinal::try_immonium('a').is_err()); + + // A bare `I` names no residue; a modified one carries a mod string. + assert!(IonSeriesOrdinal::parse("I").is_err()); + assert!(matches!( + IonSeriesOrdinal::parse("IC[Carbamidomethyl]"), + Err(IonParsingError::UnsupportedModifiedImmonium { .. }) + )); + // `p` is the only letter that takes no ordinal, and requires none. + assert!(matches!( + IonSeriesOrdinal::parse("p1"), + Err(IonParsingError::UnexpectedOrdinal { .. }) + )); + assert!(matches!( + IonSeriesOrdinal::parse("y"), + Err(IonParsingError::MissingOrdinal { .. }) + )); + assert!(matches!( + IonSeriesOrdinal::parse("q1"), + Err(IonParsingError::UnsupportedFragmentType { .. }) + )); + } +} From 501bb57d322ef2ba0fb36f9134106e1aa869f0c2 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 12:15:03 -0700 Subject: [PATCH 6/8] feat(timsquery): keep DIA-NN's neutral-loss fragments instead of dropping them The `.speclib` reader discarded every fragment carrying a neutral loss -- 152 on the pinned fixture -- because `IonAnnot` could not represent one. It can now, so the only thing missing was DIA-NN's loss-code numbering, which is undocumented. Measured it instead of guessing. A lossy fragment and its no-loss sibling differ by exactly the loss, so `(no_loss_mz - lossy_mz) * charge` is the neutral mass. On the fixture, code 1 gives 18.011 across 47 pairs and code 2 gives 17.027 across 37: water and ammonia, with nothing else within half a dalton. The two account for all 152. `loss_codes_are_water_and_ammonia` re-derives that from the fixture rather than trusting the table, and fails on any code this build does not map -- so a DIA-NN version that renumbers them stops the load instead of relabelling peaks. An unmapped code is still dropped and counted, because a wrong loss puts a real m/z on a label that collides with a different real fragment. --- rust/timsquery/src/lib.rs | 1 + rust/timsquery/src/serde/diann_speclib_io.rs | 159 +++++++++++++++++-- 2 files changed, 146 insertions(+), 14 deletions(-) diff --git a/rust/timsquery/src/lib.rs b/rust/timsquery/src/lib.rs index 84fd7bf8..65aa71aa 100644 --- a/rust/timsquery/src/lib.rs +++ b/rust/timsquery/src/lib.rs @@ -58,6 +58,7 @@ pub mod ion { IonAnnot, IonParsingError, IonSeriesOrdinal, + NeutralLoss, Series, UnknownIonCounter, }; diff --git a/rust/timsquery/src/serde/diann_speclib_io.rs b/rust/timsquery/src/serde/diann_speclib_io.rs index 34a58f8a..fa7c1e20 100644 --- a/rust/timsquery/src/serde/diann_speclib_io.rs +++ b/rust/timsquery/src/serde/diann_speclib_io.rs @@ -30,7 +30,11 @@ use super::library_file::{ TargetReadingError, TargetTable, }; -use crate::ion::IonAnnot; +use crate::ion::{ + IonAnnot, + IonSeriesOrdinal, + NeutralLoss, +}; use crate::models::{ Row, TargetCapabilities, @@ -572,9 +576,9 @@ pub struct SpeclibDecodeStats { /// Fragments flagged `ExcludeFromAssay` (`type & 0x80`). Kept (see /// `map_entry`); counted for reporting only. pub exclude_flagged: usize, - /// Fragments with a neutral loss (`loss != 0`). `IonAnnot` can represent - /// these; what is missing is the DIA-NN loss-code mapping (see #105), so - /// they are dropped and counted to keep the cost visible. + /// Fragments whose neutral-loss code `loss_from_code` does not map. Water + /// and ammonia are mapped and kept; anything else is dropped rather than + /// guessed, and counted so the cost stays visible. pub loss_dropped: usize, /// Fragments whose ion-type code or recovered series was unusable. pub unknown_ion_dropped: usize, @@ -703,6 +707,28 @@ fn append_arena(dst: &mut TargetColumns, mut src: TargetColumns Option { + match code { + 0 => Some(NeutralLoss::None), + 1 => Some(NeutralLoss::Water), + 2 => Some(NeutralLoss::Ammonia), + _ => None, + } +} + /// Map a decoded target `EntryView` directly into the columnar arena, folding /// drop counters into `stats` and appending each kept fragment's reference /// intensity to `frag_intens` IN THE SAME ORDER as the pushed fragment labels @@ -767,13 +793,18 @@ fn map_entry( if f.typ() & 0x80 != 0 { stats.exclude_flagged += 1; } - // `IonAnnot` can represent neutral losses, but the map from DIA-NN's - // loss code byte to `NeutralLoss` is not written, so these are dropped. - // The `loss_dropped` assertion in the tests below measures the cost. #105 - if f.loss() != 0 { - stats.loss_dropped += 1; - continue; - } + let loss = match loss_from_code(f.loss()) { + Some(loss) => loss, + None => { + warn!( + "speclib entry {:?}: unknown neutral-loss code {}; dropping fragment", + name, + f.loss() + ); + stats.loss_dropped += 1; + continue; + } + }; let type_char = match f.typ() & 0x7F { 1 => 'b', @@ -830,7 +861,9 @@ fn map_entry( } }; - let ion = match IonAnnot::try_new(type_char, Some(series as u8), charge, 0) { + let ion = match IonSeriesOrdinal::from_series_char(type_char, Some(series as u8)) + .and_then(|series| IonAnnot::new(series, loss, charge, 0)) + { Ok(ion) => ion, Err(e) => { warn!( @@ -1158,10 +1191,20 @@ mod tests { let (geom, _intens, stats, _eof) = parse_speclib_reader(std::io::BufReader::new(file)).unwrap(); // Reference parser: 8384 ExcludeFromAssay-flagged (kept, counted only), - // 152 neutral-loss dropped, no dc != 0 entries. + // no dc != 0 entries. assert_eq!(stats.exclude_flagged, 8384); - assert_eq!(stats.loss_dropped, 152); assert_eq!(stats.decoys_dropped, 0); + // Every neutral loss in this fixture is water or ammonia, both of which + // `loss_from_code` maps, so nothing is dropped for its loss. This was + // 152 before that mapping existed. + assert_eq!(stats.loss_dropped, 0); + // And they arrive as labelled losses rather than as bare ions. + let with_loss = geom + .frag_labels + .iter() + .filter(|l| l.loss() != NeutralLoss::None) + .count(); + assert_eq!(with_loss, 152, "every dropped loss is now a labelled one"); // Exclude-flagged fragments are present in the output, not dropped: // entry[0] keeps its flagged y9. assert!( @@ -1171,4 +1214,92 @@ mod tests { "flagged y9 must be kept" ); } + + /// Re-derives [`loss_from_code`] from the fixture instead of trusting it. + /// + /// DIA-NN's loss numbering is undocumented, so the mapping was measured: a + /// lossy fragment and its no-loss sibling differ by exactly the loss, so + /// `(base_mz - lossy_mz) * charge` is the neutral mass. If a future fixture + /// or DIA-NN version renumbers the codes, this fails here rather than + /// silently relabelling peaks. + /// + /// The tolerance is 3 mDa: the m/z values on disk are `f32`, and DIA-NN's + /// own loss constants are a shade off the monoisotopic masses (water reads + /// 18.011 against 18.0106). That is far tighter than the gap to any other + /// candidate loss, which is what makes the identification safe. + #[test] + fn loss_codes_are_water_and_ammonia() { + use std::collections::BTreeMap; + let mut buf = Vec::new(); + std::fs::File::open(fixture_path()) + .unwrap() + .read_to_end(&mut buf) + .unwrap(); + let lib = SpecLib::open(std::io::Cursor::new(&buf[..])).unwrap(); + + let mut hist: BTreeMap = BTreeMap::new(); + let mut deltas: BTreeMap> = BTreeMap::new(); + + let mut c = Cursor::new(lib.data.bytes()); + c.pos = lib.entries_start; + for _ in 0..lib.n_entries { + let entry = read_entry(&mut c, lib.version, &lib.pg_ids).unwrap(); + let frags: Vec<(u8, u8, u8, u8, f64)> = entry + .peptide + .fragments() + .map(|f| { + ( + f.typ() & 0x7F, + f.index(), + f.charge(), + f.loss(), + f.mz() as f64, + ) + }) + .collect(); + for &(_, _, _, loss, _) in &frags { + *hist.entry(loss).or_default() += 1; + } + for &(t, i, c, loss, mz) in &frags { + if loss == 0 || c == 0 { + continue; + } + if let Some(&(_, _, _, _, base)) = frags + .iter() + .find(|&&(t2, i2, c2, l2, _)| (t2, i2, c2, l2) == (t, i, c, 0)) + { + deltas.entry(loss).or_default().push((base - mz) * c as f64); + } + } + } + // Only codes this build maps appear. A new one must be measured, not + // guessed, so it fails here rather than being dropped in silence. + for &code in hist.keys() { + assert!( + loss_from_code(code).is_some(), + "unmapped loss code {code} in the fixture: measure it before mapping it" + ); + } + // The codes carrying a loss, and the neutral mass each one measures at. + let expected = [ + (1u8, NeutralLoss::Water, 18.0106_f64), + (2, NeutralLoss::Ammonia, 17.0265), + ]; + assert_eq!( + deltas.keys().copied().collect::>(), + expected.iter().map(|e| e.0).collect::>(), + "the fixture's set of lossy codes changed" + ); + for (code, loss, mass) in expected { + assert_eq!(loss_from_code(code), Some(loss)); + let observed = &deltas[&code]; + assert!(!observed.is_empty()); + for &d in observed { + assert!( + (d - mass).abs() < 3e-3, + "code {code} measures {d:.4}, not {mass:.4} -- it is not {loss:?}" + ); + } + } + } } From 967e3c2c3e488865da69d39691ec68d400a57e96 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 12:45:53 -0700 Subject: [PATCH 7/8] chore: trim PR archaeology comments --- Taskfile.yml | 4 +- rust/micromzpaf/Cargo.toml | 5 +- rust/micromzpaf/src/error.rs | 10 +- rust/micromzpaf/src/lib.rs | 157 ++---------------- rust/micromzpaf/src/loss.rs | 101 ++--------- rust/micromzpaf/src/series.rs | 54 +----- rust/timsquery/src/serde/diann_io.rs | 6 +- rust/timsquery/src/serde/diann_speclib_io.rs | 45 +---- .../src/serde/elution_group_inputs.rs | 10 +- rust/timsquery/src/serde/skyline_io.rs | 5 +- rust/timsquery/src/serde/spectronaut_io.rs | 5 +- rust/timsseek/src/ml/cv.rs | 113 +++---------- rust/timsseek/src/ml/qvalues.rs | 90 ++-------- 13 files changed, 85 insertions(+), 520 deletions(-) diff --git a/Taskfile.yml b/Taskfile.yml index 809d37a0..f4459c3d 100644 --- a/Taskfile.yml +++ b/Taskfile.yml @@ -30,8 +30,7 @@ tasks: # dashboard.rs only exists under the opt-in feature; CI only `cargo check`s it. - cargo clippy -p timsseek_cli --features dashboard --all-targets -- -D warnings - # Broken and private intra-doc links only bite the reader, so nothing else - # catches them: CI does not build docs. + # Keep workspace documentation warning-free. doc: cmds: - RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps @@ -87,4 +86,3 @@ tasks: desc: Stop local Koina server cmds: - docker stop koina-local - diff --git a/rust/micromzpaf/Cargo.toml b/rust/micromzpaf/Cargo.toml index 168e12a8..41434f29 100644 --- a/rust/micromzpaf/Cargo.toml +++ b/rust/micromzpaf/Cargo.toml @@ -7,10 +7,9 @@ license.workspace = true [dependencies] serde = { workspace = true, features = ["derive"] } thiserror = { workspace = true } -# Only to warn once when the lossy `TryFrom<&str>` drops a mass-error suffix. +# Used for the one-time warning when parsing drops a mass-error suffix. tracing = { workspace = true } [dev-dependencies] -# Only to pin that `IonAnnot` serialises as its mzPAF string rather than the -# packed word. +# Pins the string serialization format. serde_json = { workspace = true } diff --git a/rust/micromzpaf/src/error.rs b/rust/micromzpaf/src/error.rs index 052d0035..d8b1d241 100644 --- a/rust/micromzpaf/src/error.rs +++ b/rust/micromzpaf/src/error.rs @@ -10,11 +10,6 @@ use crate::{ }; /// Why an annotation could not be parsed or represented. -/// -/// Never matched outside this crate today, so the variants exist for the message -/// a user sees when a library fails to load. Each one names the offending value -/// and, where there is one, the bound it missed -- a reader surfaces this with -/// no other context. #[derive(Debug, Error)] pub enum IonParsingError { #[error("Ordinal {ordinal} out of range for series '{series}'")] @@ -41,16 +36,13 @@ pub enum IonParsingError { UnknownIonsExhausted, #[error("Could not parse '{annotation}': {context}")] ParsingError { - /// The whole annotation, not the fragment of it that failed: a reader - /// reports this with no row index, so the full text is the only handle - /// the user gets on which row broke. annotation: String, context: &'static str, }, } impl IonParsingError { - /// Shorthand for [`Self::ParsingError`], which is built at a dozen sites. + /// Create a parsing error with context. pub(crate) fn parse(annotation: &str, context: &'static str) -> Self { Self::ParsingError { annotation: annotation.to_string(), diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 67b74833..50b9af5a 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -1,14 +1,7 @@ //! Compact representation of fragment ion annotations for mass spectrometry. //! -//! A spectral library carries one annotation per fragment, so this type is -//! replicated millions of times in a loaded arena. It is a packed `u32` so that -//! the annotations an mzSpecLib-shaped library needs -- neutral losses, -//! internal fragments, immonium ions -- fit *without* growing the type. -//! -//! Bolting a loss field onto a struct of fields would have cost three bytes, -//! not one: a 5-byte key paired with an `f32` pads to a 12-byte tuple, growing -//! the inline `TinyVec` storage in timsseek's `ExpectedIntensities` from 104 to -//! 156 bytes. Packed, the tuple stays 8 bytes. +//! A packed `u32` representation of the fragment annotations used by spectral +//! libraries, including neutral losses, internal fragments and immonium ions. //! //! # Bit layout //! @@ -20,8 +13,7 @@ //! └───────┴────────────────┴──────────┴────────┴────────┴───────┘ //! ``` //! -//! The widths and their shifts are checked by `const` assertions next to the -//! constants, so this diagram cannot drift away from the code. +//! The widths and shifts are checked by `const` assertions. //! //! `payload` is reinterpreted per `kind` -- a tagged union inside the word: //! @@ -33,13 +25,10 @@ //! | precursor | unused | //! //! `unknown` is discriminant 0 and `charge` is stored biased by one, so the -//! all-zero word is the valid annotation `?0` at charge 1. That matters -//! because `IonAnnot: Default` is forced by `tinyvec::Array` and a default can -//! reach any serde path. +//! all-zero word is the valid annotation `?0` at charge 1. //! -//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned -//! (see [`ISOTOPE_MIN`]). A bit field truncates rather than wrapping loudly, so -//! every constructor range-checks -- see [`IonAnnot::try_new`]. +//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned. +//! Constructors validate both before packing. //! //! # The mzPAF subset //! @@ -111,11 +100,6 @@ use series::{ use std::fmt::Display; use std::hash::Hash; -// ── The word ───────────────────────────────────────────────────────────────── -// -// `kind` and `payload` are `series`'s to define; this module only decides where -// they sit and what surrounds them. - const KIND_SHIFT: u32 = 0; const CHARGE_SHIFT: u32 = 4; const CHARGE_BITS: u32 = 4; @@ -125,28 +109,13 @@ const LOSS_SHIFT: u32 = 12; const LOSS_BITS: u32 = 6; const PAYLOAD_SHIFT: u32 = 18; -/// Charge bounds. Observed maximum in the HUPO-PSI corpus is 3. -/// -/// Asymmetric because the field stores `charge - 1` zigzagged: that shifts the -/// whole window up by one, so `8` fits and `-8` does not. The `const` assertion -/// below is what holds these to the field rather than to this comment -- -/// `CHARGE_MIN = -8` would zigzag to 17 and truncate to 1, silently decoding as -/// charge 0. +/// Charge bounds representable by the packed field. pub const CHARGE_MIN: i8 = -7; pub const CHARGE_MAX: i8 = 8; -/// Widest isotope offset the 4-bit field holds. Observed maximum is 3. -/// -/// Isotope offsets are unsigned. mzPAF spells a negative offset `-Ni`, which -/// this crate does not parse, so allowing one in the field would let `Display` -/// emit `y5+-2i` -- text no other mzPAF reader accepts, through a type whose -/// wire format *is* that text. Rejecting it here keeps the invalid spelling -/// unreachable and spends the whole 4 bits on offsets that can round-trip. +/// Isotope offset bounds representable by the packed field. pub const ISOTOPE_MIN: i8 = 0; pub const ISOTOPE_MAX: i8 = mask(ISOTOPE_BITS) as i8; -// The layout is otherwise enforced by prose and a diagram. These make it -// self-checking, so a field cannot be widened or moved without the build -// failing. `series` asserts its own payload bounds. const _: () = assert!( KIND_BITS + CHARGE_BITS + ISOTOPE_BITS + LOSS_BITS + PAYLOAD_BITS <= 32, "the fields overflow the word" @@ -162,7 +131,6 @@ const _: () = assert!( "fields must abut" ); const _: () = assert!(PAYLOAD_SHIFT == LOSS_SHIFT + LOSS_BITS, "fields must abut"); -// The zigzag bounds are hand-derived; these are what make them checked. const _: () = assert!( zigzag_charge(CHARGE_MIN) <= mask(CHARGE_BITS) && zigzag_charge(CHARGE_MAX) <= mask(CHARGE_BITS), @@ -172,8 +140,6 @@ const _: () = assert!( ISOTOPE_MIN >= 0 && ISOTOPE_MAX as u32 <= mask(ISOTOPE_BITS), "the isotope range does not fit ISOTOPE_BITS" ); -// `pack` masks the loss discriminant, so a table grown past the field would -// truncate into a *different* loss rather than fail. const _: () = assert!( NeutralLoss::COUNT as u32 <= mask(LOSS_BITS), "the loss table has outgrown LOSS_BITS" @@ -184,8 +150,7 @@ pub(crate) const fn mask(bits: u32) -> u32 { (1u32 << bits) - 1 } -/// Zigzag: map a small signed value onto an unsigned one without losing the -/// sign bit to the field width. +/// Map a small signed value to an unsigned value without losing its sign. #[inline] const fn zigzag(v: i8) -> u32 { (((v as i32) << 1) ^ ((v as i32) >> 31)) as u32 @@ -195,11 +160,7 @@ const fn unzigzag(u: u32) -> i8 { (((u >> 1) as i32) ^ -((u & 1) as i32)) as i8 } -/// Charge is stored biased by one, so the zero field decodes to charge 1. -/// -/// Charge 0 is rejected by every constructor, so it is not a value the field -/// needs to represent -- and spending the zero word on it would make -/// `IonAnnot::default()` render an annotation that cannot be parsed back. +/// Encode charge with a one-unit bias so zero decodes to charge 1. #[inline] const fn zigzag_charge(charge: i8) -> u32 { zigzag(charge - 1) @@ -212,11 +173,6 @@ const fn unzigzag_charge(u: u32) -> i8 { /// Compact representation of fragment annotations. /// /// A packed `u32`; see the crate docs for the bit layout. -/// -/// Deliberately not `Ord`. Ordering the packed word sorts by ordinal first, -/// then loss, then isotope, then charge, then series -- an order nobody means. -/// Nothing in the workspace sorts annotations, and `KeyLike` does not require -/// it, so the trait is not offered rather than offered and meaningless. #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, Default)] pub struct IonAnnot(u32); @@ -244,13 +200,7 @@ impl<'de> Deserialize<'de> for IonAnnot { } impl IonAnnot { - /// Build an annotation from its parts. - /// - /// The one constructor. `series` has already been validated against the - /// payload width by [`IonSeriesOrdinal`]'s own constructors, so all this - /// checks is `charge` and `isotope` -- and it must, because a bit field - /// truncates silently, so an unchecked value would corrupt the annotation - /// rather than fail. + /// Build an annotation from its parts, validating the packed fields. pub fn new( series: IonSeriesOrdinal, loss: NeutralLoss, @@ -276,12 +226,7 @@ impl IonAnnot { )) } - /// A backbone / precursor / unknown annotation from its mzPAF letter. - /// - /// The shape a file reader has: a series character straight out of a column. - /// Internal fragments and immonium ions are spelled differently and go - /// through [`IonSeriesOrdinal::try_internal`] / [`IonSeriesOrdinal::try_immonium`] - /// and [`Self::new`]. + /// Build a backbone, precursor or unknown annotation from its series letter. pub fn try_new( ion_type: char, ordinal: Option, @@ -323,9 +268,6 @@ impl IonAnnot { /// Errors when the result leaves [`ISOTOPE_MIN`]..=[`ISOTOPE_MAX`], naming /// that range, since the caller cannot see the field width. pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { - // Saturating rather than checked: the field is far narrower than `i8`, - // so the range check below is the one that matters and it reports the - // bound that actually applies. let new_isotope = self.get_isotope().saturating_add(offset_neutrons); if !(ISOTOPE_MIN..=ISOTOPE_MAX).contains(&new_isotope) { return Err(IonParsingError::IsotopeOutOfRange { @@ -373,18 +315,9 @@ static MASS_ERROR_DISCARDED: std::sync::Once = std::sync::Once::new(); impl TryFrom<&str> for IonAnnot { type Error = IonParsingError; - /// Parses an annotation, **discarding any mass-error suffix**, and warns - /// once per process the first time it discards a real one. - /// - /// The suffix is a property of one observed peak, not of the ion: two peaks - /// annotated `b12` in different spectra carry different errors, so keeping - /// it on the annotation would make two `b12`s unequal and break the - /// per-precursor label uniqueness the whole crate keys on. + /// Parse an annotation, ignoring any mass-error suffix. /// - /// A caller that needs the error wants it for the *m/z*, not the label: - /// call [`split_mass_error`] first, recover the theoretical m/z with - /// [`MassError::theoretical_from_observed`], and store that. Nothing is - /// lost that way -- only the residual, which belongs to the measurement. + /// Use [`split_mass_error`] when the suffix is needed for m/z correction. fn try_from(value: &str) -> Result { let (ion, mass_error) = split_mass_error(value)?; if mass_error.is_some() { @@ -425,18 +358,7 @@ impl Display for IonAnnot { } } -/// Hands out `?1`, `?2`, ... for peaks whose annotation this crate cannot -/// represent. -/// -/// Fragment labels must be unique within a precursor -- lookup is by first -/// match, so a repeated label makes every later peak carrying it unreachable. -/// A monotonic counter makes that uniqueness structural, and returning an -/// error once the 8-bit ordinal is spent keeps the overflow from being -/// something each reader has to remember to check. -/// -/// Deliberately neither `Copy` nor `Clone`: a duplicated counter forks, and -/// each fork reissues labels the other already handed out -- exactly the bug -/// this type exists to prevent. Pass it by `&mut`. +/// Generates unique unknown-ion labels within a precursor. #[derive(Debug, Default)] pub struct UnknownIonCounter(u8); @@ -461,23 +383,12 @@ mod tests { IonAnnot::try_from(s).unwrap_or_else(|e| panic!("{s:?} must parse: {e}")) } - /// The size claim the design rests on: a loss, an internal-fragment span - /// and an immonium residue all ride along without growing the tuple that - /// timsseek stores inline. - /// - /// The predecessor was also 4 bytes, so `size_of` alone pins nothing -- - /// what this asserts is that the *added* fields cost nothing, which is only - /// meaningful together with the round-trip tests proving they are really in - /// there. + /// Keep the annotation and annotation/intensity pair compact. #[test] fn the_added_fields_cost_no_space() { assert_eq!(size_of::(), 4); - // Paired with an intensity on the scoring hot path; padding here would - // grow the inline TinyVec storage in `ExpectedIntensities` from 104 to - // 156 bytes at the current inline capacity of 13. assert_eq!(size_of::<(IonAnnot, f32)>(), 8); - // All four of these are new capacity, and none of them widened the word. let loaded = IonAnnot::new( IonSeriesOrdinal::try_internal(2, 11).expect("in range"), NeutralLoss::PhosphoricAcidWater, @@ -494,15 +405,7 @@ mod tests { assert_eq!((loaded.get_charge(), loaded.get_isotope()), (3, 2)); } - /// `IonAnnot: Default` is not optional -- `tinyvec::Array` requires - /// `Item: Default`, `TimsElutionGroup` stores labels in a `TinyVec`, and - /// timsquery's `KeyLike` propagates the bound. So a default can reach any - /// serde path, and the zero word has to mean something. - /// - /// It means `?0`: charge is stored as `zigzag(charge - 1)`, so the zero - /// field is charge 1 rather than the impossible charge 0. The default is - /// therefore a real annotation and `Serialize`/`Deserialize` are inverses - /// on it, instead of rendering a value no constructor accepts. + /// The default packed word is a valid annotation and round-trips via serde. #[test] fn the_default_annotation_is_a_real_annotation() { let d = IonAnnot::default(); @@ -532,7 +435,6 @@ mod tests { assert_eq!(annot, expected, "{input}"); assert_eq!(annot.get_charge(), charge, "{input}"); assert_eq!(annot.get_isotope(), isotope, "{input}"); - // Round-trips byte-identically when no loss is involved. assert_eq!(format!("{}", annot), input); } } @@ -577,20 +479,12 @@ mod tests { Err(IonParsingError::IsotopeOutOfRange { .. }) )); assert!(IonAnnot::try_new('y', Some(1), 0, 0).is_err()); - // The payload bounds are `series`'s to enforce; see - // `constructors_reject_what_the_payload_cannot_hold` there. This pins - // that an out-of-range span cannot reach a word through the parser. assert!(matches!( IonAnnot::try_from("m64:1"), Err(IonParsingError::OrdinalOutOfRange { .. }) )); } - /// The crate docs say negative isotope offsets are unsupported. When the - /// field held them, `Display` emitted `y5+-2i` -- not mzPAF, and reparsed - /// happily, so a library could round-trip through serde into text no other - /// mzPAF reader accepts. The field is unsigned so that spelling is - /// unreachable rather than merely undocumented. #[test] fn negative_isotopes_are_unrepresentable_not_just_unparsed() { assert!(matches!( @@ -598,11 +492,8 @@ mod tests { Err(IonParsingError::IsotopeOutOfRange { isotope: -1 }) )); assert!(ion("y5").try_with_offset_neutrons(-1).is_err()); - // mzPAF's own spelling for a negative offset is still not parsed. assert!(IonAnnot::try_from("y5-2i").is_err()); - // And the spelling that used to leak out is not accepted either. assert!(IonAnnot::try_from("y5+-2i").is_err()); - // Every representable isotope renders as something that parses back. for isotope in ISOTOPE_MIN..=ISOTOPE_MAX { let a = IonAnnot::try_new('y', Some(5), 1, isotope).expect("in range"); let text = a.to_string(); @@ -625,13 +516,9 @@ mod tests { assert_eq!(a.try_get_ordinal(), Some(5)); assert_eq!(format!("{}", a), "y5-H2O"); - // Non-canonical spelling resolves to the same annotation, and renders - // canonically -- so this pair is equal, which is the property that - // upholds per-precursor label uniqueness. assert_eq!(ion("y5-CH3SOH"), ion("y5-CH4OS")); assert_eq!(format!("{}", ion("y5-CH3SOH")), "y5-CH4OS"); - // Loss combines with charge and isotope. let b = ion("y10-NH3+i^2"); assert_eq!(b.loss(), NeutralLoss::Ammonia); assert_eq!(b.get_isotope(), 1); @@ -645,7 +532,6 @@ mod tests { a.series_ordinal(), IonSeriesOrdinal::internal { start: 2, end: 11 } ); - // An internal fragment has no ladder position. assert_eq!(a.try_get_ordinal(), None); assert_eq!(format!("{}", a), "m2:11"); @@ -663,8 +549,6 @@ mod tests { ); assert_eq!(format!("{}", a), "IA"); - // Carries an arbitrary mod string -- no field width represents it, so - // it must fail rather than silently degrade to a bare immonium. assert!(matches!( IonAnnot::try_from("IC[Carbamidomethyl]"), Err(IonParsingError::UnsupportedModifiedImmonium { .. }) @@ -676,8 +560,6 @@ mod tests { let (rest, err) = split_mass_error("y1/-0.0005").unwrap(); assert_eq!(ion(rest), ion("y1")); assert_eq!(err, Some(MassError::Da(-0.0005))); - // theoretical = observed - error; verified against a real SpectraST - // peak: y1 for C-terminal R, observed 175.1184, theoretical 175.11895. let theo = err.unwrap().theoretical_from_observed(175.1184); assert!((theo - 175.1189).abs() < 1e-9, "got {theo}"); @@ -686,15 +568,10 @@ mod tests { let theo = err.unwrap().theoretical_from_observed(700.0); assert!((theo - 699.99916).abs() < 1e-4, "got {theo}"); - // Absent suffix is not an error. assert_eq!(split_mass_error("y6").unwrap(), ("y6", None)); - // TryFrom discards it rather than failing. assert_eq!(IonAnnot::try_from("y1/-0.0005").unwrap(), ion("y1")); } - /// An unrepresentable loss must fail loudly. Parsing `y1-HCOOH` as plain - /// `y1` would put a loss peak's m/z on the `y1` label and collide with the - /// real `y1`. #[test] fn unrepresentable_loss_is_rejected_not_stripped() { assert!(matches!( diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs index 7009372b..ddd7cdaa 100644 --- a/rust/micromzpaf/src/loss.rs +++ b/rust/micromzpaf/src/loss.rs @@ -1,7 +1,6 @@ //! Neutral losses, keyed by atomic composition rather than by spelling. //! -//! Libraries write the same chemical loss different ways. Two real examples -//! from the HUPO-PSI mzSpecLib corpus: +//! Libraries may write the same chemical loss in different ways: //! //! | written | library | composition | //! |---|---|---| @@ -10,14 +9,8 @@ //! | `-NH2-CO-CH2SH` | NIST | C2H5N1O1S1 | //! | `-C2H5NOS` | SpectraST | C2H5N1O1S1 | //! -//! Keying on the string would make `y5-CH4OS` and `y5-CH3SOH` distinct labels -//! for one ion. Since fragment labels must be unique within a precursor (see -//! `ExpectedIntensities::try_from_pairs` in timsseek), that is exactly the -//! wrong direction: it hides a genuine duplicate behind two spellings. -//! -//! So parsing goes `text -> composition -> discriminant`, and the composition -//! is a *parse-time* concept only. What gets stored on an `IonAnnot` is the -//! discriminant, so there is no per-annotation cost at runtime. +//! Parsing goes `text -> composition -> discriminant`; the packed annotation +//! stores only the discriminant. use std::fmt::Display; @@ -25,14 +18,7 @@ use crate::IonParsingError; /// Atom counts for the elements that appear in peptide neutral losses. /// -/// Deliberately not a general chemical formula: these losses only ever draw -/// from C/H/N/O/S/P, and six named `u8`s make equality a 6-byte compare during -/// the parse-time table lookup. -/// -/// Named fields rather than `[u8; 6]` so that [`TABLE`] reads as chemistry -/// (`H2O` is `h: 2, o: 1`) instead of six positional numbers, where a -/// transposition would be invisible on review and would silently alias one -/// loss onto another. +/// These losses use only C/H/N/O/S/P. #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] pub(crate) struct Composition { c: u8, @@ -44,10 +30,6 @@ pub(crate) struct Composition { } /// A [`Composition`] naming only the elements it contains: `C!(h: 2, o: 1)`. -/// -/// A plain struct literal would have to spell all six counts, which is the -/// positional noise this struct exists to remove; functional update syntax -/// (`..ZERO`) is not permitted in a `const` item. macro_rules! C { ($($field:ident: $count:expr),+ $(,)?) => { Composition { $($field: $count,)+ ..Composition::ZERO } @@ -68,8 +50,6 @@ impl Composition { /// The count for one element symbol, or `None` if this crate does not /// represent that element. /// - /// The single place the symbol-to-field mapping lives, so the parser cannot - /// disagree with the struct about which letter means which count. fn count_mut(&mut self, symbol: u8) -> Option<&mut u8> { Some(match symbol { b'C' => &mut self.c, @@ -88,10 +68,7 @@ impl Composition { [self.c, self.h, self.n, self.o, self.s, self.p] } - /// Combine element-wise. Saturating: the only inputs that reach the ceiling - /// are absurd (`200H2O`), and no [`TABLE`] row holds a saturated count, so a - /// saturated result cannot alias onto a real loss. Pinned by - /// `table_compositions_are_unique`. + /// Combine element-wise, saturating each count. fn zip(self, other: Self, f: impl Fn(u8, u8) -> u8) -> Self { Self { c: f(self.c, other.c), @@ -114,9 +91,7 @@ impl Composition { /// Parse a bare formula like `H2O`, `CH4OS`, `C2H5NOS`. /// - /// Only single-letter C/H/N/O/S/P are recognized; anything else is an - /// error rather than a silent skip, so an unsupported loss surfaces as - /// "not representable" instead of being mistaken for a smaller one. + /// Only single-letter C/H/N/O/S/P elements are recognized. fn parse_formula(s: &str) -> Result { if s.is_empty() { return Err(IonParsingError::parse(s, "Empty neutral-loss formula")); @@ -149,8 +124,7 @@ impl Composition { /// Parse a full loss expression: `-` separated terms, each optionally /// prefixed by a repeat count. `2H2O`, `H2O-NH3`, `NH2-CO-CH2SH`. /// - /// Because terms are summed, ordering and multiplier spelling collapse for - /// free: `H2O-NH3` == `NH3-H2O`, and `2H2O` == `H2O-H2O`. + /// Terms are summed, so ordering and multiplier spelling are normalized. pub(crate) fn parse_expression(s: &str) -> Result { let mut total = Composition::default(); for term in s.split('-') { @@ -181,11 +155,7 @@ impl Composition { /// The neutral losses this crate can represent, as a packed discriminant. /// -/// Scoped deliberately: DIA-NN emits none, Spectronaut two (`-H2O`, `-NH3`), -/// NIST eight, plus the phospho losses that no non-phospho corpus can show. -/// SpectraST's wider combinatorics are NOT here -- an unlisted loss parses to a -/// composition that misses the table and is reported as unrepresentable, which -/// routes the peak to an unknown label rather than silently mislabelling it. +/// Inputs outside this table are reported as unrepresentable. #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] #[repr(u8)] pub enum NeutralLoss { @@ -218,17 +188,8 @@ pub enum NeutralLoss { } /// `(composition, discriminant, canonical spelling)`, indexed by discriminant. -/// -/// Rows must stay in discriminant order starting at 1 -- -/// `table_is_indexed_by_discriminant` pins that, and both `from_discriminant` -/// and `canonical` index straight into this. -/// -/// The canonical spelling is what `Display` emits, so a non-canonical input -/// (`-CH3SOH`) round-trips to the canonical form (`-CH4OS`). Round-trip tests -/// must therefore compare parsed values, not bytes. -/// -/// Compositions name their elements, so each row can be read against its own -/// spelling. +/// Entries are ordered by discriminant; canonical spelling is emitted by +/// [`Display`]. const TABLE: &[(Composition, NeutralLoss, &str)] = &[ (C!(h: 2, o: 1), NeutralLoss::Water, "H2O"), (C!(h: 3, n: 1), NeutralLoss::Ammonia, "NH3"), @@ -257,15 +218,10 @@ const TABLE: &[(Composition, NeutralLoss, &str)] = &[ ]; impl NeutralLoss { - /// How many losses the table names, so the width of the packed field can be - /// checked against it rather than assumed to be roomy. + /// Number of loss entries in the table. pub(crate) const COUNT: u8 = TABLE.len() as u8; - /// Inverse of the `#[repr(u8)]` discriminant, for unpacking out of a bit - /// field. Lives next to the enum so the two cannot drift apart. - /// - /// An unrecognized value maps to [`Self::None`]: the only way to produce - /// one is a reserved discriminant, which no constructor emits. + /// Decode a packed loss discriminant. pub(crate) fn from_discriminant(d: u8) -> Self { Self::at(d).map_or(Self::None, |(_, loss, _)| *loss) } @@ -273,20 +229,13 @@ impl NeutralLoss { /// The [`TABLE`] row a discriminant names, or `None` for [`Self::None`] and /// for the reserved values above the table. /// - /// `TABLE` is in discriminant order starting at 1, which - /// `table_is_indexed_by_discriminant` pins. Indexing it rather than - /// re-listing the variants is what keeps the discriminant, the composition - /// and the spelling from drifting apart. fn at(d: u8) -> Option<&'static (Composition, NeutralLoss, &'static str)> { TABLE.get((d as usize).checked_sub(1)?) } /// Resolve a loss expression (without the leading `-`) to a discriminant. /// - /// `None` means "parsed as a valid composition, but not one we represent" -- - /// distinct from `Err`, which means the text was not a loss expression at - /// all. Callers route the former to an unknown label and the latter to a - /// parse failure. + /// `None` means the composition is valid but not in the supported table. pub(crate) fn from_expression(s: &str) -> Result, IonParsingError> { let comp = Composition::parse_expression(s)?; Ok(TABLE @@ -314,8 +263,6 @@ impl Display for NeutralLoss { mod tests { use super::*; - /// The count parser: implicit 1, multi-digit counts, and a count past the - /// `u8` slot. Single-digit cases are covered by every TABLE row. #[test] fn atom_counts_parse_and_reject_out_of_range() { assert_eq!( @@ -339,10 +286,8 @@ mod tests { ); } - /// The two cross-library spelling collisions this module exists for. #[test] fn different_spellings_resolve_to_one_loss() { - // NIST vs SpectraST, methanesulfenic acid. assert_eq!( NeutralLoss::from_expression("CH3SOH").unwrap(), Some(NeutralLoss::Methanesulfenic) @@ -351,7 +296,6 @@ mod tests { NeutralLoss::from_expression("CH4OS").unwrap(), Some(NeutralLoss::Methanesulfenic) ); - // NIST structural notation vs SpectraST molecular notation. assert_eq!( NeutralLoss::from_expression("NH2-CO-CH2SH").unwrap(), Some(NeutralLoss::Carbamidomethylthiol) @@ -362,7 +306,6 @@ mod tests { ); } - /// Summing terms collapses ordering and multiplier spelling for free. #[test] fn ordering_and_multipliers_normalize() { assert_eq!( @@ -379,10 +322,6 @@ mod tests { ); } - /// A well-formed composition outside the table is `Ok(None)` -- "valid but - /// not representable" -- while malformed text is `Err`. Callers need to - /// tell those apart to route one to an unknown label and the other to a - /// parse failure. #[test] fn unrepresentable_is_distinct_from_malformed() { assert_eq!(NeutralLoss::from_expression("HCOOH").unwrap(), None); @@ -391,10 +330,6 @@ mod tests { assert!(NeutralLoss::from_expression("H2O-").is_err()); } - /// `from_discriminant` and `canonical` both index [`TABLE`] by discriminant, - /// so a row sitting at the wrong offset would decode as its neighbour. This - /// is the one invariant that keeps the enum, the composition and the - /// spelling in step. #[test] fn table_is_indexed_by_discriminant() { for (i, (_comp, loss, canon)) in TABLE.iter().enumerate() { @@ -414,8 +349,6 @@ mod tests { assert_eq!(NeutralLoss::None.canonical(), ""); } - /// Every table entry must survive canonical -> composition -> discriminant, - /// and back out through the bit field. #[test] fn table_round_trips_through_canonical_spelling() { for (_comp, loss, canon) in TABLE { @@ -433,18 +366,12 @@ mod tests { } } - /// Compositions must be unique: two entries sharing one would make the - /// table lookup order-dependent. #[test] fn table_compositions_are_unique() { for (i, (a, _, sa)) in TABLE.iter().enumerate() { for (b, _, sb) in TABLE.iter().skip(i + 1) { assert_ne!(a, b, "{sa} and {sb} share a composition"); } - // `scaled` and `plus` saturate, so an absurd input like `200H2O` - // lands on 255 in some slot. That is only safe to report as - // unrepresentable while no table entry holds a saturated count -- - // otherwise the saturation would alias onto a real loss. assert!( a.counts().iter().all(|&n| n < u8::MAX), "{sa} holds a saturated atom count" @@ -452,8 +379,6 @@ mod tests { } } - /// The two things the TABLE round trip cannot see: the `-` prefix Display - /// adds, and that `None` renders as nothing at all rather than "-". #[test] fn display_adds_the_prefix_and_none_renders_empty() { assert_eq!(NeutralLoss::None.to_string(), ""); diff --git a/rust/micromzpaf/src/series.rs b/rust/micromzpaf/src/series.rs index f3834ba8..025b1ae6 100644 --- a/rust/micromzpaf/src/series.rs +++ b/rust/micromzpaf/src/series.rs @@ -1,7 +1,6 @@ //! What kind of ion an annotation names, and how that fits in a packed payload. //! -//! This module owns three things that have to agree, and only agree if they sit -//! together: +//! This module owns three related details: //! //! 1. the mzPAF **spelling** of each kind (`Display` and `IonSeriesOrdinal::parse`, //! which are exact inverses), @@ -9,11 +8,8 @@ //! and `from_parts`, also inverses), //! 3. the **bounds** that layout imposes, enforced by the constructors. //! -//! Splitting any of those across a module boundary is what lets a `v` render as -//! a `w`, or an internal-fragment endpoint get validated against a width it no -//! longer has. [`IonAnnot`](crate::IonAnnot) owns the surrounding word -- -//! charge, isotope, loss, and where these two fields sit in it -- and nothing -//! else. +//! [`IonAnnot`](crate::IonAnnot) owns the surrounding word -- charge, isotope, +//! loss, and their positions. use std::fmt::Display; @@ -83,8 +79,7 @@ pub enum Series { } impl Series { - /// Every series, in discriminant order, and parallel to the private letter - /// table -- `series_letters_and_discriminants_agree` pins that pairing. + /// Every series, in discriminant order. pub const ALL: [Self; 9] = [ Self::a, Self::b, @@ -96,8 +91,7 @@ impl Series { Self::y, Self::z, ]; - /// The mzPAF letters, in discriminant order. The single place the - /// letter-to-discriminant pairing lives. + /// The mzPAF letters, in discriminant order. const CHARS: &'static [u8; 9] = b"abcdvwxyz"; /// The mzPAF letter for this series. @@ -202,9 +196,7 @@ impl IonSeriesOrdinal { /// Parse the ion-kind part of an annotation: everything left of the charge, /// isotope and loss suffixes. /// - /// The exact inverse of this type's [`Display`], which is why the two live - /// side by side -- `every_series_variant_round_trips_through_its_mzpaf_spelling` - /// pins that. + /// The inverse of this type's [`Display`]. pub(crate) fn parse(core: &str) -> Result { // Internal fragment: `m:`. Checked before the single-letter // forms because `m` is not a backbone letter, so there is no ambiguity. @@ -256,9 +248,7 @@ impl IonSeriesOrdinal { /// Split into the `kind` discriminant and its `payload`, the two fields /// [`IonAnnot`](crate::IonAnnot) packs. /// - /// Total: every arm masks to its own width, so a value built by hand rather - /// than through the constructors encodes as something decodable instead of - /// panicking on a path `Display` (and so `Serialize`) reaches. + /// Split into the fields used by the packed word. pub(crate) const fn to_parts(self) -> (u32, u32) { match self { Self::backbone { series, ordinal } => (series as u32, ordinal as u32), @@ -279,9 +269,8 @@ impl IonSeriesOrdinal { } } - /// Inverse of [`Self::to_parts`]. Total by construction: `unknown` is the - /// catch-all discriminant, so a value this build does not recognise -- only - /// reachable from a corrupted word -- decodes as an unknown ion. + /// Inverse of [`Self::to_parts`]. Unknown discriminants decode as unknown + /// ions. pub(crate) const fn from_parts(kind: u32, payload: u32) -> Self { let ordinal = payload as u8; match kind { @@ -319,10 +308,6 @@ impl Display for IonSeriesOrdinal { mod tests { use super::*; - /// One representative of every [`IonSeriesOrdinal`] case: all nine backbone - /// series (each with a distinct ordinal, so a transposition in - /// `to_parts`/`from_parts` cannot cancel out) plus the four others, at their - /// field boundaries. fn all_series() -> Vec { let mut out: Vec = Series::ALL .iter() @@ -336,22 +321,16 @@ mod tests { IonSeriesOrdinal::precursor, IonSeriesOrdinal::unknown { ordinal: 10 }, IonSeriesOrdinal::internal { start: 2, end: 11 }, - // Both endpoints at their 6-bit ceiling: the widest payload the - // 12 bits hold, and the case a narrowed field would silently clip. IonSeriesOrdinal::internal { start: INTERNAL_POS_MAX, end: INTERNAL_POS_MAX, }, IonSeriesOrdinal::immonium { residue: 'A' }, - // Highest residue index the 5-bit immonium field must hold. IonSeriesOrdinal::immonium { residue: 'Z' }, ]); out } - /// `to_parts` and `from_parts` are hand-written inverses. Without this, - /// swapping two arms (`v` encoding as `w`) mislabels a whole ion series and - /// every other test still passes. #[test] fn every_variant_round_trips_through_the_packed_parts() { for series in all_series() { @@ -364,8 +343,6 @@ mod tests { ); } - // Every backbone series must land on its own discriminant; the other - // four are singletons and are covered by the round trip above. let mut kinds: Vec = Series::ALL .iter() .map(|&series| { @@ -379,23 +356,15 @@ mod tests { assert_eq!(kinds.len(), Series::ALL.len(), "two series share a kind"); } - /// `from_parts` claims to be total, and three of the sixteen `kind` values - /// are unassigned. Renumbering the assigned ones is safe only while the - /// unassigned ones stay decodable, since a corrupted word reaches this on a - /// `Display` (and so `Serialize`) path. #[test] fn every_kind_bit_pattern_decodes_without_panicking() { for kind in 0..=mask(KIND_BITS) { for payload in [0, 1, mask(PAYLOAD_BITS)] { - // Also exercises `Display`, which is where a partial decode - // would have panicked. let _ = IonSeriesOrdinal::from_parts(kind, payload).to_string(); } } } - /// `Display` and [`IonSeriesOrdinal::parse`] are hand-written inverses, and - /// this module exists so they cannot drift apart. #[test] fn every_variant_round_trips_through_its_mzpaf_spelling() { for series in all_series() { @@ -408,8 +377,6 @@ mod tests { } } - /// The letters are the mzPAF spelling of the discriminants, and - /// `from_char`/`as_char` are hand-written inverses of each other. #[test] fn series_letters_and_discriminants_agree() { assert_eq!(Series::ALL.len(), Series::CHARS.len()); @@ -421,8 +388,6 @@ mod tests { } } - /// The bounds this module owns, checked where they are enforced rather than - /// through a packed word. #[test] fn constructors_reject_what_the_payload_cannot_hold() { assert!(IonSeriesOrdinal::try_internal(INTERNAL_POS_MAX, INTERNAL_POS_MAX).is_ok()); @@ -432,7 +397,6 @@ mod tests { )); assert!(IonSeriesOrdinal::try_immonium('A').is_ok()); assert!(IonSeriesOrdinal::try_immonium('Z').is_ok()); - // Lowercase would underflow the `- b'A'` in `to_parts`. assert!(IonSeriesOrdinal::try_immonium('a').is_err()); // A bare `I` names no residue; a modified one carries a mod string. diff --git a/rust/timsquery/src/serde/diann_io.rs b/rust/timsquery/src/serde/diann_io.rs index b979348a..02251d65 100644 --- a/rust/timsquery/src/serde/diann_io.rs +++ b/rust/timsquery/src/serde/diann_io.rs @@ -50,9 +50,7 @@ impl std::error::Error for DiannReadingError {} #[derive(Debug)] pub enum DiannPrecursorParsingError { - /// Which fragment row failed, and why. The row index is the only handle the - /// user gets on a bad row: one failure aborts the whole library load, and - /// nothing upstream knows where it happened. + /// Ion parsing failed at this fragment row. IonParsing { row: usize, source: IonParsingError, @@ -83,8 +81,6 @@ impl std::fmt::Display for DiannPrecursorParsingError { } impl DiannPrecursorParsingError { - /// `map_err` adaptor that stamps the fragment row onto an ion-parsing - /// failure. Replaces a `From` impl, which had no way to see the row. fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { move |source| Self::IonParsing { row, source } } diff --git a/rust/timsquery/src/serde/diann_speclib_io.rs b/rust/timsquery/src/serde/diann_speclib_io.rs index fa7c1e20..4659857b 100644 --- a/rust/timsquery/src/serde/diann_speclib_io.rs +++ b/rust/timsquery/src/serde/diann_speclib_io.rs @@ -8,8 +8,7 @@ //! (little-endian, no padding, no section table); it must be parsed strictly in //! order, so there is no random access across sections. //! -//! The reader is three layers. The types named here are private to this module, -//! so they are spelled rather than linked: +//! The reader is three layers: //! 1. Decode (`Cursor`, `Fragment`, `Peptide`): typed, zero-copy views //! over byte ranges of the in-memory file. No domain knowledge. //! 2. Emit (`SpecLib` + `EntryIter`): a pull parser that walks the @@ -707,19 +706,8 @@ fn append_arena(dst: &mut TargetColumns, mut src: TargetColumns Option { match code { 0 => Some(NeutralLoss::None), @@ -843,11 +831,7 @@ fn map_entry( continue; } - // `charge` is a raw wire byte. `as i8` would reinterpret anything above - // 127 as negative and hand a plausible-looking charge to the - // constructor, so it is range-checked instead. Note the sibling `typ()` - // is masked `& 0x7F` because DIA-NN sets a flag bit there; if this byte - // carries flags too, this is where it will surface (see #105). + // The wire value is unsigned; reject values that do not fit IonAnnot. let charge = match i8::try_from(f.charge()) { Ok(charge) => charge, Err(_) => { @@ -1190,15 +1174,9 @@ mod tests { let file = std::fs::File::open(fixture_path()).unwrap(); let (geom, _intens, stats, _eof) = parse_speclib_reader(std::io::BufReader::new(file)).unwrap(); - // Reference parser: 8384 ExcludeFromAssay-flagged (kept, counted only), - // no dc != 0 entries. assert_eq!(stats.exclude_flagged, 8384); assert_eq!(stats.decoys_dropped, 0); - // Every neutral loss in this fixture is water or ammonia, both of which - // `loss_from_code` maps, so nothing is dropped for its loss. This was - // 152 before that mapping existed. assert_eq!(stats.loss_dropped, 0); - // And they arrive as labelled losses rather than as bare ions. let with_loss = geom .frag_labels .iter() @@ -1215,18 +1193,6 @@ mod tests { ); } - /// Re-derives [`loss_from_code`] from the fixture instead of trusting it. - /// - /// DIA-NN's loss numbering is undocumented, so the mapping was measured: a - /// lossy fragment and its no-loss sibling differ by exactly the loss, so - /// `(base_mz - lossy_mz) * charge` is the neutral mass. If a future fixture - /// or DIA-NN version renumbers the codes, this fails here rather than - /// silently relabelling peaks. - /// - /// The tolerance is 3 mDa: the m/z values on disk are `f32`, and DIA-NN's - /// own loss constants are a shade off the monoisotopic masses (water reads - /// 18.011 against 18.0106). That is far tighter than the gap to any other - /// candidate loss, which is what makes the identification safe. #[test] fn loss_codes_are_water_and_ammonia() { use std::collections::BTreeMap; @@ -1272,15 +1238,12 @@ mod tests { } } } - // Only codes this build maps appear. A new one must be measured, not - // guessed, so it fails here rather than being dropped in silence. for &code in hist.keys() { assert!( loss_from_code(code).is_some(), "unmapped loss code {code} in the fixture: measure it before mapping it" ); } - // The codes carrying a loss, and the neutral mass each one measures at. let expected = [ (1u8, NeutralLoss::Water, 18.0106_f64), (2, NeutralLoss::Ammonia, 17.0265), diff --git a/rust/timsquery/src/serde/elution_group_inputs.rs b/rust/timsquery/src/serde/elution_group_inputs.rs index 0c05675d..53f5eee8 100644 --- a/rust/timsquery/src/serde/elution_group_inputs.rs +++ b/rust/timsquery/src/serde/elution_group_inputs.rs @@ -22,9 +22,7 @@ pub enum ElutionGroupInputError { inner: String, }, MissingFragmentLabels, - /// More fragments than the `u8` label space holds. Wrapping instead would - /// mint a duplicate label, and lookup is by first match, so every later - /// fragment sharing it would be unreachable. + /// Too many fragments for the unknown-ion label space. TooManyFragmentsToLabel { found: usize, }, @@ -77,11 +75,7 @@ impl ElutionGroupInput { }) } - /// Fill in `?1`, `?2`, ... for an input that named no fragment labels. - /// - /// Minted through [`UnknownIonCounter`] rather than an index cast, so - /// running past the label space is an error rather than a wrap into - /// duplicate labels. + /// Fill in unique unknown-ion labels for an input with no fragment labels. pub fn try_fill_labels_annot( self, ) -> Result, ElutionGroupInputError> { diff --git a/rust/timsquery/src/serde/skyline_io.rs b/rust/timsquery/src/serde/skyline_io.rs index d4630390..78683fc3 100644 --- a/rust/timsquery/src/serde/skyline_io.rs +++ b/rust/timsquery/src/serde/skyline_io.rs @@ -48,8 +48,7 @@ impl std::error::Error for SkylineReadingError {} #[derive(Debug)] pub enum SkylinePrecursorParsingError { - /// Which fragment row failed, and why. One failure aborts the whole library - /// load, and nothing upstream knows where it happened. + /// Ion parsing failed at this fragment row. IonParsing { row: usize, source: IonParsingError, @@ -69,8 +68,6 @@ impl std::fmt::Display for SkylinePrecursorParsingError { } impl SkylinePrecursorParsingError { - /// `map_err` adaptor that stamps the fragment row onto an ion-parsing - /// failure. Replaces a `From` impl, which had no way to see the row. fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { move |source| Self::IonParsing { row, source } } diff --git a/rust/timsquery/src/serde/spectronaut_io.rs b/rust/timsquery/src/serde/spectronaut_io.rs index 0b3ccf58..4b81e9ac 100644 --- a/rust/timsquery/src/serde/spectronaut_io.rs +++ b/rust/timsquery/src/serde/spectronaut_io.rs @@ -47,8 +47,7 @@ impl std::error::Error for SpectronautReadingError {} #[derive(Debug)] pub enum SpectronautPrecursorParsingError { - /// Which fragment row failed, and why. One failure aborts the whole library - /// load, and nothing upstream knows where it happened. + /// Ion parsing failed at this fragment row. IonParsing { row: usize, source: IonParsingError, @@ -70,8 +69,6 @@ impl std::fmt::Display for SpectronautPrecursorParsingError { } impl SpectronautPrecursorParsingError { - /// `map_err` adaptor that stamps the fragment row onto an ion-parsing - /// failure. Replaces a `From` impl, which had no way to see the row. fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { move |source| Self::IonParsing { row, source } } diff --git a/rust/timsseek/src/ml/cv.rs b/rust/timsseek/src/ml/cv.rs index cbd24006..f80a1abd 100644 --- a/rust/timsseek/src/ml/cv.rs +++ b/rust/timsseek/src/ml/cv.rs @@ -354,24 +354,15 @@ impl FoldDataset for StreamingDataset<'_, T> { /// A model [`CrossValidatedScorer`] can cross-fit: fit on a row-index slice, /// score another, and report per-column importance. /// -/// `fit` receives BOTH the training rows and an early-stopping (`val`) slice. -/// A model without early stopping is expected to ignore `val`; only the -/// scorer's guarantee matters, namely that neither slice is ever scored by the -/// model fitted from it. +/// `fit` receives training rows and an optional early-stopping slice. Neither +/// slice is scored by the model fitted from it. pub(crate) trait FoldModel: Sized { type Config; type Error; /// Fit one fold's model. /// - /// `fold` is the identity of the fold being fitted, and it is a PARAMETER - /// rather than config state because the caller is the only thing that knows - /// it: a config is shared across every fold, and deriving the fold from - /// `train[0]` breaks on an empty slice and is meaningless for a partition - /// that trains on several folds at once (`qvalues::crossfit_lda` trains on - /// all folds but one). Models with a stochastic initialization mix it into - /// their seed, so that folds differ from each other while a rerun of the - /// same fold does not. Deterministic models ignore it, the same way a model - /// without early stopping ignores `val`. + /// `fold` identifies the fold being fitted. Models may use it to derive a + /// fold-specific random seed. fn fit( cfg: &Self::Config, data: &D, @@ -419,9 +410,6 @@ pub(crate) trait FoldModel: Sized { /// Dataset over an already-materialized row-major matrix. `get_values` is a /// `copy_from_slice` out of the existing slab. /// -/// `pub(crate)`, matching [`PrecomputedFeatures`]: the only constructor is -/// `pub(crate)` and there is no other public inherent method, so `pub` made the -/// name visible outside the crate while leaving it unconstructable there. pub(crate) struct RowMajorDataset { features: PrecomputedFeatures, names: Vec>, @@ -544,12 +532,8 @@ impl DataBuffer { /// Gather the FEATURES of `rows` (in the given order) out of `data` into /// feature-major layout: `fold_buffer[feature_idx * nrows + sample_idx]`. /// - /// Deliberately does NOT read [`FoldDataset::is_decoy`]. The scoring path has - /// no use for the labels, and this module is where leak-freedom is enforced, - /// so "the labels are not even read here" is worth being structurally true - /// rather than true by inspection of a discarded binding. `response_buffer` - /// is cleared, so [`Self::as_matrix`] cannot hand out a stale label slice - /// from a previous gather; use [`Self::features_as_matrix`] after this. + /// Does not read labels. Clears the response buffer so a features-only gather + /// cannot expose labels from a previous gather. fn fill_features_from(&mut self, data: &D, rows: &[usize], ncols: usize) { self.fold_buffer.clear(); self.response_buffer.clear(); @@ -583,15 +567,7 @@ impl DataBuffer { /// Features + labels. Only valid after [`Self::fill_from`]. /// - /// The assertion catches a features-only gather OF A NON-EMPTY ROW SET -- - /// [`Self::fill_features_from`] clears `response_buffer`, so the labels are - /// missing while `nrows` is not, and this panics instead of handing out an - /// empty response slice. It does NOT catch the zero-row case: with - /// `nrows == 0` both sides are `0` and the two fills are indistinguishable - /// here, so an empty gather returns an empty view either way. That is the - /// shape `crossfit` produces (it passes `val = &[]`), and it is harmless - /// because that path uses [`Self::fill_from`] anyway -- but the assertion is - /// not what makes it safe. + /// Panics unless the label buffer matches the row count. fn as_matrix(&self) -> FoldView<'_> { let mat = self.features_as_matrix(); assert_eq!(self.response_buffer.len(), self.nrows); @@ -610,10 +586,8 @@ pub(crate) fn fold_weights(data: &D, rows: &[usize]) -> Vec .collect() } -/// [`FoldModel`] adapter for `forust`'s [`GradientBooster`]. A newtype only -/// because `GradientBooster` is a foreign type; it also carries the lane width -/// so `FoldModel::importance` can return a full-width, lane-indexed vector -/// (forust reports only the columns it split on). +/// [`FoldModel`] adapter for `forust`'s [`GradientBooster`]. It carries the lane +/// width so importance can be returned in full feature order. pub(crate) struct GbmFoldModel { booster: GradientBooster, ncols: usize, @@ -623,10 +597,7 @@ impl FoldModel for GbmFoldModel { type Config = GBMConfig; type Error = ForustError; - /// `fold` is DELIBERATELY IGNORED: forust seeds itself from - /// `GBMConfig::seed`, so mixing the fold in here would fight the one seed - /// the booster already has. Folds still differ, because they are fitted on - /// different rows. + /// The booster owns its seed; fold identity is not mixed into it. fn fit( cfg: &GBMConfig, data: &D, @@ -720,10 +691,8 @@ impl PrecomputedFeatures { } /// Build from an already-materialized row-major matrix - /// (`features[i*ncols + j]`) + responses, instead of walking - /// Rows MUST align with the `data` the scorer is constructed from (same - /// order). This is how the lane-matrix consumer trains GBM on a prebuilt - /// lane feature set. + /// (`features[i*ncols + j]`) and responses. Rows must align with the + /// scorer's data order. pub(crate) fn from_row_major(features: Vec, ncols: usize, responses: Vec) -> Self { assert_eq!( features.len(), @@ -752,18 +721,8 @@ pub struct FeatureStat { pub nan_ratio: f32, } -/// Per-fold feature statistics. -/// -/// `feature_stats` is in the dataset's own column order -/// (`FoldDataset::column_names`, i.e. the matrix's column order), one entry per -/// column. `feature_importance` is sorted by importance DESCENDING (top features -/// first), and carries only the columns the model reported a finite value for -- -/// see `FoldModel::importance` -- so it is generally shorter than -/// `feature_stats` and in a different order. -/// -/// Descending because the two TSV sidecars `timsseek_cli` writes -/// (`results.feature_stats.tsv` / `results.feature_importance.tsv`) emit these -/// vectors in order, one row per entry, and the importance one is read top-down. +/// Per-fold feature statistics in dataset column order. Feature importance is +/// sorted descending and includes only finite reported values. #[derive(Debug, Serialize)] pub struct FoldStats { pub fold: u8, @@ -773,20 +732,9 @@ pub struct FoldStats { pub type RescoreFeatureStats = Vec; -/// Per-fold feature means/NaN ratios + model importance, for ANY fold -/// partition. -/// -/// Partition-agnostic on purpose: `fold_rows[f]` is simply "the rows summarized -/// under fold `f`" and `models[f]` is "the model whose importance is reported -/// there". The two rescoring partitions in this crate disagree on both -- the -/// [`CrossValidatedScorer`] fits on fold `f` and scores the others, while -/// `qvalues::crossfit_lda` fits on everything BUT fold `f` and scores only fold -/// `f` -- and both are leak-free. Keeping this function ignorant of which one it -/// is handed is what lets the sidecar have one implementation without the two -/// partitions being forced to converge. -/// -/// Column names (and therefore the row width) come from `data`, so the stats -/// align with the matrix the model saw by construction. +/// Per-fold feature means, NaN ratios and model importance. `fold_rows[f]` +/// contains the rows summarized for fold `f`, and `models[f]` supplies that +/// fold's importance values. pub(crate) fn fold_feature_stats( data: &D, fold_rows: &[Vec], @@ -796,16 +744,7 @@ pub(crate) fn fold_feature_stats( let mut out: RescoreFeatureStats = Vec::with_capacity(fold_rows.len()); let mut row_buf = vec![0.0f64; names.len()]; for (fold, rows) in fold_rows.iter().enumerate() { - // --- Importance, back to the (name, gain) sidecar shape --- - // UNREPORTED (`NAN`) columns are dropped; every REPORTED column is - // emitted, `0.0` included. See the `FoldModel::importance` contract - // for why those are different things. The drop matters because the - // sidecar and the dashboard's fold-averaged gain both treat a feature - // as "reported by this fold" simply by being present, so a model's - // unmeasured columns would otherwise pad the averaging divisor with - // values it never produced. - // - // Retain finite zeroes: they are reported measurements. + // Keep finite importances, including zero; NaN means unreported. let importance: Vec<(Arc, f32)> = match models.get(fold).copied().flatten() { Some(model) => { let raw_imp = model.importance(); @@ -1472,17 +1411,7 @@ mod test { } } - /// THE `FoldModel::importance` sentinel contract, at the sidecar boundary: - /// `NAN` means "this model reports nothing for this column" and is dropped; - /// every FINITE value reaches the sidecar, `0.0` included. - /// - /// Regression guard. The boundary used to filter `!= 0.0`, which was right - /// for a tree model (where 0.0 did mean "never split on") and wrong for - /// every other model: an LDA's `|coef|` of exactly 0.0 is a measurement of a - /// dead or constant column, and silently deleting those rows removes - /// exactly what an operator reads the sidecar to find. Both halves are - /// asserted -- dropping the NAN alone would also pass if 0.0 were dropped - /// too, so the surviving zero is the load-bearing assertion. + /// NaN importances are omitted; finite values, including zero, are retained. #[test] fn importance_nan_is_unreported_but_zero_is_a_value() { let dataset = RowMajorDataset::new( @@ -1514,9 +1443,7 @@ mod test { assert_eq!(stats[0].feature_stats.len(), 3); } - /// The GBM side of the same contract: forust reports only the columns it - /// split on, and the ones it did not must come back as `NAN` (absent from - /// the sidecar), NOT as a 0.0 gain it never measured. + /// Forust reports only columns it split on; other columns are NaN. /// /// `feature_4` is constant, so no tree can split on it. The other four are /// the usual separable draws, so their presence keeps the absence diff --git a/rust/timsseek/src/ml/qvalues.rs b/rust/timsseek/src/ml/qvalues.rs index e16a08fe..01d81c54 100644 --- a/rust/timsseek/src/ml/qvalues.rs +++ b/rust/timsseek/src/ml/qvalues.rs @@ -259,48 +259,13 @@ impl CrossFit { /// Cross-fit a [`FoldModel`] over the canonical rescore fold partition and /// return each row's HELD-OUT score. /// -/// The shared statement of leak-freedom for `crossfit_lda`, `rescore_lda`, -/// and `rescore_hybrid`. -/// /// Generic over the model because the partition is independent of the fitted /// model. `what` names the model in failure logs. /// -/// # Why a row may never be scored by a model that saw it -/// These models are label-aware: they fit on the target/decoy labels. One -/// in-sample fit over all rows, scoring those same rows, lets every row's -/// discriminant peek at its own label; the target/decoy separation then looks -/// better than it is, and `assign_qval` derives the q-values from exactly that -/// separation. The result is an FDR that is wrong in the flattering direction -/// with nothing downstream to catch it. -/// -/// The hybrid needs this even more sharply: there the held-out score is not the -/// final score but one `lda_score` column fed to a -/// cross-validated GBM. An in-sample column would smuggle a row's own label -/// into a feature the GBM reads while that row is held out of its GBM fold -- so -/// the GBM's own CV cannot notice, and the leak surfaces only as an optimistic -/// FDR. -/// /// # Partition -/// Fold membership comes from [`FoldDataset::get_fold`] -- for the -/// production datasets, that is `i % N_RESCORE_FOLDS`. For -/// fold `f` the model is fitted on every row with `get_fold(i) != f` and then -/// scores only the rows with `get_fold(i) == f`, so no row contributes to the -/// model scoring it. -/// -/// This is NOT `CrossValidatedScorer`'s partition and must not be unified with -/// it: that one fits on fold `f` alone, early-stops on `f + 1`, and scores the -/// remaining `n_folds - 2` folds, so a row is scored by several models that -/// each saw a fraction of the data. Here a row is scored exactly once, by a -/// model that saw everything else. Both satisfy leak-freedom, which is the only -/// property either has to satisfy; how many rows a model sees and how many -/// models score a row are free to differ, and do. -/// -/// What the two DO have to agree on is the fold ASSIGNMENT `get_fold`, or a -/// hybrid row's cross-fit column can come from a model trained on rows the GBM -/// is holding out -- leak restored, silently. Both sides read that assignment -/// from [`FoldDataset::get_fold`], so it has one definition; the -/// `crossfit_holds_out_exactly_the_rows_the_gbm_scorer_trains_on` test pins the -/// two partitions against each other rather than against a re-typed modulo. +/// For fold `f`, the model is fitted on rows where `get_fold(i) != f` and +/// scores only rows where `get_fold(i) == f`. Fold assignment is shared with +/// the GBM scorer when this output is used as a hybrid feature. /// /// Returns an error if any fold cannot fit or score all of its held-out rows. /// Callers therefore cannot accidentally consume a partially filled score @@ -324,18 +289,10 @@ where ); for f in 0..n_folds { - // TRAIN rows = every row NOT in fold f. HELD = exactly fold f. - // Ascending in both cases, so the row order the fit reduces over is the - // dataset's own order. + // Train on every row outside fold f and score fold f. let train: Vec = (0..nrows).filter(|&i| data.get_fold(i) != f).collect(); let held: Vec = (0..nrows).filter(|&i| data.get_fold(i) == f).collect(); - // `val` is EMPTY and the partition is why: this walk trains on - // all-but-fold-`f` and scores fold `f`, so every row is already spoken - // for and there is no third slice to hand over. The LDA is closed-form - // and ignores it. `MlpFoldModel` DOES early-stop, and handles the empty - // slice by carving a deterministic inner validation set out of `train` - // -- see its `fit` for the rule; nothing about it reaches fold `f`. let model = M::fit(cfg, data, f, &train, &[]).map_err(|e| RescoreError::CrossFit { model: what, fold: f, @@ -343,7 +300,6 @@ where stage: "fit", reason: e.to_string(), })?; - // Score ONLY the held-out fold, with a model that never saw it. let preds = model .predict(data, &held) .map_err(|e| RescoreError::CrossFit { @@ -447,14 +403,9 @@ pub fn rescore(mut data: Vec) -> RescoreResult { /// is generally much cheaper and less sensitive than GBM, with neither gap /// constant across candidate counts. /// -/// The FDR machinery (`assign_qval`, target-decoy competition) is untouched -- -/// only the discriminant score source changes. -/// -/// CROSS-FIT, not a single in-sample fit: every row's score comes from an LDA -/// fitted without that row, via `crossfit_lda` -- see `crossfit` for the -/// partition and why it is mandatory. +/// Scores each row with an LDA model fitted without that row. /// -/// Returns PER-FOLD `FoldStats` (one per fold, like the GBM path): feature +/// Returns per-fold `FoldStats`: feature /// means/NaN ratios over each fold's held-out rows, `|coef|` importance from /// that fold's model. /// @@ -499,10 +450,8 @@ pub fn rescore_lda(mut data: Vec) -> RescoreResult { /// /// `data` must already be through [`canonicalize_and_shuffle`]. /// -/// The fold count is `N_RESCORE_FOLDS`, i.e. the SAME -/// `get_fold` the GBM's [`CrossValidatedScorer`] derives its partition from -/// below. That shared definition is what makes "the same fold assignment on both -/// sides" structural rather than a comment; see `crossfit`. +/// The fold count and assignment match the GBM's +/// [`CrossValidatedScorer`] partition. fn hybrid_linear_dataset(data: &[CompetedCandidate]) -> StreamingDataset<'_, CompetedCandidate> { StreamingDataset::new( data, @@ -555,16 +504,8 @@ fn hybrid_frame( ) } -/// Hybrid rescorer: cross-fit an LDA on the LINEAR lane, push its (leak-free) -/// `lda_score` as one extra column into the NONLINEAR lane, then train the GBM -/// CV on `nonlinear + lda_score` instead of the full feature frame. -/// -/// LEAK-FREEDOM: `lda_score` is cross-fit via `crossfit` -- see there for the -/// partition, why a label-aware feature fed to a CV'd GBM in particular must be -/// leak-free, and why the fold ASSIGNMENT has to match the one -/// `CrossValidatedScorer` derives its own partition from. ASSIGNMENT, not -/// partition: the two train/score splits differ deliberately and both are -/// leak-free, so "unifying" them is the refactor to not make. +/// Hybrid rescorer: cross-fit an LDA on the LINEAR lane, append its held-out +/// score to the NONLINEAR lane, then train the GBM on both lanes. /// /// Selected via the `rescore_model` config field / `--rescore-model` CLI flag /// ([`crate::ml::RescoreModel::Hybrid`]). @@ -774,8 +715,7 @@ fn write_competed_linear_row(candidate: &CompetedCandidate, out: &mut [f64]) { sink.finish(); } -/// The LINEAR-lane matrix for `data` in its CURRENT order (call AFTER any -/// shuffle, so row `i` aligns with `data[i]`). `LINEAR_NCOLS` wide. +/// The LINEAR-lane matrix for `data` in its current order. #[cfg(test)] fn build_linear_matrix(data: &[CompetedCandidate]) -> Vec { let mut out = Vec::with_capacity(data.len() * LINEAR_NCOLS); @@ -799,12 +739,8 @@ fn build_nonlinear_matrix(data: &[CompetedCandidate]) -> Vec { /// The ALL-lane matrix (linear then nonlinear, per row) -- the GBM feature set, /// `ALL_NCOLS` wide, matching [`all_feature_name_set`]'s order. /// -/// ONE pass over `rows` with ONE `Derived::compute` per row: the two lanes are -/// adjacent within a row, so there is nothing to gain from walking twice. -/// -/// Takes `(scoring, meta)` pairs rather than a row type, because both sides of -/// rescoring feed it and they agree on nothing else. See [`competed_rows`] for -/// the pre-rescore side and [`feature_frame`] for the post-rescore one. +/// Each row is projected once, with its linear features followed by its +/// nonlinear features. fn build_all_matrix<'a>( rows: impl ExactSizeIterator, ) -> Vec { From 1021be6cd86cd4b4bd30d6b13a7fe18ccca668c4 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Fri, 28 Aug 2026 12:48:29 -0700 Subject: [PATCH 8/8] style: cargo fmt --- rust/micromzpaf/src/loss.rs | 2 -- rust/timsseek/src/ml/cv.rs | 1 - 2 files changed, 3 deletions(-) diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs index ddd7cdaa..29289c84 100644 --- a/rust/micromzpaf/src/loss.rs +++ b/rust/micromzpaf/src/loss.rs @@ -49,7 +49,6 @@ impl Composition { /// The count for one element symbol, or `None` if this crate does not /// represent that element. - /// fn count_mut(&mut self, symbol: u8) -> Option<&mut u8> { Some(match symbol { b'C' => &mut self.c, @@ -228,7 +227,6 @@ impl NeutralLoss { /// The [`TABLE`] row a discriminant names, or `None` for [`Self::None`] and /// for the reserved values above the table. - /// fn at(d: u8) -> Option<&'static (Composition, NeutralLoss, &'static str)> { TABLE.get((d as usize).checked_sub(1)?) } diff --git a/rust/timsseek/src/ml/cv.rs b/rust/timsseek/src/ml/cv.rs index f80a1abd..480d6ef4 100644 --- a/rust/timsseek/src/ml/cv.rs +++ b/rust/timsseek/src/ml/cv.rs @@ -409,7 +409,6 @@ pub(crate) trait FoldModel: Sized { /// Dataset over an already-materialized row-major matrix. `get_values` is a /// `copy_from_slice` out of the existing slab. -/// pub(crate) struct RowMajorDataset { features: PrecomputedFeatures, names: Vec>,