mscore/data/
spectrum.rs

1use std::fmt;
2use std::collections::BTreeMap;
3use std::sync::Arc;
4use nalgebra::DVector;
5use std::fmt::{Display, Formatter};
6use bincode::{Decode, Encode};
7use serde::{Serialize, Deserialize};
8
9extern crate rand;
10
11use rand::distributions::{Uniform, Distribution};
12use rand::Rng;
13use statrs::distribution::Normal;
14
15/// Represents a vectorized mass spectrum.
16pub trait ToResolution {
17    fn to_resolution(&self, resolution: i32) -> Self;
18}
19
20/// Vectorized representation for Structs holding m/z values and intensities.
21pub trait Vectorized<T> {
22    fn vectorized(&self, resolution: i32) -> T;
23}
24
25/// Represents the type of spectrum.
26///
27/// # Description
28///
29/// The `SpecType` enum is used to distinguish between precursor and fragment spectra.
30///
31#[derive(Clone, PartialEq, Debug, Serialize, Deserialize, Encode, Decode)]
32pub enum MsType {
33    Precursor,
34    FragmentDda,
35    FragmentDia,
36    Unknown,
37}
38
39impl MsType {
40    /// Returns the `MsType` enum corresponding to the given integer value.
41    ///
42    /// # Arguments
43    ///
44    /// * `ms_type` - An integer value corresponding to the `MsType` enum.
45    ///
46    pub fn new(ms_type: i32) -> MsType {
47        match ms_type {
48            0 => MsType::Precursor,
49            8 => MsType::FragmentDda,
50            9 => MsType::FragmentDia,
51            _ => MsType::Unknown,
52        }
53    }
54
55    /// Returns the integer value corresponding to the `MsType` enum.
56    pub fn ms_type_numeric(&self) -> i32 {
57        match self {
58            MsType::Precursor => 0,
59            MsType::FragmentDda => 8,
60            MsType::FragmentDia => 9,
61            MsType::Unknown => -1,
62        }
63    }
64}
65
66impl Default for MsType {
67    fn default() -> Self {
68        MsType::Unknown
69    }
70}
71
72impl Display for MsType {
73    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
74        match self {
75            MsType::Precursor => write!(f, "Precursor"),
76            MsType::FragmentDda => write!(f, "FragmentDda"),
77            MsType::FragmentDia => write!(f, "FragmentDia"),
78            MsType::Unknown => write!(f, "Unknown"),
79        }
80    }
81}
82
83/// Represents a mass spectrum with associated m/z values and intensities.
84///
85/// Uses Arc<Vec<T>> for efficient cloning - clone is O(1) instead of O(n).
86#[derive(Clone, Debug, Serialize, Deserialize)]
87pub struct MzSpectrum {
88    pub mz: Arc<Vec<f64>>,
89    pub intensity: Arc<Vec<f64>>,
90}
91
92// Manual bincode implementation for Arc compatibility
93impl Encode for MzSpectrum {
94    fn encode<E: bincode::enc::Encoder>(&self, encoder: &mut E) -> Result<(), bincode::error::EncodeError> {
95        // Encode the inner Vec directly
96        bincode::Encode::encode(&*self.mz, encoder)?;
97        bincode::Encode::encode(&*self.intensity, encoder)?;
98        Ok(())
99    }
100}
101
102impl<Context> Decode<Context> for MzSpectrum {
103    fn decode<D: bincode::de::Decoder<Context = Context>>(decoder: &mut D) -> Result<Self, bincode::error::DecodeError> {
104        let mz: Vec<f64> = bincode::Decode::decode(decoder)?;
105        let intensity: Vec<f64> = bincode::Decode::decode(decoder)?;
106        Ok(MzSpectrum::new(mz, intensity))
107    }
108}
109
110impl<'de, Context> bincode::BorrowDecode<'de, Context> for MzSpectrum {
111    fn borrow_decode<D: bincode::de::BorrowDecoder<'de, Context = Context>>(decoder: &mut D) -> Result<Self, bincode::error::DecodeError> {
112        let mz: Vec<f64> = bincode::BorrowDecode::borrow_decode(decoder)?;
113        let intensity: Vec<f64> = bincode::BorrowDecode::borrow_decode(decoder)?;
114        Ok(MzSpectrum::new(mz, intensity))
115    }
116}
117
118impl MzSpectrum {
119    /// Constructs a new `MzSpectrum`.
120    ///
121    /// # Arguments
122    ///
123    /// * `mz` - A vector of m/z values.
124    /// * `intensity` - A vector of intensity values corresponding to the m/z values.
125    ///
126    /// # Panics
127    ///
128    /// Panics if the lengths of `mz` and `intensity` are not the same. (actually, it doesn't at the moment, planning on adding this later)
129    ///
130    /// # Example
131    ///
132    /// ```rust
133    /// # use mscore::data::spectrum::MzSpectrum;
134    /// let spectrum = MzSpectrum::new(vec![100.0, 200.0], vec![10.0, 20.0]);
135    /// assert_eq!(*spectrum.mz, vec![100.0, 200.0]);
136    /// assert_eq!(*spectrum.intensity, vec![10.0, 20.0]);
137    /// ```
138    pub fn new(mz: Vec<f64>, intensity: Vec<f64>) -> Self {
139        MzSpectrum {
140            mz: Arc::new(mz),
141            intensity: Arc::new(intensity),
142        }
143    }
144
145    pub fn filter_ranged(&self, mz_min: f64, mz_max: f64, intensity_min:f64, intensity_max: f64) -> Self {
146        let mut mz_vec: Vec<f64> = Vec::new();
147        let mut intensity_vec: Vec<f64> = Vec::new();
148
149        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
150            if mz_min <= *mz && *mz <= mz_max && *intensity >= intensity_min && *intensity <= intensity_max {
151                mz_vec.push(*mz);
152                intensity_vec.push(*intensity);
153            }
154        }
155        MzSpectrum::new(mz_vec, intensity_vec)
156    }
157
158    /// Splits the spectrum into a collection of windows based on m/z values.
159    ///
160    /// This function divides the spectrum into smaller spectra (windows) based on a specified window length.
161    /// Each window contains peaks from the original spectrum that fall within the m/z range of that window.
162    ///
163    /// # Arguments
164    ///
165    /// * `window_length`: The size (in terms of m/z values) of each window.
166    ///
167    /// * `overlapping`: If `true`, each window will overlap with its neighboring windows by half of the `window_length`.
168    ///   This means that a peak may belong to multiple windows. If `false`, windows do not overlap.
169    ///
170    /// * `min_peaks`: The minimum number of peaks a window must have to be retained in the result.
171    ///
172    /// * `min_intensity`: The minimum intensity value a window must have (in its highest intensity peak) to be retained in the result.
173    ///
174    /// # Returns
175    ///
176    /// A `BTreeMap` where the keys represent the window indices and the values are the spectra (`MzSpectrum`) within those windows.
177    /// Windows that do not meet the criteria of having at least `min_peaks` peaks or a highest intensity peak
178    /// greater than or equal to `min_intensity` are discarded.
179    ///
180    /// # Example
181    ///
182    /// ```rust
183    /// # use mscore::data::spectrum::MzSpectrum;
184    /// let spectrum = MzSpectrum::new(vec![100.0, 101.0, 102.5, 103.0], vec![10.0, 20.0, 30.0, 40.0]);
185    /// let windowed_spectrum = spectrum.to_windows(1.0, false, 1, 10.0);
186    /// assert!(windowed_spectrum.contains_key(&100));
187    /// assert!(windowed_spectrum.contains_key(&102));
188    /// ```
189    pub fn to_windows(&self, window_length: f64, overlapping: bool, min_peaks: usize, min_intensity: f64) -> BTreeMap<i32, MzSpectrum> {
190        // Build up vectors first, then wrap in MzSpectrum at the end
191        let mut splits: BTreeMap<i32, (Vec<f64>, Vec<f64>)> = BTreeMap::new();
192
193        for (i, &mz) in self.mz.iter().enumerate() {
194            let intensity = self.intensity[i];
195            let tmp_key = (mz / window_length).floor() as i32;
196
197            let entry = splits.entry(tmp_key).or_insert_with(|| (Vec::new(), Vec::new()));
198            entry.0.push(mz);
199            entry.1.push(intensity);
200        }
201
202        if overlapping {
203            let mut splits_offset: BTreeMap<i32, (Vec<f64>, Vec<f64>)> = BTreeMap::new();
204
205            for (i, &mmz) in self.mz.iter().enumerate() {
206                let intensity = self.intensity[i];
207                let tmp_key = -((mmz + window_length / 2.0) / window_length).floor() as i32;
208
209                let entry = splits_offset.entry(tmp_key).or_insert_with(|| (Vec::new(), Vec::new()));
210                entry.0.push(mmz);
211                entry.1.push(intensity);
212            }
213
214            for (key, (mz_vec, int_vec)) in splits_offset {
215                let entry = splits.entry(key).or_insert_with(|| (Vec::new(), Vec::new()));
216                entry.0.extend(mz_vec);
217                entry.1.extend(int_vec);
218            }
219        }
220
221        // Filter and convert to MzSpectrum
222        splits.into_iter()
223            .filter(|(_, (mz_vec, int_vec))| {
224                mz_vec.len() >= min_peaks && int_vec.iter().cloned().max_by(
225                    |a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)).unwrap_or(0.0) >= min_intensity
226            })
227            .map(|(key, (mz_vec, int_vec))| (key, MzSpectrum::new(mz_vec, int_vec)))
228            .collect()
229    }
230
231    pub fn to_centroid(&self, baseline_noise_level: i32, sigma: f64, normalize: bool) -> MzSpectrum {
232
233        let filtered = self.filter_ranged(0.0, 1e9, baseline_noise_level as f64, 1e9);
234
235        let mut cent_mz = Vec::new();
236        let mut cent_i: Vec<f64> = Vec::new();
237
238        let mut last_mz = 0.0;
239        let mut mean_mz = 0.0;
240        let mut sum_i = 0.0;
241
242        for (i, &current_mz) in filtered.mz.iter().enumerate() {
243            let current_intensity = filtered.intensity[i];
244
245            // If peak is too far away from last peak, push centroid
246            if current_mz - last_mz > sigma && mean_mz > 0.0 {
247                mean_mz /= sum_i;
248                cent_mz.push(mean_mz);
249                cent_i.push(sum_i);
250
251                // Start new centroid
252                sum_i = 0.0;
253                mean_mz = 0.0;
254            }
255
256            mean_mz += current_mz * current_intensity as f64;
257            sum_i += current_intensity;
258            last_mz = current_mz;
259        }
260
261        // Push back last remaining centroid
262        if mean_mz > 0.0 {
263            mean_mz /= sum_i;
264            cent_mz.push(mean_mz);
265            cent_i.push(sum_i);
266        }
267
268        if normalize {
269            let sum_i: f64 = cent_i.iter().sum();
270            cent_i = cent_i.iter().map(|&i| i / sum_i).collect();
271        }
272        MzSpectrum::new(cent_mz, cent_i)
273    }
274
275    pub fn from_collection(collection: Vec<MzSpectrum>) -> MzSpectrum {
276
277        let quantize = |mz: f64| -> i64 {
278            (mz * 1_000_000.0).round() as i64
279        };
280
281        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
282
283        for spectrum in collection {
284            for (mz, intensity) in spectrum.mz.iter().zip(spectrum.intensity.iter()) {
285                let key = quantize(*mz);
286                let entry = combined_map.entry(key).or_insert(0.0);
287                *entry += *intensity;
288            }
289        }
290
291        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
292        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
293
294        MzSpectrum::new(mz_combined, intensity_combined)
295    }
296
297    pub fn add_mz_noise_uniform(&self, ppm: f64, right_drag: bool) -> Self {
298        let mut rng = rand::thread_rng();
299        self.add_mz_noise_uniform_with_rng(ppm, right_drag, &mut rng)
300    }
301
302    /// Same as `add_mz_noise_uniform`, but drawing from a caller-supplied RNG so the jitter is
303    /// reproducible. `thread_rng` makes the written spectrum depend on which thread happened to
304    /// run it; a shifted m/z can cross a TOF bin boundary and change the writer's dedup.
305    pub fn add_mz_noise_uniform_with_rng<R: Rng>(&self, ppm: f64, right_drag: bool, rng: &mut R) -> Self {
306        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
307
308            let ppm_mz = match right_drag {
309                true => mz * ppm / 1e6 / 2.0,
310                false => mz * ppm / 1e6,
311            };
312
313            let dist = match right_drag {
314                true => Uniform::from(mz - (ppm_mz / 3.0)..=mz + ppm_mz),
315                false => Uniform::from(mz - ppm_mz..=mz + ppm_mz),
316            };
317
318            dist.sample(rng)
319        })
320    }
321
322    pub fn add_mz_noise_normal(&self, ppm: f64) -> Self {
323        let mut rng = rand::thread_rng();
324        self.add_mz_noise_normal_with_rng(ppm, &mut rng)
325    }
326
327    /// Reproducible counterpart of `add_mz_noise_normal`; see `add_mz_noise_uniform_with_rng`.
328    pub fn add_mz_noise_normal_with_rng<R: Rng>(&self, ppm: f64, rng: &mut R) -> Self {
329        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
330            let ppm_mz = mz * ppm / 1e6;
331            let dist = Normal::new(mz, ppm_mz / 3.0).unwrap();
332            dist.sample(rng)
333        })
334    }
335
336    fn add_mz_noise<R: Rng, F>(&self, ppm: f64, rng: &mut R, noise_fn: F) -> Self
337        where
338            F: Fn(&mut R, f64, f64) -> f64,
339    {
340        let mz: Vec<f64> = self.mz.iter().map(|&mz_value| noise_fn(rng, mz_value, ppm)).collect();
341        // Clone intensity Arc (O(1)) and wrap new mz in Arc
342        let spectrum = MzSpectrum { mz: Arc::new(mz), intensity: self.intensity.clone() };
343        // Sort the spectrum by m/z values and potentially sum up intensities at the same m/z value
344        spectrum.to_resolution(6)
345    }
346}
347
348impl ToResolution for MzSpectrum {
349    /// Bins the spectrum's m/z values to a given resolution and sums the intensities.
350    ///
351    /// # Arguments
352    ///
353    /// * `resolution` - The desired resolution in terms of decimal places. For instance, a resolution of 2
354    /// would bin m/z values to two decimal places.
355    ///
356    /// # Returns
357    ///
358    /// A new `MzSpectrum` where m/z values are binned according to the given resolution.
359    ///
360    /// # Example
361    ///
362    /// ```rust
363    /// # use mscore::data::spectrum::MzSpectrum;
364    /// # use mscore::data::spectrum::ToResolution;
365    /// let spectrum = MzSpectrum::new(vec![100.123, 100.121, 100.131], vec![10.0, 20.0, 30.0]);
366    /// let binned_spectrum_1 = spectrum.to_resolution(1);
367    /// let binned_spectrum_2 = spectrum.to_resolution(2);
368    /// // assert_eq!(*binned_spectrum_1.mz, vec![100.1]);
369    /// assert_eq!(*binned_spectrum_1.intensity, vec![60.0]);
370    /// assert_eq!(*binned_spectrum_2.mz, vec![100.12, 100.13]);
371    /// assert_eq!(*binned_spectrum_2.intensity, vec![30.0, 30.0]);
372    /// ```
373    fn to_resolution(&self, resolution: i32) -> Self {
374        let mut binned: BTreeMap<i64, f64> = BTreeMap::new();
375        let factor = 10f64.powi(resolution);
376
377        for (mz, inten) in self.mz.iter().zip(self.intensity.iter()) {
378
379            let key = (mz * factor).round() as i64;
380            let entry = binned.entry(key).or_insert(0.0);
381            *entry += *inten;
382        }
383
384        let mz: Vec<f64> = binned.keys().map(|&key| key as f64 / 10f64.powi(resolution)).collect();
385        let intensity: Vec<f64> = binned.values().cloned().collect();
386
387        MzSpectrum::new(mz, intensity)
388    }
389}
390
391impl Vectorized<MzSpectrumVectorized> for MzSpectrum {
392    /// Convert the `MzSpectrum` to a `MzSpectrumVectorized` using the given resolution for binning.
393    ///
394    /// After binning to the desired resolution, the binned m/z values are translated into integer indices.
395    fn vectorized(&self, resolution: i32) -> MzSpectrumVectorized {
396
397        let binned_spectrum = self.to_resolution(resolution);
398
399        // Translate the m/z values into integer indices
400        let indices: Vec<i32> = binned_spectrum.mz.iter().map(|&mz| (mz * 10f64.powi(resolution)).round() as i32).collect();
401
402        MzSpectrumVectorized {
403            resolution,
404            indices,
405            values: (*binned_spectrum.intensity).clone(),
406        }
407    }
408}
409
410/// Formats the `MzSpectrum` for display.
411impl Display for MzSpectrum {
412    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
413
414        let (mz, i) = self.mz.iter()
415            .zip(self.intensity.iter())
416            .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
417            .unwrap();
418
419        write!(f, "MzSpectrum(data points: {}, max  by intensity:({}, {}))", self.mz.len(), format!("{:.3}", mz), i)
420    }
421}
422
423impl std::ops::Add for MzSpectrum {
424    type Output = Self;
425    /// Combines two `MzSpectrum` instances by summing up the intensities of matching m/z values.
426    ///
427    /// # Description
428    /// Each m/z value is quantized to retain at least 6 decimals. If two spectra have m/z values
429    /// that quantize to the same integer value, their intensities are summed.
430    ///
431    /// # Example
432    /// ```
433    /// # use mscore::data::spectrum::MzSpectrum;
434    /// let spectrum1 = MzSpectrum::new(vec![100.523, 101.923], vec![10.0, 20.0]);
435    /// let spectrum2 = MzSpectrum::new(vec![101.235, 105.112], vec![15.0, 30.0]);
436    ///
437    /// let combined = spectrum1 + spectrum2;
438    ///
439    /// assert_eq!(*combined.mz, vec![100.523, 101.235, 101.923, 105.112]);
440    /// assert_eq!(*combined.intensity, vec![10.0, 15.0, 20.0, 30.0]);
441    /// ```
442    fn add(self, other: Self) -> MzSpectrum {
443        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
444
445        // Helper to quantize mz to an integer key
446        let quantize = |mz: f64| -> i64 {
447            (mz * 1_000_000.0).round() as i64
448        };
449
450        // Add the m/z and intensities from the first spectrum to the map
451        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
452            let key = quantize(*mz);
453            combined_map.insert(key, *intensity);
454        }
455
456        // Combine the second spectrum into the map
457        for (mz, intensity) in other.mz.iter().zip(other.intensity.iter()) {
458            let key = quantize(*mz);
459            let entry = combined_map.entry(key).or_insert(0.0);
460            *entry += *intensity;
461        }
462
463        // Convert the combined map back into two Vec<f64>
464        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
465        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
466
467        MzSpectrum::new(mz_combined, intensity_combined)
468    }
469}
470
471impl std::ops::Mul<f64> for MzSpectrum {
472    type Output = Self;
473    fn mul(self, scale: f64) -> Self::Output {
474        let scaled_intensities: Vec<f64> = self.intensity.iter().map(|i| scale * i).collect();
475        // Clone mz Arc (O(1)), wrap new intensities in Arc
476        Self { mz: self.mz.clone(), intensity: Arc::new(scaled_intensities) }
477    }
478}
479
480impl std::ops::Sub for MzSpectrum {
481    type Output = Self;
482    fn sub(self, other: Self) -> Self::Output {
483        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
484
485        // Helper to quantize mz to an integer key
486        let quantize = |mz: f64| -> i64 {
487            (mz * 1_000_000.0).round() as i64
488        };
489
490        // Add the m/z and intensities from the first spectrum to the map
491        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
492            let key = quantize(*mz);
493            combined_map.insert(key, *intensity);
494        }
495
496        // Combine the second spectrum into the map
497        for (mz, intensity) in other.mz.iter().zip(other.intensity.iter()) {
498            let key = quantize(*mz);
499            let entry = combined_map.entry(key).or_insert(0.0);
500            *entry -= *intensity;
501        }
502
503        // Convert the combined map back into two Vec<f64>
504        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
505        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
506
507        MzSpectrum::new(mz_combined, intensity_combined)
508    }
509}
510
511/// Represents a mass spectrum with associated m/z indices, m/z values, and intensities
512#[derive(Clone, Debug)]
513pub struct IndexedMzSpectrum {
514    pub index: Vec<i32>,
515    pub mz_spectrum: MzSpectrum,
516}
517
518// implement default (empty IndexedMzSpectrum) constructor
519impl Default for IndexedMzSpectrum {
520    fn default() -> Self {
521        IndexedMzSpectrum { index: Vec::new(), mz_spectrum: MzSpectrum::new(Vec::new(), Vec::new()) }
522    }
523}
524
525impl IndexedMzSpectrum {
526    /// Creates a new `TOFMzSpectrum` instance.
527    ///
528    /// # Arguments
529    ///
530    /// * `index` - A vector containing the mz index, e.g., time-of-flight values.
531    /// * `mz` - A vector containing the m/z values.
532    /// * `intensity` - A vector containing the intensity values.
533    ///
534    /// # Examples
535    ///
536    /// ```
537    /// use mscore::data::spectrum::IndexedMzSpectrum;
538    /// use mscore::data::spectrum::MzSpectrum;
539    ///
540    /// let spectrum = IndexedMzSpectrum::new(vec![1000, 2000], vec![100.5, 200.5], vec![50.0, 60.0]);
541    /// ```
542    pub fn new(index: Vec<i32>, mz: Vec<f64>, intensity: Vec<f64>) -> Self {
543        IndexedMzSpectrum { index, mz_spectrum: MzSpectrum::new(mz, intensity) }
544    }
545
546    /// Create an IndexedMzSpectrum from an existing MzSpectrum, sharing the Arc-wrapped vectors.
547    ///
548    /// This is more efficient than using `new()` when you already have an MzSpectrum,
549    /// as it avoids copying the mz and intensity vectors.
550    ///
551    /// # Arguments
552    ///
553    /// * `index` - A vector containing the index values.
554    /// * `mz_spectrum` - An existing MzSpectrum to share.
555    pub fn from_mz_spectrum(index: Vec<i32>, mz_spectrum: MzSpectrum) -> Self {
556        IndexedMzSpectrum { index, mz_spectrum }
557    }
558
559    /// Bins the spectrum based on a given m/z resolution, summing intensities and averaging index values
560    /// for m/z values that fall into the same bin.
561    ///
562    /// # Arguments
563    ///
564    /// * `resolution` - The desired m/z resolution for binning.
565    ///
566    /// # Examples
567    ///
568    /// ```
569    /// use mscore::data::spectrum::IndexedMzSpectrum;
570    ///
571    /// let spectrum = IndexedMzSpectrum::new(vec![1000, 2000], vec![100.42, 100.43], vec![50.0, 60.0]);
572    /// let binned_spectrum = spectrum.to_resolution(1);
573    ///
574    /// assert_eq!(*binned_spectrum.mz_spectrum.mz, vec![100.4]);
575    /// assert_eq!(*binned_spectrum.mz_spectrum.intensity, vec![110.0]);
576    /// assert_eq!(binned_spectrum.index, vec![1500]);
577    /// ```
578    pub fn to_resolution(&self, resolution: i32) -> IndexedMzSpectrum {
579
580        let mut mz_bins: BTreeMap<i64, (f64, Vec<i64>)> = BTreeMap::new();
581        let factor = 10f64.powi(resolution);
582
583        for ((mz, intensity), tof_val) in self.mz_spectrum.mz.iter().zip(self.mz_spectrum.intensity.iter()).zip(&self.index) {
584            let key = (mz * factor).round() as i64;
585            let entry = mz_bins.entry(key).or_insert((0.0, Vec::new()));
586            entry.0 += *intensity;
587            entry.1.push(*tof_val as i64);
588        }
589
590        let mz: Vec<f64> = mz_bins.keys().map(|&key| key as f64 / factor).collect();
591        let intensity: Vec<f64> = mz_bins.values().map(|(intensity, _)| *intensity).collect();
592        let tof: Vec<i32> = mz_bins.values().map(|(_, tof_vals)| {
593            let sum: i64 = tof_vals.iter().sum();
594            let count: i32 = tof_vals.len() as i32;
595            (sum as f64 / count as f64).round() as i32
596        }).collect();
597
598        IndexedMzSpectrum {index: tof, mz_spectrum: MzSpectrum::new(mz, intensity) }
599    }
600
601    /// Convert the `IndexedMzSpectrum` to a `IndexedMzVector` using the given resolution for binning.
602    ///
603    /// After binning to the desired resolution, the binned m/z values are translated into integer indices.
604    ///
605    /// # Arguments
606    ///
607    /// * `resolution` - The desired m/z resolution for binning.
608    ///
609    /// # Examples
610    ///
611    /// ```
612    /// use mscore::data::spectrum::IndexedMzSpectrum;
613    ///
614    /// let spectrum = IndexedMzSpectrum::new(vec![1000, 2000], vec![100.42, 100.43], vec![50.0, 60.0]);
615    /// let binned_spectrum = spectrum.to_resolution(1);
616    ///
617    /// assert_eq!(*binned_spectrum.mz_spectrum.mz, vec![100.4]);
618    /// assert_eq!(*binned_spectrum.mz_spectrum.intensity, vec![110.0]);
619    /// assert_eq!(binned_spectrum.index, vec![1500]);
620    /// ```
621    pub fn vectorized(&self, resolution: i32) -> IndexedMzSpectrumVectorized {
622
623        let binned_spectrum = self.to_resolution(resolution);
624
625        // Translate the m/z values into integer indices
626        let indices: Vec<i32> = binned_spectrum.mz_spectrum.mz.iter()
627            .map(|&mz| (mz * 10f64.powi(resolution)).round() as i32).collect();
628
629        IndexedMzSpectrumVectorized {
630            index: binned_spectrum.index,
631            mz_vector: MzSpectrumVectorized {
632                resolution,
633                indices,
634                values: (*binned_spectrum.mz_spectrum.intensity).clone(),
635            }
636        }
637    }
638
639    pub fn filter_ranged(&self, mz_min: f64, mz_max: f64, intensity_min:f64, intensity_max: f64) -> Self {
640        let mut mz_vec: Vec<f64> = Vec::new();
641        let mut intensity_vec: Vec<f64> = Vec::new();
642        let mut index_vec: Vec<i32> = Vec::new();
643
644        for ((&mz, &intensity), &index) in self.mz_spectrum.mz.iter().zip(self.mz_spectrum.intensity.iter()).zip(self.index.iter()) {
645            if mz_min <= mz && mz <= mz_max && intensity >= intensity_min && intensity <= intensity_max {
646                mz_vec.push(mz);
647                intensity_vec.push(intensity);
648                index_vec.push(index);
649            }
650        }
651        IndexedMzSpectrum { index: index_vec, mz_spectrum: MzSpectrum::new(mz_vec, intensity_vec) }
652    }
653}
654
655impl Display for IndexedMzSpectrum {
656    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
657        let (mz, i) = self.mz_spectrum.mz.iter()
658            .zip(self.mz_spectrum.intensity.iter())
659            .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
660            .unwrap();
661
662        write!(f, "IndexedMzSpectrum(data points: {}, max  by intensity:({}, {}))", self.mz_spectrum.mz.len(), format!("{:.3}", mz), i)
663    }
664}
665
666#[derive(Clone)]
667pub struct MzSpectrumVectorized {
668    pub resolution: i32,
669    pub indices: Vec<i32>,
670    pub values: Vec<f64>,
671}
672
673impl MzSpectrumVectorized {
674    /// Convert the `MzVector` to a dense vector with a specified maximum index.
675    ///
676    /// The resulting vector has length equal to `max_index + 1` and its values
677    /// are the intensities corresponding to each index. Indices with no associated intensity will have a value of 0.
678    ///
679    /// # Arguments
680    ///
681    /// * `max_index` - The maximum index for the dense vector.
682    
683    fn get_max_index(&self) -> usize {
684        let base: i32 = 10;
685        let max_mz: i32 = 2000;
686        let max_index: usize = (max_mz*base.pow(self.resolution as u32)) as usize;
687        max_index
688    }
689
690    pub fn to_dense(&self, max_index: Option<usize>) -> DVector<f64> {
691        let max_index = match max_index {
692            Some(max_index) => max_index,
693            None => self.get_max_index(),
694        };
695        let mut dense_intensities: DVector<f64> = DVector::<f64>::zeros(max_index + 1);
696        for (&index, &intensity) in self.indices.iter().zip(self.values.iter()) {
697            if (index as usize) <= max_index {
698                dense_intensities[index as usize] = intensity;
699            }
700        }
701        dense_intensities
702    }
703    pub fn to_dense_spectrum(&self, max_index: Option<usize>) -> MzSpectrumVectorized{
704        let max_index = match max_index {
705            Some(max_index) => max_index,
706            None => self.get_max_index(),
707        };
708        let dense_intensities: Vec<f64> = self.to_dense(Some(max_index)).data.into();
709        let dense_indices: Vec<i32> = (0..=max_index).map(|i| i as i32).collect();
710        let dense_spectrum: MzSpectrumVectorized = MzSpectrumVectorized { resolution: (self.resolution), indices: (dense_indices), values: (dense_intensities) };
711        dense_spectrum
712    }
713}
714
715#[derive(Clone)]
716pub struct IndexedMzSpectrumVectorized {
717    pub index: Vec<i32>,
718    pub mz_vector: MzSpectrumVectorized,
719}
720#[cfg(test)]
721mod mz_noise_seeded_tests {
722    use super::*;
723    use crate::simulation::noise_rng::noise_rng;
724
725    fn spectrum() -> MzSpectrum {
726        MzSpectrum::new(
727            vec![300.1234, 450.5678, 700.9012, 1100.3456],
728            vec![1000.0, 2500.0, 700.0, 150.0],
729        )
730    }
731
732    #[test]
733    fn seeded_uniform_jitter_is_reproducible() {
734        let s = spectrum();
735        let a = s.add_mz_noise_uniform_with_rng(6.5, false, &mut noise_rng(41, &[7]));
736        let b = s.add_mz_noise_uniform_with_rng(6.5, false, &mut noise_rng(41, &[7]));
737        assert_eq!(a.mz, b.mz);
738        assert_eq!(a.intensity, b.intensity);
739    }
740
741    #[test]
742    fn seeded_normal_jitter_is_reproducible() {
743        let s = spectrum();
744        let a = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
745        let b = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
746        assert_eq!(a.mz, b.mz);
747    }
748
749    #[test]
750    fn a_different_frame_key_gives_different_jitter() {
751        let s = spectrum();
752        let a = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
753        let b = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[8]));
754        assert_ne!(a.mz, b.mz, "frames must not share a noise stream");
755    }
756
757    #[test]
758    fn the_jitter_stays_within_the_requested_ppm() {
759        let s = spectrum();
760        let ppm = 6.5;
761        let out = s.add_mz_noise_uniform_with_rng(ppm, false, &mut noise_rng(41, &[1]));
762        assert_eq!(out.mz.len(), s.mz.len());
763        for (before, after) in s.mz.iter().zip(out.mz.iter()) {
764            // to_resolution(6) rounds to 1e-6, so allow that on top of the ppm window.
765            let window = before * ppm / 1e6 + 1e-6;
766            assert!((after - before).abs() <= window, "{before} -> {after} exceeds {ppm} ppm");
767        }
768    }
769}
mscore/data/spectrum.rs

mscore/data/
spectrum.rs
