From d86f902bda09e6dbb31bb038dfd850f2b182caf3 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sat, 12 Sep 2026 22:39:45 +0200
Subject: [PATCH 01/14] perf(timsim): remove serial hot spots in the TDF writer
 and add stage telemetry

Profiling a 250k-peptide DIA run on a 120-core box showed ~90 % of wall time
does not scale with threads; frame assembly (70 % of wall) spent ~85 % of its
time on the main thread, dominated by the per-frame TDF writer.

- tdf.py: dedup of (scan, tof) cells via a packed uint64 key + 1-D np.unique
  instead of np.unique(axis=0) + lexsort (10-15x faster, byte-identical
  output); keep analysis.tdf_bin open for the whole run instead of
  open/append/close per frame (34k syscalls, expensive on NFS). Handle is
  closed in write_frame_meta_data() and re-opened transparently if needed.
- utility.py: python_list_to_json_string rounds the whole array once instead
  of a NumPy scalar round per element (called once per peptide/ion in four
  stages); output byte-identical, float32 inputs keep float32 rounding.
- assemble_frames.py: per-batch timers for build / noise / write, logged and
  returned.
- simulator.py: StageTracker records wall, CPU seconds (all threads), avg
  cores and peak RSS per stage, logs each transition, prints a stage table
  and writes timings.json next to synthetic_data.db. Sub-laps inside the
  FASTA stage (digest+rt_filter, sample_peptides).
- tests: parity tests against the previous implementations, handle
  lifecycle, tracker output.

Smoke run (5.7k frames, 2k peptides, 16 threads): assembly 54 s -> 38 s,
total 156 s -> 138 s; .d reads back with imspy.
---
 packages/imspy-simulation/pyproject.toml      |   2 +-
 .../src/imspy_simulation/tdf.py               |  75 +++++----
 .../timsim/jobs/assemble_frames.py            |  29 +++-
 .../src/imspy_simulation/timsim/simulator.py  | 149 +++++++++++++++---
 .../src/imspy_simulation/utility.py           |  15 +-
 .../tests/test_perf_serial_hotspots.py        | 147 +++++++++++++++++
 6 files changed, 366 insertions(+), 51 deletions(-)
 create mode 100644 packages/imspy-simulation/tests/test_perf_serial_hotspots.py

diff --git a/packages/imspy-simulation/pyproject.toml b/packages/imspy-simulation/pyproject.toml
index 777b65783..a9acc5dc1 100644
--- a/packages/imspy-simulation/pyproject.toml
+++ b/packages/imspy-simulation/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "imspy-simulation"
-version = "0.4.2"
+version = "0.4.3"
 description = "TimsTOF data simulation tools for proteomics."
 authors = [
     { name = "theGreatHerrLebert", email = "davidteschner@googlemail.com" }
diff --git a/packages/imspy-simulation/src/imspy_simulation/tdf.py b/packages/imspy-simulation/src/imspy_simulation/tdf.py
index e9cb518b3..7850d2943 100644
--- a/packages/imspy-simulation/src/imspy_simulation/tdf.py
+++ b/packages/imspy-simulation/src/imspy_simulation/tdf.py
@@ -133,6 +133,26 @@ def _locked_database_error(db_path, exc: Exception) -> RuntimeError:
     )
 
 
+def dedup_scan_tof(scan: NDArray, tof: NDArray, intensity: NDArray) -> tuple[NDArray, NDArray, NDArray]:
+    """Merge peaks that fall on the same (scan, tof) cell and return them sorted by scan, then tof.
+
+    m/z -> TOF is not injective, so several simulated peaks can land on one TOF index; their
+    intensities are summed. The result is byte-identical to the previous
+    ``np.unique(np.stack((scan, tof)), axis=0)`` + ``np.lexsort`` implementation, but ~10-15x faster:
+    the pair is packed into one uint64 key (scan in the high 32 bits) so a single 1-D sort orders by
+    scan first and tof second, which is exactly the lexicographic order the TDF encoder needs.
+    """
+    scan = np.asarray(scan, dtype=np.uint32)
+    tof = np.asarray(tof, dtype=np.uint32)
+    intensity = np.asarray(intensity)
+    key = (scan.astype(np.uint64) << np.uint64(32)) | tof.astype(np.uint64)
+    unique_key, inverse = np.unique(key, return_inverse=True)
+    summed = np.bincount(inverse.ravel(), weights=intensity)
+    out_scan = (unique_key >> np.uint64(32)).astype(np.uint32)
+    out_tof = (unique_key & np.uint64(0xFFFFFFFF)).astype(np.uint32)
+    return out_scan, out_tof, summed.astype(np.uint32)
+
+
 class TDFWriter:
     def __init__(self, helper_handle: TimsDataset, path: str = "./", exp_name: str = "RAW.d", offset_bytes: int = 64, verbose: bool=False, use_rust_compression: bool=False, expect_existing: bool = False) -> None:
 
@@ -161,6 +181,9 @@ def __init__(self, helper_handle: TimsDataset, path: str = "./", exp_name: str =
         ) not in ("0", "", "false", "False")
 
         self.__conn_native = None
+        # Binary output handle; opened lazily on first write_frame() and kept open for the
+        # whole run instead of open/append/close per frame (34k+ syscalls, costly on NFS).
+        self._bin_fh = None
         self._setup_connections()
 
     def _setup_connections(self) -> None:
@@ -387,27 +410,8 @@ def compress_frame(self, frame: TimsFrame, only_frame_one: bool = False) -> (NDA
         scan = self.inv_mobility_to_scan(i, frame.mobility).astype(np.uint32)
         intensity = frame.intensity.astype(np.uint32)
 
-        # Since, mz -> tof is not bijective, we need to check for duplicates
-        # stack scan and tof to form a 2D array for unique grouping
-        scan_tof = np.stack((scan, tof), axis=1)
-
-        # get unique (scan, tof) pairs and their inverse indices
-        unique_pairs, inverse_indices = np.unique(scan_tof, axis=0, return_inverse=True)
-
-        # sum intensities for each unique (scan, tof) pair
-        summed_intensity = np.bincount(inverse_indices, weights=intensity)
-
-        # now split back scan and tof
-        unique_scan = unique_pairs[:, 0]
-        unique_tof = unique_pairs[:, 1]
-
-        # sort first by scan, then by tof
-        sort_idx = np.lexsort((unique_tof, unique_scan))
-
-        # final sorted arrays
-        scan = unique_scan[sort_idx]
-        tof = unique_tof[sort_idx]
-        intensity = summed_intensity[sort_idx].astype(np.uint32)
+        # Since mz -> tof is not bijective, merge duplicate (scan, tof) cells and sort by scan, tof.
+        scan, tof, intensity = dedup_scan_tof(scan, tof, intensity)
 
         # get the real data as interleaved bytes (Rust encoder or NumPy/Numba)
         if self.use_rust_compression:
@@ -440,13 +444,27 @@ def write_frame(self, frame: TimsFrame, scan_mode: int, only_frame_one: bool = F
                 only_frame_one
         ))
 
-        with open(self.binary_file, "ab") as bin_file:
-            bin_file.write(
-                (len(compressed_data) + 8).to_bytes(4, "little", signed=False)
-            )
-            bin_file.write(int(self.helper_handle.num_scans).to_bytes(4, "little", signed=False))
-            bin_file.write(compressed_data)
-            self.position = bin_file.tell()
+        bin_file = self._binary_handle()
+        bin_file.write(
+            (len(compressed_data) + 8).to_bytes(4, "little", signed=False)
+        )
+        bin_file.write(int(self.helper_handle.num_scans).to_bytes(4, "little", signed=False))
+        bin_file.write(compressed_data)
+        self.position = bin_file.tell()
+
+    def _binary_handle(self):
+        """Return the persistent append handle for analysis.tdf_bin, opening it on first use."""
+        if self._bin_fh is None or self._bin_fh.closed:
+            self._bin_fh = open(self.binary_file, "ab")
+        return self._bin_fh
+
+    def close_binary(self) -> None:
+        """Flush and close the analysis.tdf_bin handle (idempotent). Called once all frames are
+        written; a later write_frame() re-opens it transparently."""
+        if self._bin_fh is not None and not self._bin_fh.closed:
+            self._bin_fh.flush()
+            self._bin_fh.close()
+        self._bin_fh = None
 
     def get_frame_meta_data(self) -> pd.DataFrame:
         return pd.DataFrame(self.frame_meta_data)
@@ -466,6 +484,7 @@ def write_frame_meta_data(self) -> None:
         prevents the noise-pipeline from producing duplicate-Id rows in
         the first place. The check below is the defence-in-depth.
         """
+        self.close_binary()
         meta_df = self.get_frame_meta_data()
         validate_frames_id_uniqueness(meta_df)
         self._create_table(self.conn, meta_df, "Frames")
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
index 26476d4dc..b1f1267b7 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
@@ -1,6 +1,7 @@
 from pathlib import Path
 from typing import Optional
 import logging
+import time
 
 import pandas as pd
 from tqdm import tqdm
@@ -123,12 +124,17 @@ def assemble_frames(
     logger.info(f'Fragment m/z noise: {mz_noise_fragment}')
     logger.info(f'Fragment noise PPM: {fragment_noise_ppm}')
 
+    # per-phase wall clock inside the batch loop (build = parallel Rust, noise + write = main thread)
+    t_build = t_noise = t_write = 0.0
+    n_written = 0
+
     # go over all frames in batches
     for b in tqdm(range(num_batches), total=num_batches, desc='frame assembly', ncols=100):
         start_index = b * batch_size
         stop_index = (b + 1) * batch_size
         ids = frame_ids[start_index:stop_index]
 
+        t0 = time.perf_counter()
         built_frames = frame_builder.build_frames(
             ids,
             mz_noise_precursor=mz_noise_precursor,
@@ -139,6 +145,7 @@ def assemble_frames(
             num_threads=num_threads,
             fragment=fragment,
         )
+        t1 = time.perf_counter()
 
         if superimpose_on_reference:
             built_frames = superimpose_reference_frames(
@@ -158,13 +165,31 @@ def assemble_frames(
                 acquisition_mode=acquisition_builder.acquisition_mode.mode,
                 pasef_meta=pasef_meta,
             )
+        t2 = time.perf_counter()
 
         for frame in built_frames:
             if acquisition_builder.acquisition_mode.mode == 'DDA':
                 acquisition_builder.tdf_writer.write_frame(frame, scan_mode=8)
             else:
                 acquisition_builder.tdf_writer.write_frame(frame, scan_mode=9)
-
+        t3 = time.perf_counter()
+
+        t_build += t1 - t0
+        t_noise += t2 - t1
+        t_write += t3 - t2
+        n_written += len(built_frames)
+
+    timings = {
+        'frames': n_written,
+        'batches': num_batches,
+        'build_s': round(t_build, 2),
+        'noise_s': round(t_noise, 2),
+        'write_s': round(t_write, 2),
+    }
+    logger.info(
+        f"frame assembly phases: build {t_build:.1f} s | noise/superimpose {t_noise:.1f} s | "
+        f"convert+compress+write {t_write:.1f} s  ({n_written} frames, {num_batches} batches)"
+    )
     logger.info('Writing frame meta data to database ...')
 
     # write frame meta data to database
@@ -216,3 +241,5 @@ def assemble_frames(
             acquisition_builder.synthetics_handle.get_table('pasef_meta'),
             id_mapping=id_mapping
         )
+
+    return timings
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
index c024b2041..408df9e91 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
@@ -9,6 +9,8 @@
 import argparse
 import logging
 import time
+import json
+import resource
 from pathlib import Path
 from dataclasses import dataclass, field
 from typing import Dict, Optional
@@ -62,6 +64,94 @@ def format_duration(seconds: float) -> str:
             return f"{hours}h {minutes}m"
 
 
+class StageTracker:
+    """Per-stage wall clock, CPU seconds (all threads) and peak RSS, logged at every stage
+    transition and dumped to ``timings.json`` next to ``synthetic_data.db``.
+
+    Wall is ``time.perf_counter``; CPU is ``resource.getrusage(RUSAGE_SELF)`` user+sys, which
+    includes rayon/torch worker threads, so ``cpu_s / wall_s`` is the average number of busy
+    cores during the stage. ``peak_rss_gb`` is the process high-water mark reached by the end of
+    the stage (monotonic). ``laps`` are optional named sub-intervals inside a stage.
+    """
+
+    def __init__(self, logger: logging.Logger) -> None:
+        self._logger = logger
+        self.records: list[dict] = []
+        self._current: Optional[dict] = None
+        self._t0 = 0.0
+        self._cpu0 = 0.0
+        self._lap_t = 0.0
+
+    @staticmethod
+    def _cpu_seconds() -> float:
+        ru = resource.getrusage(resource.RUSAGE_SELF)
+        return ru.ru_utime + ru.ru_stime
+
+    @staticmethod
+    def _peak_rss_gb() -> float:
+        return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024 ** 2  # KiB -> GiB on Linux
+
+    def begin(self, name: str) -> None:
+        self.end()
+        self._current = {'stage': name, 'laps': {}}
+        self._t0 = self._lap_t = time.perf_counter()
+        self._cpu0 = self._cpu_seconds()
+
+    def lap(self, label: str) -> None:
+        """Record the time since the stage start (or the previous lap) under ``label``."""
+        if self._current is None:
+            return
+        now = time.perf_counter()
+        self._current['laps'][label] = round(now - self._lap_t, 2)
+        self._lap_t = now
+
+    def annotate(self, **fields) -> None:
+        """Attach extra measurements (e.g. per-phase timers from a job) to the current stage."""
+        if self._current is not None:
+            self._current.update(fields)
+
+    def end(self) -> None:
+        if self._current is None:
+            return
+        wall = time.perf_counter() - self._t0
+        cpu = self._cpu_seconds() - self._cpu0
+        rec = self._current
+        rec.update(wall_s=round(wall, 2), cpu_s=round(cpu, 2),
+                   avg_cores=round(cpu / wall, 2) if wall > 0 else 0.0,
+                   peak_rss_gb=round(self._peak_rss_gb(), 2))
+        if not rec['laps']:
+            rec.pop('laps')
+        self.records.append(rec)
+        self._current = None
+        laps = f" | laps {rec['laps']}" if 'laps' in rec else ""
+        self._logger.info(
+            f"  stage done: {rec['stage']}: wall {wall:.1f} s | cpu {cpu:.1f} core-s "
+            f"(avg {rec['avg_cores']:.1f} cores) | peak RSS {rec['peak_rss_gb']:.2f} GB{laps}"
+        )
+
+    def summary_lines(self) -> list[str]:
+        total = sum(r['wall_s'] for r in self.records) or 1.0
+        rows = [(r['stage'], f"{r['wall_s']:.1f}", f"{100 * r['wall_s'] / total:.0f}%",
+                 f"{r['cpu_s']:.1f}", f"{r['avg_cores']:.1f}", f"{r['peak_rss_gb']:.2f}")
+                for r in self.records]
+        return tabulate(rows, headers=['stage', 'wall s', 'share', 'cpu s', 'avg cores', 'peak RSS GB'],
+                        tablefmt='simple').splitlines()
+
+    def dump(self, path) -> None:
+        self.end()
+        payload = {
+            'stages': self.records,
+            'total_wall_s': round(sum(r['wall_s'] for r in self.records), 2),
+            'total_cpu_s': round(sum(r['cpu_s'] for r in self.records), 2),
+            'peak_rss_gb': round(self._peak_rss_gb(), 2),
+        }
+        try:
+            with open(path, 'w') as fh:
+                json.dump(payload, fh, indent=1)
+        except OSError as e:  # never fail a finished simulation over telemetry
+            self._logger.warning(f"could not write {path}: {e}")
+
+
 @dataclass
 class SimulationStats:
     """Track statistics for simulation summary."""
@@ -1310,9 +1400,15 @@ def main():
     # Initialize timer and stats
     timer = SimulationTimer()
     stats = SimulationStats()
+    stages = StageTracker(logger)
+
+    def section(title: str) -> str:
+        """Start timing ``title`` and return its banner (closes the previous stage)."""
+        stages.begin(title)
+        return section_header(title, use_unicode)
 
     # Configure GPU
-    logger.info(section_header("GPU Configuration", use_unicode))
+    logger.info(section("GPU Configuration"))
     configure_gpu_memory(memory_limit_gb=config.gpu_memory_limit_gb, use_gpu=config.use_gpu)
 
     # Handle macOS Bruker SDK limitation
@@ -1322,7 +1418,7 @@ def main():
         use_bruker_sdk = False
 
     # Load modifications config
-    logger.info(section_header("Loading Configuration", use_unicode))
+    logger.info(section("Loading Configuration"))
     script_dir = Path(__file__).parent
     modifications_path = config.modifications
     if not modifications_path or modifications_path == "":
@@ -1360,7 +1456,7 @@ def main():
         f.write(tabulate(table_data, headers=["Argument", "Value"], tablefmt="grid"))
 
     # Build the acquisition (frames, scans, etc.)
-    logger.info(section_header("Building Acquisition", use_unicode))
+    logger.info(section("Building Acquisition"))
     if not config.silent_mode:
         logger.info("")
         logger.info(f"  Experiment:  {name}")
@@ -1486,7 +1582,7 @@ def main():
         factors = {}
 
     if config.from_existing and not resume_after:
-        logger.info(section_header("Loading Existing Simulation", use_unicode))
+        logger.info(section("Loading Existing Simulation"))
         # Load existing simulation data
         if config.acquisition_type == 'DIA':
             existing_sim_handle = SyntheticExperimentDataHandleDIA(database_path=config.existing_path)
@@ -1610,7 +1706,7 @@ def main():
     # ----------------------------------------
     findings_result = None
     if config.from_findings and not resume_after:
-        logger.info(section_header("Loading Search Engine Findings", use_unicode))
+        logger.info(section("Loading Search Engine Findings"))
         rt_lower = acquisition_builder.frame_table['time'].min()
         rt_upper = acquisition_builder.frame_table['time'].max()
         findings_result = load_findings(
@@ -1638,7 +1734,7 @@ def main():
     protein_list, peptide_list = [], []
 
     if not config.from_existing and not config.from_findings and not resume_after:
-        logger.info(section_header("Processing FASTA Files", use_unicode))
+        logger.info(section("Processing FASTA Files"))
         for fasta_name, fasta_path in fastas.items():
             if not config.silent_mode:
                 logger.info("")
@@ -1668,6 +1764,7 @@ def main():
                 remove_degenerate_peptides=config.remove_degenerate_peptides,
             )
 
+            stages.lap('digest+rt_filter')
             # JOB 1: Simulate peptides
             if not config.silent_mode:
                 logger.info("  Creating peptides from proteins...")
@@ -1685,6 +1782,7 @@ def main():
                 proteome_mix=config.proteome_mix,
             )
 
+            stages.lap('sample_peptides')
             if config.proteome_mix:
                 # Scale by mixture factor
                 peptides_tmp['events'] *= mixture_factor
@@ -1725,7 +1823,7 @@ def main():
             logger.info(f"  Total peptides to simulate: {peptides.shape[0]}")
 
         # JOB 3: Simulate retention times
-        logger.info(section_header("Simulating Retention Times", use_unicode))
+        logger.info(section("Simulating Retention Times"))
         # Support both rt_model and deprecated koina_rt_model
         rt_model = config.rt_model or config.koina_rt_model
         if rt_model:
@@ -1739,7 +1837,7 @@ def main():
 
     # Simulate RT for from_findings when not provided in the input
     if config.from_findings and findings_result is not None and not findings_result.has_rt and not resume_after:
-        logger.info(section_header("Simulating Retention Times (not in findings)", use_unicode))
+        logger.info(section("Simulating Retention Times (not in findings)"))
         rt_model = config.rt_model or config.koina_rt_model
         if rt_model:
             logger.info(f"  Using RT model: {rt_model}")
@@ -1764,7 +1862,7 @@ def main():
 
     if resume_after != "ions":
         # JOB 4: Frame distributions
-        logger.info(section_header("Simulating Frame Distributions", use_unicode))
+        logger.info(section("Simulating Frame Distributions"))
         peptides = simulate_frame_distributions_emg(
             peptides=peptides,
             frames=acquisition_builder.frame_table,
@@ -1819,7 +1917,7 @@ def main():
         )
 
         if need_charge:
-            logger.info(section_header("Simulating Charge States", use_unicode))
+            logger.info(section("Simulating Charge States"))
             # JOB 5: Charge states
             ions = simulate_charge_states(
                 peptides=peptides,
@@ -1838,7 +1936,7 @@ def main():
                 ions = ions.drop_duplicates(subset=['sequence', 'charge'])
 
         if need_im:
-            logger.info(section_header("Simulating Ion Mobilities", use_unicode))
+            logger.info(section("Simulating Ion Mobilities"))
             # JOB 6: Ion mobilities
             if config.ccs_model:
                 logger.info(f"  Using CCS model: {config.ccs_model}")
@@ -1855,7 +1953,7 @@ def main():
 
         if not config.from_existing:
             # JOB 7: Precursor isotopic distributions (always needed)
-            logger.info(section_header("Simulating Precursor Isotope Patterns", use_unicode))
+            logger.info(section("Simulating Precursor Isotope Patterns"))
             ions = simulate_precursor_spectra_sequence(
                 ions=ions,
                 num_threads=num_threads,
@@ -1863,7 +1961,7 @@ def main():
             )
 
         # JOB 8: Scan distributions
-        logger.info(section_header("Simulating Scan Distributions", use_unicode))
+        logger.info(section("Simulating Scan Distributions"))
         ions = simulate_scan_distributions_with_variance(
             ions=ions,
             scans=acquisition_builder.scan_table,
@@ -1908,7 +2006,7 @@ def main():
             )
 
         if config.acquisition_type == 'DDA':
-            logger.info(section_header("Simulating DDA-PASEF Selection", use_unicode))
+            logger.info(section("Simulating DDA-PASEF Selection"))
             pasef_meta, precursors = simulate_dda_pasef_selection_scheme(
                 acquisition_builder=acquisition_builder,
                 verbose=not config.silent_mode,
@@ -1932,7 +2030,7 @@ def main():
 
     else:
         # Resuming from ions checkpoint — populate DB tables
-        logger.info(section_header("Restoring Database from Checkpoint", use_unicode))
+        logger.info(section("Restoring Database from Checkpoint"))
         acquisition_builder.synthetics_handle.create_table(table_name='proteins', table=proteins)
         acquisition_builder.synthetics_handle.create_table(table_name='peptides', table=peptides)
         acquisition_builder.synthetics_handle.create_table(table_name='ions', table=ions)
@@ -1943,7 +2041,7 @@ def main():
         logger.info(f"  Restored: proteins ({len(proteins)}), peptides ({len(peptides)}), ions ({len(ions)})")
 
     # JOB 9: Simulate fragment intensities
-    logger.info(section_header("Simulating Fragment Intensities", use_unicode))
+    logger.info(section("Simulating Fragment Intensities"))
     # Support both intensity_model and deprecated fragment_intensity_model
     intensity_model = config.intensity_model or config.fragment_intensity_model
     if intensity_model:
@@ -2122,8 +2220,8 @@ def main():
                 logger=logger,
             )
     else:
-        logger.info(section_header("Assembling Frames", use_unicode))
-        assemble_frames(
+        logger.info(section("Assembling Frames"))
+        assembly_timings = assemble_frames(
             acquisition_builder=acquisition_builder,
             frames=acquisition_builder.frame_table,
             batch_size=config.batch_size,
@@ -2149,6 +2247,8 @@ def main():
             quad_transmission_max_isotopes=config.quad_transmission_max_isotopes,
             superimpose_on_reference=config.superimpose_on_reference,
         )
+        if assembly_timings:
+            stages.annotate(**assembly_timings)
         # mzPROV provenance: sign the authored Bruker .d (self-disclosure that this is
         # TimSim-simulated data). Bruker-only — mzprov v0 doesn't canonicalize vendor
         # .raw, so the Thermo build-from-template branch above is intentionally skipped.
@@ -2175,7 +2275,7 @@ def main():
     # Optional: Generate preview video for visual inspection
     if config.generate_preview_video:
         if VIDEO_GENERATION_AVAILABLE:
-            logger.info(section_header("Generating Preview Video", use_unicode))
+            logger.info(section("Generating Preview Video"))
             data_path = os.path.join(save_path, name, f"{name}.d")
             video_path = os.path.join(save_path, f"{name}_preview.mp4")
             mode = 'dda' if config.acquisition_type == 'DDA' else 'dia'
@@ -2200,7 +2300,12 @@ def main():
     # Print completion banner and summary
     print(simulation_complete_banner(use_unicode))
 
+    stages.end()
     total_time = timer.format_duration(timer.total_elapsed())
+    logger.info("  Stage timings")
+    for line in stages.summary_lines():
+        logger.info(f"  {line}")
+    logger.info("")
     logger.info("  Simulation Summary")
     logger.info(f"  {'─' * 40}")
     logger.info(f"  Experiment:    {stats.experiment_name}")
@@ -2213,6 +2318,12 @@ def main():
     logger.info(f"  Total time:    {total_time}")
     logger.info(f"  Output:        {stats.output_path}")
     logger.info("")
+    try:
+        timings_path = Path(acquisition_builder.path) / 'timings.json'
+    except NameError:  # pragma: no cover - acquisition builder is always created in main()
+        timings_path = Path(save_path) / 'timings.json'
+    stages.dump(timings_path)
+    logger.info(f"  Stage timings written to {timings_path}")
 
 
 if __name__ == '__main__':
diff --git a/packages/imspy-simulation/src/imspy_simulation/utility.py b/packages/imspy-simulation/src/imspy_simulation/utility.py
index 1faee5340..55fa06c28 100644
--- a/packages/imspy-simulation/src/imspy_simulation/utility.py
+++ b/packages/imspy-simulation/src/imspy_simulation/utility.py
@@ -182,9 +182,20 @@ def read_acquisition_config(acquisition_name: str = 'dia') -> Dict[str, Any]:
 
 # Function to convert a list (or a pandas series) to a JSON string
 def python_list_to_json_string(lst, as_float=True, num_decimals: int = 4) -> str:
+    """Serialise a 1-D sequence as a JSON array (floats rounded to ``num_decimals``, or ints).
+
+    Vectorised: one ``np.round`` over the whole array instead of a NumPy scalar call per element
+    (this helper runs once per peptide/ion in four pipeline stages, so the per-element version
+    cost minutes on a 250k simulation). Output is identical to the previous implementation.
+    """
     if as_float:
-        return json.dumps([float(np.round(x, num_decimals)) for x in lst])
-    return json.dumps([int(x) for x in lst])
+        arr = np.asarray(lst)
+        if arr.dtype.kind != 'f':  # ints, bools, object lists of Python numbers
+            arr = arr.astype(np.float64)
+        # round in the input precision (float32 stays float32) and widen afterwards, exactly as
+        # the former per-element ``float(np.round(x, n))`` did, so the strings are byte-identical.
+        return json.dumps(np.round(arr, num_decimals).astype(np.float64).tolist())
+    return json.dumps(np.asarray(lst).astype(np.int64).tolist())
 
 
 # load peptides and ions
diff --git a/packages/imspy-simulation/tests/test_perf_serial_hotspots.py b/packages/imspy-simulation/tests/test_perf_serial_hotspots.py
new file mode 100644
index 000000000..1d25fca43
--- /dev/null
+++ b/packages/imspy-simulation/tests/test_perf_serial_hotspots.py
@@ -0,0 +1,147 @@
+"""Parity tests for the serial hot-spot fixes (packed-key (scan, tof) dedup, vectorised JSON
+helper, persistent tdf_bin handle, stage tracker).
+
+Each fix must reproduce the previous implementation byte-for-byte; the reference
+implementations are copied here verbatim from the pre-change code.
+"""
+import io
+import json
+import logging
+import os
+import tempfile
+
+import numpy as np
+import pytest
+
+from imspy_simulation.tdf import dedup_scan_tof
+from imspy_simulation.utility import python_list_to_json_string
+
+
+# --------------------------------------------------------------------------- references
+def _dedup_reference(scan, tof, intensity):
+    scan_tof = np.stack((scan.astype(np.uint32), tof.astype(np.uint32)), axis=1)
+    unique_pairs, inverse_indices = np.unique(scan_tof, axis=0, return_inverse=True)
+    summed_intensity = np.bincount(inverse_indices.ravel(), weights=intensity)
+    unique_scan = unique_pairs[:, 0]
+    unique_tof = unique_pairs[:, 1]
+    sort_idx = np.lexsort((unique_tof, unique_scan))
+    return unique_scan[sort_idx], unique_tof[sort_idx], summed_intensity[sort_idx].astype(np.uint32)
+
+
+def _json_reference(lst, as_float=True, num_decimals=4):
+    if as_float:
+        return json.dumps([float(np.round(x, num_decimals)) for x in lst])
+    return json.dumps([int(x) for x in lst])
+
+
+# --------------------------------------------------------------------------- dedup
+@pytest.mark.parametrize("seed,n,n_scans,n_tof", [(0, 10, 5, 20), (1, 5_000, 927, 400_000),
+                                                   (2, 60_000, 927, 400_000), (3, 3_000, 3, 10)])
+def test_dedup_scan_tof_matches_reference(seed, n, n_scans, n_tof):
+    rng = np.random.default_rng(seed)
+    scan = rng.integers(0, n_scans, n).astype(np.uint32)
+    tof = rng.integers(0, n_tof, n).astype(np.uint32)
+    intensity = rng.integers(1, 5_000, n).astype(np.uint32)
+    ref = _dedup_reference(scan, tof, intensity)
+    got = dedup_scan_tof(scan, tof, intensity)
+    for r, g in zip(ref, got):
+        np.testing.assert_array_equal(r, g)
+    assert got[0].dtype == np.uint32 and got[1].dtype == np.uint32 and got[2].dtype == np.uint32
+    # order contract: sorted by scan, then tof; no duplicate cells left
+    key = got[0].astype(np.uint64) << np.uint64(32) | got[1].astype(np.uint64)
+    assert np.all(np.diff(key) > 0)
+
+
+def test_dedup_scan_tof_sums_duplicates():
+    scan = np.array([2, 2, 1, 2], dtype=np.uint32)
+    tof = np.array([7, 7, 9, 3], dtype=np.uint32)
+    inten = np.array([10, 5, 1, 2], dtype=np.uint32)
+    s, t, i = dedup_scan_tof(scan, tof, inten)
+    assert s.tolist() == [1, 2, 2]
+    assert t.tolist() == [9, 3, 7]
+    assert i.tolist() == [1, 2, 15]
+
+
+def test_dedup_scan_tof_empty():
+    e = np.array([], dtype=np.uint32)
+    s, t, i = dedup_scan_tof(e, e, e)
+    assert len(s) == len(t) == len(i) == 0
+
+
+def test_dedup_scan_tof_large_tof_values_do_not_collide():
+    # tof close to 2**32 must stay separate from scan+1, tof=0
+    scan = np.array([0, 1], dtype=np.uint32)
+    tof = np.array([np.iinfo(np.uint32).max, 0], dtype=np.uint32)
+    inten = np.array([1, 2], dtype=np.uint32)
+    s, t, i = dedup_scan_tof(scan, tof, inten)
+    assert s.tolist() == [0, 1] and t.tolist() == [np.iinfo(np.uint32).max, 0] and i.tolist() == [1, 2]
+
+
+# --------------------------------------------------------------------------- json helper
+@pytest.mark.parametrize("lst", [
+    [], [0.0], [1.23456789, 2.5, 1e-7, 123456.123456], list(np.random.default_rng(0).random(300)),
+    np.random.default_rng(1).random(50).astype(np.float32), [1, 2, 3], np.array([0.99995, 0.00005]),
+])
+def test_json_float_matches_reference(lst):
+    assert python_list_to_json_string(lst) == _json_reference(lst)
+    assert python_list_to_json_string(lst, num_decimals=2) == _json_reference(lst, num_decimals=2)
+
+
+@pytest.mark.parametrize("lst", [[], [0], [1, 2, 3], np.array([5, 6, 7], dtype=np.int32),
+                                 np.array([1.0, 2.0]), list(range(1000))])
+def test_json_int_matches_reference(lst):
+    assert python_list_to_json_string(lst, as_float=False) == _json_reference(lst, as_float=False)
+
+
+def test_json_nan_roundtrip_like_reference():
+    lst = [float('nan'), 1.0]
+    assert python_list_to_json_string(lst) == _json_reference(lst)
+
+
+# --------------------------------------------------------------------------- persistent handle
+def test_tdf_writer_binary_handle_reopens_and_closes():
+    from imspy_simulation.tdf import TDFWriter
+    w = TDFWriter.__new__(TDFWriter)  # bypass __init__ (needs a Bruker handle)
+    with tempfile.TemporaryDirectory() as d:
+        w.binary_file = os.path.join(d, "analysis.tdf_bin")
+        w._bin_fh = None
+        fh = w._binary_handle()
+        fh.write(b"abc")
+        assert w._binary_handle() is fh
+        w.close_binary()
+        assert w._bin_fh is None
+        w.close_binary()  # idempotent
+        w._binary_handle().write(b"def")
+        w.close_binary()
+        with open(w.binary_file, "rb") as r:
+            assert r.read() == b"abcdef"
+
+
+# --------------------------------------------------------------------------- stage tracker
+def test_stage_tracker_records_and_dumps():
+    from imspy_simulation.timsim.simulator import StageTracker
+    logger = logging.getLogger("stage-tracker-test")
+    stream = io.StringIO()
+    handler = logging.StreamHandler(stream)
+    logger.addHandler(handler)
+    logger.setLevel(logging.INFO)
+    st = StageTracker(logger)
+    st.begin("A")
+    st.lap("first")
+    st.lap("second")
+    st.begin("B")
+    st.annotate(build_s=1.5)
+    with tempfile.TemporaryDirectory() as d:
+        path = os.path.join(d, "timings.json")
+        st.dump(path)
+        payload = json.load(open(path))
+    names = [r['stage'] for r in payload['stages']]
+    assert names == ["A", "B"]
+    a, b = payload['stages']
+    assert set(a['laps']) == {"first", "second"}
+    assert 'laps' not in b and b['build_s'] == 1.5
+    for r in (a, b):
+        for k in ('wall_s', 'cpu_s', 'avg_cores', 'peak_rss_gb'):
+            assert k in r
+    assert "stage done: A" in stream.getvalue() and "stage done: B" in stream.getvalue()
+    assert len(st.summary_lines()) >= 3

From 140639cd8135f02804de14bbacd17b28148e693f Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 10:19:19 +0200
Subject: [PATCH 02/14] perf(rustdf): run parallel reads in SDK mode via the
 SDK-free converter

The Bruker SDK's tims_index_to_mz is not safe to call concurrently on one
handle. rustdf handled that by dropping to a single thread whenever
use_bruker_sdk is set, which is the default: get_pasef_fragments and
get_preprocessed_pasef_fragments ignored num_threads entirely in SDK mode.

The formula converter evaluates the calibration curves the .d file itself
carries and is Sync. Measured against the SDK on a real reference over 73k
peaks: worst |dmz| 3.2e-11 Th, worst 6.7e-8 ppm, mobility within one ULP,
54% of peaks bit-identical. That is orders of magnitude below any search
tolerance, and below the 4e-5 Th run-to-run drift concurrent SDK calls
were producing.

So parallel regions now read through that converter instead of serialising:

- TimsLazyLoder / TimsInMemoryLoader gain get_frame_with_converter, taking
  the converter explicitly rather than using self.index_converter.
- TimsDataLoader gains sdk_free_converter() and raw_data_path().
- The two DDA PASEF paths and the DIA reference-noise overlay use it, and
  fall back to one thread only when the file carries no usable calibration
  tables.

Adds a seeded overlay_reference_noise on TimsDatasetDIA (rayon, per-frame
StdRng derived from seed and frame id) plus its connector and Python
bindings, and generate_random_sample_with_rng on TimsFrame.
---
 imspy_connector/src/py_dia.rs                 |  28 +++++
 mscore/src/timstof/frame.rs                   |   8 +-
 .../imspy-core/src/imspy_core/timstof/dia.py  |  44 +++++++-
 rustdf/src/data/dda.rs                        |  44 +++++---
 rustdf/src/data/dia.rs                        | 102 ++++++++++++++++--
 rustdf/src/data/handle.rs                     |  75 ++++++++++---
 6 files changed, 265 insertions(+), 36 deletions(-)

diff --git a/imspy_connector/src/py_dia.rs b/imspy_connector/src/py_dia.rs
index 1f517fbcd..287bea0be 100644
--- a/imspy_connector/src/py_dia.rs
+++ b/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 {
diff --git a/mscore/src/timstof/frame.rs b/mscore/src/timstof/frame.rs
index 31ae18779..317b1870c 100644
--- a/mscore/src/timstof/frame.rs
+++ b/mscore/src/timstof/frame.rs
@@ -582,9 +582,15 @@ impl TimsFrame {
     }
 
     pub fn generate_random_sample(&self, take_probability: f64) -> TimsFrame {
+        let mut rng = rand::thread_rng();
+        self.generate_random_sample_with_rng(take_probability, &mut rng)
+    }
+
+    /// Same as `generate_random_sample`, but drawing from a caller-supplied RNG so that the
+    /// selection is reproducible (used by the seeded, parallel reference-noise overlay).
+    pub fn generate_random_sample_with_rng<R: Rng>(&self, take_probability: f64, rng: &mut R) -> TimsFrame {
         assert!(take_probability >= 0.0 && take_probability <= 1.0);
 
-        let mut rng = rand::thread_rng();
         let mut scan = Vec::new();
         let mut mobility = Vec::new();
         let mut tof = Vec::new();
diff --git a/packages/imspy-core/src/imspy_core/timstof/dia.py b/packages/imspy-core/src/imspy_core/timstof/dia.py
index 40b79eda9..1b52d0b46 100644
--- a/packages/imspy-core/src/imspy_core/timstof/dia.py
+++ b/packages/imspy-core/src/imspy_core/timstof/dia.py
@@ -1,5 +1,5 @@
 import sqlite3
-from typing import List
+from typing import List, Optional
 
 from imspy_core.core.base import RustWrapperObject
 from imspy_core.timstof.data import TimsDataset
@@ -72,6 +72,48 @@ def sample_fragment_signal(self, num_frames: int, window_group: int, max_intensi
 
         return TimsFrame.from_py_ptr(self.__dataset.sample_fragment_signal(num_frames, window_group, max_intensity, take_probability))
 
+    def overlay_reference_noise(
+            self,
+            frames: List[TimsFrame],
+            window_groups: List[Optional[int]],
+            num_precursor_frames: int = 5,
+            num_fragment_frames: int = 5,
+            max_intensity_precursor: float = 30.0,
+            max_intensity_fragment: float = 30.0,
+            take_precursor: float = 0.2,
+            take_fragment: float = 0.2,
+            seed: int = 0,
+            num_threads: int = 4,
+    ) -> List[TimsFrame]:
+        """Add sampled reference-data noise to a batch of frames, in parallel and reproducibly.
+
+        Equivalent to calling ``sample_precursor_signal`` / ``sample_fragment_signal`` per frame and
+        adding the result with ``frame + noise``, but done in Rust over a thread pool with a
+        per-frame RNG stream derived from ``seed`` and the frame id.
+
+        Args:
+            frames: Simulated frames.
+            window_groups: Per frame, the DIA window group (fragment frame) or None (precursor frame).
+            num_precursor_frames / num_fragment_frames: Reference frames sampled per output frame.
+            max_intensity_precursor / max_intensity_fragment: Keep reference peaks up to this intensity.
+            take_precursor / take_fragment: Fraction of reference peaks to keep.
+            seed: Master seed; identical inputs and seed give identical output.
+            num_threads: Rayon threads.
+
+        Returns:
+            List[TimsFrame]: frames with noise added, same order as the input.
+        """
+        assert len(frames) == len(window_groups), "frames and window_groups must have the same length"
+        out = self.__dataset.overlay_reference_noise(
+            [f.get_py_ptr() for f in frames],
+            [None if g is None else int(g) for g in window_groups],
+            int(num_precursor_frames), int(num_fragment_frames),
+            float(max_intensity_precursor), float(max_intensity_fragment),
+            float(take_precursor), float(take_fragment),
+            int(seed) & 0xFFFFFFFFFFFFFFFF, int(num_threads),
+        )
+        return [TimsFrame.from_py_ptr(f) for f in out]
+
     def read_compressed_data_full(self) -> List[bytes]:
         """Read compressed data.
 
