From b65c32168c347590032aeb2466f8d7ae16373ad1 Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Sun, 13 Sep 2026 12:24:31 -0700
Subject: [PATCH 1/6] Keep BiCV's held-out cell blocks lazy for duck-typed
 backends

_bicv_trial evaluated a rank's fit by row-restricting X.X with a boolean
cell mask (X_mat[train_cell_mask]/X_mat[test_cell_mask]) and handing the
result to rmatmul/calc_W, deliberately avoiding materializing the raw data
("reaches the raw data through products... rather than materialising a
block of it"). That holds for a plain ndarray or scipy-sparse X_mat, but
not for a vsparse normalized view (e.g. BAL-Pf2's lazy-normalized-view
AnnData): bracket indexing such a view is documented to always eagerly
build a dense ndarray for the selection. At BAL-Pf2's real scale (1.3M
cells), a single train/test split (~50% of all cells) would materialize
on the order of tens of GB, once or twice per BiCV trial, across every
rank/repeat in a sweep -- never surfaced before because an unrelated
vsparse memory issue always killed these runs earlier in the pipeline.

Adds _restrict_rows(X_mat, mask), which uses vsparse's new
select(recalculate=False) (meyer-lab/vsparse#50) to stay a genuinely lazy
view for a duck-typed backend, falling back to ordinary indexing for
plain dense/sparse X_mat (unaffected, still cheap for those). Also adds a
third branch to _test_block_moments for the same duck-typed case: selects
the cell subset lazily, then streams it in bounded row chunks rather than
ever materializing the whole subset as one dense block.

Temporarily pins vsparse to its (as yet unmerged/unreleased)
select-without-recalculate branch in pyproject.toml -- drop once
meyer-lab/vsparse#50 merges and releases, reverting to the plain PyPI
version constraint.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 pyproject.toml                      |   8 ++
 scrise/rank_selection.py            |  60 +++++++++++--
 scrise/tests/test_rank_selection.py | 133 ++++++++++++++++++++++++++++
 uv.lock                             |  10 +--
 4 files changed, 197 insertions(+), 14 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 0852ae0c..00d17ecc 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,6 +39,14 @@ fbuild = "analysis.figures.common:genFigure"
 requires = ["uv_build>=0.12.0,<0.13"]
 build-backend = "uv_build"
 
+# TEMPORARY: pin vsparse to an unreleased PR branch for select(recalculate=False)
+# (meyer-lab/vsparse#50), needed by rank_selection's BiCV evaluation to stay a
+# lazy view over a large test-cell block instead of eagerly materializing it.
+# Drop this once #50 merges and a release picks it up, reverting to the plain
+# PyPI version constraint above.
+[tool.uv.sources]
+vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "select-without-recalculate" }
+
 [tool.uv.build-backend]
 module-name = ["scrise"]
 module-root = ""
