mscore/chemistry/
sum_formula.rs

1use std::collections::HashMap;
2use crate::algorithm::isotope::generate_isotope_distribution;
3use crate::chemistry::constants::MASS_PROTON;
4use crate::chemistry::elements::atomic_weights_mono_isotopic;
5use crate::data::spectrum::MzSpectrum;
6
7pub struct SumFormula {
8    pub formula: String,
9    pub elements: HashMap<String, i32>,
10}
11
12impl SumFormula {
13    pub fn new(formula: &str) -> Self {
14        let elements = parse_formula(formula).unwrap();
15        SumFormula {
16            formula: formula.to_string(),
17            elements,
18        }
19    }
20    pub fn monoisotopic_weight(&self) -> f64 {
21        let atomic_weights = atomic_weights_mono_isotopic();
22
23        // Iterate elements in a deterministic (sorted) order. `elements` is a
24        // HashMap whose iteration order is randomized per instance, and
25        // floating-point addition is not associative, so an unsorted fold can
26        // return masses that differ in their least-significant bits between
27        // calls. Sorting pins one canonical accumulation order, matching what
28        // `generate_isotope_distribution` already does for the same map.
29        let mut elements: Vec<(&String, i32)> =
30            self.elements.iter().map(|(k, v)| (k, *v)).collect();
31        elements.sort_unstable_by(|a, b| a.0.cmp(b.0));
32
33        elements.into_iter().fold(0.0, |acc, (element, count)| {
34            acc + atomic_weights[element.as_str()] * count as f64
35        })
36    }
37
38    pub fn isotope_distribution(&self, charge: i32) -> MzSpectrum {
39        let distribution = generate_isotope_distribution(&self.elements, 1e-3, 1e-9, 200);
40        let intensity = distribution.iter().map(|(_, i)| *i).collect();
41        let mz = distribution.iter().map(|(m, _)| (*m + charge as f64 * MASS_PROTON) / charge as f64).collect();
42        MzSpectrum::new(mz, intensity)
43    }
44}
45
46fn parse_formula(formula: &str) -> Result<HashMap<String, i32>, String> {
47    let atomic_weights = atomic_weights_mono_isotopic();
48    let mut element_counts = HashMap::new();
49    let mut current_element = String::new();
50    let mut current_count = String::new();
51    let mut chars = formula.chars().peekable();
52
53    while let Some(c) = chars.next() {
54        if c.is_ascii_uppercase() {
55            if !current_element.is_empty() {
56                let count = current_count.parse::<i32>().unwrap_or(1);
57                if atomic_weights.contains_key(current_element.as_str()) {
58                    *element_counts.entry(current_element.clone()).or_insert(0) += count;
59                } else {
60                    return Err(format!("Unknown element: {}", current_element));
61                }
62            }
63            current_element = c.to_string();
64            current_count = String::new();
65        } else if c.is_ascii_digit() {
66            current_count.push(c);
67        } else if c.is_ascii_lowercase() {
68            current_element.push(c);
69        }
70
71        if chars.peek().map_or(true, |next_c| next_c.is_ascii_uppercase()) {
72            let count = current_count.parse::<i32>().unwrap_or(1);
73            if atomic_weights.contains_key(current_element.as_str()) {
74                *element_counts.entry(current_element.clone()).or_insert(0) += count;
75            } else {
76                return Err(format!("Unknown element: {}", current_element));
77            }
78            current_element = String::new();
79            current_count = String::new();
80        }
81    }
82
83    Ok(element_counts)
84}
85
86#[cfg(test)]
87mod tests {
88    use super::*;
89    use std::collections::BTreeSet;
90
91    // A formula whose element masses are genuinely order-sensitive: summing its
92    // five elements in different orders yields two distinct f64 bit patterns
93    // (a 66/54 split over the 120 permutations), so a randomized HashMap order
94    // reliably produces both within this many draws. Small formulas such as
95    // C6H12O6 round identically under every permutation and cannot catch the bug.
96    const ORDER_SENSITIVE_FORMULA: &str = "C100H150N20O30S2";
97
98    #[test]
99    fn monoisotopic_weight_is_bitwise_stable() {
100        let mut observed_bits = BTreeSet::new();
101        for _ in 0..1024 {
102            let formula = SumFormula::new(ORDER_SENSITIVE_FORMULA);
103            observed_bits.insert(formula.monoisotopic_weight().to_bits());
104        }
105
106        assert_eq!(
107            observed_bits.len(),
108            1,
109            "identical sum formula produced multiple exact monoisotopic weights: {observed_bits:?}"
110        );
111    }
112
113    #[test]
114    fn monoisotopic_weight_matches_expected_mass() {
115        let formula = SumFormula::new("C6H12O6");
116        let quantized = (formula.monoisotopic_weight() * 1e6).round() as i64;
117        assert_eq!(quantized, 180063388);
118    }
119}
mscore/chemistry/sum_formula.rs

mscore/chemistry/
sum_formula.rs
