Backport PR #3853: fix: scanpy benchmark uses pip for igraph + setup_cache (#3859)

Author: ilan-gold

.github/workflows/benchmark.yml

@@ -48,7 +48,7 @@ jobs:
           key: benchmark-state-${{ hashFiles('benchmarks/**') }}
 
       - name: Install dependencies
-        run: pip install 'asv>=0.6.4'
+        run: pip install 'asv>=0.6.4' py-rattler
 
       - name: Configure ASV
         working-directory: ${{ env.ASV_DIR }}

benchmarks/asv.conf.json

@@ -43,7 +43,7 @@
     // If missing or the empty string, the tool will be automatically
     // determined by looking for tools on the PATH environment
     // variable.
-    "environment_type": "conda",
+    "environment_type": "rattler",
 
     // timeout in seconds for installing any dependencies in environment
     // defaults to 10 min
@@ -78,13 +78,12 @@
         "natsort": [""],
         "pandas": [""],
         "memory_profiler": [""],
-        "zarr": ["2.18.4"],
+        "zarr": [""],
         "pytest": [""],
-        "scanpy": [""],
-        "python-igraph": [""],
+        "pip+igraph": [""], // https://github.com/airspeed-velocity/asv/issues/1554
         // "psutil": [""]
         "pooch": [""],
-        "scikit-image": [""],
+        "scikit-image": [""], // https://github.com/conda-forge/scikit-misc-feedstock/pull/29
         // "scikit-misc": [""],
     },
 

benchmarks/benchmarks/preprocessing_counts.py

@@ -5,73 +5,64 @@
 
 from __future__ import annotations
 
+from itertools import product
 from typing import TYPE_CHECKING
 
+import anndata as ad
+
 import scanpy as sc
 
 from ._utils import get_count_dataset
 
 if TYPE_CHECKING:
-    from anndata import AnnData
-
     from ._utils import Dataset, KeyCount
 
-# setup variables
-
-adata: AnnData
-batch_key: str | None
-
-
-def setup(dataset: Dataset, layer: KeyCount, *_):
-    """Set up global variables before each benchmark."""
-    global adata, batch_key
-    adata, batch_key = get_count_dataset(dataset, layer=layer)
-    assert "log1p" not in adata.uns
-
 
 # ASV suite
+class PreprocessingCountsSuite:  # noqa: D101
+    params: tuple[list[Dataset], list[KeyCount]] = (
+        ["pbmc68k_reduced", "pbmc3k"],
+        ["counts", "counts-off-axis"],
+    )
+    param_names = ("dataset", "layer")
 
-params: tuple[list[Dataset], list[KeyCount]] = (
-    ["pbmc68k_reduced", "pbmc3k"],
-    ["counts", "counts-off-axis"],
-)
-param_names = ["dataset", "layer"]
-
-
-def time_filter_cells(*_):
-    sc.pp.filter_cells(adata, min_genes=100)
-
-
-def peakmem_filter_cells(*_):
-    sc.pp.filter_cells(adata, min_genes=100)
-
-
-def time_filter_genes(*_):
-    sc.pp.filter_genes(adata, min_cells=3)
-
+    def setup_cache(self) -> None:
+        """Without this caching, asv was running several processes which meant the data was repeatedly downloaded."""
+        for dataset, layer in product(*self.params):
+            adata, batch_key = get_count_dataset(dataset, layer=layer)
+            assert "lop1p" not in adata.uns
+            adata.uns["batch_key"] = batch_key
+            adata.write_h5ad(f"{dataset}_{layer}.h5ad")
 
-def peakmem_filter_genes(*_):
-    sc.pp.filter_genes(adata, min_cells=3)
+    def setup(self, dataset, layer) -> None:
+        self.adata = ad.read_h5ad(f"{dataset}_{layer}.h5ad")
 
+    def time_filter_cells(self, *_) -> None:
+        sc.pp.filter_cells(self.adata, min_genes=100)
 
-def time_scrublet(*_):
-    sc.pp.scrublet(adata, batch_key=batch_key)
+    def peakmem_filter_cells(self, *_) -> None:
+        sc.pp.filter_cells(self.adata, min_genes=100)
 
