Files
cellxgene/backend/common/compute/estimate_distribution.py
T
Bruce Martin 1ea2b7fe80 fix for incorrect stats computation in diff exp t-test (#2318)
* 2211 fixes

* lint

* lint

* add missing test and bug found by test

* change terminology for count distribution

* update scanpy requirement

* update scanpy requirement
2021-07-23 11:36:26 -07:00

63 lines
1.9 KiB
Python

import numba
import concurrent.futures
import numpy as np
from scipy import sparse
from backend.common.constants import XApproxDistribution
@numba.njit(fastmath=True, error_model="numpy", nogil=True)
def min_max(arr):
"""Return (min, max) values for the ndarray."""
n = arr.size
odd = n % 2
if not odd:
n -= 1
max_val = min_val = arr[0]
i = 1
while i < n:
x = arr[i]
y = arr[i + 1]
if x > y:
x, y = y, x
min_val = min(x, min_val)
max_val = max(y, max_val)
i += 2
if not odd:
x = arr[n]
min_val = min(x, min_val)
max_val = max(x, max_val)
return min_val, max_val
def estimate_approximate_distribution(X) -> XApproxDistribution:
"""
Estimate the distribution (normal, count) of the X matrix.
Currently this is based upon the assumption that scRNA-seq data is
exponentially distributed in its raw (count) form, and when logged,
any (max-min) range in excess of 24 is implies tens of millions of
observations of a single feature and so is extremely unlikely.
"""
if sparse.isspmatrix_csc(X) or sparse.isspmatrix_csr(X):
Xdata = X.data
elif type(X) is np.ndarray:
Xdata = X.reshape(
X.size,
)
else:
raise TypeError(f"Unsupported matrix type: {str(type(X))}")
CHUNKSIZE = 1 << 24
if Xdata.size > CHUNKSIZE:
min_val = max_val = Xdata[0]
with concurrent.futures.ThreadPoolExecutor() as tp:
for (_min, _max) in tp.map(min_max, [Xdata[i : i + CHUNKSIZE] for i in range(0, Xdata.size, CHUNKSIZE)]):
min_val = min(_min, min_val)
max_val = max(_max, max_val)
else:
min_val, max_val = min_max(Xdata)
excess_range = (max_val - min_val) > 24
return XApproxDistribution.COUNT if excess_range else XApproxDistribution.NORMAL