diff --git a/CHANGELOG.md b/CHANGELOG.md index be097b9..9ae66cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -177,6 +177,63 @@ than a number. Both are recorded in every run's `manifest.json`. cancellation flag was written and never read. The waiter is now the only writer: it reads the intent after reaping the engine and publishes `cancelled`, `done` when the engine had already finished, or `failed`; until then the run shows "Stopping" (#14). +- Code review F, whole-repository review (`docs/31_code_review_2026-09-08_full.md`, + F1 to F10): + - `prescan` read the infinite-bounds sentinel that `rt-im-train` writes for "calibration + unavailable, search the whole gradient" as "cannot be screened" and dropped the + candidate, so a run with no confident seeds discarded the entire library and exited 0 + with a zero-row survivors table. An unbounded window now screens over the whole + gradient, a candidate with no window row is treated the same, both are counted, and + screening every candidate away is an error (F1). + - A present-but-wrong-typed `is_transferred` was swallowed as "no transfers", silently + removing every match-between-runs identification from `peptides.tsv` and + `proteins.tsv` while the parquet still carried them. Present columns are read in their + declared type and a mismatch is an error; only an absent column falls back (F2). + - `sidecar::resolve_script` tried the working directory before the directory beside the + binary, the ordering `python::resolve_script_dir` was hardened against, so a `scripts/` + directory inside an untrusted dataset could have its worker executed. An absolute + directory is taken as given, then the executable's directory, then `/scripts`, and + the working directory last (F3). + - `Loess::predict` indexed before the start of its grid for a non-finite query, so one + library row with a null `predicted_irt` could abort or misread memory at rt-im-train. + It returns NaN, and rt-im-train treats a non-finite library iRT as "no calibrated RT" + and counts those rows (F4). + - The rescorer's in-memory TSV backend standardised with median/IQR while the parquet and + streaming backends used mean/std, so the same pool scored differently depending on + `rescore.handoff` and on the 4 GB streaming threshold. All three use mean/std, which + leaves the shipped parquet default and every published benchmark unchanged. The + `MUMDIA_NN_FOLD_KEYS` companion is length-checked instead of being sliced short, which + used to leave the tail rows unscored at a fabricated mid-rank score, and the estimate + that picks the backend counts feature columns by name (F5). + - `refuse_output_over_input` was wired into two of eighteen stages, so + `compete --features f.parquet --out f.parquet` replaced the widest artifact of the run + with the competed subset at exit 0. It now guards every output of `search-seed`, + `rt-im-train`, `extract`, `features`, `compete`, `rescore`, `quant` and `audit` (F6). + - The LOESS boundary extrapolation slope introduced in the previous package was the + pointwise local slope at the sparsest, most one-sided point of the fit: unbounded, free + to be negative, and multiplying an unbounded distance. It is the secant of the fitted + curve over its end decile, clamped non-negative and to at most four times the global + slope, and the test uses noisy anchors rather than a noiseless quadratic (F7). + Measured on HYE B01 against the previous behaviour, same library and settings: 48,533 + stripped peptides at 1%, 53,127 PSM-q 1% targets, 6,519 protein groups and 1,961,800 + extracted rows in both arms, identical to the row. The two extrapolations agree + wherever the anchors are dense and differ only outside the anchor range. + - A desktop stop arriving between the reap and the end of `publish_exit` could pass a + recycled process id to the tree kill. The waiter retires the id the instant `wait` + returns, before it reads the output directory (F8). + - The conversion lock added in the previous package spun without pause on an undeletable + stale lock, mistook clock skew and a peer's partial file for evidence about its own + holder, could be held by two processes at once, and left every interrupted conversion's + partial mzML behind for ever. Take-overs are bounded and paced, the holder is + identified by a token it reads back, a future modification time counts as fresh, the + partial-file probe matches this destination only, and abandoned partials are swept + under the lock (F9). + - Dropping an unpredicted candidate with everything sharing its pair key also removed + positional isomers that predicted correctly, bounded only by the library being emptied. + The direct misses and the collateral are counted separately and exceeding 2% of the + library is an error naming the sidecar. The key stays position-free deliberately: a + positional key would stop matching a reverse decoy to its target, trading a sensitivity + defect for an FDR one (F10). - Code review E, follow-up (`docs/30_code_review_2026-09-08.md`, R1 to R9): - Enabling DeepLC fine-tuning with its own defaults was rejected at load, because the documented automatic batch size is `finetune_batch = 0` and the new validation demanded diff --git a/desktop/src-tauri/src/run.rs b/desktop/src-tauri/src/run.rs index 67f9129..991d582 100644 --- a/desktop/src-tauri/src/run.rs +++ b/desktop/src-tauri/src/run.rs @@ -275,6 +275,22 @@ impl Run { matches!(self.snapshot().status.as_str(), "running" | "starting") } + /// Forget the process id, under the lock a stop takes before it kills. + /// + /// Called the instant `wait` returns, before anything else. The pid is free for the + /// operating system to reuse from that moment, and a stop landing later would + /// otherwise pass it to `kill_tree`, which on Windows terminates whatever now owns it + /// and its whole tree (docs/31 F8). Retiring it inside `publish_exit` was too late: + /// that function scans the output directory and reads the result reports first, so + /// the window was as long as that disk work. It is now a few instructions, and a stop + /// that reaches the lock inside it still finds the pid this run really owns, because + /// the reap has only just returned. + fn retire_pid(&self) { + if let Ok(mut p) = self.pid.lock() { + *p = None; + } + } + /// Apply `f` only while the run is still active; returns whether it was. fn set_if_active(&self, f: F) -> bool { if let Ok(mut s) = self.snapshot.lock() { @@ -298,12 +314,11 @@ impl Run { /// engine finished before the kill landed and its outputs are complete, and /// calling them cancelled would hide a finished result. fn publish_exit(&self, outcome: std::io::Result, out_dir: &Path) { - // Retire the pid first. This waits for a stop that is still killing (it holds the - // same lock), so nothing below overlaps a kill, and a later stop finds nothing to - // signal (docs/30 R3). - if let Ok(mut p) = self.pid.lock() { - *p = None; - } + // Idempotent: the waiter retires the pid the moment `wait` returns (docs/31 F8), + // and this call is what makes `publish_exit` safe to reach from a test or any + // other path. Taking the lock here also waits for a stop that is still killing, + // so nothing below overlaps a kill. + self.retire_pid(); let cancelled = self.cancelled.load(Ordering::SeqCst); if cancelled { // The only sweep: after the reap, before the release, inside this run's @@ -748,6 +763,8 @@ pub fn start(id: String, req: Request) -> Result, String> { let out_dir = PathBuf::from(&req.out_dir); std::thread::spawn(move || { let outcome = child.wait(); + // Before anything else: the pid is reusable from here (docs/31 F8). + run.retire_pid(); run.publish_exit(outcome, &out_dir); }); } @@ -988,6 +1005,22 @@ mod tests { let _ = std::fs::remove_dir_all(&dir); } + #[test] + fn a_stop_after_the_reap_has_no_pid_to_kill() { + // docs/31 F8: the pid is reusable the moment `wait` returns, so the waiter retires + // it there rather than after reading the output directory. A stop arriving in + // between must find nothing, not a recycled pid. + let (run, dir) = running("reaped"); + *run.pid.lock().unwrap() = Some(4242); + run.retire_pid(); + assert_eq!(*run.pid.lock().unwrap(), None); + // The run is still active, so cancel proceeds and simply has nothing to signal. + run.cancel(); + assert!(run.snapshot().cancel_requested); + assert_eq!(*run.pid.lock().unwrap(), None); + let _ = std::fs::remove_dir_all(&dir); + } + #[test] fn cancel_records_the_intent_without_publishing_a_terminal_status() { // Until the engine is reaped the run is still running, whatever the button diff --git a/docs/06_predict_frag_index_matchers.md b/docs/06_predict_frag_index_matchers.md index 3701a1f..17c721d 100644 --- a/docs/06_predict_frag_index_matchers.md +++ b/docs/06_predict_frag_index_matchers.md @@ -429,7 +429,7 @@ m/z (`Library::local_frag_index`, `index.rs:325`). | `RtPredictor` / `FragmentPredictor` | `predict.rs:13` / `predict.rs:19` | predictor traits (predict + `identity`); implemented only by the native structs | | `NativeRt` | `predict.rs:25` | additive retention-coefficient model + `sqrt(len)` + `0.01*mod` term, `identity` `native-rt-v1` | | `NativeFrag` | `predict.rs:73` | heuristic b/y intensity model (y=1.0, b=0.75, mid-seq positional, charge-2 x0.5), max-normalized, `identity` `native-frag-v1` | -| `resolve_script` | `sidecar.rs:20` | locate a worker script (CWD, exe dir/dir, exe dir/scripts, else CWD-relative) | +| `resolve_script` | `sidecar.rs:20` | locate a worker script (an absolute directory as given, else exe dir/dir, exe dir/scripts, and the working directory LAST; docs/31 F3) | | `run_ms2pip` | `sidecar.rs:42` | MS2PIP client; in `id`/`peptidoform`/`charge`, out `id`/`ion_type`/`ordinal`/`intensity`; returns `cid -> (ion_byte, ordinal) -> intensity` | | `run_deeplc` | `sidecar.rs:81` | DeepLC client; in `id`/`peptidoform`, out `id`/`predicted_rt`; returns `id -> predicted_rt` | | `run_deeplc_finetune` | `sidecar.rs:111` | DeepLC multitask fine-tune; `deeplc_finetune.py ` + epoch/patience/q-train/batch flags (called by `run`, not predict-frag) | diff --git a/docs/08_rt_im_train.md b/docs/08_rt_im_train.md index 4fa1227..9f02f7c 100644 --- a/docs/08_rt_im_train.md +++ b/docs/08_rt_im_train.md @@ -208,8 +208,10 @@ The math: contribute. If the weighted system is degenerate, meaning `sw < 1e-12` (all weights vanished) or `sw*swxx - swx^2` near zero (calibrate.rs:146-149), it falls back to the global line. -- **`Loess::predict`** (calibrate.rs:87-102) uses the linear fallback for `x` at or - outside the grid ends, and also when the grid has fewer than 2 nodes +- **`Loess::predict`** (calibrate.rs:87-102) returns NaN for a non-finite `x` (a null + library iRT reads as NaN, and interpolating it used to index before the start of the + grid; docs/31 F4), extrapolates from the nearer grid end for `x` outside the range, + uses the linear fallback when the grid has fewer than 2 nodes (calibrate.rs:89-94), and otherwise linearly interpolates between the two bracketing grid nodes found by `partition_point` (calibrate.rs:95-101). When the two bracketing nodes are within `1e-12` in `x` it returns the lower node's `y` @@ -514,9 +516,14 @@ cheaper option with equal RT residuals on the one run measured here. | `candidate_window` | rt_im_train.rs:65-70 | Builds `(cal, lo, hi)`; returns `(NaN, -inf, +inf)` when calibrated RT or width is absent. | | `rt_im_train::run` | rt_im_train.rs:72-354 | The stage: join iRT, select anchors, fit, window, apply, write. | | `linear_fit` | calibrate.rs:6-28 | OLS `y = slope*x + intercept` with degenerate-case guards. | -| `Loess` | calibrate.rs:31-37 | Grid-based local-linear smoother; carries the boundary local slopes for extrapolation and the global line as the degenerate fallback. | +| `Loess` | calibrate.rs:31-37 | Grid-based local-linear smoother; carries the two boundary extrapolation slopes and the global line as the degenerate fallback. | | `Loess::fit` | calibrate.rs:42-83 | Sorts anchors, builds a `grid_n`-point local-linear grid, `k = clamp(ceil(span*n),3,n)`. | -| `Loess::predict` | calibrate.rs:87-102 | Grid interpolation inside range; outside it, the boundary grid value continued with the boundary local slope, so the map is continuous at both ends. Until docs/29 #10 it switched to the global line there, which on `y = 200 + 10x^2` (span 0.3) jumped from 193.4 to 38.3 at `x = 0` and from 1173.5 to 1018.4 at the top: about 155 s discontinuities for gradient-edge peptides. Measured on HYE B01 with the imported iRT and `native_tda`: 45,946 stripped peptides at 1% before, 45,957 after, decoy fraction unchanged; 1.9% of candidates got a different window, almost all with iRT above the anchor range, which the global line had placed past the end of the run. With the DeepLC 4.1.1 re-predicted precursors (`w_rt` 414 s): 48,533 in both arms, 0.2% of windows moved. | +| `Loess::predict` | calibrate.rs:87-102 | NaN for a non-finite query (docs/31 F4). Grid interpolation inside range; outside it, the boundary grid value continued with the boundary extrapolation slope, so the map is continuous at both ends. Until docs/29 #10 it switched to the global line there, which on `y = 200 + 10x^2` (span 0.3) jumped from 193.4 to 38.3 at `x = 0` and from 1173.5 to 1018.4 at the top: about 155 s discontinuities for gradient-edge peptides. The extrapolation slope is the +secant of the fitted curve over its end decile, clamped non-negative and to at most four +times the global slope: the pointwise local slope it first used comes from the sparsest, +most one-sided window in the fit, and on noisy anchors it was free to be negative (which +inverts the iRT-to-RT map) or several times the global slope, multiplying a distance that +is unbounded by construction (docs/31 F7). Measured on HYE B01 with the imported iRT and `native_tda`: 45,946 stripped peptides at 1% before, 45,957 after, decoy fraction unchanged; 1.9% of candidates got a different window, almost all with iRT above the anchor range, which the global line had placed past the end of the run. With the DeepLC 4.1.1 re-predicted precursors (`w_rt` 414 s): 48,533 in both arms, 0.2% of windows moved. | | `local_linear` | calibrate.rs:107-153 | Tricubic-weighted local least squares at one point. | | `percentile` | calibrate.rs:156-164 | Nearest-rank percentile: sorts a copy, `rank = round(p.clamp(0,1)*(len-1))`. Not interpolated. Empty input returns 0.0. | | `CalibrationMethod` | config.rs:56-61 | Enum `{ Loess, Linear, None }`; default `Loess`. `None` is rejected at load. | @@ -621,9 +628,9 @@ though the enum variant still exists. - **`slope`/`intercept` are emitted in `cal.json` whenever calibration is available, including under LOESS** (rt_im_train.rs:128, 286-287, 313-314). Under LOESS they are the degenerate fallback (fewer than four anchors, or a local window without - spread), not the extrapolation model: since docs/29 #10 the map continues the - boundary local fit outside the anchor range. Do not read them as the calibration - when `method == "loess"`. They are serialized as `null` only when calibration is + spread), not the extrapolation model: since docs/29 #10 the map continues the fitted + curve outside the anchor range. Do not read them as the calibration when + `method == "loess"`. They are serialized as `null` only when calibration is unavailable (`n_train < 2`), since the fit is not computed in that case. - **`CalibrationMethod::None` still exists but is rejected** at config load (config.rs:1336-1342). The stage would otherwise fall through to the linear path diff --git a/docs/24_config_reference.md b/docs/24_config_reference.md index 10899bf..4bcbfb1 100644 --- a/docs/24_config_reference.md +++ b/docs/24_config_reference.md @@ -713,49 +713,49 @@ listed with the file it is in. | `MUMDIA_MOKAPOT_WORKERS` | sidecar | `"3"` | `scripts/mokapot_worker.py:120` | | `MUMDIA_MSCONVERT` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/raw.rs:193` | | `MUMDIA_NN_ALPHA` | sidecar | `"1e-4"` | `scripts/mokapot_worker.py:90` | -| `MUMDIA_NN_BATCH` | sidecar | `4096` | `scripts/nn_rescore_worker.py:280` | -| `MUMDIA_NN_CHUNK` | sidecar | `250000` | `scripts/nn_rescore_worker.py:284` | -| `MUMDIA_NN_DEVICE` | sidecar | `"auto"` | `scripts/nn_rescore_worker.py:292` | -| `MUMDIA_NN_DROPOUT` | sidecar | `0.3` | `scripts/nn_rescore_worker.py:277` | -| `MUMDIA_NN_EARLY_STOP` | sidecar | `1` | `scripts/nn_rescore_worker.py:285` | -| `MUMDIA_NN_EARLY_STOP_TOL` | sidecar | `0.01` | `scripts/nn_rescore_worker.py:286` | -| `MUMDIA_NN_EPOCHS` | sidecar | `25` | `scripts/nn_rescore_worker.py:275` | -| `MUMDIA_NN_FEATURES` | sidecar | `""` | `scripts/nn_rescore_worker.py:368` | -| `MUMDIA_NN_FOLDS` | sidecar | `3` | `scripts/nn_rescore_worker.py:263` | -| `MUMDIA_NN_FOLD_KEYS` | sidecar | `""` | `scripts/nn_rescore_worker.py:407` | -| `MUMDIA_NN_HIDDEN` | sidecar | `"128,64"` in nn_rescore_worker.py; `"128,64,64,32"` in mokapot_worker.py | `scripts/mokapot_worker.py:83`, `scripts/nn_rescore_worker.py:276` | -| `MUMDIA_NN_INIT_SAMPLE` | sidecar | `300000` | `scripts/nn_rescore_worker.py:537` | -| `MUMDIA_NN_INIT_TOPK` | sidecar | `0` | `scripts/nn_rescore_worker.py:659` | -| `MUMDIA_NN_ITERS` | sidecar | `5` | `scripts/nn_rescore_worker.py:274` | -| `MUMDIA_NN_LR` | sidecar | `1e-3` | `scripts/nn_rescore_worker.py:278` | -| `MUMDIA_NN_MARGIN_FRAC` | sidecar | `0.5` | `scripts/nn_rescore_worker.py:273` | +| `MUMDIA_NN_BATCH` | sidecar | `4096` | `scripts/nn_rescore_worker.py:296` | +| `MUMDIA_NN_CHUNK` | sidecar | `250000` | `scripts/nn_rescore_worker.py:300` | +| `MUMDIA_NN_DEVICE` | sidecar | `"auto"` | `scripts/nn_rescore_worker.py:308` | +| `MUMDIA_NN_DROPOUT` | sidecar | `0.3` | `scripts/nn_rescore_worker.py:293` | +| `MUMDIA_NN_EARLY_STOP` | sidecar | `1` | `scripts/nn_rescore_worker.py:301` | +| `MUMDIA_NN_EARLY_STOP_TOL` | sidecar | `0.01` | `scripts/nn_rescore_worker.py:302` | +| `MUMDIA_NN_EPOCHS` | sidecar | `25` | `scripts/nn_rescore_worker.py:291` | +| `MUMDIA_NN_FEATURES` | sidecar | `""` | `scripts/nn_rescore_worker.py:387` | +| `MUMDIA_NN_FOLDS` | sidecar | `3` | `scripts/nn_rescore_worker.py:279` | +| `MUMDIA_NN_FOLD_KEYS` | sidecar | `""` | `scripts/nn_rescore_worker.py:426` | +| `MUMDIA_NN_HIDDEN` | sidecar | `"128,64"` in nn_rescore_worker.py; `"128,64,64,32"` in mokapot_worker.py | `scripts/mokapot_worker.py:83`, `scripts/nn_rescore_worker.py:292` | +| `MUMDIA_NN_INIT_SAMPLE` | sidecar | `300000` | `scripts/nn_rescore_worker.py:564` | +| `MUMDIA_NN_INIT_TOPK` | sidecar | `0` | `scripts/nn_rescore_worker.py:686` | +| `MUMDIA_NN_ITERS` | sidecar | `5` | `scripts/nn_rescore_worker.py:290` | +| `MUMDIA_NN_LR` | sidecar | `1e-3` | `scripts/nn_rescore_worker.py:294` | +| `MUMDIA_NN_MARGIN_FRAC` | sidecar | `0.5` | `scripts/nn_rescore_worker.py:289` | | `MUMDIA_NN_MAX_ITER` | sidecar | `"200"` | `scripts/mokapot_worker.py:91` | -| `MUMDIA_NN_NEG_RATIO` | sidecar | `0.0` | `scripts/nn_rescore_worker.py:267` | -| `MUMDIA_NN_NEG_SELECT` | sidecar | `"random"` | `scripts/nn_rescore_worker.py:268` | -| `MUMDIA_NN_PREGATHER_GB` | sidecar | `8` | `scripts/nn_rescore_worker.py:287` | -| `MUMDIA_NN_SEED` | sidecar | `0` | `scripts/nn_rescore_worker.py:283` | -| `MUMDIA_NN_SEEDS` | sidecar | `1` | `scripts/nn_rescore_worker.py:282` | +| `MUMDIA_NN_NEG_RATIO` | sidecar | `0.0` | `scripts/nn_rescore_worker.py:283` | +| `MUMDIA_NN_NEG_SELECT` | sidecar | `"random"` | `scripts/nn_rescore_worker.py:284` | +| `MUMDIA_NN_PREGATHER_GB` | sidecar | `8` | `scripts/nn_rescore_worker.py:303` | +| `MUMDIA_NN_SEED` | sidecar | `0` | `scripts/nn_rescore_worker.py:299` | +| `MUMDIA_NN_SEEDS` | sidecar | `1` | `scripts/nn_rescore_worker.py:298` | | `MUMDIA_NN_SOLVER` | sidecar | `"adam"` | `scripts/mokapot_worker.py:89` | -| `MUMDIA_NN_STREAM` | sidecar | `"auto"` | `scripts/nn_rescore_worker.py:329` | -| `MUMDIA_NN_STREAM_GB` | sidecar | `4` | `scripts/nn_rescore_worker.py:338` | -| `MUMDIA_NN_THREADS` | both | `16` | `rust/mumdia/crates/mumdia/src/main.rs:93`, `scripts/nn_rescore_worker.py:313`, `scripts/nn_rescore_worker.py:314` | -| `MUMDIA_NN_TRAIN_FDR` | sidecar | `0.01` | `scripts/nn_rescore_worker.py:281` | -| `MUMDIA_NN_TRAIN_SUB` | sidecar | `0.0` | `scripts/nn_rescore_worker.py:264` | -| `MUMDIA_NN_WARM_EPOCHS` | sidecar | `0` | `scripts/nn_rescore_worker.py:266` | -| `MUMDIA_NN_WARM_START` | sidecar | `0` | `scripts/nn_rescore_worker.py:265` | -| `MUMDIA_NN_WD` | sidecar | `1e-4` | `scripts/nn_rescore_worker.py:279` | +| `MUMDIA_NN_STREAM` | sidecar | `"auto"` | `scripts/nn_rescore_worker.py:345` | +| `MUMDIA_NN_STREAM_GB` | sidecar | `4` | `scripts/nn_rescore_worker.py:357` | +| `MUMDIA_NN_THREADS` | both | `16` | `rust/mumdia/crates/mumdia/src/main.rs:93`, `scripts/nn_rescore_worker.py:329`, `scripts/nn_rescore_worker.py:330` | +| `MUMDIA_NN_TRAIN_FDR` | sidecar | `0.01` | `scripts/nn_rescore_worker.py:297` | +| `MUMDIA_NN_TRAIN_SUB` | sidecar | `0.0` | `scripts/nn_rescore_worker.py:280` | +| `MUMDIA_NN_WARM_EPOCHS` | sidecar | `0` | `scripts/nn_rescore_worker.py:282` | +| `MUMDIA_NN_WARM_START` | sidecar | `0` | `scripts/nn_rescore_worker.py:281` | +| `MUMDIA_NN_WD` | sidecar | `1e-4` | `scripts/nn_rescore_worker.py:295` | | `MUMDIA_PYTHON` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/python.rs:198` | | `MUMDIA_PYTHON_DEEPLC` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/python.rs:197` | | `MUMDIA_PYTHON_MBR` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/python.rs:197` | | `MUMDIA_PYTHON_MS2PIP` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/python.rs:197` | | `MUMDIA_PYTHON_RESCORE` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/python.rs:197` | -| `MUMDIA_RESCORE_MODEL` | both | `"nn"` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:345`, `scripts/mokapot_worker.py:181`, `scripts/mokapot_worker.py:37` | +| `MUMDIA_RESCORE_MODEL` | both | `"nn"` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:353`, `scripts/mokapot_worker.py:181`, `scripts/mokapot_worker.py:37` | | `MUMDIA_THERMO_PARSER` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/raw.rs:171` | | `MUMDIA_XGB_DEPTH` | sidecar | `"6"` | `scripts/mokapot_worker.py:63` | | `MUMDIA_XGB_JOBS` | sidecar | `"0"` | `scripts/mokapot_worker.py:68` | | `MUMDIA_XGB_LR` | sidecar | `"0.1"` | `scripts/mokapot_worker.py:64` | | `MUMDIA_XGB_TREES` | sidecar | `"200"` | `scripts/mokapot_worker.py:62` | -| `OMP_NUM_THREADS` | both | `16` | `rust/mumdia/crates/mumdia/src/main.rs:93`, `scripts/nn_rescore_worker.py:316`, `scripts/nn_rescore_worker.py:317` | +| `OMP_NUM_THREADS` | both | `16` | `rust/mumdia/crates/mumdia/src/main.rs:93`, `scripts/nn_rescore_worker.py:332`, `scripts/nn_rescore_worker.py:333` | | `PATH` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/raw.rs:293` | | `ProgramFiles` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/raw.rs:201` | | `ProgramFiles(x86)` | engine | none (unset means off) | `rust/mumdia/crates/mumdia/src/raw.rs:202` | @@ -774,23 +774,23 @@ one exception noted in its own help text: it sets `MUMDIA_NN_THREADS` and |---|---|---|---| | `KMP_DUPLICATE_LIB_OK` | sidecar | `"TRUE"` | `scripts/deeplc_finetune.py:28` | | `MKL_NUM_THREADS` | sidecar | `"1"` | `scripts/deeplc_finetune.py:31` | -| `MUMDIA_NN_FOLDS` | engine | `p.cfg.folds.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1347` | -| `MUMDIA_NN_FOLD_KEYS` | engine | `&foldkeys` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1374` | -| `MUMDIA_NN_ITERS` | engine | `p.cfg.num_iter.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1348` | -| `MUMDIA_NN_MARGIN_FRAC` | engine | `p.cfg.train_margin_frac.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1372` | -| `MUMDIA_NN_NEG_RATIO` | engine | `p.cfg.train_neg_ratio.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1353` | -| `MUMDIA_NN_NEG_SELECT` | engine | `match p.cfg.train_neg_select { mumdia_core::config::NegSelect::Random => "random", mumdia_core::config::NegSelect::Margin => "margin", mumdia_core::config::NegSelect::Hybrid => "hybrid", }` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1354` | -| `MUMDIA_NN_SEEDS` | engine | `p.cfg.seeds.max(1).to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1373` | +| `MUMDIA_NN_FOLDS` | engine | `p.cfg.folds.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1355` | +| `MUMDIA_NN_FOLD_KEYS` | engine | `&foldkeys` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1382` | +| `MUMDIA_NN_ITERS` | engine | `p.cfg.num_iter.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1356` | +| `MUMDIA_NN_MARGIN_FRAC` | engine | `p.cfg.train_margin_frac.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1380` | +| `MUMDIA_NN_NEG_RATIO` | engine | `p.cfg.train_neg_ratio.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1361` | +| `MUMDIA_NN_NEG_SELECT` | engine | `match p.cfg.train_neg_select { mumdia_core::config::NegSelect::Random => "random", mumdia_core::config::NegSelect::Margin => "margin", mumdia_core::config::NegSelect::Hybrid => "hybrid", }` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1362` | +| `MUMDIA_NN_SEEDS` | engine | `p.cfg.seeds.max(1).to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1381` | | `MUMDIA_NN_THREADS` | engine | `n.to_string()` | `rust/mumdia/crates/mumdia/src/main.rs:94` | -| `MUMDIA_NN_TRAIN_FDR` | engine | `p.cfg.train_fdr.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1349` | -| `MUMDIA_NN_TRAIN_SUB` | engine | `p.cfg.train_subsample.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1362` | -| `MUMDIA_NN_WARM_EPOCHS` | engine | `p.cfg.train_warm_epochs.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1371` | -| `MUMDIA_NN_WARM_START` | engine | `if p.cfg.train_warm_epochs > 0 { "1" } else { "0" }` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1363` | +| `MUMDIA_NN_TRAIN_FDR` | engine | `p.cfg.train_fdr.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1357` | +| `MUMDIA_NN_TRAIN_SUB` | engine | `p.cfg.train_subsample.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1370` | +| `MUMDIA_NN_WARM_EPOCHS` | engine | `p.cfg.train_warm_epochs.to_string()` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1379` | +| `MUMDIA_NN_WARM_START` | engine | `if p.cfg.train_warm_epochs > 0 { "1" } else { "0" }` | `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1371` | | `NUMEXPR_NUM_THREADS` | sidecar | `"1"` | `scripts/deeplc_finetune.py:32` | | `OMP_NUM_THREADS` | both | `"1"` in deeplc_finetune.py; `n.to_string()` in main.rs | `rust/mumdia/crates/mumdia/src/main.rs:94`, `scripts/deeplc_finetune.py:29` | | `OPENBLAS_NUM_THREADS` | sidecar | `"1"` | `scripts/deeplc_finetune.py:30` | -| `PYTHONIOENCODING` | engine | `"utf-8"` | `rust/mumdia/crates/mumdia/src/sidecar.rs:409` | -| `PYTHONUTF8` | engine | `"1"` | `rust/mumdia/crates/mumdia/src/sidecar.rs:409`, `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1103`, `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1341` | +| `PYTHONIOENCODING` | engine | `"utf-8"` | `rust/mumdia/crates/mumdia/src/sidecar.rs:425` | +| `PYTHONUTF8` | engine | `"1"` | `rust/mumdia/crates/mumdia/src/sidecar.rs:425`, `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1111`, `rust/mumdia/crates/mumdia/src/stages/rescore.rs:1349` | ## Unresolved by the generator @@ -806,8 +806,8 @@ Every field whose struct has an `impl Default` resolved from the source. 2 environment read(s) whose name is not a literal: -- `rust/mumdia/crates/mumdia/src/stages/extract.rs:2685: env read via closure of `&mut flushed`` -- `rust/mumdia/crates/mumdia/src/stages/extract.rs:2703: env read via closure of `&mut cand_hits`` +- `rust/mumdia/crates/mumdia/src/stages/extract.rs:2693: env read via closure of `&mut flushed`` +- `rust/mumdia/crates/mumdia/src/stages/extract.rs:2711: env read via closure of `&mut cand_hits`` ## Coverage diff --git a/docs/31_code_review_2026-09-08_full.md b/docs/31_code_review_2026-09-08_full.md new file mode 100644 index 0000000..d5ed834 --- /dev/null +++ b/docs/31_code_review_2026-09-08_full.md @@ -0,0 +1,201 @@ +# MuMDIA full code review — 2026-09-08 + +Baseline: `21b8ac7`, the head of `review/e-followup` (work packages A to E). Whole +repository rather than a diff: engine, core, IO, the Python workers and library helpers, +the desktop application, the tests, the benchmark utilities and the delivery +configuration. Twelve parallel finders and one verification sweep; the maintainer +re-verified eleven findings against the source before recording them here. + +Twenty-one findings. Six are confirmed correctness or security defects in code that +predates this review series, five are regressions introduced by packages A, D and E, and +ten are second-tier robustness, coverage and delivery problems. Package F below closes +the first eleven. + +**Verdict at the time of review: do not tag a release from this tree.** Two of the +findings produce a plausible wrong answer rather than an error, and one is a code-execution +path through an untrusted working directory. The 2026-09-08 HYE release check passing does +not contradict that: it exercised the library path with calibration available and no +match-between-runs transfer, so it could not reach any of them. + +## Findings + +Priority P1 means fix before relying on the affected workflow. *Reproduced* means an +executed probe demonstrated the failure; *source-confirmed* means the trigger and the +consequence follow from the implementation without a full run. + +| | Priority | Finding | Origin | +|---|---|---|---| +| F1 | P1 | prescan reads the no-calibration sentinel as "cannot be screened" and discards the entire library at exit 0 | pre-existing | +| F2 | P1 | a malformed `is_transferred` silently deletes every transfer from both TSV reports | pre-existing | +| F3 | P1 | `sidecar::resolve_script` tries the working directory first, executing a planted worker | pre-existing | +| F4 | P1 | `Loess::predict` indexes out of bounds when the query is NaN | pre-existing | +| F5 | P1 | three rescorer backends, two standardisations, and an unchecked fold-key companion | pre-existing | +| F6 | P1 | the output-over-input guard is wired into 2 of 18 stages | pre-existing | +| F7 | P2 | the LOESS boundary extrapolation slope is unbounded and may be negative | package D (#65) | +| F8 | P2 | a desktop stop can kill a recycled process id after the reap | package E (#66) | +| F9 | P2 | the conversion lock spins, mis-detects staleness three ways, and leaks partial files | package E (#66) | +| F10 | P2 | a position-free pair key deletes positional isomers that predicted correctly | package A (#62) | +| F11 | P2 | ten further robustness, coverage and delivery findings, listed at the end | mixed | + +**F1. prescan discards the library on the documented no-calibration path. Reproduced.** + +`rt_im_train::candidate_window` returns `(NaN, -inf, +inf)` whenever calibration is +unavailable and its doc comment calls the infinite bounds recall-safe. prescan's guard +rejected exactly that: `if !lo.is_finite() || !hi.is_finite() { return None; }`. This is the +FASTA/MS2PIP failure CLAUDE.md already records, where the seed search finds no confident +PSMs and every window row is infinite; prescan then wrote a zero-row survivors table with a +NaN target/decoy ratio and exited 0. The single-label bail could not see it, being gated on +`surv.len() > 1000`. A candidate simply absent from `run_windows` was dropped the same way, +uncounted. + +**F2. A malformed transfer flag deletes every transfer from the reports. Reproduced.** + +`mbr_worker.py` writes `is_transferred` through pandas, so a nullable boolean dtype, a null +or int8 0/1 reaches the report stage. `match t.bool(...) { Err(_) => vec![false; n] }` then +degraded acceptance to the q threshold alone: every transferred identification vanished +from `peptides.tsv` and `proteins.tsv`, indistinguishable from an MBR run that transferred +nothing, while the parquet still carried them. Four sites. The same absent-versus-malformed +mistake was removed from `quant` in package A and from `audit` in package D, and left in +the stage that decides what users read. + +**F3. Script resolution prefers the working directory. Source-confirmed.** + +`python::resolve_script_dir` was reordered in an earlier pass to config-relative, +exe-relative, working-directory-last, with a comment explaining the attack: the shipped +default is the relative `"scripts"`, which both example configurations carry, so unpacking +a dataset archive and running with an example config could execute a worker the archive +contained. That resolver only claims a directory holding `mbr_worker.py` or +`deeplc_worker.py`, so a directory containing any of the other ten workers passed through +unresolved, and `sidecar::resolve_script`, used at all eleven call sites, still tried the +working directory first. + +**F4. A NaN query indexes before the start of the grid. Reproduced.** + +With `x = NaN` both boundary comparisons are false and `partition_point` returns 0, so +`grid_x[j - 1]` underflows: a panic in debug, an out-of-bounds read in release. One library +row with a null `predicted_irt` is enough, because the parquet reader maps a null f32 to +NaN. The anchor loop in `rt_im_train` filters for finite values; the whole-library +application does not, and `align` and `search_seed` share the exposure. + +**F5. The rescorer's three backends do not agree. Reproduced.** + +In-memory TSV standardised with median/IQR while in-memory parquet and the streaming +memmap used mean/std, so the same pool scored differently depending on `rescore.handoff` +and on which side of `MUMDIA_NN_STREAM_GB` the matrix fell, with nothing logged. The +docs/28 comparison cited as identical identifications was measured on a 4.03 GB matrix, so +both arms streamed and the divergence was masked. Separately, the `MUMDIA_NN_FOLD_KEYS` +companion was sliced without a length check: numpy returns a short array rather than +raising, the tail rows land in no fold, are never scored, keep the zero initialiser and +emerge from the final rank-normalisation as a plausible tied mid-rank score, which +satisfies the caller's completeness contract so `rescore.strict` does not catch it. + +**F6. The output-over-input guard is almost unwired. Reproduced.** + +`refuse_output_over_input` had three references in the workspace: its definition and two +call sites. `mumdia compete --features features.parquet --out features.parquet` opens the +table footer-only, streams the surviving rows, and publishes by rename after the read has +finished, so nothing errors and the widest artifact of the run is replaced by the competed +subset at exit 0. Recovering it means re-running extract and features, the tallest stage. +The same shape applies to `rescore`, `features`, `quant`, `extract`, `search_seed`, +`rt-im-train` and `audit`. + +**F7 to F10** are regressions from packages A, D and E and are described in the fix table +below. + +## Package F + +One stacked pull request on `review/e-followup`, with a test per finding. + +| | Change | Test | +|---|---|---| +| F1 | an unbounded window means "screen over the whole gradient", the recall-safe reading of the sentinel; candidates with no window row are treated the same; both are counted and warned about; screening every candidate away is now an error naming the likely causes | `prescan` counts and bail | +| F2 | `transfer_columns` reads present columns in their declared type and errors on a type mismatch; only an absent column falls back | `a_malformed_transfer_column_is_an_error_not_a_silent_loss_of_every_transfer`, `a_table_that_never_saw_mbr_reports_no_transfers_without_complaint` | +| F3 | absolute directories as given, then the executable's directory, then `/scripts`, and the working directory last | `the_shipped_directory_beside_the_binary_wins_over_the_working_directory` and two more | +| F4 | `Loess::predict` returns NaN for a non-finite query; `rt_im_train` treats a non-finite library iRT as "no calibrated RT" and counts it | `a_nan_query_returns_nan_instead_of_indexing_out_of_bounds` | +| F5 | one mean/std standardisation in all three backends, which leaves the shipped parquet default and every published benchmark unchanged; the fold-key companion is length-checked; the backend-size estimate counts feature columns by name | `a_short_fold_key_file_is_refused_rather_than_leaving_rows_unfolded` | +| F6 | the guard is wired into `search_seed`, `rt-im-train`, `extract`, `features`, `compete`, `rescore`, `quant` and `audit`, for every output each writes | `writing_the_output_over_the_input_is_refused` | +| F7 | the extrapolation slope is the secant of the fitted curve over its end decile, clamped non-negative and to at most four times the global least-squares slope | `boundary_extrapolation_stays_monotone_and_bounded_under_noise`, on noisy anchors rather than a noiseless quadratic | +| F8 | the waiter retires the process id the instant `wait` returns, before it reads the output directory | `a_stop_after_the_reap_has_no_pid_to_kill` | +| F9 | bounded take-over attempts with a pause and a deadline; a token written and read back, so only the real holder owns or removes a lock; a future modification time counts as fresh; the partial-file probe matches this destination only; abandoned partials of this destination are swept under the lock | five tests in `raw::tests` | +| F10 | the direct misses and the rows dropped for sharing a pair key are counted separately, and exceeding 2% of the library is an error naming the sidecar | `a_large_unpredicted_fraction_is_a_failure_not_a_warning`, `a_positional_isomer_still_shares_its_pair_key_by_design` | + +### Two changes deliberately not made + +**F10 does not make the pair key position-aware.** `M[Oxidation]PEPTIDEMK` and +`MPEPTIDEM[Oxidation]K` share a base peptide, a charge and a modification multiset, so a +miss on either drops both. Adding position would be worse rather than better: a reverse +decoy carries its modifications at mirrored positions, so a positional key would stop +matching a target to its decoy, and a target could be dropped while its decoy stayed. That +is an FDR defect in exchange for a sensitivity one. The collateral is counted and bounded +instead, and a test pins the trade so a future change has to argue with it. + +**F9 does not add a recipe-keyed conversion cache.** Unique temporary names, an owned lock +and the staleness rules cover the reproduced failure. Two different conversion recipes +aimed at one destination name remain a documented limitation of writing beside the input. + +### F7 merge check + +HYE B01 on doxy, the same arm as the docs/29 #10 pair 2: the DeepLC 4.1.1 re-predicted +precursor table, `rt_im_train.library_irt = library` so no DeepLC draw, `native_tda`, 32 +threads, one binary changed. The comparison arm is the package-D head, whose pointwise +boundary slope this replaces. + +| | package D (pointwise slope) | package F (clamped end-decile secant) | +|---|---|---| +| stripped peptides at `peptide_q_value` 1% | 48,533 | 48,533 | +| PSM-q 1% targets | 53,127 | 53,127 | +| PSM-q 1% decoy fraction | 0.010 | 0.010 | +| protein groups at `pg_q_value` 1% | 6,519 | 6,519 | +| extract accepted rows | 1,961,800 | 1,961,800 | +| `w_rt`, in-sample residual median | 414 s, 78.3 s | 414 s, 78.3 s | +| wall (32 threads) | 25:46 | 24:38 | + +Identical on every count, which is the expected result and the reason the change is safe: +both extrapolations are continuous at the boundary and agree closely wherever the anchors +are dense, and they differ only for queries outside the anchor range, where this benchmark +has almost nothing. The guard exists for the case the counts cannot show, a sparse or noisy +boundary window producing a negative or several-times-global slope, which the unit test +covers directly on noisy anchors. + +## F11: the second tier, not in package F + +Recorded for a later pass, in rough order of value: + +- `matchers/binning.rs`: `LogBins::new` has no bound on the derived bin count while + `validate()` only requires a positive tolerance, so `frag_tol_ppm: 2e-5` asks for 547 GB + and dies on an allocation failure rather than on a message naming the setting. +- `digest.rs`: `collision_safe_decoy` XORs `fnv1a(pep)` and the identical term cancels + inside `make_decoy`, so the scramble seed has no peptide dependence and every peptide of + one length gets the same positional permutation. The shipped `reverse` default only + correlates retried decoys; `strategy = scramble` gets a structured null. +- `ci/smoke.sh`: the NaN-retention-time regression compares two reports through process + substitution with no existence check, so both files missing is a passing `diff`. +- `sbom.cdx.json`: the committed SBOM, which ships in every release archive, references two + components it does not define, so it fails CycloneDX validation and `--check` regenerates + the same document. +- `tests/python/test_predictor_workers.py`: the DeepLC import-order test asserts on two + Windows-only strings and no CI job is both Windows and DeepLC-installed, so it cannot + fail. +- `features/similarity.rs`: `rank_overlap_top3` and `frac_top3_predicted_observed` divide by + a hardcoded 3.0 rather than the number of positive predicted fragments, so they encode how + many fragments were predicted; `mass_uncertainty.rs` computes the same quantity correctly. +- `config.rs`: the `max_feature_matrix_gib` doc comment, copied into the CI-checked + configuration reference, still describes the f64 layout the ceiling no longer uses. +- `run.rs`: the audit is hardcoded at q 0.01 against the report's `quant.q_threshold`. +- `build.rs` cannot distinguish a clean tree from a failed `git status`, so a dirty tree can + be stamped clean. +- `check_doc_refs.py` skips `../..`-prefixed links, `make_fixture_mzml.py` truncates rather + than spans the m/z range when planting peptides, `cargo install tauri-cli --version "^2"` + is a range under a comment calling it pinned, and one workflow pins + `softprops/action-gh-release` at two different revisions. + +## Verified and sound + +Recorded because a review that only lists defects invites the same checks again: all +seventeen documented defaults match the code; no `HashMap` iteration feeds a float +reduction or an output ordering; every configuration struct denies unknown fields and all +five shipped configurations parse; the four DeepLC 4.1.1 enforcement points agree with the +Rust constant and both Python literals match it; neither sidecar leaks a target/decoy label +into the feature matrix; every worker's argument surface matches the arguments its Rust +caller builds; the experiment report headers match the writer in order; and +`split_by_source` runs before per-run quantification. diff --git a/docs/README.md b/docs/README.md index 3ef92df..c577b25 100644 --- a/docs/README.md +++ b/docs/README.md @@ -101,3 +101,4 @@ measurement locally, so the same finding cannot drift into two versions. | [28_feature_selection_analysis.md](28_feature_selection_analysis.md) | Feature-selection analysis for `nn_torch` rescoring: 43 of 387 Extended features are dead by construction, ~120 chosen multivariately reproduce the full set on two HYE runs (real worker confirmed), training time is flat in the feature count, and the parquet handoff is the larger memory lever. | | [29_code_review_2026-09-07.md](29_code_review_2026-09-07.md) | Repository-wide code review at `ab7049e`: 21 findings with reproductions (pooled `source` type, streaming NULLs, entrapment in-sample scores, artifact publication, output ownership, MBR confidence, Met excision, LOESS endpoints, provenance, validation), the agreed four work packages A-D and the deferrals. The decision record the fixes cite by finding number. | | [30_code_review_2026-09-08.md](30_code_review_2026-09-08.md) | Follow-up review at `d7874f2` after work packages A-D: nine remaining findings (fine-tune batch regression, Windows run-name aliases, cancellation cleanup lifetime, concurrent conversions, remaining numeric domains, library rewrite coverage, audit flag consistency, overlapping reservations, debug-binary stack), the status of the original 21, and the implementation status of package E. | +| [31_code_review_2026-09-08_full.md](31_code_review_2026-09-08_full.md) | Whole-repository review at `21b8ac7`, after work packages A-E: 21 findings (prescan's no-calibration sentinel, the malformed transfer flag, script-resolution order, a NaN LOESS query, three rescorer standardisations, the unwired output-over-input guard, and five regressions from A/D/E), work package F, and the second tier left for later. | diff --git a/rust/mumdia/crates/mumdia/src/calibrate.rs b/rust/mumdia/crates/mumdia/src/calibrate.rs index 4111157..b9902ea 100644 --- a/rust/mumdia/crates/mumdia/src/calibrate.rs +++ b/rust/mumdia/crates/mumdia/src/calibrate.rs @@ -45,8 +45,19 @@ pub struct Loess { /// extrapolation model. slope: f64, intercept: f64, - /// Slope of the local line at the first and at the last grid point; extrapolation - /// continues these from the boundary values, so `predict` is continuous there. + /// Extrapolation slope below the first and above the last grid point; the boundary + /// grid VALUE plus this slope keeps `predict` continuous at both ends. + /// + /// The secant of the fitted curve over its end decile, not the pointwise local slope + /// at the boundary. The pointwise slope comes from the one place where the tricubic + /// window is one-sided and the anchors are sparsest, so it is the noisiest number the + /// fit produces, unbounded, and free to be negative: on 55 anchors with 120 s of + /// scatter it ranged 6.9 to 28.4 against a global slope near 24, and a peptidoform + /// one iRT decade outside the anchors landed 1077 s from where the previous global + /// line put it, against an `w_rt` near 180 s (docs/31 F7). A negative slope inverts + /// the iRT-to-RT map, which the global line it replaced could not do. The secant + /// averages the same curve over many grid points, and the clamp below keeps the + /// result monotone and within a factor of the global least-squares slope. lo_slope: f64, hi_slope: f64, } @@ -86,19 +97,13 @@ impl Loess { let gn = grid_n.max(2); let mut grid_x = Vec::with_capacity(gn); let mut grid_y = Vec::with_capacity(gn); - let (mut lo_slope, mut hi_slope) = (slope, slope); for g in 0..gn { let x = lo + (hi - lo) * g as f64 / (gn - 1) as f64; - let (y, b) = local_linear(&sx, &sy, x, k, slope, intercept); + let (y, _) = local_linear(&sx, &sy, x, k, slope, intercept); grid_x.push(x); grid_y.push(y); - if g == 0 { - lo_slope = b; - } - if g == gn - 1 { - hi_slope = b; - } } + let (lo_slope, hi_slope) = boundary_slopes(&grid_x, &grid_y, slope); Loess { grid_x, grid_y, @@ -114,6 +119,15 @@ impl Loess { /// slope), so the prediction is continuous across the boundary; see the type's /// documentation for what the previous global-line switch did. pub fn predict(&self, x: f64) -> f64 { + // NaN in, NaN out. With x = NaN both boundary comparisons are false and + // `partition_point` returns 0, so the interpolation below indexed `grid_x[-1]`: + // a panic in debug and an out-of-bounds read in release, reachable from one + // library row whose `predicted_irt` is null, because the parquet reader maps a + // null f32 to NaN (docs/31 F4). Callers treat a NaN prediction as "no + // calibration for this row", which is the documented recall-safe sentinel. + if !x.is_finite() { + return f64::NAN; + } let g = &self.grid_x; if g.len() < 2 { return self.slope * x + self.intercept; @@ -135,6 +149,45 @@ impl Loess { } } +/// Extrapolation slopes for the two ends of a fitted grid. +/// +/// Each is the secant of the fitted curve over its end decile, clamped to be +/// non-negative (an inverted iRT-to-RT map is never right) and, when the global +/// least-squares slope is positive, to at most four times it. Both guards exist because +/// this number multiplies a distance that has no bound: it is applied to queries outside +/// the anchor range entirely. See the field documentation on [`Loess`]. +fn boundary_slopes(grid_x: &[f64], grid_y: &[f64], global_slope: f64) -> (f64, f64) { + let n = grid_x.len(); + if n < 2 { + return (global_slope, global_slope); + } + let step = (n / 10).max(1).min(n - 1); + let secant = |a: usize, b: usize| { + let dx = grid_x[b] - grid_x[a]; + if dx.abs() < 1e-12 { + global_slope + } else { + (grid_y[b] - grid_y[a]) / dx + } + }; + let cap = if global_slope.is_finite() && global_slope > 0.0 { + 4.0 * global_slope + } else { + f64::INFINITY + }; + let clamp = |v: f64| { + if !v.is_finite() { + return if global_slope.is_finite() { + global_slope.max(0.0) + } else { + 0.0 + }; + } + v.clamp(0.0, cap) + }; + (clamp(secant(0, step)), clamp(secant(n - 1 - step, n - 1))) +} + /// Weighted local linear regression at `x0` over the `k` nearest points /// (tricubic weights): the fitted value at `x0` and the slope of the local line. /// Falls back to the global line (value and slope) if degenerate. @@ -257,6 +310,61 @@ mod tests { ); } + #[test] + fn a_nan_query_returns_nan_instead_of_indexing_out_of_bounds() { + // docs/31 F4: one library row with a null predicted_irt reaches this as NaN. + let xs: Vec = (0..40).map(|i| i as f64).collect(); + let ys: Vec = xs.iter().map(|x| 2.0 * x + 1.0).collect(); + let lo = Loess::fit(&xs, &ys, 0.3, 20); + assert!(lo.predict(f64::NAN).is_nan()); + assert!(lo.predict(f64::INFINITY).is_nan()); + assert!(lo.predict(f64::NEG_INFINITY).is_nan()); + // A degenerate fit takes the linear path and must not panic either. + let tiny = Loess::fit(&[1.0, 2.0], &[1.0, 2.0], 0.5, 2); + assert!(tiny.predict(f64::NAN).is_nan()); + } + + #[test] + fn boundary_extrapolation_stays_monotone_and_bounded_under_noise() { + // docs/31 F7: the pointwise boundary slope came from the sparsest, most one-sided + // local window in the fit. On noisy anchors it was free to be negative (inverting + // the iRT-to-RT map) or several times the global slope, and it multiplies a + // distance that is unbounded because it applies outside the anchor range. + // Deterministic pseudo-noise so the assertion is reproducible. + let mut seed = 0x243f_6a88_85a3_08d3u64; + let mut noise = || { + seed = seed + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((seed >> 33) as f64 / (1u64 << 31) as f64 - 1.0) * 120.0 + }; + let xs: Vec = (0..55).map(|i| i as f64 * 2.0).collect(); + let ys: Vec = xs.iter().map(|x| 500.0 + 24.0 * x + noise()).collect(); + let lo = Loess::fit(&xs, &ys, 0.3, 200); + let (slope, _) = linear_fit(&xs, &ys); + for s in [lo.lo_slope, lo.hi_slope] { + assert!(s >= 0.0, "an inverted extrapolation slope: {s}"); + assert!( + s <= 4.0 * slope, + "an unbounded extrapolation slope: {s} vs {slope}" + ); + } + // Monotone outside the range, and continuous at both ends. + let x0 = xs[0]; + let xn = xs[xs.len() - 1]; + assert!(lo.predict(x0 - 30.0) <= lo.predict(x0)); + assert!(lo.predict(xn + 30.0) >= lo.predict(xn)); + assert!((lo.predict(x0) - lo.predict(x0 + 1e-6)).abs() < 1.0); + assert!((lo.predict(xn) - lo.predict(xn - 1e-6)).abs() < 1.0); + // And the error one decade outside the anchors stays within a search window. + let far = lo.predict(x0 - 30.0); + let linear = slope * (x0 - 30.0) + lo.intercept; + assert!( + (far - linear).abs() < 600.0, + "extrapolation {far} is implausibly far from the global line {linear}" + ); + } + #[test] fn percentile_basic() { let v: Vec = (0..=100).map(|i| i as f64).collect(); diff --git a/rust/mumdia/crates/mumdia/src/raw.rs b/rust/mumdia/crates/mumdia/src/raw.rs index 9a186dc..5dda7cb 100644 --- a/rust/mumdia/crates/mumdia/src/raw.rs +++ b/rust/mumdia/crates/mumdia/src/raw.rs @@ -539,6 +539,10 @@ fn unique_tag() -> String { /// into the same destination. Dropping the guard releases the lock. struct ConvertLock { path: PathBuf, + /// Written into the lock file and verified after creation, so a lock this process + /// took cannot be one another process took a moment earlier, and `Drop` cannot + /// delete a lock that is no longer ours. + token: String, } impl ConvertLock { @@ -553,9 +557,22 @@ impl ConvertLock { /// Take the lock for `out`, waiting for a holder to finish first. `Ok(None)` means the /// holder finished and left a usable conversion at `out`, which the caller reuses. + /// + /// Three things this has to survive, all of which the first version did not + /// (docs/31 F9). A stale lock that cannot be deleted is bounded rather than retried + /// forever with no pause. Two waiters that both judge one lock stale cannot both come + /// away holding it, because each writes a token and reads it back, and only the + /// process whose token survives owns the lock. And waiting itself has a deadline, so + /// a lock nothing will ever release fails the run with a message instead of hanging. fn acquire(out: &Path, src: &Path, reuse: bool) -> Result> { + const RETRY: std::time::Duration = std::time::Duration::from_secs(2); + const MAX_TAKEOVERS: u32 = 8; + const MAX_WAIT: std::time::Duration = std::time::Duration::from_secs(12 * 3600); let path = Self::path_for(out); + let token = format!("{}:{}", std::process::id(), unique_tag()); + let waiting_since = std::time::Instant::now(); let mut announced = false; + let mut takeovers = 0u32; loop { match std::fs::OpenOptions::new() .write(true) @@ -564,17 +581,47 @@ impl ConvertLock { { Ok(mut f) => { use std::io::Write; - let _ = writeln!(f, "{}", std::process::id()); - return Ok(Some(ConvertLock { path })); + f.write_all(token.as_bytes()) + .and_then(|()| f.sync_all()) + .with_context(|| { + format!("writing the conversion lock {}", path.display()) + })?; + drop(f); + // Read it back. If another waiter removed this file and created its + // own between the two calls, the token no longer matches and the lock + // is theirs; wait and try again rather than converting alongside them. + let held = std::fs::read_to_string(&path) + .map(|t| t.trim() == token) + .unwrap_or(false); + if held { + return Ok(Some(ConvertLock { path, token })); + } + std::thread::sleep(RETRY); + continue; } Err(e) if e.kind() == std::io::ErrorKind::AlreadyExists => { - if lock_is_stale(&path) { + if lock_is_stale(&path, out) { + takeovers += 1; + if takeovers > MAX_TAKEOVERS { + bail!( + "the conversion lock {} looks abandoned but cannot be \ + removed after {MAX_TAKEOVERS} attempts. Delete it by hand \ + once nothing is converting {}, then rerun.", + path.display(), + out.display() + ); + } warn!( lock = %path.display(), + attempt = takeovers, "convert: removing a stale conversion lock; its holder stopped \ writing" ); + // A failed removal used to `continue` straight back into + // `create_new` with no pause, pinning a core and writing one + // warning per iteration for as long as the run lasted. let _ = std::fs::remove_file(&path); + std::thread::sleep(RETRY); continue; } if !announced { @@ -585,7 +632,16 @@ impl ConvertLock { ); announced = true; } - std::thread::sleep(std::time::Duration::from_secs(2)); + if waiting_since.elapsed() > MAX_WAIT { + bail!( + "waited {} hours for the conversion lock {} and it is still \ + held and still fresh. Check for another MuMDIA converting {}.", + MAX_WAIT.as_secs() / 3600, + path.display(), + out.display() + ); + } + std::thread::sleep(RETRY); if !path.exists() && reuse && out.is_file() && is_newer_than(out, src) { return Ok(None); } @@ -602,30 +658,64 @@ impl ConvertLock { impl Drop for ConvertLock { fn drop(&mut self) { - let _ = std::fs::remove_file(&self.path); + // Only if it is still ours: a lock taken over after this one was judged stale + // belongs to the taker, and removing it would leave two converters writing one + // destination, which is the situation the lock exists to prevent. + let ours = std::fs::read_to_string(&self.path) + .map(|t| t.trim() == self.token) + .unwrap_or(false); + if ours { + let _ = std::fs::remove_file(&self.path); + } } } -/// A lock whose holder has stopped: neither the lock nor any partial conversion file -/// beside it has been written for `STALE_AFTER`. Converters write their output -/// continuously, so a live conversion keeps a `.partial-` file fresh; a crashed or killed -/// holder leaves both untouched. -fn lock_is_stale(lock: &Path) -> bool { - const STALE_AFTER: std::time::Duration = std::time::Duration::from_secs(15 * 60); - let fresh = |p: &Path| { - std::fs::metadata(p) - .and_then(|m| m.modified()) - .ok() - .and_then(|t| t.elapsed().ok()) - .is_some_and(|age| age < STALE_AFTER) - }; - if fresh(lock) { +/// How long a lock and its partial output may go untouched before the holder counts as +/// gone. +const STALE_AFTER: std::time::Duration = std::time::Duration::from_secs(15 * 60); + +/// Was `p` written within [`STALE_AFTER`]? +/// +/// A modification time in the future counts as fresh. It is an unreadable clock, not +/// evidence that nobody is writing: `SystemTime::elapsed` returns `Err` for a future +/// timestamp, and treating that as "not fresh" made ordinary clock skew on a network +/// share declare a live lock stale one second after it was taken (docs/31 F9). +fn written_recently(p: &Path) -> bool { + match std::fs::metadata(p).and_then(|m| m.modified()) { + Ok(t) => match t.elapsed() { + Ok(age) => age < STALE_AFTER, + Err(_) => true, + }, + Err(_) => false, + } +} + +/// The `.partial-` prefix of the temporary files belonging to one destination. +fn partial_prefix(out: &Path) -> String { + let name = out.file_name().unwrap_or_default().to_string_lossy(); + let stem = name.strip_suffix(".mzML").unwrap_or(&name); + format!("{stem}.partial-") +} + +/// A lock whose holder has stopped: neither the lock nor a partial conversion file OF +/// THIS DESTINATION has been written for [`STALE_AFTER`]. Converters write continuously, +/// so a live conversion keeps its partial file fresh; a crashed or killed holder leaves +/// both untouched. +/// +/// Matching this destination's prefix rather than any `.partial-` name is what makes the +/// probe mean anything with several conversions in one directory: it used to accept any +/// partial file, so one live conversion kept a dead peer's lock fresh indefinitely, and +/// with `experiment.parallel_runs > 1` that is the normal case (docs/31 F9). +fn lock_is_stale(lock: &Path, out: &Path) -> bool { + if written_recently(lock) { return false; } let dir = lock.parent().unwrap_or(Path::new(".")); + let prefix = partial_prefix(out); if let Ok(rd) = std::fs::read_dir(dir) { for e in rd.flatten() { - if e.file_name().to_string_lossy().contains(".partial-") && fresh(&e.path()) { + let name = e.file_name().to_string_lossy().into_owned(); + if name.starts_with(&prefix) && written_recently(&e.path()) { return false; } } @@ -633,6 +723,41 @@ fn lock_is_stale(lock: &Path) -> bool { true } +/// Remove this destination's abandoned partial conversions. +/// +/// Called with the lock held, so anything matching that is older than [`STALE_AFTER`] +/// belongs to a conversion that died. Giving every attempt a unique temporary name fixed +/// one bug and created this one: the pre-conversion `remove_file` then named a file that +/// cannot exist yet, so a killed conversion left its multi-gigabyte partial mzML beside +/// the input for ever and nothing anywhere deleted it (docs/31 F9). +fn sweep_stale_partials(out_dir: &Path, out: &Path) { + let prefix = partial_prefix(out); + let Ok(rd) = std::fs::read_dir(out_dir) else { + return; + }; + for e in rd.flatten() { + let name = e.file_name().to_string_lossy().into_owned(); + if !name.starts_with(&prefix) || !name.ends_with(".mzML") { + continue; + } + let path = e.path(); + if written_recently(&path) { + continue; + } + match std::fs::remove_file(&path) { + Ok(()) => warn!( + partial = %path.display(), + "convert: removed an abandoned partial conversion" + ), + Err(e) => warn!( + partial = %path.display(), + error = %e, + "convert: could not remove an abandoned partial conversion" + ), + } + } +} + pub fn ensure_mzml( input: &str, cfg: &ConvertConfig, @@ -762,9 +887,11 @@ msconvert was not usable either: {e}" // `reuse_converted` accepted it on every later run: a partial acquisition searched // silently, for ever, presenting as unexplained low identification counts with // nothing in the interface to reveal it. Neither failure path removed it. + // Under the lock, so every abandoned partial of THIS destination is dead by + // definition. The unique name below cannot collide with a live one. + sweep_stale_partials(&out_dir, &out); let tmp_name = partial_name(&out_name, &unique_tag()); let tmp = out_dir.join(&tmp_name); - let _ = std::fs::remove_file(&tmp); let args: Vec = if is_thermo { // `-f 2` is indexed mzML, which is what msconvert produces by default and so @@ -989,27 +1116,121 @@ mod tests { let _ = std::fs::remove_dir_all(&d); } + /// Backdate a file's modification time by `secs`. + fn age(path: &Path, secs: u64) { + let t = std::time::SystemTime::now() - std::time::Duration::from_secs(secs); + std::fs::OpenOptions::new() + .write(true) + .open(path) + .unwrap() + .set_modified(t) + .unwrap(); + } + #[test] fn a_stale_lock_is_broken_and_a_fresh_one_is_honoured() { let d = tmp("stale"); + let out = d.join("run.mzML"); let lock = d.join("run.mzML.converting"); std::fs::write(&lock, b"1").unwrap(); - assert!(!lock_is_stale(&lock), "a lock written a moment ago is live"); - // Age the lock past the staleness window; nothing partial is being written. - let old = std::time::SystemTime::now() - std::time::Duration::from_secs(20 * 60); + assert!( + !lock_is_stale(&lock, &out), + "a lock written a moment ago is live" + ); + age(&lock, 20 * 60); + assert!( + lock_is_stale(&lock, &out), + "an old lock with no live partial file is stale" + ); + // A partial file of THIS destination, still being written, keeps an old lock alive. + std::fs::write(d.join("run.partial-99-0.mzML"), b"...").unwrap(); + assert!(!lock_is_stale(&lock, &out)); + let _ = std::fs::remove_dir_all(&d); + } + + #[test] + fn another_destinations_partial_does_not_keep_this_lock_alive() { + // docs/31 F9: the probe accepted any `.partial-` name, so one live conversion in a + // shared directory kept a dead peer's lock fresh for ever. With + // `experiment.parallel_runs > 1` that is the ordinary case. + let d = tmp("peer"); + let out = d.join("dead.mzML"); + let lock = d.join("dead.mzML.converting"); + std::fs::write(&lock, b"1").unwrap(); + age(&lock, 20 * 60); + std::fs::write(d.join("alive.partial-1-0.mzML"), b"...").unwrap(); + assert!( + lock_is_stale(&lock, &out), + "a peer's live partial file must not vouch for this destination" + ); + std::fs::write(d.join("dead.partial-1-0.mzML"), b"...").unwrap(); + assert!(!lock_is_stale(&lock, &out), "its own partial file does"); + let _ = std::fs::remove_dir_all(&d); + } + + #[test] + fn a_future_modification_time_reads_as_fresh_not_stale() { + // Clock skew on a network share, not evidence that nobody is writing. + let d = tmp("skew"); + let f = d.join("run.mzML.converting"); + std::fs::write(&f, b"1").unwrap(); + let future = std::time::SystemTime::now() + std::time::Duration::from_secs(3600); std::fs::OpenOptions::new() .write(true) - .open(&lock) + .open(&f) .unwrap() - .set_modified(old) + .set_modified(future) .unwrap(); + assert!(written_recently(&f)); + assert!(!lock_is_stale(&f, &d.join("run.mzML"))); + let _ = std::fs::remove_dir_all(&d); + } + + #[test] + fn abandoned_partials_are_swept_and_live_ones_are_not() { + let d = tmp("sweep"); + let out = d.join("run.mzML"); + let dead = d.join("run.partial-1-0.mzML"); + let live = d.join("run.partial-2-0.mzML"); + let peer = d.join("other.partial-3-0.mzML"); + let real = d.join("run.mzML"); + for f in [&dead, &live, &peer, &real] { + std::fs::write(f, b"x").unwrap(); + } + age(&dead, 20 * 60); + age(&peer, 20 * 60); + sweep_stale_partials(&d, &out); assert!( - lock_is_stale(&lock), - "an old lock with no live partial file is stale" + !dead.is_file(), + "an abandoned partial of this destination goes" + ); + assert!(live.is_file(), "a partial still being written stays"); + assert!( + peer.is_file(), + "another destination's partial is not ours to remove" + ); + assert!(real.is_file(), "the converted mzML itself is untouched"); + let _ = std::fs::remove_dir_all(&d); + } + + #[test] + fn a_lock_is_only_removed_by_the_holder_that_still_owns_it() { + // docs/31 F9: two waiters that both judged a lock stale each removed it and each + // created their own, so the second unlinked the first's and both believed they + // held it; the first's Drop then removed the second's. + let d = tmp("token"); + let out = d.join("run.mzML"); + std::fs::write(d.join("run.raw"), b"raw").unwrap(); + let first = ConvertLock::acquire(&out, &d.join("run.raw"), true) + .unwrap() + .expect("first holder"); + // Simulate a take-over: another process replaces the lock with its own token. + std::fs::write(ConvertLock::path_for(&out), b"9999:other").unwrap(); + drop(first); + assert!( + ConvertLock::path_for(&out).is_file(), + "dropping a lock we no longer own must not remove the new holder's" ); - // A partial file still being written keeps even an old lock alive. - std::fs::write(d.join("run.partial-99-0.mzML"), b"...").unwrap(); - assert!(!lock_is_stale(&lock)); let _ = std::fs::remove_dir_all(&d); } diff --git a/rust/mumdia/crates/mumdia/src/sidecar.rs b/rust/mumdia/crates/mumdia/src/sidecar.rs index 168abaf..da5420f 100644 --- a/rust/mumdia/crates/mumdia/src/sidecar.rs +++ b/rust/mumdia/crates/mumdia/src/sidecar.rs @@ -17,27 +17,43 @@ pub type FragmentIntensityMap = HashMap>; /// Resolve a sidecar worker script path so a deployed binary finds its workers /// regardless of the working directory: try the configured dir relative to the -/// CWD, then relative to the binary's own directory, then `/scripts`. -/// Falls back to the CWD-relative path (so the eventual error names it) if none -/// of those exist. +/// the binary's own directory, then `/scripts`, and the current working +/// directory LAST. Falls back to the directory-relative path (so the eventual error +/// names it) if none of those exist. +/// +/// The working directory used to be tried first. `python::resolve_script_dir` was +/// reordered away from exactly that and documents why at length: the shipped default is +/// the relative `"scripts"`, which both sidecar example configs carry, so unpacking a +/// dataset archive, `cd`-ing into it and running with an example config executed any +/// worker the archive happened to contain. That resolver only claims a directory holding +/// `mbr_worker.py` or `deeplc_worker.py`, so a directory with any of the other ten +/// workers reached this function still relative, and this function ran it (docs/31 F3). +/// An absolute directory is taken as given: naming one is how a user is unambiguous. pub fn resolve_script(dir: &str, worker: &str) -> String { - let cwd_rel = format!("{dir}/{worker}"); - if std::path::Path::new(&cwd_rel).exists() { - return cwd_rel; + let exe_dir = std::env::current_exe() + .ok() + .and_then(|e| e.parent().map(|p| p.to_path_buf())); + resolve_script_in(dir, worker, exe_dir.as_deref()) +} + +/// [`resolve_script`] with the executable's directory supplied, so the ordering can be +/// tested without a second binary. +fn resolve_script_in(dir: &str, worker: &str, exe_dir: Option<&std::path::Path>) -> String { + let dir_rel = format!("{dir}/{worker}"); + if std::path::Path::new(dir).is_absolute() { + return dir_rel; } - if let Ok(exe) = std::env::current_exe() { - if let Some(base) = exe.parent() { - for cand in [ - base.join(dir).join(worker), - base.join("scripts").join(worker), - ] { - if cand.exists() { - return cand.to_string_lossy().into_owned(); - } + if let Some(base) = exe_dir { + for cand in [ + base.join(dir).join(worker), + base.join("scripts").join(worker), + ] { + if cand.exists() { + return cand.to_string_lossy().into_owned(); } } } - cwd_rel + dir_rel } /// MS2PIP: predict singly-charged b/y intensities per (peptidoform, charge). @@ -440,3 +456,67 @@ fn run_worker(python: &str, script: &str, args: &[&str], utf8: bool) -> Result<( } Ok(()) } + +#[cfg(test)] +mod resolve_tests { + use super::resolve_script_in; + use std::path::Path; + + fn scratch(name: &str) -> std::path::PathBuf { + let d = std::env::temp_dir().join(format!("mumdia_resolve_{}_{name}", std::process::id())); + let _ = std::fs::remove_dir_all(&d); + std::fs::create_dir_all(&d).unwrap(); + d + } + + #[test] + fn the_shipped_directory_beside_the_binary_wins_over_the_working_directory() { + // docs/31 F3. Both locations hold a worker of the same name; the one that ships + // with the binary must win, because the working directory can be an untrusted + // dataset the user merely unpacked and `cd`-ed into. + let exe = scratch("exe"); + std::fs::create_dir_all(exe.join("scripts")).unwrap(); + std::fs::write(exe.join("scripts").join("ms2pip_worker.py"), b"# shipped").unwrap(); + let got = resolve_script_in("scripts", "ms2pip_worker.py", Some(&exe)); + assert_eq!( + got, + exe.join("scripts") + .join("ms2pip_worker.py") + .to_string_lossy() + ); + assert_ne!(got, "scripts/ms2pip_worker.py"); + let _ = std::fs::remove_dir_all(&exe); + } + + #[test] + fn an_absolute_directory_is_taken_as_given() { + let abs = if cfg!(windows) { + "C:/opt/mumdia/scripts" + } else { + "/opt/mumdia/scripts" + }; + let exe = scratch("abs"); + std::fs::create_dir_all(exe.join("scripts")).unwrap(); + std::fs::write(exe.join("scripts").join("mbr_worker.py"), b"# shipped").unwrap(); + assert_eq!( + resolve_script_in(abs, "mbr_worker.py", Some(&exe)), + format!("{abs}/mbr_worker.py"), + "naming a directory outright is how a user is unambiguous" + ); + let _ = std::fs::remove_dir_all(&exe); + } + + #[test] + fn nothing_beside_the_binary_falls_back_to_the_relative_path_for_the_error() { + let exe = scratch("empty"); + assert_eq!( + resolve_script_in("scripts", "deeplc_worker.py", Some(&exe)), + "scripts/deeplc_worker.py" + ); + assert_eq!( + resolve_script_in("scripts", "deeplc_worker.py", None::<&Path>), + "scripts/deeplc_worker.py" + ); + let _ = std::fs::remove_dir_all(&exe); + } +} diff --git a/rust/mumdia/crates/mumdia/src/stages/audit.rs b/rust/mumdia/crates/mumdia/src/stages/audit.rs index 8b90198..c8aec13 100644 --- a/rust/mumdia/crates/mumdia/src/stages/audit.rs +++ b/rust/mumdia/crates/mumdia/src/stages/audit.rs @@ -64,6 +64,18 @@ fn load_extract_reasons(psms_path: &str) -> HashMap { pub fn run(p: AuditParams) -> Result { let t0 = Instant::now(); + // `--out` must not be one of this stage's own inputs: every input is read + // before the output is published, so writing over one replaces it and exits 0 + // (docs/31 F6). The shared guard existed and was wired into two stages. + mumdia_io::refuse_output_over_input( + p.out, + &[ + ("--lib-precursors", p.library_precursors), + ("--psms", p.psms), + ("--competed", p.competed), + ("--psms-scored", p.scored), + ], + )?; // Search space = all library precursors. let lib = TableFile::open(p.library_precursors) diff --git a/rust/mumdia/crates/mumdia/src/stages/compete.rs b/rust/mumdia/crates/mumdia/src/stages/compete.rs index 4c22d5f..3441ccf 100644 --- a/rust/mumdia/crates/mumdia/src/stages/compete.rs +++ b/rust/mumdia/crates/mumdia/src/stages/compete.rs @@ -32,6 +32,10 @@ pub struct CompeteParams<'a> { pub fn run(p: CompeteParams) -> Result { let t0 = Instant::now(); + // `--out` must not be one of this stage's own inputs: every input is read + // before the output is published, so writing over one replaces it and exits 0 + // (docs/31 F6). The shared guard existed and was wired into two stages. + mumdia_io::refuse_output_over_input(p.out, &[("--features", p.features)])?; // Footer-only open. The key columns below stream one at a time and the feature columns // (hundreds of them) are never materialised: the previous path read the whole features // table into Arrow and then copied every column into an owned Vec, so compete held two @@ -539,6 +543,44 @@ fn resolve_competition( mod tests { use super::*; + #[test] + fn writing_the_output_over_the_input_is_refused() { + // docs/31 F6: the features table is opened footer-only and streamed, so the read + // completes before `AtomicPath::publish` renames over it. Nothing errored: the + // widest artifact of the run, hundreds of columns and gigabytes on a real library, + // was replaced by the competed subset at exit 0, recoverable only by re-running + // extract and features. + let dir = std::env::temp_dir().join(format!("mumdia_compete_guard_{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let features = dir.join("features.parquet"); + mumdia_io::table::write_table( + features.to_str().unwrap(), + vec![mumdia_io::table::Col::U32( + "candidate_id".into(), + vec![1, 2], + )], + ) + .unwrap(); + let cfg = mumdia_core::config::CompeteConfig::default(); + let e = run(CompeteParams { + features: features.to_str().unwrap(), + out: features.to_str().unwrap(), + cfg: &cfg, + config_hash: "h", + }) + .unwrap_err() + .to_string(); + assert!(e.contains("its own input"), "{e}"); + assert!( + mumdia_io::table::TableFile::open(features.to_str().unwrap()) + .unwrap() + .nrows + == 2, + "the input must still be there" + ); + let _ = std::fs::remove_dir_all(&dir); + } + fn one_group(members: Vec) -> HashMap<(u32, u8, i64, i32), Vec> { let mut g = HashMap::new(); g.insert((0u32, 0u8, 0i64, 0i32), members); diff --git a/rust/mumdia/crates/mumdia/src/stages/extract.rs b/rust/mumdia/crates/mumdia/src/stages/extract.rs index f57bdef..abef7bd 100644 --- a/rust/mumdia/crates/mumdia/src/stages/extract.rs +++ b/rust/mumdia/crates/mumdia/src/stages/extract.rs @@ -1430,6 +1430,14 @@ fn extract_twopass_windows( pub fn run(p: ExtractParams) -> Result<(u64, u64)> { let t0 = Instant::now(); + // Neither output may be one of the inputs (docs/31 F6). + let inputs = [ + ("--ms2", p.ms2), + ("--lib-precursors", p.library_precursors), + ("--lib-fragments", p.library_fragments), + ]; + mumdia_io::refuse_output_over_input(p.out_psms, &inputs)?; + mumdia_io::refuse_output_over_input(p.out_chrom, &inputs)?; // Skip the bucketed page_search index when the fragindex backend is selected (the // default): it is never read on that path and costs a full sort plus several full // copies of every library fragment. diff --git a/rust/mumdia/crates/mumdia/src/stages/features.rs b/rust/mumdia/crates/mumdia/src/stages/features.rs index 40f0d2e..7a678ed 100644 --- a/rust/mumdia/crates/mumdia/src/stages/features.rs +++ b/rust/mumdia/crates/mumdia/src/stages/features.rs @@ -1217,6 +1217,15 @@ fn confident_global_bounds( } pub fn run(p: FeaturesParams) -> Result { + // Neither output may be one of the inputs (docs/31 F6). + let mut inputs = vec![("--psms", p.psms), ("--chromatograms", p.chromatograms)]; + if let Some(seed) = p.seed { + inputs.push(("--seed-psms", seed)); + } + mumdia_io::refuse_output_over_input(p.out, &inputs)?; + if !p.out_pin.is_empty() { + mumdia_io::refuse_output_over_input(p.out_pin, &inputs)?; + } run_with_chunk_rows(p, CHUNK_CHROM_ROWS) } diff --git a/rust/mumdia/crates/mumdia/src/stages/predict_frag.rs b/rust/mumdia/crates/mumdia/src/stages/predict_frag.rs index 8e885f1..1ec1253 100644 --- a/rust/mumdia/crates/mumdia/src/stages/predict_frag.rs +++ b/rust/mumdia/crates/mumdia/src/stages/predict_frag.rs @@ -139,10 +139,28 @@ pub fn run(p: PredictFragParams) -> Result<(u64, u64)> { // the FDR rests on is kept. Substituting a value was the previous behaviour, iRT 0.0 // for DeepLC and the native heuristic for MS2PIP, and it produced a plausible, finite // library whose rows came from an unrecorded mixture of predictors (docs/29 #17). + let n_before_drop = raws.len() as u64; let (n_dropped_rows, n_dropped_pairs) = drop_unpredicted(&mut raws, &rt_missing, &frag_missing); + // The rows the predictors actually missed, and the rows that went with them for + // sharing a pair key. Reported separately because the second number is the cost of + // the position-free key (see `rows_to_drop`) and the first is the sidecar's own miss + // rate; one warning conflating them hid a worker returning nothing for a large batch. + let n_direct = { + let mut d: Vec = rt_missing + .iter() + .chain(frag_missing.iter()) + .copied() + .collect(); + d.sort_unstable(); + d.dedup(); + d.len() as u64 + }; + let n_collateral = n_dropped_rows.saturating_sub(n_direct); if n_dropped_rows > 0 { tracing::warn!( candidates_dropped = n_dropped_rows, + unpredicted = n_direct, + dropped_with_their_pair = n_collateral, pairs_dropped = n_dropped_pairs, without_irt = rt_missing.len(), without_intensities = frag_missing.len(), @@ -150,6 +168,17 @@ pub fn run(p: PredictFragParams) -> Result<(u64, u64)> { instead of receiving a substitute value; the counts are in the library report" ); } + if let Some(frac) = dropped_fraction_exceeded(n_dropped_rows, n_before_drop) { + bail!( + "predict-frag: {n_dropped_rows} of {n_before_drop} candidates ({:.2}%) have no \ + prediction ({n_direct} the predictor missed, {n_collateral} dropped with their \ + pair) and would silently leave the library. That is a failing sidecar rather \ + than a few unsupported peptidoforms; the ceiling is {:.0}%. Check the worker's \ + output above, then rerun.", + frac * 100.0, + MAX_DROPPED_FRACTION * 100.0 + ); + } if raws.is_empty() && (n_dropped_rows > 0 || n_parse_err > 0) { bail!( "predict-frag: no candidate received a prediction ({n_dropped_rows} dropped for \ @@ -295,6 +324,10 @@ pub fn run(p: PredictFragParams) -> Result<(u64, u64)> { stats.insert("candidates".to_string(), json!(n_prec)); stats.insert("fragments".to_string(), json!(n_frag)); stats.insert("parse_errors".to_string(), json!(n_parse_err)); + stats.insert( + "candidates_dropped_with_their_pair".to_string(), + json!(n_collateral), + ); stats.insert( "candidates_dropped_unpredicted".to_string(), json!(n_dropped_rows), @@ -496,8 +529,34 @@ fn pair_key(base_peptide_id: u32, charge: i32, peptidoform: &str) -> (u32, i32, (base_peptide_id, charge, mod_signature(peptidoform)) } +/// The largest fraction of the library that may be dropped for want of a prediction +/// before the build is treated as a failed sidecar rather than a few unsupported +/// peptidoforms. +/// +/// The only previous guard fired when the library was emptied entirely, so a worker that +/// returned nothing for 5% of a 9.8M-peptidoform batch removed those rows, plus every +/// positional isomer sharing their signature, behind one warning (docs/31 F10). +const MAX_DROPPED_FRACTION: f64 = 0.02; + +/// The dropped fraction when it exceeds [`MAX_DROPPED_FRACTION`], else `None`. +fn dropped_fraction_exceeded(dropped: u64, total: u64) -> Option { + if total == 0 || dropped == 0 { + return None; + } + let frac = dropped as f64 / total as f64; + (frac > MAX_DROPPED_FRACTION).then_some(frac) +} + /// Which rows to drop so that every row sharing a pair key with an unpredicted row goes /// with it. `keys[i]` is the pair key of row `i`; `missing` lists unpredicted rows. +/// +/// The key is deliberately position-free, and that over-deletes: `M[Oxidation]PEPTIDEMK` +/// and `MPEPTIDEM[Oxidation]K` share a base peptide, a charge and a modification +/// multiset, so a miss on either removes both. Making the key position-aware would be +/// worse, not better: a reverse decoy carries its modifications at mirrored positions, so +/// a positional key would stop matching a target to its decoy and a target could be +/// dropped while its decoy stayed, which is an FDR defect rather than a sensitivity one. +/// The collateral is counted separately and bounded by [`MAX_DROPPED_FRACTION`] instead. fn rows_to_drop(keys: &[(u32, i32, String)], missing: &[usize]) -> Vec { let doomed: std::collections::HashSet<&(u32, i32, String)> = missing.iter().map(|&i| &keys[i]).collect(); @@ -653,7 +712,10 @@ fn fragment_cardinality(cid: &[u32], mz: &[f64]) -> Vec { #[cfg(test)] mod cardinality_tests { - use super::{fragment_cardinality, mod_signature, ms2pip_values, pair_key, rows_to_drop}; + use super::{ + dropped_fraction_exceeded, fragment_cardinality, mod_signature, ms2pip_values, pair_key, + rows_to_drop, MAX_DROPPED_FRACTION, + }; use std::collections::HashMap; #[test] @@ -702,6 +764,42 @@ mod cardinality_tests { assert_eq!(v2, vec![1.0, 0.0]); } + #[test] + fn a_large_unpredicted_fraction_is_a_failure_not_a_warning() { + // docs/31 F10: the only previous guard fired when the library was emptied, so a + // worker returning nothing for a large batch removed those rows plus every + // positional isomer sharing their signature behind one warning. + assert_eq!(dropped_fraction_exceeded(0, 1000), None); + assert_eq!(dropped_fraction_exceeded(0, 0), None); + assert_eq!( + dropped_fraction_exceeded(5, 1000), + None, + "a few misses are fine" + ); + let frac = dropped_fraction_exceeded(500_000, 9_800_000).expect("5% must be refused"); + assert!(frac > MAX_DROPPED_FRACTION); + // The boundary itself is accepted; only exceeding it is refused. + assert_eq!(dropped_fraction_exceeded(20, 1000), None); + assert!(dropped_fraction_exceeded(21, 1000).is_some()); + } + + #[test] + fn a_positional_isomer_still_shares_its_pair_key_by_design() { + // Deliberate over-deletion, documented on `rows_to_drop`: a positional key would + // stop matching a reverse decoy to its target, which is an FDR defect. The test + // pins the trade so a future change has to argue with it. + let a = pair_key(7, 2, "M[Oxidation]PEPTIDEMK"); + let b = pair_key(7, 2, "MPEPTIDEM[Oxidation]K"); + assert_eq!(a, b); + let keys = vec![a.clone(), b, pair_key(7, 3, "M[Oxidation]PEPTIDEMK")]; + let drop = rows_to_drop(&keys, &[0]); + assert_eq!( + drop, + vec![true, true, false], + "the other charge is a different precursor" + ); + } + #[test] fn ms2pip_values_use_one_scale_when_the_model_emitted_charge_2() { let keys = [(b'b', 2, 1), (b'y', 3, 1), (b'y', 3, 2), (b'y', 4, 2)]; diff --git a/rust/mumdia/crates/mumdia/src/stages/prescan.rs b/rust/mumdia/crates/mumdia/src/stages/prescan.rs index b24e9e2..ae42387 100644 --- a/rust/mumdia/crates/mumdia/src/stages/prescan.rs +++ b/rust/mumdia/crates/mumdia/src/stages/prescan.rs @@ -32,7 +32,7 @@ use mumdia_io::report::ArtifactReport; use mumdia_io::table::{write_table, Col, TableFile}; use rayon::prelude::*; use serde_json::json; -use tracing::info; +use tracing::{info, warn}; pub struct PrescanParams<'a> { /// `spectra_ms2` for this run. @@ -290,6 +290,12 @@ pub fn run(p: PrescanParams) -> Result { for (k, v) in per_spectrum { obs.entry(k).or_default().extend(v); } + // The full observed RT-bin range, for candidates whose RT window is unbounded (see + // the screening loop). Empty when no spectrum produced a tag, in which case nothing + // can survive on any path. + let (bin_min, bin_max) = obs.keys().fold((i64::MAX, i64::MIN), |(lo, hi), &(_, b)| { + (lo.min(b), hi.max(b)) + }); info!( cells = obs.len(), spectra = ms2.nrows, @@ -330,6 +336,8 @@ pub fn run(p: PrescanParams) -> Result { let slack = p.cfg.rt_slack_s; let t1 = Instant::now(); + // Candidates screened over the whole gradient because their RT window is unbounded. + let unbounded = std::sync::atomic::AtomicU64::new(0); // Screening is independent per candidate, which is what makes this worth doing in Rust: the // equivalent Python loop was single-threaded and ~40% of the per-file wall clock. let mut surv: Vec<(u32, &str)> = (0..lib.nrows) @@ -345,12 +353,26 @@ pub fn run(p: PrescanParams) -> Result { return None; } let c = cid[i] as usize; - let (lo, hi) = (*lo_by.get(c)?, *hi_by.get(c)?); - if !lo.is_finite() || !hi.is_finite() { - return None; - } - let b0 = ((lo - slack) / bin).floor() as i64; - let b1 = ((hi + slack) / bin).floor() as i64; + // An unbounded RT window is "search the whole gradient", not "cannot be + // screened". `rt_im_train::candidate_window` writes (NaN, -inf, +inf) whenever + // calibration is unavailable and documents the infinite bounds as recall-safe; + // this guard used to read them as a screening failure and drop the candidate, + // so a run with no confident seeds -- the documented FASTA/MS2PIP failure -- + // discarded the ENTIRE library and exited 0 (docs/31 F1). A candidate with no + // run_windows row at all is the same case: unknown RT, not absent evidence. + let (lo, hi) = match (lo_by.get(c), hi_by.get(c)) { + (Some(&l), Some(&h)) => (l, h), + _ => (f64::NEG_INFINITY, f64::INFINITY), + }; + let (b0, b1) = if lo.is_finite() && hi.is_finite() { + ( + ((lo - slack) / bin).floor() as i64, + ((hi + slack) / bin).floor() as i64, + ) + } else { + unbounded.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + (bin_min, bin_max) + }; let m = pmz[i]; for w in 0..win.nrows { if !(w_lo[w] <= m && m < w_hi[w]) { @@ -378,15 +400,39 @@ pub fn run(p: PrescanParams) -> Result { } else { f64::NAN }; + let n_unbounded = unbounded.load(std::sync::atomic::Ordering::Relaxed); info!( screened = lib.nrows, survivors = surv.len(), targets = n_t, decoys = n_d, target_decoy_ratio = ratio, + rt_unbounded = n_unbounded, elapsed_ms = t1.elapsed().as_millis() as u64, "prescan: screened candidates" ); + if n_unbounded > 0 { + warn!( + candidates = n_unbounded, + of = lib.nrows, + "prescan: these candidates have no bounded RT window (no calibration, or no \ + run_windows row) and were screened over the whole gradient, which is the \ + recall-safe reading of the sentinel but a weaker screen" + ); + } + // Screening everything away is not a result. Before docs/31 F1 this was the silent + // outcome of an uncalibrated run: every candidate dropped, a zero-row survivors table, + // exit 0. The single-label bail below cannot see it, because it is gated on a + // survivor count. + if lib.nrows > 0 && surv.is_empty() { + anyhow::bail!( + "prescan screened {} candidates and none survived; a search would have no \ + candidates at all. Check that the observed tag index is not empty (spectra, \ + prescan.tol_da, prescan.top_peaks) and that prescan.anchor_mods matches the \ + modifications in this library", + lib.nrows + ); + } // A one-sided survival means the screen has become label-dependent and the modification's // null is gone. Fail loudly rather than emit a library whose modified q-values cannot be // estimated. diff --git a/rust/mumdia/crates/mumdia/src/stages/quant.rs b/rust/mumdia/crates/mumdia/src/stages/quant.rs index 30a45c2..c99823d 100644 --- a/rust/mumdia/crates/mumdia/src/stages/quant.rs +++ b/rust/mumdia/crates/mumdia/src/stages/quant.rs @@ -454,6 +454,17 @@ fn rollup_protein_bases( pub fn run(p: QuantParams) -> Result<(u64, u64)> { let t0 = Instant::now(); + // No output may be one of the inputs (docs/31 F6). + let inputs = [ + ("--psms-scored", p.psms_scored), + ("--chromatograms", p.chromatograms), + ]; + for out in [Some(p.out_peptide), Some(p.out_protein), p.out_fragment] + .into_iter() + .flatten() + { + mumdia_io::refuse_output_over_input(out, &inputs)?; + } // Identified target PSMs below the peptide q threshold. let ps = TableFile::open(p.psms_scored)?; diff --git a/rust/mumdia/crates/mumdia/src/stages/report.rs b/rust/mumdia/crates/mumdia/src/stages/report.rs index 091aeb3..1eff48e 100644 --- a/rust/mumdia/crates/mumdia/src/stages/report.rs +++ b/rust/mumdia/crates/mumdia/src/stages/report.rs @@ -10,7 +10,7 @@ use std::collections::{HashMap, HashSet}; use std::io::Write; -use anyhow::Result; +use anyhow::{Context, Result}; use mumdia_io::table::TableFile; pub struct ReportParams<'a> { @@ -47,6 +47,41 @@ fn qcell(q: f64) -> String { } } +/// The match-between-runs acceptance basis of every row: the flag and the q it was +/// accepted at, or all-false and all-NaN when the table has never been through +/// `mumdia mbr`. +/// +/// Absent means absent; present means readable. A present column of the wrong type used +/// to be swallowed by an `Err(_) => vec![false; n]` arm, so a boolean written as a +/// nullable dtype or as int8 (which `mbr_worker.py` writes through pandas) silently +/// removed EVERY transfer from both TSVs at exit 0, indistinguishable from an MBR run +/// that transferred nothing while the parquet still showed them (docs/31 F2). This is the +/// absent-versus-malformed rule `quant` and `audit` already follow. +fn transfer_columns(t: &TableFile, n: usize) -> Result<(Vec, Vec)> { + let flag = if t.has_column("is_transferred") { + t.bool("is_transferred").with_context(|| { + "report: `is_transferred` is present but not a boolean column; it is written by \ + `mumdia mbr` and decides which rows are reported" + })? + } else { + vec![false; n] + }; + let q = if t.has_column("transfer_q") { + t.f64("transfer_q") + .with_context(|| "report: `transfer_q` is present but not a float column")? + } else { + vec![f64::NAN; n] + }; + if flag.len() != n || q.len() != n { + anyhow::bail!( + "report: the transfer columns have {} and {} rows against {n} scored rows", + flag.len(), + q.len() + ); + } + Ok((flag, q)) +} + /// Write peptides.tsv + proteins.tsv from a scored PSM table. Returns /// (n_peptides, n_protein_groups) at the FDR threshold. pub fn run(p: ReportParams) -> Result<(u64, u64)> { @@ -65,10 +100,7 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { // this stage filters on, so without this a `mumdia mbr` followed by `mumdia report` // showed no transfers at all. Same contract as `quant`: a transfer has already // passed `mbr.q_transfer`, and a decoy is still never reported. - let is_transferred: Vec = match t.bool("is_transferred") { - Ok(v) => v, - Err(_) => vec![false; n], - }; + let (is_transferred, transfer_q) = transfer_columns(&t, n)?; let accepted = |i: usize, q: &[f64]| label[i] == "target" && (is_transferred[i] || q[i] <= p.q_threshold); // The acceptance basis is exported with every row (`is_transferred`, `transfer_q`), @@ -77,10 +109,6 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { // already passed `mbr.q_transfer` (docs/29 #19). A transfer is a peptide-level // acceptance and not protein-group confidence; a group admitted through one carries // the flag so a reader can tell. - let transfer_q: Vec = match t.f64("transfer_q") { - Ok(v) => v, - Err(_) => vec![f64::NAN; n], - }; let transfer_cells = |i: usize| -> String { if is_transferred[i] { let q = if transfer_q[i].is_nan() { @@ -257,10 +285,7 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { let source = t.u32("source")?; let run_q = t.f64("run_psm_q")?; let n = t.nrows; - let is_transferred: Vec = match t.bool("is_transferred") { - Ok(v) => v, - Err(_) => vec![false; n], - }; + let (is_transferred, transfer_q) = transfer_columns(&t, n)?; let accepted = |i: usize, q: &[f64]| label[i] == "target" && (is_transferred[i] || q[i] <= p.q_threshold); // The acceptance basis is exported with every row (`is_transferred`, `transfer_q`), @@ -269,10 +294,6 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { // already passed `mbr.q_transfer` (docs/29 #19). A transfer is a peptide-level // acceptance and not protein-group confidence; a group admitted through one carries // the flag so a reader can tell. - let transfer_q: Vec = match t.f64("transfer_q") { - Ok(v) => v, - Err(_) => vec![f64::NAN; n], - }; let transfer_cells = |i: usize| -> String { if is_transferred[i] { let q = if transfer_q[i].is_nan() { @@ -769,6 +790,69 @@ mod tests { path } + #[test] + fn a_malformed_transfer_column_is_an_error_not_a_silent_loss_of_every_transfer() { + // docs/31 F2: `mbr_worker.py` writes `is_transferred` through pandas, so a + // nullable dtype or int8 0/1 reaches this stage. The old `Err(_) => vec![false; n]` + // arm turned that into "no transfers", removing every match-between-runs + // identification from both TSVs at exit 0 while the parquet still showed them. + let path = tmp("scored_int_flag.parquet"); + write_table( + &path, + vec![ + Col::Str("peptidoform".into(), vec!["PEPTIDEK".into()]), + Col::I32("charge".into(), vec![2]), + Col::Str("protein".into(), vec!["P1".into()]), + Col::Str("label".into(), vec!["target".into()]), + Col::F64("peptide_q_value".into(), vec![0.5]), + Col::Str("protein_group".into(), vec!["PG1".into()]), + Col::F64("pg_q_value".into(), vec![0.5]), + Col::F64("score".into(), vec![1.0]), + // The flag as an integer, which is what a pandas round-trip can produce. + Col::I32("is_transferred".into(), vec![1]), + Col::F64("transfer_q".into(), vec![0.004]), + ], + ) + .unwrap(); + let out = tmp("report_badflag"); + std::fs::create_dir_all(&out).unwrap(); + let e = run(ReportParams { + scored: &path, + peptide_quant: None, + protein_quant: None, + out_peptides: &format!("{out}/peptides.tsv"), + out_proteins: &format!("{out}/proteins.tsv"), + q_threshold: 0.01, + }) + .unwrap_err() + .to_string(); + assert!(e.contains("is_transferred"), "{e}"); + } + + #[test] + fn a_table_that_never_saw_mbr_reports_no_transfers_without_complaint() { + // The absent case stays absent: `mumdia report` on a plain scored table must not + // start demanding columns only `mumdia mbr` writes. + let scored = scored_table(); + let out = tmp("report_nombr"); + std::fs::create_dir_all(&out).unwrap(); + let (n_pep, _) = run(ReportParams { + scored: &scored, + peptide_quant: None, + protein_quant: None, + out_peptides: &format!("{out}/peptides.tsv"), + out_proteins: &format!("{out}/proteins.tsv"), + q_threshold: 0.01, + }) + .unwrap(); + assert_eq!(n_pep, 2); + let text = std::fs::read_to_string(format!("{out}/peptides.tsv")).unwrap(); + for line in text.lines().skip(1) { + let c: Vec<&str> = line.split('\t').collect(); + assert_eq!((c[7], c[8]), ("false", "")); + } + } + #[test] fn report_exports_the_transfer_basis_and_a_tighter_threshold_does_not_revoke_it() { // docs/29 #19: the acceptance rule admits any target flagged `is_transferred` diff --git a/rust/mumdia/crates/mumdia/src/stages/rescore.rs b/rust/mumdia/crates/mumdia/src/stages/rescore.rs index 9f7f835..351b876 100644 --- a/rust/mumdia/crates/mumdia/src/stages/rescore.rs +++ b/rust/mumdia/crates/mumdia/src/stages/rescore.rs @@ -112,6 +112,14 @@ pub fn run(p: RescoreParams) -> Result { if p.competed.is_empty() { anyhow::bail!("rescore requires at least one competed input"); } + // `--out` must not be one of the competed tables: they are read before the scored + // table is published, so writing over one replaces it and exits 0 (docs/31 F6). + let competed_inputs: Vec<(&str, &str)> = p + .competed + .iter() + .map(|c| ("--competed", c.as_str())) + .collect(); + mumdia_io::refuse_output_over_input(p.out, &competed_inputs)?; if p.cfg.folds < 2 { anyhow::bail!("rescore.folds must be >= 2 for out-of-fold scoring"); } diff --git a/rust/mumdia/crates/mumdia/src/stages/rt_im_train.rs b/rust/mumdia/crates/mumdia/src/stages/rt_im_train.rs index f248bff..250cd01 100644 --- a/rust/mumdia/crates/mumdia/src/stages/rt_im_train.rs +++ b/rust/mumdia/crates/mumdia/src/stages/rt_im_train.rs @@ -99,6 +99,23 @@ fn candidate_window(calibrated_rt: Option, width: Option) -> (f64, f64 pub fn run(p: RtImTrainParams) -> Result { let t0 = Instant::now(); + // `--out` must not be one of this stage's own inputs: every input is read + // before the output is published, so writing over one replaces it and exits 0 + // (docs/31 F6). The shared guard existed and was wired into two stages. + mumdia_io::refuse_output_over_input( + p.out_windows, + &[ + ("--seed-psms", p.seed_psms), + ("--lib-precursors", p.library_precursors), + ], + )?; + mumdia_io::refuse_output_over_input( + p.out_cal, + &[ + ("--seed-psms", p.seed_psms), + ("--lib-precursors", p.library_precursors), + ], + )?; let holdout_frac = p.cfg.window_holdout_frac; if !(0.0..=0.9).contains(&holdout_frac) { @@ -357,8 +374,17 @@ pub fn run(p: RtImTrainParams) -> Result { Vec::with_capacity(n), Vec::with_capacity(n), ); + let mut n_nonfinite_irt = 0u64; for i in 0..n { - let calibrated_rt = calibration_available.then(|| predict(irt[i] as f64)); + // A row whose library iRT is not finite has no calibrated RT, which is the + // documented "search the whole gradient" sentinel rather than an arithmetic + // accident. The parquet reader maps a null f32 to NaN, so one failed prediction + // in an imported library reaches here (docs/31 F4). + let usable_irt = (irt[i] as f64).is_finite(); + if !usable_irt { + n_nonfinite_irt += 1; + } + let calibrated_rt = (calibration_available && usable_irt).then(|| predict(irt[i] as f64)); let width = calibrated_rt.map(|cal| match &adaptive { Some((rt_min, span, widths)) => { let nb = widths.len(); @@ -457,6 +483,10 @@ pub fn run(p: RtImTrainParams) -> Result { stats.insert("n_train".to_string(), json!(n_train)); stats.insert("w_rt".to_string(), json!(w_rt)); stats.insert("calibration_status".to_string(), json!(status)); + stats.insert( + "candidates_without_finite_irt".to_string(), + json!(n_nonfinite_irt), + ); ArtifactReport { logical_name: artifact::RUN_WINDOWS.0.to_string(), schema_name: artifact::RUN_WINDOWS.0.to_string(), @@ -471,10 +501,18 @@ pub fn run(p: RtImTrainParams) -> Result { } .write_for(p.out_windows)?; + if n_nonfinite_irt > 0 { + tracing::warn!( + candidates = n_nonfinite_irt, + of = rows, + "rt-im-train: these candidates have no finite library iRT, so they get the unbounded RT window rather than a calibrated one; a null predicted_irt reads as NaN (docs/31 F4)" + ); + } info!( rows, w_rt = ?w_rt, status, + without_finite_irt = n_nonfinite_irt, elapsed_ms = elapsed, "rt-im-train: done" ); diff --git a/rust/mumdia/crates/mumdia/src/stages/search_seed.rs b/rust/mumdia/crates/mumdia/src/stages/search_seed.rs index 710ef78..1c3880c 100644 --- a/rust/mumdia/crates/mumdia/src/stages/search_seed.rs +++ b/rust/mumdia/crates/mumdia/src/stages/search_seed.rs @@ -44,6 +44,17 @@ struct Best { pub fn run(p: SearchSeedParams) -> Result { let t0 = Instant::now(); + // `--out` must not be one of this stage's own inputs: every input is read + // before the output is published, so writing over one replaces it and exits 0 + // (docs/31 F6). The shared guard existed and was wired into two stages. + mumdia_io::refuse_output_over_input( + p.out, + &[ + ("--ms2", p.ms2), + ("--lib-precursors", p.library_precursors), + ("--lib-fragments", p.library_fragments), + ], + )?; // See extract: the bucketed index is dead weight on the fragindex path. let build_bucketed = !matches!(p.cfg.matcher, MatcherKind::Fragindex); let lib = Library::load_with( diff --git a/scripts/nn_rescore_worker.py b/scripts/nn_rescore_worker.py index d8650b6..fa9cce1 100644 --- a/scripts/nn_rescore_worker.py +++ b/scripts/nn_rescore_worker.py @@ -18,10 +18,12 @@ MEMORY (multi-run / large PINs): two feature backends behind one accessor. - in-memory (default for PINs <= MUMDIA_NN_STREAM_GB, 4 GB): the full standardised - feature matrix is held in RAM (median/IQR standardisation). + feature matrix is held in RAM. Mean/std standardisation, the same as the streaming + backend: one transform for every backend and every handoff, so a score does not + depend on which one ran (docs/31 F5). - streaming memmap (large PINs, or MUMDIA_NN_STREAM=1): the PIN is read ONCE in - chunks into a disk-backed float32 memmap (mean/std standardisation accumulated - in the same pass); training and scoring then draw MINIBATCHES indexed into the + chunks into a disk-backed float32 memmap (the same mean/std, accumulated in the + same pass); training and scoring then draw MINIBATCHES indexed into the memmap, so peak RAM is one batch + per-row metadata, NOT the whole matrix. This is what makes combining many runs into one rescoring tractable: the full PIN never lives in RAM at once. @@ -165,7 +167,21 @@ def folds_for(peptides, fold_keys, folds, off=0): serve the chunked streaming backend as well as the two whole-table ones. """ if fold_keys is not None: - return (fold_keys[off:off + len(peptides)] % folds).astype(np.int16) + want = len(peptides) + got = fold_keys[off:off + want] + # Numpy slicing past the end returns a SHORT array rather than raising, and a + # short `fold` puts the tail rows in no fold at all: `np.where(fold == f)` cannot + # index them, they are never scored, and they keep the zero initialiser, which the + # final rank-normalisation turns into a plausible tied mid-rank score. The Rust + # caller's completeness contract is satisfied by that, so `rescore.strict` does not + # catch it either (docs/31 F5). + if len(got) != want: + raise SystemExit( + "MUMDIA_NN_FOLD_KEYS has %d rows but the PIN needs at least %d " + "(rows %d..%d); the companion table does not belong to this PIN, and " + "folding on a truncated one would leave the tail unscored." + % (len(fold_keys), off + want, off, off + want)) + return (got % folds).astype(np.int16) return np.array([peptide_fold(x, folds) for x in peptides], np.int16) @@ -333,7 +349,10 @@ def main(): # times smaller than the equivalent text, so the raw file size would understate the # memory a full read actually needs. _md = pq.read_metadata(pin_path) - _nf_guess = max(1, _md.num_columns - 3) # minus SpecId / Label / Peptide + # Count the feature columns by name. Subtracting a hardcoded 3 undercounted the + # non-feature columns, of which NON_FEATURE lists 7, so the estimate that picks the + # backend was biased upward and could stream a matrix that fits (docs/31 F5). + _nf_guess = max(1, sum(1 for c in pq.read_schema(pin_path).names if c not in NON_FEATURE)) filesize = int(_md.num_rows) * _nf_guess * 4 stream_gb = env_f("MUMDIA_NN_STREAM_GB", 4) stream = stream_env in ("1", "on", "true") or ( @@ -475,10 +494,18 @@ def main(): pin[feat_cols].to_numpy(np.float32), nan=0.0, posinf=0.0, neginf=0.0 ) del pin - med = np.median(X, axis=0) - iqr = np.subtract(*np.percentile(X, [75, 25], axis=0)) - iqr[iqr == 0] = 1.0 - Xs = np.clip((X - med) / iqr, -8, 8).astype(np.float32) + # Mean/std, the same transform as the parquet and streaming backends. + # + # This path used median/IQR while the other two used mean/std, so the same PSM + # pool produced different scores depending on `rescore.handoff` and on which + # side of MUMDIA_NN_STREAM_GB the matrix fell: a 3.99 GB PIN was standardised + # one way and a 4.01 GB PIN the other, with nothing in the log to say which + # (docs/31 F5). Converging on mean/std leaves the shipped default (parquet) + # and every published benchmark unchanged, and only moves this legacy path. + mean = X.mean(axis=0, dtype=np.float64).astype(np.float32) + std = X.std(axis=0, dtype=np.float64).astype(np.float32) + std[std == 0] = 1.0 + Xs = np.clip((X - mean) / std, -8, 8).astype(np.float32) del X n = len(y) get = lambda idx: Xs[idx] diff --git a/tests/python/test_nn_rescore_worker.py b/tests/python/test_nn_rescore_worker.py index 79afd07..5fa3130 100644 --- a/tests/python/test_nn_rescore_worker.py +++ b/tests/python/test_nn_rescore_worker.py @@ -237,6 +237,32 @@ def test_explicit_fold_keys_pair_a_target_with_its_decoy(): assert hashed[0] != hashed[1] +def test_a_short_fold_key_file_is_refused_rather_than_leaving_rows_unfolded(): + """A fold-key companion shorter than the PIN must stop the run. + + Numpy slicing past the end returns a SHORT array instead of raising, and a short + `fold` puts the tail rows in no fold at all: `np.where(fold == f)` cannot reach + them, they are never scored, and they keep the zero initialiser, which the final + rank-normalisation turns into a plausible tied mid-rank score. The Rust caller's + completeness contract is satisfied by that, so `rescore.strict` does not catch it + either (docs/31 F5). + """ + w = _import_worker() + np = pytest.importorskip("numpy") + + keys = np.arange(6, dtype=np.uint32) + peptides = [f"p{i}" for i in range(10)] + with pytest.raises(SystemExit) as exc: + w.folds_for(peptides, keys, 3) + assert "MUMDIA_NN_FOLD_KEYS" in str(exc.value) + # A companion that covers the rows is unaffected, at an offset too. + full = np.arange(10, dtype=np.uint32) + assert len(w.folds_for(peptides, full, 3)) == 10 + assert len(w.folds_for(peptides[7:], full, 3, off=7)) == 3 + with pytest.raises(SystemExit): + w.folds_for(peptides[7:], keys, 3, off=7) + + def test_fold_keys_respect_the_streaming_row_offset(): """The chunked backend passes a flat row offset; the keys must be sliced by it.