perf: remove the serial hot spots, and make a run reproducible from its config

imspy_connector/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "imspy-connector"
-version = "0.4.2"
+version = "0.4.4"
 edition = "2021"
 
 [lib]

imspy_connector/pyproject.toml

@@ -7,7 +7,7 @@ name = "imspy_connector"
 dependencies = [
     "opentims-bruker-bridge>=1.1.0",
 ]
-version = "0.4.2"
+version = "0.4.4"
 requires-python = ">=3.11"
 classifiers = [
     "Programming Language :: Rust",

imspy_connector/src/py_dataset.rs

@@ -5,12 +5,26 @@ use rustdf::data::utility::{zstd_compress, zstd_decompress, reconstruct_compress
 
 use crate::py_tims_frame::{PyTimsFrame};
 use crate::py_tims_slice::PyTimsSlice;
-use numpy::{IntoPyArray, PyArray1};
+use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
+use pyo3::exceptions::PyRuntimeError;
+use rustdf::data::handle::IndexConverter;
+use rustdf::data::utility::build_compressed_frames as rust_build_compressed_frames;
 use pyo3::types::PyList;
 use pyo3::{PyResult, Python};
 use rustdf::data::acquisition::AcquisitionMode;
 use rustdf::data::handle::TimsData;
 
+/// Borrow each numpy array as a contiguous slice, keeping the arrays' lifetime.
+fn readonly_slices<'a>(arrays: &'a [PyReadonlyArray1<'_, f64>]) -> PyResult<Vec<&'a [f64]>> {
+    arrays
+        .iter()
+        .map(|a| {
+            a.as_slice()
+                .map_err(|_| PyRuntimeError::new_err("build_compressed_frames: arrays must be contiguous"))
+        })
+        .collect()
+}
+
 #[pyclass]
 pub struct PyTimsDataset {
     inner: TimsDataset,
@@ -119,6 +133,60 @@ impl PyTimsDataset {
         self.inner.loader.get_index_converter().inverse_mobility_to_scan(frame_id, &inverse_mobility_values)
     }
 
+    /// Run the whole per-frame TDF writer pipeline for a batch of frames in parallel:
+    /// m/z -> TOF, 1/K0 -> scan, `(scan, tof)` dedup, interleave and zstd.
+    ///
+    /// Returns `(num_peaks, max_intensity, summed_intensity, compressed_bytes)` per frame, in input
+    /// order; the statistics are taken after the dedup, which is what the `Frames` row records.
+    ///
+    /// Raises if the dataset was opened with the Bruker SDK and the file carries no calibration
+    /// tables the SDK-free converter can use — the SDK cannot be called from several threads at
+    /// once, so the caller must fall back to the per-frame path in that case.
+    #[pyo3(signature = (frame_ids, mz, mobility, intensity, max_scans, compression_level=0, num_threads=4))]
+    pub fn build_compressed_frames(
+        &self,
+        py: Python<'_>,
+        frame_ids: Vec<u32>,
+        mz: Vec<PyReadonlyArray1<'_, f64>>,
+        mobility: Vec<PyReadonlyArray1<'_, f64>>,
+        intensity: Vec<PyReadonlyArray1<'_, f64>>,
+        max_scans: u32,
+        compression_level: i32,
+        num_threads: usize,
+    ) -> PyResult<Vec<(u32, u32, u64, Vec<u8>)>> {
+        let sdk_free = self.inner.loader.sdk_free_converter();
+        let converter: Option<&(dyn IndexConverter + Sync)> = match sdk_free.as_ref() {
+            Some(c) => Some(c),
+            None => self.inner.loader.sync_index_converter(),
+        };
+        let Some(converter) = converter else {
+            return Err(PyRuntimeError::new_err(
+                "UNAVAILABLE: build_compressed_frames needs a converter that reproduces the \
+                 Bruker SDK exactly; this .d either has no usable calibration tables or uses \
+                 MzCalibration ModelType 2, which is only accurate to ~2 ppm. Use the per-frame \
+                 writer path.",
+            ));
+        };
+
+        let mz_slices = readonly_slices(&mz)?;
+        let im_slices = readonly_slices(&mobility)?;
+        let it_slices = readonly_slices(&intensity)?;
+
+        let frames = py
+            .detach(|| {
+                rust_build_compressed_frames(
+                    converter, &frame_ids, &mz_slices, &im_slices, &it_slices,
+                    max_scans, compression_level, num_threads,
+                )
+            })
+            .map_err(PyRuntimeError::new_err)?;
+
+        Ok(frames
+            .into_iter()
+            .map(|f| (f.num_peaks, f.max_intensity, f.summed_intensity, f.data))
+            .collect())
+    }
+
     #[staticmethod]
     pub fn compress_bytes_zstd(bytes: Vec<u8>, compression_level: i32) -> Vec<u8> {
         let result = zstd_compress(&bytes, compression_level).unwrap();

imspy_connector/src/py_dia.rs

@@ -2226,6 +2226,34 @@ impl PyTimsDatasetDIA {
         PyTimsFrame::from_inner(self.inner.sample_fragment_signal(num_frames, window_group, max_intensity, take_probability))
     }
 
+    /// Batched, parallel, seeded reference-noise overlay (see rustdf `overlay_reference_noise`).
+    /// `window_groups[i]` is the DIA window group of frame i, or None for a precursor frame.
+    #[pyo3(signature = (frames, window_groups, num_precursor_frames, num_fragment_frames, max_intensity_precursor, max_intensity_fragment, take_precursor, take_fragment, seed, num_threads))]
+    pub fn overlay_reference_noise(
+        &self,
+        frames: Vec<PyTimsFrame>,
+        window_groups: Vec<Option<u32>>,
+        num_precursor_frames: usize,
+        num_fragment_frames: usize,
+        max_intensity_precursor: f64,
+        max_intensity_fragment: f64,
+        take_precursor: f64,
+        take_fragment: f64,
+        seed: u64,
+        num_threads: usize,
+    ) -> Vec<PyTimsFrame> {
+        let inner: Vec<mscore::timstof::frame::TimsFrame> = frames.into_iter().map(|f| f.inner).collect();
+        self.inner
+            .overlay_reference_noise(
+                inner, window_groups, num_precursor_frames, num_fragment_frames,
+                max_intensity_precursor, max_intensity_fragment, take_precursor, take_fragment,
+                seed, num_threads,
+            )
+            .into_iter()
+            .map(PyTimsFrame::from_inner)
+            .collect()
+    }
+
     /// Convenience helper for Python-side planner
     #[getter]
     pub fn max_global_num_scans(&self) -> usize {

imspy_connector/src/py_peptide.rs

@@ -1,5 +1,7 @@
 use std::collections::{HashMap};
 use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+use rayon::prelude::*;
 
 use mscore::data::peptide::{FragmentType, PeptideSequence, PeptideProductIon,
                             PeptideProductIonSeries, PeptideProductIonSeriesCollection, PeptideIon};
@@ -345,12 +347,36 @@ impl PyPeptideProductIon {
     }
 }
 
+/// Monoisotopic mass for a batch of UniMod sequences, in parallel.
+///
+/// Constructing one `PyPeptideSequence` per peptide just to read `.mono_isotopic_mass` costs about
+/// 200 us each — 250 000 peptides took 48 s in a single Python loop, which was the largest remaining
+/// serial cost in the simulation. The work is per-sequence independent, so do the whole batch over a
+/// rayon pool and hand back one array.
+#[pyfunction]
+#[pyo3(signature = (sequences, num_threads=4))]
+pub fn mono_isotopic_masses(py: Python<'_>, sequences: Vec<String>, num_threads: usize) -> PyResult<Vec<f64>> {
+    let pool = rayon::ThreadPoolBuilder::new()
+        .num_threads(num_threads.max(1))
+        .build()
+        .map_err(|e| pyo3::exceptions::PyRuntimeError::new_err(format!("thread pool: {e}")))?;
+    Ok(py.detach(|| {
+        pool.install(|| {
+            sequences
+                .par_iter()
+                .map(|s| PeptideSequence::new(s.clone(), None).mono_isotopic_mass())
+                .collect()
+        })
+    }))
+}
+
 #[pymodule]
 pub fn py_peptide(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<PyPeptideSequence>()?;
     m.add_class::<PyPeptideIon>()?;
     m.add_class::<PyPeptideProductIon>()?;
     m.add_class::<PyPeptideProductIonSeries>()?;
     m.add_class::<PyPeptideProductIonSeriesCollection>()?;
+    m.add_function(wrap_pyfunction!(mono_isotopic_masses, m)?)?;
     Ok(())
 }

imspy_connector/src/py_simulation.rs

@@ -202,6 +202,12 @@ pub struct PyTimsTofSyntheticsPrecursorFrameBuilder {
 
 #[pymethods]
 impl PyTimsTofSyntheticsPrecursorFrameBuilder {
+
+    /// Seed the m/z jitter so a run's noise is reproducible. Without it the jitter comes from
+    /// `thread_rng` and the written spectra depend on thread scheduling.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.inner.set_noise_seed(seed);
+    }
     #[new]
     pub fn new(db_path: &str) -> Self {
         let path = std::path::Path::new(db_path);
@@ -238,6 +244,12 @@ pub struct PyTimsTofSyntheticsFrameBuilderDIA {
 
 #[pymethods]
 impl PyTimsTofSyntheticsFrameBuilderDIA {
+
+    /// Seed the m/z jitter so a run's noise is reproducible. Without it the jitter comes from
+    /// `thread_rng` and the written spectra depend on thread scheduling.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.inner.set_noise_seed(seed);
+    }
     #[new]
     #[pyo3(signature = (db_path, with_annotations, num_threads, isotope_config=None, projection_mode=None, target_p=0.999, frame_step_size=0.001, scan_step_size=0.0001, n_steps=1000, remove_epsilon=1e-4))]
     pub fn new(
@@ -337,6 +349,12 @@ pub struct PyTimsTofSyntheticsFrameBuilderDDA {
 
 #[pymethods]
 impl PyTimsTofSyntheticsFrameBuilderDDA {
+
+    /// Seed the m/z jitter so a run's noise is reproducible. Without it the jitter comes from
+    /// `thread_rng` and the written spectra depend on thread scheduling.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.inner.set_noise_seed(seed);
+    }
     #[new]
     #[pyo3(signature = (db_path, with_annotations, num_threads, isotope_config=None))]
     pub fn new(
@@ -425,6 +443,12 @@ pub struct PyTimsTofLazyFrameBuilderDIA {
 
 #[pymethods]
 impl PyTimsTofLazyFrameBuilderDIA {
+
+    /// Seed the m/z jitter so a run's noise is reproducible. Without it the jitter comes from
+    /// `thread_rng` and the written spectra depend on thread scheduling.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.inner.set_noise_seed(seed);
+    }
     #[new]
     #[pyo3(signature = (db_path, num_threads=4, projection_mode=None, target_p=0.999, frame_step_size=0.001, scan_step_size=0.0001, n_steps=1000, remove_epsilon=1e-4))]
     pub fn new(

mscore/src/algorithm/peptide.rs

@@ -11,6 +11,7 @@ use rayon::ThreadPoolBuilder;
 use regex::Regex;
 use statrs::distribution::{Binomial, Discrete};
 use std::collections::HashMap;
+use std::sync::LazyLock;
 
 /// calculate the monoisotopic mass of a peptide sequence
 ///
@@ -33,10 +34,20 @@ use std::collections::HashMap;
 /// let mass_quantized = (mass * 1e6).round() as i32;
 /// assert_eq!(mass_quantized, 936418877);
 /// ```
+/// Tables and patterns hoisted out of the per-peptide path.
+///
+/// `calculate_peptide_mono_isotopic_mass` rebuilt the 2 073-entry UniMod table, rebuilt the amino
+/// acid table and recompiled the regex on **every** call. That put one mass calculation at ~190 us,
+/// and 250 000 peptides at ~48 s — the largest remaining serial cost in the simulation. Built once.
+static AA_MASSES: LazyLock<HashMap<&'static str, f64>> = LazyLock::new(amino_acid_masses);
+static UNIMOD_MASSES: LazyLock<HashMap<u32, f64>> = LazyLock::new(unimod_modifications_mass_numerical);
+static UNIMOD_PATTERN: LazyLock<Regex> =
+    LazyLock::new(|| Regex::new(r"\[UNIMOD:(\d+)]").unwrap());
+
 pub fn calculate_peptide_mono_isotopic_mass(peptide_sequence: &PeptideSequence) -> f64 {
-    let amino_acid_masses = amino_acid_masses();
-    let modifications_mz_numerical = unimod_modifications_mass_numerical();
-    let pattern = Regex::new(r"\[UNIMOD:(\d+)]").unwrap();
+    let amino_acid_masses = &*AA_MASSES;
+    let modifications_mz_numerical = &*UNIMOD_MASSES;
+    let pattern = &*UNIMOD_PATTERN;
 
     let sequence = peptide_sequence.sequence.as_str();
 
@@ -102,7 +113,7 @@ pub fn calculate_peptide_product_ion_mono_isotopic_mass(sequence: &str, kind: Fr
         return 0.0;
     }
 
-    let amino_acid_masses = amino_acid_masses();
+    let amino_acid_masses = &*AA_MASSES;
 
     // Add up raw amino acid masses and potential modifications
     let mass_sequence: f64 = sequence