Add a per-row/per-column statistics sidecar to VCS arrays

README.md

@@ -143,6 +143,15 @@ va = vsparse.load_packed("archived.h5ad")
 norm = va.X.normalized("parafac2")  # lazy view; norm @ B works without densifying
 ```
 
+### Statistics sidecar
+
+Every `VCSCArray`/`VCSRArray` carries per-row and per-column summary statistics as `arr.stats` (a `vsparse.VCSStats`): `row_sums`/`col_sums` (total counts), `row_nnz`/`col_nnz` (n-genes/n-cells by counts), and `col_depth_sums` (per-gene sums of depth-normalized counts). They are computed on first access and written next to `X`/`raw_X` by `write_h5ad`/`write_zarr` whether or not they were ever computed or a normalized view was ever made, so a file read back has them without touching the data. Every recipe's per-cell and per-gene scales (`a`/`b`) are derived from them, so only centering/scaling recipes still need a pass over the data. Files without a usable sidecar (older versions, other tools) load as before, and the statistics are recomputed when first needed.
+
+```python
+va = vsparse.VCSCAnnData.read_h5ad("archived.h5ad")
+keep = va.X.stats.row_sums >= 500  # QC filter from the sidecar alone
+```
+
 ### CUDA (`to_gpu`)
 
 A normalized view moves onto an NVIDIA GPU with `.to_gpu()`, keeping its value-compressed layout there rather than expanding to one float per nonzero — so a matrix that only fits in device memory compressed still fits. Everything floating-point on the device is float32.

docs/api.md

@@ -13,6 +13,9 @@
    :show-inheritance:
    :inherited-members:
 
+.. autoclass:: vsparse.VCSStats
+   :members:
+
 .. autofunction:: vsparse.from_anndata
 
 .. autofunction:: vsparse.to_layer

docs/usage.md

@@ -127,6 +127,23 @@ va = vsparse.load_packed("archived.ivcsr.h5ad")
 norm = va.X.normalized("parafac2")
 ```
 
+### Statistics sidecar
+
+Each `VCSCArray`/`VCSRArray` carries a {class}`vsparse.VCSStats` as `arr.stats`:
+per-row `row_sums`/`row_nnz` (total counts, n-genes-by-counts), per-column
+`col_sums`/`col_nnz` (total counts, n-cells-by-counts), and `col_depth_sums`, the
+per-column sum of depth-normalized values `sum_i x[i, j] / row_sums[i]`. They are
+computed on first access (one pass over the data) and kept.
+
+`write_h5ad`/`write_zarr` always store them, as a versioned `stats` subgroup of
+`X`/`raw_X`, and readers (including `load_packed`) attach them, so QC filters on
+these columns need no read of `X`. Normalization takes every recipe's
+per-cell and per-gene scales (`a`/`b` in `y = (g(x * a * b) - c) * s`) from the
+sidecar; the centering/scaling terms `c`/`s` depend on the nonlinear `g`, so they
+still take one pass per recipe. A file without a usable sidecar -- written by an
+earlier version, or with an unrecognized sidecar version -- reads as before, and the
+statistics are recomputed when first needed.
+
 ## CUDA
 
 A normalized view can be moved onto an NVIDIA GPU with `.to_gpu()`, which needs

src/vsparse/__init__.py

@@ -12,6 +12,7 @@
 from vsparse._cuda import CudaNormalizedView, cuda_is_available
 from vsparse._norm_common import RECIPES, Recipe
 from vsparse._rapid_load import load_packed