+    def time_filter_genes(self, *_) -> None:
+        sc.pp.filter_genes(self.adata, min_cells=3)
 
-def peakmem_scrublet(*_):
-    sc.pp.scrublet(adata, batch_key=batch_key)
+    def peakmem_filter_genes(self, *_) -> None:
+        sc.pp.filter_genes(self.adata, min_cells=3)
 
+    def time_scrublet(self, *_) -> None:
+        sc.pp.scrublet(self.adata, batch_key=self.adata.uns["batch_key"])
 
-# Can’t do seurat v3 yet: https://github.com/conda-forge/scikit-misc-feedstock/issues/17
-"""
-def time_hvg_seurat_v3(*_):
-    # seurat v3 runs on counts
-    sc.pp.highly_variable_genes(adata, flavor="seurat_v3_paper")
+    def peakmem_scrublet(self, *_) -> None:
+        sc.pp.scrublet(self.adata, batch_key=self.adata.uns["batch_key"])
 
+    # sciki-misc does not exit on osx-arm64
+    # https://github.com/conda-forge/scikit-misc-feedstock/pull/29
+    # def time_hvg_seurat_v3(self, *_):
+    #     # seurat v3 runs on counts
+    #     sc.pp.highly_variable_genes(self.adata, flavor="seurat_v3_paper")
 
-def peakmem_hvg_seurat_v3(*_):
-    sc.pp.highly_variable_genes(adata, flavor="seurat_v3_paper")
-"""
+    # def peakmem_hvg_seurat_v3(self, *_):
+    #     sc.pp.highly_variable_genes(self.adata, flavor="seurat_v3_paper")
 
 
 class FastSuite:
@@ -83,28 +74,38 @@ class FastSuite:
     )
     param_names = ("dataset", "layer")
 
-    def time_calculate_qc_metrics(self, *_):
+    def setup_cache(self) -> None:
+        """Without this caching, asv was running several processes which meant the data was repeatedly downloaded."""
+        for dataset, layer in product(*self.params):
+            adata, _ = get_count_dataset(dataset, layer=layer)
+            assert "lop1p" not in adata.uns
+            adata.write_h5ad(f"{dataset}_{layer}.h5ad")
+
+    def setup(self, dataset, layer) -> None:
+        self.adata = ad.read_h5ad(f"{dataset}_{layer}.h5ad")
+
+    def time_calculate_qc_metrics(self, *_) -> None:
         sc.pp.calculate_qc_metrics(
-            adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
+            self.adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
         )
 
-    def peakmem_calculate_qc_metrics(self, *_):
+    def peakmem_calculate_qc_metrics(self, *_) -> None:
         sc.pp.calculate_qc_metrics(
-            adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
+            self.adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
         )
 
-    def time_normalize_total(self, *_):
-        sc.pp.normalize_total(adata, target_sum=1e4)
+    def time_normalize_total(self, *_) -> None:
+        sc.pp.normalize_total(self.adata, target_sum=1e4)
 
-    def peakmem_normalize_total(self, *_):
-        sc.pp.normalize_total(adata, target_sum=1e4)
+    def peakmem_normalize_total(self, *_) -> None:
+        sc.pp.normalize_total(self.adata, target_sum=1e4)
 
-    def time_log1p(self, *_):
-        # TODO: This would fail: assert "log1p" not in adata.uns, "ASV bug?"
+    def time_log1p(self, *_) -> None:
+        # TODO: This would fail: assert "log1p" not in self.adata.uns, "ASV bug?"
         # https://github.com/scverse/scanpy/issues/3052
-        adata.uns.pop("log1p", None)
-        sc.pp.log1p(adata)
+        self.adata.uns.pop("log1p", None)
+        sc.pp.log1p(self.adata)
 
-    def peakmem_log1p(self, *_):
-        adata.uns.pop("log1p", None)
-        sc.pp.log1p(adata)
+    def peakmem_log1p(self, *_) -> None:
+        self.adata.uns.pop("log1p", None)
+        sc.pp.log1p(self.adata)

benchmarks/benchmarks/preprocessing_log.py

@@ -5,89 +5,70 @@
 
 from __future__ import annotations
 
+from itertools import product
 from typing import TYPE_CHECKING
 
