From fef5a9572ca804d458c72a67d72749f407d763c9 Mon Sep 17 00:00:00 2001
From: Vladimir Shitov <vladimirshitov98@gmail.com>
Date: Sun, 25 Oct 2026 07:17:38 +0100
Subject: [PATCH 1/4] correlation: rank CSC matrices, so CSR predictions get
 the right Spearman scores

dynutils::spearman_rank_sparse() ranks along the compressed dimension, which is a
column only for CSC. ss_opm writes its predictions as CSR, so its per-cell and
per-gene Spearman correlations were computed on ranks along the wrong axis on every
dataset. Convert the solution and the prediction to CSC after reading them, and add
a component test that scores a prediction as CSR and as CSC against dense ranks.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  6 ++
 src/metrics/correlation/config.vsh.yaml       |  4 ++
 src/metrics/correlation/script.R              | 11 ++--
 .../correlation/test_prediction_layout.py     | 57 +++++++++++++++++++
 4 files changed, 73 insertions(+), 5 deletions(-)
 create mode 100644 src/metrics/correlation/test_prediction_layout.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 583e59ca..3310492d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,9 @@
+# task_predict_modality 0.2.1
+
+## BUG FIXES
+
+* `correlation`: Convert the solution and the prediction to CSC before correlating them. `dynutils::spearman_rank_sparse()` ranks along the compressed dimension, which is a column only for CSC, so a CSR matrix (read as a `dgRMatrix`) was ranked per feature for the per-cell Spearman correlation and per cell for the per-gene one. In 0.2.0 this hit `mean_spearman_per_cell` and `mean_spearman_per_gene` of `ss_opm`, which writes CSR, on every dataset, and of every method on GEX -> ATAC, whose solution holds the ATAC layer as CSR (`random_predict` also writes CSR there). All other scores are unchanged. A new component test scores a prediction as CSR and as CSC and checks both against dense ranks.
+
 # task_predict_modality 0.2.0
 
 Adds eight methods, selects the ATAC target peaks by `hvg_score`, and fixes the bugs surfaced by the first full benchmark runs.
diff --git a/src/metrics/correlation/config.vsh.yaml b/src/metrics/correlation/config.vsh.yaml
index 4c1db788..06980942 100644
--- a/src/metrics/correlation/config.vsh.yaml
+++ b/src/metrics/correlation/config.vsh.yaml
@@ -73,6 +73,10 @@ info:
 resources:
   - type: r_script
     path: script.R
+# __merge__ appends to the api's test_resources
+test_resources:
+  - type: python_script
+    path: test_prediction_layout.py
 engines:
   - type: docker
     image: openproblems/base_r:1
