Re-Added to_numpy()

Author: ilaykav

benchmarks/benchmarks/preprocessing_log.py

@@ -9,10 +9,11 @@
 from typing import TYPE_CHECKING
 
 import anndata as ad
+import numpy as np
 
 import scanpy as sc
 
-from ._utils import get_dataset, param_skipper
+from ._utils import get_dataset, param_skipper, pbmc3k
 
 if TYPE_CHECKING:
     from ._utils import Dataset, KeyX
@@ -72,3 +73,26 @@ def time_scale(self, *_) -> None:
 
     def peakmem_scale(self, *_) -> None:
         sc.pp.scale(self.adata, max_value=10)
+
+
+class CombatSuite:
+    """Benchmark combat batch correction."""
+
+    def setup_cache(self) -> None:
+        adata = pbmc3k()
+        sc.pp.highly_variable_genes(adata, n_top_genes=500)
+        adata = adata[:, adata.var["highly_variable"]].copy()
+        sc.pp.scale(adata, max_value=10)
+        # assign cells to 3 batches deterministically
+        rng = np.random.default_rng(0)
+        adata.obs["batch"] = rng.choice(["A", "B", "C"], size=adata.n_obs)
+        adata.write_h5ad("adata_combat.h5ad")
+
+    def setup(self) -> None:
+        self.adata = ad.read_h5ad("adata_combat.h5ad")
+
+    def time_combat(self) -> None:
+        sc.pp.combat(self.adata, key="batch")
+
+    def peakmem_combat(self) -> None:
+        sc.pp.combat(self.adata, key="batch")

benchmarks/benchmarks/tools.py

@@ -9,7 +9,7 @@
 
 import scanpy as sc
 
-from ._utils import pbmc3k, pbmc68k_reduced
+from ._utils import pbmc68k_reduced
 
 
 class ToolsSuite:  # noqa: D101
@@ -44,28 +44,3 @@ def time_rank_genes_groups(self) -> None:
 
     def peakmem_rank_genes_groups(self) -> None:
         sc.tl.rank_genes_groups(self.adata, "bulk_labels", method="wilcoxon")
-
-
-class CombatSuite:
-    """Benchmark combat batch correction."""
-
-    def setup_cache(self) -> None:
-        import numpy as np
-
-        adata = pbmc3k()
-        sc.pp.highly_variable_genes(adata, n_top_genes=500)
-        adata = adata[:, adata.var["highly_variable"]].copy()
-        sc.pp.scale(adata, max_value=10)
-        # assign cells to 3 batches deterministically
-        np.random.seed(0)
-        adata.obs["batch"] = np.random.choice(["A", "B", "C"], size=adata.n_obs)
-        adata.write_h5ad("adata_combat.h5ad")
-
-    def setup(self) -> None:
-        self.adata = ad.read_h5ad("adata_combat.h5ad")
-
-    def time_combat(self) -> None:
-        sc.pp.combat(self.adata, key="batch")
-
-    def peakmem_combat(self) -> None:
-        sc.pp.combat(self.adata, key="batch")

src/scanpy/preprocessing/_combat.py

@@ -105,21 +105,17 @@ def _standardize_data(
 
     design = _design_matrix(model, batch_key, batch_levels)
 
-    # use numpty .values extration only once to avoid pandas overhead
-    design_arr = design.values
     # compute pooled variance estimator
-    b_hat = np.dot(
-        np.dot(la.inv(np.dot(design_arr.T, design_arr)), design_arr.T), data.values.T
-    )
+    b_hat = np.dot(np.dot(la.inv(np.dot(design.T, design)), design.T), data.T)
     grand_mean = np.dot((n_batches / n_array).T, b_hat[:n_batch, :])
-    var_pooled = (data.values - np.dot(design_arr, b_hat).T) ** 2
+    var_pooled = np.asarray((data - np.dot(design, b_hat).T) ** 2)
     var_pooled = np.mean(var_pooled, axis=1, keepdims=True)
 
     # Compute the means
     if np.sum(var_pooled == 0) > 0:
         print(f"Found {np.sum(var_pooled == 0)} genes with zero variance.")
-    stand_mean = grand_mean[:, np.newaxis]
-    tmp = design_arr.copy()
+    stand_mean = np.asarray(grand_mean)[:, np.newaxis]
+    tmp = np.array(design.copy())
     tmp[:, :n_batch] = 0
     stand_mean = stand_mean + np.dot(tmp, b_hat).T
 
@@ -128,7 +124,7 @@ def _standardize_data(
     s_data = np.where(
         var_pooled == 0,
         0,
-        (data.values - stand_mean) / np.sqrt(var_pooled),
+        (np.asarray(data) - stand_mean) / np.sqrt(var_pooled),
     )
     s_data = pd.DataFrame(s_data, index=data.index, columns=data.columns)
 
@@ -272,27 +268,24 @@ def combat(  # noqa: PLR0915
 
     # we now apply the parametric adjustment to the standardized data from above
     # loop over all batches in the data
-    bayesdata_arr = bayesdata.to_numpy(copy=True)
-    batch_design_arr = batch_design.values
     for j, batch_idxs in enumerate(batch_info.values()):
         # we basically subtract the additive batch effect, rescale by the ratio
         # of multiplicative batch effect to pooled variance and add the overall gene
         # wise mean
         dsq = np.sqrt(delta_star[j, :])
-        numer = (
-            bayesdata_arr[:, batch_idxs]
-            - np.dot(batch_design_arr[batch_idxs], gamma_star).T
+        numer = np.array(
+            bayesdata.iloc[:, batch_idxs]
+            - np.dot(batch_design.iloc[batch_idxs], gamma_star).T
         )
-        bayesdata_arr[:, batch_idxs] = numer / dsq[:, np.newaxis]
+        bayesdata.iloc[:, batch_idxs] = numer / dsq[:, np.newaxis]
 
-    bayesdata_arr = bayesdata_arr * np.sqrt(var_pooled) + stand_mean
+    bayesdata = bayesdata * np.sqrt(var_pooled) + stand_mean
 
     # put back into the adata object or return
-    x = bayesdata.to_numpy().transpose()
     if inplace:
-        adata.X = bayesdata_arr.T
+        adata.X = bayesdata.to_numpy().transpose()
         return None
-    return bayesdata_arr.T
+    return bayesdata.to_numpy().transpose()
 
 
 def _it_sol(