Fix/correlation csr guanlab pbmc multiome

CHANGELOG.md

@@ -1,3 +1,17 @@
+# 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.
+
+* `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.

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

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)

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

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"])

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

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'])

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

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'])

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
@@ -48,4 +48,4 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [hightime, highmem, highcpu]
+      label: [veryhightime, veryhighmem, highcpu]

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
 
 
@@ -92,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']
@@ -101,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

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

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

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

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"]]

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

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)

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

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'])

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

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

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

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