From 76523d4f32956579dc61002de427719cb482e5dd Mon Sep 17 00:00:00 2001 From: bmccandless Date: Tue, 2 Jun 2020 08:23:52 -0700 Subject: [PATCH] sparse column shift encoding. (#1502) Many of our matrices are log normalized, which tends to eliminate the number of non zero values (if there were any). This prevents the matrix from being stored as a sparse matrix. The solution here is to use a simple transformation to make it sparse again. The most common value from each column is subtracted from that column. These values that were subtracted are saved in an array called X_col_shift. The cellxgene code needs to understand how to undo the transformation when operating over the X matrix. - added script to create a synthetic dataset for testing - added a script to convert an existing CXG dataset to a sparse CXG dataset --- server/compute/diffexp_cxg.py | 6 ++ server/converters/cxgtool.py | 116 ++++++++++++++++++++++++----- server/converters/to_sparse.py | 62 +++++++++++++++ server/data_common/data_adaptor.py | 3 +- server/data_cxg/cxg_adaptor.py | 29 +++++--- server/test/__init__.py | 3 +- server/test/create_test_matrix.py | 44 +++++++++++ server/test/test_diffexp.py | 97 ++++++++++++++++-------- 8 files changed, 297 insertions(+), 63 deletions(-) create mode 100644 server/converters/to_sparse.py create mode 100644 server/test/create_test_matrix.py diff --git a/server/compute/diffexp_cxg.py b/server/compute/diffexp_cxg.py index 087d9c61..2a84abda 100644 --- a/server/compute/diffexp_cxg.py +++ b/server/compute/diffexp_cxg.py @@ -107,6 +107,12 @@ def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01): future.cancel() raise ComputeError(str(e)) + if is_sparse: + if adaptor.has_array("X_col_shift"): + X_col_shift = adaptor.open_array("X_col_shift")[:] + meanA += X_col_shift + meanB += X_col_shift + r = diffexp_ttest_from_mean_var( meanA.astype(dtype), varA.astype(dtype), diff --git a/server/converters/cxgtool.py b/server/converters/cxgtool.py index 3290740c..504e99b6 100644 --- a/server/converters/cxgtool.py +++ b/server/converters/cxgtool.py @@ -8,8 +8,11 @@ The organization of the TileDB structure is: ├─ obs TileDB array containing cell (row) attributes, one attribute per │ dataframe column, shape (n_obs,) ├─ var TileDB array containing gene (column) attributes, with one attribute per - │ dataframe column, shape (n_obs,) + │ dataframe column, shape (n_var,) ├─ X Main count matrix as a 2D TileDB array, single unnamed numeric attribute + ├─ X_col_shift TilebDB Array used in column shift encoding, shape (n_var,), dtype = X.dtype. + │ Single unnamed numeric attribute. If this array is sparse, and X_col_shift exists, + │ then all values in the i'th column were subtracted by X_col_shift[i]. ├─ emb TileDB group, storing optional embeddings (group may be empty) │ └─ TileDB Array, single anon attribute, ND numeric array, shape (n_obs, N) └─ cxg_group_metadata Empty array used only to stash metadata about the overall object. @@ -70,6 +73,7 @@ import argparse import numpy as np from os.path import splitext, basename import json +from scipy.stats import mode from server.common.colors import convert_anndata_category_colors_to_cxg_category_colors from server.common.errors import ColorFormatException @@ -205,7 +209,7 @@ def write_cxg( log(1, "\t...embeddings created") # X matrix - save_X(container, adata, ctx, sparse_threshold) + save_X(container, adata.X, ctx, sparse_threshold) log(1, "\t...X created") @@ -437,72 +441,146 @@ def create_X(X_name, shape, is_sparse): tiledb.DenseArray.create(X_name, schema) -def check_sparse_X(adata, sparse_threshold): - shape = adata.X.shape +def evaluate_for_sparse_encoding(xdata, sparse_threshold): + """ + This function determines if the X matrix has a sparsity below the sparse_threshold. + This function also returns the number of non-zeros encountered and number + of elements evaluated. This function may return before evaluating the whole X matrix + if it can be determined that X is not sparse enough. + """ + shape = xdata.shape stride = min(int(np.power(10, np.around(np.log10(1e9 / shape[1])))), 10_000) nnz = 0 maxnnz = int(shape[0] * shape[1] * sparse_threshold / 100) for row in range(0, shape[0], stride): lim = min(row + stride, shape[0]) - a = adata.X[row:lim, :] + a = xdata[row:lim, :] if type(a) is not np.ndarray: a = a.toarray() nnz += np.count_nonzero(a) if nnz > maxnnz: - return (nnz, lim * shape[1]) - log(2, "\t...rows", row, "to", lim, "nnz", nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[1]))) + return (False, nnz, lim * shape[1]) + log(2, "\t...rows", lim, "of", shape[0], "nnz", nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[1]))) - return (nnz, shape[0] * shape[1]) + is_sparse = (100.0 * nnz / (shape[0] * shape[1])) < sparse_threshold + return (is_sparse, nnz, shape[0] * shape[1]) -def save_X(container, adata, ctx, sparse_threshold): +def evaluate_for_sparse_column_shift_encoding(xdata, sparse_threshold): + """Column shift encoding works by taking the most common value in each column, then + subtracting that value from each element of the column. If each column mostly contains + its most common value, then the resulting matrix can be very sparse. + + This function determines if column shift encoding can be used to transform + the X matrix into a sparse matrix with a sparsity below the sparse_threshold. + If so, return the col_shift array that stores this encoding. + This function also returns the number of non-zeros encountered and number + of elements evaluated. This function may return before evaluating the whole X matrix + if it can be determined that X cannot benefit from column shift encoding. + """ + shape = xdata.shape + stride = max(1, 128_000_000 // shape[0]) + col_shift = np.zeros(shape[1]) + nnz = 0 + maxnnz = int(shape[0] * shape[1] * sparse_threshold / 100) + for col in range(0, shape[1], stride): + lim = min(col + stride, shape[1]) + a = xdata[:, col:lim] + if type(a) is not np.ndarray: + a = a.toarray() + m = mode(a) + col_shift[col:lim] = m.mode + nnz += shape[0] * (lim - col) - np.sum(m.count) + if nnz > maxnnz: + return (None, nnz, shape[0] * lim) + log(2, "\t...cols", lim, "of", shape[1], "nnz", + nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[0]))) + + is_sparse = (100.0 * nnz / (shape[0] * shape[1])) < sparse_threshold + return (col_shift if is_sparse else None, nnz, shape[0] * shape[1]) + + +def save_X(container, xdata, ctx, sparse_threshold, expect_sparse=False): # Save X count matrix X_name = f"{container}/X" - shape = adata.X.shape + shape = xdata.shape + log(1, "\t...shape:", str(shape)) + + col_shift = None if sparse_threshold == 100: is_sparse = True elif sparse_threshold == 0: is_sparse = False else: - nnz, nelem = check_sparse_X(adata, sparse_threshold) + is_sparse, nnz, nelem = evaluate_for_sparse_encoding(xdata, sparse_threshold) percent = 100.0 * nnz / nelem if nelem != shape[0] * shape[1]: - log(1, "\t...shape:", str(shape), "non-zeros percent (estimate): %6.2f" % percent) + log(1, "\t...sparse=", is_sparse, "non-zeros percent (estimate): %6.2f" % percent) else: - log(1, "\t...shape:", str(shape), "non-zeros:", nnz, "percent: %6.2f" % percent) + log(1, "\t...sparse=", is_sparse, "non-zeros:", nnz, "percent: %6.2f" % percent) is_sparse = percent < sparse_threshold + if not is_sparse: + col_shift, nnz, nelem = evaluate_for_sparse_column_shift_encoding(xdata, sparse_threshold) + is_sparse = col_shift is not None + percent = 100.0 * nnz / nelem + if nelem != shape[0] * shape[1]: + log(1, "\t...sparse=", is_sparse, "col shift non-zeros percent (estimate): %6.2f" % percent) + else: + log(1, "\t...sparse=", is_sparse, "col shift non-zeros:", nnz, "percent: %6.2f" % percent) + + if expect_sparse is True and is_sparse is False: + return False create_X(X_name, shape, is_sparse) - stride = min(int(np.power(10, np.around(np.log10(1e9 / shape[1])))), 10_000) if is_sparse: - log(1, "\t...output X as sparse matrix") + if col_shift is not None: + log(1, "\t...output X as sparse matrix with column shift encoding") + X_col_shift_name = f"{container}/X_col_shift" + filters = tiledb.FilterList([tiledb.ZstdFilter()]) + attrs = [tiledb.Attr(dtype=np.float32, filters=filters)] + domain = tiledb.Domain(tiledb.Dim(domain=(0, shape[1] - 1), tile=min(shape[1], 5000), dtype=np.uint32)) + schema = tiledb.ArraySchema(domain=domain, attrs=attrs) + tiledb.DenseArray.create(X_col_shift_name, schema) + with tiledb.DenseArray(X_col_shift_name, mode="w", ctx=ctx) as X_col_shift: + X_col_shift[:] = col_shift + tiledb.consolidate(X_col_shift_name, ctx=ctx) + else: + log(1, "\t...output X as sparse matrix") + with tiledb.SparseArray(X_name, mode="w", ctx=ctx) as X: + nnz = 0 for row in range(0, shape[0], stride): lim = min(row + stride, shape[0]) - a = adata.X[row:lim, :] + a = xdata[row:lim, :] if type(a) is not np.ndarray: a = a.toarray() + if col_shift is not None: + a = a - col_shift indices = np.nonzero(a) trow = indices[0] + row + nnz += indices[0].shape[0] X[trow, indices[1]] = a[indices[0], indices[1]] - log(2, "\t...rows", row, "to", lim) - tiledb.consolidate(X_name, ctx=ctx) + log(2, "\t...rows", lim, "of", shape[0], "nnz", nnz, "sparse", nnz / (lim * shape[1])) else: log(1, "\t...output X as dense matrix") with tiledb.DenseArray(X_name, mode="w", ctx=ctx) as X: for row in range(0, shape[0], stride): lim = min(row + stride, shape[0]) - a = adata.X[row:lim, :] + a = xdata[row:lim, :] if type(a) is not np.ndarray: a = a.toarray() X[row:lim, :] = a log(2, "\t...rows", row, "to", lim) tiledb.consolidate(X_name, ctx=ctx) + if hasattr(tiledb, "vacuum"): + tiledb.vacuum(X_name) + + return is_sparse def save_metadata(container, metadata_dict): diff --git a/server/converters/to_sparse.py b/server/converters/to_sparse.py new file mode 100644 index 00000000..ef704c4d --- /dev/null +++ b/server/converters/to_sparse.py @@ -0,0 +1,62 @@ +""" +Script to create a sparse dataset in CXG format based on an input dataset in CXG format. +The input dataset is not modified. +""" +import os +import shutil +import tiledb +import argparse +import sys +import server.converters.cxgtool as cxgtool + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("input", help="input cxg directory") + parser.add_argument("output", help="output cxg directory") + parser.add_argument("--overwrite", action="store_true", help="replace output cxg directory") + parser.add_argument("--verbose", "-v", action="count", default=0, help="verbose output") + parser.add_argument( + "--sparse-threshold", + "-s", + type=float, + default=5.0, # default is 5% non-zero values + help="The X array will be sparse if the percent of non-zeros falls below this value", + ) + args = parser.parse_args() + + if os.path.exists(args.output): + print("output dir exists:", args.output) + if args.overwrite: + print("output dir removed:", args.output) + shutil.rmtree(args.output) + else: + print("use the overwrite option to remove the output directory") + sys.exit(1) + + if not os.path.isdir(args.input): + print("input is not a directory", args.input) + sys.exit(1) + + shutil.copytree(args.input, args.output, + ignore=shutil.ignore_patterns("X", "X_col_shift")) + + ctx = tiledb.Ctx( + { + "sm.num_reader_threads": 32, + "sm.num_writer_threads": 32, + "sm.consolidation.buffer_size": 1 * 1024 * 1024 * 1024, + } + ) + + with tiledb.DenseArray(os.path.join(args.input, "X"), mode="r", ctx=ctx) as X_in: + is_sparse = cxgtool.save_X(args.output, X_in, ctx, args.sparse_threshold, expect_sparse=True) + + if is_sparse is False: + print("The array is not sparse, cleaning up, abort.") + shutil.rmtree(args.output) + sys.exit(1) + + +if __name__ == "__main__": + main() diff --git a/server/data_common/data_adaptor.py b/server/data_common/data_adaptor.py index cd542cc0..617addad 100644 --- a/server/data_common/data_adaptor.py +++ b/server/data_common/data_adaptor.py @@ -241,8 +241,7 @@ class DataAdaptor(metaclass=ABCMeta): duplicate_columns = list(set(labels_df.columns) & set(obs_columns)) if len(duplicate_columns) > 0: raise KeyError( - "Labels file may not contain column names which overlap " - f"with h5ad obs columns {duplicate_columns}" + "Labels file may not contain column names which overlap " f"with h5ad obs columns {duplicate_columns}" ) # labels must have same count as obs annotations diff --git a/server/data_cxg/cxg_adaptor.py b/server/data_cxg/cxg_adaptor.py index 968be4eb..a5a75813 100644 --- a/server/data_cxg/cxg_adaptor.py +++ b/server/data_cxg/cxg_adaptor.py @@ -133,6 +133,10 @@ class CxgAdaptor(DataAdaptor): return False return True + def has_array(self, name): + a_type = tiledb.object_type(path_join(self.url, name), ctx=self.tiledb_ctx) + return a_type == "array" + def _validate_and_initialize(self): """ remember, preload_validation() has already been called, so @@ -147,13 +151,7 @@ class CxgAdaptor(DataAdaptor): * version 0.1 -- metadata attache to cxg_group_metadata array. Same as 0, except it adds group metadata. """ - a_type = tiledb.object_type(path_join(self.url, "cxg_group_metadata"), ctx=self.tiledb_ctx) - if a_type is None: - # version 0 - cxg_version = "0.0" - title = None - about = None - elif a_type == "array": + if self.has_array("cxg_group_metadata"): # version >0 gmd = self.open_array("cxg_group_metadata") cxg_version = gmd.meta["cxg_version"] @@ -161,6 +159,11 @@ class CxgAdaptor(DataAdaptor): cxg_properties = json.loads(gmd.meta["cxg_properties"]) title = cxg_properties.get("title", None) about = cxg_properties.get("about", None) + else: + # version 0 + cxg_version = "0.0" + title = None + about = None if cxg_version not in ["0.0", "0.1"]: raise DatasetAccessError(f"cxg matrix is not valid: {self.url}") @@ -251,10 +254,18 @@ class CxgAdaptor(DataAdaptor): data = X[:, :] else: data = X.multi_index[obs_items, var_items] - nrows, obsindices = self.__remap_indices(X.shape[0], obs_mask, data["obs"]) - ncols, varindices = self.__remap_indices(X.shape[1], var_mask, data["var"]) + + nrows, obsindices = self.__remap_indices(X.shape[0], obs_mask, data.get("coords", data)["obs"]) + ncols, varindices = self.__remap_indices(X.shape[1], var_mask, data.get("coords", data)["var"]) densedata = np.zeros((nrows, ncols), dtype=self.get_X_array_dtype()) densedata[obsindices, varindices] = data[""] + if self.has_array("X_col_shift"): + X_col_shift = self.open_array("X_col_shift") + if var_items == slice(None): + densedata += X_col_shift[:] + else: + densedata += X_col_shift.multi_index[var_items][""] + return densedata else: diff --git a/server/test/__init__.py b/server/test/__init__.py index 38d471c7..5b1b36be 100644 --- a/server/test/__init__.py +++ b/server/test/__init__.py @@ -60,7 +60,7 @@ def skip_if(condition, reason: str): return decorator -def app_config(data_locator, backed=False): +def app_config(data_locator, backed=False, extra={}): args = { "embeddings__names": ["umap", "tsne", "pca"], "presentation__max_categories": 100, @@ -74,6 +74,7 @@ def app_config(data_locator, backed=False): } config = AppConfig() config.update(**args) + config.update(**extra) config.complete_config() return config diff --git a/server/test/create_test_matrix.py b/server/test/create_test_matrix.py new file mode 100644 index 00000000..476b6387 --- /dev/null +++ b/server/test/create_test_matrix.py @@ -0,0 +1,44 @@ +import anndata +import argparse +import random +import scipy +import numpy as np + + +def main(): + parser = argparse.ArgumentParser("A command to generate test h5ad files") + parser.add_argument("output", help="Name of the output file") + parser.add_argument("nobs", type=int, help="Number of observations (rows)") + parser.add_argument("nvar", type=int, help="Number of variables (columns)") + parser.add_argument("-n", "--nnz-percent", type=float, default=100, help="percent of non-zeros") + parser.add_argument("-c", "--col-shift", action="store_true", help="add a random value to each column") + parser.add_argument("--seed", type=int, default=None, help="add a random value to each column") + + args = parser.parse_args() + create_test_h5ad(args.output, args.nobs, args.nvar, args.nnz_percent, args.col_shift, args.seed) + + +def create_test_h5ad(outfile, nobs, nvar, nnz_percent=100, apply_col_shift=False, seed=None): + random.seed(seed) + np.random.seed(seed) + x = create_X_array(nobs, nvar, nnz_percent, apply_col_shift) + obsm = {"X_random": np.random.rand(nobs, 2).astype(np.float32)} + adata = anndata.AnnData(x, obsm=obsm) + adata.write(outfile) + + +def create_X_array(nobs, nvar, nnz_percent, apply_col_shift): + if nnz_percent < 100: + array = scipy.sparse.random(nobs, nvar, nnz_percent * 0.01, dtype=np.float32, format="csc") + else: + array = np.random.rand(nobs, nvar).astype(np.float32) + + if apply_col_shift: + col_shift = np.random.rand((nvar)) + array += col_shift + + return array + + +if __name__ == "__main__": + main() diff --git a/server/test/test_diffexp.py b/server/test/test_diffexp.py index a5d222f7..c69365ac 100644 --- a/server/test/test_diffexp.py +++ b/server/test/test_diffexp.py @@ -4,9 +4,10 @@ from server.test import PROJECT_ROOT, app_config import server.compute.diffexp_cxg as diffexp_cxg import server.compute.diffexp_generic as diffexp_generic from server.converters.cxgtool import write_cxg +from server.test.create_test_matrix import create_test_h5ad +from server.data_common.fbs.matrix import encode_matrix_fbs, decode_matrix_fbs import numpy as np import tempfile -import anndata import os @@ -14,8 +15,8 @@ class DiffExpTest(unittest.TestCase): """Tests the diffexp returns the expected results for one test case, using different adaptor types and different algorithms.""" - def load_dataset(self, path): - config = app_config(path) + def load_dataset(self, path, extra={}): + config = app_config(path, extra=extra) loader = MatrixDataLoader(path) adaptor = loader.open(config) return adaptor @@ -28,6 +29,14 @@ class DiffExpTest(unittest.TestCase): mask[sel] = True return mask + def compare_diffexp_results(self, results, expects): + self.assertEqual(len(results), len(expects)) + for result, expect in zip(results, expects): + self.assertEqual(result[0], expect[0]) + self.assertTrue(np.isclose(result[1], expect[1], 1e-6, 1e-4)) + self.assertTrue(np.isclose(result[2], expect[2], 1e-6, 1e-4)) + self.assertTrue(np.isclose(result[3], expect[3], 1e-6, 1e-4)) + def check_1_10_2_10(self, results): """Checks the results for a specific set of rows selections""" expects = [ @@ -42,12 +51,7 @@ class DiffExpTest(unittest.TestCase): [1575, 1.0317602, 0.007830310753043345, 1.0], [576, 0.97873515, 0.008272092578813124, 1.0], ] - self.assertEqual(len(results), len(expects)) - for result, expect in zip(results, expects): - self.assertEqual(result[0], expect[0]) - self.assertTrue(np.isclose(result[1], expect[1], 1e-6, 1e-6)) - self.assertTrue(np.isclose(result[2], expect[2], 1e-6, 1e-6)) - self.assertTrue(np.isclose(result[3], expect[3], 1e-6, 1e-6)) + self.compare_diffexp_results(results, expects) def get_X_col(self, adaptor, cols): varmask = np.zeros(adaptor.get_shape()[1], dtype=bool) @@ -86,33 +90,62 @@ class DiffExpTest(unittest.TestCase): self.check_1_10_2_10(results) def test_cxg_sparse(self): + self.sparse_diffexp(False) + + def test_cxg_sparse_col_shift(self): + self.sparse_diffexp(True) + + def sparse_diffexp(self, apply_col_shift): with tempfile.TemporaryDirectory() as dirname: + # create a sparse matrix + h5adfile = os.path.join(dirname, "sparse.h5ad") + create_test_h5ad(h5adfile, 2000, 2000, 10, apply_col_shift) + adaptor_anndata = self.load_dataset(h5adfile, extra=dict(embeddings__names=[])) + adata = adaptor_anndata.data + sparsename = os.path.join(dirname, "sparse.cxg") - densename = os.path.join(dirname, "dense.cxg") - source_h5ad = anndata.read_h5ad(f"{PROJECT_ROOT}/example-dataset/pbmc3k.h5ad") - # create a cxg sparse array - write_cxg(adata=source_h5ad, container=sparsename, title="pbmc3k", sparse_threshold=100) - write_cxg(adata=source_h5ad, container=densename, title="pbmc3k", sparse_threshold=0) + write_cxg(adata=adata, container=sparsename, title="sparse", sparse_threshold=11) adaptor_sparse = self.load_dataset(sparsename) + assert adaptor_sparse.open_array("X").schema.sparse + assert adaptor_sparse.has_array("X_col_shift") == apply_col_shift + + densename = os.path.join(dirname, "dense.cxg") + write_cxg(adata=adata, container=densename, title="dense", sparse_threshold=0) adaptor_dense = self.load_dataset(densename) + assert not adaptor_dense.open_array("X").schema.sparse + assert not adaptor_dense.has_array("X_col_shift") - col_results = [] - for adaptor in (adaptor_sparse, adaptor_dense): - maskA = self.get_mask(adaptor, 1, 10) - maskB = self.get_mask(adaptor, 2, 10) - diffexp_results = diffexp_cxg.diffexp_ttest(adaptor, maskA, maskB, 10) - self.check_1_10_2_10(diffexp_results) - topcols = [x[0] for x in diffexp_results] - cols = self.get_X_col(adaptor, topcols) - assert cols.shape[0] == adaptor.get_shape()[0] - assert cols.shape[1] == len(diffexp_results) - col_results.append(cols) + maskA = self.get_mask(adaptor_anndata, 1, 10) + maskB = self.get_mask(adaptor_anndata, 2, 10) - x = adaptor.get_X_array() - print(x) + diffexp_results_anndata = diffexp_generic.diffexp_ttest(adaptor_anndata, maskA, maskB, 10) + diffexp_results_sparse = diffexp_cxg.diffexp_ttest(adaptor_sparse, maskA, maskB, 10) + diffexp_results_dense = diffexp_cxg.diffexp_ttest(adaptor_dense, maskA, maskB, 10) - for row in range(col_results[0].shape[0]): - for col in range(col_results[0].shape[1]): - sval = col_results[0][row][col] - dval = col_results[1][row][col] - self.assertTrue(np.isclose(sval, dval, 1e-6, 1e-6)) + self.compare_diffexp_results(diffexp_results_anndata, diffexp_results_sparse) + self.compare_diffexp_results(diffexp_results_anndata, diffexp_results_dense) + + topcols = np.array([x[0] for x in diffexp_results_anndata]) + cols_anndata = self.get_X_col(adaptor_anndata, topcols) + cols_sparse = self.get_X_col(adaptor_sparse, topcols) + cols_dense = self.get_X_col(adaptor_dense, topcols) + assert cols_anndata.shape[0] == adaptor_sparse.get_shape()[0] + assert cols_anndata.shape[1] == len(diffexp_results_anndata) + + def convert(mat, cols): + return decode_matrix_fbs(encode_matrix_fbs(mat, col_idx=cols)).to_numpy() + + cols_anndata = convert(cols_anndata, topcols) + cols_sparse = convert(cols_sparse, topcols) + cols_dense = convert(cols_dense, topcols) + + x = adaptor_sparse.get_X_array() + assert x.shape == adaptor_sparse.get_shape() + + for row in range(cols_anndata.shape[0]): + for col in range(cols_anndata.shape[1]): + vanndata = cols_anndata[row][col] + vsparse = cols_sparse[row][col] + vdense = cols_dense[row][col] + self.assertTrue(np.isclose(vanndata, vsparse, 1e-6, 1e-6)) + self.assertTrue(np.isclose(vanndata, vdense, 1e-6, 1e-6))