diff --git a/src/metrics/correlation/script.R b/src/metrics/correlation/script.R
index 977a6fdd..f285aa88 100644
--- a/src/metrics/correlation/script.R
+++ b/src/metrics/correlation/script.R
@@ -30,15 +30,16 @@ expect_true(
 
 cat("Computing correlation metrics\n")
 # Wrangle data
-tv <- ad_sol$layers[["normalized"]]
-pv <- ad_pred$layers[["normalized"]]
+# dynutils' spearman_rank_sparse() ranks along the compressed dimension, so a CSR (dgRMatrix) prediction
+# would be ranked along the wrong axis
+tv <- methods::as(ad_sol$layers[["normalized"]], "CsparseMatrix")
+pv <- methods::as(ad_pred$layers[["normalized"]], "CsparseMatrix")
 
 # score non-finite predictions as zero, like the zero-variance case below
-pv_values <- if (methods::is(pv, "sparseMatrix")) pv@x else pv
-non_finite <- !is.finite(pv_values)
+non_finite <- !is.finite(pv@x)
 if (any(non_finite)) {
   cat("Prediction contains", sum(non_finite), "non-finite values, scoring them as 0\n")
-  if (methods::is(pv, "sparseMatrix")) pv@x[non_finite] <- 0 else pv[non_finite] <- 0
+  pv@x[non_finite] <- 0
 }
 
 # precompute sds
diff --git a/src/metrics/correlation/test_prediction_layout.py b/src/metrics/correlation/test_prediction_layout.py
new file mode 100644
index 00000000..7d694bc0
--- /dev/null
+++ b/src/metrics/correlation/test_prediction_layout.py
@@ -0,0 +1,57 @@
+import subprocess
+
+import anndata as ad
+import numpy as np
+import scipy.sparse
+from scipy.stats import rankdata
+
+## VIASH START
+meta = {
+    "executable": "target/executable/metrics/correlation/correlation",
+    "resources_dir": ".",
+}
+## VIASH END
+
+# A prediction must get the same scores whether its layer is stored as CSR or CSC, and its per-cell and per-gene
+# Spearman correlations must match those of dense ranks. dynutils' spearman_rank_sparse() ranks along the compressed
+# dimension, so a CSR prediction used to be ranked along the wrong axis.
+input_dir = f"{meta['resources_dir']}/resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/swap"
+solution = ad.read_h5ad(f"{input_dir}/test_mod2.h5ad")
+prediction = ad.read_h5ad(f"{input_dir}/prediction.h5ad")
+
+
+def score(layout):
+    prediction.layers["normalized"] = scipy.sparse.csr_matrix(prediction.layers["normalized"]).asformat(layout)
+    prediction.write_h5ad(f"prediction_{layout}.h5ad")
+    subprocess.run(
+        [
+            meta["executable"],
+            "--input_test_mod2", f"{input_dir}/test_mod2.h5ad",
+            "--input_prediction", f"prediction_{layout}.h5ad",
+            "--output", f"scores_{layout}.h5ad",
+        ],
+        check=True,
+    )
+    return ad.read_h5ad(f"scores_{layout}.h5ad")
+
+
+def paired_spearman(truth, predicted, axis):
+    """Spearman correlation of matching rows (axis=1) or columns (axis=0); constant vectors score 0."""
+    truth_ranks, predicted_ranks = (rankdata(values, axis=axis) for values in (truth, predicted))
+    truth_ranks -= truth_ranks.mean(axis=axis, keepdims=True)
+    predicted_ranks -= predicted_ranks.mean(axis=axis, keepdims=True)
+    norms = np.sqrt((truth_ranks**2).sum(axis=axis) * (predicted_ranks**2).sum(axis=axis))
+    covariance = (truth_ranks * predicted_ranks).sum(axis=axis)
+    return np.divide(covariance, norms, out=np.zeros_like(covariance), where=norms > 0)
+
+
+scores = {layout: score(layout) for layout in ["csr", "csc"]}
+truth = solution.layers["normalized"].toarray()
+predicted = prediction.layers["normalized"].toarray()
+expected_spearman = {"obs": paired_spearman(truth, predicted, axis=1), "var": paired_spearman(truth, predicted, axis=0)}
+
+for layout, layout_scores in scores.items():
+    print(f"{layout}:", dict(zip(layout_scores.uns["metric_ids"], np.round(layout_scores.uns["metric_values"], 4))))
+    np.testing.assert_allclose(layout_scores.uns["metric_values"], scores["csc"].uns["metric_values"], rtol=0, atol=1e-10)
+    for annotation, expected in expected_spearman.items():
+        np.testing.assert_allclose(getattr(layout_scores, annotation)["spearman"], expected, rtol=0, atol=1e-8)

From 0d84b482bcf95abf2f07bf033e30c754feb1b372 Mon Sep 17 00:00:00 2001
From: Vladimir Shitov <vladimirshitov98@gmail.com>
Date: Sun, 25 Oct 2026 07:17:38 +0100
Subject: [PATCH 2/4] guanlab_dengkw_pm: solve the kernel ridge regression with
 numpy's 64-bit-integer LAPACK

scipy's 32-bit-integer OpenBLAS Cholesky reports the kernel matrix as not positive
definite once a half of the batches holds more than ~30-46k cells, so the method
failed on pbmc_multiome ATAC -> GEX in 0.2.0. Also ask for veryhighmem (118 GB peak)
and veryhightime (GEX -> ATAC hit hightime's 8 hours).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  2 ++
 src/methods/guanlab_dengkw_pm/config.vsh.yaml |  2 +-
 src/methods/guanlab_dengkw_pm/script.py       | 21 +++++++------------
 3 files changed, 11 insertions(+), 14 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3310492d..39ee6825 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,8 @@
 
 * `correlation`: Convert the solution and the prediction to CSC before correlating them. `dynutils::spearman_rank_sparse()` ranks along the compressed dimension, which is a column only for CSC, so a CSR matrix (read as a `dgRMatrix`) was ranked per feature for the per-cell Spearman correlation and per cell for the per-gene one. In 0.2.0 this hit `mean_spearman_per_cell` and `mean_spearman_per_gene` of `ss_opm`, which writes CSR, on every dataset, and of every method on GEX -> ATAC, whose solution holds the ATAC layer as CSR (`random_predict` also writes CSR there). All other scores are unchanged. A new component test scores a prediction as CSR and as CSC and checks both against dense ranks.
 
+* `guanlab_dengkw_pm`: Solve the kernel ridge regression with `numpy.linalg.solve()` instead of `scipy.linalg.cho_factor()`/`cho_solve()`. The Cholesky of the 32-bit-integer OpenBLAS that scipy's wheels link reports the kernel matrix as not positive definite once a half of the batches holds more than ~30-46k cells (`LinAlgError: Internal potrf return info = [49313]` on `pbmc_multiome` ATAC -> GEX in 0.2.0, where the halves hold 56k-74k cells); numpy's wheels link a 64-bit-integer OpenBLAS whose LU solve handles these systems. Also ask for `veryhighmem` and `veryhightime` rather than `highmem` and `hightime`: ATAC -> GEX peaks at 112-118 GB, and GEX -> ATAC, which 0.2.0 stopped at `hightime`'s 8 hours, takes 6.5-8 hours on 30 cores.
+
 # task_predict_modality 0.2.0
 
 Adds eight methods, selects the ATAC target peaks by `hvg_score`, and fixes the bugs surfaced by the first full benchmark runs.
diff --git a/src/methods/guanlab_dengkw_pm/config.vsh.yaml b/src/methods/guanlab_dengkw_pm/config.vsh.yaml
index 4f1fc942..4bad9a8d 100644
--- a/src/methods/guanlab_dengkw_pm/config.vsh.yaml
+++ b/src/methods/guanlab_dengkw_pm/config.vsh.yaml
@@ -48,4 +48,4 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [hightime, highmem, highcpu]
+      label: [veryhightime, veryhighmem, highcpu]
diff --git a/src/methods/guanlab_dengkw_pm/script.py b/src/methods/guanlab_dengkw_pm/script.py
index 1a0bdd81..98649d97 100644
--- a/src/methods/guanlab_dengkw_pm/script.py
+++ b/src/methods/guanlab_dengkw_pm/script.py
@@ -1,7 +1,6 @@
 import sys
 import anndata as ad
 import numpy as np
-from scipy.linalg import cho_factor, cho_solve
 from scipy.sparse import csc_matrix
 from sklearn.decomposition import TruncatedSVD
 from sklearn.gaussian_process.kernels import RBF
@@ -32,22 +31,18 @@
 # Workaround for sklearn's KernelRidge: this computes
 # KernelRidge(alpha=alpha, kernel=kernel).fit(X, y).predict(X_test), but solves the system itself.
 #
-# KernelRidge solves (K + alpha * I) w = y with scipy.linalg.solve(assume_a="pos"). In this image
-# (scipy 1.17 and 1.18 with OpenBLAS) that call raises a MemoryError or segfaults once a half of the
-# batches holds more than ~30-40k cells, no matter how much memory the job gets (200 GB was not enough
-# for 40k cells, which need ~13 GB). So the method failed on every dataset except bmmc_multiome.
+# KernelRidge solves (K + alpha * I) w = y with scipy.linalg.solve(assume_a="pos"). The Cholesky of the
+# 32-bit-integer OpenBLAS that scipy's wheels link breaks down on large systems: solve() raised a MemoryError
+# or segfaulted, and cho_factor() reports K + alpha * I as not positive definite ("potrf return info") for
+# every RBF system of 46k-60k cells tried (30k still works) and for every half of the batches of pbmc_multiome
+# (56k-74k cells). numpy's wheels link a 64-bit-integer OpenBLAS (scipy-openblas64), whose LU solve handles
+# them; where scipy's Cholesky works, the coefficients agree to 6e-13.
 #
-# K + alpha * I is symmetric positive definite, so solve() does a Cholesky factorisation anyway.
-# cho_factor()/cho_solve() run that same factorisation directly and handle 46k cells in about a minute;
-# the predictions are identical to KernelRidge's. They also factorise K in place, whereas KernelRidge
-# copies it first, which saves one n x n float64 matrix of peak memory.
-#
-# Switch back to KernelRidge once scipy.linalg.solve() handles matrices of this size again.
+# Switch back to KernelRidge once scipy ships a 64-bit-integer LAPACK.
 def kernel_ridge_predict(kernel, alpha, X, y, X_test):
     K = kernel(X)
     K.flat[::K.shape[0] + 1] += alpha
-    # K is symmetric, so its transpose is a Fortran-ordered view LAPACK can factorise in place
-    dual_coef = cho_solve(cho_factor(K.T, overwrite_a=True), y)
+    dual_coef = np.linalg.solve(K, y)
     return kernel(X_test, X) @ dual_coef
 
 

From a79b960d92f97fda3b3799cc5cdf37ecfa32016a Mon Sep 17 00:00:00 2001
From: Vladimir Shitov <vladimirshitov98@gmail.com>
Date: Sun, 25 Oct 2026 07:18:07 +0100
Subject: [PATCH 3/4] simple_mlp: keep the inputs sparse and densify one batch
 at a time

split() densified whole training and validation folds (143 GB for the ATAC input of
pbmc_multiome), so ATAC -> GEX failed on it in every run since run_2026-08-14.
utils.SparseRows serves one batch at a time as dense tensors through a BatchSampler,
drawing the same batches as DataLoader(shuffle=True) did; the unused counts layers
are dropped after reading.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  2 +
 src/methods/simple_mlp/resources/utils.py     | 41 +++++++++++++++++--
 .../simple_mlp/simple_mlp_predict/script.py   | 10 ++---
 .../simple_mlp/simple_mlp_train/script.py     | 22 ++++------
 4 files changed, 51 insertions(+), 24 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 39ee6825..62809422 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,8 @@
 
 * `guanlab_dengkw_pm`: Solve the kernel ridge regression with `numpy.linalg.solve()` instead of `scipy.linalg.cho_factor()`/`cho_solve()`. The Cholesky of the 32-bit-integer OpenBLAS that scipy's wheels link reports the kernel matrix as not positive definite once a half of the batches holds more than ~30-46k cells (`LinAlgError: Internal potrf return info = [49313]` on `pbmc_multiome` ATAC -> GEX in 0.2.0, where the halves hold 56k-74k cells); numpy's wheels link a 64-bit-integer OpenBLAS whose LU solve handles these systems. Also ask for `veryhighmem` and `veryhightime` rather than `highmem` and `hightime`: ATAC -> GEX peaks at 112-118 GB, and GEX -> ATAC, which 0.2.0 stopped at `hightime`'s 8 hours, takes 6.5-8 hours on 30 cores.
 
+* `simple_mlp_train`, `simple_mlp_predict`: Keep the inputs and targets sparse and densify one batch at a time (`utils.SparseRows` + `utils.batch_loader()`), and drop the unused `counts` layers after reading. `split()` densified whole folds, 143 GB for the ATAC input of `pbmc_multiome`, so ATAC -> GEX failed on it in every run since run_2026-08-14, 0.2.0 included. Batches are drawn as before.
+
 # task_predict_modality 0.2.0
 
 Adds eight methods, selects the ATAC target peaks by `hvg_score`, and fixes the bugs surfaced by the first full benchmark runs.
diff --git a/src/methods/simple_mlp/resources/utils.py b/src/methods/simple_mlp/resources/utils.py
index f37e54c2..47ad338c 100644
--- a/src/methods/simple_mlp/resources/utils.py
+++ b/src/methods/simple_mlp/resources/utils.py
@@ -1,6 +1,9 @@
 import numpy as np
+import scipy.sparse
+import torch
 import yaml
 from collections import namedtuple
+from torch.utils.data import BatchSampler, DataLoader, Dataset, RandomSampler, SequentialSampler
 
 
 def to_site_donor(data):
@@ -26,17 +29,47 @@ def split(tr1, tr2, fold):
 
     maskr = ~mask
 
-    Xt = tr1[mask].layers["normalized"].toarray()
-    X = tr1[maskr].layers["normalized"].toarray()
+    # kept sparse: densified, the NeurIPS 2022 ATAC input alone needs 143 GB; see batch_loader
+    Xt = float32_rows(tr1.layers["normalized"], mask)
+    X = float32_rows(tr1.layers["normalized"], maskr)
 
-    yt = tr2[mask].layers["normalized"].toarray()
-    y = tr2[maskr].layers["normalized"].toarray()
+    yt = float32_rows(tr2.layers["normalized"], mask)
+    y = float32_rows(tr2.layers["normalized"], maskr)
 
     print(f"{X.shape}, {y.shape}, {Xt.shape}, {yt.shape}")
 
     return X,y,Xt,yt
 
 
+def float32_rows(matrix, rows):
+    """The selected rows of a sparse matrix as float32 CSR, the precision the MLP trains in."""
+    return scipy.sparse.csr_matrix(matrix[rows], dtype=np.float32)
+
+
+class SparseRows(Dataset):
+    """Rows of one or more sparse matrices, served as dense float32 tensors a whole batch at a time."""
+
+    def __init__(self, *matrices):
+        self.matrices = [scipy.sparse.csr_matrix(matrix, dtype=np.float32) for matrix in matrices]
+
+    def __len__(self):
+        return self.matrices[0].shape[0]
+
+    def __getitem__(self, rows):
+        return tuple(torch.from_numpy(matrix[rows].toarray()) for matrix in self.matrices)
+
+
+def batch_loader(dataset, batch_size, shuffle, drop_last, num_workers):
+    """DataLoader that hands `dataset` the row indices of a whole batch, so only one batch is densified at a time."""
+    rows = RandomSampler(dataset) if shuffle else SequentialSampler(dataset)
+    return DataLoader(
+        dataset,
+        sampler=BatchSampler(rows, batch_size=batch_size, drop_last=drop_last),
+        batch_size=None,
+        num_workers=num_workers,
+    )
+
+
 def load_yaml(path):
     with open(path) as f:
         x = yaml.safe_load(f)
diff --git a/src/methods/simple_mlp/simple_mlp_predict/script.py b/src/methods/simple_mlp/simple_mlp_predict/script.py
index aaf31647..1a8b8b5e 100644
--- a/src/methods/simple_mlp/simple_mlp_predict/script.py
+++ b/src/methods/simple_mlp/simple_mlp_predict/script.py
@@ -4,7 +4,6 @@
 from scipy.sparse import csc_matrix
 import anndata as ad
 import torch
-from torch.utils.data import TensorDataset,DataLoader
 
 ## VIASH START
 par = {
@@ -68,10 +67,7 @@ def _predict(model,dl):
     ymean = torch.from_numpy(ymean).float()
     yaml_path=f"{resources_dir}/yaml/mlp_{task}.yaml"
     config = utils.load_yaml(yaml_path)
-    X = input_test_mod1.layers["normalized"].toarray()
-    X = torch.from_numpy(X).float()
-
-    te_ds = TensorDataset(X)
+    te_ds = utils.SparseRows(input_test_mod1.layers["normalized"])
 
     yp = 0
     for fold in folds:
@@ -81,12 +77,12 @@ def _predict(model,dl):
         ckpt = glob(load_path)[0]
         model_inf = MLP.load_from_checkpoint(
             ckpt,
-            in_dim=X.shape[1],
+            in_dim=input_test_mod1.n_vars,
             out_dim=input_train_mod2.n_vars,
             ymean=ymean,
             config=config
         )
-        te_loader = DataLoader(
+        te_loader = utils.batch_loader(
             te_ds,
             batch_size=config.batch_size,
             num_workers=0,
diff --git a/src/methods/simple_mlp/simple_mlp_train/script.py b/src/methods/simple_mlp/simple_mlp_train/script.py
index 7f90f4e0..7f35e36e 100644
--- a/src/methods/simple_mlp/simple_mlp_train/script.py
+++ b/src/methods/simple_mlp/simple_mlp_train/script.py
@@ -8,7 +8,6 @@
 
 import torch
 import pytorch_lightning as pl
-from torch.utils.data import TensorDataset, DataLoader
 from pytorch_lightning.callbacks import ModelCheckpoint
 from pytorch_lightning.loggers import TensorBoardLogger,WandbLogger
 
@@ -37,24 +36,18 @@
 
 def _train(X, y, Xt, yt, logger, config, num_workers):
 
-    X = torch.from_numpy(X).float()
-    y = torch.from_numpy(y).float()
-    ymean = torch.mean(y, dim=0, keepdim=True)
+    ymean = torch.from_numpy(np.asarray(y.mean(axis=0), dtype=np.float32))
 
-    tr_ds = TensorDataset(X,y)
-    tr_loader = DataLoader(
-        tr_ds,
+    tr_loader = utils.batch_loader(
+        utils.SparseRows(X, y),
         batch_size=config.batch_size,
         num_workers=num_workers,
         shuffle=True,
         drop_last=True
     )
 
-    Xt = torch.from_numpy(Xt).float()
-    yt = torch.from_numpy(yt).float()
-    te_ds = TensorDataset(Xt,yt)
-    te_loader = DataLoader(
-        te_ds,
+    te_loader = utils.batch_loader(
+        utils.SparseRows(Xt, yt),
         batch_size=config.batch_size,
         num_workers=num_workers,
         shuffle=False,
@@ -83,7 +76,7 @@ def _train(X, y, Xt, yt, logger, config, num_workers):
     yp = trainer.predict(net, te_loader, ckpt_path='best')
     yp = torch.cat(yp, dim=0)
 
-    score = ((yp-yt)**2).mean()**0.5
+    score = ((yp - torch.from_numpy(yt.toarray()).float())**2).mean()**0.5
     print(f"VALID RMSE {score:.3f}")
     del trainer
     return score,yp.detach().numpy()
@@ -92,6 +85,9 @@ def _train(X, y, Xt, yt, logger, config, num_workers):
 
 input_train_mod1 = ad.read_h5ad(par['input_train_mod1'])
 input_train_mod2 = ad.read_h5ad(par['input_train_mod2'])
+# only the normalized layers are used, and the counts take as much memory again
+for adata in (input_train_mod1, input_train_mod2):
+    del adata.layers["counts"]
 
 mod_1 = input_train_mod1.uns["modality"]
 mod_2 = input_train_mod2.uns["modality"]

From cbc458904cdc9501c7e1d65618952ccf80939004 Mon Sep 17 00:00:00 2001
From: Vladimir Shitov <vladimirshitov98@gmail.com>
Date: Sun, 25 Oct 2026 09:52:28 +0100
Subject: [PATCH 4/4] Seed every stochastic method and control

Rerunning the benchmark did not reproduce its predictions: guanlab_dengkw_pm's and
knnr_py's truncated SVDs, the lmds landmarks of lm and knnr_r, random_predict, the
scvi-tools and pynndescent steps of cellmapper, and the training of simple_mlp, novel,
babel and senkin_tmp's networks were unseeded. src/utils/seed.py seeds python, numpy,
torch and tensorflow and makes them use deterministic GPU algorithms (also for the
already seeded ss_opm, scipenn and scbutterfly); the R components call set.seed().
Methods without one get a --seed argument (default 0).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  4 +++
 .../random_predict/config.vsh.yaml            |  5 ++++
 src/control_methods/random_predict/script.R   |  5 +++-
 src/methods/babel/babel_train/config.vsh.yaml |  5 ++++
 src/methods/babel/babel_train/script.py       |  4 +++
 src/methods/cellmapper_linear/config.vsh.yaml |  5 ++++
 src/methods/cellmapper_linear/script.py       | 11 +++++++-
 src/methods/cellmapper_scvi/config.vsh.yaml   |  5 ++++
 src/methods/cellmapper_scvi/script.py         |  7 ++++-
 src/methods/guanlab_dengkw_pm/config.vsh.yaml |  2 +-
 src/methods/guanlab_dengkw_pm/script.py       |  4 +--
 src/methods/knnr_py/config.vsh.yaml           |  4 +++
 src/methods/knnr_py/script.py                 |  3 ++-
 src/methods/knnr_r/config.vsh.yaml            |  4 +++
 src/methods/knnr_r/script.R                   |  5 +++-
 src/methods/lm/config.vsh.yaml                |  4 +++
 src/methods/lm/script.R                       |  5 +++-
 src/methods/novel/novel_train/config.vsh.yaml |  1 +
 src/methods/novel/novel_train/script.py       |  3 +++
 src/methods/scbutterfly/butterfly_common.py   |  6 +++++
 .../scbutterfly_predict/config.vsh.yaml       |  1 +
 .../scbutterfly_train/config.vsh.yaml         |  1 +
 src/methods/scipenn/config.vsh.yaml           |  1 +
 src/methods/scipenn/script.py                 |  4 +--
 .../senkin_tmp_train/config.vsh.yaml          |  1 +
 .../senkin_tmp/senkin_tmp_train/script.py     |  3 +++
 .../simple_mlp_train/config.vsh.yaml          |  5 ++++
 .../simple_mlp/simple_mlp_train/script.py     |  7 ++++-
 .../ss_opm/ss_opm_train/config.vsh.yaml       |  1 +
 src/methods/ss_opm/ss_opm_train/script.py     |  4 +--
 src/utils/seed.py                             | 26 +++++++++++++++++++
 31 files changed, 132 insertions(+), 14 deletions(-)
 create mode 100644 src/utils/seed.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 62809422..953f115d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # task_predict_modality 0.2.1
 
+## MAJOR CHANGES
+
+* Seed every stochastic method and control, so that rerunning the benchmark reproduces its predictions. Unseeded in 0.2.0 were the truncated SVDs of `guanlab_dengkw_pm` (its `np.random.seed()` came after them) and `knnr_py`, the `lmds` landmarks of `lm` and `knnr_r`, the draws of `random_predict`, the scvi-tools model and pynndescent search of `cellmapper_scvi` and `cellmapper_linear`, and the training of `simple_mlp`, `novel`, `babel` and the `senkin_tmp` networks. The new `src/utils/seed.py` seeds python, numpy, torch and tensorflow and makes them use deterministic GPU algorithms, which `ss_opm`, `scipenn` and `scbutterfly`, already seeded, now use too; the R components call `set.seed()`. Methods without one get a `--seed` argument (default 0).
+
 ## BUG FIXES
 
 * `correlation`: Convert the solution and the prediction to CSC before correlating them. `dynutils::spearman_rank_sparse()` ranks along the compressed dimension, which is a column only for CSC, so a CSR matrix (read as a `dgRMatrix`) was ranked per feature for the per-cell Spearman correlation and per cell for the per-gene one. In 0.2.0 this hit `mean_spearman_per_cell` and `mean_spearman_per_gene` of `ss_opm`, which writes CSR, on every dataset, and of every method on GEX -> ATAC, whose solution holds the ATAC layer as CSR (`random_predict` also writes CSR there). All other scores are unchanged. A new component test scores a prediction as CSR and as CSC and checks both against dense ranks.
diff --git a/src/control_methods/random_predict/config.vsh.yaml b/src/control_methods/random_predict/config.vsh.yaml
index 04d36dcd..135b895d 100644
--- a/src/control_methods/random_predict/config.vsh.yaml
+++ b/src/control_methods/random_predict/config.vsh.yaml
@@ -3,6 +3,11 @@ name: random_predict
 label: Random predictions
 summary: Returns random training profiles.
 description: Returns random training profiles.
+arguments:
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for drawing the random training cells.
 resources:
   - type: r_script
     path: script.R
diff --git a/src/control_methods/random_predict/script.R b/src/control_methods/random_predict/script.R
index 94b9d4d4..0aac9b18 100644
--- a/src/control_methods/random_predict/script.R
+++ b/src/control_methods/random_predict/script.R
@@ -7,11 +7,14 @@ par <- list(
   input_train_mod1 = "resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/train_mod1.h5ad",
   input_test_mod1 = "resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/test_mod1.h5ad",
   input_train_mod2 = "resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/train_mod2.h5ad",
-  output = "output.h5ad"
+  output = "output.h5ad",
+  seed = 0L
 )
 meta <- list(name = "foo")
 ## VIASH END
 
+set.seed(par$seed)
+
 cat("Reading h5ad files\n")
 input_train_mod2 <- anndata::read_h5ad(par$input_train_mod2)
 input_test_mod1 <- anndata::read_h5ad(par$input_test_mod1)
diff --git a/src/methods/babel/babel_train/config.vsh.yaml b/src/methods/babel/babel_train/config.vsh.yaml
index b5c2c531..9c886ad2 100644
--- a/src/methods/babel/babel_train/config.vsh.yaml
+++ b/src/methods/babel/babel_train/config.vsh.yaml
@@ -17,6 +17,7 @@ resources:
   - path: ../chrom_utils.py
     dest: chrom_utils.py
   - path: /src/utils/exit_codes.py
+  - path: /src/utils/seed.py
 test_resources:
   - path: /resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
     dest: resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
@@ -47,6 +48,10 @@ arguments:
     type: integer
     default: 20
     description: Epochs without validation loss improvement before stopping early.
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the weight initialisation, the validation split and the order of the training batches.
 engines:
   - type: docker
     image: openproblems/base_pytorch_nvidia:1
diff --git a/src/methods/babel/babel_train/script.py b/src/methods/babel/babel_train/script.py
index eca8bac1..5059baab 100644
--- a/src/methods/babel/babel_train/script.py
+++ b/src/methods/babel/babel_train/script.py
@@ -26,6 +26,7 @@
     "batch_size": 64,
     "loss2_weight": 3.0,
     "early_stopping_patience": 20,
+    "seed": 0,
 }
 meta = {"name": "babel", "resources_dir": "src/methods/babel"}
 ## VIASH END
@@ -36,6 +37,7 @@
 from exit_codes import exit_non_applicable
 from model import AssymSplicedAutoEncoder
 from losses import QuadLoss
+from seed import seed_everything
 
 
 def _as_csr(X):
@@ -131,6 +133,8 @@ def get_loss(self, y_pred, y_true, X=None, training=False):
         return self.criterion_(preds11, preds12, preds21, preds22, target1, target2_bin)
 
 
+seed_everything(par["seed"])
+
 logger.info("Reading input files...")
 adata_mod1_train = ad.read_h5ad(par["input_train_mod1"])
 adata_mod2_train = ad.read_h5ad(par["input_train_mod2"])
diff --git a/src/methods/cellmapper_linear/config.vsh.yaml b/src/methods/cellmapper_linear/config.vsh.yaml
index 33e03ddf..c432a65b 100644
--- a/src/methods/cellmapper_linear/config.vsh.yaml
+++ b/src/methods/cellmapper_linear/config.vsh.yaml
@@ -57,9 +57,14 @@ arguments:
     type: "integer"
     default: 50
     description: Number of neighbors to consider for k-NN graph construction.
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for the approximate nearest neighbor search (pynndescent) and the fallback representation.
 resources:
   - type: python_script
     path: script.py
+  - path: /src/utils/seed.py
 engines:
   - type: docker
     image: openproblems/base_python:1
diff --git a/src/methods/cellmapper_linear/script.py b/src/methods/cellmapper_linear/script.py
index 7059f7fa..3a8518e8 100644
--- a/src/methods/cellmapper_linear/script.py
+++ b/src/methods/cellmapper_linear/script.py
@@ -1,3 +1,5 @@
+import sys
+
 import anndata as ad
 import cellmapper as cm
 from scipy.sparse import csc_matrix
@@ -13,13 +15,20 @@
     'fallback_representation': 'joint_pca',  # or None for fallback_representation
     'n_neighbors': 30,
     'kernel_method': 'gauss',
-    'mask_var': "hvg"  # variable to mask for fallback representation
+    'mask_var': "hvg",  # variable to mask for fallback representation
+    'seed': 0,
 }
 meta = {
   'name': 'cellmapper_linear',
+  'resources_dir': 'src/utils',
 }
 ## VIASH END
 
+sys.path.append(meta['resources_dir'])
+from seed import seed_everything
+
+seed_everything(par['seed'])
+
 print('Reading input files', flush=True)
 input_train_mod1 = ad.read_h5ad(par['input_train_mod1'])
 input_train_mod2 = ad.read_h5ad(par['input_train_mod2'])
diff --git a/src/methods/cellmapper_scvi/config.vsh.yaml b/src/methods/cellmapper_scvi/config.vsh.yaml
index aff93911..1b896e5f 100644
--- a/src/methods/cellmapper_scvi/config.vsh.yaml
+++ b/src/methods/cellmapper_scvi/config.vsh.yaml
@@ -62,11 +62,16 @@ arguments:
     type: boolean
     default: false
     description: Whether to plot the UMAP embedding of the latent space (for diagnoscic purposes)
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for the scvi-tools model and the approximate nearest neighbor search (pynndescent).
 resources:
   - type: python_script
     path: script.py
   - path: utils.py
     dest: utils.py 
+  - path: /src/utils/seed.py
 engines:
   - type: docker
     image: openproblems/base_pytorch_nvidia:1
diff --git a/src/methods/cellmapper_scvi/script.py b/src/methods/cellmapper_scvi/script.py
index ef873109..7493a46c 100644
--- a/src/methods/cellmapper_scvi/script.py
+++ b/src/methods/cellmapper_scvi/script.py
@@ -16,7 +16,7 @@
     'use_hvg': False,
     'adt_normalization': 'clr',  # Normalization method for ADT data
     'plot_umap': True,
-
+    'seed': 0,
 }
 meta = {
   'name': 'cellmapper_scvi',
@@ -25,8 +25,13 @@
 ## VIASH END
 
 sys.path.append(meta['resources_dir'])
+import scvi
+from seed import seed_everything
 from utils import get_representation
 
+seed_everything(par['seed'])
+scvi.settings.seed = par['seed']
+
 print('Reading input files', flush=True)
 input_train_mod1 = ad.read_h5ad(par['input_train_mod1'])
 input_train_mod2 = ad.read_h5ad(par['input_train_mod2'])
diff --git a/src/methods/guanlab_dengkw_pm/config.vsh.yaml b/src/methods/guanlab_dengkw_pm/config.vsh.yaml
index 4bad9a8d..ea6e18f4 100644
--- a/src/methods/guanlab_dengkw_pm/config.vsh.yaml
+++ b/src/methods/guanlab_dengkw_pm/config.vsh.yaml
@@ -29,7 +29,7 @@ arguments:
   - name: "--seed"
     type: "integer"
     default: 1000
-    description: Seed for the batch reshuffling.
+    description: Seed for the truncated SVDs and the batch reshuffling.
 resources:
   - type: python_script
     path: script.py
diff --git a/src/methods/guanlab_dengkw_pm/script.py b/src/methods/guanlab_dengkw_pm/script.py
index 98649d97..2417be5f 100644
--- a/src/methods/guanlab_dengkw_pm/script.py
+++ b/src/methods/guanlab_dengkw_pm/script.py
@@ -87,7 +87,7 @@ def kernel_ridge_predict(kernel, alpha, X, y, X_test):
 print('Models using the Truncated SVD to reduce the dimension', flush=True)
 
 if n_mod1 is not None and n_mod1 < input_train.n_vars:
-    embedder_mod1 = TruncatedSVD(n_components=n_mod1)
+    embedder_mod1 = TruncatedSVD(n_components=n_mod1, random_state=par['seed'])
     mod1_pca = embedder_mod1.fit_transform(input_train.layers["normalized"]).astype(np.float32)
     train_matrix = mod1_pca[input_train.obs['group'] == 'train']
     test_matrix = mod1_pca[input_train.obs['group'] == 'test']
@@ -96,7 +96,7 @@ def kernel_ridge_predict(kernel, alpha, X, y, X_test):
     test_matrix = input_test_mod1.to_df(layer="normalized").values.astype(np.float32)
   
 if n_mod2 is not None and n_mod2 < input_train_mod2.n_vars:
-    embedder_mod2 = TruncatedSVD(n_components=n_mod2)
+    embedder_mod2 = TruncatedSVD(n_components=n_mod2, random_state=par['seed'])
     train_gs = embedder_mod2.fit_transform(input_train_mod2.layers["normalized"]).astype(np.float32)
 else:
     # no dimensionality reduction, so the KRR predicts the mod2 features directly
diff --git a/src/methods/knnr_py/config.vsh.yaml b/src/methods/knnr_py/config.vsh.yaml
index db59bc7c..942161a4 100644
--- a/src/methods/knnr_py/config.vsh.yaml
+++ b/src/methods/knnr_py/config.vsh.yaml
@@ -23,6 +23,10 @@ arguments:
     type: "integer"
     default: 100
     description: Number of neighbors to use.
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for the truncated SVD.
 resources:
   - type: python_script
     path: script.py
diff --git a/src/methods/knnr_py/script.py b/src/methods/knnr_py/script.py
index b5029132..e9d6ad1c 100644
--- a/src/methods/knnr_py/script.py
+++ b/src/methods/knnr_py/script.py
@@ -12,6 +12,7 @@
     'output': 'output.h5ad',
     'n_pcs': 4,
     'n_neighbors': 5,
+    'seed': 0,
 }
 meta = { 'name': 'foo' }
 ## VIASH END
@@ -31,7 +32,7 @@
 )
 
 print('Performing dimensionality reduction on modality 1 values...', flush=True)
-embedder = TruncatedSVD(n_components=par['n_pcs'])
+embedder = TruncatedSVD(n_components=par['n_pcs'], random_state=par['seed'])
 X = embedder.fit_transform(input_train.layers["normalized"])
 
 # split dimred back up
diff --git a/src/methods/knnr_r/config.vsh.yaml b/src/methods/knnr_r/config.vsh.yaml
index d8d3d22e..3427cff2 100644
--- a/src/methods/knnr_r/config.vsh.yaml
+++ b/src/methods/knnr_r/config.vsh.yaml
@@ -23,6 +23,10 @@ arguments:
     type: "integer"
     default: 20
     description: Number of neighbors to use in the knn regression.
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for the landmarks of the landmark MDS.
 resources:
   - type: r_script
     path: script.R
diff --git a/src/methods/knnr_r/script.R b/src/methods/knnr_r/script.R
index 0c9589b3..4813b492 100644
--- a/src/methods/knnr_r/script.R
+++ b/src/methods/knnr_r/script.R
@@ -11,11 +11,14 @@ par <- list(
   output = "output.h5ad",
   n_pcs = 4L,
   n_neighbors = 3,
-  distance_method = "pearson"
+  distance_method = "pearson",
+  seed = 0L
 )
 meta <- list(name = "foo")
 ## VIASH END
 
+set.seed(par$seed)
+
 cat("Reading mod1 h5ad files\n")
 input_train_mod1 <- anndata::read_h5ad(par$input_train_mod1)
 dataset_id <- input_train_mod1$uns[["dataset_id"]]
diff --git a/src/methods/lm/config.vsh.yaml b/src/methods/lm/config.vsh.yaml
index c452ed1e..fbc783d4 100644
--- a/src/methods/lm/config.vsh.yaml
+++ b/src/methods/lm/config.vsh.yaml
@@ -19,6 +19,10 @@ arguments:
     type: "integer"
     default: 50
     description: Number of principal components to use.
+  - name: "--seed"
+    type: "integer"
+    default: 0
+    description: Seed for the landmarks of the landmark MDS.
 resources:
   - type: r_script
     path: script.R
diff --git a/src/methods/lm/script.R b/src/methods/lm/script.R
index e46efc91..90245a1e 100644
--- a/src/methods/lm/script.R
+++ b/src/methods/lm/script.R
@@ -9,11 +9,14 @@ par <- list(
   input_test_mod1 = paste0(path, "test_mod1.h5ad"),
   input_train_mod2 = paste0(path, "train_mod2.h5ad"),
   output = "output.h5ad",
-  n_pcs = 4L
+  n_pcs = 4L,
+  seed = 0L
 )
 meta <- list(name = "foo")
 ## VIASH END
 
+set.seed(par$seed)
+
 cat("Reading mod1 files\n")
 input_train_mod1 <- anndata::read_h5ad(par$input_train_mod1)
 input_test_mod1 <- anndata::read_h5ad(par$input_test_mod1)
diff --git a/src/methods/novel/novel_train/config.vsh.yaml b/src/methods/novel/novel_train/config.vsh.yaml
index bd584d6d..0359d6da 100644
--- a/src/methods/novel/novel_train/config.vsh.yaml
+++ b/src/methods/novel/novel_train/config.vsh.yaml
@@ -15,6 +15,7 @@ resources:
   - path: script.py
     type: python_script
   - path: ../helper_functions.py
+  - path: /src/utils/seed.py
 engines:
   - type: docker
     image: openproblems/base_pytorch_nvidia:1
diff --git a/src/methods/novel/novel_train/script.py b/src/methods/novel/novel_train/script.py
index 91e35885..b0621123 100644
--- a/src/methods/novel/novel_train/script.py
+++ b/src/methods/novel/novel_train/script.py
@@ -35,6 +35,9 @@
 sys.path.append(meta['resources_dir'])
 from helper_functions import train_and_valid, lsiTransformer, ModalityMatchingDataset
 from helper_functions import ModelRegressionAtac2Gex, ModelRegressionAdt2Gex, ModelRegressionGex2Adt, ModelRegressionGex2Atac
+from seed import seed_everything
+
+seed_everything(par['seed'])
 
 print('Load data', flush=True)
 input_train_mod1 = ad.read_h5ad(par['input_train_mod1'])
diff --git a/src/methods/scbutterfly/butterfly_common.py b/src/methods/scbutterfly/butterfly_common.py
index 1f0b6e62..d97856e1 100644
--- a/src/methods/scbutterfly/butterfly_common.py
+++ b/src/methods/scbutterfly/butterfly_common.py
@@ -20,9 +20,13 @@
 from scipy.sparse import csr_matrix, issparse
 
 import chrom_utils
+from seed import seed_everything
 
 logger = logging.getLogger(__name__)
 
+# scButterfly seeds python, numpy and torch with this itself (Butterfly(), train_model())
+SEED = 19193
+
 
 def apply_runtime_patches():
     """Patch scButterfly/torch quirks. Call before importing scButterfly.
@@ -386,6 +390,8 @@ def build_butterfly(train_mod1, train_mod2, test_mod1, n_top_genes, Butterfly,
     sort_index, chrom_list = chrom_utils.sorted_chrom_order(ATAC_data)
     ATAC_data = chrom_utils.apply_sort(ATAC_data, sort_index)
 
+    # adds deterministic GPU algorithms to the seeding scButterfly does
+    seed_everything(SEED)
     butterfly = Butterfly()
     butterfly.load_data(RNA_data, ATAC_data, train_id_final, test_id, validation_id)
     del RNA_data, ATAC_data  # load_data keeps its own copies
diff --git a/src/methods/scbutterfly/scbutterfly_predict/config.vsh.yaml b/src/methods/scbutterfly/scbutterfly_predict/config.vsh.yaml
index c4289f5f..1b0f6132 100644
--- a/src/methods/scbutterfly/scbutterfly_predict/config.vsh.yaml
+++ b/src/methods/scbutterfly/scbutterfly_predict/config.vsh.yaml
@@ -14,6 +14,7 @@ resources:
   - path: ../chrom_utils.py
   - path: /src/utils/exit_codes.py
   - path: /src/utils/cell_scale.py
+  - path: /src/utils/seed.py
 test_resources:
   - path: /resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
     dest: resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
diff --git a/src/methods/scbutterfly/scbutterfly_train/config.vsh.yaml b/src/methods/scbutterfly/scbutterfly_train/config.vsh.yaml
index f3542921..943e3343 100644
--- a/src/methods/scbutterfly/scbutterfly_train/config.vsh.yaml
+++ b/src/methods/scbutterfly/scbutterfly_train/config.vsh.yaml
@@ -18,6 +18,7 @@ resources:
   - path: ../chrom_utils.py
   - path: /src/utils/exit_codes.py
   - path: /src/utils/cell_scale.py
+  - path: /src/utils/seed.py
 test_resources:
   - path: /resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
     dest: resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap
diff --git a/src/methods/scipenn/config.vsh.yaml b/src/methods/scipenn/config.vsh.yaml
index 8b4d170e..07e8a840 100644
--- a/src/methods/scipenn/config.vsh.yaml
+++ b/src/methods/scipenn/config.vsh.yaml
@@ -45,6 +45,7 @@ resources:
   - type: python_script
     path: script.py
   - path: /src/utils/exit_codes.py
+  - path: /src/utils/seed.py
 engines:
   # sciPENN is installed --no-deps: it pins a spurious numba<=0.50.0 (never imported) and a
   # 2021-era stack, but runs on the base image's torch/scanpy/anndata. The base already ships
diff --git a/src/methods/scipenn/script.py b/src/methods/scipenn/script.py
index 566acffd..ad21832c 100644
--- a/src/methods/scipenn/script.py
+++ b/src/methods/scipenn/script.py
@@ -29,6 +29,7 @@
 
 sys.path.append(meta["resources_dir"])
 from exit_codes import exit_non_applicable
+from seed import seed_everything
 
 logger.info("Reading input files...")
 rna_train = ad.read_h5ad(par["input_train_mod1"])   # mod1 = RNA (GEX)
@@ -46,8 +47,7 @@
         f"got mod1={mod1!r}, mod2={mod2!r}."
     )
 
-np.random.seed(par["seed"])
-torch.manual_seed(par["seed"])
+seed_everything(par["seed"])
 
 # Feed sciPENN the task's log_cp10k layer (its internal normalize_total/log1p are
 # disabled below) so predictions live in the same space as the ground truth. sciPENN
diff --git a/src/methods/senkin_tmp/senkin_tmp_train/config.vsh.yaml b/src/methods/senkin_tmp/senkin_tmp_train/config.vsh.yaml
index 87f24bbd..df0bd346 100644
--- a/src/methods/senkin_tmp/senkin_tmp_train/config.vsh.yaml
+++ b/src/methods/senkin_tmp/senkin_tmp_train/config.vsh.yaml
@@ -14,6 +14,7 @@ resources:
     type: python_script
   - path: /src/utils/exit_codes.py
   - path: /src/utils/cell_scale.py
+  - path: /src/utils/seed.py
 arguments:
   - name: "--n_folds"
     type: integer
diff --git a/src/methods/senkin_tmp/senkin_tmp_train/script.py b/src/methods/senkin_tmp/senkin_tmp_train/script.py
index 8404b4ab..70949f85 100644
--- a/src/methods/senkin_tmp/senkin_tmp_train/script.py
+++ b/src/methods/senkin_tmp/senkin_tmp_train/script.py
@@ -45,12 +45,15 @@
 sys.path.append(meta["resources_dir"])
 from cell_scale import CellScale
 from exit_codes import exit_non_applicable
+from seed import seed_everything
 
 # The original solution corrected batch effects per day and computed gene-protein correlations per donor and day.
 # The benchmark datasets only carry a generic `batch` column (e.g. "s1d1" in NeurIPS 2021, "day_donor" in
 # NeurIPS 2022), which is used for both purposes here. Nothing dataset specific is assumed about its format.
 BATCH_KEY = "batch"
 SEED = 42
+# the neural networks (tensorflow, imported by senkin_tmp_cite_pred) shuffle with numpy and were not seeded
+seed_everything(SEED)
 
 # ---------------------------------------------------------------------------
 # Load data
diff --git a/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml b/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
index e10fd7bd..4376a020 100644
--- a/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
+++ b/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
@@ -7,10 +7,15 @@ arguments:
     description: Number of training epochs. Defaults to the value in the bundled model config.
     info:
       test_default: 2
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the weight initialisation and the order of the training batches.
 resources:
   - type: python_script
     path: script.py
   - path: ../resources/
+  - path: /src/utils/seed.py
 engines:
   - type: docker
     image: openproblems/base_pytorch_nvidia:1
diff --git a/src/methods/simple_mlp/simple_mlp_train/script.py b/src/methods/simple_mlp/simple_mlp_train/script.py
index 7f35e36e..c2fb49cf 100644
--- a/src/methods/simple_mlp/simple_mlp_train/script.py
+++ b/src/methods/simple_mlp/simple_mlp_train/script.py
@@ -19,7 +19,8 @@
     'input_train_mod2': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap/train_mod2.h5ad',
     'input_test_mod1': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/swap/test_mod1.h5ad',
     'output': 'output/model',
-    'n_epochs': None
+    'n_epochs': None,
+    'seed': 0
 }
 meta = {
     'resources_dir': 'src/methods/simple_mlp',
@@ -31,7 +32,9 @@
 
 import sys
 sys.path.append(resources_dir)
+sys.path.append(meta['resources_dir'])
 from models import MLP
+from seed import seed_everything
 import utils
 
 def _train(X, y, Xt, yt, logger, config, num_workers):
@@ -83,6 +86,8 @@ def _train(X, y, Xt, yt, logger, config, num_workers):
 
 
 
+seed_everything(par['seed'])
+
 input_train_mod1 = ad.read_h5ad(par['input_train_mod1'])
 input_train_mod2 = ad.read_h5ad(par['input_train_mod2'])
 # only the normalized layers are used, and the counts take as much memory again
diff --git a/src/methods/ss_opm/ss_opm_train/config.vsh.yaml b/src/methods/ss_opm/ss_opm_train/config.vsh.yaml
index 2efde543..c3ca44c4 100644
--- a/src/methods/ss_opm/ss_opm_train/config.vsh.yaml
+++ b/src/methods/ss_opm/ss_opm_train/config.vsh.yaml
@@ -53,6 +53,7 @@ resources:
     path: script.py
   - path: ../ss_opm_common.py
   - path: /src/utils/cell_scale.py
+  - path: /src/utils/seed.py
 engines:
   - type: docker
     image: openproblems/base_pytorch_nvidia:1
diff --git a/src/methods/ss_opm/ss_opm_train/script.py b/src/methods/ss_opm/ss_opm_train/script.py
index bb8431aa..377dd4d7 100644
--- a/src/methods/ss_opm/ss_opm_train/script.py
+++ b/src/methods/ss_opm/ss_opm_train/script.py
@@ -28,6 +28,7 @@
 
 sys.path.append(meta["resources_dir"])
 from cell_scale import CellScale  # noqa: E402
+from seed import seed_everything  # noqa: E402
 from ss_opm_common import (  # noqa: E402
     ROW_BLOCK,
     BatchSingularVectors,
@@ -55,7 +56,6 @@
 from ss_opm.model.encoder_decoder.encoder_decoder import EncoderDecoder  # noqa: E402
 from ss_opm.pre_post_processing.pre_post_processing import PrePostProcessing  # noqa: E402
 from ss_opm.utility.row_normalize import row_normalize  # noqa: E402
-from ss_opm.utility.set_seed import set_seed  # noqa: E402
 
 device = "cuda" if torch.cuda.is_available() else "cpu"
 print(f"Using device: {device}", flush=True)
@@ -80,7 +80,7 @@ def _float32_train_step_forward(self, batch, training_length_ratio):
     os.environ.setdefault("OMP_NUM_THREADS", str(max(1, min(int(meta["cpus"]), 8))))
 
 SEED = 42
-set_seed(SEED)
+seed_everything(SEED)  # the original's set_seed(), plus deterministic GPU algorithms
 
 # ---- Load data ----
 print("Loading data...", flush=True)
diff --git a/src/utils/seed.py b/src/utils/seed.py
new file mode 100644
index 00000000..10e79cd5
--- /dev/null
+++ b/src/utils/seed.py
@@ -0,0 +1,26 @@
+import os
+import random
+import sys
+
+import numpy as np
+
+
+def seed_everything(seed):
+    """Seed python's and numpy's random number generators, and those of torch and tensorflow when the method has
+    imported them (so call this after the imports), and make them pick deterministic GPU algorithms, so that a method
+    gives the same predictions every time it runs on the same data."""
+    random.seed(seed)
+    np.random.seed(seed)
+    # cuBLAS is only deterministic with a fixed workspace, read when CUDA initialises
+    os.environ.setdefault("CUBLAS_WORKSPACE_CONFIG", ":4096:8")
+    if "torch" in sys.modules:
+        torch = sys.modules["torch"]
+        torch.manual_seed(seed)  # seeds every CUDA device too
+        torch.backends.cudnn.deterministic = True
+        torch.backends.cudnn.benchmark = False
+        # warn rather than fail on the few operations that have no deterministic implementation
+        torch.use_deterministic_algorithms(True, warn_only=True)
+    if "tensorflow" in sys.modules:
+        tensorflow = sys.modules["tensorflow"]
+        tensorflow.keras.utils.set_random_seed(seed)
+        tensorflow.config.experimental.enable_op_determinism()