fixed new xarray free handling in ttest and rank test

Author: davidsebfischer

diffxpy/stats/stats.py

@@ -1,8 +1,8 @@
-from typing import Union
-
 import numpy as np
 import numpy.linalg
+import scipy.sparse
 import scipy.stats
+from typing import Union
 
 
 def likelihood_ratio_test(
@@ -39,19 +39,18 @@ def likelihood_ratio_test(
 
 
 def mann_whitney_u_test(
-        x0: np.ndarray,
-        x1: np.ndarray,
+        x0: Union[np.ndarray, scipy.sparse.csr_matrix],
+        x1: Union[np.ndarray, scipy.sparse.csr_matrix]
 ):
     """
-    Perform Wilcoxon rank sum test (Mann-Whitney U test) along second axis
-    (for each gene).
+    Perform Wilcoxon rank sum test (Mann-Whitney U test) along second axis, ie. for each gene.
 
     The Wilcoxon rank sum test is a non-parameteric test
     to compare two groups of observations.
 
-    :param x0: np.array (observations x genes)
+    :param x0: (observations x genes)
         Observations in first group by gene
-    :param x1:  np.array (observations x genes)
+    :param x1: (observations x genes)
         Observations in second group by gene.
     """
     axis = 1
@@ -66,8 +65,8 @@ def mann_whitney_u_test(
 
     pvals = np.asarray([
         scipy.stats.mannwhitneyu(
-            x=x0[:, i].flatten(),
-            y=x1[:, i].flatten(),
+            x=np.asarray(x0[:, i].todense()).flatten() if isinstance(x0, scipy.sparse.csr_matrix) else x0[:, i],
+            y=np.asarray(x1[:, i].todense()).flatten() if isinstance(x0, scipy.sparse.csr_matrix) else x1[:, i],
             use_continuity=True,
             alternative="two-sided"
         ).pvalue for i in range(x0.shape[1])

diffxpy/testing/det.py

@@ -3,6 +3,7 @@
 from typing import Union, Dict, Tuple, List, Set
 import pandas as pd
 from random import sample
+import scipy.sparse
 
 import numpy as np
 import patsy
@@ -521,7 +522,7 @@ def _ave(self):
         :return: np.ndarray
         """
 
-        return np.mean(self.full_estim.x, axis=0)
+        return np.asarray(np.mean(self.full_estim.x, axis=0)).flatten()
 
     def _log_fold_change(self, factors: Union[Dict, Tuple, Set, List], base=np.e):
         """
@@ -715,16 +716,16 @@ def __init__(
         self._store_ols = None
 
         try:
-            if model_estim._error_codes is not None:
-                self._error_codes = model_estim._error_codes
+            if self.model_estim.error_codes is not None:
+                self._error_codes = self.model_estim.error_codes
             else:
                 self._error_codes = None
         except Exception as e:
             self._error_codes = None
 
         try:
-            if model_estim._niter is not None:
-                self._niter = model_estim._niter
+            if self.model_estim.niter is not None:
+                self._niter = self.model_estim.niter
             else:
                 self._niter = None
         except Exception as e:
@@ -776,7 +777,7 @@ def _ave(self):
 
         :return: np.ndarray
         """
-        return self.x.mean(axis=0)
+        return np.asarray(self.x.mean(axis=0)).flatten()
 
     def _test(self):
         """
@@ -1530,8 +1531,8 @@ def __init__(
         x0, x1 = split_x(data, grouping)
 
         # Only compute p-values for genes with non-zero observations and non-zero group-wise variance.
-        mean_x0 = x0.mean(axis=0).astype(dtype=np.float)
-        mean_x1 = x1.mean(axis=0).astype(dtype=np.float)
+        mean_x0 = np.asarray(np.mean(x0, axis=0)).flatten().astype(dtype=np.float)
+        mean_x1 = np.asarray(np.mean(x1, axis=0)).flatten().astype(dtype=np.float)
         # Avoid unnecessary mean computation:
         self._mean = np.average(
             a=np.vstack([mean_x0, mean_x1]),
@@ -1541,8 +1542,13 @@ def __init__(
             returned=False
         )
         self._ave_nonzero = self._mean != 0  # omit all-zero features
-        var_x0 = np.asarray(x0.var(axis=0)).flatten().astype(dtype=np.float)
-        var_x1 = np.asarray(x1.var(axis=0)).flatten().astype(dtype=np.float)
+        if isinstance(x0, scipy.sparse.csr_matrix):
+            # Efficient analytic expression of variance without densification.
+            var_x0 = np.asarray(np.mean(x0.power(2), axis=0)).flatten().astype(dtype=np.float) - np.square(mean_x0)
+            var_x1 = np.asarray(np.mean(x1.power(2), axis=0)).flatten().astype(dtype=np.float) - np.square(mean_x1)
+        else:
+            var_x0 = np.asarray(np.var(x0, axis=0)).flatten().astype(dtype=np.float)
+            var_x1 = np.asarray(np.var(x1, axis=0)).flatten().astype(dtype=np.float)
         self._var_geq_zero = np.logical_or(
             var_x0 > 0,
             var_x1 > 0
@@ -1649,8 +1655,8 @@ def __init__(
 
         x0, x1 = split_x(data, grouping)
 
-        mean_x0 = x0.mean(axis=0).astype(dtype=np.float)
-        mean_x1 = x1.mean(axis=0).astype(dtype=np.float)
+        mean_x0 = np.asarray(np.mean(x0, axis=0)).flatten().astype(dtype=np.float)
+        mean_x1 = np.asarray(np.mean(x1, axis=0)).flatten().astype(dtype=np.float)
         # Avoid unnecessary mean computation:
         self._mean = np.average(
             a=np.vstack([mean_x0, mean_x1]),
@@ -1659,19 +1665,24 @@ def __init__(
             axis=0,
             returned=False
         )
-        var_x0 = np.asarray(x0.var(axis=0)).flatten().astype(dtype=np.float)
-        var_x1 = np.asarray(x1.var(axis=0)).flatten().astype(dtype=np.float)
+        if isinstance(x0, scipy.sparse.csr_matrix):
+            # Efficient analytic expression of variance without densification.
+            var_x0 = np.asarray(np.mean(x0.power(2), axis=0)).flatten().astype(dtype=np.float) - np.square(mean_x0)
+            var_x1 = np.asarray(np.mean(x1.power(2), axis=0)).flatten().astype(dtype=np.float) - np.square(mean_x1)
+        else:
+            var_x0 = np.asarray(np.var(x0, axis=0)).flatten().astype(dtype=np.float)
+            var_x1 = np.asarray(np.var(x1, axis=0)).flatten().astype(dtype=np.float)
         self._var_geq_zero = np.logical_or(
             var_x0 > 0,
             var_x1 > 0
         )
         idx_run = np.where(np.logical_and(self._mean != 0, self._var_geq_zero))[0]
 
-        # TODO: can this be done on sparse?
+        # TODO: can this be done directly on sparse?
         pval = np.zeros([data.shape[1]]) + np.nan
         pval[idx_run] = stats.mann_whitney_u_test(
-            x0=np.asarray(x0[:, idx_run]),
-            x1=np.asarray(x1[:, idx_run])
+            x0=x0[:, idx_run],
+            x1=x1[:, idx_run]
         )
         pval[np.where(np.logical_and(
             np.logical_and(mean_x0 == mean_x1, self._mean > 0),