diff --git a/scrise/rank_selection.py b/scrise/rank_selection.py
index feb9a4cf..a53b5a2a 100644
--- a/scrise/rank_selection.py
+++ b/scrise/rank_selection.py
@@ -105,6 +105,31 @@ def _holdout_scale(
 # temporaries to a few hundred MB regardless of how large the matrix is.
 _MOMENT_CHUNK_NNZ = 50_000_000
 
+# Row-chunk size (in bytes of the dense block materialized per chunk) when
+# streaming a duck-typed backend (e.g. a vsparse normalized view) that has no
+# CSR internals to stream directly, but does support a lazy, stats-preserving
+# `select()`. Keeps a chunk's dense materialization bounded regardless of how
+# large the selected block is.
+_MOMENT_CHUNK_BUDGET_BYTES = 64 << 20
+
+
+def _restrict_rows(X_mat: Any, mask: np.ndarray) -> Any:
+    """Row-restrict ``X_mat`` by a boolean mask, staying lazy where possible.
+
+    For a duck-typed backend that supports ``select()`` (e.g. a vsparse
+    normalized view), keeps the view's *existing* statistics fixed rather
+    than renormalizing the selected rows on their own -- correct here, since
+    the caller is evaluating a model fit against a slice, not treating that
+    slice as its own dataset -- and, unlike bracket indexing on such a view,
+    never eagerly materializes the (potentially huge) selection as a dense
+    array. Plain dense/scipy-sparse ``X_mat`` falls back to ordinary
+    indexing, which is already cheap for those.
+    """
+    select = getattr(X_mat, "select", None)
+    if select is not None:
+        return select(mask, recalculate=False)
+    return X_mat[mask]
+
 
 def _test_block_moments(
     X_mat: Any, cell_mask: np.ndarray, gene_idx: np.ndarray
@@ -116,8 +141,8 @@ def _test_block_moments(
     sums = np.zeros(n_genes)
     squares = np.zeros(n_genes)
 
-    if not sps.issparse(X_mat):
-        block = np.asarray(X_mat)[cell_mask]
+    if isinstance(X_mat, np.ndarray):
+        block = X_mat[cell_mask]
         # Both moments accumulate in float64. A float32 block summed in its own
         # dtype carries ~1e-4 relative error over a few hundred thousand rows,
         # which would reach the reported R2X through `ss_tot`.
@@ -126,6 +151,25 @@ def _test_block_moments(
             np.sum(block.astype(np.float64) ** 2, axis=0)[gene_idx],
         )
 
+    if not sps.issparse(X_mat):
+        # A duck-typed backend (e.g. a vsparse normalized view): select the
+        # cell subset lazily (keeping the view's existing statistics fixed,
+        # rather than renormalizing this subset on its own -- see
+        # `_restrict_rows`) and stream it in row chunks, rather than ever
+        # materializing the whole (potentially huge) subset as one dense
+        # block.
+        sub = X_mat.select(cell_mask, recalculate=False)
+        n_sub_rows = sub.shape[0]
+        chunk_rows = max(1, _MOMENT_CHUNK_BUDGET_BYTES // (n_genes * 8))
+        start = 0
+        while start < n_sub_rows:
+            stop = min(n_sub_rows, start + chunk_rows)
+            block = np.asarray(sub[start:stop], dtype=np.float64)
+            sums += block.sum(axis=0)
+            squares += np.sum(block**2, axis=0)
+            start = stop
+        return sums[gene_idx], squares[gene_idx]
+
     mat = X_mat.tocsr() if X_mat.format != "csr" else X_mat
     n_rows = mat.shape[0]
     start = 0
@@ -229,7 +273,8 @@ def _bicv_trial(
     cond_train = cond_idx[train_cell_mask]
     Z = _cell_loadings(P_train, B, A, cond_train, n_cond)
     ZtY = np.asarray(
-        rmatmul(np.ascontiguousarray(Z.T), X_mat[train_cell_mask]), dtype=np.float64
+        rmatmul(np.ascontiguousarray(Z.T), _restrict_rows(X_mat, train_cell_mask)),
+        dtype=np.float64,
     )[:, test_gene_idx]
     ZtY -= np.outer(Z.sum(axis=0), means_test_genes)
     # `lstsq` on the rank x rank system keeps the minimum-norm
@@ -240,7 +285,8 @@ def _bicv_trial(
     C_full = np.zeros((n_genes, C.shape[1]))
     C_full[train_gene_mask] = C
     cond_test = cond_idx[test_cell_mask]
-    W_test = calc_W(X_mat[test_cell_mask], means, C_full)
+    X_test = _restrict_rows(X_mat, test_cell_mask)
+    W_test = calc_W(X_test, means, C_full)
     cond_slices_test = condition_slices(cond_test, n_cond)
     P_test, _ = project_data(W_test, [A, B, C], cond_slices_test)
 
@@ -248,9 +294,9 @@ def _bicv_trial(
     # held-out loadings need rescaling for the held-out slice's size.
     L = _cell_loadings(P_test, B, A, cond_test, n_cond)
     L *= _holdout_scale(cond_train, cond_test, n_cond)
-    LtY = np.asarray(
-        rmatmul(np.ascontiguousarray(L.T), X_mat[test_cell_mask]), dtype=np.float64
-    )[:, test_gene_idx]
+    LtY = np.asarray(rmatmul(np.ascontiguousarray(L.T), X_test), dtype=np.float64)[
+        :, test_gene_idx
+    ]
     LtY -= np.outer(L.sum(axis=0), means_test_genes)
 
     col_sums, col_squares = _test_block_moments(X_mat, test_cell_mask, test_gene_idx)
diff --git a/scrise/tests/test_rank_selection.py b/scrise/tests/test_rank_selection.py
index 3155ea59..16da0e25 100644
--- a/scrise/tests/test_rank_selection.py
+++ b/scrise/tests/test_rank_selection.py
@@ -507,6 +507,139 @@ def test_block_moments_stream_across_more_than_one_row_block(monkeypatch):
     np.testing.assert_allclose(squares, np.sum(block**2, axis=0), rtol=1e-12)
 
 
+# -- duck-typed backend: a vsparse normalized view ---------------------------
+#
+# `X.X` for a real vsparse-backed dataset (e.g. BAL-Pf2's lazy-normalized-view
+# AnnData) is neither a plain ndarray nor scipy-sparse -- it's a
+# VCSRArrayNormalized/VCSCArrayNormalized. Bracket indexing it with a large
+# boolean mask (the old `X_mat[cell_mask]`) eagerly materializes the whole
+# selection as dense, which is fine at test scale but would be tens of GB at
+# BAL-Pf2's real scale. `_restrict_rows`/`_test_block_moments`'s duck-typed
+# branch exist to avoid that; these tests exercise them directly.
+
+
+def _normalized_view(rng, shape=(120, 30), vsparse_cls=None):
+    import vsparse
+
+    vsparse_cls = vsparse_cls or vsparse.VCSRArray
+    dense = rng.random(shape)
+    dense[dense < 0.6] = 0.0
+    mat = (
+        sps.csr_array(dense)
+        if vsparse_cls is vsparse.VCSRArray
+        else sps.csc_array(dense)
+    )
+    return dense, vsparse_cls.from_scipy(mat).normalized()
+
+
+@pytest.mark.parametrize("fmt", ["csr", "csc"])
+def test_restrict_rows_matches_bracket_indexing_for_a_normalized_view(fmt):
+    import vsparse
+
+    rng = np.random.default_rng(2)
+    vsparse_cls = vsparse.VCSRArray if fmt == "csr" else vsparse.VCSCArray
+    _dense, nv = _normalized_view(rng, vsparse_cls=vsparse_cls)
+    mask = rng.random(nv.shape[0]) < 0.6
+
+    from ..rank_selection import _restrict_rows
+
+    restricted = _restrict_rows(nv, mask)
+    assert isinstance(restricted, type(nv))  # stayed a lazy view, not densified
+    np.testing.assert_allclose(
+        np.asarray(restricted.toarray()), np.asarray(nv[mask, :])
+    )
+
+
+def test_restrict_rows_is_a_no_op_passthrough_for_plain_arrays():
+    from ..rank_selection import _restrict_rows
+
+    rng = np.random.default_rng(3)
+    dense = rng.random((20, 5))
+    mask = rng.random(20) < 0.5
+
+    np.testing.assert_allclose(_restrict_rows(dense, mask), dense[mask])
+    np.testing.assert_allclose(
+        _restrict_rows(sps.csr_array(dense), mask).toarray(), dense[mask]
+    )
+
+
+def test_block_moments_match_a_normalized_view_reference():
+    """`_test_block_moments`'s duck-typed branch matches a dense reference."""
+    import vsparse
+
+    rng = np.random.default_rng(4)
+    dense, nv = _normalized_view(rng, vsparse_cls=vsparse.VCSRArray)
+    mask = rng.random(dense.shape[0]) < 0.7
+    gene_idx = np.sort(rng.choice(dense.shape[1], size=11, replace=False))
+
+    sums, squares = _test_block_moments(nv, mask, gene_idx)
+    block = np.asarray(nv[mask, :])[:, gene_idx]
+    np.testing.assert_allclose(sums, block.sum(axis=0), rtol=1e-10)
+    np.testing.assert_allclose(squares, np.sum(block**2, axis=0), rtol=1e-10)
+
+
+def test_block_moments_stream_a_normalized_view_across_more_than_one_chunk(monkeypatch):
+    """Force several row chunks so the normalized-view streaming path is exercised."""
+    import vsparse
+
+    import scrise.rank_selection as rs
+
+    monkeypatch.setattr(rs, "_MOMENT_CHUNK_BUDGET_BYTES", 64)  # forces 1-row chunks
+    rng = np.random.default_rng(5)
+    dense, nv = _normalized_view(rng, shape=(50, 12), vsparse_cls=vsparse.VCSRArray)
+    mask = rng.random(dense.shape[0]) < 0.8
+    gene_idx = np.arange(12)
+
+    sums, squares = rs._test_block_moments(nv, mask, gene_idx)
+    block = np.asarray(nv[mask, :])
+    np.testing.assert_allclose(sums, block.sum(axis=0), rtol=1e-10)
+    np.testing.assert_allclose(squares, np.sum(block**2, axis=0), rtol=1e-10)
+
+
+def test_block_moments_streaming_peak_is_bounded_by_the_chunk_budget_not_the_block():
+    """The row-chunk streaming loop's own peak tracks the chunk size, not the
+    selected block's total size.
+
+    ``select(recalculate=False)`` on a boolean mask does its own one-time
+    ``O(nnz)`` structural rebuild (the same cost bracket indexing or any
+    other full-array boolean selection pays in vsparse today) -- that part
+    isn't what this streaming code is trying to bound, so it's built once,
+    outside the measurement, isolating the chunk loop itself.
+    """
+    import tracemalloc
+
+    import vsparse
+
+    rng = np.random.default_rng(6)
+    dense, nv = _normalized_view(rng, shape=(4000, 400), vsparse_cls=vsparse.VCSRArray)
+    mask = np.ones(dense.shape[0], dtype=bool)  # select everything
+    sub = nv.select(mask, recalculate=False)
+    n_rows = sub.shape[0]
+
+    def _stream(chunk_rows: int) -> int:
+        tracemalloc.start()
+        try:
+            tracemalloc.reset_peak()
+            start = 0
+            while start < n_rows:
+                stop = min(n_rows, start + chunk_rows)
+                block = np.asarray(sub[start:stop], dtype=np.float64)
+                block.sum(axis=0)
+                start = stop
+            return tracemalloc.get_traced_memory()[1]
+        finally:
+            tracemalloc.stop()
+
+    _stream(50)  # warm up any lazy imports/JIT before measuring
+    small_chunk_peak = _stream(50)
+    large_chunk_peak = _stream(n_rows)  # the whole block in one "chunk"
+
+    # A peak that scaled with the full block regardless of chunk size would
+    # make these roughly equal; bounded streaming keeps the small-chunk peak
+    # well under the whole-block one.
+    assert small_chunk_peak < large_chunk_peak / 4
+
+
 # Held-out slice scaling
 
 
diff --git a/uv.lock b/uv.lock
index bdfe0ae6..cc8cd0ea 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2177,7 +2177,7 @@ requires-dist = [
     { name = "tensorly", specifier = ">=0.9.0" },
     { name = "tlviz", specifier = ">=0.1.1" },
     { name = "tqdm", specifier = ">=4.66.1" },
-    { name = "vsparse", specifier = ">=0.2.0" },
+    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=select-without-recalculate" },
 ]
 provides-extras = ["gpu"]
 
@@ -2471,8 +2471,8 @@ wheels = [
 
 [[package]]
 name = "vsparse"
-version = "0.2.0"
-source = { registry = "https://pypi.org/simple" }
+version = "0.3.0"
+source = { git = "https://github.com/meyer-lab/vsparse?branch=select-without-recalculate#87a9a13c9d0cc84f59608b5595289ba5f9147f78" }
 dependencies = [
     { name = "anndata" },
     { name = "hdf5plugin" },
@@ -2480,10 +2480,6 @@ dependencies = [
     { name = "numpy" },
     { name = "scipy" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/3d/a5/78db8f12cb1dfe2e0a2f7d637efef572ffb75578d01e34f2c1d2a4327747/vsparse-0.2.0.tar.gz", hash = "sha256:aba6221ec50ce3cbe6bc149d69216202c66b2361c9b7a7ba09e3f3bc21a52cff", size = 37836, upload-time = "2026-09-04T12:50:17.51Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/28/14/f2a51a5aa3a3d131429ed7d238d9a13d32dd10493b3019e69e9fd4399420/vsparse-0.2.0-py3-none-any.whl", hash = "sha256:567035afa9a51c5edfc3103a9f98f95568f4ba5e2a9477f73031b703a376d277", size = 45084, upload-time = "2026-09-04T12:50:16.137Z" },
-]
 
 [[package]]
 name = "watchdog"

From 2b5e886683d13d4fa54c54bf38a682550d9a9863 Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Sun, 13 Sep 2026 12:26:03 -0700
Subject: [PATCH 2/6] Point the temporary vsparse pin at the combined testing
 branch

Avoids a conflicting-git-ref resolution error for downstream consumers
(e.g. BAL-Pf2) that pin vsparse's combined bal-pf2-gpu-testing branch
(carrying both #49's matmul kernel and #50's select(recalculate=False))
rather than #50's branch alone.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 pyproject.toml | 15 +++++++++------
 uv.lock        |  4 ++--
 2 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 00d17ecc..e675e0b2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,13 +39,16 @@ fbuild = "analysis.figures.common:genFigure"
 requires = ["uv_build>=0.12.0,<0.13"]
 build-backend = "uv_build"
 
-# TEMPORARY: pin vsparse to an unreleased PR branch for select(recalculate=False)
-# (meyer-lab/vsparse#50), needed by rank_selection's BiCV evaluation to stay a
-# lazy view over a large test-cell block instead of eagerly materializing it.
-# Drop this once #50 merges and a release picks it up, reverting to the plain
-# PyPI version constraint above.
+# TEMPORARY: pin vsparse to an unreleased combined testing branch (carrying
+# both meyer-lab/vsparse#49's matmul kernel and meyer-lab/vsparse#50's
+# select(recalculate=False), which rank_selection's BiCV evaluation needs to
+# stay a lazy view over a large test-cell block instead of eagerly
+# materializing it) rather than #50's own branch directly, so downstream
+# consumers pinning the same combined branch (e.g. BAL-Pf2) don't hit a
+# conflicting-git-ref resolution error. Drop this once #49/#50 merge and a
+# release picks them up, reverting to the plain PyPI version constraint above.
 [tool.uv.sources]
-vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "select-without-recalculate" }
+vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "bal-pf2-gpu-testing" }
 
 [tool.uv.build-backend]
 module-name = ["scrise"]
diff --git a/uv.lock b/uv.lock
index cc8cd0ea..4b126c9e 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2177,7 +2177,7 @@ requires-dist = [
     { name = "tensorly", specifier = ">=0.9.0" },
     { name = "tlviz", specifier = ">=0.1.1" },
     { name = "tqdm", specifier = ">=4.66.1" },
-    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=select-without-recalculate" },
+    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=bal-pf2-gpu-testing" },
 ]
 provides-extras = ["gpu"]
 
@@ -2472,7 +2472,7 @@ wheels = [
 [[package]]
 name = "vsparse"
 version = "0.3.0"
-source = { git = "https://github.com/meyer-lab/vsparse?branch=select-without-recalculate#87a9a13c9d0cc84f59608b5595289ba5f9147f78" }
+source = { git = "https://github.com/meyer-lab/vsparse?branch=bal-pf2-gpu-testing#0beb1af8dfd155a2a21fa61feb824b1c11be65a3" }
 dependencies = [
     { name = "anndata" },
     { name = "hdf5plugin" },

From 460594c72aa6c7b7e181c536c75506d0bd317f92 Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Sun, 13 Sep 2026 13:26:42 -0700
Subject: [PATCH 3/6] Compress a BiCV trial's held-out split once, reused
 across every rank

bicv() previously called run_parafac2 once per (rank, repeat) pair, and
run_parafac2 unconditionally recompressed its input via compress_dataset
whenever compression_kwarg was given -- so a 100-rank x 3-repeat sweep
ran ~300 full CANDELINC compressions of freshly-split raw data, each an
O(nnz) pass, even though CompressedData is explicitly designed to support
fitting any rank <= its own L_g/L_c ("enabling fast rank sweeps without
touching the raw data again" -- parafac2.compress.CompressedData) and
parafac2_nd already accepts a precomputed CompressedData directly.

Each BiCV trial's train/test split is independent of rank (only the
downstream fit is), and the in-sample fit's full dataset doesn't change
across ranks either -- so both were needlessly recompressing the same
data at every rank.

Adds _fit_at_ranks(X_in, ranks, ...), which -- when compression_kwarg is
given -- compresses X_in exactly once via compress_dataset, sized for
max(ranks) via compress_dataset's own rank parameter, then calls
parafac2_nd(compressed, rank=r) directly for every r in ranks, reusing
that one compressed representation. Without compression_kwarg, behavior
is unchanged: each rank still goes through run_parafac2's own per-call
compression shortcut, since there's no CompressedData object to hoist out
in that path.

_bicv_trial now takes `ranks: Sequence[int]` instead of a single `rank`
and returns one result dict per rank, computing the once-per-trial split
and its downstream references (test-block moments, restricted train/test
views) exactly once and reusing them for every rank's held-out scoring,
rather than recomputing per rank. bicv()'s main loop restructures to
match: one held-out split per repeat, evaluated at every rank, instead of
one independent split per (rank, repeat) pair.

New tests confirm compress_dataset is called exactly (1 + n_repeats)
times for an n-rank sweep (not n * (1 + n_repeats)), each sized for the
largest requested rank, and that the compression_kwarg-less path is
unaffected.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014MHUrbvYnPr1naj92oYTab
---
 scrise/rank_selection.py            | 301 +++++++++++++++++++---------
 scrise/tests/test_rank_selection.py |  83 +++++++-
 2 files changed, 284 insertions(+), 100 deletions(-)

diff --git a/scrise/rank_selection.py b/scrise/rank_selection.py
index a53b5a2a..2621d8c4 100644
--- a/scrise/rank_selection.py
+++ b/scrise/rank_selection.py
@@ -18,6 +18,8 @@
 import numpy as np
 import pandas as pd
 import scipy.sparse as sps
+from parafac2.compress import compress_dataset
+from parafac2.parafac2 import parafac2_nd
 from parafac2.utils import calc_W, condition_slices, project_data, rmatmul
 from tqdm import tqdm
 
@@ -195,9 +197,89 @@ def _test_block_moments(
     return sums[gene_idx], squares[gene_idx]
 
 
+def _fit_at_ranks(
+    X_in: anndata.AnnData,
+    ranks: Sequence[int],
+    tolerance: float,
+    max_iter: int,
+    compress: int | tuple[int, int | None] | str | bool | None,
+    compression_kwarg: dict[str, Any] | None,
+    parafac2_kwarg: dict[str, Any] | None,
+    rng: np.random.Generator,
+) -> list[tuple[tuple, float]]:
+    """Fit PARAFAC2 against ``X_in`` at every rank in ``ranks``.
+
+    When ``compression_kwarg`` is given, ``X_in`` is compressed *once*,
+    explicitly, sized for the *largest* rank in ``ranks`` (CANDELINC
+    compression to ``L`` dimensions is valid for fitting any rank ``<= L``,
+    not just the rank it was sized for -- see
+    :func:`parafac2.compress.compress_dataset`/:class:`~parafac2.compress.CompressedData`),
+    and every rank's fit reuses that same compressed representation. Compression
+    is the expensive, ``O(nnz)`` raw-data pass; the fit itself
+    (:func:`~parafac2.parafac2.parafac2_nd` against an already-compressed
+    ``CompressedData``) only touches the small dense compressed cores, so
+    this turns what used to be one full compression *per rank* into one
+    compression for the whole set.
+
+    Without ``compression_kwarg``, each rank's fit instead goes through
+    :func:`~scrise._pf2_utils.run_parafac2` exactly as before -- unchanged,
+    since ``parafac2_nd``'s own internal per-call compression shortcut
+    doesn't hand back a reusable ``CompressedData`` object to hoist out of
+    the loop.
+
+    Returns a list of ``((weights, (A, B, C), projections), r2x)`` tuples,
+    one per rank, in the same order as ``ranks``.
+    """
+    if not compression_kwarg:
+        return [
+            run_parafac2(
+                X_in,
+                rank=rank,
+                random_state=int(rng.integers(np.iinfo(np.int32).max)),
+                tol=tolerance,
+                n_iter_max=max_iter,
+                compress=compress,
+                compression_kwarg=None,
+                parafac2_kwarg=parafac2_kwarg,
+            )
+            for rank in ranks
+        ]
+
+    if not compress:
+        raise ValueError("compression_kwarg requires compress to be set.")
+
+    parafac2_kwarg = dict(parafac2_kwarg) if parafac2_kwarg else {}
+    normalize_slices = parafac2_kwarg.pop("normalize_slices", False)
+    backend = parafac2_kwarg.pop("backend", None)
+
+    compressed = compress_dataset(
+        X_in,
+        L=compress,
+        rank=max(ranks),
+        random_state=int(rng.integers(np.iinfo(np.int32).max)),
+        normalize_slices=normalize_slices,
+        backend=backend,
+        **compression_kwarg,
+    )
+    return [
+        parafac2_nd(
+            compressed,
+            rank=rank,
+            random_state=int(rng.integers(np.iinfo(np.int32).max)),
+            tol=tolerance,
+            n_iter_max=max_iter,
+            normalize_slices=normalize_slices,
+            backend=backend,
+            compress=None,  # already compressed above
+            **parafac2_kwarg,
+        )
+        for rank in ranks
+    ]
+
+
 def _bicv_trial(
     X: anndata.AnnData,
-    rank: int,
+    ranks: Sequence[int],
     held_out_cell_frac: float,
     held_out_gene_frac: float,
     seed: int,
@@ -206,12 +288,14 @@ def _bicv_trial(
     compress: int | tuple[int, int | None] | str | bool | None = "auto",
     compression_kwarg: dict[str, Any] | None = None,
     parafac2_kwarg: dict[str, Any] | None = None,
-) -> dict[str, float]:
-    """Run a single bi-cross-validation trial and return its scores and shape.
+) -> list[dict[str, float]]:
+    """Run a single bi-cross-validation trial, evaluated at every rank in ``ranks``.
 
-    Splits cells (stratified by condition) and genes into train/test blocks,
-    fits PARAFAC2 on the train-cell x train-gene block, then predicts the
-    held-out test-cell x test-gene block:
+    Splits cells (stratified by condition) and genes into train/test blocks
+    *once*, fits PARAFAC2 on the train-cell x train-gene block at each rank
+    (compressing that train block once, not once per rank -- see
+    :func:`_fit_at_ranks`), then predicts the held-out test-cell x test-gene
+    block from each rank's fit:
 
     - Gene loadings for the held-out genes are estimated by regressing the
       train cells' expression of those genes onto the fitted (train-cell)
@@ -222,12 +306,16 @@ def _bicv_trial(
 
     R2X is then computed by reconstructing the held-out block from these
     estimates and comparing against the (mean-centered) observed values.
+    That reference block (and its moments) depends only on the split, not on
+    the rank, so it too is computed once and reused across every rank below.
 
     Returns
     -------
-    dict[str, float]
-        ``BiCV R2X`` (the held-out block), ``Train Block R2X`` (in-sample, on
-        the block the model was fit to), the four block sizes, and ``Seed``.
+    list[dict[str, float]]
+        One dict per rank (in the same order as ``ranks``), each with
+        ``Rank``, ``BiCV R2X`` (the held-out block), ``Train Block R2X``
+        (in-sample, on the block the model was fit to), the four block
+        sizes, and ``Seed``.
     """
     rng = np.random.default_rng(seed)
     cond_idx = X.obs["condition_unique_idxs"].to_numpy().astype(int)
@@ -239,17 +327,16 @@ def _bicv_trial(
     train_gene_mask, test_gene_mask = _split_genes(X.n_vars, held_out_gene_frac, rng)
 
     X_train = X[train_cell_mask][:, train_gene_mask].copy()
-    (weights, (A, B, C), P_train), train_block_r2x = run_parafac2(
+    fits = _fit_at_ranks(
         X_train,
-        rank=rank,
-        random_state=int(rng.integers(np.iinfo(np.int32).max)),
-        tol=tolerance,
-        n_iter_max=max_iter,
-        compress=compress,
-        compression_kwarg=compression_kwarg,
-        parafac2_kwarg=parafac2_kwarg,
+        ranks,
+        tolerance,
+        max_iter,
+        compress,
+        compression_kwarg,
+        parafac2_kwarg,
+        rng,
     )
-    A = A * weights
 
     # Everything below reaches the raw data through products against the raw
     # matrix rather than materialising a block of it. Cells are restricted by
@@ -267,38 +354,17 @@ def _bicv_trial(
     n_genes = X.n_vars
     test_gene_idx = np.flatnonzero(test_gene_mask)
     means_test_genes = means[test_gene_mask]
-
-    # Estimate gene loadings for the held-out genes from the train cells.
-    #   Z^T (X[train, test] - 1 mu^T) = (Z^T X[train])[:, test] - (Z^T 1) mu^T
     cond_train = cond_idx[train_cell_mask]
-    Z = _cell_loadings(P_train, B, A, cond_train, n_cond)
-    ZtY = np.asarray(
-        rmatmul(np.ascontiguousarray(Z.T), _restrict_rows(X_mat, train_cell_mask)),
-        dtype=np.float64,
-    )[:, test_gene_idx]
-    ZtY -= np.outer(Z.sum(axis=0), means_test_genes)
-    # `lstsq` on the rank x rank system keeps the minimum-norm
-    # behaviour when the fit is rank deficient.
-    C_test = np.linalg.lstsq(Z.T @ Z, ZtY, rcond=None)[0].T
-
-    # Estimate projections for the held-out cells from the train genes.
-    C_full = np.zeros((n_genes, C.shape[1]))
-    C_full[train_gene_mask] = C
     cond_test = cond_idx[test_cell_mask]
-    X_test = _restrict_rows(X_mat, test_cell_mask)
-    W_test = calc_W(X_test, means, C_full)
     cond_slices_test = condition_slices(cond_test, n_cond)
-    P_test, _ = project_data(W_test, [A, B, C], cond_slices_test)
-
-    # Score the held-out block. `A` carries the training slice's energy, so the
-    # held-out loadings need rescaling for the held-out slice's size.
-    L = _cell_loadings(P_test, B, A, cond_test, n_cond)
-    L *= _holdout_scale(cond_train, cond_test, n_cond)
-    LtY = np.asarray(rmatmul(np.ascontiguousarray(L.T), X_test), dtype=np.float64)[
-        :, test_gene_idx
-    ]
-    LtY -= np.outer(L.sum(axis=0), means_test_genes)
-
+    scale = _holdout_scale(cond_train, cond_test, n_cond)
+
+    # These depend only on the (fixed, once-per-trial) split, not on rank, so
+    # -- unlike compression/fitting above -- they were already shared across
+    # ranks even before this change; now made explicit and computed once here
+    # rather than once per rank.
+    X_train_lazy = _restrict_rows(X_mat, train_cell_mask)
+    X_test_lazy = _restrict_rows(X_mat, test_cell_mask)
     col_sums, col_squares = _test_block_moments(X_mat, test_cell_mask, test_gene_idx)
     n_test_cells = int(test_cell_mask.sum())
     ss_tot = float(
@@ -308,19 +374,57 @@ def _bicv_trial(
             + n_test_cells * means_test_genes**2
         )
     )
-    cross = float(np.sum(C_test.T * LtY))
-    ss_fit = float(np.sum((L.T @ L) * (C_test.T @ C_test)))
-    ss_res = ss_tot - 2.0 * cross + ss_fit
 
-    return {
-        "BiCV R2X": 1.0 - ss_res / ss_tot,
-        "Train Block R2X": float(train_block_r2x),
-        "NTrainGenes": int(train_gene_mask.sum()),
-        "NTestGenes": int(test_gene_mask.sum()),
-        "NTrainCells": int(train_cell_mask.sum()),
-        "NTestCells": int(test_cell_mask.sum()),
-        "Seed": int(seed),
-    }
+    results = []
+    for rank, ((weights, (A, B, C), P_train), train_block_r2x) in zip(
+        ranks, fits, strict=True
+    ):
+        A = A * weights
+
+        # Estimate gene loadings for the held-out genes from the train cells.
+        #   Z^T (X[train, test] - 1 mu^T) = (Z^T X[train])[:, test] - (Z^T 1) mu^T
+        Z = _cell_loadings(P_train, B, A, cond_train, n_cond)
+        ZtY = np.asarray(
+            rmatmul(np.ascontiguousarray(Z.T), X_train_lazy), dtype=np.float64
+        )[:, test_gene_idx]
+        ZtY -= np.outer(Z.sum(axis=0), means_test_genes)
+        # `lstsq` on the rank x rank system keeps the minimum-norm
+        # behaviour when the fit is rank deficient.
+        C_test = np.linalg.lstsq(Z.T @ Z, ZtY, rcond=None)[0].T
+
+        # Estimate projections for the held-out cells from the train genes.
+        C_full = np.zeros((n_genes, C.shape[1]))
+        C_full[train_gene_mask] = C
+        W_test = calc_W(X_test_lazy, means, C_full)
+        P_test, _ = project_data(W_test, [A, B, C], cond_slices_test)
+
+        # Score the held-out block. `A` carries the training slice's energy,
+        # so the held-out loadings need rescaling for the held-out slice's
+        # size.
+        L = _cell_loadings(P_test, B, A, cond_test, n_cond)
+        L = L * scale
+        LtY = np.asarray(
+            rmatmul(np.ascontiguousarray(L.T), X_test_lazy), dtype=np.float64
+        )[:, test_gene_idx]
+        LtY -= np.outer(L.sum(axis=0), means_test_genes)
+
+        cross = float(np.sum(C_test.T * LtY))
+        ss_fit = float(np.sum((L.T @ L) * (C_test.T @ C_test)))
+        ss_res = ss_tot - 2.0 * cross + ss_fit
+
+        results.append(
+            {
+                "Rank": rank,
+                "BiCV R2X": 1.0 - ss_res / ss_tot,
+                "Train Block R2X": float(train_block_r2x),
+                "NTrainGenes": int(train_gene_mask.sum()),
+                "NTestGenes": int(test_gene_mask.sum()),
+                "NTrainCells": int(train_cell_mask.sum()),
+                "NTestCells": int(test_cell_mask.sum()),
+                "Seed": int(seed),
+            }
+        )
+    return results
 
 
 def _resolve_bicv_inputs(
@@ -394,11 +498,22 @@ def bicv(
 
     For each candidate rank, computes both the ordinary in-sample fit R2X
     (using the full dataset, as in :func:`scrise.factorization.rise_pca_r2x`)
-    and the BiCV R2X (repeated ``n_repeats`` times with independent random
-    cell/gene splits). The fit R2X increases monotonically with rank; the
-    BiCV R2X penalizes overfitting and typically peaks near the rank that
-    best generalizes to held-out data. Plot both with
-    :func:`scrise.plotting.plot_bicv_r2x` to select a rank.
+    and the BiCV R2X (``n_repeats`` independent random cell/gene splits,
+    each evaluated at every rank -- see :func:`_bicv_trial`). The fit R2X
+    increases monotonically with rank; the BiCV R2X penalizes overfitting and
+    typically peaks near the rank that best generalizes to held-out data.
+    Plot both with :func:`scrise.plotting.plot_bicv_r2x` to select a rank.
+
+    Each of the ``n_repeats`` splits (and the full dataset, for the in-sample
+    fit) is compressed once -- sized for the *largest* rank requested -- and
+    every rank's PARAFAC2 fit reuses that same compression, rather than
+    recompressing per rank (see :func:`_fit_at_ranks`). Compression is by far
+    the most expensive, ``O(nnz)`` step against the raw data; fitting a
+    smaller rank from an already-compressed representation only touches the
+    small dense compressed cores. This applies whenever ``compression_kwarg``
+    is given (as it must be to control e.g. ``n_power_iter``); without it,
+    each rank still goes through its own call into ``parafac2_nd``'s internal
+    compression shortcut, unchanged.
 
     Parameters
     ----------
@@ -409,8 +524,9 @@ def bicv(
     ranks : sequence of int
         Candidate rank values to evaluate (e.g., [5, 10, 15, 20, 25, 30]).
     n_repeats : int, optional (default: 3)
-        Number of independent random cell/gene splits per rank. Higher
-        values give a less noisy BiCV estimate but take longer.
+        Number of independent random cell/gene splits, each evaluated at
+        every rank in ``ranks``. Higher values give a less noisy BiCV
+        estimate but take longer.
     held_out_cell_frac : float, optional (default: 0.5)
         Fraction of cells held out per condition in each BiCV trial.
     held_out_gene_frac : float, optional (default: 0.5)
@@ -470,35 +586,38 @@ def bicv(
 
     rng = np.random.default_rng(random_state)
     rows = []
-    for rank in tqdm(ranks, desc="BiCV rank selection"):
-        _, fit_r2x = run_parafac2(
-            X,
-            rank=rank,
-            random_state=int(rng.integers(np.iinfo(np.int32).max)),
-            tol=tolerance,
-            n_iter_max=max_iter,
-            compress=compress,
-            compression_kwarg=compression_kwarg,
-            parafac2_kwarg=parafac2_kwarg,
-        )
+
+    # The full (unsplit) dataset is the same for every rank, so -- like the
+    # per-trial train blocks below -- it's compressed once (sized for the
+    # largest rank) and reused, rather than once per rank.
+    fit_results = _fit_at_ranks(
+        X, ranks, tolerance, max_iter, compress, compression_kwarg, parafac2_kwarg, rng
+    )
+    for rank, (_, fit_r2x) in zip(ranks, fit_results, strict=True):
         # The full-data fit has no split, so the per-trial columns are absent
         # here rather than zero.
         rows.append({"Rank": rank, "Repeat": 0, "Metric": "Fit R2X", "R2X": fit_r2x})
 
-        for repeat in range(n_repeats):
-            trial_seed = int(rng.integers(np.iinfo(np.int32).max))
-            trial = _bicv_trial(
-                X,
-                rank,
-                held_out_cell_frac,
-                held_out_gene_frac,
-                trial_seed,
-                tolerance,
-                max_iter,
-                compress,
-                compression_kwarg,
-                parafac2_kwarg,
-            )
+    # One held-out split per repeat, evaluated at every rank (see
+    # `_bicv_trial`), rather than one independent split per (rank, repeat)
+    # pair -- the split, and the compression of its train block, no longer
+    # need to be redone for each rank.
+    for repeat in tqdm(range(n_repeats), desc="BiCV repeats"):
+        trial_seed = int(rng.integers(np.iinfo(np.int32).max))
+        trial_rows = _bicv_trial(
+            X,
+            ranks,
+            held_out_cell_frac,
+            held_out_gene_frac,
+            trial_seed,
+            tolerance,
+            max_iter,
+            compress,
+            compression_kwarg,
+            parafac2_kwarg,
+        )
+        for trial in trial_rows:
+            rank = trial.pop("Rank")
             rows.append(
                 {
                     "Rank": rank,
diff --git a/scrise/tests/test_rank_selection.py b/scrise/tests/test_rank_selection.py
index 16da0e25..6e1a0474 100644
--- a/scrise/tests/test_rank_selection.py
+++ b/scrise/tests/test_rank_selection.py
@@ -157,6 +157,71 @@ def test_bicv_parafac2_kwarg_and_compression_kwarg():
         )
 
 
+def test_bicv_compresses_once_per_trial_not_once_per_rank(monkeypatch):
+    """The whole point of the redesign: one compression per held-out split
+    (and one for the in-sample fit), reused across every rank -- not one
+    compression per (rank, repeat) pair as before.
+    """
+    import scrise.rank_selection as rs
+
+    X = _make_test_data()
+    ranks = [2, 4, 6]
+    n_repeats = 2
+
+    calls: list[int] = []
+    real_compress_dataset = rs.compress_dataset
+
+    def _counting_compress_dataset(*args, **kwargs):
+        calls.append(kwargs["rank"])
+        return real_compress_dataset(*args, **kwargs)
+
+    monkeypatch.setattr(rs, "compress_dataset", _counting_compress_dataset)
+
+    bicv(
+        X,
+        ranks,
+        n_repeats=n_repeats,
+        random_state=0,
+        max_iter=20,
+        compress="auto",
+        compression_kwarg={"n_power_iter": 1},
+    )
+
+    # One compression for the in-sample fit, plus one per repeat -- not
+    # one per (rank, repeat) pair (which would be len(ranks) * (1 + n_repeats)).
+    assert len(calls) == 1 + n_repeats
+    # Every compression is sized for the *largest* requested rank.
+    assert calls == [max(ranks)] * len(calls)
+
+
+def test_bicv_without_compression_kwarg_still_compresses_per_rank(monkeypatch):
+    """Without compression_kwarg there's no explicit compress_dataset call to
+    hoist -- parafac2_nd's own internal shortcut still runs once per rank,
+    exactly as before this change.
+    """
+    import scrise.rank_selection as rs
+
+    X = _make_test_data()
+    ranks = [2, 4]
+    n_repeats = 1
+
+    calls: list[int] = []
+    real_run_parafac2 = rs.run_parafac2
+
+    def _counting_run_parafac2(*args, **kwargs):
+        calls.append(kwargs["rank"])
+        return real_run_parafac2(*args, **kwargs)
+
+    monkeypatch.setattr(rs, "run_parafac2", _counting_run_parafac2)
+
+    bicv(X, ranks, n_repeats=n_repeats, random_state=0, max_iter=20, compress="auto")
+
+    # One run_parafac2 call per rank for the in-sample fit, plus one per
+    # (rank, repeat) pair for BiCV trials -- unchanged from before.
+    assert len(calls) == len(ranks) + len(ranks) * n_repeats
+    assert sorted(calls) == sorted(ranks * (1 + n_repeats))
+
+
 def test_bicv_adata_alias():
     X = _make_test_data()
     results = bicv(adata=X, ranks=[2], n_repeats=1, random_state=0, max_iter=10)
@@ -325,13 +390,13 @@ def test_reported_seed_reproduces_its_own_trial():
     target = rows.iloc[1]
     replayed = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.5,
         held_out_gene_frac=0.5,
         seed=int(target["Seed"]),
         tolerance=1e-6,
         max_iter=60,
-    )
+    )[0]
     assert replayed["BiCV R2X"] == pytest.approx(target["R2X"], rel=1e-9)
     assert replayed["Train Block R2X"] == pytest.approx(
         target["Train Block R2X"], rel=1e-9
@@ -415,13 +480,13 @@ def test_streamed_scoring_matches_the_dense_formulation(sparse, seed):
 
     got = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.2,
         held_out_gene_frac=0.2,
         seed=seed,
         tolerance=1e-6,
         max_iter=60,
-    )["BiCV R2X"]
+    )[0]["BiCV R2X"]
     want = _dense_reference_trial(X, rank=3, seed=seed)
     assert got == pytest.approx(want, rel=1e-6, abs=1e-9)
 
@@ -436,13 +501,13 @@ def test_streamed_scoring_matches_on_interleaved_conditions():
 
     got = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.2,
         held_out_gene_frac=0.2,
         seed=1,
         tolerance=1e-6,
         max_iter=60,
-    )["BiCV R2X"]
+    )[0]["BiCV R2X"]
     want = _dense_reference_trial(X, rank=3, seed=1)
     assert got == pytest.approx(want, rel=1e-6, abs=1e-9)
 
@@ -705,7 +770,7 @@ def test_exactly_low_rank_data_scores_near_one_at_its_own_rank():
     held-out score still said it was worse than predicting the mean.
     """
     adata = _exact_low_rank(rank_true=4)
-    trial = _bicv_trial(adata, 4, 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=300)
+    trial = _bicv_trial(adata, [4], 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=300)[0]
 
     assert trial["Train Block R2X"] > 0.99
     assert trial["BiCV R2X"] > 0.95
@@ -719,7 +784,7 @@ def test_score_barely_moves_with_the_held_out_fraction(frac):
     same data at the same rank, because the error was sqrt(n_train/n_test).
     """
     adata = _exact_low_rank(rank_true=4)
-    trial = _bicv_trial(adata, 4, frac, frac, seed=0, tolerance=1e-8, max_iter=300)
+    trial = _bicv_trial(adata, [4], frac, frac, seed=0, tolerance=1e-8, max_iter=300)[0]
     assert trial["BiCV R2X"] > 0.9
 
 
@@ -743,7 +808,7 @@ def test_pure_noise_does_not_score_positive():
     adata.obs["condition_unique_idxs"] = cond
     adata.var["means"] = np.zeros(60)
 
-    trial = _bicv_trial(adata, 5, 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=200)
+    trial = _bicv_trial(adata, [5], 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=200)[0]
     assert trial["BiCV R2X"] < 0.02
 
 

From 6a3794924a9f18e9259baf7289d0a322810dd47e Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Mon, 14 Sep 2026 08:51:59 -0700
Subject: [PATCH 4/6] Point the vsparse pin at the new to-scipy-sparse-dtype
 branch

vsparse#49/#50 (the fixes this pin originally existed for) have merged to
vsparse's main. Repointing to meyer-lab/vsparse@to-scipy-sparse-dtype
(main plus vsparse#51's still-open dtype argument for to_scipy_sparse,
which RISE doesn't itself need) rather than main directly, so downstream
consumers pinning the same branch for #51 (e.g. BAL-Pf2) don't hit a
conflicting-git-ref resolution error against RISE's own vsparse pin.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014MHUrbvYnPr1naj92oYTab
---
 pyproject.toml | 18 ++++++++++--------
 uv.lock        |  4 ++--
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index e675e0b2..65b2cee1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,16 +39,18 @@ fbuild = "analysis.figures.common:genFigure"
 requires = ["uv_build>=0.12.0,<0.13"]
 build-backend = "uv_build"
 
-# TEMPORARY: pin vsparse to an unreleased combined testing branch (carrying
-# both meyer-lab/vsparse#49's matmul kernel and meyer-lab/vsparse#50's
-# select(recalculate=False), which rank_selection's BiCV evaluation needs to
+# TEMPORARY: meyer-lab/vsparse#49's matmul kernel and #50's
+# select(recalculate=False) (which rank_selection's BiCV evaluation needs to
 # stay a lazy view over a large test-cell block instead of eagerly
-# materializing it) rather than #50's own branch directly, so downstream
-# consumers pinning the same combined branch (e.g. BAL-Pf2) don't hit a
-# conflicting-git-ref resolution error. Drop this once #49/#50 merge and a
-# release picks them up, reverting to the plain PyPI version constraint above.
+# materializing it) have both merged to vsparse's main. Pinned to
+# to-scipy-sparse-dtype (main plus vsparse#51's still-open dtype argument
+# for to_scipy_sparse, which RISE itself doesn't need) rather than main
+# directly, so downstream consumers pinning the same branch (e.g. BAL-Pf2,
+# which does need #51) don't hit a conflicting-git-ref resolution error.
+# Drop this once #51 merges and a release picks it up, reverting to the
+# plain PyPI version constraint above.
 [tool.uv.sources]
-vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "bal-pf2-gpu-testing" }
+vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "to-scipy-sparse-dtype" }
 
 [tool.uv.build-backend]
 module-name = ["scrise"]
diff --git a/uv.lock b/uv.lock
index 4b126c9e..e34984ba 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2177,7 +2177,7 @@ requires-dist = [
     { name = "tensorly", specifier = ">=0.9.0" },
     { name = "tlviz", specifier = ">=0.1.1" },
     { name = "tqdm", specifier = ">=4.66.1" },
-    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=bal-pf2-gpu-testing" },
+    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=to-scipy-sparse-dtype" },
 ]
 provides-extras = ["gpu"]
 
@@ -2472,7 +2472,7 @@ wheels = [
 [[package]]
 name = "vsparse"
 version = "0.3.0"
-source = { git = "https://github.com/meyer-lab/vsparse?branch=bal-pf2-gpu-testing#0beb1af8dfd155a2a21fa61feb824b1c11be65a3" }
+source = { git = "https://github.com/meyer-lab/vsparse?branch=to-scipy-sparse-dtype#671273895466be14474b9cc64b87ef0ee650c573" }
 dependencies = [
     { name = "anndata" },
     { name = "hdf5plugin" },

From 58d427c06caf80154fc969557e962f9bdd5a27b6 Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Mon, 14 Sep 2026 10:33:25 -0700
Subject: [PATCH 5/6] Pin parafac2 to at least the latest PyPI release (1.7.0)

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 pyproject.toml |  2 +-
 uv.lock        | 90 +++++++++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 87 insertions(+), 5 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index afbdc101..4a21ebd1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -13,7 +13,7 @@ dependencies = [
     "pandas>=3.0.0",
     "pyarrow>=19.0",
     "tensorly>=0.9.0",
-    "parafac2>=1.6.0",
+    "parafac2>=1.7.0",
     "anndata>=0.13",
     "pacmap>=0.9",
     "tqdm>=4.66.1",
diff --git a/uv.lock b/uv.lock
index e34984ba..56d445b1 100644
--- a/uv.lock
+++ b/uv.lock
@@ -271,6 +271,48 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/ec/82/32e3bd191d498e64f6f911ad55d14006a0861e54869d2d32452326399e65/coverage-7.15.2-py3-none-any.whl", hash = "sha256:eb6bcae8d1a9d305351ecb108232441d11c5cfe9de840a04388ba5d2db8d735c", size = 213375, upload-time = "2026-07-15T18:56:17.305Z" },
 ]
 
+[[package]]
+name = "cuda-bindings"
+version = "13.4.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-pathfinder" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/38/cf/165b4d449f94956c2a60930cf5dfeb27132ead60a7e7f2c37819df1cba07/cuda_bindings-13.4.1-cp313-cp313-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:b601c0cbf0dffb648f68e56a60b320738a20210293f33896a1964a6438cc65f1", size = 6313772, upload-time = "2026-09-10T01:17:03.969Z" },
+    { url = "https://files.pythonhosted.org/packages/d9/f9/cf021d1560541caa1f35f3e7e311d2678dbacb4fa6a4573b63470fe1ae00/cuda_bindings-13.4.1-cp313-cp313-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0e2c357698588b06ebd65811ee2013b6650dd0a10d16d899924aabad0d606d76", size = 6924300, upload-time = "2026-09-10T01:17:06.23Z" },
+    { url = "https://files.pythonhosted.org/packages/66/93/d41e182cde0f1dd3c6c6aac018cb550ffedf55ba2a67c427320e91048592/cuda_bindings-13.4.1-cp313-cp313-win_amd64.whl", hash = "sha256:11bce22fe70a19eb6e727ae7126940245160e0223b3f1cb3e209a313449f75ea", size = 6289385, upload-time = "2026-09-10T01:17:08.372Z" },
+    { url = "https://files.pythonhosted.org/packages/d0/cd/e2100b3d5e2b2e023a603a6463df4582c73d057daed9e7cf6802cd3b8455/cuda_bindings-13.4.1-cp313-cp313-win_arm64.whl", hash = "sha256:73d86523ee2f749cf312fdada7be83160c84ddf4845232769f044695ec1c0e68", size = 5293378, upload-time = "2026-09-10T01:17:10.356Z" },
+    { url = "https://files.pythonhosted.org/packages/58/17/74346b49114779920929ec0ea1361f0a0357262d6020cfdd017f670da638/cuda_bindings-13.4.1-cp314-cp314-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:d7df2dddb81feb15787e8c4a13b7aa3f4c23eabeec7746037e94dbbed065fc6c", size = 6406821, upload-time = "2026-09-10T01:17:12.344Z" },
+    { url = "https://files.pythonhosted.org/packages/f6/77/2f9a38be7399a34e3703b6ce1be60bb2c56611d324750f8d544ff90b0473/cuda_bindings-13.4.1-cp314-cp314-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:62df23df11074e9833bf348bbcf0b8eec2fcbded4f305c6fbaa3ed067433e97d", size = 6972525, upload-time = "2026-09-10T01:17:15.098Z" },
+    { url = "https://files.pythonhosted.org/packages/fb/88/28195fa9f0c7aebf9978b5508e463289091a574240c4c50f9035b1d96764/cuda_bindings-13.4.1-cp314-cp314-win_amd64.whl", hash = "sha256:fb66fa7ea3f257a0451d54edc4b5a4a052eb29c32ba15b104f465ecb24907ffb", size = 6434166, upload-time = "2026-09-10T01:17:17.097Z" },
+    { url = "https://files.pythonhosted.org/packages/2d/f5/41388671cbe81477aa30fbccbf195ca6f207eeaca5b1268282f1aba827cd/cuda_bindings-13.4.1-cp314-cp314-win_arm64.whl", hash = "sha256:506457cc214f3dac300044431a36b3f608f7a8c495c5bb40b9008d162ff1d5d3", size = 5450466, upload-time = "2026-09-10T01:17:19.099Z" },
+    { url = "https://files.pythonhosted.org/packages/7d/de/41197aebf91f6c5f82b35e06e3b4bbe08edddb335c4d6e53c1f1fe542e0f/cuda_bindings-13.4.1-cp314-cp314t-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:0947a6c491a622b076e9afcbc0444eefd9cbe001630558bbcba8d0a90dab7fe7", size = 6243282, upload-time = "2026-09-10T01:17:21.158Z" },
+    { url = "https://files.pythonhosted.org/packages/be/d6/1b697092f53cfd4d721fbbe13e7d66dc72a9d508156af9be4c86607e4373/cuda_bindings-13.4.1-cp314-cp314t-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:51e730e80d997b6037566033a79aae3995e1981654054536f512e5c091679c70", size = 6803245, upload-time = "2026-09-10T01:17:23.309Z" },
+    { url = "https://files.pythonhosted.org/packages/f4/19/c4e0cb61555f553bc01a533f558b2ba775c88b0df5ac3dfba70125948712/cuda_bindings-13.4.1-cp314-cp314t-win_amd64.whl", hash = "sha256:2a95d8e231d219c20beed1c3c97f58603375a3f2d1883d076a6584190606c9b0", size = 7331806, upload-time = "2026-09-10T01:17:25.641Z" },
+    { url = "https://files.pythonhosted.org/packages/15/dd/80fb0ea15ed15ac1b2b27d2c1d0eef4e725f37effbc4946ddd0795cfc8e3/cuda_bindings-13.4.1-cp314-cp314t-win_arm64.whl", hash = "sha256:76c10f8ad5cf55676518b2cb3d1e3d867826389b29d8418b402d23ba4df7d655", size = 5421917, upload-time = "2026-09-10T01:17:27.82Z" },
+]
+
+[[package]]
+name = "cuda-core"
+version = "1.2.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-pathfinder" },
+    { name = "numpy" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/c5/88/145fdb860aa0cafbe197af627567f2f5347d33e61651b38a7361487eb64c/cuda_core-1.2.0-cp313-cp313-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:2a77d40f613e4d151f4466818d1baab271237122c7240450fe66fe8a071ff86f", size = 6154425, upload-time = "2026-09-03T00:12:36.418Z" },
+    { url = "https://files.pythonhosted.org/packages/aa/bf/ccfe842dc59003a721a37103c1964bb797b4b3362d07f82b6bbc71926aec/cuda_core-1.2.0-cp313-cp313-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0dddba476d990c62a17b58b2c198a3b2768d575987577febda8f3e88140df192", size = 6562184, upload-time = "2026-09-03T00:12:38.441Z" },
+    { url = "https://files.pythonhosted.org/packages/9f/fd/1ebc50d239f61211f8732a7b6e8a35d363e2828b5eeca1930dc3a89925f0/cuda_core-1.2.0-cp313-cp313-win_amd64.whl", hash = "sha256:f19593699a49a52d08300d44c1c84db8629dc01702b36336401c0ba87fe5067f", size = 5923555, upload-time = "2026-09-03T00:12:40.511Z" },
+    { url = "https://files.pythonhosted.org/packages/e0/7c/fd915c7f8cdb954744c1fda75fefa75a39e73508b6aae8ec7369800c8491/cuda_core-1.2.0-cp314-cp314-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:4206c19f3c787c6dd4fdce3625e842686dc73ba665d46d16d6c57659a94db0ac", size = 6233228, upload-time = "2026-09-03T00:12:42.122Z" },
+    { url = "https://files.pythonhosted.org/packages/76/69/2e255b5db595188ba12911fcd9efa45ccf71b02cd96809bc950923dc76d4/cuda_core-1.2.0-cp314-cp314-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:e138be12af795c69c1cb61e38562ea70e5a38616b2753e2a265231f1941394ef", size = 6595561, upload-time = "2026-09-03T00:12:44.001Z" },
+    { url = "https://files.pythonhosted.org/packages/16/c5/7d8a1d90071f5aac2a3bf11f1829535a60519a673809b5c3bc3cdb30df54/cuda_core-1.2.0-cp314-cp314-win_amd64.whl", hash = "sha256:419e884cbe26d85936c8ffea27865a6f2589677f34a81eecb0fe1ec3715c06a0", size = 6051744, upload-time = "2026-09-03T00:12:45.754Z" },
+    { url = "https://files.pythonhosted.org/packages/e7/f3/22b67d9b035fb22d055a39ca21affc4246008da7e72aa024a8fd771b9990/cuda_core-1.2.0-cp314-cp314t-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:c4228d27af1993f8516b51db6c74846ebd1f772bb34132ee721da61369b7febd", size = 6408313, upload-time = "2026-09-03T00:12:48.128Z" },
+    { url = "https://files.pythonhosted.org/packages/a1/82/5edfcbff353cb4dedd13415e5793ca710b302a26f9ee9144d27e9a680bc1/cuda_core-1.2.0-cp314-cp314t-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:ccfd9905e43dc5ef0f88e26a3e3f1509d84c0306710dec351c5ab1db0da5752e", size = 6683167, upload-time = "2026-09-03T00:12:50.178Z" },
+    { url = "https://files.pythonhosted.org/packages/38/eb/30258a7fe759902f1603b960aaedd843666179e686d6ad47d04ba9914a1c/cuda_core-1.2.0-cp314-cp314t-win_amd64.whl", hash = "sha256:29ec78a8a45ad029e8d31b0ee67fe658964efc97a0176eb9d7ad354677116923", size = 6972170, upload-time = "2026-09-03T00:12:52.247Z" },
+]
+
 [[package]]
 name = "cuda-pathfinder"
 version = "1.6.0"
@@ -1345,6 +1387,29 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/dd/7e/2eecb277d8a98184d881fb98a738363fd4f14577a4d2d7f8264266e82623/nvidia_nvjitlink_cu12-12.9.86-py3-none-win_amd64.whl", hash = "sha256:cc6fcec260ca843c10e34c936921a1c426b351753587fdd638e8cff7b16bb9db", size = 35584936, upload-time = "2025-06-05T20:16:08.525Z" },
 ]
 
+[[package]]
+name = "nvmath-python"
+version = "1.0.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-bindings" },
+    { name = "cuda-core" },
+    { name = "cuda-pathfinder" },
+    { name = "numpy" },
+    { name = "pywin32", marker = "sys_platform == 'win32'" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/3b/43/1bf0b6d6e97ad3ea5040a026c5aff03f1c60fceb977bb943c2a6bbc25294/nvmath_python-1.0.0-cp313-cp313-manylinux_2_28_aarch64.whl", hash = "sha256:f550e28e87c253ab8747c559c1c4faeb15b79f7c2684bb23e8d4c42df4f5852d", size = 5056945, upload-time = "2026-07-09T20:16:00.903Z" },
+    { url = "https://files.pythonhosted.org/packages/83/f2/7ac589040aeaa92bf72506b480f5a6b9065cb160e10513db07420f0a0eff/nvmath_python-1.0.0-cp313-cp313-manylinux_2_28_x86_64.whl", hash = "sha256:f1e1bc567e56e7ff72515ac34c5af328478e7b8226c6f92b6b2f324925bb38a8", size = 5379834, upload-time = "2026-07-09T20:16:24.375Z" },
+    { url = "https://files.pythonhosted.org/packages/53/52/1cd2fb1fdc2ce3c29341a192ca47209593c7a4e8e1f79f2512c2980ed03f/nvmath_python-1.0.0-cp313-cp313-win_amd64.whl", hash = "sha256:d3a4085a9f98faa223fbb13d6244f8db4d61d48f3f5785e84f6c115b6ee7e8a4", size = 4087396, upload-time = "2026-07-09T20:16:47.375Z" },
+    { url = "https://files.pythonhosted.org/packages/7a/a2/512f5d0a9e2bdfe0b52c7afe7dbdd3c140fa18ab41c1c863f59eee8808df/nvmath_python-1.0.0-cp314-cp314-manylinux_2_28_aarch64.whl", hash = "sha256:7f78fb221df17fbb8e1303d536f5b945122c10febd74b37dea9cb8e44cc586a3", size = 5094362, upload-time = "2026-07-09T20:17:09.604Z" },
+    { url = "https://files.pythonhosted.org/packages/4c/74/a07afd16866ee9f3afcf60e7d5dd1ad3ddc5fb9d56998fce0629eb2671a9/nvmath_python-1.0.0-cp314-cp314-manylinux_2_28_x86_64.whl", hash = "sha256:b6ebb9ae1a2d16dc0b7431d0e6c3c72aa7b11f3e3e83d4dc4f046b086fff8fb8", size = 5395750, upload-time = "2026-07-09T20:17:33.616Z" },
+    { url = "https://files.pythonhosted.org/packages/43/da/d94c5904b0c7878be75994ba541b7937c9ecef50dc633fa10d149885cf4f/nvmath_python-1.0.0-cp314-cp314-win_amd64.whl", hash = "sha256:89807bd2e8a8ef79bf73868638c83fe5f94082f7d8d7b87e08cc94d072959ec4", size = 4204235, upload-time = "2026-07-09T20:17:57.785Z" },
+    { url = "https://files.pythonhosted.org/packages/85/20/f4627d20fef8a2df6dabe8758717d254535529104e6a17c80078276b3ba3/nvmath_python-1.0.0-cp314-cp314t-manylinux_2_28_aarch64.whl", hash = "sha256:e572183b26b850e9e7fd11abbab2b3cc74ab5e13afb04025049e293fee84b4ad", size = 5009618, upload-time = "2026-07-09T20:18:19.401Z" },
+    { url = "https://files.pythonhosted.org/packages/a2/bc/cc88e62b71c37e1b569680507e9eca5f8ea27a8a31a58db3aa6f8250b9ed/nvmath_python-1.0.0-cp314-cp314t-manylinux_2_28_x86_64.whl", hash = "sha256:96fe75363269113d933c9e158da5ce2d585adfab1b7bfb8fdda1a16a202cba52", size = 5292397, upload-time = "2026-07-09T20:18:39.616Z" },
+    { url = "https://files.pythonhosted.org/packages/35/65/958fb26ab46f11b4e8503792f9dc86ef3a18f705206f71181793ca6b3be1/nvmath_python-1.0.0-cp314-cp314t-win_amd64.whl", hash = "sha256:b6ab77f61fca32a77c50075b78f366c7dc05733f46f744068644d63a9530b611", size = 4818972, upload-time = "2026-07-09T20:19:01.313Z" },
+]
+
 [[package]]
 name = "packaging"
 version = "26.2"
@@ -1417,7 +1482,7 @@ wheels = [
 
 [[package]]
 name = "parafac2"
-version = "1.6.0"
+version = "1.7.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "anndata" },
@@ -1426,15 +1491,16 @@ dependencies = [
     { name = "tensorly" },
     { name = "tqdm" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/5b/11/f69fc560bf17fbdcf269fdb63d9106724ded27dee395e74afb40e061e93f/parafac2-1.6.0.tar.gz", hash = "sha256:b7e7817a28fc59ba1f1d80469c65d74d4a6044bf9f383364eb07cef1f4d3516e", size = 28064, upload-time = "2026-09-02T01:02:14.823Z" }
+sdist = { url = "https://files.pythonhosted.org/packages/40/55/35d771bcd672150a84993ea89764da418d7685df1925c17728e9fc35edfe/parafac2-1.7.0.tar.gz", hash = "sha256:074907a94c64ec0fce4c22cdf73396bb511f7cabd04af353f831302f14e4b1ee", size = 43927, upload-time = "2026-09-14T17:07:34.24Z" }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/06/58/654144ff025927b36895e515b7b01083959afa7ada96909acf49d364a08f/parafac2-1.6.0-py3-none-any.whl", hash = "sha256:dd5b1575bc7c1a779f794ec1c0d2e9c5dc0c32feb8ad6e46e60b0c592b57223a", size = 33469, upload-time = "2026-09-02T01:02:13.391Z" },
+    { url = "https://files.pythonhosted.org/packages/72/8e/634f3f244c9cb7f6b277f1240d09710e4f5fe3c62843242a21239d3c0d74/parafac2-1.7.0-py3-none-any.whl", hash = "sha256:c49ecc4eb4d6e971bce65cbb8ae3f9979ae0450a49b492026bce64633661dd74", size = 53172, upload-time = "2026-09-14T17:07:32.935Z" },
 ]
 
 [package.optional-dependencies]
 gpu = [
     { name = "cupy-cuda12x", extra = ["ctk"], marker = "sys_platform != 'darwin'" },
     { name = "mlx", marker = "platform_machine == 'arm64' and sys_platform == 'darwin'" },
+    { name = "nvmath-python", marker = "sys_platform != 'darwin'" },
 ]
 
 [[package]]
@@ -1828,6 +1894,22 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/0b/d7/1959b9648791274998a9c3526f6d0ec8fd2233e4d4acce81bbae76b44b2a/python_dotenv-1.2.2-py3-none-any.whl", hash = "sha256:1d8214789a24de455a8b8bd8ae6fe3c6b69a5e3d64aa8a8e5d68e694bbcb285a", size = 22101, upload-time = "2026-03-01T16:00:25.09Z" },
 ]
 
+[[package]]
+name = "pywin32"
+version = "312"
+source = { registry = "https://pypi.org/simple" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/2d/41/12fbfd7f36ed2146d8bc9de96c2741296bf0d490b98508496cff322e274c/pywin32-312-cp313-cp313-win32.whl", hash = "sha256:7a27df850933d16a8eabfbaeb73d52b273e2da667f80d70b01a89d1f6828d02c", size = 6370184, upload-time = "2026-06-04T07:49:36.253Z" },
+    { url = "https://files.pythonhosted.org/packages/ba/db/36a78e3403099d31d9746d13fdcde5accc43c1155f375a34d15983a479a7/pywin32-312-cp313-cp313-win_amd64.whl", hash = "sha256:c53e878d15a1c44788082bfe712a905433473aa38f86375b7cf8b45e3acbaaf9", size = 6914298, upload-time = "2026-06-04T07:49:38.876Z" },
+    { url = "https://files.pythonhosted.org/packages/84/37/c1697194092b76de9ed47ca124323f02c57ffc8a45c06f88a3d5acaf01eb/pywin32-312-cp313-cp313-win_arm64.whl", hash = "sha256:59aba5d5940842075343a5ddc6b11f1cdf0d1567fe745290359dfbcc7c2eb831", size = 6727640, upload-time = "2026-06-04T07:49:41.083Z" },
+    { url = "https://files.pythonhosted.org/packages/fc/2b/1f3cded5822fd49c02f40544cbb5f58c7cfd6b1694869fd476cb6170ee97/pywin32-312-cp314-cp314-win32.whl", hash = "sha256:a77a90fbb6881238d2ca9c6fd797b25817f3768fe78d214a90137ff055a75f5b", size = 6468928, upload-time = "2026-06-04T07:49:43.188Z" },
+    { url = "https://files.pythonhosted.org/packages/21/82/3bf86d2e2808902013132e1ce905a7da0da53790f3836c64bf44d55e24f3/pywin32-312-cp314-cp314-win_amd64.whl", hash = "sha256:a4dd3a848290ef724347b19f301045831d8e802fa4464f491b98b1e0a081432e", size = 7024157, upload-time = "2026-06-04T07:49:45.34Z" },
+    { url = "https://files.pythonhosted.org/packages/a4/0e/73f6d6800b4f27655abd9e9f6aaeaefcddb2b946e4674efa2bab184a7f7b/pywin32-312-cp314-cp314-win_arm64.whl", hash = "sha256:9fce94568364e0155e6dfb781ac5d95903be8baf28670632beab1b523f300daa", size = 6839598, upload-time = "2026-06-04T07:49:47.613Z" },
+    { url = "https://files.pythonhosted.org/packages/eb/61/caa39686032d2ebdd04ff0ab5cbe163126c0066d98e00c9018646e42393b/pywin32-312-cp315-cp315-win32.whl", hash = "sha256:5c1fbe4a937a73ae9297384a3da38518cbc694c68ad8a809b2e19acd350f03ed", size = 6471159, upload-time = "2026-06-04T07:49:50.035Z" },
+    { url = "https://files.pythonhosted.org/packages/0f/cd/7e1de64a4a6f69c04214169657ccab0d93a670ea50e35eb8f489d7378249/pywin32-312-cp315-cp315-win_amd64.whl", hash = "sha256:c2f03a0f73f804a13c2735b99392b0cd426bb4f2c4d0178e5ac966a0f21618d5", size = 7025293, upload-time = "2026-06-04T07:49:54.857Z" },
+    { url = "https://files.pythonhosted.org/packages/23/ed/4532e9388e65fa16b46776ef47ad631a64eda1631884488af707666350ed/pywin32-312-cp315-cp315-win_arm64.whl", hash = "sha256:a8597d28f267b39074aef51fa593530082b39cbe5a074226096857b1fed2dfb9", size = 6840337, upload-time = "2026-06-04T07:49:57.531Z" },
+]
+
 [[package]]
 name = "pyyaml"
 version = "6.0.3"
@@ -2168,7 +2250,7 @@ requires-dist = [
     { name = "numpy", specifier = ">=2.2" },
     { name = "pacmap", specifier = ">=0.9" },
     { name = "pandas", specifier = ">=3.0.0" },
-    { name = "parafac2", specifier = ">=1.6.0" },
+    { name = "parafac2", specifier = ">=1.7.0" },
     { name = "parafac2", extras = ["gpu"], marker = "extra == 'gpu'" },
     { name = "pyarrow", specifier = ">=19.0" },
     { name = "scikit-learn", specifier = ">=1.6" },

From 374b6448f0e7750d09faab35fd63b025ae58d708 Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Mon, 14 Sep 2026 10:36:07 -0700
Subject: [PATCH 6/6] Remove the temporary vsparse git pin now that 0.4.0 is on
 PyPI

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 pyproject.toml | 13 -------------
 uv.lock        | 10 +++++++---
 2 files changed, 7 insertions(+), 16 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 4a21ebd1..ea329e0b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,19 +39,6 @@ fbuild = "analysis.figures.common:genFigure"
 requires = ["uv_build>=0.12.0,<0.13"]
 build-backend = "uv_build"
 
-# TEMPORARY: meyer-lab/vsparse#49's matmul kernel and #50's
-# select(recalculate=False) (which rank_selection's BiCV evaluation needs to
-# stay a lazy view over a large test-cell block instead of eagerly
-# materializing it) have both merged to vsparse's main. Pinned to
-# to-scipy-sparse-dtype (main plus vsparse#51's still-open dtype argument
-# for to_scipy_sparse, which RISE itself doesn't need) rather than main
-# directly, so downstream consumers pinning the same branch (e.g. BAL-Pf2,
-# which does need #51) don't hit a conflicting-git-ref resolution error.
-# Drop this once #51 merges and a release picks it up, reverting to the
-# plain PyPI version constraint above.
-[tool.uv.sources]
-vsparse = { git = "https://github.com/meyer-lab/vsparse", branch = "to-scipy-sparse-dtype" }
-
 [tool.uv.build-backend]
 module-name = ["scrise"]
 module-root = ""
diff --git a/uv.lock b/uv.lock
index 56d445b1..c3015e90 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2259,7 +2259,7 @@ requires-dist = [
     { name = "tensorly", specifier = ">=0.9.0" },
     { name = "tlviz", specifier = ">=0.1.1" },
     { name = "tqdm", specifier = ">=4.66.1" },
-    { name = "vsparse", git = "https://github.com/meyer-lab/vsparse?branch=to-scipy-sparse-dtype" },
+    { name = "vsparse", specifier = ">=0.4.0" },
 ]
 provides-extras = ["gpu"]
 
@@ -2553,8 +2553,8 @@ wheels = [
 
 [[package]]
 name = "vsparse"
-version = "0.3.0"
-source = { git = "https://github.com/meyer-lab/vsparse?branch=to-scipy-sparse-dtype#671273895466be14474b9cc64b87ef0ee650c573" }
+version = "0.4.0"
+source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "anndata" },
     { name = "hdf5plugin" },
@@ -2562,6 +2562,10 @@ dependencies = [
     { name = "numpy" },
     { name = "scipy" },
 ]
+sdist = { url = "https://files.pythonhosted.org/packages/90/36/c4e280f14af925b71c223777389bfa3ad472f4540375880bfe9b50ac6225/vsparse-0.4.0.tar.gz", hash = "sha256:72faa01087fcdd3fd72467638303c53838ff226685b0126cf5c060fc54626742", size = 47657, upload-time = "2026-09-14T17:12:39.176Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/ff/4e/3a167ca192bcb10623cbf403e76798d600e249c7c26f0bd28f24d0ed15b6/vsparse-0.4.0-py3-none-any.whl", hash = "sha256:a294b647dd8fac06746a8d8ffd1a778a31d60a6387d50f5915a406a37d768f2c", size = 54961, upload-time = "2026-09-14T17:12:38.011Z" },
+]
 
 [[package]]
 name = "watchdog"