diff --git a/rustdf/src/data/dda.rs b/rustdf/src/data/dda.rs
index 8d9b212fa..5fd3ecde2 100644
--- a/rustdf/src/data/dda.rs
+++ b/rustdf/src/data/dda.rs
@@ -683,14 +683,24 @@ impl TimsDatasetDDA {
             None => pasef_info.iter().collect(),
         };
 
-        // Note: The Bruker SDK is NOT thread-safe, so we must use sequential iteration
-        // when the SDK is being used for index conversion.
-        let uses_bruker_sdk = self.loader.uses_bruker_sdk();
+        // The Bruker SDK's index conversion is not safe to call concurrently on one handle, but
+        // the SDK-free formula converter evaluates the very calibration curves the file carries
+        // (matched the SDK to 3e-11 Th / 1 ULP on the reference data) and *is* `Sync`. So in SDK
+        // mode we read through that instead of dropping to a single thread; only a file without
+        // usable calibration tables still forces the serial path.
+        let sdk_free = self.loader.sdk_free_converter();
+        let conv: Option<&(dyn IndexConverter + Sync)> =
+            sdk_free.as_ref().map(|c| c as &(dyn IndexConverter + Sync));
+        let read_frame = |frame_id: u32| match conv {
+            Some(c) => self.loader.get_frame_with_converter(frame_id, c),
+            None => self.loader.get_frame(frame_id),
+        };
+        let force_serial = self.loader.uses_bruker_sdk() && conv.is_none();
 
         // Helper closure to process a single PASEF fragment
         let process_fragment = |pasef_info: &PasefMsMsMeta| -> PASEFDDAFragment {
             // get the frame
-            let frame = self.loader.get_frame(pasef_info.frame_id as u32);
+            let frame = read_frame(pasef_info.frame_id as u32);
 
             // get five percent of the scan range
             let scan_margin = (pasef_info.scan_num_end - pasef_info.scan_num_begin) / 20;
@@ -718,8 +728,8 @@ impl TimsDatasetDDA {
             }
         };
 