+import anndata as ad
+
 import scanpy as sc
-from scanpy.preprocessing._utils import _get_mean_var
 
 from ._utils import get_dataset, param_skipper
 
 if TYPE_CHECKING:
-    from anndata import AnnData
-
     from ._utils import Dataset, KeyX
 
-# setup variables
-
-
-adata: AnnData
-batch_key: str | None
-
-
-def setup(dataset: Dataset, layer: KeyX, *_):
-    """Set up global variables before each benchmark."""
-    global adata, batch_key
-    adata, batch_key = get_dataset(dataset, layer=layer)
-
 
 # ASV suite
 
 params: tuple[list[Dataset], list[KeyX]] = (
     ["pbmc68k_reduced", "pbmc3k"],
     [None, "off-axis"],
 )
-param_names = ["dataset", "layer"]
-
+param_names = ("dataset", "layer")
 skip_when = param_skipper(param_names, params)
 
 
-def time_pca(*_):
-    sc.pp.pca(adata, svd_solver="arpack")
-
-
-def peakmem_pca(*_):
-    sc.pp.pca(adata, svd_solver="arpack")
-
-
-def time_highly_variable_genes(*_):
-    # the default flavor runs on log-transformed data
-    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
-
-
-def peakmem_highly_variable_genes(*_):
-    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
-
-
-# regress_out is very slow for this dataset
-@skip_when(dataset={"pbmc3k"})
-def time_regress_out(*_):
-    sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
-
+class PreprocessingSuite:  # noqa: D101
+    params = params
+    param_names = param_names
 
-@skip_when(dataset={"pbmc3k"})
-def peakmem_regress_out(*_):
-    sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
+    def setup_cache(self) -> None:
+        """Without this caching, asv was running several processes which meant the data was repeatedly downloaded."""
+        for dataset, layer in product(*self.params):
+            adata, _ = get_dataset(dataset, layer=layer)
+            adata.write_h5ad(f"{dataset}_{layer}.h5ad")
 
+    def setup(self, dataset, layer) -> None:
+        self.adata = ad.read_h5ad(f"{dataset}_{layer}.h5ad")
 
-def time_scale(*_):
-    sc.pp.scale(adata, max_value=10)
+    def time_pca(self, *_) -> None:
+        sc.pp.pca(self.adata, svd_solver="arpack")
 
+    def peakmem_pca(self, *_) -> None:
+        sc.pp.pca(self.adata, svd_solver="arpack")
 
-def peakmem_scale(*_):
-    sc.pp.scale(adata, max_value=10)
+    def time_highly_variable_genes(self, *_) -> None:
+        # the default flavor runs on log-transformed data
+        sc.pp.highly_variable_genes(
+            self.adata, min_mean=0.0125, max_mean=3, min_disp=0.5
+        )
 
+    def peakmem_highly_variable_genes(self, *_) -> None:
+        sc.pp.highly_variable_genes(
+            self.adata, min_mean=0.0125, max_mean=3, min_disp=0.5
+        )
 
-class FastSuite:
-    """Suite for fast preprocessing operations."""
+    # regress_out is very slow for this dataset
+    @skip_when(dataset={"pbmc3k"})
+    def time_regress_out(self, *_) -> None:
+        sc.pp.regress_out(self.adata, ["total_counts", "pct_counts_mt"])
 
-    params: tuple[list[Dataset], list[KeyX]] = (
-        ["pbmc3k", "pbmc68k_reduced", "bmmc", "lung93k"],
-        [None, "off-axis"],
-    )
-    param_names = ("dataset", "layer")
+    @skip_when(dataset={"pbmc3k"})
+    def peakmem_regress_out(self, *_) -> None:
+        sc.pp.regress_out(self.adata, ["total_counts", "pct_counts_mt"])
 
-    def time_mean_var(self, *_):
-        _get_mean_var(adata.X)
+    def time_scale(self, *_) -> None:
+        sc.pp.scale(self.adata, max_value=10)
 
-    def peakmem_mean_var(self, *_):
-        _get_mean_var(adata.X)
+    def peakmem_scale(self, *_) -> None:
+        sc.pp.scale(self.adata, max_value=10)