+from vsparse._stats import VCSStats
 from vsparse._vcs_norm import VCSCArrayNormalized, VCSRArrayNormalized
 
 __all__ = [
@@ -23,6 +24,7 @@
     "VCSCArrayNormalized",
     "VCSRArray",
     "VCSRArrayNormalized",
+    "VCSStats",
     "cuda_is_available",
     "from_anndata",
     "load_packed",

src/vsparse/_base.py

@@ -17,7 +17,7 @@
 import numpy as np
 import scipy.sparse as sp
 
-from vsparse import _construct, _ops
+from vsparse import _construct, _ops, _stats
 from vsparse._indexutils import is_full_slice as _is_full_slice
 from vsparse._indexutils import normalize_major_idx as _normalize_major_idx
 from vsparse._indexutils import smallest_index_dtype as _smallest_index_dtype
@@ -42,7 +42,7 @@ class _VCSBase:
     # __rmatmul__/__rmul__ instead of trying to broadcast us as an ndarray.
     __array_ufunc__ = None
 
-    __slots__ = ("_norm_cache", "indices", "major_ptr", "shape", "value_ptr", "values")
+    __slots__ = ("_norm_cache", "_stats", "indices", "major_ptr", "shape", "value_ptr", "values")
 
     def __init__(
         self,
@@ -81,6 +81,7 @@ def __init__(
         self.value_ptr = value_ptr
         self.indices = indices
         self._norm_cache = _NormCache()
+        self._stats: _stats.VCSStats | None = None
 
     # -- axis bookkeeping ------------------------------------------------
 
@@ -117,14 +118,29 @@ def __repr__(self) -> str:  # pragma: no cover - cosmetic
             f"<{cls} shape={self.shape} dtype={self.dtype} nnz={self.nnz} n_unique={self.n_unique}>"
         )
 
+    @property
+    def stats(self) -> _stats.VCSStats:
+        """Per-row/per-column totals and nonzero counts -- see :mod:`vsparse._stats`.
+
+        Computed on first access (``O(nnz)``) and kept, or restored from the
+        sidecar of a file written by :meth:`~vsparse.VCSCAnnData.write_h5ad`/
+        :meth:`~vsparse.VCSCAnnData.write_zarr`. Normalization reads its
+        per-cell/per-gene scales (``a``/``b``) from here.
+        """
+        if self._stats is None:
+            self._stats = _stats.compute_stats(self)
+        return self._stats
+
     def copy(self):
-        return type(self)(
+        out = type(self)(
             self.shape,
             self.major_ptr.copy(),
             self.values.copy(),
             self.value_ptr.copy(),
             self.indices.copy(),
         )
+        out._stats = self._stats  # same data, so the same (read-only) statistics
+        return out
 
     # -- construction / conversion ---------------------------------------
 

src/vsparse/_construct.py

@@ -53,7 +53,9 @@ def _compress(major_ptr_in, minor_indices, data, n_major):
 
         seg_vals = data[start:end]
         seg_minor = minor_indices[start:end]
-        order = np.argsort(seg_vals)
+        # Stable, so each value group keeps its minor indices in input order
+        # (sorted, for canonical scipy input) -- see vsparse._stats.
+        order = np.argsort(seg_vals, kind="mergesort")
         sorted_vals = seg_vals[order]
         sorted_minor = seg_minor[order]
 

src/vsparse/_io.py

@@ -6,6 +6,12 @@
 :class:`~vsparse.VCSCAnnData` -- and survive a round trip through HDF5 or zarr
 without any manual (de)serialization.
 
+Both formats store the array's statistics sidecar (see :mod:`vsparse._stats`)
+as a ``stats`` subgroup next to the array data. Readers attach it when it's
+present and valid, and otherwise leave the array to recompute it on demand, so
+files without one -- written by earlier versions, or by other tools -- still
+read normally.
+
 This relies on ``anndata._io.specs.registry``, an internal (not yet public)
 extension mechanism -- the same one anndata itself uses to support
 dask/awkward/cupy arrays. It is not guaranteed stable across anndata
@@ -23,7 +29,7 @@
 import numpy as np
 import zarr
 
-from vsparse import _ivcsc
+from vsparse import _ivcsc, _stats
 from vsparse._base import VCSCArray, VCSRArray, _VCSBase
 from vsparse._compression import read_array
 from vsparse._indexutils import smallest_index_dtype
@@ -55,20 +61,31 @@ def _write_vcs(
     g.attrs["shape"] = v.shape
     for name in _ARRAY_KEYS:
         _writer.write_elem(g, name, getattr(v, name), dataset_kwargs=dataset_kwargs)
+    _stats.write_stats(g, v.stats, write_elem=_writer.write_elem, dataset_kwargs=dataset_kwargs)
+
+
+def _read_shape(elem: GroupStorageType) -> tuple[int, int]:
+    shape_vals = [int(s) for s in np.asarray(elem.attrs["shape"]).tolist()]
+    return (shape_vals[0], shape_vals[1])
+
+
+def _with_stats(arr: _VCSBase, elem: GroupStorageType) -> _VCSBase:
+    """Attach the sidecar stored under ``elem``, if there's a usable one -- see :mod:`vsparse._stats`."""
+    arr._stats = _stats.read_stats(elem, arr.shape)
+    return arr
 
 
 def _make_read_vcs(cls: type[_VCSBase]):
     def _read(elem: GroupStorageType, *, _reader: Reader) -> _VCSBase:
-        shape_vals = [int(s) for s in np.asarray(elem.attrs["shape"]).tolist()]
-        shape = (shape_vals[0], shape_vals[1])
         arrays = {name: read_array(elem[name]) for name in _ARRAY_KEYS}
-        return cls(
-            shape,
+        arr = cls(
+            _read_shape(elem),
             arrays["major_ptr"],
             arrays["values"],
             arrays["value_ptr"],
             arrays["indices"],
         )
+        return _with_stats(arr, elem)
 
     return _read
 
@@ -88,6 +105,9 @@ def write_ivcs_elem(
     the registry, e.g. via ``anndata.io.read_elem``) always reconstructs a
     normal VCSCArray/VCSRArray with a plain ``indices`` array.
 
+    ``v``'s statistics sidecar (:attr:`~vsparse.VCSCArray.stats`) is written
+    alongside, computing it first if ``v`` doesn't carry it yet.
+
     Unlike the VCSC/VCSR codec registered below, this is not reachable
     through anndata's generic ``write_elem`` type dispatch (which can only
     route by the Python type of ``v``, and both formats share the same
@@ -106,21 +126,21 @@ def write_ivcs_elem(
         ad.io.write_elem(g, name, getattr(v, name), dataset_kwargs=dataset_kwargs)
     packed = _ivcsc.pack_indices(v.value_ptr, v.indices)
     ad.io.write_elem(g, "packed_indices", packed, dataset_kwargs=dataset_kwargs)
+    _stats.write_stats(g, v.stats, write_elem=ad.io.write_elem, dataset_kwargs=dataset_kwargs)
     g.attrs["encoding-type"] = "ivcsc" if v._format == "csc" else "ivcsr"
     g.attrs["encoding-version"] = _IVCS_SPEC_VERSION
 
 
 def _make_read_ivcs(cls: type[_VCSBase]):
     def _read(elem: GroupStorageType, *, _reader: Reader) -> _VCSBase:
-        shape_vals = [int(s) for s in np.asarray(elem.attrs["shape"]).tolist()]
-        shape = (shape_vals[0], shape_vals[1])
         major_ptr = read_array(elem["major_ptr"])
         values = read_array(elem["values"])
         value_ptr = read_array(elem["value_ptr"])
         packed = read_array(elem["packed_indices"])
         dtype = np.dtype(cast(str, elem.attrs["indices_dtype"]))
         indices = _ivcsc.unpack_indices(value_ptr, packed, dtype)
-        return cls(shape, major_ptr, values, value_ptr, indices)
+        arr = cls(_read_shape(elem), major_ptr, values, value_ptr, indices)
+        return _with_stats(arr, elem)
 
     return _read