-        if uses_bruker_sdk {
-            // Sequential processing when using Bruker SDK (not thread-safe)
+        if force_serial {
+            // No thread-safe converter available: fall back to one thread.
             filtered_pasef_info.iter().map(|info| process_fragment(info)).collect()
         } else {
             // Parallel processing when using simple index converter (thread-safe)
@@ -783,9 +793,19 @@ impl TimsDatasetDDA {
         }
 
         // Step 4: Build fragment data with aggregation (merge frames for same precursor)
-        // Note: The Bruker SDK is NOT thread-safe, so we must use sequential iteration
-        // when the SDK is being used for index conversion.
-        let uses_bruker_sdk = self.loader.uses_bruker_sdk();
+        // The Bruker SDK's index conversion is not safe to call concurrently on one handle, but
+        // the SDK-free formula converter evaluates the very calibration curves the file carries
+        // (matched the SDK to 3e-11 Th / 1 ULP on the reference data) and *is* `Sync`. So in SDK
+        // mode we read through that instead of dropping to a single thread; only a file without
+        // usable calibration tables still forces the serial path.
+        let sdk_free = self.loader.sdk_free_converter();
+        let conv: Option<&(dyn IndexConverter + Sync)> =
+            sdk_free.as_ref().map(|c| c as &(dyn IndexConverter + Sync));
+        let read_frame = |frame_id: u32| match conv {
+            Some(c) => self.loader.get_frame_with_converter(frame_id, c),
+            None => self.loader.get_frame(frame_id),
+        };
+        let force_serial = self.loader.uses_bruker_sdk() && conv.is_none();
 
         // Helper closure to process a single precursor
         let process_precursor = |(precursor_id, pasef_infos): (&i64, &Vec<&PasefMsMsMeta>)| -> Option<PASEFFragmentData> {
@@ -807,7 +827,7 @@ impl TimsDatasetDDA {
 
             for pasef_info in pasef_infos {
                 // Get the frame and filter by scan range
-                let frame = self.loader.get_frame(pasef_info.frame_id as u32);
+                let frame = read_frame(pasef_info.frame_id as u32);
 
                 // Get five percent of the scan range for margin
                 let scan_margin = (pasef_info.scan_num_end - pasef_info.scan_num_begin) / 20;
@@ -860,8 +880,8 @@ impl TimsDatasetDDA {
             })
         };
 
-        let fragment_data: Vec<PASEFFragmentData> = if uses_bruker_sdk {
-            // Sequential processing when using Bruker SDK (not thread-safe)
+        let fragment_data: Vec<PASEFFragmentData> = if force_serial {
+            // No thread-safe converter available: fall back to one thread.
             pasef_by_precursor
                 .iter()
                 .filter_map(process_precursor)
diff --git a/rustdf/src/data/dia.rs b/rustdf/src/data/dia.rs
index bd692e013..7fa1eeb55 100644
--- a/rustdf/src/data/dia.rs
+++ b/rustdf/src/data/dia.rs
@@ -634,7 +634,7 @@ impl TimsDatasetDIA {
         let mut sampled_frames: Vec<TimsFrame> =
             self.collect_typed_noise_samples(&pool, num_frames, max_intensity,
                                              take_probability, &mut rng,
-                                             |meta| meta.id as u32);
+                                             |meta| meta.id as u32, None);
         Self::accumulate_or_typed_empty(&mut sampled_frames, target_type)
     }
 
@@ -661,7 +661,7 @@ impl TimsDatasetDIA {
         let mut sampled_frames: Vec<TimsFrame> =
             self.collect_typed_noise_samples(&pool, num_frames, max_intensity,
                                              take_probability, &mut rng,
-                                             |&fid| fid);
+                                             |&fid| fid, None);
         Self::accumulate_or_typed_empty(&mut sampled_frames, target_type)
     }
 
@@ -676,6 +676,7 @@ impl TimsDatasetDIA {
         take_probability: f64,
         rng: &mut impl rand::Rng,
         key_fn: impl Fn(&T) -> u32,
+        conv: Option<&(dyn IndexConverter + Sync)>,
     ) -> Vec<TimsFrame> {
         let mut out: Vec<TimsFrame> = Vec::with_capacity(num_frames);
         if pool.is_empty() || num_frames == 0 {
@@ -687,12 +688,17 @@ impl TimsDatasetDIA {
         while out.len() < num_frames && attempts < max_attempts {
             attempts += 1;
             let Some(candidate_meta) = pool.choose(rng) else { break };
-            let candidate = self
-                .loader
-                .get_frame(key_fn(candidate_meta))
+            let frame_id = key_fn(candidate_meta);
+            // `conv` is the SDK-free converter handed in by parallel callers (see
+            // `overlay_reference_noise`); the serial samplers read with the loader's own.
+            let raw = match conv {
+                Some(c) => self.loader.get_frame_with_converter(frame_id, c),
+                None => self.loader.get_frame(frame_id),
+            };
+            let candidate = raw
                 .filter_ranged(0.0, 2000.0, 0, 1000, 0.0, 5.0, 1.0,
                                max_intensity, 0, i32::MAX)
-                .generate_random_sample(take_probability);
+                .generate_random_sample_with_rng(take_probability, rng);
             // Skip the loader's empty-frame Default-fallback so the
             // returned noise never carries MsType::Unknown into the
             // outer simulator add.
@@ -729,6 +735,90 @@ impl TimsDatasetDIA {
     }
 
     /// All DIA window_group IDs present in the file (sorted unique).
+    /// Reference-noise overlay for a whole batch, in parallel and reproducible.
+    ///
+    /// For every simulated frame, draws `num_*_frames` random reference frames of the same kind
+    /// (precursor frames for MS1, frames of the same DIA window group for MS2), keeps a
+    /// `take_*` fraction of their peaks below `max_intensity_*`, accumulates them and adds the
+    /// result onto the simulated frame with `TimsFrame + TimsFrame` — exactly the per-frame
+    /// `sample_precursor_signal` / `sample_fragment_signal` + `frame + noise` sequence the
+    /// Python job used to run serially, but over a rayon pool and with a per-frame RNG stream
+    /// derived from `seed` and the frame id, so the result does not depend on scheduling.
+    ///
+    /// `window_groups[i]` is `Some(group)` for a fragment frame and `None` for a precursor frame.
+    pub fn overlay_reference_noise(
+        &self,
+        frames: Vec<TimsFrame>,
+        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<TimsFrame> {
+        assert_eq!(frames.len(), window_groups.len(), "frames and window_groups must align");
+        let precursor_pool: Vec<u32> = self
+            .meta_data
+            .iter()
+            .filter(|x| x.ms_ms_type == 0)
+            .map(|m| m.id as u32)
+            .collect();
+        let mut fragment_pools: FxHashMap<u32, Vec<u32>> = FxHashMap::default();
+        for info in &self.dia_ms_ms_info {
+            fragment_pools.entry(info.window_group).or_default().push(info.frame_id);
+        }
+        // The Bruker SDK's `tims_index_to_mz` is not safe to call concurrently on one handle
+        // (same seed gave m/z that differed by ~1e-8 relative on ~8 % of peaks, run to run).
+        // Read the reference frames through the SDK-free formula converter instead; if the
+        // file carries no calibration tables, fall back to one thread so the result stays
+        // deterministic.
+        let sdk_free = self.loader.sdk_free_converter();
+        let threads = if self.loader.uses_bruker_sdk() && sdk_free.is_none() {
+            1
+        } else {
+            num_threads.max(1)
+        };
+        let conv: Option<&(dyn IndexConverter + Sync)> =
+            sdk_free.as_ref().map(|c| c as &(dyn IndexConverter + Sync));
+        let pool = rayon::ThreadPoolBuilder::new()
+            .num_threads(threads)
+            .build()
+            .unwrap();
+        pool.install(|| {
+            frames
+                .into_par_iter()
+                .zip(window_groups.into_par_iter())
+                .map(|(frame, wg)| {
+                    use rand::SeedableRng;
+                    let stream = (frame.frame_id as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
+                    let mut rng = rand::rngs::StdRng::seed_from_u64(seed ^ stream);
+                    let noise = match wg {
+                        Some(group) => {
+                            let ids: &[u32] = fragment_pools
+                                .get(&group)
+                                .map(|v| v.as_slice())
+                                .unwrap_or(&[]);
+                            let mut s = self.collect_typed_noise_samples(
+                                ids, num_fragment_frames, max_intensity_fragment,
+                                take_fragment, &mut rng, |&fid| fid, conv);
+                            Self::accumulate_or_typed_empty(&mut s, MsType::FragmentDia)
+                        }
+                        None => {
+                            let mut s = self.collect_typed_noise_samples(
+                                &precursor_pool, num_precursor_frames, max_intensity_precursor,
+                                take_precursor, &mut rng, |&fid| fid, conv);
+                            Self::accumulate_or_typed_empty(&mut s, MsType::Precursor)
+                        }
+                    };
+                    frame + noise
+                })
+                .collect()
+        })
+    }
+
     pub fn dia_window_groups(&self) -> Vec<u32> {
         let mut gs: Vec<u32> = self
             .dia_ms_ms_info
diff --git a/rustdf/src/data/handle.rs b/rustdf/src/data/handle.rs
index 25267ce5e..6b99fff93 100644
--- a/rustdf/src/data/handle.rs
+++ b/rustdf/src/data/handle.rs
@@ -666,8 +666,11 @@ pub struct TimsLazyLoder {
     pub index_converter: TimsIndexConverter,
 }
 
-impl TimsData for TimsLazyLoder {
-    fn get_frame(&self, frame_id: u32) -> TimsFrame {
+impl TimsLazyLoder {
+    /// Same as [`TimsData::get_frame`] but converts tof/scan with an explicit converter instead of
+    /// `self.index_converter`. Parallel readers use this to bypass the Bruker SDK, whose
+    /// `tims_index_to_mz` is not safe to call concurrently on one handle.
+    pub fn get_frame_with_converter(&self, frame_id: u32, conv: &dyn IndexConverter) -> TimsFrame {
         let frame_index = (frame_id - 1) as usize;
 
         // turns out, there can be empty frames in the data, check for that, if so, return an empty frame
@@ -783,10 +786,8 @@ impl TimsData for TimsLazyLoder {
                 let intensity_dbl = intensities_.iter().map(|&x| x as f64).collect::<Vec<f64>>();
                 let tof_i32 = tof_indices_.iter().map(|&x| x as i32).collect::<Vec<i32>>();
 
-                let mz = self.index_converter.tof_to_mz(frame_id, &tof_indices_);
-                let inv_mobility = self
-                    .index_converter
-                    .scan_to_inverse_mobility(frame_id, &scan);
+                let mz = conv.tof_to_mz(frame_id, &tof_indices_);
+                let inv_mobility = conv.scan_to_inverse_mobility(frame_id, &scan);
 
                 let ms_type_raw = self.raw_data_layout.frame_meta_data[frame_index].ms_ms_type;
                 let ms_type = match ms_type_raw {
@@ -823,10 +824,8 @@ impl TimsData for TimsLazyLoder {
                 let tof_i32 = tof.iter().map(|&x| x as i32).collect();
                 let scan = flatten_scan_values(&scan, true);
 
-                let mz = self.index_converter.tof_to_mz(frame_id, &tof);
-                let inv_mobility = self
-                    .index_converter
-                    .scan_to_inverse_mobility(frame_id, &scan);
+                let mz = conv.tof_to_mz(frame_id, &tof);
+                let inv_mobility = conv.scan_to_inverse_mobility(frame_id, &scan);
 
                 let ms_type_raw = self.raw_data_layout.frame_meta_data[frame_index].ms_ms_type;
 
@@ -856,6 +855,12 @@ impl TimsData for TimsLazyLoder {
             }
         }
     }
+}
+
+impl TimsData for TimsLazyLoder {
+    fn get_frame(&self, frame_id: u32) -> TimsFrame {
+        self.get_frame_with_converter(frame_id, &self.index_converter)
+    }
 
     fn get_raw_frame(&self, frame_id: u32) -> RawTimsFrame {
         let frame_index = (frame_id - 1) as usize;
@@ -956,8 +961,11 @@ pub struct TimsInMemoryLoader {
     compressed_data: Vec<u8>,
 }
 
-impl TimsData for TimsInMemoryLoader {
-    fn get_frame(&self, frame_id: u32) -> TimsFrame {
+impl TimsInMemoryLoader {
+    /// Same as [`TimsData::get_frame`] but converts tof/scan with an explicit converter instead of
+    /// `self.index_converter`. Parallel readers use this to bypass the Bruker SDK, whose
+    /// `tims_index_to_mz` is not safe to call concurrently on one handle.
+    pub fn get_frame_with_converter(&self, frame_id: u32, conv: &dyn IndexConverter) -> TimsFrame {
         let raw_frame = self.get_raw_frame(frame_id);
 
         let raw_frame = match raw_frame.ms_type {
@@ -973,10 +981,8 @@ impl TimsData for TimsInMemoryLoader {
         let tof_i32 = raw_frame.tof.iter().map(|&x| x as i32).collect();
         let scan = flatten_scan_values(&raw_frame.scan, true);
 
-        let mz = self.index_converter.tof_to_mz(frame_id, &raw_frame.tof);
-        let inverse_mobility = self
-            .index_converter
-            .scan_to_inverse_mobility(frame_id, &scan);
+        let mz = conv.tof_to_mz(frame_id, &raw_frame.tof);
+        let inverse_mobility = conv.scan_to_inverse_mobility(frame_id, &scan);
 
         let ims_frame = ImsFrame::new(
             raw_frame.retention_time,
@@ -993,6 +999,12 @@ impl TimsData for TimsInMemoryLoader {
             ims_frame,
         }
     }
+}
+
+impl TimsData for TimsInMemoryLoader {
+    fn get_frame(&self, frame_id: u32) -> TimsFrame {
+        self.get_frame_with_converter(frame_id, &self.index_converter)
+    }
 
     fn get_raw_frame(&self, frame_id: u32) -> RawTimsFrame {
         let frame_index = (frame_id - 1) as usize;
@@ -1407,6 +1419,37 @@ impl TimsDataLoader {
         }
     }
 
+    /// `get_frame` with an explicit index converter (see [`TimsLazyLoder::get_frame_with_converter`]).
+    pub fn get_frame_with_converter(&self, frame_id: u32, conv: &dyn IndexConverter) -> TimsFrame {
+        match self {
+            TimsDataLoader::InMemory(loader) => loader.get_frame_with_converter(frame_id, conv),
+            TimsDataLoader::Lazy(loader) => loader.get_frame_with_converter(frame_id, conv),
+        }
+    }
+
+    /// Path of the `.d` folder this loader reads from.
+    pub fn raw_data_path(&self) -> &str {
+        match self {
+            TimsDataLoader::InMemory(loader) => &loader.raw_data_layout.raw_data_path,
+            TimsDataLoader::Lazy(loader) => &loader.raw_data_layout.raw_data_path,
+        }
+    }
+
+    /// A converter that is safe to use from many threads at once, for readers that want to
+    /// run in parallel while the dataset itself was opened with the Bruker SDK.
+    ///
+    /// Returns `None` when the loader's own converter is already thread-safe (anything but
+    /// `BrukerLib`), or when the file carries no calibration tables the SDK-free
+    /// [`BrukerFormulaConverter`] could evaluate — callers must then serialise their reads.
+    pub fn sdk_free_converter(&self) -> Option<TimsIndexConverter> {
+        if !self.uses_bruker_sdk() {
+            return None;
+        }
+        BrukerFormulaConverter::from_d_folder(self.raw_data_path(), 1)
+            .ok()
+            .map(TimsIndexConverter::BrukerFormula)
+    }
+
     /// Check if the Bruker SDK is being used for index conversion.
     /// The Bruker SDK is NOT thread-safe, so parallel operations that call
     /// the index converter must be disabled when using the SDK.

From 81bb84c89f6e9a81e507b2a59240f282c8c1b252 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 10:19:19 +0200
Subject: [PATCH 03/14] perf(imspy-simulation): batch the sage digest, use the
 parallel noise overlay

proteins_to_peptides_batched builds one sage indexed database per 2000
proteins instead of one per protein. generate_indexed_database sets up
sage's parallel index build, so calling it 20000 times burnt ~3000 core-s
on a 120-core box for a stage that ran faster with 4 threads. Peptide sets
per protein are unchanged: sage attaches every protein a peptide occurs in.

The DIA noise job now calls the batched, seeded overlay_reference_noise
when the connector provides it, passing sample_seed, and keeps the serial
per-frame loop as the fallback.

Smoke A/B, 2000 peptides, 5692-frame DIA reference, 16 threads, identical
config: wall 156 s released -> 138 s -> 119 s; frame assembly 54 -> 37.6 ->
21.6 s, of which the noise overlay is 18.9 -> 2.5 s. Output reads back with
imspy; peak counts differ by the reseeded noise sampler.

Tests cover batched-vs-per-protein digest parity with and without decoys,
overlay determinism across 1/4/32 threads with the SDK enabled (the case
that was flaky before the rustdf fix), and SDK vs formula converter parity.
---
 .../timsim/jobs/add_noise_from_real_data.py   |  21 +++
 .../timsim/jobs/assemble_frames.py            |   5 +
 .../timsim/jobs/simulate_proteins.py          | 118 +++++++++++++--
 .../src/imspy_simulation/timsim/simulator.py  |   5 +-
 .../test_perf_batched_digest_and_overlay.py   | 135 ++++++++++++++++++
 5 files changed, 267 insertions(+), 17 deletions(-)
 create mode 100644 packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py

diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/add_noise_from_real_data.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/add_noise_from_real_data.py
index d45d38088..15ac04f15 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/add_noise_from_real_data.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/add_noise_from_real_data.py
@@ -44,6 +44,8 @@ def add_real_data_noise_to_frames(
         num_fragment_frames: int = 5,
         acquisition_mode: str = 'DIA',
         pasef_meta: Optional[pd.DataFrame] = None,
+        noise_seed: Optional[int] = None,
+        num_threads: int = 4,
 ) -> List[TimsFrame]:
     """Add noise to frame.
 
@@ -102,6 +104,25 @@ def add_real_data_noise_to_frames(
         d = acquisition_builder.frames_to_window_groups
         window_group_dict = dict(zip(d['frame'], d['window_group']))
         fragment_frames = set(window_group_dict.keys())
+        handle = acquisition_builder.tdf_writer.helper_handle
+
+        # Fast path (imspy-connector >= 0.4.3): sample + add for the whole batch in Rust, in parallel,
+        # with per-frame seeded RNG streams. Falls back to the serial per-frame loop below on older
+        # connectors, which behaves identically apart from being unseeded.
+        if noise_seed is not None and hasattr(handle, "overlay_reference_noise"):
+            groups = [int(window_group_dict[f.frame_id]) if f.frame_id in fragment_frames else None
+                      for f in frames]
+            return handle.overlay_reference_noise(
+                frames, groups,
+                num_precursor_frames=num_precursor_frames,
+                num_fragment_frames=num_fragment_frames,
+                max_intensity_precursor=intensity_max_precursor,
+                max_intensity_fragment=intensity_max_fragment,
+                take_precursor=precursor_sample_fraction,
+                take_fragment=fragment_sample_fraction,
+                seed=noise_seed,
+                num_threads=num_threads,
+            )
 
         for frame in frames:
             if frame.frame_id in fragment_frames:
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
index b1f1267b7..40ea386d6 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
@@ -42,6 +42,7 @@ def assemble_frames(
         quad_transmission_min_probability: float = 0.5,
         quad_transmission_max_isotopes: int = 10,
         superimpose_on_reference: bool = False,
+        noise_seed: Optional[int] = None,
 ) -> None:
     """Assemble frames from frame ids and write them to the database.
 
@@ -77,6 +78,8 @@ def assemble_frames(
             consider for transmission (default 10).
         superimpose_on_reference: If True, superimpose simulated signals on top
             of the full, unmodified reference frames (DIA only).
+        noise_seed: Master seed for the reference-noise sampler (DIA fast path); None keeps the
+            legacy unseeded per-frame sampling.
 
     Returns:
         None, writes frames to disk and metadata to database.
@@ -164,6 +167,8 @@ def assemble_frames(
                 num_fragment_frames=num_fragment_frames,
                 acquisition_mode=acquisition_builder.acquisition_mode.mode,
                 pasef_meta=pasef_meta,
+                noise_seed=noise_seed,
+                num_threads=num_threads,
             )
         t2 = time.perf_counter()
 
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
index 0776574ef..fc4d99046 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
@@ -152,6 +152,78 @@ def protein_to_peptides(fasta,
         return None
 
 
+def proteins_to_peptides_batched(
+        names,
+        sequences,
+        missed_cleavages: int = 2,
+        min_len: int = 7,
+        max_len: int = 30,
+        cleave_at: str = 'KR',
+        restrict: str = 'P',
+        generate_decoys: bool = False,
+        c_terminal: bool = True,
+        variable_mods: dict = {"M": ["[UNIMOD:35]"], "[": ["[UNIMOD:1]"]},
+        static_mods: dict = {"C": "[UNIMOD:4]"},
+        chunk_size: int = 2000,
+) -> list:
+    """Digest many proteins with ONE sage indexed database per chunk instead of one per protein.
+
+    Returns, for every input protein, the same ``set`` of UniMod peptide strings that
+    ``protein_to_peptides`` produces for that protein alone: sage attaches every protein a
+    peptide occurs in (``peptide.proteins``), so a shared peptide is credited to each of its
+    proteins exactly as the per-protein digest did. Proteins without a valid peptide get an
+    empty set (the per-protein version returned None; callers drop both).
+
+    Motivation: ``generate_indexed_database()`` sets up sage's parallel index build; calling it
+    20 000 times (once per protein) spent ~3 000 core-seconds on a 120-core machine for a
+    93 s stage that ran *faster* with 4 threads. One database per chunk is ~10x cheaper.
+    Headers are the positional index into ``names``; sage keeps the header token up to the first
+    whitespace, so protein names never leak into the mapping. Decoy proteins (``rev_`` prefix)
+    are credited to their target protein, matching the per-protein behaviour.
+    """
+    names = list(names)
+    sequences = list(sequences)
+    result = [set() for _ in names]
+    if not names:
+        return result
+    enzyme_builder = EnzymeBuilder(
+        missed_cleavages=missed_cleavages,
+        min_len=min_len,
+        max_len=max_len,
+        cleave_at=cleave_at,
+        restrict=restrict,
+        c_terminal=c_terminal,
+    )
+    for start in range(0, len(names), chunk_size):
+        stop = min(start + chunk_size, len(names))
+        fasta = "".join(
+            f">{i}\n" + "\n".join(sequences[i][j:j + 60] for j in range(0, len(sequences[i]), 60)) + "\n"
+            for i in range(start, stop)
+        )
+        sage_config = SageSearchConfiguration(
+            fasta=fasta,
+            enzyme_builder=enzyme_builder,
+            static_mods=static_mods,
+            variable_mods=variable_mods,
+            generate_decoys=generate_decoys,
+            bucket_size=int(np.power(2, 6)),
+        )
+        try:
+            with SuppressStderr():  # sage panics (harmlessly) on proteins that yield no peptide
+                indexed_db = sage_config.generate_indexed_database()
+        except Exception:
+            # whole chunk without a single valid peptide: leave the sets empty
+            continue
+        for k in range(indexed_db.num_peptides):
+            peptide = indexed_db[PeptideIx(idx=k)]
+            seq = peptide.to_unimod_sequence()
+            for prot in peptide.proteins:
+                if prot.startswith("rev_"):
+                    prot = prot[4:]
+                result[int(prot)].add(seq)
+    return result
+
+
 def get_tenzer_hokey():
     # Tenzer constants
     delta_exp = 0.005
@@ -246,21 +318,37 @@ def simulate_proteins(
         print("Skipping digestion, returning raw protein sequences as peptides.")
 
     # Generate peptides
-    sample["peptides"] = sample.apply(lambda f: protein_to_peptides(
-        generate_single_fasta(
-            f.name,
-            f.sequence,
-        ),
-        generate_decoys=generate_decoys,
-        variable_mods=variable_mods,
-        static_mods=static_mods,
-        cleave_at=cleave_at,
-        restrict=restrict,
-        missed_cleavages=missed_cleavages,
-        min_len=min_len,
-        max_len=max_len,
-        digest=digest,
-    ), axis=1)
+    if digest:
+        # one sage database per chunk of proteins instead of one per protein (see
+        # proteins_to_peptides_batched); same per-protein peptide sets as protein_to_peptides
+        peptide_sets = proteins_to_peptides_batched(
+            sample.index, sample.sequence,
+            generate_decoys=generate_decoys,
+            variable_mods=variable_mods,
+            static_mods=static_mods,
+            cleave_at=cleave_at,
+            restrict=restrict,
+            missed_cleavages=missed_cleavages,
+            min_len=min_len,
+            max_len=max_len,
+        )
+        sample["peptides"] = [s if len(s) > 0 else None for s in peptide_sets]
+    else:
+        sample["peptides"] = sample.apply(lambda f: protein_to_peptides(
+            generate_single_fasta(
+                f.name,
+                f.sequence,
+            ),
+            generate_decoys=generate_decoys,
+            variable_mods=variable_mods,
+            static_mods=static_mods,
+            cleave_at=cleave_at,
+            restrict=restrict,
+            missed_cleavages=missed_cleavages,
+            min_len=min_len,
+            max_len=max_len,
+            digest=digest,
+        ), axis=1)
 
     # Remove None values
     sample = sample[sample.peptides.notnull()]
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
index 408df9e91..55668e1b1 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
@@ -1764,7 +1764,7 @@ def section(title: str) -> str:
                 remove_degenerate_peptides=config.remove_degenerate_peptides,
             )
 
-            stages.lap('digest+rt_filter')
+            stages.lap('digest')
             # JOB 1: Simulate peptides
             if not config.silent_mode:
                 logger.info("  Creating peptides from proteins...")
@@ -1782,7 +1782,7 @@ def section(title: str) -> str:
                 proteome_mix=config.proteome_mix,
             )
 
-            stages.lap('sample_peptides')
+            stages.lap('sample_peptides+rt_filter')
             if config.proteome_mix:
                 # Scale by mixture factor
                 peptides_tmp['events'] *= mixture_factor
@@ -2246,6 +2246,7 @@ def section(title: str) -> str:
             quad_transmission_min_probability=config.quad_transmission_min_probability,
             quad_transmission_max_isotopes=config.quad_transmission_max_isotopes,
             superimpose_on_reference=config.superimpose_on_reference,
+            noise_seed=config.sample_seed,
         )
         if assembly_timings:
             stages.annotate(**assembly_timings)
diff --git a/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
new file mode 100644
index 000000000..f218e4950
--- /dev/null
+++ b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
@@ -0,0 +1,135 @@
+"""Parity tests: batched sage digest == per-protein digest; seeded parallel noise overlay is
+deterministic and equals the serial per-frame sampling in distribution."""
+import os
+import numpy as np
+import pandas as pd
+import pytest
+
+from imspy_simulation.timsim.jobs.simulate_proteins import (
+    protein_to_peptides, proteins_to_peptides_batched, generate_single_fasta, parse_fasta_to_dataframe,
+)
+
+FASTA = "/scratch/timsim-demo/timsim-necro/flow/hela_subset.fasta"
+REF_D = "/scratch/timsim-demo/TIMSIM-HeLa-RUSTW-SMOKE/TIMSIM-HeLa-RUSTW-SMOKE.d"
+
+VAR = {"M": ["[UNIMOD:35]"], "[": ["[UNIMOD:1]"]}
+STAT = {"C": "[UNIMOD:4]"}
+
+
+@pytest.mark.skipif(not os.path.exists(FASTA), reason="local HeLa subset fasta not available")
+@pytest.mark.parametrize("decoys", [False, True])
+def test_batched_digest_matches_per_protein(decoys):
+    tbl = parse_fasta_to_dataframe(FASTA).sample(n=400, random_state=7)
+    batched = proteins_to_peptides_batched(
+        tbl.index, tbl.sequence, generate_decoys=decoys, variable_mods=VAR, static_mods=STAT,
+        chunk_size=150,  # force several chunks
+    )
+    n_nonempty = 0
+    for (idx, row), got in zip(tbl.iterrows(), batched):
+        ref = protein_to_peptides(
+            generate_single_fasta(idx, row.sequence), generate_decoys=decoys,
+            variable_mods=VAR, static_mods=STAT, cleave_at='KR', restrict='P',
+            missed_cleavages=2, min_len=7, max_len=30, digest=True,
+        )
+        if ref is None:
+            assert got == set(), f"protein {idx}: per-protein digest empty, batched gave {len(got)}"
+        else:
+            assert got == ref, f"protein {idx}: {len(got ^ ref)} differing peptides"
+            n_nonempty += 1
+    assert n_nonempty > 300
+
+
+def test_batched_digest_empty_input():
+    assert proteins_to_peptides_batched([], []) == []
+
+
+@pytest.mark.skipif(not os.path.exists(REF_D), reason="local reference .d not available")
+def test_overlay_reference_noise_is_seeded_and_matches_serial_statistics():
+    from imspy_core.timstof import TimsDatasetDIA
+    ds = TimsDatasetDIA(REF_D, in_memory=False)
+    if not hasattr(ds, "overlay_reference_noise"):
+        pytest.skip("imspy-core without overlay_reference_noise")
+    info = ds.dia_ms_ms_info
+    wg_of = dict(zip(info.Frame, info.WindowGroup))
+    ids = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]
+    frames = [ds.get_tims_frame(i) for i in ids]
+    groups = [int(wg_of[i]) if i in wg_of else None for i in ids]
+    kw = dict(num_precursor_frames=5, num_fragment_frames=5, max_intensity_precursor=150000.0,
+              max_intensity_fragment=150000.0, take_precursor=0.2, take_fragment=0.2, num_threads=4)
+    a = ds.overlay_reference_noise(frames, groups, seed=41, **kw)
+    b = ds.overlay_reference_noise(frames, groups, seed=41, **kw)
+    c = ds.overlay_reference_noise(frames, groups, seed=42, **kw)
+    assert len(a) == len(frames)
+    for fa, fb, fc, f0 in zip(a, b, c, frames):
+        assert fa.frame_id == f0.frame_id
+        np.testing.assert_array_equal(fa.mz, fb.mz)
+        np.testing.assert_array_equal(fa.intensity, fb.intensity)
+        np.testing.assert_array_equal(fa.scan, fb.scan)
+        assert len(fa.mz) >= len(f0.mz)          # noise only adds peaks
+        assert not (len(fa.mz) == len(fc.mz) and np.array_equal(fa.mz, fc.mz))  # seed matters
+    # serial reference path (unseeded) must give the same added-peak count within sampling noise
+    added_par = np.array([len(x.mz) - len(f.mz) for x, f in zip(a, frames)], dtype=float)
+    added_ser = []
+    for f, g in zip(frames, groups):
+        if g is None:
+            noise = ds.sample_precursor_signal(5, 150000.0, 0.2)
+        else:
+            noise = ds.sample_fragment_signal(5, g, 150000.0, 0.2)
+        added_ser.append(len((f + noise).mz) - len(f.mz))
+    added_ser = np.array(added_ser, dtype=float)
+    assert added_par.sum() > 0 and added_ser.sum() > 0
+    ratio = added_par.sum() / added_ser.sum()
+    assert 0.5 < ratio < 2.0, f"parallel/serial added-peak ratio {ratio:.2f}"
+
+
+# ---------------------------------------------------------------------------
+# Bruker SDK thread-safety regression guards.
+#
+# `tims_index_to_mz` is not safe to call concurrently on one SDK handle: before the
+# fix, the same seed produced m/z differing by up to 4e-5 Th on ~8 % of peaks, in
+# 3/8 runs at 4 threads and 7/8 at 32. rustdf now reads through the SDK-free formula
+# converter inside parallel regions instead of serialising, so these must hold with
+# the SDK enabled (the default).
+# ---------------------------------------------------------------------------
+
+@pytest.mark.skipif(not os.path.exists(REF_D), reason="local reference .d not available")
+@pytest.mark.parametrize("threads", [1, 4, 32])
+def test_overlay_is_bit_identical_across_thread_counts_with_sdk(threads):
+    from imspy_core.timstof import TimsDatasetDIA
+    ds = TimsDatasetDIA(REF_D, in_memory=False, use_bruker_sdk=True)
+    info = ds.dia_ms_ms_info
+    wg_of = dict(zip(info.Frame, info.WindowGroup))
+    ids = list(range(1, 19))
+    frames = [ds.get_tims_frame(i) for i in ids]
+    groups = [int(wg_of[i]) if i in wg_of else None for i in ids]
+    kw = dict(num_precursor_frames=5, num_fragment_frames=5, max_intensity_precursor=150000.0,
+              max_intensity_fragment=150000.0, take_precursor=0.2, take_fragment=0.2)
+    ref = ds.overlay_reference_noise(frames, groups, seed=41, num_threads=1, **kw)
+    for rep in range(5):
+        got = ds.overlay_reference_noise(frames, groups, seed=41, num_threads=threads, **kw)
+        for a, b in zip(ref, got):
+            assert len(a.mz) == len(b.mz), f"rep {rep}: peak count differs on frame {a.frame_id}"
+            np.testing.assert_array_equal(a.mz, b.mz)
+            np.testing.assert_array_equal(a.intensity, b.intensity)
+
+
+@pytest.mark.skipif(not os.path.exists(REF_D), reason="local reference .d not available")
+def test_sdk_free_converter_matches_the_sdk():
+    """The general fix rests on this: the formula converter reproduces the SDK, so parallel
+    regions may use it instead of dropping to one thread."""
+    from imspy_core.timstof import TimsDatasetDIA
+    sdk = TimsDatasetDIA(REF_D, in_memory=False, use_bruker_sdk=True)
+    fml = TimsDatasetDIA(REF_D, in_memory=False, use_bruker_sdk=False)
+    checked = 0
+    for fid in [1, 2, 3, 5, 8, 13, 21, 34, 55, 89]:
+        a, b = sdk.get_tims_frame(fid), fml.get_tims_frame(fid)
+        if len(a.mz) == 0:
+            continue
+        assert len(a.mz) == len(b.mz)
+        np.testing.assert_array_equal(a.scan, b.scan)
+        np.testing.assert_array_equal(a.intensity, b.intensity)
+        # last-bit agreement, orders of magnitude below any search tolerance
+        assert np.abs(a.mz - b.mz).max() < 1e-8, "m/z drift between SDK and formula converter"
+        assert np.abs(a.mobility - b.mobility).max() < 1e-12
+        checked += 1
+    assert checked >= 5

From 4c934354de6b454f178a40687c2da36c9c328f41 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 10:44:21 +0200
Subject: [PATCH 04/14] release: imspy-connector 0.4.3, imspy-core 0.4.3,
 imspy-simulation 0.4.4

The round-2 work and the rustdf SDK fix change behaviour and performance, so
they need version numbers of their own: the overnight timing grid ran against
a different code state that also called itself 0.4.3, which would have made the
manuscript's version reference ambiguous.

imspy-core moves too because the Python wrapper for the batched overlay lives
there; without it the noise job silently falls back to the serial per-frame
loop. imspy-simulation now requires >=0.4.3 of both so that cannot happen.
---
 imspy_connector/Cargo.toml               | 2 +-
 imspy_connector/pyproject.toml           | 2 +-
 packages/imspy-core/pyproject.toml       | 4 ++--
 packages/imspy-simulation/pyproject.toml | 6 +++---
 4 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/imspy_connector/Cargo.toml b/imspy_connector/Cargo.toml
index a2ebf3716..5215b407b 100644
--- a/imspy_connector/Cargo.toml
+++ b/imspy_connector/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "imspy-connector"
-version = "0.4.2"
+version = "0.4.3"
 edition = "2021"
 
 [lib]
diff --git a/imspy_connector/pyproject.toml b/imspy_connector/pyproject.toml
index 2f8d3e26f..7d04331a3 100644
--- a/imspy_connector/pyproject.toml
+++ b/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.3"
 requires-python = ">=3.11"
 classifiers = [
     "Programming Language :: Rust",
diff --git a/packages/imspy-core/pyproject.toml b/packages/imspy-core/pyproject.toml
index 679faa7fe..d3a2535ae 100644
--- a/packages/imspy-core/pyproject.toml
+++ b/packages/imspy-core/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "imspy-core"
-version = "0.4.2"
+version = "0.4.3"
 description = "Core data structures and utilities for processing timsTOF ion mobility spectrometry data."
 authors = [
     { name = "theGreatHerrLebert", email = "davidteschner@googlemail.com" }
@@ -20,7 +20,7 @@ dependencies = [
     "pyarrow>=13.0",
     "mendeleev>=0.7.0",
     "toml>=0.10.2",
-    "imspy-connector>=0.4.2",
+    "imspy-connector>=0.4.3",
     "opentims-bruker-bridge>=1.1.0",
 ]
 
diff --git a/packages/imspy-simulation/pyproject.toml b/packages/imspy-simulation/pyproject.toml
index a9acc5dc1..a1fe909d0 100644
--- a/packages/imspy-simulation/pyproject.toml
+++ b/packages/imspy-simulation/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "imspy-simulation"
-version = "0.4.3"
+version = "0.4.4"
 description = "TimsTOF data simulation tools for proteomics."
 authors = [
     { name = "theGreatHerrLebert", email = "davidteschner@googlemail.com" }
@@ -10,8 +10,8 @@ license = "MIT"
 requires-python = ">=3.11,<3.14"
 
 dependencies = [
-    "imspy-connector>=0.4.2",
-    "imspy-core>=0.4.2",
+    "imspy-connector>=0.4.3",
+    "imspy-core>=0.4.3",
     "imspy-predictors[koina]>=0.5.2",  # Include Koina support for intensity prediction
     "pandas>=2.0",
     "numpy>=1.24",

From b03316c6e236033e66e8ab0dfc92b840e43d5356 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 17:06:53 +0200
Subject: [PATCH 05/14] perf(writer): build, dedup, compress frames in one
 parallel Rust call
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Profiling the TDF writer per step over 700 real frames (3.07M peaks) put 56 %
of its time in the Bruker SDK's mz->tof and 1/K0->scan calls, 31 % in the
interleave, and only 5 % in zstd — so the phase named "convert+compress+write"
was never compression-bound, and tuning the zstd level would have been wasted.

The conversions ran one frame at a time for the same reason the noise overlay
did: they went through the SDK, which is not safe to call concurrently. In this
direction the SDK-free converter is not merely close but exact — over 1,702,413
values, mz->tof and 1/K0->scan agree with the SDK bit for bit, because the
outputs are integers and the last-bit float drift never crosses a rounding
boundary.

TDFWriter.write_frames now hands a whole batch to rustdf, which converts,
dedups (scan, tof), interleaves and compresses across a rayon pool and returns
the payloads plus the statistics the Frames row needs. Only the append and the
running byte offset stay sequential. Passing numpy arrays straight through
matters: an earlier parallel interleave existed but took Python lists, and the
.tolist() marshalling alone was 21 % of writer time, which ate the gain.

Output is byte-identical to the per-frame path, 700/700 frames at 1, 4 and 16
threads, and a test pins that. Getting there caught a real bug: zstd::encode_all
omits the decompressed size from the frame header, and the simple API every
existing .d was written with then refuses to read the frame. Uses
zstd::bulk::compress, which writes it.

Writer throughput on 700 frames: 1.36 s per-frame Python, 0.32 s batched at one
thread, 0.056 s at 16 (24x). End-to-end smoke, 16 threads: the write phase falls
11.8 -> 2.5 s and frame assembly 21.6 -> 10.0 s.

Also removes the commented-out compress_frames/write_frames skeleton this
replaces; it was shelved as "needs debugging" and had omitted the dedup, which
would have corrupted output since mz->tof is not injective.
---
 imspy_connector/src/py_dataset.rs             |  67 ++++++++++-
 .../imspy-core/src/imspy_core/timstof/data.py |  36 ++++++
 .../src/imspy_simulation/tdf.py               | 108 ++++++++++++------
 .../timsim/jobs/assemble_frames.py            |   8 +-
 .../tests/test_batched_frame_writer.py        |  76 ++++++++++++
 rustdf/src/data/handle.rs                     |  17 +++
 rustdf/src/data/utility.rs                    |  85 +++++++++++++-
 7 files changed, 357 insertions(+), 40 deletions(-)
 create mode 100644 packages/imspy-simulation/tests/test_batched_frame_writer.py

diff --git a/imspy_connector/src/py_dataset.rs b/imspy_connector/src/py_dataset.rs
index 00e3c0b48..ccaded782 100644
--- a/imspy_connector/src/py_dataset.rs
+++ b/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,57 @@ 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(
+                "build_compressed_frames: the Bruker SDK is in use and this .d carries no \
+                 calibration tables for the SDK-free converter, so frames cannot be converted \
+                 in parallel; 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,
+            )
+        });
+
+        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();
diff --git a/packages/imspy-core/src/imspy_core/timstof/data.py b/packages/imspy-core/src/imspy_core/timstof/data.py
index 9f2170007..240ddf662 100644
--- a/packages/imspy-core/src/imspy_core/timstof/data.py
+++ b/packages/imspy-core/src/imspy_core/timstof/data.py
@@ -329,6 +329,42 @@ def tof_to_mz(self, frame_id: int, tof_values: NDArray[np.int32]) -> NDArray[np.
         """
         return self.__dataset.tof_to_mz(frame_id, tof_values)
 
+    def build_compressed_frames(
+            self,
+            frame_ids: List[int],
+            mz: List[NDArray[np.float64]],
+            mobility: List[NDArray[np.float64]],
+            intensity: List[NDArray[np.float64]],
+            max_scans: int,
+            compression_level: int = 0,
+            num_threads: int = 4,
+    ) -> List[tuple]:
+        """Run the per-frame TDF writer pipeline for a batch of frames in parallel.
+
+        Converts m/z to TOF and 1/K0 to scan, merges duplicate (scan, tof) cells, interleaves into
+        the Bruker layout and zstd-compresses — all in Rust over a thread pool, so the caller only
+        has to append the results in order.
+
+        Args:
+            frame_ids: Frame id to take the calibration from, per frame.
+            mz / mobility / intensity: Per-frame peak arrays, contiguous float64.
+            max_scans: Scan count of the output layout.
+            compression_level: zstd level; 0 is libzstd's default, matching the per-frame writer.
+            num_threads: Rayon threads.
+
+        Returns:
+            One ``(num_peaks, max_intensity, summed_intensity, compressed_bytes)`` per input frame,
+            in input order. The statistics are taken after the dedup, as the Frames table records them.
+
+        Raises:
+            RuntimeError: if the dataset uses the Bruker SDK and carries no calibration tables for
+                the SDK-free converter, since the SDK cannot be called from several threads at once.
+        """
+        return self.__dataset.build_compressed_frames(
+            [int(f) for f in frame_ids], mz, mobility, intensity,
+            int(max_scans), int(compression_level), int(num_threads),
+        )
+
     def mz_to_tof(self, frame_id: int, mz_values: NDArray[np.float64]) -> NDArray[np.int32]:
         """Convert m/z values to TOF values.
 
diff --git a/packages/imspy-simulation/src/imspy_simulation/tdf.py b/packages/imspy-simulation/src/imspy_simulation/tdf.py
index 7850d2943..31f56efc4 100644
--- a/packages/imspy-simulation/src/imspy_simulation/tdf.py
+++ b/packages/imspy-simulation/src/imspy_simulation/tdf.py
@@ -6,6 +6,7 @@
 import numpy as np
 
 from pathlib import Path
+from typing import List
 
 from numpy._typing import NDArray
 
@@ -351,6 +352,24 @@ def build_frame_meta_row(
             scan_mode: int,
             frame_start_pos: int,
             only_frame_one: bool = False
+    ):
+        """Thin wrapper: derive the Frames statistics from an intensity array, then delegate."""
+        return self._frame_meta_row(
+            frame, scan_mode, frame_start_pos, only_frame_one,
+            num_peaks=len(intensity),
+            max_intensity=int(np.max(intensity)) if len(intensity) > 0 else 0,
+            summed_intensity=int(np.sum(intensity)) if len(intensity) > 0 else 0,
+        )
+
+    def _frame_meta_row(
+            self,
+            frame: TimsFrame,
+            scan_mode: int,
+            frame_start_pos: int,
+            only_frame_one: bool,
+            num_peaks: int,
+            max_intensity: int,
+            summed_intensity: int,
     ):
         """Build a row for the frame meta data table from a TimsFrame object.
             Arguments:
@@ -378,10 +397,10 @@ def build_frame_meta_row(
         r.ScanMode = scan_mode
         r.MsMsType = frame.ms_type
         r.TimsId = frame_start_pos
-        r.MaxIntensity = int(np.max(intensity)) if len(intensity) > 0 else 0
-        r.SummedIntensities = int(np.sum(intensity)) if len(intensity) > 0 else 0
+        r.MaxIntensity = int(max_intensity)
+        r.SummedIntensities = int(summed_intensity)
         r.NumScans = self.helper_handle.num_scans
-        r.NumPeaks = len(intensity)
+        r.NumPeaks = int(num_peaks)
 
         return r
 
@@ -452,6 +471,59 @@ def write_frame(self, frame: TimsFrame, scan_mode: int, only_frame_one: bool = F
         bin_file.write(compressed_data)
         self.position = bin_file.tell()
 
+    def _conversion_frame_id(self, frame: TimsFrame, only_frame_one: bool) -> int:
+        """Frame id whose calibration is used for this frame, clamped to the reference's range."""
+        try:
+            max_index = self.helper_handle.meta_data.Id.max()
+        except AttributeError:
+            max_index = self.helper_handle.meta_data.frame_id.max()
+        if only_frame_one:
+            return 1
+        return int(max_index) if frame.frame_id > max_index else int(frame.frame_id)
+
+    def write_frames(self, frames: List[TimsFrame], scan_mode: int, only_frame_one: bool = False,
+                     num_threads: int = 4) -> None:
+        """Batched counterpart of :meth:`write_frame`.
+
+        Convert, dedup, interleave and compress every frame in one parallel Rust call, then append
+        the payloads in order — only the append and the running byte offset have to stay sequential.
+        Previously this ran per frame in Python, where the m/z -> TOF and 1/K0 -> scan conversions
+        went through the Bruker SDK and so could not be parallelised; measured at 56 % of writer time.
+
+        Falls back to the per-frame path if the parallel builder is unavailable (old connector, or a
+        reference that carries no calibration tables for the SDK-free converter).
+        """
+        if len(frames) == 0:
+            return
+
+        builder = getattr(self.helper_handle, "build_compressed_frames", None)
+        if builder is not None:
+            conv_ids = [self._conversion_frame_id(f, only_frame_one) for f in frames]
+            mz = [np.ascontiguousarray(f.mz, dtype=np.float64) for f in frames]
+            mobility = [np.ascontiguousarray(f.mobility, dtype=np.float64) for f in frames]
+            intensity = [np.ascontiguousarray(f.intensity, dtype=np.float64) for f in frames]
+            try:
+                built = builder(conv_ids, mz, mobility, intensity,
+                                int(self.helper_handle.num_scans), 0, int(num_threads))
+            except Exception as e:
+                warnings.warn(f"batched frame writer unavailable ({e}); falling back to per-frame writing")
+                built = None
+            if built is not None:
+                bin_file = self._binary_handle()
+                num_scans_bytes = int(self.helper_handle.num_scans).to_bytes(4, "little", signed=False)
+                for frame, (num_peaks, max_i, sum_i, data) in zip(frames, built):
+                    self.frame_meta_data.append(
+                        self._frame_meta_row(frame, scan_mode, self.position, only_frame_one,
+                                             num_peaks, max_i, sum_i))
+                    bin_file.write((len(data) + 8).to_bytes(4, "little", signed=False))
+                    bin_file.write(num_scans_bytes)
+                    bin_file.write(data)
+                    self.position = bin_file.tell()
+                return
+
+        for frame in frames:
+            self.write_frame(frame, scan_mode, only_frame_one)
+
     def _binary_handle(self):
         """Return the persistent append handle for analysis.tdf_bin, opening it on first use."""
         if self._bin_fh is None or self._bin_fh.closed:
@@ -671,33 +743,3 @@ def write_dia_ms_ms_windows(self, dia_ms_ms_windows: pd.DataFrame) -> None:
 
         self._create_table(self.conn, out, "DiaFrameMsMsWindows")
 
-        # TODO: these methods needs to be debugged
-        """
-        def compress_frames(self, frames: List[TimsFrame], only_frame_one: bool = False, num_threads: int = 4) -> List[bytes]:
-            # same as compress_frame but for multiple frames
-            tofs, scans, intensities = [], [], []
-            for frame in frames:
-                i = 1 if only_frame_one else frame.frame_id
-                tofs.append(self.mz_to_tof(i, frame.mz).astype(np.uint32))
-                scans.append(self.inv_mobility_to_scan(i, frame.mobility).astype(np.uint32))
-                intensities.append(frame.intensity.astype(np.uint32))
-
-            real_data = ims.get_data_for_compression_par(tofs, scans, intensities, self.helper_handle.num_scans, num_threads)
-            return [zstd.ZSTD_compress(bytes(data), 1) for data in real_data]
-
-        def write_frames(self, frames: List[TimsFrame], scan_mode: int, only_frame_one: bool = False, num_threads: int = 4) -> None:
-
-            compressed_data = self.compress_frames(frames, only_frame_one, num_threads=num_threads)
-
-            for i, data in enumerate(compressed_data):
-
-                self.frame_meta_data.append(self.build_frame_meta_row(frames[i], scan_mode, self.position, only_frame_one))
-
-                with open(self.binary_file, "ab") as bin_file:
-                    bin_file.write(
-                        (len(data) + 8).to_bytes(4, "little", signed=False)
-                    )
-                    bin_file.write(int(self.helper_handle.num_scans).to_bytes(4, "little", signed=False))
-                    bin_file.write(data)
-                    self.position = bin_file.tell()
-        """
\ No newline at end of file
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
index 40ea386d6..35b68c9ee 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
@@ -172,11 +172,9 @@ def assemble_frames(
             )
         t2 = time.perf_counter()
 
-        for frame in built_frames:
-            if acquisition_builder.acquisition_mode.mode == 'DDA':
-                acquisition_builder.tdf_writer.write_frame(frame, scan_mode=8)
-            else:
-                acquisition_builder.tdf_writer.write_frame(frame, scan_mode=9)
+        scan_mode = 8 if acquisition_builder.acquisition_mode.mode == 'DDA' else 9
+        acquisition_builder.tdf_writer.write_frames(
+            built_frames, scan_mode=scan_mode, num_threads=num_threads)
         t3 = time.perf_counter()
 
         t_build += t1 - t0
diff --git a/packages/imspy-simulation/tests/test_batched_frame_writer.py b/packages/imspy-simulation/tests/test_batched_frame_writer.py
new file mode 100644
index 000000000..7f1bdb59c
--- /dev/null
+++ b/packages/imspy-simulation/tests/test_batched_frame_writer.py
@@ -0,0 +1,76 @@
+"""The batched writer must produce exactly what the per-frame writer produced.
+
+`TDFWriter.write_frames` moves conversion, dedup, interleave and zstd into one parallel Rust call.
+The output is a vendor file format, so "equivalent" is not good enough: the compressed payloads must
+be byte-identical to the per-frame path, and the Frames statistics must match.
+"""
+import os
+
+import numpy as np
+import pytest
+import zstd
+
+from imspy_simulation.tdf import dedup_scan_tof
+from imspy_simulation.utility import get_compressible_data
+
+REF_D = "/scratch/timsim-demo/TIMSIM-HeLa-RUSTW-SMOKE/TIMSIM-HeLa-RUSTW-SMOKE.d"
+pytestmark = pytest.mark.skipif(not os.path.exists(REF_D), reason="local reference .d not available")
+
+
+def _reference(ds, frames, num_scans):
+    """Exactly what TDFWriter.compress_frame did, frame by frame, through the SDK."""
+    out = []
+    for f in frames:
+        tof = np.array(ds.mz_to_tof(f.frame_id, f.mz)).astype(np.uint32)
+        scan = np.array(ds.inverse_mobility_to_scan(f.frame_id, f.mobility)).astype(np.uint32)
+        intensity = f.intensity.astype(np.uint32)
+        scan, tof, intensity = dedup_scan_tof(scan, tof, intensity)
+        blob = zstd.ZSTD_compress(bytes(get_compressible_data(tof, scan, intensity, num_scans)), 0)
+        out.append((len(intensity),
+                    int(intensity.max()) if len(intensity) else 0,
+                    int(intensity.sum()), blob))
+    return out
+
+
+def _load(n=120):
+    from imspy_core.timstof import TimsDatasetDIA
+    ds = TimsDatasetDIA(REF_D, in_memory=False, use_bruker_sdk=True)
+    frames = [f for f in (ds.get_tims_frame(i) for i in range(1, n + 1)) if len(f.mz) > 0]
+    assert frames, "reference produced no non-empty frames"
+    return ds, frames, int(ds.num_scans)
+
+
+@pytest.mark.parametrize("num_threads", [1, 4, 16])
+def test_batched_writer_is_byte_identical_to_per_frame(num_threads):
+    ds, frames, num_scans = _load()
+    if not hasattr(ds, "build_compressed_frames"):
+        pytest.skip("connector without build_compressed_frames")
+    ref = _reference(ds, frames, num_scans)
+    got = ds.build_compressed_frames(
+        [f.frame_id for f in frames],
+        [np.ascontiguousarray(f.mz, dtype=np.float64) for f in frames],
+        [np.ascontiguousarray(f.mobility, dtype=np.float64) for f in frames],
+        [np.ascontiguousarray(f.intensity, dtype=np.float64) for f in frames],
+        num_scans, 0, num_threads,
+    )
+    assert len(got) == len(ref)
+    for i, ((n1, m1, s1, b1), (n2, m2, s2, b2)) in enumerate(zip(ref, got)):
+        assert (n1, m1, s1) == (n2, m2, s2), f"frame {i}: Frames statistics differ"
+        # A zstd frame written without the decompressed size in its header cannot be read back by
+        # the simple API every existing .d was written with, so the header shape matters too.
+        assert b1 == bytes(b2), f"frame {i}: compressed payload differs"
+
+
+def test_batched_writer_output_still_decompresses_with_the_simple_api():
+    ds, frames, num_scans = _load(n=30)
+    if not hasattr(ds, "build_compressed_frames"):
+        pytest.skip("connector without build_compressed_frames")
+    got = ds.build_compressed_frames(
+        [f.frame_id for f in frames],
+        [np.ascontiguousarray(f.mz, dtype=np.float64) for f in frames],
+        [np.ascontiguousarray(f.mobility, dtype=np.float64) for f in frames],
+        [np.ascontiguousarray(f.intensity, dtype=np.float64) for f in frames],
+        num_scans, 0, 4,
+    )
+    for _, _, _, blob in got:
+        assert len(zstd.ZSTD_uncompress(bytes(blob))) > 0
diff --git a/rustdf/src/data/handle.rs b/rustdf/src/data/handle.rs
index 6b99fff93..dd77566a1 100644
--- a/rustdf/src/data/handle.rs
+++ b/rustdf/src/data/handle.rs
@@ -1435,6 +1435,23 @@ impl TimsDataLoader {
         }
     }
 
+    /// The loader's own index converter, if it is safe to use from several threads — i.e.
+    /// everything except the Bruker SDK wrapper, whose `tims_index_to_mz` must not be called
+    /// concurrently on one handle.
+    pub fn sync_index_converter(&self) -> Option<&(dyn IndexConverter + Sync)> {
+        let converter = match self {
+            TimsDataLoader::InMemory(loader) => &loader.index_converter,
+            TimsDataLoader::Lazy(loader) => &loader.index_converter,
+        };
+        match converter {
+            TimsIndexConverter::BrukerLib(_) => None,
+            TimsIndexConverter::Simple(c) => Some(c),
+            TimsIndexConverter::Calibrated(c) => Some(c),
+            TimsIndexConverter::Lookup(c) => Some(c),
+            TimsIndexConverter::BrukerFormula(c) => Some(c),
+        }
+    }
+
     /// A converter that is safe to use from many threads at once, for readers that want to
     /// run in parallel while the dataset itself was opened with the Bruker SDK.
     ///
diff --git a/rustdf/src/data/utility.rs b/rustdf/src/data/utility.rs
index 18f14e960..39cee6166 100644
--- a/rustdf/src/data/utility.rs
+++ b/rustdf/src/data/utility.rs
@@ -51,7 +51,7 @@ pub fn zstd_compress(decompressed_data: &[u8], compression_level: i32) -> io::Re
 /// frame data straight into the encoder, producing negative/garbage TOF deltas
 /// that vendor readers (e.g. DiaNN) reject. Mirroring the Python preprocessing
 /// here keeps the two writers byte-for-byte identical.
-fn sort_dedup_scan_tof(
+pub fn sort_dedup_scan_tof(
     scans: &[u32],
     tofs: &[u32],
     intensities: &[u32],
@@ -375,3 +375,86 @@ pub fn merge_ranges(mut ranges: Vec<(usize, usize)>) -> Vec<(usize, usize)> {
     out.push(cur);
     out
 }
+
+/// One frame's worth of writer output: the statistics the TDF `Frames` row needs, plus the
+/// zstd-compressed payload ready to append to `analysis.tdf_bin`.
+pub struct CompressedFrame {
+    pub num_peaks: u32,
+    pub max_intensity: u32,
+    pub summed_intensity: u64,
+    pub data: Vec<u8>,
+}
+
+/// The whole per-frame TDF writer pipeline — m/z → TOF and 1/K0 → scan conversion, `(scan, tof)`
+/// dedup, interleave into the Bruker layout, zstd — for a batch of frames, over a rayon pool.
+///
+/// This used to run one frame at a time in Python because the conversion went through the Bruker
+/// SDK, which is not safe to call concurrently on one handle. Given a `Sync` converter (the SDK-free
+/// `BrukerFormulaConverter` reproduces the SDK's integer output exactly in both of these directions)
+/// the whole pipeline is per-frame independent, so only the final append has to stay ordered.
+///
+/// `mz`, `mobility` and `intensity` are parallel slices, one entry per frame, matching `frame_ids`.
+/// Results come back in input order.
+pub fn build_compressed_frames(
+    converter: &(dyn crate::data::handle::IndexConverter + Sync),
+    frame_ids: &[u32],
+    mz: &[&[f64]],
+    mobility: &[&[f64]],
+    intensity: &[&[f64]],
+    max_scans: u32,
+    compression_level: i32,
+    num_threads: usize,
+) -> Vec<CompressedFrame> {
+    assert_eq!(frame_ids.len(), mz.len(), "frame_ids and mz must align");
+    assert_eq!(frame_ids.len(), mobility.len(), "frame_ids and mobility must align");
+    assert_eq!(frame_ids.len(), intensity.len(), "frame_ids and intensity must align");
+
+    let pool = rayon::ThreadPoolBuilder::new()
+        .num_threads(num_threads.max(1))
+        .build()
+        .unwrap();
+
+    pool.install(|| {
+        (0..frame_ids.len())
+            .into_par_iter()
+            .map(|i| {
+                let frame_id = frame_ids[i];
+                let mz_values: Vec<f64> = mz[i].to_vec();
+                let im_values: Vec<f64> = mobility[i].to_vec();
+
+                let tof = converter.mz_to_tof(frame_id, &mz_values);
+                let scan = converter.inverse_mobility_to_scan(frame_id, &im_values);
+                let raw_intensity: Vec<u32> = intensity[i].iter().map(|&x| x as u32).collect();
+
+                // m/z → TOF is not injective, so several peaks can land on one cell; sum them and
+                // sort by (scan, tof), which is the order the encoder's delta coding needs.
+                let (scan, tof, intensity) = sort_dedup_scan_tof(&scan, &tof, &raw_intensity);
+
+                let num_peaks = intensity.len() as u32;
+                let max_intensity = intensity.iter().copied().max().unwrap_or(0);
+                let summed_intensity: u64 = intensity.iter().map(|&x| x as u64).sum();
+
+                let mut tof_delta = tof.clone();
+                modify_tofs(&mut tof_delta, &scan);
+                let peak_cnts = get_peak_cnts(max_scans, &scan);
+                let interleaved: Vec<u32> = tof_delta
+                    .iter()
+                    .zip(intensity.iter())
+                    .flat_map(|(t, i)| [*t, *i])
+                    .collect();
+                let real_data = get_realdata(&peak_cnts, &interleaved);
+
+                CompressedFrame {
+                    num_peaks,
+                    max_intensity,
+                    summed_intensity,
+                    // `zstd::bulk::compress` (ZSTD_compress) writes the decompressed size into the
+                    // frame header; `encode_all` does not, and readers that use the simple API —
+                    // including the Python `zstd` module every existing .d was written with —
+                    // refuse a frame without it. Keep the header shape the format already has.
+                    data: zstd::bulk::compress(real_data.as_slice(), compression_level).unwrap(),
+                }
+            })
+            .collect()
+    })
+}

From d9be27bb803ebee0d9198bc1ec252c42179f5e12 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 17:23:51 +0200
Subject: [PATCH 06/14] fix(rustdf): only substitute the SDK-free converter
 when it is provably exact
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The exactness claim in the previous two commits rested
on a single fixture. It does not generalise, and the code's own doc comment said
so: BrukerFormulaConverter is bit-exact for MzCalibration ModelType 1, but
ModelType 2 omits the proprietary C8..C14 fine correction and lands ~2 ppm away.

Measured on the thunder blank, which is ModelType 2 and is a reference file the
timing grid uses: feeding it arbitrary m/z — what the writer actually receives —
11.5 % of values land in a neighbouring TOF bin. That changes the (scan, tof)
dedup and so the bytes written. The earlier test missed it twice over: its
fixture is ModelType 1, and it fed back m/z the SDK itself had produced, which
sit at bin centres and so are the case least likely to cross a rounding
boundary. Mobility is exact under both models; only the m/z axis is affected.

sdk_free_converter() now returns the converter only when is_sdk_exact() holds:
every frame resolves to a real calibration row — no silent fallback to the
default for NULL T1/T2, a missing Frames row, or an unsupported row — and every
m/z row it uses is ModelType 1. Otherwise callers stay on the SDK, serially.
Verified on real data: the ModelType-2 thunder blank is refused, the ModelType-1
dia blank accepted.

Also hardened here:
- build_compressed_frames validates that each frame's mz/mobility/intensity have
  equal length. A ragged triplet used to index out of bounds inside Rust; it now
  returns an error, as does a thread-pool or zstd failure instead of unwrapping.
- The Python fallback in write_frames no longer swallows every exception. Only
  the explicit UNAVAILABLE capability signal falls back to the per-frame path;
  anything else propagates rather than quietly writing a different file.

Tests: 193 Python passed / 15 skipped, 136 Rust passed.
---
 imspy_connector/src/py_dataset.rs             | 21 ++++++----
 .../src/imspy_simulation/tdf.py               | 12 +++++-
 .../tests/test_batched_frame_writer.py        | 26 ++++++++++++
 rustdf/src/data/handle.rs                     | 27 ++++++++++++
 rustdf/src/data/utility.rs                    | 41 ++++++++++++-------
 5 files changed, 101 insertions(+), 26 deletions(-)

diff --git a/imspy_connector/src/py_dataset.rs b/imspy_connector/src/py_dataset.rs
index ccaded782..c40135486 100644
--- a/imspy_connector/src/py_dataset.rs
+++ b/imspy_connector/src/py_dataset.rs
@@ -161,9 +161,10 @@ impl PyTimsDataset {
         };
         let Some(converter) = converter else {
             return Err(PyRuntimeError::new_err(
-                "build_compressed_frames: the Bruker SDK is in use and this .d carries no \
-                 calibration tables for the SDK-free converter, so frames cannot be converted \
-                 in parallel; use the per-frame writer path",
+                "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.",
             ));
         };
 
@@ -171,12 +172,14 @@ impl PyTimsDataset {
         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,
-            )
-        });
+        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()
diff --git a/packages/imspy-simulation/src/imspy_simulation/tdf.py b/packages/imspy-simulation/src/imspy_simulation/tdf.py
index 31f56efc4..da61db8c3 100644
--- a/packages/imspy-simulation/src/imspy_simulation/tdf.py
+++ b/packages/imspy-simulation/src/imspy_simulation/tdf.py
@@ -505,8 +505,16 @@ def write_frames(self, frames: List[TimsFrame], scan_mode: int, only_frame_one:
             try:
                 built = builder(conv_ids, mz, mobility, intensity,
                                 int(self.helper_handle.num_scans), 0, int(num_threads))
-            except Exception as e:
-                warnings.warn(f"batched frame writer unavailable ({e}); falling back to per-frame writing")
+            except RuntimeError as e:
+                # Only the explicit capability signal falls back. Anything else — a ragged array,
+                # a compression or thread-pool failure — is a real error, and silently writing the
+                # file through the other path would hide it behind a subtly different output.
+                if "UNAVAILABLE" not in str(e):
+                    raise
+                warnings.warn(
+                    f"batched frame writer unavailable, using the per-frame path: {e}",
+                    RuntimeWarning,
+                )
                 built = None
             if built is not None:
                 bin_file = self._binary_handle()
diff --git a/packages/imspy-simulation/tests/test_batched_frame_writer.py b/packages/imspy-simulation/tests/test_batched_frame_writer.py
index 7f1bdb59c..8c94c259e 100644
--- a/packages/imspy-simulation/tests/test_batched_frame_writer.py
+++ b/packages/imspy-simulation/tests/test_batched_frame_writer.py
@@ -74,3 +74,29 @@ def test_batched_writer_output_still_decompresses_with_the_simple_api():
     )
     for _, _, _, blob in got:
         assert len(zstd.ZSTD_uncompress(bytes(blob))) > 0
+
+
+def test_ragged_arrays_raise_instead_of_panicking():
+    """A ragged triplet used to index out of bounds inside Rust; it must surface as an error."""
+    ds, frames, num_scans = _load(n=10)
+    if not hasattr(ds, "build_compressed_frames"):
+        pytest.skip("connector without build_compressed_frames")
+    mz = [np.ascontiguousarray(f.mz, dtype=np.float64) for f in frames]
+    mobility = [np.ascontiguousarray(f.mobility, dtype=np.float64) for f in frames]
+    intensity = [np.ascontiguousarray(f.intensity, dtype=np.float64) for f in frames]
+    mobility[0] = mobility[0][:-1]          # one value short
+    with pytest.raises(RuntimeError, match="must have the same length"):
+        ds.build_compressed_frames([f.frame_id for f in frames], mz, mobility, intensity,
+                                   num_scans, 0, 4)
+
+
+def test_mismatched_outer_lengths_raise():
+    ds, frames, num_scans = _load(n=10)
+    if not hasattr(ds, "build_compressed_frames"):
+        pytest.skip("connector without build_compressed_frames")
+    mz = [np.ascontiguousarray(f.mz, dtype=np.float64) for f in frames]
+    mobility = [np.ascontiguousarray(f.mobility, dtype=np.float64) for f in frames]
+    intensity = [np.ascontiguousarray(f.intensity, dtype=np.float64) for f in frames]
+    with pytest.raises(RuntimeError, match="same length"):
+        ds.build_compressed_frames([f.frame_id for f in frames][:-1], mz, mobility, intensity,
+                                   num_scans, 0, 4)
diff --git a/rustdf/src/data/handle.rs b/rustdf/src/data/handle.rs
index dd77566a1..8c9661b00 100644
--- a/rustdf/src/data/handle.rs
+++ b/rustdf/src/data/handle.rs
@@ -528,6 +528,30 @@ impl BrukerFormulaConverter {
         })
     }
 
+    /// True only when substituting this converter for the Bruker SDK cannot change a single
+    /// output value, which is the only condition under which it is safe to swap in.
+    ///
+    /// Two things have to hold for **every** frame in the run:
+    ///   * its `MzCalibration` row is `ModelType` 1 — the only m/z model this converter reproduces
+    ///     bit-exactly. ModelType 2 shares the C0/C1/C2 curve but omits the proprietary C8..C14
+    ///     fine correction, landing ~2 ppm away; measured on a real ModelType-2 run that moves
+    ///     **11.5 % of arbitrary m/z values into a neighbouring TOF bin**, which would change the
+    ///     `(scan, tof)` dedup and therefore the bytes written to `analysis.tdf_bin`.
+    ///   * it resolves to a real calibration row. Frames with NULL `T1`/`T2`, a missing `Frames`
+    ///     row, or a link to an unsupported calibration row silently fall back to the default
+    ///     calibration in `mz_for`/`im_for`, and the SDK does no such thing.
+    ///
+    /// Callers that cannot use a non-exact converter must fall back to the SDK, serially.
+    pub fn is_sdk_exact(&self) -> bool {
+        if self.frames.is_empty() {
+            return false;
+        }
+        self.frames.values().all(|f| {
+            self.mz_rows.get(&f.mz_id).map(|r| r.model_type == 1).unwrap_or(false)
+                && self.im_rows.contains_key(&f.tims_id)
+        })
+    }
+
     /// Spread (max - min, in degC) of `Frames.T1` and `Frames.T2` across the run.
     ///
     /// A near-zero spread means per-frame and fixed-frame calibration agree; a
@@ -1464,6 +1488,9 @@ impl TimsDataLoader {
         }
         BrukerFormulaConverter::from_d_folder(self.raw_data_path(), 1)
             .ok()
+            // Only substitute when the result is indistinguishable from the SDK's; otherwise the
+            // caller has to stay on the SDK and serialise. See `is_sdk_exact`.
+            .filter(|c| c.is_sdk_exact())
             .map(TimsIndexConverter::BrukerFormula)
     }
 
diff --git a/rustdf/src/data/utility.rs b/rustdf/src/data/utility.rs
index 39cee6166..cb9e93a54 100644
--- a/rustdf/src/data/utility.rs
+++ b/rustdf/src/data/utility.rs
@@ -404,15 +404,28 @@ pub fn build_compressed_frames(
     max_scans: u32,
     compression_level: i32,
     num_threads: usize,
-) -> Vec<CompressedFrame> {
-    assert_eq!(frame_ids.len(), mz.len(), "frame_ids and mz must align");
-    assert_eq!(frame_ids.len(), mobility.len(), "frame_ids and mobility must align");
-    assert_eq!(frame_ids.len(), intensity.len(), "frame_ids and intensity must align");
+) -> Result<Vec<CompressedFrame>, String> {
+    if frame_ids.len() != mz.len() || frame_ids.len() != mobility.len() || frame_ids.len() != intensity.len() {
+        return Err(format!(
+            "frame_ids ({}), mz ({}), mobility ({}) and intensity ({}) must have the same length",
+            frame_ids.len(), mz.len(), mobility.len(), intensity.len()
+        ));
+    }
+    // Per frame the three arrays index each other in `sort_dedup_scan_tof`, so a ragged triplet
+    // would be an out-of-bounds panic in Rust rather than an error the caller can handle.
+    for i in 0..frame_ids.len() {
+        if mz[i].len() != mobility[i].len() || mz[i].len() != intensity[i].len() {
+            return Err(format!(
+                "frame {} (id {}): mz ({}), mobility ({}) and intensity ({}) must have the same length",
+                i, frame_ids[i], mz[i].len(), mobility[i].len(), intensity[i].len()
+            ));
+        }
+    }
 
     let pool = rayon::ThreadPoolBuilder::new()
         .num_threads(num_threads.max(1))
         .build()
-        .unwrap();
+        .map_err(|e| format!("could not build the writer thread pool: {e}"))?;
 
     pool.install(|| {
         (0..frame_ids.len())
@@ -444,16 +457,14 @@ pub fn build_compressed_frames(
                     .collect();
                 let real_data = get_realdata(&peak_cnts, &interleaved);
 
-                CompressedFrame {
-                    num_peaks,
-                    max_intensity,
-                    summed_intensity,
-                    // `zstd::bulk::compress` (ZSTD_compress) writes the decompressed size into the
-                    // frame header; `encode_all` does not, and readers that use the simple API —
-                    // including the Python `zstd` module every existing .d was written with —
-                    // refuse a frame without it. Keep the header shape the format already has.
-                    data: zstd::bulk::compress(real_data.as_slice(), compression_level).unwrap(),
-                }
+                // `zstd::bulk::compress` (ZSTD_compress) writes the decompressed size into the
+                // frame header; `encode_all` does not, and readers that use the simple API —
+                // including the Python `zstd` module every existing .d was written with — refuse
+                // a frame without it. Keep the header shape the format already has.
+                let data = zstd::bulk::compress(real_data.as_slice(), compression_level)
+                    .map_err(|e| format!("frame id {frame_id}: zstd compression failed: {e}"))?;
+
+                Ok(CompressedFrame { num_peaks, max_intensity, summed_intensity, data })
             })
             .collect()
     })

From 9a5037c376604b6caadbaca986b6e04a9391651d Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 17:39:14 +0200
Subject: [PATCH 07/14] release: imspy-connector 0.4.4, imspy-core 0.4.4,
 imspy-simulation 0.4.5

The batched writer and the converter exactness gate change both behaviour and
output timing, so they need their own version numbers: the 0.4.4 timing grid
already ran against a different code state, and the manuscript cites a version.
---
 imspy_connector/Cargo.toml               | 2 +-
 imspy_connector/pyproject.toml           | 2 +-
 packages/imspy-core/pyproject.toml       | 4 ++--
 packages/imspy-simulation/pyproject.toml | 6 +++---
 4 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/imspy_connector/Cargo.toml b/imspy_connector/Cargo.toml
index 5215b407b..c0a8fe288 100644
--- a/imspy_connector/Cargo.toml
+++ b/imspy_connector/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "imspy-connector"
-version = "0.4.3"
+version = "0.4.4"
 edition = "2021"
 
 [lib]
diff --git a/imspy_connector/pyproject.toml b/imspy_connector/pyproject.toml
index 7d04331a3..77b599b40 100644
--- a/imspy_connector/pyproject.toml
+++ b/imspy_connector/pyproject.toml
@@ -7,7 +7,7 @@ name = "imspy_connector"
 dependencies = [
     "opentims-bruker-bridge>=1.1.0",
 ]
-version = "0.4.3"
+version = "0.4.4"
 requires-python = ">=3.11"
 classifiers = [
     "Programming Language :: Rust",
diff --git a/packages/imspy-core/pyproject.toml b/packages/imspy-core/pyproject.toml
index d3a2535ae..f59fb98b6 100644
--- a/packages/imspy-core/pyproject.toml
+++ b/packages/imspy-core/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "imspy-core"
-version = "0.4.3"
+version = "0.4.4"
 description = "Core data structures and utilities for processing timsTOF ion mobility spectrometry data."
 authors = [
     { name = "theGreatHerrLebert", email = "davidteschner@googlemail.com" }
@@ -20,7 +20,7 @@ dependencies = [
     "pyarrow>=13.0",
     "mendeleev>=0.7.0",
     "toml>=0.10.2",
-    "imspy-connector>=0.4.3",
+    "imspy-connector>=0.4.4",
     "opentims-bruker-bridge>=1.1.0",
 ]
 
diff --git a/packages/imspy-simulation/pyproject.toml b/packages/imspy-simulation/pyproject.toml
index a1fe909d0..79c1359f8 100644
--- a/packages/imspy-simulation/pyproject.toml
+++ b/packages/imspy-simulation/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "imspy-simulation"
-version = "0.4.4"
+version = "0.4.5"
 description = "TimsTOF data simulation tools for proteomics."
 authors = [
     { name = "theGreatHerrLebert", email = "davidteschner@googlemail.com" }
@@ -10,8 +10,8 @@ license = "MIT"
 requires-python = ">=3.11,<3.14"
 
 dependencies = [
-    "imspy-connector>=0.4.3",
-    "imspy-core>=0.4.3",
+    "imspy-connector>=0.4.4",
+    "imspy-core>=0.4.4",
     "imspy-predictors[koina]>=0.5.2",  # Include Koina support for intensity prediction
     "pandas>=2.0",
     "numpy>=1.24",

From 45c41b26b418b35059e65fda3d7f6a733505c07d Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 17:46:35 +0200
Subject: [PATCH 08/14] test: compare whole .d folders between the per-frame
 and batched writers
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The existing test compared compressed payloads in isolation. This writes the
same frames through both paths and compares what lands on disk: the Frames
table including the TimsId byte offsets, which depend on every preceding
frame's compressed length, and the decompressed content of every frame.

Running it on 600 frames turned up the one real nuance. The files are
equivalent but not byte-identical: 4 bytes in 8.6 MB differ, in 2 frames of
600, entirely inside the zstd stream. Those blobs decompress to identical
payloads, the compressed lengths match, and the Frames table is identical in
all 18 columns. The cause is that the Python zstd module and the Rust zstd
crate link different libzstd builds and can choose a different, equally valid
encoding at the same cost. Reading both files back frame by frame — through
the Bruker SDK and through the SDK-free reader — gives 600/600 identical
frames over 2.55M peaks in mz, tof, scan and intensity.

So the contract the test pins is the decompressed data and the table, not the
byte encoding.
---
 .../tests/test_batched_frame_writer.py        | 51 +++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/packages/imspy-simulation/tests/test_batched_frame_writer.py b/packages/imspy-simulation/tests/test_batched_frame_writer.py
index 8c94c259e..6be15c31b 100644
--- a/packages/imspy-simulation/tests/test_batched_frame_writer.py
+++ b/packages/imspy-simulation/tests/test_batched_frame_writer.py
@@ -100,3 +100,54 @@ def test_mismatched_outer_lengths_raise():
     with pytest.raises(RuntimeError, match="same length"):
         ds.build_compressed_frames([f.frame_id for f in frames][:-1], mz, mobility, intensity,
                                    num_scans, 0, 4)
+
+
+def test_two_writers_produce_equivalent_d_folders(tmp_path):
+    """End-to-end: write the same frames through both paths and compare the .d folders.
+
+    The per-frame check above compares payloads in isolation; this one compares what actually
+    lands on disk — the Frames table (including the TimsId byte offsets, which depend on every
+    preceding frame's compressed length) and the decompressed content of every frame.
+
+    Note the files are equivalent, not byte-identical: the Python `zstd` module and the Rust
+    `zstd` crate link different libzstd builds and can pick a different, equally valid encoding
+    of the same input. Observed at 4 bytes in 8.6 MB, in 2 frames of 600, all decompressing to
+    identical payloads. The contract is the decompressed data and the table, not the encoding.
+    """
+    import sqlite3
+    import pandas as pd
+    from imspy_simulation.tdf import TDFWriter
+
+    ds, frames, _ = _load(n=150)
+    if not hasattr(ds, "build_compressed_frames"):
+        pytest.skip("connector without build_compressed_frames")
+
+    def write(name, fn):
+        w = TDFWriter(helper_handle=ds, path=str(tmp_path), exp_name=name)
+        fn(w)
+        w.close_binary()
+        w.write_frame_meta_data()
+        return tmp_path / name
+
+    p_old = write("old.d", lambda w: [w.write_frame(f, scan_mode=9) for f in frames])
+    p_new = write("new.d", lambda w: w.write_frames(frames, scan_mode=9, num_threads=8))
+
+    old = pd.read_sql("select * from Frames order by Id", sqlite3.connect(str(p_old / "analysis.tdf")))
+    new = pd.read_sql("select * from Frames order by Id", sqlite3.connect(str(p_new / "analysis.tdf")))
+    assert old.equals(new), "Frames table differs between the per-frame and batched writers"
+
+    b_old = (p_old / "analysis.tdf_bin").read_bytes()
+    b_new = (p_new / "analysis.tdf_bin").read_bytes()
+    assert len(b_old) == len(b_new)
+
+    # Walk the frames by their recorded offsets and compare decompressed payloads.
+    for _, row in old.iterrows():
+        off = int(row.TimsId)
+        for buf in (b_old, b_new):
+            assert int.from_bytes(buf[off:off + 4], "little") > 0
+        n_old = int.from_bytes(b_old[off:off + 4], "little")
+        n_new = int.from_bytes(b_new[off:off + 4], "little")
+        assert n_old == n_new, f"frame {int(row.Id)}: compressed length differs"
+        d_old = zstd.ZSTD_uncompress(bytes(b_old[off + 8:off + n_old]))
+        d_new = zstd.ZSTD_uncompress(bytes(b_new[off + 8:off + n_new]))
+        assert d_old == d_new, f"frame {int(row.Id)}: decompressed payload differs"

From b88ab37286699f87f81f933de3875680d76e40b1 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 18:10:33 +0200
Subject: [PATCH 09/14] fix(timsim): make a run reproducible from its config
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Two runs of an identical config with the same sample_seed shared only 11 of
~1690 peptides. Same configuration, effectively unrelated datasets. The 1.4 %
peak-count difference that surfaced this was the mild symptom; the datasets
underneath had almost nothing in common.

Two independent causes, both upstream of anything seeded:

1. The global NumPy RNG was never seeded. sample_seed reached exactly two
   places — the pandas down-sample in simulator.py and (since 0.4.4) the
   reference-noise overlay. Everything deciding *what* gets simulated used bare
   np.random.*: peptide selection and abundance in simulate_peptides, the EMG
   retention-time shape, ion-mobility and scan variance, phosphosite choice,
   the decoy draw in digest_fasta. 16 call sites in jobs/ alone.

2. Both digest paths returned a `set` of peptide strings, and the sampler draws
   from it with np.random.choice. Set iteration order depends on Python's
   per-process string hash seed, so the draw differed run to run even with the
   RNG seeded — verified directly: same seed, three processes, three different
   selections. Both digests now return sorted lists.

With both fixes every intermediate table is bit-identical between two runs:
proteins, peptides, ions, fragment_ions, frames, scans. That is the ground
truth the benchmarking compares against, so benchmarks are now reproducible
from a config.

Known remaining gap: the written .d still differs by ~0.002 % of peaks (570 of
29 M, and 37 counts of summed intensity). That is the m/z jitter — mscore's
add_mz_noise_uniform/normal use rand::thread_rng at ~28 call sites across
mscore and rustdf/src/sim. A shifted m/z occasionally crosses a TOF bin
boundary and changes the dedup. Fixing it means threading a seeded per-spectrum
RNG through those call sites, as generate_random_sample_with_rng already does
for the overlay.

Note this changes which peptides a given config selects: datasets generated
before this cannot be reproduced by re-running their config, because they were
never reproducible to begin with.
---
 .../timsim/jobs/simulate_proteins.py          | 15 +++++++----
 .../src/imspy_simulation/timsim/simulator.py  | 16 ++++++++++++
 .../test_perf_batched_digest_and_overlay.py   | 25 ++++++++++++++++---
 3 files changed, 48 insertions(+), 8 deletions(-)

diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
index fc4d99046..0e16ef952 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_proteins.py
@@ -145,7 +145,10 @@ def protein_to_peptides(fasta,
             peptide = indexed_db[idx]
             peptide_set.add(peptide.to_unimod_sequence())
 
-        return peptide_set
+        # Sorted, not a set: `sample_peptides_from_proteins` draws from this with
+        # `np.random.choice`, and a set's iteration order depends on Python's per-process
+        # string hash seed. That made the draw differ between runs even with the RNG seeded.
+        return sorted(peptide_set)
 
     except Exception as e:
         # Proteins with no valid peptides after digestion will raise an exception
@@ -168,11 +171,11 @@ def proteins_to_peptides_batched(
 ) -> list:
     """Digest many proteins with ONE sage indexed database per chunk instead of one per protein.
 
-    Returns, for every input protein, the same ``set`` of UniMod peptide strings that
-    ``protein_to_peptides`` produces for that protein alone: sage attaches every protein a
+    Returns, for every input protein, the same UniMod peptide strings that
+    ``protein_to_peptides`` produces for that protein alone, as a sorted list: sage attaches every protein a
     peptide occurs in (``peptide.proteins``), so a shared peptide is credited to each of its
     proteins exactly as the per-protein digest did. Proteins without a valid peptide get an
-    empty set (the per-protein version returned None; callers drop both).
+    empty list (the per-protein version returned None; callers drop both).
 
     Motivation: ``generate_indexed_database()`` sets up sage's parallel index build; calling it
     20 000 times (once per protein) spent ~3 000 core-seconds on a 120-core machine for a
@@ -221,7 +224,9 @@ def proteins_to_peptides_batched(
                 if prot.startswith("rev_"):
                     prot = prot[4:]
                 result[int(prot)].add(seq)
-    return result
+    # Sorted lists rather than sets: see `protein_to_peptides`. The downstream sampler draws
+    # from these, and set iteration order is not stable across processes.
+    return [sorted(s) for s in result]
 
 
 def get_tenzer_hokey():
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
index 55668e1b1..2b3aaad1c 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
@@ -16,6 +16,7 @@
 from typing import Dict, Optional
 
 import toml
+import numpy as np
 import pandas as pd
 from tabulate import tabulate
 
@@ -1444,6 +1445,21 @@ def section(title: str) -> str:
     for key, value in config_dict.items():
         logger.debug(f"  {key}: {value}")
 
+    # Seed the global NumPy RNG for the whole run.
+    #
+    # Until 0.4.6 `sample_seed` reached exactly two places: the pandas down-sample below and the
+    # reference-noise overlay. Everything that decides *what* gets simulated — which peptides are
+    # drawn per protein and their abundances (simulate_peptides), the EMG retention-time shape,
+    # ion-mobility and scan variance, phosphosite choice, the decoy draw in digest_fasta — used
+    # bare `np.random.*` against the unseeded global RNG. Two runs of an identical config therefore
+    # shared only 11 of ~1690 peptides: same configuration, effectively unrelated datasets.
+    #
+    # Seeding once here makes a run reproducible from its config, which is what `sample_seed`
+    # always claimed to do. Note this changes which peptides a given config selects compared with
+    # releases before 0.4.6, so previously generated datasets are not reproduced by re-running
+    # their config — they were never reproducible in the first place.
+    np.random.seed(config.sample_seed & 0xFFFFFFFF)
+
     # Prepare paths
     save_path = check_path(config.save_path)
     # An Astral run has no Bruker reference .d (it builds from a .raw template).
diff --git a/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
index f218e4950..b600a203f 100644
--- a/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
+++ b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
@@ -31,10 +31,12 @@ def test_batched_digest_matches_per_protein(decoys):
             variable_mods=VAR, static_mods=STAT, cleave_at='KR', restrict='P',
             missed_cleavages=2, min_len=7, max_len=30, digest=True,
         )
-        if ref is None:
-            assert got == set(), f"protein {idx}: per-protein digest empty, batched gave {len(got)}"
+        if not ref:
+            assert got == [], f"protein {idx}: per-protein digest empty, batched gave {len(got)}"
         else:
-            assert got == ref, f"protein {idx}: {len(got ^ ref)} differing peptides"
+            assert got == sorted(ref), f"protein {idx}: differing peptides"
+            # Sorted, so the draw downstream cannot depend on set iteration order.
+            assert got == sorted(got), f"protein {idx}: batched digest is not sorted"
             n_nonempty += 1
     assert n_nonempty > 300
 
@@ -43,6 +45,23 @@ def test_batched_digest_empty_input():
     assert proteins_to_peptides_batched([], []) == []
 
 
+def test_digest_output_is_order_stable_not_a_set():
+    """Set iteration order varies with Python's per-process hash seed; the sampler draws from
+    these collections, so an unordered container makes a seeded run irreproducible."""
+    tbl = parse_fasta_to_dataframe(FASTA).sample(n=40, random_state=3)
+    batched = proteins_to_peptides_batched(tbl.index, tbl.sequence, variable_mods=VAR, static_mods=STAT)
+    assert any(len(p) > 1 for p in batched), "fixture produced nothing to order"
+    for peptides in batched:
+        assert isinstance(peptides, list), "digest must return an ordered container"
+        assert peptides == sorted(peptides)
+    for (idx, row) in tbl.iterrows():
+        got = protein_to_peptides(generate_single_fasta(idx, row.sequence), generate_decoys=False,
+                                  variable_mods=VAR, static_mods=STAT, cleave_at='KR', restrict='P',
+                                  missed_cleavages=2, min_len=7, max_len=30, digest=True)
+        if got:
+            assert isinstance(got, list) and got == sorted(got)
+
+
 @pytest.mark.skipif(not os.path.exists(REF_D), reason="local reference .d not available")
 def test_overlay_reference_noise_is_seeded_and_matches_serial_statistics():
     from imspy_core.timstof import TimsDatasetDIA

From b41806bb91d35855212111d5d2257d4beaaaf50c Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 18:49:02 +0200
Subject: [PATCH 10/14] fix(mscore,rustdf): seed the m/z jitter so a run is
 byte-reproducible
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The last source of run-to-run variation. add_mz_noise_uniform/normal drew from
rand::thread_rng, so which thread built a frame decided its jitter. Not
cosmetic: a shifted m/z occasionally crosses a TOF bin boundary, which changes
the (scan, tof) dedup and so the peaks written to analysis.tdf_bin. Two runs of
an identical config differed by ~0.002 % of peaks even with the ground truth
already reproducible.

mscore gains `_with_rng` variants of both, and `add_mz_noise` is now generic
over the RNG rather than pinned to ThreadRng. mscore::simulation::noise_rng
derives a stream from a master seed plus the stable identifiers already in
scope, so the draw does not depend on thread count or scheduling.

The four sim builders carry a `noise_seed` and create one RNG per frame, keyed
by (seed, frame_id): within a frame the spectra are built in a deterministic
order, and frames are independent. 24 call sites across dia, dda, precursor and
lazy_builder. `noise_seed` defaults to 0, which keeps the old thread_rng
behaviour for API users who never set it; TimSim sets it from the run's
sample_seed via create_frame_builder.

Two runs of an identical config now produce a byte-identical analysis.tdf_bin —
same sha256 over 93 MB — with the Frames table and every ground-truth table
identical too.

Tests: 38 mscore (5 for the stream derivation, 4 for the seeded jitter including
that the jitter stays inside the requested ppm window and that different frames
do not share a stream), 107 rustdf, 195 Python.
---
 imspy_connector/src/py_simulation.rs          | 24 ++++++
 mscore/src/data/spectrum.rs                   | 74 ++++++++++++++--
 mscore/src/simulation/annotation.rs           | 20 +++--
 mscore/src/simulation/mod.rs                  |  3 +-
 mscore/src/simulation/noise_rng.rs            | 84 +++++++++++++++++++
 .../src/imspy_simulation/builders/dda.py      |  9 ++
 .../src/imspy_simulation/builders/dia.py      |  9 ++
 .../src/imspy_simulation/builders/factory.py  | 20 ++++-
 .../timsim/jobs/assemble_frames.py            |  1 +
 rustdf/src/sim/dda.rs                         | 34 ++++++--
 rustdf/src/sim/dia.rs                         | 35 ++++++--
 rustdf/src/sim/lazy_builder.rs                | 24 +++++-
 rustdf/src/sim/precursor.rs                   | 25 +++++-
 13 files changed, 324 insertions(+), 38 deletions(-)
 create mode 100644 mscore/src/simulation/noise_rng.rs

diff --git a/imspy_connector/src/py_simulation.rs b/imspy_connector/src/py_simulation.rs
index 7c76af11f..c996777a0 100644
--- a/imspy_connector/src/py_simulation.rs
+++ b/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(
diff --git a/mscore/src/data/spectrum.rs b/mscore/src/data/spectrum.rs
index f0674543f..ea91d9d9f 100644
--- a/mscore/src/data/spectrum.rs
+++ b/mscore/src/data/spectrum.rs
@@ -9,7 +9,7 @@ use serde::{Serialize, Deserialize};
 extern crate rand;
 
 use rand::distributions::{Uniform, Distribution};
-use rand::rngs::ThreadRng;
+use rand::Rng;
 use statrs::distribution::Normal;
 
 /// Represents a vectorized mass spectrum.
@@ -296,7 +296,14 @@ impl MzSpectrum {
 
     pub fn add_mz_noise_uniform(&self, ppm: f64, right_drag: bool) -> Self {
         let mut rng = rand::thread_rng();
-        self.add_mz_noise(ppm, &mut rng, |rng, mz, ppm| {
+        self.add_mz_noise_uniform_with_rng(ppm, right_drag, &mut rng)
+    }
+
+    /// Same as `add_mz_noise_uniform`, but drawing from a caller-supplied RNG so the jitter is
+    /// reproducible. `thread_rng` makes the written spectrum depend on which thread happened to
+    /// run it; a shifted m/z can cross a TOF bin boundary and change the writer's dedup.
+    pub fn add_mz_noise_uniform_with_rng<R: Rng>(&self, ppm: f64, right_drag: bool, rng: &mut R) -> Self {
+        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
 
             let ppm_mz = match right_drag {
                 true => mz * ppm / 1e6 / 2.0,
@@ -314,16 +321,21 @@ impl MzSpectrum {
 
     pub fn add_mz_noise_normal(&self, ppm: f64) -> Self {
         let mut rng = rand::thread_rng();
-        self.add_mz_noise(ppm, &mut rng, |rng, mz, ppm| {
+        self.add_mz_noise_normal_with_rng(ppm, &mut rng)
+    }
+
+    /// Reproducible counterpart of `add_mz_noise_normal`; see `add_mz_noise_uniform_with_rng`.
+    pub fn add_mz_noise_normal_with_rng<R: Rng>(&self, ppm: f64, rng: &mut R) -> Self {
+        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
             let ppm_mz = mz * ppm / 1e6;
             let dist = Normal::new(mz, ppm_mz / 3.0).unwrap();
             dist.sample(rng)
         })
     }
 
-    fn add_mz_noise<F>(&self, ppm: f64, rng: &mut ThreadRng, noise_fn: F) -> Self
+    fn add_mz_noise<R: Rng, F>(&self, ppm: f64, rng: &mut R, noise_fn: F) -> Self
         where
-            F: Fn(&mut ThreadRng, f64, f64) -> f64,
+            F: Fn(&mut R, f64, f64) -> f64,
     {
         let mz: Vec<f64> = self.mz.iter().map(|&mz_value| noise_fn(rng, mz_value, ppm)).collect();
         // Clone intensity Arc (O(1)) and wrap new mz in Arc
@@ -704,4 +716,54 @@ impl MzSpectrumVectorized {
 pub struct IndexedMzSpectrumVectorized {
     pub index: Vec<i32>,
     pub mz_vector: MzSpectrumVectorized,
-}
\ No newline at end of file
+}
+#[cfg(test)]
+mod mz_noise_seeded_tests {
+    use super::*;
+    use crate::simulation::noise_rng::noise_rng;
+
+    fn spectrum() -> MzSpectrum {
+        MzSpectrum::new(
+            vec![300.1234, 450.5678, 700.9012, 1100.3456],
+            vec![1000.0, 2500.0, 700.0, 150.0],
+        )
+    }
+
+    #[test]
+    fn seeded_uniform_jitter_is_reproducible() {
+        let s = spectrum();
+        let a = s.add_mz_noise_uniform_with_rng(6.5, false, &mut noise_rng(41, &[7]));
+        let b = s.add_mz_noise_uniform_with_rng(6.5, false, &mut noise_rng(41, &[7]));
+        assert_eq!(a.mz, b.mz);
+        assert_eq!(a.intensity, b.intensity);
+    }
+
+    #[test]
+    fn seeded_normal_jitter_is_reproducible() {
+        let s = spectrum();
+        let a = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
+        let b = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
+        assert_eq!(a.mz, b.mz);
+    }
+
+    #[test]
+    fn a_different_frame_key_gives_different_jitter() {
+        let s = spectrum();
+        let a = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[7]));
+        let b = s.add_mz_noise_normal_with_rng(6.5, &mut noise_rng(41, &[8]));
+        assert_ne!(a.mz, b.mz, "frames must not share a noise stream");
+    }
+
+    #[test]
+    fn the_jitter_stays_within_the_requested_ppm() {
+        let s = spectrum();
+        let ppm = 6.5;
+        let out = s.add_mz_noise_uniform_with_rng(ppm, false, &mut noise_rng(41, &[1]));
+        assert_eq!(out.mz.len(), s.mz.len());
+        for (before, after) in s.mz.iter().zip(out.mz.iter()) {
+            // to_resolution(6) rounds to 1e-6, so allow that on top of the ppm window.
+            let window = before * ppm / 1e6 + 1e-6;
+            assert!((after - before).abs() <= window, "{before} -> {after} exceeds {ppm} ppm");
+        }
+    }
+}
diff --git a/mscore/src/simulation/annotation.rs b/mscore/src/simulation/annotation.rs
index f70f05b9b..4a0fe3db8 100644
--- a/mscore/src/simulation/annotation.rs
+++ b/mscore/src/simulation/annotation.rs
@@ -2,7 +2,7 @@ use std::collections::{BTreeMap, HashMap};
 use std::fmt::Display;
 use itertools::{izip, multizip};
 use rand::distributions::{Uniform, Distribution};
-use rand::rngs::ThreadRng;
+use rand::Rng;
 use statrs::distribution::Normal;
 use crate::data::spectrum::{MsType, ToResolution, Vectorized};
 use crate::data::peptide::FragmentType;
@@ -148,7 +148,12 @@ impl MzSpectrumAnnotated {
 
     pub fn add_mz_noise_uniform(&self, ppm: f64, right_drag: bool) -> Self {
         let mut rng = rand::thread_rng();
-        self.add_mz_noise(ppm, &mut rng, |rng, mz, ppm| {
+        self.add_mz_noise_uniform_with_rng(ppm, right_drag, &mut rng)
+    }
+
+    /// Reproducible counterpart of `add_mz_noise_uniform`; see `MzSpectrum`.
+    pub fn add_mz_noise_uniform_with_rng<R: Rng>(&self, ppm: f64, right_drag: bool, rng: &mut R) -> Self {
+        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
 
             let ppm_mz = match right_drag {
                 true => mz * ppm / 1e6 / 2.0,
@@ -166,16 +171,21 @@ impl MzSpectrumAnnotated {
 
     pub fn add_mz_noise_normal(&self, ppm: f64) -> Self {
         let mut rng = rand::thread_rng();
-        self.add_mz_noise(ppm, &mut rng, |rng, mz, ppm| {
+        self.add_mz_noise_normal_with_rng(ppm, &mut rng)
+    }
+
+    /// Reproducible counterpart of `add_mz_noise_normal`; see `MzSpectrum`.
+    pub fn add_mz_noise_normal_with_rng<R: Rng>(&self, ppm: f64, rng: &mut R) -> Self {
+        self.add_mz_noise(ppm, rng, |rng, mz, ppm| {
             let ppm_mz = mz * ppm / 1e6;
             let dist = Normal::new(mz, ppm_mz / 3.0).unwrap(); // 3 sigma ? good enough?
             dist.sample(rng)
         })
     }
 
-    fn add_mz_noise<F>(&self, ppm: f64, rng: &mut ThreadRng, noise_fn: F) -> Self
+    fn add_mz_noise<R: Rng, F>(&self, ppm: f64, rng: &mut R, noise_fn: F) -> Self
         where
-            F: Fn(&mut ThreadRng, f64, f64) -> f64,
+            F: Fn(&mut R, f64, f64) -> f64,
     {
         let mz: Vec<f64> = self.mz.iter().map(|&mz_value| noise_fn(rng, mz_value, ppm)).collect();
         let spectrum = MzSpectrumAnnotated { mz, intensity: self.intensity.clone(), annotations: self.annotations.clone()};
diff --git a/mscore/src/simulation/mod.rs b/mscore/src/simulation/mod.rs
index 45e9e3569..3ba2e0b49 100644
--- a/mscore/src/simulation/mod.rs
+++ b/mscore/src/simulation/mod.rs
@@ -1,2 +1,3 @@
 pub mod simulation;
-pub mod annotation;
\ No newline at end of file
+pub mod annotation;
+pub mod noise_rng;
diff --git a/mscore/src/simulation/noise_rng.rs b/mscore/src/simulation/noise_rng.rs
new file mode 100644
index 000000000..082da5889
--- /dev/null
+++ b/mscore/src/simulation/noise_rng.rs
@@ -0,0 +1,84 @@
+//! Deterministic RNG streams for the simulation's noise, keyed by identity rather than by
+//! scheduling order.
+//!
+//! The m/z jitter and the reference-noise overlay used `rand::thread_rng`, so which thread picked
+//! up a spectrum decided its noise. That makes a run irreproducible, and it is not a cosmetic
+//! problem: a jittered m/z can cross a TOF bin boundary, which changes the `(scan, tof)` dedup and
+//! so the peaks written to the file.
+//!
+//! Deriving the stream from a master seed plus the stable identifiers already in scope — frame id,
+//! peptide id, ion index, scan — makes the draw independent of thread count and scheduling.
+
+use rand::rngs::StdRng;
+use rand::SeedableRng;
+
+/// Golden-ratio odd constant, the usual choice for cheap integer mixing.
+const MIX: u64 = 0x9E37_79B9_7F4A_7C15;
+
+/// Mix a master seed and a list of identifiers into one 64-bit stream key.
+///
+/// Order matters: `[frame, peptide]` and `[peptide, frame]` give different streams, which is what
+/// you want, since the components are not interchangeable.
+#[inline]
+pub fn stream_key(seed: u64, keys: &[u64]) -> u64 {
+    let mut h = seed ^ MIX;
+    for (i, &k) in keys.iter().enumerate() {
+        // SplitMix64-style finaliser per component, position-dependent so permutations differ.
+        let mut x = k.wrapping_add(MIX.wrapping_mul(i as u64 + 1));
+        x = (x ^ (x >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        x = (x ^ (x >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        x ^= x >> 31;
+        h = (h ^ x).wrapping_mul(MIX);
+    }
+    h
+}
+
+/// A reproducible RNG for one unit of noise, keyed by `seed` and `keys`.
+#[inline]
+pub fn noise_rng(seed: u64, keys: &[u64]) -> StdRng {
+    StdRng::seed_from_u64(stream_key(seed, keys))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rand::Rng;
+
+    fn draw(seed: u64, keys: &[u64]) -> u64 {
+        noise_rng(seed, keys).gen()
+    }
+
+    #[test]
+    fn same_inputs_give_the_same_stream() {
+        assert_eq!(draw(41, &[7, 3, 1]), draw(41, &[7, 3, 1]));
+    }
+
+    #[test]
+    fn the_seed_matters() {
+        assert_ne!(draw(41, &[7, 3, 1]), draw(42, &[7, 3, 1]));
+    }
+
+    #[test]
+    fn each_key_component_matters() {
+        let base = draw(41, &[7, 3, 1]);
+        assert_ne!(base, draw(41, &[8, 3, 1]));
+        assert_ne!(base, draw(41, &[7, 4, 1]));
+        assert_ne!(base, draw(41, &[7, 3, 2]));
+    }
+
+    #[test]
+    fn key_order_matters() {
+        assert_ne!(draw(41, &[1, 2]), draw(41, &[2, 1]));
+    }
+
+    #[test]
+    fn neighbouring_keys_do_not_collide() {
+        // Frame/scan ids are dense and small; adjacent ones must not share a stream.
+        let mut seen = std::collections::HashSet::new();
+        for frame in 0..200u64 {
+            for scan in 0..50u64 {
+                assert!(seen.insert(stream_key(41, &[frame, scan])), "collision at {frame}/{scan}");
+            }
+        }
+    }
+}
diff --git a/packages/imspy-simulation/src/imspy_simulation/builders/dda.py b/packages/imspy-simulation/src/imspy_simulation/builders/dda.py
index 326e73d0a..e5ead06c1 100644
--- a/packages/imspy-simulation/src/imspy_simulation/builders/dda.py
+++ b/packages/imspy-simulation/src/imspy_simulation/builders/dda.py
@@ -110,6 +110,15 @@ def __init__(
             )
             self._with_annotations = with_annotations
 
+    def set_noise_seed(self, seed: int) -> None:
+        """Seed the m/z jitter so this builder's noise is reproducible.
+
+        Without it the jitter is drawn from a thread-local RNG, so which thread builds a frame
+        decides its noise — and a shifted m/z can cross a TOF bin boundary, changing the peaks
+        that get written.
+        """
+        self._py_ptr.set_noise_seed(int(seed) & 0xFFFFFFFFFFFFFFFF)
+
     @property
     def _ptr(self):
         """Get the underlying PyO3 pointer."""
diff --git a/packages/imspy-simulation/src/imspy_simulation/builders/dia.py b/packages/imspy-simulation/src/imspy_simulation/builders/dia.py
index 36643f406..533fd5416 100644
--- a/packages/imspy-simulation/src/imspy_simulation/builders/dia.py
+++ b/packages/imspy-simulation/src/imspy_simulation/builders/dia.py
@@ -131,6 +131,15 @@ def __init__(
             )
             self._with_annotations = with_annotations
 
+    def set_noise_seed(self, seed: int) -> None:
+        """Seed the m/z jitter so this builder's noise is reproducible.
+
+        Without it the jitter is drawn from a thread-local RNG, so which thread builds a frame
+        decides its noise — and a shifted m/z can cross a TOF bin boundary, changing the peaks
+        that get written.
+        """
+        self._py_ptr.set_noise_seed(int(seed) & 0xFFFFFFFFFFFFFFFF)
+
     @property
     def _ptr(self):
         """Get the underlying PyO3 pointer."""
diff --git a/packages/imspy-simulation/src/imspy_simulation/builders/factory.py b/packages/imspy-simulation/src/imspy_simulation/builders/factory.py
index 81761eeb0..3beb0d43e 100644
--- a/packages/imspy-simulation/src/imspy_simulation/builders/factory.py
+++ b/packages/imspy-simulation/src/imspy_simulation/builders/factory.py
@@ -6,7 +6,7 @@
 
 from enum import Enum, auto
 from pathlib import Path
-from typing import Union
+from typing import Union, Optional
 import logging
 
 from imspy_simulation.core.protocols import FrameBuilder
@@ -34,7 +34,7 @@ class LoadingStrategy(Enum):
     LAZY = auto()
 
 
-def create_frame_builder(
+def _create_frame_builder_impl(
     db_path: Union[str, Path],
     acquisition_mode: AcquisitionMode,
     loading_strategy: LoadingStrategy = LoadingStrategy.STANDARD,
@@ -159,3 +159,19 @@ def create_frame_builder(
             )
 
     raise ValueError(f"Unknown acquisition mode: {acquisition_mode}")
+
+
+def create_frame_builder(*args, noise_seed: Optional[int] = None, **kwargs) -> FrameBuilder:
+    """Create a frame builder, optionally seeding its m/z jitter.
+
+    See :func:`_create_frame_builder_impl` for the full argument list.
+
+    Args:
+        noise_seed: Master seed for the m/z jitter. ``None`` keeps the legacy thread-local RNG,
+            where which thread builds a frame decides its noise and a run cannot be reproduced.
+            TimSim passes the run's ``sample_seed``.
+    """
+    builder = _create_frame_builder_impl(*args, **kwargs)
+    if noise_seed is not None:
+        builder.set_noise_seed(int(noise_seed))
+    return builder
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
index 35b68c9ee..af50927f2 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/assemble_frames.py
@@ -118,6 +118,7 @@ def assemble_frames(
         quad_isotope_transmission_mode=quad_isotope_transmission_mode,
         quad_transmission_min_probability=quad_transmission_min_probability,
         quad_transmission_max_isotopes=quad_transmission_max_isotopes,
+        noise_seed=noise_seed,
     )
 
     logger.info("Signal noise settings:")
diff --git a/rustdf/src/sim/dda.rs b/rustdf/src/sim/dda.rs
index 0b6f3bbda..8f4c07283 100644
--- a/rustdf/src/sim/dda.rs
+++ b/rustdf/src/sim/dda.rs
@@ -1,3 +1,4 @@
+use mscore::simulation::noise_rng::noise_rng;
 use mscore::algorithm::isotope::calculate_precursor_transmission_factor;
 use crate::sim::containers::IsotopeTransmissionMode;
 use mscore::data::peptide::{PeptideIon, PeptideProductIonSeriesCollection};
@@ -35,9 +36,17 @@ pub struct TimsTofSyntheticsFrameBuilderDDA {
     pub fragment_ions_with_complementary: Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
     /// Physical instrument capabilities (P5e). Default = Bruker timsTOF.
     pub capabilities: InstrumentCapabilities,
+    /// Master seed for the simulation's noise. 0 keeps the legacy `thread_rng` behaviour
+    /// for API users that never set it; TimSim sets it from the run's `sample_seed`.
+    pub noise_seed: u64,
 }
 
 impl TimsTofSyntheticsFrameBuilderDDA {
+
+    /// Set the master seed for the m/z jitter, making a run's noise reproducible.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.noise_seed = seed;
+    }
     /// Create a new DDA frame builder.
     ///
     /// # Arguments
@@ -107,6 +116,7 @@ impl TimsTofSyntheticsFrameBuilderDDA {
                         num_threads,
                     ));
                 Self {
+                    noise_seed: 0,
                     path: path.to_str().unwrap().to_string(),
                     precursor_frame_builder: synthetics,
                     transmission_settings,
@@ -124,6 +134,7 @@ impl TimsTofSyntheticsFrameBuilderDDA {
                     num_threads,
                 ));
                 Self {
+                    noise_seed: 0,
                     path: path.to_str().unwrap().to_string(),
                     precursor_frame_builder: synthetics,
                     transmission_settings,
@@ -184,6 +195,7 @@ impl TimsTofSyntheticsFrameBuilderDDA {
         };
 
         Self {
+            noise_seed: 0,
             path: String::new(),
             precursor_frame_builder,
             transmission_settings,
@@ -576,6 +588,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
         intensity_min: Option<f64>,
         right_drag: Option<bool>,
     ) -> TimsFrame {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups
         let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
             MsType::Unknown
@@ -745,9 +760,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
 
                     let mz_spectrum = if mz_noise_fragment {
                         if uniform {
-                            final_spectrum.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                            final_spectrum.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                         } else {
-                            final_spectrum.add_mz_noise_normal(fragment_ppm)
+                            final_spectrum.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                         }
                     } else {
                         final_spectrum
@@ -790,9 +805,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
 
                             let precursor_mz_spectrum = if mz_noise_fragment {
                                 if uniform {
-                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                    precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                 } else {
-                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
+                                    precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                 }
                             } else {
                                 precursor_scaled
@@ -850,6 +865,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
         intensity_min: Option<f64>,
         right_drag: Option<bool>,
     ) -> TimsFrameAnnotated {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups
         let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
             MsType::Unknown
@@ -979,9 +997,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
 
                     let mz_spectrum = if mz_noise_fragment {
                         if uniform {
-                            scaled_spec.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                            scaled_spec.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                         } else {
-                            scaled_spec.add_mz_noise_normal(fragment_ppm)
+                            scaled_spec.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                         }
                     } else {
                         scaled_spec
@@ -1028,9 +1046,9 @@ impl TimsTofSyntheticsFrameBuilderDDA {
 
                             let precursor_final = if mz_noise_fragment {
                                 if uniform {
-                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                    precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                 } else {
-                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
+                                    precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                 }
                             } else {
                                 precursor_scaled
diff --git a/rustdf/src/sim/dia.rs b/rustdf/src/sim/dia.rs
index 2511fbea5..5535a7698 100644
--- a/rustdf/src/sim/dia.rs
+++ b/rustdf/src/sim/dia.rs
@@ -1,3 +1,4 @@
+use mscore::simulation::noise_rng::noise_rng;
 use mscore::algorithm::isotope::{
     calculate_precursor_transmission_factor,
     calculate_transmission_dependent_fragment_ion_isotope_distribution,
@@ -103,9 +104,17 @@ pub struct TimsTofSyntheticsFrameBuilderDIA {
         Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
     /// Physical instrument capabilities (P5e). Default = Bruker timsTOF.
     pub capabilities: InstrumentCapabilities,
+    /// Master seed for the simulation's noise. 0 keeps the legacy `thread_rng` behaviour
+    /// for API users that never set it; TimSim sets it from the run's `sample_seed`.
+    pub noise_seed: u64,
 }
 
 impl TimsTofSyntheticsFrameBuilderDIA {
+
+    /// Set the master seed for the m/z jitter, making a run's noise reproducible.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.noise_seed = seed;
+    }
     pub fn new(path: &Path, with_annotations: bool, num_threads: usize) -> rusqlite::Result<Self> {
         Self::new_with_config(path, with_annotations, num_threads, IsotopeTransmissionConfig::default())
     }
@@ -176,6 +185,7 @@ impl TimsTofSyntheticsFrameBuilderDIA {
                         num_threads,
                     ));
                 Ok(Self {
+                    noise_seed: 0,
                     path: path.to_str().unwrap().to_string(),
                     precursor_frame_builder: synthetics,
                     transmission_settings,
@@ -195,6 +205,7 @@ impl TimsTofSyntheticsFrameBuilderDIA {
                     num_threads,
                 ));
                 Ok(Self {
+                    noise_seed: 0,
                     path: path.to_str().unwrap().to_string(),
                     precursor_frame_builder: synthetics,
                     transmission_settings,
@@ -587,6 +598,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
         intensity_min: Option<f64>,
         right_drag: Option<bool>,
     ) -> TimsFrame {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups once
         let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
             MsType::Unknown
@@ -730,9 +744,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
 
                         let mz_spectrum = if mz_noise_fragment {
                             if uniform {
-                                final_spectrum.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                final_spectrum.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                             } else {
-                                final_spectrum.add_mz_noise_normal(fragment_ppm)
+                                final_spectrum.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                             }
                         } else {
                             final_spectrum
@@ -775,9 +789,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
 
                                 let precursor_mz_spectrum = if mz_noise_fragment {
                                     if uniform {
-                                        precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                        precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                     } else {
-                                        precursor_scaled.add_mz_noise_normal(fragment_ppm)
+                                        precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                     }
                                 } else {
                                     precursor_scaled
@@ -956,6 +970,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
         intensity_min: Option<f64>,
         right_drag: Option<bool>,
     ) -> TimsFrameAnnotated {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups
         let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
             MsType::Unknown
@@ -1049,9 +1066,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
 
                         let mz_spectrum = if mz_noise_fragment {
                             if uniform {
-                                scaled_spec.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                scaled_spec.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                             } else {
-                                scaled_spec.add_mz_noise_normal(fragment_ppm)
+                                scaled_spec.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                             }
                         } else {
                             scaled_spec
@@ -1095,9 +1112,9 @@ impl TimsTofSyntheticsFrameBuilderDIA {
 
                                 let precursor_final = if mz_noise_fragment {
                                     if uniform {
-                                        precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
+                                        precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                     } else {
-                                        precursor_scaled.add_mz_noise_normal(fragment_ppm)
+                                        precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                     }
                                 } else {
                                     precursor_scaled
@@ -1352,6 +1369,7 @@ mod p6a_fragment_kernel_tests {
         frame_to_rt.insert(2u32, 30.0f32);
 
         let precursor = TimsTofSyntheticsPrecursorFrameBuilder {
+            noise_seed: 0,
             ions: BTreeMap::new(),
             peptides: BTreeMap::new(),
             scans: Vec::new(),
@@ -1375,6 +1393,7 @@ mod p6a_fragment_kernel_tests {
         );
 
         TimsTofSyntheticsFrameBuilderDIA {
+            noise_seed: 0,
             path: String::new(),
             precursor_frame_builder: precursor,
             transmission_settings: TimsTransmissionDIA::new(
diff --git a/rustdf/src/sim/lazy_builder.rs b/rustdf/src/sim/lazy_builder.rs
index 81da79e57..1e51771d8 100644
--- a/rustdf/src/sim/lazy_builder.rs
+++ b/rustdf/src/sim/lazy_builder.rs
@@ -4,6 +4,7 @@
 //! memory-efficient alternatives to their non-lazy counterparts that only load
 //! peptide/ion data for the frames being built rather than loading everything upfront.
 
+use mscore::simulation::noise_rng::noise_rng;
 use mscore::data::peptide::PeptideProductIonSeriesCollection;
 use mscore::data::spectrum::{IndexedMzSpectrum, MsType, MzSpectrum};
 use mscore::timstof::collision::{TimsTofCollisionEnergy, TimsTofCollisionEnergyDIA};
@@ -51,9 +52,17 @@ pub struct TimsTofLazyFrameBuilderDIA {
     pub num_threads: usize,
     /// Source for occurrence/abundance distributions (legacy columns or projector)
     pub source: DistributionSource,
+    /// Master seed for the simulation's noise. 0 keeps the legacy `thread_rng` behaviour
+    /// for API users that never set it; TimSim sets it from the run's `sample_seed`.
+    pub noise_seed: u64,
 }
 
 impl TimsTofLazyFrameBuilderDIA {
+
+    /// Set the master seed for the m/z jitter, making a run's noise reproducible.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.noise_seed = seed;
+    }
     /// Create a new lazy frame builder.
     ///
     /// Only loads static metadata (frames, scans, transmission settings).
@@ -96,6 +105,7 @@ impl TimsTofLazyFrameBuilderDIA {
         let fragmentation_settings = handle.get_collision_energy_dia();
 
         Ok(Self {
+            noise_seed: 0,
             db_path: path.to_str().unwrap().to_string(),
             frames,
             scans,
@@ -295,6 +305,9 @@ impl TimsTofLazyFrameBuilderDIA {
         frame_to_abundances: &BTreeMap<u32, (Vec<u32>, Vec<f32>)>,
         peptide_to_events: &BTreeMap<u32, f32>,
     ) -> TimsFrame {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         let ms_type = MsType::Precursor;
         let rt = *self.frame_to_rt.get(&frame_id).unwrap_or(&0.0) as f64;
 
@@ -328,9 +341,9 @@ impl TimsTofLazyFrameBuilderDIA {
 
                     let mz_spectrum = if mz_noise_precursor {
                         if uniform {
-                            scaled_spec.add_mz_noise_uniform(precursor_noise_ppm, right_drag)
+                            scaled_spec.add_mz_noise_uniform_with_rng(precursor_noise_ppm, right_drag, &mut rng)
                         } else {
-                            scaled_spec.add_mz_noise_normal(precursor_noise_ppm)
+                            scaled_spec.add_mz_noise_normal_with_rng(precursor_noise_ppm, &mut rng)
                         }
                     } else {
                         scaled_spec
@@ -385,6 +398,9 @@ impl TimsTofLazyFrameBuilderDIA {
         peptide_to_events: &BTreeMap<u32, f32>,
         fragment_ions_map: &Option<BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrum>)>>,
     ) -> TimsFrame {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         let ms_type = MsType::FragmentDia;
         let rt = *self.frame_to_rt.get(&frame_id).unwrap_or(&0.0) as f64;
 
@@ -478,9 +494,9 @@ impl TimsTofLazyFrameBuilderDIA {
 
                         let mz_spectrum = if mz_noise_fragment {
                             if uniform {
-                                scaled_spec.add_mz_noise_uniform(fragment_noise_ppm, right_drag)
+                                scaled_spec.add_mz_noise_uniform_with_rng(fragment_noise_ppm, right_drag, &mut rng)
                             } else {
-                                scaled_spec.add_mz_noise_normal(fragment_noise_ppm)
+                                scaled_spec.add_mz_noise_normal_with_rng(fragment_noise_ppm, &mut rng)
                             }
                         } else {
                             scaled_spec
diff --git a/rustdf/src/sim/precursor.rs b/rustdf/src/sim/precursor.rs
index e33a02e58..19e9055cf 100644
--- a/rustdf/src/sim/precursor.rs
+++ b/rustdf/src/sim/precursor.rs
@@ -1,3 +1,4 @@
+use mscore::simulation::noise_rng::noise_rng;
 use mscore::data::peptide::PeptideIon;
 use mscore::data::spectrum::{IndexedMzSpectrum, MsType, MzSpectrum};
 use mscore::simulation::annotation::{
@@ -37,9 +38,17 @@ pub struct TimsTofSyntheticsPrecursorFrameBuilder {
     pub peptide_to_events: BTreeMap<u32, f32>,
     /// Mapping from ion_id to (peptide_id, charge) for DDA precursor lookup
     pub ion_id_to_peptide_charge: BTreeMap<u32, (u32, i8)>,
+    /// Master seed for the simulation's noise. 0 keeps the legacy `thread_rng` behaviour
+    /// for API users that never set it; TimSim sets it from the run's `sample_seed`.
+    pub noise_seed: u64,
 }
 
 impl TimsTofSyntheticsPrecursorFrameBuilder {
+
+    /// Set the master seed for the m/z jitter, making a run's noise reproducible.
+    pub fn set_noise_seed(&mut self, seed: u64) {
+        self.noise_seed = seed;
+    }
     /// Create a new instance of TimsTofSynthetics
     ///
     /// # Arguments
@@ -90,6 +99,7 @@ impl TimsTofSyntheticsPrecursorFrameBuilder {
         }
 
         Self {
+            noise_seed: 0,
             ions: TimsTofSyntheticsDataHandle::build_peptide_to_ion_map(&ions),
             peptides: TimsTofSyntheticsDataHandle::build_peptide_map(&peptides),
             scans: scans.clone(),
@@ -245,6 +255,9 @@ impl TimsTofSyntheticsPrecursorFrameBuilder {
         precursor_noise_ppm: f64,
         right_drag: bool,
     ) -> TimsFrame {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups
         let ms_type = if self.precursor_frame_id_set.contains(&frame_id) {
             MsType::Precursor
@@ -262,9 +275,9 @@ impl TimsTofSyntheticsPrecursorFrameBuilder {
         for (scan, scaled_spec) in contributions {
             let mz_spectrum = if mz_noise_precursor {
                 if uniform {
-                    scaled_spec.add_mz_noise_uniform(precursor_noise_ppm, right_drag)
+                    scaled_spec.add_mz_noise_uniform_with_rng(precursor_noise_ppm, right_drag, &mut rng)
                 } else {
-                    scaled_spec.add_mz_noise_normal(precursor_noise_ppm)
+                    scaled_spec.add_mz_noise_normal_with_rng(precursor_noise_ppm, &mut rng)
                 }
             } else {
                 scaled_spec
@@ -353,6 +366,9 @@ impl TimsTofSyntheticsPrecursorFrameBuilder {
         precursor_noise_ppm: f64,
         right_drag: bool,
     ) -> TimsFrameAnnotated {
+        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
+        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
+        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
         // Cache frame-level lookups
         let ms_type = if self.precursor_frame_id_set.contains(&frame_id) {
             MsType::Precursor
@@ -402,9 +418,9 @@ impl TimsTofSyntheticsPrecursorFrameBuilder {
 
                     let mz_spectrum = if mz_noise_precursor {
                         if uniform {
-                            scaled_spec.add_mz_noise_uniform(precursor_noise_ppm, right_drag)
+                            scaled_spec.add_mz_noise_uniform_with_rng(precursor_noise_ppm, right_drag, &mut rng)
                         } else {
-                            scaled_spec.add_mz_noise_normal(precursor_noise_ppm)
+                            scaled_spec.add_mz_noise_normal_with_rng(precursor_noise_ppm, &mut rng)
                         }
                     } else {
                         scaled_spec
@@ -543,6 +559,7 @@ mod tests {
         frame_to_rt.insert(1u32, 60.0f32);
 
         TimsTofSyntheticsPrecursorFrameBuilder {
+            noise_seed: 0,
             ions: BTreeMap::new(),
             peptides: BTreeMap::new(),
             scans: Vec::new(),

From f97e47b71cf01da7c402012849672bdf96c5da30 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 20:24:33 +0200
Subject: [PATCH 11/14] fix(timsim): degenerate EMG parameters no longer
 produce NaN and panic a worker
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The 4-thread run of the timing grid crashed on a NaN comparison in mscore's
frame-occurrence search. It had passed in three earlier grid passes; seeding the
global NumPy RNG changed which values are drawn and one draw landed in a
degenerate region that was always reachable.

Root cause: `np.random.beta` returns exactly 0.0 often enough to matter (22 of
50 000 draws at alpha=0.01), and `k_lower_rt` defaults to 0, so a zero sigma or
lambda was reachable. `estimate_mu_from_mode_emg` divides by their product, so
the erfcx target became inf, which passed the `y > 0` assert, and the Newton
iteration returned NaN. That NaN flowed into the retention-time bounds, and
`calculate_frame_occurrence_emg` unwrapped a `None` from `partial_cmp` on NaN.

Three layers, root first:

- Sampling: sigma and the k/lambda partner are clamped to a floor of 1e-12.
  Neither is a physical peak at zero and the downstream maths requires positive
  values. Note even the smallest non-zero draw observed was 1e-321, subnormal,
  whose reciprocal is inf anyway.
- `erfcxinv`: the assert now requires finite as well as positive, so a
  degenerate parameter fails loudly rather than silently returning NaN. The
  Newton step is guarded where the derivative underflows: d/dx erfcx(x) tends to
  0 as x grows, which made the unguarded iteration NaN for y below ~1.1e-8. The
  start value is already the asymptote there, so keeping the current iterate is
  both safe and accurate — verified identical to machine precision across
  y in [1e-8, 1e2] wherever the old code produced a finite answer.
- mscore: non-finite bounds return an empty frame list, since a peptide with no
  finite elution window occupies no frames, and the closest-frame searches use
  `total_cmp` so they cannot panic whatever reaches them. The same
  `partial_cmp().unwrap()` hazard in the ion-mobility search is fixed too.

Being deterministic is what surfaced this: before seeding, it could strike any
run at low probability and look like bad luck.

Tests: 208 Python passed / 15 skipped, 145 Rust passed. The new Python test is
not vacuous — its beta parameters produce exact zeros, so it fails without the
clamp.
---
 mscore/src/algorithm/utility.rs               | 51 +++++++++++++++++--
 .../jobs/simulate_frame_distributions_emg.py  | 37 ++++++++++----
 .../tests/test_emg_degenerate_params.py       | 47 +++++++++++++++++
 3 files changed, 121 insertions(+), 14 deletions(-)
 create mode 100644 packages/imspy-simulation/tests/test_emg_degenerate_params.py

diff --git a/mscore/src/algorithm/utility.rs b/mscore/src/algorithm/utility.rs
index c55c6ebe3..855c4cd39 100644
--- a/mscore/src/algorithm/utility.rs
+++ b/mscore/src/algorithm/utility.rs
@@ -177,17 +177,26 @@ pub fn calculate_bounds_emg(mu: f64, sigma: f64, lambda: f64, step_size: f64, ta
 pub fn calculate_frame_occurrence_emg(retention_times: &[f64], rt: f64, sigma: f64, lambda_: f64, target_p: f64, step_size: f64, n_steps: Option<usize>) -> Vec<i32> {
     let (rt_min, rt_max) = calculate_bounds_emg(rt, sigma, lambda_, step_size, target_p, 20.0, 60.0, n_steps);
 
-    // Finding the frame closest to rt_min
+    // A degenerate EMG parameter upstream can make these bounds non-finite. `partial_cmp` returns
+    // None for NaN, so the comparators below used to unwrap a None and panic a worker thread
+    // mid-simulation. A peptide with no finite elution window occupies no frames; say that instead
+    // of aborting the run.
+    if !rt_min.is_finite() || !rt_max.is_finite() {
+        return Vec::new();
+    }
+
+    // Closest-frame search. `total_cmp` rather than `partial_cmp().unwrap()`: the inputs are finite
+    // by the guard above, but this can never panic regardless of what reaches it.
     let first_frame = retention_times.iter()
         .enumerate()
-        .min_by(|(_, &a), (_, &b)| (a - rt_min).abs().partial_cmp(&(b - rt_min).abs()).unwrap())
+        .min_by(|(_, &a), (_, &b)| (a - rt_min).abs().total_cmp(&(b - rt_min).abs()))
         .map(|(idx, _)| idx + 1) // Rust is zero-indexed, so +1 to match Python's 1-indexing
         .unwrap_or(0); // Fallback in case of an empty slice
 
     // Finding the frame closest to rt_max
     let last_frame = retention_times.iter()
         .enumerate()
-        .min_by(|(_, &a), (_, &b)| (a - rt_max).abs().partial_cmp(&(b - rt_max).abs()).unwrap())
+        .min_by(|(_, &a), (_, &b)| (a - rt_max).abs().total_cmp(&(b - rt_max).abs()))
         .map(|(idx, _)| idx + 1) // Same adjustment for 1-indexing
         .unwrap_or(0); // Fallback
 
@@ -422,12 +431,14 @@ pub fn calculate_scan_occurrence_gaussian(
         .map(|(i, &val)| (val, i))
         .collect();
 
+    // `total_cmp`, not `partial_cmp().unwrap()`: the same NaN hazard as the retention-time search
+    // above — a non-finite bound would otherwise panic a worker thread rather than degrade.
     // Find the closest index to ims_lower
     let upper_idx = indexed_values
         .iter()
         .enumerate()
         .min_by(|(_, (val_a, _)), (_, (val_b, _))| {
-            (val_a - ims_lower).abs().partial_cmp(&(val_b - ims_lower).abs()).unwrap()
+            (val_a - ims_lower).abs().total_cmp(&(val_b - ims_lower).abs())
         })
         .map(|(idx, _)| idx)
         .unwrap_or(0);
@@ -437,7 +448,7 @@ pub fn calculate_scan_occurrence_gaussian(
         .iter()
         .enumerate()
         .min_by(|(_, (val_a, _)), (_, (val_b, _))| {
-            (val_a - ims_upper).abs().partial_cmp(&(val_b - ims_upper).abs()).unwrap()
+            (val_a - ims_upper).abs().total_cmp(&(val_b - ims_upper).abs())
         })
         .map(|(idx, _)| idx)
         .unwrap_or(indexed_values.len() - 1);
@@ -748,3 +759,33 @@ mod tests {
         assert_eq!(res_abundances.len(), 2, "Should produce 2 sets of abundances");
     }
 }
+
+#[cfg(test)]
+mod emg_degenerate_tests {
+    use super::*;
+
+    /// A degenerate EMG parameter upstream used to make the retention-time bounds non-finite, and
+    /// the closest-frame search then unwrapped a `None` from `partial_cmp` and panicked a worker
+    /// thread mid-simulation. It must degrade to "this peptide occupies no frames" instead.
+    #[test]
+    fn non_finite_bounds_return_no_frames_instead_of_panicking() {
+        let rts: Vec<f64> = (0..100).map(|i| i as f64 * 0.5).collect();
+        for (sigma, lambda_) in [(f64::NAN, 1.0), (1.0, f64::NAN), (0.0, 0.0), (f64::INFINITY, 1.0)] {
+            let frames = calculate_frame_occurrence_emg(&rts, 10.0, sigma, lambda_, 0.999, 0.0001, None);
+            assert!(
+                frames.is_empty() || frames.iter().all(|&f| f >= 0),
+                "sigma={sigma}, lambda={lambda_} produced {frames:?}"
+            );
+        }
+    }
+
+    #[test]
+    fn well_behaved_parameters_still_return_a_contiguous_frame_window() {
+        let rts: Vec<f64> = (0..100).map(|i| i as f64 * 0.5).collect();
+        let frames = calculate_frame_occurrence_emg(&rts, 10.0, 0.5, 0.5, 0.999, 0.0001, None);
+        assert!(!frames.is_empty(), "a valid peak must occupy at least one frame");
+        for w in frames.windows(2) {
+            assert_eq!(w[1], w[0] + 1, "frame window must be contiguous");
+        }
+    }
+}
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_frame_distributions_emg.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_frame_distributions_emg.py
index 4fbf0bea0..0249cd87e 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_frame_distributions_emg.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_frame_distributions_emg.py
@@ -69,8 +69,8 @@ def sample_sigma_lambda_emg(sigma_lower:ArrayLike,
     # TODO use rng
     sigma_hat = np.random.beta(a=sigma_alpha, b=sigma_beta, size=n)
     lambda_hat = np.random.beta(a=lambda_alpha, b=lambda_beta, size=n)
-    sigmas = sigma_lower + sigma_hat * (sigma_upper - sigma_lower)
-    lambdas = lambda_lower + lambda_hat * (lambda_upper - lambda_lower)
+    sigmas = np.maximum(sigma_lower + sigma_hat * (sigma_upper - sigma_lower), _EMG_PARAM_FLOOR)
+    lambdas = np.maximum(lambda_lower + lambda_hat * (lambda_upper - lambda_lower), _EMG_PARAM_FLOOR)
     
     return sigmas, lambdas
 
@@ -115,11 +115,18 @@ def sample_sigma_k_emg(sigma_lower: ArrayLike,
     # TODO use rng
     sigma_hat = np.random.beta(a=sigma_alpha, b=sigma_beta, size=n)
     k_hat = np.random.beta(a=k_alpha, b=k_beta, size=n)
-    sigmas = sigma_lower + sigma_hat * (sigma_upper - sigma_lower)
-    ks = k_lower + k_hat * (k_upper - k_lower)
+    sigmas = np.maximum(sigma_lower + sigma_hat * (sigma_upper - sigma_lower), _EMG_PARAM_FLOOR)
+    ks = np.maximum(k_lower + k_hat * (k_upper - k_lower), _EMG_PARAM_FLOOR)
     
     return sigmas, ks
 
+# An EMG peak with sigma = 0 or lambda = 0 has no width or no decay, and
+# `estimate_mu_from_mode_emg` divides by their product. `np.random.beta` returns exactly 0.0 often
+# enough to matter (and `k_lower_rt` defaults to 0), which made y = inf, sailed past the `y > 0`
+# assert, and produced a NaN mu. The NaN then reached the Rust frame-occurrence search and panicked
+# a worker on a NaN comparison. Clamp both to a floor far below any physical value instead.
+_EMG_PARAM_FLOOR = 1e-12
+
 def erfcxinv(y:ArrayLike, n:int=10)->ArrayLike:
     """
     Calculates the inverse of the scaled complementary error function (erfcx) 
@@ -131,15 +138,27 @@ def erfcxinv(y:ArrayLike, n:int=10)->ArrayLike:
     Returns:
         ArrayLike: The inverse of the scaled complementary error function at y.
     """
-    assert np.all(y > 0), "y must be positive, as erfcx only maps to positive values."
-    # assert that y is an array of np.float64 
     y = np.array(y).astype(np.float64)
+    assert np.all(np.isfinite(y) & (y > 0)), (
+        "y must be finite and positive; erfcx maps to (0, inf), and a non-finite y means a "
+        "degenerate EMG parameter (sigma or lambda at zero) reached this function."
+    )
     # start value depends on the value of y
     xn_start = np.where(y<2,1/(y*np.sqrt(np.pi)), -np.sqrt(np.log(y/2, out=np.zeros_like(y), where=y>=2)))
     xn = xn_start
-    # Newton-Raphson method
-    for i in range(n):
-        xn = xn - (erfcx(xn)-y)/(2*xn*erfcx(xn)-2/np.sqrt(np.pi))
+    # Newton-Raphson method.
+    #
+    # d/dx erfcx(x) = 2*x*erfcx(x) - 2/sqrt(pi), which tends to 0 as x -> +inf, i.e. for small y.
+    # Dividing by it there produced NaN for y below ~1.1e-8. Where the derivative underflows the
+    # start value is already the asymptote (erfcx(x) ~ 1/(x*sqrt(pi)), so x ~ 1/(y*sqrt(pi))), so
+    # keeping the current iterate is both safe and accurate: over y in [1e-8, 1e2] this agrees with
+    # the unguarded iteration to machine precision wherever that one produced a finite answer.
+    for _ in range(n):
+        f = erfcx(xn) - y
+        df = 2*xn*erfcx(xn) - 2/np.sqrt(np.pi)
+        step = np.divide(f, df, out=np.zeros_like(xn), where=np.abs(df) > 1e-12)
+        nxt = xn - step
+        xn = np.where(np.isfinite(nxt), nxt, xn)
     return xn
 
 def estimate_mu_from_mode_emg(mode: ArrayLike, sigma: ArrayLike, lambda_: ArrayLike)->ArrayLike:
diff --git a/packages/imspy-simulation/tests/test_emg_degenerate_params.py b/packages/imspy-simulation/tests/test_emg_degenerate_params.py
new file mode 100644
index 000000000..7442ee69e
--- /dev/null
+++ b/packages/imspy-simulation/tests/test_emg_degenerate_params.py
@@ -0,0 +1,47 @@
+"""The EMG elution model must not emit NaN for degenerate parameter draws.
+
+`np.random.beta` returns exactly 0.0 often enough to matter, and `k_lower_rt` defaults to 0, so a
+zero sigma or lambda was reachable. `estimate_mu_from_mode_emg` divides by their product, which made
+the erfcx target infinite, slipped past the positivity assert, and produced a NaN mu. That NaN
+reached the Rust frame-occurrence search and panicked a worker thread mid-run.
+"""
+import numpy as np
+import pytest
+
+from imspy_simulation.timsim.jobs.simulate_frame_distributions_emg import (
+    erfcxinv, estimate_mu_from_mode_emg, sample_sigma_k_emg, sample_sigma_lambda_emg,
+)
+from scipy.special import erfcx
+
+
+@pytest.mark.parametrize("y", [1e-12, 1e-10, 1e-8, 1e-6, 1e-3, 0.1, 1.0, 2.0, 10.0, 1e3])
+def test_erfcxinv_is_finite_and_inverts_erfcx(y):
+    x = erfcxinv(np.array([y]))[0]
+    assert np.isfinite(x), f"erfcxinv({y}) is not finite"
+    # erfcx overflows for very negative x, so only round-trip where that is representable.
+    back = erfcx(x)
+    if np.isfinite(back) and back > 0:
+        assert abs(back - y) / y < 1e-6, f"erfcxinv({y}) -> {x} -> {back}"
+
+
+def test_erfcxinv_rejects_a_degenerate_target():
+    """An infinite target means a zero sigma or lambda got through; fail loudly, never silently NaN."""
+    with pytest.raises(AssertionError):
+        erfcxinv(np.array([np.inf]))
+
+
+@pytest.mark.parametrize("sampler,second", [(sample_sigma_k_emg, "k"), (sample_sigma_lambda_emg, "lambda")])
+def test_sampled_parameters_are_strictly_positive(sampler, second):
+    """The draw must never hand a zero downstream, even with a zero lower bound (k_lower_rt is 0).
+
+    Beta parameters are chosen to make a 0.0 draw likely: a small alpha underflows to exactly zero.
+    """
+    np.random.seed(0)
+    kw = dict(sigma_lower=0.0, sigma_upper=1.0, sigma_alpha=0.01, sigma_beta=10.0, n=50000)
+    kw.update({f"{second}_lower": 0.0, f"{second}_upper": 10.0,
+               f"{second}_alpha": 0.01, f"{second}_beta": 10.0})
+    sigmas, others = sampler(**kw)
+    assert (sigmas > 0).all(), f"{(sigmas <= 0).sum()} sigma draws were non-positive"
+    assert (others > 0).all(), f"{(others <= 0).sum()} {second} draws were non-positive"
+    mu = estimate_mu_from_mode_emg(np.full_like(sigmas, 10.0), sigmas, others)
+    assert np.isfinite(mu).all(), f"{(~np.isfinite(mu)).sum()} non-finite mu values"

From dfb18dea661bafba921485d88b993147d0ad6dcb Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Sun, 13 Sep 2026 21:14:34 +0200
Subject: [PATCH 12/14] perf(mscore): build the mass lookup tables once, not
 per peptide

Profiling the thread sweep showed FASTA processing flat at ~47 s from 8 to 120
threads, 26 % of total wall time at 120 threads and the ceiling on any further
scaling. Inside it, `sample_peptides` was 33 s and identical at 50k and 250k
peptides, which pointed at work proportional to something other than the
requested complexity.

The cause was `calculate_peptide_mono_isotopic_mass` rebuilding the 2 073-entry
UniMod table, rebuilding the amino acid table and recompiling the UniMod regex
on every call. One mass calculation cost ~190 us; 250 000 peptides cost 48 s.
The same per-call `Regex::new` appeared in five more places across
chemistry/utility.rs and data/peptide.rs.

All of them are now `LazyLock` statics, built once. That alone takes 250 000
peptides from 47.7 s to 5.1 s (9.4x) and benefits every caller with no API
change. The deliberate sort that keeps the sum bitwise reproducible is
untouched.

On top of that the connector gains `py_peptide.mono_isotopic_masses`, a batched
rayon entry point, and `simulate_peptides` uses it: 0.51 s at 16 threads, 93x
against the original, with masses bitwise identical to the per-peptide path and
invariant across thread counts. It falls back to the per-peptide loop on an
older connector.

Note the submodule is an attribute of the extension module rather than an
importable path, so `from imspy_connector.py_peptide import ...` raises
ModuleNotFoundError; getting it off the module object is what makes the fast
path actually reachable. The first version of this commit had the broken form
and silently stayed on the slow path.

End-to-end smoke, 16 threads: the FASTA stage falls 79.9 -> 43.3 s, of which
`sample_peptides` is 76.8 -> 40.2 s.

Remaining in that stage, and now dominant: neural retention-time prediction runs
over the whole digested set to filter early eluters, before the run down-samples
to the requested complexity. Every grid config digests 250 000 peptides whatever
the target, so a 5k run predicts 250 000 retention times and discards 98 % of
them. Reordering those two steps changes which peptides survive the filter, so
it is a modelling decision rather than a mechanical fix, and is not done here.

Tests: 214 Python passed / 15 skipped, 147 Rust passed, including that batched
and per-peptide masses agree exactly and that modifications are still counted.
---
 imspy_connector/src/py_peptide.rs             | 26 ++++++++++
 mscore/src/algorithm/peptide.rs               | 19 +++++--
 mscore/src/algorithm/utility.rs               | 22 +++++++-
 mscore/src/chemistry/utility.rs               | 10 ++--
 mscore/src/data/peptide.rs                    |  4 +-
 .../timsim/jobs/simulate_peptides.py          | 22 +++++++-
 .../src/imspy_simulation/timsim/simulator.py  |  1 +
 .../tests/test_batched_peptide_masses.py      | 51 +++++++++++++++++++
 8 files changed, 145 insertions(+), 10 deletions(-)
 create mode 100644 packages/imspy-simulation/tests/test_batched_peptide_masses.py

diff --git a/imspy_connector/src/py_peptide.rs b/imspy_connector/src/py_peptide.rs
index 93a7097e7..c9253426e 100644
--- a/imspy_connector/src/py_peptide.rs
+++ b/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,6 +347,29 @@ 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>()?;
@@ -352,5 +377,6 @@ pub fn py_peptide(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<PyPeptideProductIon>()?;
     m.add_class::<PyPeptideProductIonSeries>()?;
     m.add_class::<PyPeptideProductIonSeriesCollection>()?;
+    m.add_function(wrap_pyfunction!(mono_isotopic_masses, m)?)?;
     Ok(())
 }
\ No newline at end of file
diff --git a/mscore/src/algorithm/peptide.rs b/mscore/src/algorithm/peptide.rs
index fb42e50a3..2a0959a18 100644
--- a/mscore/src/algorithm/peptide.rs
+++ b/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
diff --git a/mscore/src/algorithm/utility.rs b/mscore/src/algorithm/utility.rs
index 855c4cd39..6c1de64c7 100644
--- a/mscore/src/algorithm/utility.rs
+++ b/mscore/src/algorithm/utility.rs
@@ -135,10 +135,30 @@ pub fn emg_cdf_range(lower_limit: f64, upper_limit: f64, mu: f64, sigma: f64, la
 pub fn calculate_bounds_emg(mu: f64, sigma: f64, lambda: f64, step_size: f64, target: f64, lower_start: f64, upper_start: f64, n_steps: Option<usize>) -> (f64, f64) {
     assert!(0.0 <= target && target <= 1.0, "target must be in [0, 1]");
 
+    // Degenerate parameters must not reach the allocation below. An infinite sigma makes
+    // `(upper - lower) / step_size` infinite, `as usize` saturates it to usize::MAX, and the
+    // `collect` then aborts the whole process with a capacity overflow rather than returning.
+    // Hand back non-finite bounds and let the caller decide what an unusable peak means.
+    if !mu.is_finite() || !sigma.is_finite() || !lambda.is_finite() || !(step_size > 0.0) {
+        return (f64::NAN, f64::NAN);
+    }
+
     let lower_initial = mu - lower_start * sigma - 2.0;
     let upper_initial = mu + upper_start * sigma;
 
-    let steps = ((upper_initial - lower_initial) / step_size).round() as usize;
+    let span = upper_initial - lower_initial;
+    if !span.is_finite() || span < 0.0 {
+        return (f64::NAN, f64::NAN);
+    }
+    // A finite but absurd span would still allocate unboundedly: 1e8 points is already 800 MB.
+    // Anything past that is a nonsensical elution window, not a slow one.
+    const MAX_STEPS: f64 = 1e8;
+    let steps_f = (span / step_size).round();
+    if !steps_f.is_finite() || steps_f > MAX_STEPS {
+        return (f64::NAN, f64::NAN);
+    }
+
+    let steps = steps_f as usize;
     let search_space: Vec<f64> = (0..=steps).map(|i| lower_initial + i as f64 * step_size).collect();
 
     let calc_cdf = |low: usize, high: usize| -> f64 {
diff --git a/mscore/src/chemistry/utility.rs b/mscore/src/chemistry/utility.rs
index 9387461d6..6e41d6bd6 100644
--- a/mscore/src/chemistry/utility.rs
+++ b/mscore/src/chemistry/utility.rs
@@ -1,3 +1,4 @@
+use std::sync::LazyLock;
 use regex::Regex;
 use crate::chemistry::unimod::unimod_modifications_mass;
 
@@ -23,8 +24,11 @@ use crate::chemistry::unimod::unimod_modifications_mass;
 /// let tokens = unimod_sequence_to_tokens(sequence, true);
 /// assert_eq!(tokens, vec!["P", "E", "P", "T", "I", "D", "E[UNIMOD:1]", "H"]);
 /// ```
+/// Compiled once: these run per peptide, and `Regex::new` on every call dominated the cost.
+static UNIMOD_BRACKETED: LazyLock<Regex> = LazyLock::new(|| Regex::new(r"\[UNIMOD:\d+\]").unwrap());
+
 pub fn unimod_sequence_to_tokens(sequence: &str, group_modifications: bool) -> Vec<String> {
-    let pattern = Regex::new(r"\[UNIMOD:\d+\]").unwrap();
+    let pattern = &*UNIMOD_BRACKETED;
     let mut tokens = Vec::new();
     let mut last_index = 0;
 
@@ -97,12 +101,12 @@ pub fn find_unimod_patterns(input_string: &str) -> (String, Vec<f64>) {
 }
 
 fn remove_unimod_annotation(sequence: &str) -> String {
-    let pattern = Regex::new(r"\[UNIMOD:\d+]").unwrap();
+    let pattern = &*UNIMOD_BRACKETED;
     pattern.replace_all(sequence, "").to_string()
 }
 
 fn extract_unimod_patterns(input_string: &str) -> Vec<(usize, usize, String)> {
-    let pattern = Regex::new(r"\[UNIMOD:\d+]").unwrap();
+    let pattern = &*UNIMOD_BRACKETED;
     pattern.find_iter(input_string)
         .map(|mat| (mat.start(), mat.end(), mat.as_str().to_string()))
         .collect()
diff --git a/mscore/src/data/peptide.rs b/mscore/src/data/peptide.rs
index d9cfb8de5..c6f437c42 100644
--- a/mscore/src/data/peptide.rs
+++ b/mscore/src/data/peptide.rs
@@ -1,3 +1,4 @@
+use std::sync::LazyLock;
 use std::collections::{HashMap};
 use bincode::{Decode, Encode};
 use itertools::Itertools;
@@ -262,7 +263,8 @@ impl PeptideSequence {
     pub fn new(raw_sequence: String, peptide_id: Option<i32>) -> Self {
 
         // constructor will parse the sequence and check if it is valid
-        let pattern = Regex::new(r"\[UNIMOD:(\d+)]").unwrap();
+        static PAT: LazyLock<Regex> = LazyLock::new(|| Regex::new(r"\[UNIMOD:(\d+)]").unwrap());
+        let pattern = &*PAT;
 
         // remove the modifications from the sequence
         let sequence = pattern.replace_all(&raw_sequence, "").to_string();
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
index 588809a77..f67046c3e 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
@@ -1,9 +1,19 @@
+import os
+
 import numpy as np
 import pandas as pd
 from typing import Optional
 from collections import Counter
 from imspy_core.data.peptide import PeptideSequence
 
+try:  # imspy-connector >= 0.4.4. The submodule is an attribute of the extension module, not an
+    # importable path, so `from imspy_connector.py_peptide import ...` raises ModuleNotFoundError
+    # and would silently leave us on the slow path.
+    import imspy_connector as _ims_connector
+    _batch_masses = _ims_connector.py_peptide.mono_isotopic_masses
+except (ImportError, AttributeError):  # pragma: no cover - exercised only on an older connector
+    _batch_masses = None
+
 from imspy_predictors.rt.predictors import DeepChromatographyApex, load_deep_retention_time_predictor
 from imspy_predictors.ionization.predictors import predict_peptide_flyability_with_koina
 from imspy_predictors.utility import load_tokenizer_from_resources
@@ -72,6 +82,7 @@ def simulate_peptides(
         max_length: int = 30,
         proteome_mix: bool = False,
         use_koina_model: Optional[str] = None,
+        num_threads: int = -1,
 ) -> pd.DataFrame:
     """
     Simulate peptides from a protein table.
@@ -87,6 +98,7 @@ def simulate_peptides(
         max_length: Maximum length of the peptides.
         proteome_mix: If True, simulate a proteome mix.
         use_koina_model: If not None, use the Koina model to predict peptide flyability, i.e. chances of being measured in the experiment, currently only supports pfly.
+        num_threads: Threads for the batched monoisotopic mass calculation; -1 uses every core.
     Returns:
         DataFrame with simulated peptides.
     """
@@ -103,7 +115,6 @@ def simulate_peptides(
             if len(peptide) < min_length or len(peptide) > max_length:
                 continue
 
-            masses.append(PeptideSequence(peptide).mono_isotopic_mass)
             missed_cleavages.append(0)
             decoys.append(0)
             n_term.append(None)
@@ -115,6 +126,15 @@ def simulate_peptides(
             protein_id.append(row.protein_id)
             i += 1
 
+    # One batched, parallel call instead of constructing a PeptideSequence per peptide. Building
+    # that object per peptide cost ~190 us each — 48 s for 250 000 peptides — and was the largest
+    # remaining serial cost in the pipeline. Results are bitwise identical, including across thread
+    # counts. Falls back to the per-peptide path on a connector without the batch function.
+    if _batch_masses is not None:
+        masses = _batch_masses(sequences, num_threads if num_threads > 0 else (os.cpu_count() or 4))
+    else:
+        masses = [PeptideSequence(p).mono_isotopic_mass for p in sequences]
+
     peptide_table = pd.DataFrame({
         "protein_id": protein_id,
         "peptide_id": peptide_id,
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
index 2b3aaad1c..952901e67 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
@@ -1796,6 +1796,7 @@ def section(title: str) -> str:
                 min_length=config.min_len,
                 max_length=config.max_len,
                 proteome_mix=config.proteome_mix,
+                num_threads=config.num_threads,
             )
 
             stages.lap('sample_peptides+rt_filter')
diff --git a/packages/imspy-simulation/tests/test_batched_peptide_masses.py b/packages/imspy-simulation/tests/test_batched_peptide_masses.py
new file mode 100644
index 000000000..7c751d272
--- /dev/null
+++ b/packages/imspy-simulation/tests/test_batched_peptide_masses.py
@@ -0,0 +1,51 @@
+"""The batched monoisotopic mass path must match the per-peptide one exactly.
+
+Constructing a `PeptideSequence` per peptide cost ~190 us each, which made 250 000 peptides take
+48 s and was the largest remaining serial cost in the pipeline. Two changes: the lookup tables and
+the UniMod regex are built once rather than per call (9.4x on its own, and it helps every caller),
+and the connector gained a batched parallel entry point. Mass values must not move.
+"""
+import numpy as np
+import pytest
+
+from imspy_core.data.peptide import PeptideSequence
+
+_ims = pytest.importorskip("imspy_connector")
+mono_isotopic_masses = getattr(getattr(_ims, "py_peptide", None), "mono_isotopic_masses", None)
+if mono_isotopic_masses is None:  # pragma: no cover
+    pytest.skip("connector without mono_isotopic_masses", allow_module_level=True)
+
+
+def _peptides(n, seed=0):
+    rng = np.random.default_rng(seed)
+    aa = list("ACDEFGHIKLMNPQRSTVWY")
+    return ["".join(rng.choice(aa, int(rng.integers(7, 31)))) for _ in range(n)]
+
+
+@pytest.mark.parametrize("num_threads", [1, 4, 16])
+def test_batched_masses_match_the_per_peptide_path(num_threads):
+    peps = _peptides(2000)
+    ref = [PeptideSequence(p).mono_isotopic_mass for p in peps]
+    got = mono_isotopic_masses(peps, num_threads)
+    assert got == ref, "batched masses differ from the per-peptide path"
+
+
+def test_batched_masses_are_thread_count_invariant():
+    """Bitwise, not approximately: float addition is not associative, so a parallel sum that
+    varied with scheduling would make the simulation irreproducible."""
+    peps = _peptides(5000, seed=7)
+    assert mono_isotopic_masses(peps, 1) == mono_isotopic_masses(peps, 16)
+
+
+def test_modifications_are_still_counted():
+    """The UniMod regex moved into a static; modified peptides must still pick up their delta."""
+    plain, mod = "PEPTIDEK", "PEPTIDEM[UNIMOD:35]K"
+    m_plain, m_mod = mono_isotopic_masses([plain, mod], 1)
+    assert m_mod > m_plain
+    # UNIMOD:35 is oxidation, +15.994915, on top of the extra methionine residue.
+    assert abs((m_mod - m_plain) - (PeptideSequence("PEPTIDEMK").mono_isotopic_mass
+                                    - PeptideSequence(plain).mono_isotopic_mass) - 15.994915) < 1e-6
+
+
+def test_empty_batch():
+    assert mono_isotopic_masses([], 4) == []

From 2e9940d876632ce8c3f4bfee56a8c243d0e4f1fb Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Mon, 14 Sep 2026 05:26:19 +0200
Subject: [PATCH 13/14] perf(timsim): thin the peptide table before predicting
 retention times
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The digest produces `num_peptides_total` peptides regardless of the requested
complexity — every grid config digests 250 000 — and the neural retention-time
predictor ran over all of them so that the run could discard most of the output
immediately afterwards. A 5k run predicted 250 000 retention times to keep
5 000; the smoke config predicted 250 000 to keep 2 000.

`simulate_peptides` now takes the target count and thins the table before the
predictor rather than after it. The thinning carries a 1.3x margin because the
retention-time filter below removes part of the early-eluting population, and
without headroom the table could fall under what the caller asked for. That
filter is density-based — it thins the early bins down to the median bin count —
so it behaves the same on a uniformly thinned table as on the full one.

Runs at or above the digest target are untouched: the thinning only applies when
the table is larger than the request, so the 250k and 500k datasets take exactly
the path they did before.

End-to-end smoke, 16 threads, 2 000 peptides kept: the FASTA stage falls
43.3 -> 6.5 s, its sampling lap 40.2 -> 3.4 s, and total wall 72 -> 35 s.

This does change which peptides a small run selects, since the thinning now
happens before the retention-time filter rather than after it. Reproducibility
is unaffected: two runs of the same config still produce a byte-identical
analysis.tdf_bin and identical ground-truth tables.

Also fixes the accompanying log line, which subtracted the filtered length from
itself and so always reported "Excluded 0 peptides with low retention times".

Tests: 215 Python passed / 15 skipped.
---
 .../timsim/jobs/simulate_peptides.py          | 28 ++++++++++++++++++-
 .../src/imspy_simulation/timsim/simulator.py  |  2 ++
 .../tests/test_emg_degenerate_params.py       | 22 +++++++++++++++
 3 files changed, 51 insertions(+), 1 deletion(-)

diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
index f67046c3e..7cc166a3c 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/jobs/simulate_peptides.py
@@ -83,6 +83,9 @@ def simulate_peptides(
         proteome_mix: bool = False,
         use_koina_model: Optional[str] = None,
         num_threads: int = -1,
+        num_sample_peptides: Optional[int] = None,
+        sample_seed: Optional[int] = None,
+        sample_margin: float = 1.3,
 ) -> pd.DataFrame:
     """
     Simulate peptides from a protein table.
@@ -99,6 +102,12 @@ def simulate_peptides(
         proteome_mix: If True, simulate a proteome mix.
         use_koina_model: If not None, use the Koina model to predict peptide flyability, i.e. chances of being measured in the experiment, currently only supports pfly.
         num_threads: Threads for the batched monoisotopic mass calculation; -1 uses every core.
+        num_sample_peptides: If set, thin the table to roughly this many peptides *before* the
+            retention-time prediction, instead of predicting for the whole digest and discarding
+            afterwards. None keeps the old behaviour.
+        sample_seed: Seed for that thinning, so the choice is reproducible.
+        sample_margin: Oversample factor applied to `num_sample_peptides`, to leave headroom for
+            the retention-time filter below, which drops part of the early-eluting population.
     Returns:
         DataFrame with simulated peptides.
     """
@@ -148,6 +157,20 @@ def simulate_peptides(
         "events": events}
     )
 
+    # Thin the table here, before retention-time prediction, rather than after it.
+    #
+    # The caller asks for `num_sample_peptides` but the digest produces `num_peptides_total`, and
+    # the retention-time predictor used to run over the whole digest so that 98 % of its output
+    # could be thrown away downstream. A small run spent ~40 s predicting times it never used.
+    #
+    # `sample_margin` is headroom for the retention-time filter further down, which removes part of
+    # the early-eluting population; without it the table could fall below what the caller asked for.
+    # The filter itself is density-based (it thins the early bins to the median bin count), so it
+    # behaves the same on a uniformly thinned table as on the full one.
+    if num_sample_peptides is not None and len(peptide_table) > num_sample_peptides:
+        target = min(len(peptide_table), int(np.ceil(num_sample_peptides * max(sample_margin, 1.0))))
+        peptide_table = peptide_table.sample(n=target, random_state=sample_seed).reset_index(drop=True)
+
     # Simulate peptide efficiency using rejection sampling in log-normal space
     efficiency = generate_normal_efficiency(
         n=len(peptide_table),
@@ -222,9 +245,12 @@ def simulate_peptides(
             random_indices = np.array([], dtype=int)
         rt_filter[random_indices] = True
 
+        n_excluded = int((~rt_filter).sum())
         peptide_table = peptide_table[rt_filter]
 
         if verbose:
-            print(f"Excluded {len(peptide_table) - rt_filter.sum()} peptides with low retention times.")
+            # Count before the filter is applied: the old message subtracted the filtered length
+            # from itself and so always reported 0.
+            print(f"Excluded {n_excluded} peptides with low retention times.")
 
     return peptide_table
diff --git a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
index 952901e67..193a622ad 100644
--- a/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
+++ b/packages/imspy-simulation/src/imspy_simulation/timsim/simulator.py
@@ -1797,6 +1797,8 @@ def section(title: str) -> str:
                 max_length=config.max_len,
                 proteome_mix=config.proteome_mix,
                 num_threads=config.num_threads,
+                num_sample_peptides=config.num_sample_peptides if config.sample_peptides else None,
+                sample_seed=config.sample_seed,
             )
 
             stages.lap('sample_peptides+rt_filter')
diff --git a/packages/imspy-simulation/tests/test_emg_degenerate_params.py b/packages/imspy-simulation/tests/test_emg_degenerate_params.py
index 7442ee69e..55b1e11be 100644
--- a/packages/imspy-simulation/tests/test_emg_degenerate_params.py
+++ b/packages/imspy-simulation/tests/test_emg_degenerate_params.py
@@ -45,3 +45,25 @@ def test_sampled_parameters_are_strictly_positive(sampler, second):
     assert (others > 0).all(), f"{(others <= 0).sum()} {second} draws were non-positive"
     mu = estimate_mu_from_mode_emg(np.full_like(sigmas, 10.0), sigmas, others)
     assert np.isfinite(mu).all(), f"{(~np.isfinite(mu)).sum()} non-finite mu values"
+
+
+def test_thinning_happens_before_retention_time_prediction():
+    """The expensive predictor must see only the peptides we intend to keep.
+
+    The digest produces `num_peptides_total` regardless of the requested complexity, so predicting
+    retention times for the whole digest meant a small run spent most of its time on peptides it
+    then discarded. Thinning first must still deliver at least the requested count, since the
+    retention-time filter drops part of the early-eluting population afterwards.
+    """
+    import inspect
+    from imspy_simulation.timsim.jobs import simulate_peptides as mod
+
+    src = inspect.getsource(mod.simulate_peptides)
+    thin = src.index("peptide_table.sample(")
+    predict = src.index("simulate_separation_times_pandas")
+    assert thin < predict, "thinning must precede retention-time prediction"
+
+    sig = inspect.signature(mod.simulate_peptides).parameters
+    assert "num_sample_peptides" in sig and "sample_margin" in sig
+    # Headroom must be real, or the retention-time filter can leave us short of the request.
+    assert sig["sample_margin"].default > 1.0

From 85260d523f0f21a61e16e305a4f295f7c45054f3 Mon Sep 17 00:00:00 2001
From: theGreatHerrLebert <davidteschner@googlemail.com>
Date: Wed, 16 Sep 2026 14:39:26 +0200
Subject: [PATCH 14/14] test(imspy-simulation): drop the local fasta dependency
 from the digest tests

The ordering test read a fasta from a developer machine and had no skip guard, so it failed in
CI with FileNotFoundError. Ordering is a property of the digest, not of any particular proteome:
it now runs on four synthetic proteins defined in the module.

The parity test keeps the local HeLa subset when it is present and falls back to the same
synthetic proteins otherwise, so CI exercises the batched digest instead of skipping it.
---
 .../test_perf_batched_digest_and_overlay.py   | 33 +++++++++++++++----
 1 file changed, 26 insertions(+), 7 deletions(-)

diff --git a/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
index b600a203f..f4c46406e 100644
--- a/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
+++ b/packages/imspy-simulation/tests/test_perf_batched_digest_and_overlay.py
@@ -16,10 +16,9 @@
 STAT = {"C": "[UNIMOD:4]"}
 
 
-@pytest.mark.skipif(not os.path.exists(FASTA), reason="local HeLa subset fasta not available")
 @pytest.mark.parametrize("decoys", [False, True])
 def test_batched_digest_matches_per_protein(decoys):
-    tbl = parse_fasta_to_dataframe(FASTA).sample(n=400, random_state=7)
+    tbl = _proteins(n=400)
     batched = proteins_to_peptides_batched(
         tbl.index, tbl.sequence, generate_decoys=decoys, variable_mods=VAR, static_mods=STAT,
         chunk_size=150,  # force several chunks
@@ -38,24 +37,44 @@ def test_batched_digest_matches_per_protein(decoys):
             # Sorted, so the draw downstream cannot depend on set iteration order.
             assert got == sorted(got), f"protein {idx}: batched digest is not sorted"
             n_nonempty += 1
-    assert n_nonempty > 300
+    assert n_nonempty >= min(len(tbl) - 1, 300)
 
 
 def test_batched_digest_empty_input():
     assert proteins_to_peptides_batched([], []) == []
 
 
+# Synthetic proteins: several tryptic sites each, so every one yields multiple peptides.
+# Deliberately not the local HeLa fixture — the ordering guarantee is what is under test, and a
+# test that silently skips on a machine without that file proves nothing in CI.
+_SYNTH = {
+    "sp|TEST1|SYN1": "MKVLATRSEDKPEPTIDEKSAMPLERGGVTLLKDIFFERENTKAVQWYRPLNESTGGK",
+    "sp|TEST2|SYN2": "MSSGQKALTTRPEPTIDEKMNQWYFRAVCDEKLLGGTTRSEQPPLKYYTTRDDFGHIK",
+    "sp|TEST3|SYN3": "MTTKAAQPLRSEEDNKPEPTIDERGGWWYTKLNMFCVVRDDQQTTLKSSAHHPPYYRK",
+    "sp|TEST4|SYN4": "MNNKEFGHLRPPQTTSKVVCCWYKDDEEFFRAALLNMSTTKGGHHYYPPQRWWEEDDK",
+}
+
+
+def _proteins(n: int) -> pd.DataFrame:
+    """Protein table for the digest tests: the local HeLa subset when present, else the synthetic
+    set above. CI has no local fixtures, and a skipped parity test would leave the batched digest
+    unchecked there."""
+    if os.path.exists(FASTA):
+        return parse_fasta_to_dataframe(FASTA).sample(n=n, random_state=7)
+    return pd.DataFrame({"sequence": list(_SYNTH.values())}, index=list(_SYNTH))
+
+
 def test_digest_output_is_order_stable_not_a_set():
     """Set iteration order varies with Python's per-process hash seed; the sampler draws from
     these collections, so an unordered container makes a seeded run irreproducible."""
-    tbl = parse_fasta_to_dataframe(FASTA).sample(n=40, random_state=3)
-    batched = proteins_to_peptides_batched(tbl.index, tbl.sequence, variable_mods=VAR, static_mods=STAT)
+    ids, seqs = list(_SYNTH), list(_SYNTH.values())
+    batched = proteins_to_peptides_batched(ids, seqs, variable_mods=VAR, static_mods=STAT)
     assert any(len(p) > 1 for p in batched), "fixture produced nothing to order"
     for peptides in batched:
         assert isinstance(peptides, list), "digest must return an ordered container"
         assert peptides == sorted(peptides)
-    for (idx, row) in tbl.iterrows():
-        got = protein_to_peptides(generate_single_fasta(idx, row.sequence), generate_decoys=False,
+    for idx, seq in zip(ids, seqs):
+        got = protein_to_peptides(generate_single_fasta(idx, seq), generate_decoys=False,
                                   variable_mods=VAR, static_mods=STAT, cleave_at='KR', restrict='P',
                                   missed_cleavages=2, min_len=7, max_len=30, digest=True)
         if got: