diff --git a/CHANGELOG.md b/CHANGELOG.md
index 583e59ca..4943d5de 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,17 @@
+# task_predict_modality 0.2.1
+
+## MAJOR CHANGES
+
+* Add `--seed` to all methods, controls and `run_benchmark`, so that reruns reproduce the predictions (PR #83).
+
+## BUG FIXES
+
+* `correlation`: Convert matrices to CSC, since dynutils < 1.0.13 gave wrong Spearman scores on CSR (PR #83).
+
+* `guanlab_dengkw_pm`: Replace `scipy.linalg.cho_factor()`, which failed on `pbmc_multiome`, with `numpy.linalg.solve()`, and ask for more memory and time (PR #83).
+
+* `simple_mlp`: Densify one batch at a time, which fixes running out of memory on `pbmc_multiome` (PR #83).
+
 # 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/README.md b/README.md
index c07e9896..49049fde 100644
--- a/README.md
+++ b/README.md
@@ -415,6 +415,7 @@ Arguments:
 | `--input_train_mod2` | `file` | The mod2 expression values of the train cells. |
 | `--input_test_mod1` | `file` | The mod1 expression values of the test cells. |
 | `--output` | `file` | (*Output*) A prediction of the mod2 expression values of the test cells. |
+| `--seed` | `integer` | (*Optional*) Seed for the random number generators of the method. Default: `0`. |
 
 </div>
 
@@ -432,6 +433,7 @@ Arguments:
 | `--input_train_mod2` | `file` | The mod2 expression values of the train cells. |
 | `--input_test_mod1` | `file` | (*Optional*) The mod1 expression values of the test cells. |
 | `--output` | `file` | (*Output*) A pretrained model for predicting the expression of one modality from another. |
+| `--seed` | `integer` | (*Optional*) Seed for the random number generators used for training the model. Default: `0`. |
 
 </div>
 
@@ -450,6 +452,7 @@ Arguments:
 | `--input_test_mod1` | `file` | The mod1 expression values of the test cells. |
 | `--input_test_mod2` | `file` | The ground-truth mod2 expression values of the test cells. |
 | `--output` | `file` | (*Output*) A prediction of the mod2 expression values of the test cells. |
+| `--seed` | `integer` | (*Optional*) Seed for the random number generators of the control method, if it uses any. Default: `0`. |
 
 </div>
 
diff --git a/src/api/comp_control_method.yaml b/src/api/comp_control_method.yaml
index 7fc6f201..11921d4c 100644
--- a/src/api/comp_control_method.yaml
+++ b/src/api/comp_control_method.yaml
@@ -32,6 +32,10 @@ arguments:
     __merge__: file_prediction.yaml
     direction: output
     required: true
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the random number generators of the control method, if it uses any.
 test_resources:
   - type: python_script
     path: /common/component_tests/check_config.py
diff --git a/src/api/comp_method.yaml b/src/api/comp_method.yaml
index 9414b6ea..8451064f 100644
--- a/src/api/comp_method.yaml
+++ b/src/api/comp_method.yaml
@@ -23,6 +23,10 @@ arguments:
     __merge__: file_prediction.yaml
     direction: output
     required: true
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the random number generators of the method.
 test_resources:
   - type: python_script
     path: /common/component_tests/check_config.py
diff --git a/src/api/comp_method_train.yaml b/src/api/comp_method_train.yaml
index 3e3a1c8f..e8565062 100644
--- a/src/api/comp_method_train.yaml
+++ b/src/api/comp_method_train.yaml
@@ -23,6 +23,10 @@ arguments:
     __merge__: file_pretrained_model.yaml
     direction: output
     required: true
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the random number generators used for training the model.
 test_resources:
   - type: python_script
     path: /common/component_tests/run_and_check_output.py
diff --git a/src/api/wf_method.yaml b/src/api/wf_method.yaml
index 38fd735a..47476a83 100644
--- a/src/api/wf_method.yaml
+++ b/src/api/wf_method.yaml
@@ -25,6 +25,10 @@ arguments:
     __merge__: file_prediction.yaml
     direction: output
     required: true
+  - name: "--seed"
+    type: integer
+    default: 0
+    description: Seed for the random number generators of the method.
 test_resources:
   - type: python_script
     path: /common/component_tests/check_config.py
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/main.nf b/src/methods/babel/babel/main.nf
index 9ed90555..f5fe0bd6 100644
--- a/src/methods/babel/babel/main.nf
+++ b/src/methods/babel/babel/main.nf
@@ -3,7 +3,7 @@ workflow run_wf {
   main:
   output_ch = input_ch
     | babel_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1"],
+      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1", "seed"],
       toState: ["input_model": "output"]
     )
     | babel_predict.run(
diff --git a/src/methods/babel/babel_train/config.vsh.yaml b/src/methods/babel/babel_train/config.vsh.yaml
index b5c2c531..9c0617a0 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
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..b7368b06 100644
--- a/src/methods/cellmapper_linear/config.vsh.yaml
+++ b/src/methods/cellmapper_linear/config.vsh.yaml
@@ -60,6 +60,7 @@ arguments:
 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..f5befaac 100644
--- a/src/methods/cellmapper_scvi/config.vsh.yaml
+++ b/src/methods/cellmapper_scvi/config.vsh.yaml
@@ -67,6 +67,7 @@ resources:
     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 4f1fc942..d69daec4 100644
--- a/src/methods/guanlab_dengkw_pm/config.vsh.yaml
+++ b/src/methods/guanlab_dengkw_pm/config.vsh.yaml
@@ -26,10 +26,6 @@ arguments:
       kernel ridge model, so the consensus averages `2 * n_repeats` predictions.
     info:
       test_default: 1
-  - name: "--seed"
-    type: "integer"
-    default: 1000
-    description: Seed for the batch reshuffling.
 resources:
   - type: python_script
     path: script.py
@@ -48,4 +44,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..4def211a 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
@@ -14,7 +13,7 @@
     'input_test_mod1': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/test_mod1.h5ad',
     'output': 'output.h5ad',
     'n_repeats': 5,
-    'seed': 1000
+    'seed': 0
 }
 meta = {
     'name': 'guanlab_dengkw_pm',
@@ -32,22 +31,23 @@
 # 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"), which fails on large systems.
+# Tried on the halves of the batches of pbmc_multiome (56k-74k cells), or at the same size:
+# - scipy.linalg.solve(assume_a="pos"): raises a MemoryError or segfaults, already from ~30-40k cells
+# - scipy.linalg.cho_factor(): reports K + alpha * I as not positive definite ("potrf return info")
+# - numpy.linalg.cholesky(): segfaults
+# - numpy.linalg.solve(): works, 270 s for 66k cells on 30 threads. 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.
+# The Cholesky failures depend on the CPU: on the Cascade Lake workers of de.NBI, OpenBLAS picks its AVX-512 kernel
+# (SkylakeX). cho_factor() succeeds at 66k cells with the Haswell kernel, on AMD CPUs or with
+# OPENBLAS_CORETYPE=Haswell, but numpy.linalg.solve() doesn't need that set on every backend.
 #
-# Switch back to KernelRidge once scipy.linalg.solve() handles matrices of this size again.
+# Switch back to KernelRidge once OpenBLAS's Cholesky handles these systems on every CPU.
 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
 
 
@@ -92,7 +92,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']
@@ -101,7 +101,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/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/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/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/main.nf b/src/methods/novel/novel/main.nf
index 0d973afb..209c45fa 100644
--- a/src/methods/novel/novel/main.nf
+++ b/src/methods/novel/novel/main.nf
@@ -3,7 +3,7 @@ workflow run_wf {
   main:
   output_ch = input_ch
     | novel_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2"],
+      fromState: ["input_train_mod1", "input_train_mod2", "seed"],
       toState: ["input_model": "output"]
     )
     | novel_predict.run(
diff --git a/src/methods/novel/novel_train/config.vsh.yaml b/src/methods/novel/novel_train/config.vsh.yaml
index bd584d6d..302b2893 100644
--- a/src/methods/novel/novel_train/config.vsh.yaml
+++ b/src/methods/novel/novel_train/config.vsh.yaml
@@ -7,14 +7,11 @@ arguments:
     description: Number of training epochs.
     info:
       test_default: 2
-  - name: "--seed"
-    type: "integer"
-    default: 1
-    description: "Seed for choosing the internal validation batches, when the dataset does not contain the batches used in the competition."
 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..2cc3f53c 100644
--- a/src/methods/novel/novel_train/script.py
+++ b/src/methods/novel/novel_train/script.py
@@ -25,7 +25,7 @@
   'input_train_mod2': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/train_mod2.h5ad',
   'output': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/models/novel',
   'n_epochs': 100,
-  'seed': 1
+  'seed': 0
 }
 meta = {
   'resources_dir': 'src/methods/novel',
@@ -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..a4084827 100644
--- a/src/methods/scbutterfly/butterfly_common.py
+++ b/src/methods/scbutterfly/butterfly_common.py
@@ -20,6 +20,7 @@
 from scipy.sparse import csr_matrix, issparse
 
 import chrom_utils
+from seed import seed_everything
 
 logger = logging.getLogger(__name__)
 
@@ -169,7 +170,7 @@ class _SparseBatchLoader:
     Batches hold the same values in the same layout (RNA columns, then ATAC), and come
     back on the device as float32, so the training loop's ``.cuda().to(torch.float32)``
     is a no-op. Shuffling draws a ``torch.randperm``, like ``DataLoader``, so it follows
-    the seed scButterfly sets.
+    the seed.
     """
 
     def __init__(self, dataset, batch_size, shuffle=False, num_workers=0, drop_last=False):
@@ -322,7 +323,7 @@ def validation_cells(n_train):
     return sorted(shuffled[:max(1, int(0.1 * n_train))])
 
 
-def build_butterfly(train_mod1, train_mod2, test_mod1, n_top_genes, Butterfly,
+def build_butterfly(train_mod1, train_mod2, test_mod1, n_top_genes, Butterfly, seed,
                     model_all_target_features=True):
     """Build + preprocess + construct a Butterfly for the given data.
 
@@ -337,7 +338,7 @@ def build_butterfly(train_mod1, train_mod2, test_mod1, n_top_genes, Butterfly,
     Returns a dict with the constructed ``butterfly``, the resolved ``direction``,
     ``test_id`` and ``validation_id`` (row indices of the test and held-out training
     cells), ``chrom_list``, and the preprocessed and the original target var names.
-    Deterministic given the same inputs, so train and predict produce identical
+    Deterministic given the same inputs and seed, so train and predict produce identical
     architecture/preprocessing.
     """
     direction, mod1, mod2 = detect_direction(train_mod1, train_mod2)
@@ -387,6 +388,8 @@ def build_butterfly(train_mod1, train_mod2, test_mod1, n_top_genes, Butterfly,
     ATAC_data = chrom_utils.apply_sort(ATAC_data, sort_index)
 
     butterfly = Butterfly()
+    # Butterfly() seeds python, numpy and torch with 19193
+    seed_everything(seed)
     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
     butterfly.data_preprocessing(
diff --git a/src/methods/scbutterfly/scbutterfly/main.nf b/src/methods/scbutterfly/scbutterfly/main.nf
index ceb9e913..4ca0b5f7 100644
--- a/src/methods/scbutterfly/scbutterfly/main.nf
+++ b/src/methods/scbutterfly/scbutterfly/main.nf
@@ -3,7 +3,7 @@ workflow run_wf {
   main:
   output_ch = input_ch
     | scbutterfly_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1"],
+      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1", "seed"],
       toState: ["input_model": "output"]
     )
     | scbutterfly_predict.run(
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_predict/script.py b/src/methods/scbutterfly/scbutterfly_predict/script.py
index c2983dc0..1c25ee61 100644
--- a/src/methods/scbutterfly/scbutterfly_predict/script.py
+++ b/src/methods/scbutterfly/scbutterfly_predict/script.py
@@ -49,6 +49,8 @@
 built = butterfly_common.build_butterfly(
     train_mod1, train_mod2, test_mod1,
     n_top_genes=metadata["n_top_genes"], Butterfly=Butterfly,
+    # Models trained before --seed was added used scButterfly's own seed.
+    seed=metadata.get("seed", 19193),
     # Models trained before every target feature was modelled lack the key.
     model_all_target_features=metadata.get("model_all_target_features", False),
 )
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/scbutterfly/scbutterfly_train/script.py b/src/methods/scbutterfly/scbutterfly_train/script.py
index b9425829..98eec193 100644
--- a/src/methods/scbutterfly/scbutterfly_train/script.py
+++ b/src/methods/scbutterfly/scbutterfly_train/script.py
@@ -20,6 +20,7 @@
     "patience": 50,
     "batch_size": 64,
     "n_top_genes": 3000,
+    "seed": 0,
 }
 meta = {"name": "scbutterfly", "resources_dir": "src/methods/scbutterfly"}
 ## VIASH END
@@ -66,7 +67,7 @@
 logger.info("Building scButterfly model...")
 built = butterfly_common.build_butterfly(
     train_mod1, train_mod2, test_mod1,
-    n_top_genes=par["n_top_genes"], Butterfly=Butterfly,
+    n_top_genes=par["n_top_genes"], Butterfly=Butterfly, seed=par["seed"],
 )
 butterfly = built["butterfly"]
 
@@ -78,6 +79,7 @@
     patience=par["patience"],
     batch_size=par["batch_size"],
     output_path=par["output"],
+    seed=par["seed"],
 )
 
 # Spread of the per-cell scale, from the held-out cells predicted with the weights
@@ -100,6 +102,7 @@
     "n_top_genes": par["n_top_genes"],
     "model_all_target_features": True,
     "batch_size": par["batch_size"],
+    "seed": par["seed"],
     "target_var_names": list(train_mod2.var_names),
     "dataset_id": train_mod1.uns.get("dataset_id", ""),
 }
diff --git a/src/methods/scipenn/config.vsh.yaml b/src/methods/scipenn/config.vsh.yaml
index 8b4d170e..70212eac 100644
--- a/src/methods/scipenn/config.vsh.yaml
+++ b/src/methods/scipenn/config.vsh.yaml
@@ -37,14 +37,11 @@ arguments:
     type: integer
     description: "Training batch size."
     default: 128
-  - name: "--seed"
-    type: integer
-    description: "Seed for the train/validation split, shuffling and weight initialization."
-    default: 0
 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/main.nf b/src/methods/senkin_tmp/senkin_tmp/main.nf
index d527c816..36f0baff 100644
--- a/src/methods/senkin_tmp/senkin_tmp/main.nf
+++ b/src/methods/senkin_tmp/senkin_tmp/main.nf
@@ -3,7 +3,7 @@ workflow run_wf {
   main:
   output_ch = input_ch
     | senkin_tmp_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1"],
+      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1", "seed"],
       toState: ["input_model": "output"]
     )
     | senkin_tmp_predict.run(
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..967f4219 100644
--- a/src/methods/senkin_tmp/senkin_tmp_train/script.py
+++ b/src/methods/senkin_tmp/senkin_tmp_train/script.py
@@ -38,6 +38,7 @@
     "lgbm_learning_rate": 0.1,
     "lgbm_max_bin": 63,
     "lgbm_n_jobs": -1,
+    "seed": 0,
 }
 meta = {"name": "senkin_tmp", "resources_dir": "src/methods/senkin_tmp/senkin_tmp_train", "cpus": None, "memory_gb": None}
 ## VIASH END
@@ -45,12 +46,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
+SEED = par["seed"]
+# the neural networks (tensorflow, imported by senkin_tmp_cite_pred) shuffle with numpy
+seed_everything(SEED)
 
 # ---------------------------------------------------------------------------
 # Load data
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/main.nf b/src/methods/simple_mlp/simple_mlp/main.nf
index 997a1d33..dfba2315 100644
--- a/src/methods/simple_mlp/simple_mlp/main.nf
+++ b/src/methods/simple_mlp/simple_mlp/main.nf
@@ -4,7 +4,7 @@ workflow run_wf {
   output_ch = input_ch
 
     | simple_mlp_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2"],
+      fromState: ["input_train_mod1", "input_train_mod2", "seed"],
       toState: ["input_model": "output"]
     )
 
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/config.vsh.yaml b/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
index e10fd7bd..0635d00a 100644
--- a/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
+++ b/src/methods/simple_mlp/simple_mlp_train/config.vsh.yaml
@@ -11,6 +11,7 @@ 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 7f90f4e0..c2fb49cf 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
 
@@ -20,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',
@@ -32,29 +32,25 @@
 
 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):
 
-    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,15 +79,20 @@ 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()
 
 
 
+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
+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"]
diff --git a/src/methods/ss_opm/ss_opm/main.nf b/src/methods/ss_opm/ss_opm/main.nf
index 9f30a8b4..681f28a1 100644
--- a/src/methods/ss_opm/ss_opm/main.nf
+++ b/src/methods/ss_opm/ss_opm/main.nf
@@ -3,7 +3,7 @@ workflow run_wf {
   main:
   output_ch = input_ch
     | ss_opm_train.run(
-      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1"],
+      fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1", "seed"],
       toState: ["input_model": "output"]
     )
     | ss_opm_predict.run(
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..a57a70c1 100644
--- a/src/methods/ss_opm/ss_opm_train/script.py
+++ b/src/methods/ss_opm/ss_opm_train/script.py
@@ -22,12 +22,14 @@
     "n_epochs": 40,
     "burnin_length_epoch": 10,
     "n_rescaling_cells": 5000,
+    "seed": 0,
 }
 meta = {"name": "ss_opm_train", "resources_dir": "src/methods/ss_opm", "cpus": None}
 ## VIASH END
 
 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 +57,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)
@@ -79,8 +80,8 @@ def _float32_train_step_forward(self, batch, training_length_ratio):
 if meta.get("cpus"):
     os.environ.setdefault("OMP_NUM_THREADS", str(max(1, min(int(meta["cpus"]), 8))))
 
-SEED = 42
-set_seed(SEED)
+SEED = par["seed"]
+seed_everything(SEED)  # the original's set_seed(), plus deterministic GPU algorithms
 
 # ---- Load data ----
 print("Loading data...", flush=True)
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..e2d04725 100644
--- a/src/metrics/correlation/script.R
+++ b/src/metrics/correlation/script.R
@@ -30,15 +30,15 @@ expect_true(
 
 cat("Computing correlation metrics\n")
 # Wrangle data
-tv <- ad_sol$layers[["normalized"]]
-pv <- ad_pred$layers[["normalized"]]
+# dynutils < 1.0.13 ranks a CSR matrix (dgRMatrix) 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..e3940e90
--- /dev/null
+++ b/src/metrics/correlation/test_prediction_layout.py
@@ -0,0 +1,56 @@
+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 < 1.0.13 ranks a CSR prediction 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)
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()
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 9ee417f2..440e2bcd 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -47,6 +47,12 @@ argument_groups:
         required: true
         direction: output
         default: task_info.yaml
+  - name: Method arguments
+    arguments:
+      - name: "--seed"
+        type: integer
+        default: 0
+        description: Seed for all of the methods and control methods. A `seed` in a parameter set takes precedence.
   - name: Method filtering
     description: |
       Use these arguments to filter methods by name. By default, all methods are
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index e6ce777a..efd41698 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -126,7 +126,8 @@ workflow run_wf {
         def new_args = [
           input_train_mod1: state.input_train_mod1,
           input_train_mod2: state.input_train_mod2,
-          input_test_mod1: state.input_test_mod1
+          input_test_mod1: state.input_test_mod1,
+          seed: state.seed
         ]
         if (comp.config.info.type == "control_method") {
           new_args.input_test_mod2 = state.input_test_mod2