rustdf/sim/
containers.rs

1use mscore::chemistry::formulas::{
2    ccs_to_one_over_reduced_mobility, one_over_reduced_mobility_to_ccs,
3};
4use mscore::data::peptide::PeptideSequence;
5use mscore::data::spectrum::{MsType, MzSpectrum};
6use rand::distributions::{Distribution, Uniform};
7use serde::{Deserialize, Serialize};
8
9#[derive(Serialize, Deserialize, Debug, Clone)]
10pub struct SignalDistribution {
11    pub mean: f32,
12    pub variance: f32,
13    pub error: f32,
14    pub occurrence: Vec<u32>,
15    pub abundance: Vec<f32>,
16}
17
18impl SignalDistribution {
19    pub fn new(
20        mean: f32,
21        variance: f32,
22        error: f32,
23        occurrence: Vec<u32>,
24        abundance: Vec<f32>,
25    ) -> Self {
26        SignalDistribution {
27            mean,
28            variance,
29            error,
30            occurrence,
31            abundance,
32        }
33    }
34
35    pub fn add_noise(&self, noise_level: f32) -> Vec<f32> {
36        let mut rng = rand::thread_rng();
37        let noise_dist = Uniform::new(0.0, noise_level);
38
39        let noise: Vec<f32> = self
40            .abundance
41            .iter()
42            .map(|_| noise_dist.sample(&mut rng))
43            .collect();
44        let noise_relative: Vec<f32> = self
45            .abundance
46            .iter()
47            .zip(noise.iter())
48            .map(|(&abu, &noise)| abu * noise)
49            .collect();
50        let noised_signal: Vec<f32> = self
51            .abundance
52            .iter()
53            .zip(noise_relative.iter())
54            .map(|(&abu, &noise_rel)| abu + noise_rel)
55            .collect();
56
57        let sum_noised_signal: f32 = noised_signal.iter().sum();
58        let sum_rt_abu: f32 = self.abundance.iter().sum();
59
60        noised_signal
61            .iter()
62            .map(|&x| (x / sum_noised_signal) * sum_rt_abu)
63            .collect()
64    }
65}
66
67#[derive(Debug, Clone)]
68pub struct PeptidesSim {
69    pub protein_id: u32,
70    pub peptide_id: u32,
71    pub sequence: PeptideSequence,
72    pub proteins: String,
73    pub decoy: bool,
74    pub missed_cleavages: i8,
75    pub n_term: Option<bool>,
76    pub c_term: Option<bool>,
77    pub mono_isotopic_mass: f32,
78    pub retention_time: f32,
79    pub events: f32,
80    pub frame_start: u32,
81    pub frame_end: u32,
82    pub frame_distribution: SignalDistribution,
83}
84
85impl PeptidesSim {
86    pub fn new(
87        protein_id: u32,
88        peptide_id: u32,
89        sequence: String,
90        proteins: String,
91        decoy: bool,
92        missed_cleavages: i8,
93        n_term: Option<bool>,
94        c_term: Option<bool>,
95        mono_isotopic_mass: f32,
96        retention_time: f32,
97        events: f32,
98        frame_start: u32,
99        frame_end: u32,
100        frame_occurrence: Vec<u32>,
101        frame_abundance: Vec<f32>,
102    ) -> Self {
103        PeptidesSim {
104            protein_id,
105            peptide_id,
106            sequence: PeptideSequence::new(sequence, Some(peptide_id as i32)),
107            proteins,
108            decoy,
109            missed_cleavages,
110            n_term,
111            c_term,
112            mono_isotopic_mass,
113            retention_time,
114            events,
115            frame_start,
116            frame_end,
117            frame_distribution: SignalDistribution::new(
118                0.0,
119                0.0,
120                0.0,
121                frame_occurrence,
122                frame_abundance,
123            ),
124        }
125    }
126}
127
128#[derive(Debug, Clone)]
129pub struct WindowGroupSettingsSim {
130    pub window_group: u32,
131    pub scan_start: u32,
132    pub scan_end: u32,
133    pub isolation_mz: f32,
134    pub isolation_width: f32,
135    pub collision_energy: f32,
136}
137
138impl WindowGroupSettingsSim {
139    pub fn new(
140        window_group: u32,
141        scan_start: u32,
142        scan_end: u32,
143        isolation_mz: f32,
144        isolation_width: f32,
145        collision_energy: f32,
146    ) -> Self {
147        WindowGroupSettingsSim {
148            window_group,
149            scan_start,
150            scan_end,
151            isolation_mz,
152            isolation_width,
153            collision_energy,
154        }
155    }
156}
157
158#[derive(Debug, Clone)]
159pub struct FrameToWindowGroupSim {
160    pub frame_id: u32,
161    pub window_group: u32,
162}
163
164impl FrameToWindowGroupSim {
165    pub fn new(frame_id: u32, window_group: u32) -> Self {
166        FrameToWindowGroupSim {
167            frame_id,
168            window_group,
169        }
170    }
171}
172
173#[derive(Debug, Clone)]
174pub struct IonSim {
175    pub ion_id: u32,
176    pub peptide_id: u32,
177    pub sequence: String,
178    pub charge: i8,
179    pub relative_abundance: f32,
180    pub mobility: f32,
181    pub simulated_spectrum: MzSpectrum,
182    pub scan_distribution: SignalDistribution,
183}
184
185impl IonSim {
186    pub fn new(
187        ion_id: u32,
188        peptide_id: u32,
189        sequence: String,
190        charge: i8,
191        relative_abundance: f32,
192        mobility: f32,
193        simulated_spectrum: MzSpectrum,
194        scan_occurrence: Vec<u32>,
195        scan_abundance: Vec<f32>,
196    ) -> Self {
197        IonSim {
198            ion_id,
199            peptide_id,
200            sequence,
201            charge,
202            relative_abundance,
203            mobility,
204            simulated_spectrum,
205            scan_distribution: SignalDistribution::new(
206                0.0,
207                0.0,
208                0.0,
209                scan_occurrence,
210                scan_abundance,
211            ),
212        }
213    }
214}
215
216// --------------------------------------------------------------------------- //
217// Instrument-dispatch P1: parallel scalar-native entities.
218//
219// These mirror PeptidesSim / IonSim but hold ONLY the vendor-neutral scalar
220// physics (the trunk / "ionized sample" of INSTRUMENT_DISPATCH.md) — no
221// device-sampled occurrence/abundance vectors. They are additive: the legacy
222// SignalDistribution-bearing entities and their readers are untouched.
223//
224// Mobility ownership (plan §2.3): the trunk stores CCS (intrinsic); 1/K0 is
225// derived per instrument under that instrument's mobility environment. For
226// legacy DBs that persisted only 1/K0, the scalar reader converts 1/K0 -> CCS
227// under a declared reference `MobilityEnv` (see handle.rs).
228// --------------------------------------------------------------------------- //
229
230/// Drift-gas mobility environment for the CCS <-> 1/K0 conversion. Defaults
231/// match the timsTOF constants used by `mscore::chemistry::formulas` (N2,
232/// 31.85 °C, 273.15 K offset).
233#[derive(Debug, Clone, Copy, PartialEq)]
234pub struct MobilityEnv {
235    pub gas_mass: f64,
236    pub temp_c: f64,
237    pub t_diff: f64,
238}
239
240impl Default for MobilityEnv {
241    fn default() -> Self {
242        MobilityEnv { gas_mass: 28.013, temp_c: 31.85, t_diff: 273.15 }
243    }
244}
245
246impl MobilityEnv {
247    /// CCS for an ion observed at `one_over_k0` (legacy 1/K0 -> trunk CCS).
248    /// `charge` is clamped to >= 1 (the conversion is only defined for real
249    /// ions; the pipeline never produces charge < 1, asserted in debug).
250    pub fn ccs_from_inv_mobility(&self, one_over_k0: f64, mz: f64, charge: i8) -> f64 {
251        debug_assert!(charge >= 1, "ion charge must be >= 1, got {charge}");
252        one_over_reduced_mobility_to_ccs(
253            one_over_k0,
254            mz,
255            charge.max(1) as u32,
256            self.gas_mass,
257            self.temp_c,
258            self.t_diff,
259        )
260    }
261    /// 1/K0 for an ion of `ccs` under this environment (trunk CCS -> device 1/K0).
262    /// `charge` is clamped to >= 1 (see `ccs_from_inv_mobility`).
263    pub fn inv_mobility_from_ccs(&self, ccs: f64, mz: f64, charge: i8) -> f64 {
264        debug_assert!(charge >= 1, "ion charge must be >= 1, got {charge}");
265        ccs_to_one_over_reduced_mobility(
266            ccs,
267            mz,
268            charge.max(1) as u32,
269            self.gas_mass,
270            self.temp_c,
271            self.t_diff,
272        )
273    }
274}
275
276/// Scalar-native peptide: trunk physics with no frame-occurrence vectors.
277#[derive(Debug, Clone)]
278pub struct PeptideScalar {
279    pub protein_id: u32,
280    pub peptide_id: u32,
281    pub sequence: PeptideSequence,
282    pub proteins: String,
283    pub decoy: bool,
284    pub missed_cleavages: i8,
285    pub n_term: Option<bool>,
286    pub c_term: Option<bool>,
287    pub mono_isotopic_mass: f32,
288    /// Predicted RT apex (seconds) from the GRU predictor — provenance.
289    pub retention_time: f32,
290    /// EMG location parameter (`mu`, seconds) — the value the time projection
291    /// integrates around. NOT equal to `retention_time`: the legacy pipeline
292    /// derives it via `estimate_mu_from_mode_emg(rt_apex, sigma, lambda)` and
293    /// stores it as `rt_mu`. Falls back to `retention_time` when absent.
294    ///
295    /// Held as f64 (the DB column is REAL): the Python column writer computed the
296    /// occurrence/abundance distributions from these f64 params, so the projector
297    /// must read them at full precision to byte-reproduce the columns (an f32
298    /// round-trip straddles the 4-decimal rounding / remove_epsilon boundary).
299    pub rt_mu: f64,
300    pub rt_sigma: f64,
301    pub rt_lambda: f64,
302    pub events: f32,
303    /// Reserved per-analyte condition override (NULL -> the run's single row).
304    pub condition_id: Option<i64>,
305}
306
307/// Scalar-native ion: trunk physics with no scan-occurrence vectors. Stores CCS
308/// as canonical; 1/K0 is derived per instrument via `inv_mobility`.
309#[derive(Debug, Clone)]
310pub struct IonScalar {
311    pub ion_id: u32,
312    pub peptide_id: u32,
313    pub sequence: String,
314    pub charge: i8,
315    pub relative_abundance: f32,
316    pub mz: f64,
317    pub ccs: f64,
318    /// Legacy mobility spread (conformer width) in **1/K0 units** as stored — not
319    /// CCS-space (renamed from a misleading `ccs_std`). Transforming the spread
320    /// into CCS space is deferred; until then read it as a 1/K0 std.
321    ///
322    /// Held as f64 (DB column is REAL) so the scan-distribution projection reads
323    /// it at the precision the column writer used (see PeptideScalar::rt_mu).
324    pub inv_mobility_std: f64,
325    /// Isotope composition (m/z + relative intensity), pre-detector.
326    pub simulated_spectrum: MzSpectrum,
327    pub condition_id: Option<i64>,
328}
329
330impl IonScalar {
331    /// Derive this ion's 1/K0 under the given instrument mobility environment.
332    pub fn inv_mobility(&self, env: &MobilityEnv) -> f64 {
333        env.inv_mobility_from_ccs(self.ccs, self.mz, self.charge)
334    }
335}
336
337#[derive(Debug, Clone)]
338pub struct ScansSim {
339    pub scan: u32,
340    pub mobility: f32,
341}
342
343impl ScansSim {
344    pub fn new(scan: u32, mobility: f32) -> Self {
345        ScansSim { scan, mobility }
346    }
347}
348
349#[derive(Debug, Clone)]
350pub struct FramesSim {
351    pub frame_id: u32,
352    pub time: f32,
353    pub ms_type: i64,
354}
355
356impl FramesSim {
357    pub fn new(frame_id: u32, time: f32, ms_type: i64) -> Self {
358        FramesSim {
359            frame_id,
360            time,
361            ms_type,
362        }
363    }
364    pub fn parse_ms_type(&self) -> MsType {
365        match self.ms_type {
366            0 => MsType::Precursor,
367            8 => MsType::FragmentDda,
368            9 => MsType::FragmentDia,
369            _ => MsType::Unknown,
370        }
371    }
372}
373
374pub struct FragmentIonSim {
375    pub peptide_id: u32,
376    pub ion_id: u32,
377    pub collision_energy: f64,
378    pub charge: i8,
379    pub indices: Vec<u32>,
380    pub values: Vec<f64>,
381}
382
383impl FragmentIonSim {
384    pub fn new(
385        peptide_id: u32,
386        ion_id: u32,
387        collision_energy: f64,
388        charge: i8,
389        indices: Vec<u32>,
390        values: Vec<f64>,
391    ) -> Self {
392        FragmentIonSim {
393            peptide_id,
394            ion_id,
395            charge,
396            collision_energy,
397            indices,
398            values,
399        }
400    }
401
402    pub fn to_dense(&self, length: usize) -> Vec<f64> {
403        let mut dense = vec![0.0; length];
404        for (i, &idx) in self.indices.iter().enumerate() {
405            dense[idx as usize] = self.values[i];
406        }
407        dense
408    }
409}
410
411/// Mode for quad-selection dependent isotope transmission calculation.
412#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
413pub enum IsotopeTransmissionMode {
414    /// Disabled - no transmission-dependent calculation
415    None,
416    /// Precursor-based scaling - calculate transmission factor from precursor isotope
417    /// distribution and apply uniform scaling to all fragment intensities.
418    /// This is computationally efficient and captures the main intensity reduction effect.
419    PrecursorScaling,
420    /// Per-fragment calculation - calculate transmission-dependent isotope distribution
421    /// for each individual fragment ion based on its complementary fragment.
422    /// This is more accurate but computationally expensive.
423    /// Implements the algorithm from OpenMS's CoarseIsotopePatternGenerator.
424    PerFragment,
425}
426
427impl Default for IsotopeTransmissionMode {
428    fn default() -> Self {
429        Self::None
430    }
431}
432
433/// Configuration for quad-selection dependent isotope transmission.
434///
435/// When enabled, fragment ion isotope distributions and/or intensities are adjusted
436/// based on which precursor isotopes were transmitted through the quadrupole isolation
437/// window.
438#[derive(Debug, Clone, Serialize, Deserialize)]
439pub struct IsotopeTransmissionConfig {
440    /// Mode for transmission-dependent calculation
441    pub mode: IsotopeTransmissionMode,
442    /// Minimum probability threshold for isotope transmission
443    pub min_probability: f64,
444    /// Maximum number of isotope peaks to consider
445    pub max_isotopes: usize,
446    /// Minimum fraction of precursor ions that survive fragmentation intact (0.0-1.0)
447    pub precursor_survival_min: f64,
448    /// Maximum fraction of precursor ions that survive fragmentation intact (0.0-1.0)
449    pub precursor_survival_max: f64,
450}
451
452impl Default for IsotopeTransmissionConfig {
453    fn default() -> Self {
454        Self {
455            mode: IsotopeTransmissionMode::None,
456            min_probability: 0.5,
457            max_isotopes: 10,
458            precursor_survival_min: 0.0,
459            precursor_survival_max: 0.0,
460        }
461    }
462}
463
464impl IsotopeTransmissionConfig {
465    pub fn new(
466        mode: IsotopeTransmissionMode,
467        min_probability: f64,
468        max_isotopes: usize,
469        precursor_survival_min: f64,
470        precursor_survival_max: f64,
471    ) -> Self {
472        Self {
473            mode,
474            min_probability,
475            max_isotopes,
476            precursor_survival_min,
477            precursor_survival_max,
478        }
479    }
480
481    /// Create config with precursor scaling mode
482    pub fn precursor_scaling(min_probability: f64) -> Self {
483        Self {
484            mode: IsotopeTransmissionMode::PrecursorScaling,
485            min_probability,
486            max_isotopes: 10,
487            precursor_survival_min: 0.0,
488            precursor_survival_max: 0.0,
489        }
490    }
491
492    /// Create config with per-fragment mode
493    pub fn per_fragment(min_probability: f64, max_isotopes: usize) -> Self {
494        Self {
495            mode: IsotopeTransmissionMode::PerFragment,
496            min_probability,
497            max_isotopes,
498            precursor_survival_min: 0.0,
499            precursor_survival_max: 0.0,
500        }
501    }
502
503    /// Check if precursor survival is enabled
504    pub fn has_precursor_survival(&self) -> bool {
505        self.precursor_survival_max > 0.0
506    }
507
508    /// Check if any transmission mode is enabled
509    pub fn is_enabled(&self) -> bool {
510        self.mode != IsotopeTransmissionMode::None
511    }
512
513    /// Apply instrument-capability gating (P5e). Mobility-dependent quadrupole
514    /// isotope transmission (PrecursorScaling / PerFragment) is a Bruker timsTOF
515    /// behaviour; an instrument that lacks it (e.g. a no-IMS Astral —
516    /// `has_quad_isotope_transmission = false`) must NOT apply that scaling, so
517    /// force the mode to `None`. Bruker (flag true) is returned unchanged, so the
518    /// rendered output is byte-identical. (The m/z-isolation vs scan/mobility
519    /// transmission split is gated by `has_tims_mobility` and lands in P6 with the
520    /// Thermo acquisition windows — see the instrument-dispatch plan.)
521    pub fn gated_by(&self, capabilities: crate::sim::scheme::InstrumentCapabilities) -> Self {
522        if capabilities.has_quad_isotope_transmission {
523            self.clone()
524        } else {
525            let mut gated = self.clone();
526            gated.mode = IsotopeTransmissionMode::None;
527            gated
528        }
529    }
530}
531
532#[cfg(test)]
533mod scalar_entity_tests {
534    use super::*;
535    use crate::sim::scheme::InstrumentCapabilities;
536
537    #[test]
538    fn isotope_config_gated_by_capabilities() {
539        let cfg = IsotopeTransmissionConfig {
540            mode: IsotopeTransmissionMode::PerFragment,
541            min_probability: 0.5,
542            max_isotopes: 10,
543            precursor_survival_min: 0.0,
544            precursor_survival_max: 0.0,
545        };
546        // Bruker (default): quad isotope transmission present -> unchanged.
547        let bruker = cfg.gated_by(InstrumentCapabilities::default());
548        assert_eq!(bruker.mode, IsotopeTransmissionMode::PerFragment);
549        assert!(bruker.is_enabled());
550        // No-quad-isotope instrument (e.g. Astral): mode forced to None.
551        let astral = cfg.gated_by(InstrumentCapabilities {
552            has_tims_mobility: false,
553            has_quad_isotope_transmission: false,
554        });
555        assert_eq!(astral.mode, IsotopeTransmissionMode::None);
556        assert!(!astral.is_enabled());
557    }
558
559    #[test]
560    fn mobility_env_ccs_inv_mobility_round_trips() {
561        // 1/K0 -> CCS -> 1/K0 must be identity under the same environment
562        // (the legacy-1/K0 -> trunk-CCS migration must lose nothing).
563        let env = MobilityEnv::default();
564        let (mz, charge, one_over_k0) = (1000.0_f64, 2_i8, 0.85_f64);
565        let ccs = env.ccs_from_inv_mobility(one_over_k0, mz, charge);
566        let back = env.inv_mobility_from_ccs(ccs, mz, charge);
567        assert!((back - one_over_k0).abs() < 1e-9, "round-trip drift: {back} vs {one_over_k0}");
568    }
569
570    #[test]
571    fn ion_scalar_derives_inv_mobility_per_env() {
572        let warm = MobilityEnv { gas_mass: 28.013, temp_c: 40.0, t_diff: 273.15 };
573        let cold = MobilityEnv { gas_mass: 28.013, temp_c: 20.0, t_diff: 273.15 };
574        let ion = IonScalar {
575            ion_id: 1,
576            peptide_id: 1,
577            sequence: "PEPTIDEK".to_string(),
578            charge: 2,
579            relative_abundance: 1.0,
580            mz: 500.0,
581            ccs: 350.0,
582            inv_mobility_std: 0.0,
583            simulated_spectrum: MzSpectrum::new(vec![500.0], vec![1.0]),
584            condition_id: None,
585        };
586        // Same CCS yields different 1/K0 in different drift environments — the
587        // whole point of storing CCS in the trunk, not 1/K0.
588        assert!(ion.inv_mobility(&warm) != ion.inv_mobility(&cold));
589    }
590}
rustdf/sim/containers.rs

rustdf/sim/
containers.rs
