From 3fd8765023200494336cf3a6258837eb4cf3985f Mon Sep 17 00:00:00 2001
From: Aaron Meyer <ameyer@ucla.edu>
Date: Wed, 2 Sep 2026 09:40:48 -0700
Subject: [PATCH 1/2] Add property-based, contract, regression, and plotting
 tests

Fills the gaps identified in the testing-strategy review: Hypothesis-driven
invariant tests for the core factorization/OPQ primitives, explicit
error-path tests for documented AnnData-key/parameter contracts, a golden
regression fixture pinning pf2() output on fixed synthetic data, and smoke
tests for the previously ~0-30%-covered plotting modules. Adds hypothesis
as a dev dependency and a shared synthetic-data conftest, and gates CI on
85% coverage (currently 89%) via tool.coverage.report.fail_under.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 pyproject.toml                          |   5 +
 scrise/tests/conftest.py                | 135 +++++++++++++++
 scrise/tests/test_contracts.py          | 148 ++++++++++++++++
 scrise/tests/test_golden_regression.py  |  75 ++++++++
 scrise/tests/test_invariants.py         | 216 ++++++++++++++++++++++++
 scrise/tests/test_plotting_factors.py   | 114 +++++++++++++
 scrise/tests/test_plotting_general.py   | 148 ++++++++++++++++
 scrise/tests/test_plotting_pacmap.py    |  89 ++++++++++
 scrise/tests/test_plotting_stability.py |  99 +++++++++++
 uv.lock                                 |  67 ++++++++
 10 files changed, 1096 insertions(+)
 create mode 100644 scrise/tests/conftest.py
 create mode 100644 scrise/tests/test_contracts.py
 create mode 100644 scrise/tests/test_golden_regression.py
 create mode 100644 scrise/tests/test_invariants.py
 create mode 100644 scrise/tests/test_plotting_factors.py
 create mode 100644 scrise/tests/test_plotting_general.py
 create mode 100644 scrise/tests/test_plotting_pacmap.py
 create mode 100644 scrise/tests/test_plotting_stability.py

diff --git a/pyproject.toml b/pyproject.toml
index e6f06b3a..e21d6f7d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -68,6 +68,7 @@ dev = [
     { include-group = "analysis" },
     "pytest>=9.0",
     "pytest-cov>=7.0",
+    "hypothesis>=6.140",
     "ty",
     "ruff>=0.16",
 ]
@@ -84,3 +85,7 @@ filterwarnings = [
     "ignore::DeprecationWarning",
     "ignore::PendingDeprecationWarning:seaborn",
 ]
+
+[tool.coverage.report]
+fail_under = 85
+omit = ["scrise/tests/*"]
diff --git a/scrise/tests/conftest.py b/scrise/tests/conftest.py
new file mode 100644
index 00000000..063252ff
--- /dev/null
+++ b/scrise/tests/conftest.py
@@ -0,0 +1,135 @@
+"""
+Shared fixtures and synthetic-data factories for the scrise test suite.
+"""
+
+from collections.abc import Mapping, Sequence
+from typing import Any, cast
+
+import anndata
+import matplotlib
+import numpy as np
+import pandas as pd
+import pytest
+
+matplotlib.use("Agg")
+
+import matplotlib.pyplot as plt  # noqa: E402
+
+
+@pytest.fixture(autouse=True)
+def _close_figures():
+    """Prevent matplotlib figures from accumulating across tests."""
+    yield
+    plt.close("all")
+
+
+def make_synthetic_pf2_data(
+    n_cond: int = 5,
+    n_genes: int = 40,
+    rank: int = 3,
+    seed: int = 0,
+    cells_per_cond: tuple[int, int] = (60, 90),
+) -> anndata.AnnData:
+    """Build a synthetic AnnData with a known low-rank PARAFAC2-like structure.
+
+    Mirrors the generator originally written for ``test_rank_selection.py``
+    so every test module shares one synthetic-data factory instead of
+    hand-rolling its own. The data is a noisy realization of
+    ``(Z_i @ B) * A[i] @ C.T`` per condition ``i``, which is exactly the
+    structure PARAFAC2 assumes, so a real ``pf2()`` fit on this data
+    recovers ``A``/``B``/``C`` up to the usual permutation/sign/scale
+    ambiguity.
+    """
+    rng = np.random.default_rng(seed)
+    B = rng.normal(size=(rank, rank))
+    C = rng.normal(size=(n_genes, rank))
+    A = rng.normal(size=(n_cond, rank))
+
+    X_list = []
+    cond_idx = []
+    for i in range(n_cond):
+        n_cells = int(rng.integers(*cells_per_cond))
+        Z = rng.normal(size=(n_cells, rank))
+        signal = (Z @ B) * A[i] @ C.T
+        noise = rng.normal(scale=0.2, size=signal.shape)
+        X_list.append(signal + noise)
+        cond_idx += [i] * n_cells
+
+    X = np.concatenate(X_list, axis=0).astype(np.float32)
+    n_cells_total = X.shape[0]
+
+    obs = pd.DataFrame(
+        {
+            "condition_unique_idxs": pd.Categorical(cond_idx),
+            "Condition": pd.Categorical([f"cond_{i}" for i in cond_idx]),
+        }
+    )
+    var = pd.DataFrame(
+        {"gene_name": [f"gene_{j}" for j in range(n_genes)]},
+        index=[f"gene_{j}" for j in range(n_genes)],
+    )
+
+    adata = anndata.AnnData(X=X, obs=obs, var=var)
+    adata.var["means"] = np.zeros(n_genes)
+    adata.obs_names = [f"cell_{i}" for i in range(n_cells_total)]
+    return adata
+
+
+def make_mock_factored_adata(
+    n_cells: int = 40,
+    n_genes: int = 25,
+    n_conditions: int = 6,
+    rank: int = 3,
+    seed: int = 0,
+    with_embedding: bool = False,
+) -> anndata.AnnData:
+    """Build an AnnData already populated with (random, not fitted) RISE
+    factors -- i.e. what ``pf2()`` would have attached -- for testing
+    downstream consumers (reordering, plotting, export) without paying for
+    an actual PARAFAC2 fit."""
+    rng = np.random.default_rng(seed)
+
+    A = rng.normal(size=(n_conditions, rank)).astype(np.float32)
+    C = rng.normal(size=(n_genes, rank)).astype(np.float32)
+    B = rng.normal(size=(rank, rank)).astype(np.float32)
+    weights = rng.random(rank).astype(np.float32)
+    projections, _ = np.linalg.qr(rng.normal(size=(n_cells, rank)))
+    projections = projections.astype(np.float32)
+
+    cond_idx = [i % n_conditions for i in range(n_cells)]
+    obs = pd.DataFrame(
+        {
+            "Condition": pd.Categorical([f"cond_{i}" for i in cond_idx]),
+            "condition_unique_idxs": pd.Categorical(cond_idx),
+            "Cell Type": pd.Categorical([f"type_{i % 3}" for i in range(n_cells)]),
+        }
+    )
+    var = pd.DataFrame(
+        {"means": np.zeros(n_genes)},
+        index=[f"gene_{j}" for j in range(n_genes)],
+    )
+
+    obsm = {"projections": projections, "weighted_projections": projections @ B}
+    if with_embedding:
+        obsm["X_pf2_PaCMAP"] = rng.normal(size=(n_cells, 2)).astype(np.float32)
+
+    adata = anndata.AnnData(
+        X=rng.normal(size=(n_cells, n_genes)).astype(np.float32),
+        obs=obs,
+        var=var,
+        uns={"Pf2_A": A, "Pf2_B": B, "Pf2_weights": weights},
+        varm=cast(Mapping[str, Sequence[Any]], {"Pf2_C": C}),
+        obsm=cast(Mapping[str, Sequence[Any]], obsm),
+    )
+    adata.obs_names = [f"cell_{i}" for i in range(n_cells)]
+    return adata
+
+
+@pytest.fixture
+def synthetic_pf2_adata():
+    return make_synthetic_pf2_data()
+
+
+@pytest.fixture
+def mock_factored_adata():
+    return make_mock_factored_adata()
diff --git a/scrise/tests/test_contracts.py b/scrise/tests/test_contracts.py
new file mode 100644
index 00000000..9f5acee9
--- /dev/null
+++ b/scrise/tests/test_contracts.py
@@ -0,0 +1,148 @@
+"""
+Contract / error-path tests.
+
+Every public function in `scrise.factorization` and `scrise.rank_selection`
+documents required AnnData keys and value ranges in its docstring. These
+tests pin down what actually happens when a caller violates that contract
+(missing keys, wrong dtypes, out-of-range parameters) so a future change
+can't silently turn a clear error into a confusing downstream crash (or
+worse, a silently wrong result) without a test noticing.
+"""
+
+from collections.abc import Mapping, Sequence
+from typing import Any, cast
+
+import anndata
+import numpy as np
+import pandas as pd
+import pytest
+
+from ..factorization import (
+    correct_conditions,
+    export_factors,
+    order_components_by_energy,
+    pf2,
+)
+from ..rank_selection import bicv
+from .conftest import make_synthetic_pf2_data
+
+
+def test_pf2_missing_condition_idxs_raises():
+    """pf2() requires X.obs['condition_unique_idxs']; without it the
+    caller gets a clear KeyError rather than a cryptic failure deep inside
+    the PARAFAC2 solver."""
+    X = anndata.AnnData(X=np.random.rand(10, 5).astype(np.float32))
+    with pytest.raises(KeyError, match="condition_unique_idxs"):
+        pf2(X, rank=2, doEmbedding=False, compress=None)
+
+
+def test_correct_conditions_missing_condition_idxs_raises():
+    X = anndata.AnnData(X=np.random.rand(10, 5).astype(np.float32))
+    with pytest.raises(KeyError, match="condition_unique_idxs"):
+        correct_conditions(X)
+
+
+def test_correct_conditions_none_X_raises_typeerror():
+    """correct_conditions has an explicit guard for X.X is None (e.g. a
+    factors-only AnnData produced by export_factors) -- it should fail
+    clearly rather than crash inside `.sum()`."""
+    obs = pd.DataFrame({"condition_unique_idxs": [0, 0, 1, 1]})
+    X = anndata.AnnData(X=None, obs=obs, shape=(4, 3))
+    X.uns["Pf2_A"] = np.ones((2, 2))
+    with pytest.raises(TypeError, match="X.X must not be None"):
+        correct_conditions(X)
+
+
+@pytest.mark.parametrize("missing_key", ["Pf2_A", "Pf2_B", "Pf2_weights"])
+def test_order_components_by_energy_missing_uns_key_raises(missing_key):
+    keys = {"Pf2_A", "Pf2_B", "Pf2_weights"}
+    rank, n_genes, n_conditions = 2, 5, 3
+    uns = {
+        "Pf2_A": np.ones((n_conditions, rank)),
+        "Pf2_B": np.ones((rank, rank)),
+        "Pf2_weights": np.ones(rank),
+    }
+    del uns[missing_key]
+    X = anndata.AnnData(
+        X=np.zeros((6, n_genes), dtype=np.float32),
+        uns=uns,
+        varm=cast(Mapping[str, Sequence[Any]], {"Pf2_C": np.ones((n_genes, rank))}),
+    )
+    with pytest.raises(KeyError, match=missing_key):
+        order_components_by_energy(X)
+    assert missing_key in keys  # sanity: parametrization matches the real keys
+
+
+@pytest.mark.parametrize(
+    "missing",
+    ["uns:Pf2_A", "uns:Pf2_B", "uns:Pf2_weights", "varm:Pf2_C", "obsm:projections"],
+)
+def test_export_factors_missing_any_required_field_raises_keyerror(tmp_path, missing):
+    n_cells, n_genes, rank = 10, 6, 2
+    fields = {
+        "uns:Pf2_A": ("uns", "Pf2_A", np.ones((3, rank))),
+        "uns:Pf2_B": ("uns", "Pf2_B", np.ones((rank, rank))),
+        "uns:Pf2_weights": ("uns", "Pf2_weights", np.ones(rank)),
+        "varm:Pf2_C": ("varm", "Pf2_C", np.ones((n_genes, rank))),
+        "obsm:projections": ("obsm", "projections", np.ones((n_cells, rank))),
+    }
+    uns, varm, obsm = {}, {}, {}
+    dest = {"uns": uns, "varm": varm, "obsm": obsm}
+    for key, (section, name, value) in fields.items():
+        if key == missing:
+            continue
+        dest[section][name] = value
+
+    X = anndata.AnnData(
+        X=np.zeros((n_cells, n_genes), dtype=np.float32), uns=uns, varm=varm, obsm=obsm
+    )
+    with pytest.raises(KeyError):
+        export_factors(X, str(tmp_path / "out.h5ad"))
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        {"held_out_cell_frac": 0.0},
+        {"held_out_cell_frac": 1.0},
+        {"held_out_cell_frac": -0.1},
+        {"held_out_gene_frac": 0.0},
+        {"held_out_gene_frac": 1.5},
+        {"n_repeats": 0},
+        {"n_repeats": -1},
+    ],
+)
+def test_bicv_rejects_invalid_arguments(kwargs):
+    X = make_synthetic_pf2_data(n_cond=4, n_genes=20, rank=2, seed=0)
+    with pytest.raises(ValueError):
+        bicv(X, [2], **kwargs)
+
+
+def test_bicv_rejects_rank_exceeding_max_feasible_rank():
+    X = make_synthetic_pf2_data(n_cond=3, n_genes=10, rank=2, seed=0)
+    with pytest.raises(ValueError, match="exceeds the maximum feasible rank"):
+        bicv(X, [10_000])
+
+
+def test_export_factors_output_directory_is_created(tmp_path):
+    """export_factors should create any missing parent directories for the
+    output path rather than failing with FileNotFoundError."""
+    n_cells, n_genes, rank = 10, 6, 2
+    X = anndata.AnnData(
+        X=np.zeros((n_cells, n_genes), dtype=np.float32),
+        uns={
+            "Pf2_A": np.random.rand(3, rank),
+            "Pf2_B": np.random.rand(rank, rank),
+            "Pf2_weights": np.random.rand(rank),
+        },
+        varm=cast(
+            Mapping[str, Sequence[Any]], {"Pf2_C": np.random.rand(n_genes, rank)}
+        ),
+        obsm=cast(
+            Mapping[str, Sequence[Any]],
+            {"projections": np.random.rand(n_cells, rank).astype(np.float32)},
+        ),
+    )
+    out_path = tmp_path / "nested" / "dir" / "out.h5ad"
+    export_factors(X, str(out_path))
+    assert out_path.exists()
diff --git a/scrise/tests/test_golden_regression.py b/scrise/tests/test_golden_regression.py
new file mode 100644
index 00000000..14621de7
--- /dev/null
+++ b/scrise/tests/test_golden_regression.py
@@ -0,0 +1,75 @@
+"""
+Golden regression test for the end-to-end pf2() fit.
+
+The invariant/contract tests elsewhere check *properties* that must hold
+for any input; this test instead pins the actual numbers produced by a
+fixed synthetic dataset + fixed random_state, so a silent change in
+behavior -- e.g. from a `parafac2` dependency bump changing its solver's
+convergence path, or an accidental change to `pf2`'s defaults -- shows up
+as a diff here even if it doesn't happen to violate any general property.
+
+The expected values below were captured from the current `pf2()` output on
+the fixture defined in this file; if they need to be regenerated (a
+deliberate algorithmic change, a new `parafac2` release that legitimately
+changes results), rerun this file's `if __name__ == "__main__"` block and
+paste the printed values back in.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+from parafac2.parafac2 import parafac2_nd
+
+from .. import pf2
+from .conftest import make_synthetic_pf2_data
+
+
+def _fixture():
+    X = make_synthetic_pf2_data(
+        n_cond=4, n_genes=25, rank=3, seed=7, cells_per_cond=(20, 30)
+    )
+    X.obs["condition_unique_idxs"] = pd.Categorical(
+        X.obs["condition_unique_idxs"].astype(int)
+    )
+    return X
+
+
+def test_pf2_golden_values_on_fixed_synthetic_fixture():
+    X = _fixture()
+
+    # Reconstruction quality (R2X) is deterministic for a fixed seed/data and
+    # is the single most sensitive summary of "did the fit change at all".
+    # parafac2_nd returns it directly, so use that rather than recomputing
+    # by hand (pf2() reorders/signs the factors before returning them, so a
+    # manual reconstruction from its output must exactly replicate that
+    # bookkeeping to agree -- fragile for what this test wants to check).
+    _, r2x = parafac2_nd(
+        X, rank=3, random_state=1, tol=1e-6, n_iter_max=100, compress=None
+    )
+    assert r2x == pytest.approx(0.9814689334613979, abs=1e-3)
+
+    result = pf2(
+        X, rank=3, random_state=1, doEmbedding=False, compress=None, max_iter=100
+    )
+    C = np.array(result.varm["Pf2_C"])
+    weights = np.array(result.uns["Pf2_weights"])
+
+    # Component weights, order, and which gene dominates each component are
+    # a much coarser (and thus more robust-to-numerical-noise) fingerprint
+    # of the fit than raw factor values.
+    np.testing.assert_allclose(weights, [22.2080, 56.6194, 24.3079], rtol=1e-2)
+
+    top_gene_per_component = np.argmax(np.abs(C), axis=0)
+    np.testing.assert_array_equal(top_gene_per_component, [12, 13, 5])
+
+
+if __name__ == "__main__":
+    X = _fixture()
+    result = pf2(
+        X, rank=3, random_state=1, doEmbedding=False, compress=None, max_iter=100
+    )
+    A = np.array(result.uns["Pf2_A"])
+    C = np.array(result.varm["Pf2_C"])
+    weights = np.array(result.uns["Pf2_weights"])
+    print("weights:", weights.tolist())
+    print("top_gene_per_component:", np.argmax(np.abs(C), axis=0).tolist())
diff --git a/scrise/tests/test_invariants.py b/scrise/tests/test_invariants.py
new file mode 100644
index 00000000..cba20083
--- /dev/null
+++ b/scrise/tests/test_invariants.py
@@ -0,0 +1,216 @@
+"""
+Property-based invariant tests for the core factorization/compression
+primitives.
+
+These tests fuzz shapes, ranks, and random seeds with Hypothesis rather than
+relying only on a handful of hand-picked examples, to surface edge cases
+(rank == 1, a single condition, near-degenerate columns, ...) that
+example-based tests tend to miss. Each test asserts a mathematical property
+that must hold for *any* valid input, not just the specific arrays checked
+elsewhere in the suite.
+"""
+
+import anndata
+import numpy as np
+import pytest
+from hypothesis import HealthCheck, given, settings
+from hypothesis import strategies as st
+from hypothesis.extra.numpy import arrays
+
+from ..factorization import (
+    canonical_component_signs,
+    correct_conditions,
+    match_components_across_ranks,
+    order_components_by_energy,
+)
+from ..opq import OPQQuantizer
+
+# Bounded, finite, non-degenerate floats: avoids the near-zero-norm columns
+# that would make sign/cosine-similarity assertions numerically meaningless.
+_finite_floats = st.floats(
+    min_value=-100.0,
+    max_value=100.0,
+    allow_nan=False,
+    allow_infinity=False,
+    width=32,
+)
+
+_slow_settings = settings(
+    max_examples=15, deadline=None, suppress_health_check=[HealthCheck.too_slow]
+)
+
+
+def _nonzero_column_matrix(n_rows: int, n_cols: int):
+    """A matrix strategy in which every column has norm bounded away from 0,
+    generated by adding a fixed unit spike to each column of a random
+    matrix (cheaper and more reliable than a `.filter()` on norm)."""
+    return arrays(
+        dtype=np.float64, shape=(n_rows, n_cols), elements=_finite_floats
+    ).map(lambda C: C + np.eye(n_rows, n_cols) * 10.0)
+
+
+@given(
+    n_rows=st.integers(min_value=1, max_value=12),
+    n_cols=st.integers(min_value=1, max_value=6),
+    data=st.data(),
+)
+@_slow_settings
+def test_canonical_component_signs_idempotent_and_positive_max(n_rows, n_cols, data):
+    """Signing is idempotent, and always leaves the largest-magnitude entry
+    of each column positive -- the defining property from the docstring."""
+    C = data.draw(_nonzero_column_matrix(n_rows, n_cols))
+
+    signs = canonical_component_signs(C)
+    assert set(np.unique(signs)).issubset({-1.0, 1.0})
+
+    C_signed = C * signs
+    max_idx = np.argmax(np.abs(C_signed), axis=0)
+    assert np.all(C_signed[max_idx, np.arange(n_cols)] >= 0)
+
+    # Applying the canonical sign to an already-canonicalized matrix is a
+    # no-op (all signs come back +1).
+    signs_again = canonical_component_signs(C_signed)
+    np.testing.assert_array_equal(signs_again, np.ones(n_cols))
+
+
+def _mock_adata(A, B, C, weights, projections):
+    n_cells = projections.shape[0]
+    n_genes = C.shape[0]
+    n_conditions = A.shape[0]
+    obs = {"Condition": [f"cond_{i % n_conditions}" for i in range(n_cells)]}
+    var = {"gene_name": [f"gene_{j}" for j in range(n_genes)]}
+    return anndata.AnnData(
+        obs=obs,
+        var=var,
+        uns={"Pf2_A": A, "Pf2_B": B, "Pf2_weights": weights},
+        varm={"Pf2_C": C},
+        obsm={"projections": projections},
+    )
+
+
+@given(
+    rank=st.integers(min_value=1, max_value=5),
+    n_conditions=st.integers(min_value=1, max_value=6),
+    n_genes=st.integers(min_value=2, max_value=15),
+    n_cells=st.integers(min_value=3, max_value=20),
+    seed=st.integers(min_value=0, max_value=2**31 - 1),
+)
+@_slow_settings
+def test_order_components_by_energy_is_a_permutation(
+    rank, n_conditions, n_genes, n_cells, seed
+):
+    """For any shape/seed, reordering by energy must: (1) leave energy sorted
+    descending, and (2) be a genuine permutation of the original columns --
+    weighted_projections after reordering must equal the original
+    weighted_projections with columns permuted (order_components_by_energy
+    keeps B as the unflipped sign reference, so this holds exactly)."""
+    rng = np.random.default_rng(seed)
+    A = rng.normal(size=(n_conditions, rank))
+    C = rng.normal(size=(n_genes, rank)) + np.eye(n_genes, rank) * 10.0
+    B = rng.normal(size=(rank, rank))
+    weights = rng.random(rank)
+    if rank <= n_cells:
+        projections, _ = np.linalg.qr(rng.normal(size=(n_cells, rank)))
+    else:
+        projections = rng.normal(size=(n_cells, rank))
+
+    before_wp = projections @ B
+    adata = _mock_adata(A.copy(), B.copy(), C.copy(), weights.copy(), projections)
+
+    ordered = order_components_by_energy(adata)
+
+    new_A = np.array(ordered.uns["Pf2_A"])
+    new_C = np.array(ordered.varm["Pf2_C"])
+    new_energy = np.linalg.norm(new_A, axis=0) * np.linalg.norm(new_C, axis=0)
+    assert np.all(np.diff(new_energy) <= 1e-6)
+
+    energy = np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+    order = np.argsort(energy)[::-1]
+    np.testing.assert_allclose(
+        np.array(ordered.obsm["weighted_projections"]), before_wp[:, order], atol=1e-5
+    )
+    np.testing.assert_allclose(new_A.shape, A.shape)
+
+
+@given(
+    n_genes=st.integers(min_value=2, max_value=15),
+    rank=st.integers(min_value=1, max_value=6),
+    seed=st.integers(min_value=0, max_value=2**31 - 1),
+)
+@_slow_settings
+def test_match_components_across_ranks_self_match_is_identity(n_genes, rank, seed):
+    """Matching a (sign-canonicalized) gene factor matrix against an
+    identical copy of itself must recover the identity permutation with
+    every component matched -- the base case any correct matching
+    implementation has to satisfy."""
+    rng = np.random.default_rng(seed)
+    C = rng.normal(size=(n_genes, rank)) + np.eye(n_genes, rank) * 10.0
+    C = C * canonical_component_signs(C)
+
+    matched_pairs, unmatched_high = match_components_across_ranks(C, C.copy())
+
+    assert matched_pairs.shape[0] == rank
+    assert unmatched_high.size == 0
+    np.testing.assert_array_equal(np.sort(matched_pairs[:, 0]), np.arange(rank))
+    np.testing.assert_array_equal(np.sort(matched_pairs[:, 1]), np.arange(rank))
+
+
+@given(
+    n_conditions=st.integers(min_value=2, max_value=8),
+    rank=st.integers(min_value=1, max_value=5),
+    scale=st.floats(min_value=0.1, max_value=50.0, allow_nan=False),
+    seed=st.integers(min_value=0, max_value=2**31 - 1),
+)
+@_slow_settings
+def test_correct_conditions_is_invariant_to_uniform_rescaling(
+    n_conditions, rank, scale, seed
+):
+    """correct_conditions fits an intercept-including linear regression of
+    the (geometric-mean) condition factor magnitude against read counts,
+    then divides by the fit. Ordinary least squares is equivariant under
+    scaling its target by a positive constant c (both the fitted slope and
+    intercept scale by c), so the c introduced by rescaling Pf2_A cancels
+    against the c picked up by the fitted correction term: the corrected
+    output must be unchanged, for *any* underlying count data."""
+    rng = np.random.default_rng(seed)
+    n_cells_per_cond = 5
+    n_cells = n_conditions * n_cells_per_cond
+    n_genes = 8
+
+    cond_idx = np.repeat(np.arange(n_conditions), n_cells_per_cond)
+    counts_X = rng.integers(1, 50, size=(n_cells, n_genes)).astype(np.float64)
+    A = rng.uniform(0.1, 10.0, size=(n_conditions, rank))
+
+    def _make(A_):
+        adata = anndata.AnnData(X=counts_X.copy())
+        adata.obs["condition_unique_idxs"] = cond_idx
+        adata.uns["Pf2_A"] = A_
+        return adata
+
+    baseline = correct_conditions(_make(A.copy()))
+    scaled = correct_conditions(_make(A.copy() * scale))
+
+    np.testing.assert_allclose(scaled, baseline, rtol=1e-4, atol=1e-6)
+
+
+@pytest.mark.parametrize("D", [8, 16])
+def test_opq_reconstruction_r2_nondecreasing_in_num_subquantizers(D):
+    """Fidelity (R^2) must not get worse as the projection matrix is split
+    into more, finer-grained sub-quantizers -- more codebooks can only add
+    representational capacity."""
+    rng = np.random.default_rng(0)
+    N = 300
+    P = rng.normal(size=(N, D)).astype(np.float32)
+    P, _ = np.linalg.qr(P)
+
+    prev_r2 = -np.inf
+    for M in sorted({1, 2, 4, D}):
+        if M > D:
+            continue
+        q = OPQQuantizer(M=M, random_state=0)
+        _, _, r2 = q.fit_transform(P)
+        # A small numerical slack: OPQ's alternating optimization is a
+        # heuristic, not an exact solver, so strict monotonicity can be
+        # violated by noise near saturation (r2 close to 1).
+        assert r2 >= prev_r2 - 0.05
+        prev_r2 = r2
diff --git a/scrise/tests/test_plotting_factors.py b/scrise/tests/test_plotting_factors.py
new file mode 100644
index 00000000..3d20d4f5
--- /dev/null
+++ b/scrise/tests/test_plotting_factors.py
@@ -0,0 +1,114 @@
+"""
+Smoke + data-mapping tests for scrise.plotting.factors and
+scrise.plotting.rank_selection.
+
+We don't assert on rendered pixels. Instead we check two things a plotting
+bug is actually likely to break: (1) the function runs without raising on
+both a normal input and small edge cases (rank 1, a single row), and (2)
+where cheaply checkable, that the *right numbers* ended up on the axes
+(e.g. plot_bicv_r2x puts the correct R2X values on the y-axis) rather than
+just "some heatmap was drawn".
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import pytest
+
+from ..plotting.factors import (
+    plot_condition_factors,
+    plot_eigenstate_factors,
+    plot_gene_factors,
+    reorder_table,
+)
+from ..plotting.rank_selection import plot_bicv_r2x
+from ..rank_selection import bicv
+from .conftest import make_mock_factored_adata, make_synthetic_pf2_data
+
+
+@pytest.mark.parametrize("rank,n_cond", [(3, 6), (1, 4)])
+def test_plot_condition_factors_smoke(rank, n_cond):
+    adata = make_mock_factored_adata(n_conditions=n_cond, rank=rank)
+    # Force all-positive values so the default log_transform=True doesn't
+    # hit log(negative).
+    adata.uns["Pf2_A"] = np.abs(adata.uns["Pf2_A"]) + 0.1
+
+    fig, ax = plt.subplots()
+    plot_condition_factors(adata, ax)
+    assert ax.get_xlabel() == "Component"
+    assert len(ax.get_yticklabels()) == n_cond
+
+
+def test_plot_condition_factors_with_group_labels_and_legend():
+    n_cond, rank = 8, 3
+    adata = make_mock_factored_adata(n_conditions=n_cond, rank=rank)
+    adata.uns["Pf2_A"] = np.abs(adata.uns["Pf2_A"]) + 0.1
+
+    yt = pd.Series(np.unique(adata.obs["Condition"]))
+    group_labels = pd.Series(
+        ["groupA" if i % 2 == 0 else "groupB" for i in range(len(yt))]
+    )
+
+    fig, ax = plt.subplots()
+    plot_condition_factors(adata, ax, cond_group_labels=group_labels, group_cond=True)
+    # A legend patch per unique group should have been added.
+    legend = ax.get_legend()
+    assert legend is not None
+    assert len(legend.get_texts()) == group_labels.nunique()
+
+
+def test_plot_eigenstate_factors_smoke():
+    adata = make_mock_factored_adata(rank=4)
+    fig, ax = plt.subplots()
+    plot_eigenstate_factors(adata, ax)
+    assert ax.get_xlabel() == "Component"
+    assert len(ax.get_xticklabels()) == 4
+
+
+@pytest.mark.parametrize("trim", [True, False])
+def test_plot_gene_factors_smoke(trim):
+    adata = make_mock_factored_adata(n_genes=30, rank=3)
+    fig, ax = plt.subplots()
+    plot_gene_factors(adata, ax, trim=trim)
+    assert ax.get_xlabel() == "Component"
+
+
+def test_plot_gene_factors_trim_reduces_or_keeps_gene_count():
+    adata = make_mock_factored_adata(n_genes=30, rank=3)
+    fig1, ax1 = plt.subplots()
+    plot_gene_factors(adata, ax1, trim=False, weight=0.08)
+    n_all = len(ax1.get_yticklabels())
+
+    fig2, ax2 = plt.subplots()
+    plot_gene_factors(adata, ax2, trim=True, weight=0.08)
+    n_trimmed = len(ax2.get_yticklabels())
+
+    assert n_trimmed <= n_all
+
+
+def test_reorder_table_groups_rows_by_dominant_component():
+    """Rows should come back grouped by which column has their largest
+    magnitude entry."""
+    projs = np.array(
+        [
+            [0.1, 5.0, 0.0],  # dominant in col 1
+            [3.0, 0.0, 0.1],  # dominant in col 0
+            [0.0, 0.2, 4.0],  # dominant in col 2
+            [1.0, 0.0, 0.1],  # dominant in col 0, smaller than row 1
+        ]
+    )
+    ind = reorder_table(projs)
+    dominant = np.argmax(np.abs(projs[ind]), axis=1)
+    assert np.all(np.diff(dominant) >= 0)
+
+
+def test_plot_bicv_r2x_smoke_and_axis_labels():
+    X = make_synthetic_pf2_data(n_cond=4, n_genes=20, rank=2, seed=0)
+    results = bicv(X, [2, 4], n_repeats=1, random_state=0, max_iter=30)
+
+    fig, ax = plt.subplots()
+    plot_bicv_r2x(results, ax)
+
+    assert ax.get_xlabel() == "Rank"
+    assert ax.get_ylabel() == "R2X"
+    assert ax.get_legend() is not None
diff --git a/scrise/tests/test_plotting_general.py b/scrise/tests/test_plotting_general.py
new file mode 100644
index 00000000..360c8758
--- /dev/null
+++ b/scrise/tests/test_plotting_general.py
@@ -0,0 +1,148 @@
+"""
+Smoke + data-shaping tests for scrise.plotting.general.
+
+Several functions here (cell_count_perc_df, avegene_per_status, ...) are
+really data-reshaping helpers that happen to feed a plot; those are tested
+on their returned DataFrame contents directly rather than the plot.
+"""
+
+import anndata
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import pytest
+
+from ..plotting.general import (
+    avegene_per_status,
+    cell_count_perc_df,
+    cell_count_perc_lupus_df,
+    gene_plot_cells,
+    plot_avegene_per_celltype,
+    plot_cell_gene_corr,
+    plot_r2x,
+    rotate_xaxis,
+    rotate_yaxis,
+)
+from .conftest import make_synthetic_pf2_data
+
+
+def _two_gene_adata(n_cells=40, seed=0):
+    rng = np.random.default_rng(seed)
+    X = rng.normal(size=(n_cells, 2)).astype(np.float32)
+    obs = pd.DataFrame(
+        {
+            "hue": pd.Categorical([f"grp_{i % 3}" for i in range(n_cells)]),
+            "Cell Type": pd.Categorical([f"type_{i % 2}" for i in range(n_cells)]),
+        }
+    )
+    var = pd.DataFrame({"means": [0.0, 0.0]}, index=["gene_a", "gene_b"])
+    return anndata.AnnData(X=X, obs=obs, var=var)
+
+
+def test_rotate_xaxis_and_yaxis_set_tick_rotation():
+    fig, ax = plt.subplots()
+    rotate_xaxis(ax, rotation=45)
+    rotate_yaxis(ax, rotation=30)
+    # matplotlib doesn't expose rotation as a simple getter pre-draw; the
+    # smoke check is that these don't raise and the axes object still works.
+    ax.figure.canvas.draw()
+
+
+def test_cell_count_perc_df_percentages_sum_to_100_per_condition():
+    n_cells = 60
+    obs = pd.DataFrame(
+        {
+            "Cell Type": np.random.choice(["A", "B", "C"], size=n_cells),
+            "Condition": np.random.choice(["cond_0", "cond_1"], size=n_cells),
+        }
+    )
+    X = anndata.AnnData(X=np.zeros((n_cells, 2), dtype=np.float32), obs=obs)
+
+    df = cell_count_perc_df(X)
+    totals = df.groupby("Condition")["Cell Type Percentage"].sum()
+    np.testing.assert_allclose(totals.to_numpy(), 100.0, rtol=1e-6)
+    assert df["Cell Count"].sum() == n_cells
+
+
+def test_cell_count_perc_lupus_df_attaches_metadata_columns():
+    n_cells = 40
+    obs = pd.DataFrame(
+        {
+            "Cell Type": np.random.choice(["A", "B"], size=n_cells),
+            "Condition": np.repeat(["cond_0", "cond_1"], n_cells // 2),
+            "SLE_status": np.repeat(["healthy", "SLE"], n_cells // 2),
+            "Processing_Cohort": np.repeat([1, 2], n_cells // 2),
+            "condition_unique_idxs": np.repeat([0, 1], n_cells // 2),
+        }
+    )
+    X = anndata.AnnData(X=np.zeros((n_cells, 2), dtype=np.float32), obs=obs)
+
+    df = cell_count_perc_lupus_df(X)
+    assert {"SLE_status", "Processing_Cohort", "condition_unique_idxs"}.issubset(
+        df.columns
+    )
+    cond0_status = df.loc[df["Condition"] == "cond_0", "SLE_status"].unique()
+    assert list(cond0_status) == ["healthy"]
+
+
+def test_avegene_per_status_returns_expected_columns():
+    n_cells = 30
+    X = _two_gene_adata(n_cells)
+    X.obs["SLE_status"] = np.random.choice(["healthy", "SLE"], size=n_cells)
+    X.obs["Condition"] = np.random.choice(["cond_0", "cond_1"], size=n_cells)
+
+    df = avegene_per_status(X[:, "gene_a"], "gene_a")
+    assert {
+        "Status",
+        "Cell Type",
+        "Gene",
+        "Condition",
+        "Average Gene Expression",
+    } <= set(df.columns)
+    assert set(df["Gene"]) == {"gene_a"}
+
+
+def test_plot_avegene_per_celltype_smoke():
+    n_cells = 30
+    X = _two_gene_adata(n_cells)
+    X.obs["Condition"] = np.random.choice(["cond_0", "cond_1"], size=n_cells)
+
+    fig, ax = plt.subplots()
+    plot_avegene_per_celltype(X, ["gene_a", "gene_b"], ax)
+
+
+def test_gene_plot_cells_smoke_and_shape_assertion():
+    X = _two_gene_adata(20)
+    fig, ax = plt.subplots()
+    gene_plot_cells(X, hue="hue", ax=ax)
+
+    # gene_plot_cells asserts exactly 2 genes/columns are present.
+    X_three_genes = anndata.AnnData(
+        X=np.zeros((5, 3), dtype=np.float32),
+        var=pd.DataFrame({"means": [0.0, 0.0, 0.0]}, index=["a", "b", "c"]),
+        obs=pd.DataFrame({"hue": ["x"] * 5, "Cell Type": ["t"] * 5}),
+    )
+    fig2, ax2 = plt.subplots()
+    with pytest.raises(AssertionError):
+        gene_plot_cells(X_three_genes, hue="hue", ax=ax2)
+
+
+def test_plot_cell_gene_corr_smoke_with_missing_pivot_columns():
+    """When the requested (gene, cell-type) combination isn't present after
+    pivoting, plot_cell_gene_corr should fall back to an empty frame rather
+    than raising a KeyError."""
+    X = _two_gene_adata(20)
+    fig, ax = plt.subplots()
+    plot_cell_gene_corr(
+        X, hue="hue", cells=["type_0", "type_1"], ax=ax, unique=["grp_0"]
+    )
+
+
+def test_plot_r2x_smoke_and_axes_labels():
+    X = make_synthetic_pf2_data(n_cond=4, n_genes=15, rank=2, seed=0)
+    fig, ax = plt.subplots()
+    plot_r2x(X, np.array([1, 2, 3]), ax, compress=None)
+
+    assert ax.get_xlabel() == "Number of Components"
+    assert ax.get_ylabel() == "Variance Explained"
+    assert ax.get_legend() is not None
diff --git a/scrise/tests/test_plotting_pacmap.py b/scrise/tests/test_plotting_pacmap.py
new file mode 100644
index 00000000..3199f745
--- /dev/null
+++ b/scrise/tests/test_plotting_pacmap.py
@@ -0,0 +1,89 @@
+"""
+Smoke tests for scrise.plotting.pacmap.
+
+PaCMAP itself is not exercised here -- these functions only ever consume a
+precomputed 2D embedding (X.obsm["X_pf2_PaCMAP"]), so we supply one
+directly rather than paying for an actual PaCMAP fit.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+
+from ..plotting.pacmap import (
+    _get_canvas,
+    _to_hex,
+    assign_labels,
+    plot_gene_pacmap,
+    plot_labels_pacmap,
+    plot_wp_pacmap,
+)
+from .conftest import make_mock_factored_adata
+
+
+def test_get_canvas_bounds_pad_outward_from_data_range():
+    points = np.array([[0.0, 0.0], [10.0, 20.0]])
+    canvas = _get_canvas(points)
+    assert canvas.x_range[0] <= 0.0
+    assert canvas.x_range[1] >= 10.0
+    assert canvas.y_range[0] <= 0.0
+    assert canvas.y_range[1] >= 20.0
+
+
+def test_to_hex_returns_hex_strings():
+    colors = _to_hex(plt.get_cmap("tab10")(np.linspace(0, 1, 4)))
+    assert len(colors) == 4
+    assert all(c.startswith("#") for c in colors)
+
+
+def test_assign_labels_sets_pacmap_axis_labels():
+    fig, ax = plt.subplots()
+    out = assign_labels(ax)
+    assert out is ax
+    assert ax.get_xlabel() == "PaCMAP1"
+    assert ax.get_ylabel() == "PaCMAP2"
+    assert list(ax.get_xticks()) == []
+    assert list(ax.get_yticks()) == []
+
+
+def test_plot_gene_pacmap_smoke():
+    adata = make_mock_factored_adata(n_cells=50, n_genes=10, with_embedding=True)
+    gene = adata.var_names[0]
+    fig, ax = plt.subplots()
+    plot_gene_pacmap(gene, adata, ax)
+    assert ax.get_title() == gene
+
+
+@pytest.mark.parametrize("cmp", [1, 3])
+def test_plot_wp_pacmap_smoke(cmp):
+    adata = make_mock_factored_adata(n_cells=50, rank=3, with_embedding=True)
+    fig, ax = plt.subplots()
+    plot_wp_pacmap(adata, cmp=cmp, ax=ax)
+    assert ax.get_title() == f"Cmp. {cmp}"
+
+
+def test_plot_labels_pacmap_smoke_and_legend_matches_categories():
+    adata = make_mock_factored_adata(n_cells=60, with_embedding=True)
+    fig, ax = plt.subplots()
+    plot_labels_pacmap(adata, labelType="Cell Type", ax=ax)
+
+    legend = ax.get_legend()
+    assert legend is not None
+    n_categories = adata.obs["Cell Type"].nunique()
+    assert len(legend.get_texts()) == n_categories
+
+
+def test_plot_labels_pacmap_condition_filter_collapses_to_other():
+    adata = make_mock_factored_adata(n_cells=60, with_embedding=True)
+    present = adata.obs["Cell Type"].unique().tolist()
+    keep = [present[0]]
+
+    fig, ax = plt.subplots()
+    plot_labels_pacmap(adata, labelType="Cell Type", ax=ax, condition=keep)
+
+    legend = ax.get_legend()
+    assert legend is not None
+    legend_labels = {t.get_text() for t in legend.get_texts()}
+    # Everything not in `keep` should be collapsed into a single "Other"
+    # category, so at most len(keep) + 1 categories should appear.
+    assert legend_labels <= set(keep) | {"Other"}
diff --git a/scrise/tests/test_plotting_stability.py b/scrise/tests/test_plotting_stability.py
new file mode 100644
index 00000000..4a3c14f3
--- /dev/null
+++ b/scrise/tests/test_plotting_stability.py
@@ -0,0 +1,99 @@
+"""
+Tests for scrise.plotting.stability (Factor Match Score and its plots).
+
+calculateFMS and resample are pure/cheap enough to test directly and
+precisely. plot_fms_diff_ranks / plot_fms_percent_drop each refit pf2()
+several times, so they're only smoke-tested, on a tiny synthetic dataset
+and the smallest run/rank counts that still exercise the real code path.
+"""
+
+import anndata
+import matplotlib.pyplot as plt
+import numpy as np
+
+from ..plotting.stability import calculateFMS, plot_fms_diff_ranks, resample
+from .conftest import make_synthetic_pf2_data
+
+
+def _mock_factors_adata(A, B, C, weights):
+    return anndata.AnnData(
+        X=np.zeros((2, C.shape[0]), dtype=np.float32),
+        uns={"Pf2_A": A, "Pf2_B": B, "Pf2_weights": weights},
+        varm={"Pf2_C": C},
+    )
+
+
+def test_calculate_fms_identical_decompositions_score_near_one():
+    rng = np.random.default_rng(0)
+    n_cond, n_genes, rank = 5, 20, 3
+    A = rng.normal(size=(n_cond, rank))
+    B = rng.normal(size=(rank, rank))
+    C = rng.normal(size=(n_genes, rank))
+    weights = rng.random(rank)
+
+    X = _mock_factors_adata(A.copy(), B.copy(), C.copy(), weights.copy())
+    Y = _mock_factors_adata(A.copy(), B.copy(), C.copy(), weights.copy())
+
+    assert np.isclose(calculateFMS(X, Y), 1.0, atol=1e-6)
+
+
+def test_calculate_fms_permuted_components_still_scores_near_one():
+    """FMS should be invariant to a permutation of component order -- it's
+    matching components, not comparing them positionally."""
+    rng = np.random.default_rng(1)
+    n_cond, n_genes, rank = 5, 20, 3
+    A = rng.normal(size=(n_cond, rank))
+    B = rng.normal(size=(rank, rank))
+    C = rng.normal(size=(n_genes, rank))
+    weights = rng.random(rank)
+
+    perm = np.array([2, 0, 1])
+    X = _mock_factors_adata(A.copy(), B.copy(), C.copy(), weights.copy())
+    Y = _mock_factors_adata(
+        A[:, perm].copy(), B[:, perm].copy(), C[:, perm].copy(), weights[perm].copy()
+    )
+
+    assert np.isclose(calculateFMS(X, Y), 1.0, atol=1e-6)
+
+
+def test_calculate_fms_unrelated_decompositions_scores_low():
+    rng = np.random.default_rng(2)
+    n_cond, n_genes, rank = 5, 20, 3
+    A1, B1, C1 = (
+        rng.normal(size=(n_cond, rank)),
+        rng.normal(size=(rank, rank)),
+        rng.normal(size=(n_genes, rank)),
+    )
+    A2, B2, C2 = (
+        rng.normal(size=(n_cond, rank)),
+        rng.normal(size=(rank, rank)),
+        rng.normal(size=(n_genes, rank)),
+    )
+    weights = rng.random(rank)
+
+    X = _mock_factors_adata(A1, B1, C1, weights.copy())
+    Y = _mock_factors_adata(A2, B2, C2, weights.copy())
+
+    assert calculateFMS(X, Y) < 0.5
+
+
+def test_resample_preserves_shape_and_samples_with_replacement():
+    n_cells, n_genes = 50, 10
+    X = anndata.AnnData(
+        X=np.arange(n_cells * n_genes, dtype=np.float32).reshape(n_cells, n_genes)
+    )
+    np.random.seed(0)
+    resampled = resample(X)
+    assert resampled.shape == X.shape
+    # With replacement over 50 draws, near-certain to hit a duplicate row.
+    assert len(np.unique(np.asarray(resampled.X)[:, 0])) < n_cells
+
+
+def test_plot_fms_diff_ranks_smoke():
+    X = make_synthetic_pf2_data(n_cond=4, n_genes=15, rank=2, seed=0)
+    fig, ax = plt.subplots()
+    plot_fms_diff_ranks(X, ax, ranksList=[2], runs=1, compress=None)
+
+    assert ax.get_xlabel() == "Component"
+    assert ax.get_ylabel() == "FMS"
+    assert ax.get_ylim() == (0.0, 1.0)
diff --git a/uv.lock b/uv.lock
index c5ac43e6..01b3904b 100644
--- a/uv.lock
+++ b/uv.lock
@@ -613,6 +613,71 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/4e/14/32cf2aae083c74b95678875e11d25d0cb9e51a33d27f4218eb9aeacfcd70/hdf5plugin-7.0.0-py3-none-win_amd64.whl", hash = "sha256:2e052af8d7848e8bac92646584617503a08bb9b466cfa810a49ecd93e89b7ffa", size = 3523827, upload-time = "2026-06-25T20:59:37.758Z" },
 ]
 
+[[package]]
+name = "hypothesis"
+version = "6.167.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "sortedcontainers" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/c2/c9/8cee74c1390b2932406faaab76980f18946f258fa5a8afca17189b3bc655/hypothesis-6.167.1.tar.gz", hash = "sha256:62eefcb4d2791423626e9901c3027a6e0c5ffda2ac0b44b3c7e797ab9d2d5a4c", size = 505849, upload-time = "2026-08-30T19:53:09.05Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/53/4d/3592ca336deafbd3e9b0f47dc4c727aa32d30e765ef6370da8ecd590d388/hypothesis-6.167.1-cp310-abi3-macosx_10_12_x86_64.whl", hash = "sha256:d28118fd70e4e15ff9c308a98b312b544b6145ae45aaa3b566328c1fdee8058f", size = 785476, upload-time = "2026-08-30T19:51:02.154Z" },
+    { url = "https://files.pythonhosted.org/packages/18/bf/e33c431148994cbcb3332c6df94b833ecfb4aa6a8e51ea4b83da55ddd581/hypothesis-6.167.1-cp310-abi3-macosx_11_0_arm64.whl", hash = "sha256:e517be7f82a0a917758cc489a88b826b5371f56381fd94b4a8a09ce82d8de406", size = 781033, upload-time = "2026-08-30T19:51:27.314Z" },
+    { url = "https://files.pythonhosted.org/packages/94/a3/e0de9a82c7e790a1def0801076e0ef43110f98e95ed54a3554877d0cb66d/hypothesis-6.167.1-cp310-abi3-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:26f8cec74c4fad7aeb0852cb34c2134b16db05d878ad3946a53337dace7016f4", size = 1117814, upload-time = "2026-08-30T19:53:04.109Z" },
+    { url = "https://files.pythonhosted.org/packages/71/a4/8dd6bdc909324d1c39da1c86d65f75512ae049c159952af4cfe8feb5f8d4/hypothesis-6.167.1-cp310-abi3-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:fd202d02d129197a5e771f8a11c7d30559927284c23ec3a8bd4f37a7955964d1", size = 1141639, upload-time = "2026-08-30T19:51:50.399Z" },
+    { url = "https://files.pythonhosted.org/packages/fa/bd/c13ed6145c360d0770415efd7d5a7e63c29905aeef52ab88004fe7e7f924/hypothesis-6.167.1-cp310-abi3-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:8b1e393ab01b71f683ba2a783785871cc6b81a6e41017780c64a5bc0b99759ae", size = 1143334, upload-time = "2026-08-30T19:50:38.045Z" },
+    { url = "https://files.pythonhosted.org/packages/c0/a8/060d79ed8504b54ced9ad16f33d674b1b98a9debe9733c02709d7dd5c71c/hypothesis-6.167.1-cp310-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:8c385c7741893404306f9e5559ab3835432e85e7c153e25f854c502c410bbcbb", size = 1163345, upload-time = "2026-08-30T19:53:06.583Z" },
+    { url = "https://files.pythonhosted.org/packages/cb/f7/6e68e2b705f729a6f7b4f41030022b1a5264c5434d3bcd917233d6801c6a/hypothesis-6.167.1-cp310-abi3-manylinux_2_31_riscv64.whl", hash = "sha256:94920ca1fae70c26b0bd3fabbeef9437ffc17a39fe85696fb9a86187d92f6dba", size = 1123029, upload-time = "2026-08-30T19:52:03.134Z" },
+    { url = "https://files.pythonhosted.org/packages/a5/c0/d274fe37ed5ecd5ad8ed555edc1f5e2abc8e1c3be3d5404b7edd5cc353a8/hypothesis-6.167.1-cp310-abi3-manylinux_2_5_i686.manylinux1_i686.whl", hash = "sha256:b8d90ded2ffdc7e56b5e571993f384b52fade0a7b424e614f999cc2491789970", size = 1154003, upload-time = "2026-08-30T19:50:55.053Z" },
+    { url = "https://files.pythonhosted.org/packages/b4/2f/2b5bb386f43fc965eb86fd69fcb2bd62c08cb6d7c6708a40dc39b3b97440/hypothesis-6.167.1-cp310-abi3-musllinux_1_2_aarch64.whl", hash = "sha256:40cd5de7dd252942a08480639f5850594b1aca4a463e8a7f15e1fb6c2c3760c1", size = 1293729, upload-time = "2026-08-30T19:52:59.481Z" },
+    { url = "https://files.pythonhosted.org/packages/ac/32/22436b072d79011fe81abb933edcd2476057c7b971588c5f3caf07519a88/hypothesis-6.167.1-cp310-abi3-musllinux_1_2_armv7l.whl", hash = "sha256:b9c33f921ddc7fea93660eca408b25fe755516e22ec7ab21cb9951031f1cd608", size = 1419248, upload-time = "2026-08-30T19:50:52.903Z" },
+    { url = "https://files.pythonhosted.org/packages/9e/3e/abf39faaff0f78112112a82316a5c9fe472574480c1ecee526734775b812/hypothesis-6.167.1-cp310-abi3-musllinux_1_2_ppc64le.whl", hash = "sha256:495989cf0a5ee03f7f9598ee9efeaabf15fd861ec52b5a9d6435849453e17e5d", size = 1274903, upload-time = "2026-08-30T19:52:27.25Z" },
+    { url = "https://files.pythonhosted.org/packages/29/83/69b89ed5692ba3aad117facfb9ce099633a22c35acc3d64829b72253ec8c/hypothesis-6.167.1-cp310-abi3-musllinux_1_2_riscv64.whl", hash = "sha256:bc73c46ce8ff93b0eb220f2b75adcbd9fcc9112078a74d3522d2532ad8069bad", size = 1294185, upload-time = "2026-08-30T19:50:35.038Z" },
+    { url = "https://files.pythonhosted.org/packages/2b/1e/55dfcbe45c72df0a5c5b86a6b7c9365121acab69c2cc060bd55336a48c8f/hypothesis-6.167.1-cp310-abi3-musllinux_1_2_x86_64.whl", hash = "sha256:36e83e1d7e97aaacbf6cd778e14a841344f848a674b20dfe4fe997546a6a2151", size = 1330013, upload-time = "2026-08-30T19:52:54.841Z" },
+    { url = "https://files.pythonhosted.org/packages/3b/a6/0a36cead4ccff58bedb5d1aa2f894f7a580c317424b4fa788c6b22232b2d/hypothesis-6.167.1-cp310-abi3-win32.whl", hash = "sha256:fb4d87454d2459c2ccb541a4c61c92ce13058b91305ed3304695a409a1d886e4", size = 671942, upload-time = "2026-08-30T19:51:32.732Z" },
+    { url = "https://files.pythonhosted.org/packages/ec/5b/360285ed42109f5ef48d98ca9ffcf71d1477130c0ff1f1a8d535c8507259/hypothesis-6.167.1-cp310-abi3-win_amd64.whl", hash = "sha256:5e35f98b427bf438a946203426b485dd5b62485f3d5a69a0e0862870a545e518", size = 678637, upload-time = "2026-08-30T19:50:33.573Z" },
+    { url = "https://files.pythonhosted.org/packages/79/2d/cc084c1a8bfa296048ec0461f0fa11731abe0097f5c2310c8d39d94c8dc4/hypothesis-6.167.1-cp310-abi3-win_arm64.whl", hash = "sha256:dd6a0808a2eb8b5b1ac06bca4244eee18ed2c0e7b105599e1662203d164317b5", size = 676657, upload-time = "2026-08-30T19:51:34.494Z" },
+    { url = "https://files.pythonhosted.org/packages/4b/44/dca7b211804f60c789aced2792b1e7803ccd8b70b79041cbb92788df5d19/hypothesis-6.167.1-cp313-cp313-macosx_10_12_x86_64.whl", hash = "sha256:6478d19a7887731cc2afaa1ec15f62811c9ceb6fd18e5b7563e0a18399a9528f", size = 786947, upload-time = "2026-08-30T19:51:29.165Z" },
+    { url = "https://files.pythonhosted.org/packages/b9/6a/3cffa138492c9e3d5f98f4ff8b467273dc87af6ca3c18084272d106bde10/hypothesis-6.167.1-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:8f13167a4b81c93e7e051d1f02790814a6495fb79cacf3fb89560a796a2f7d00", size = 778584, upload-time = "2026-08-30T19:52:25.091Z" },
+    { url = "https://files.pythonhosted.org/packages/7a/7d/e8039791aaca3b21557bc520a71cdb88751892f66fd1a0a459b59872e463/hypothesis-6.167.1-cp313-cp313-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5ef7dd225f7df7d74d1c5a905592cd8b4cd348e6be639b189a43def8b0b5dd79", size = 1116749, upload-time = "2026-08-30T19:52:38.178Z" },
+    { url = "https://files.pythonhosted.org/packages/ee/b8/f9b8d93bd6178870f0daa868ca99915f6d9df1f99dc7291e9ce2743a6dc5/hypothesis-6.167.1-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7d7585429f2263d3ceeb3474bae3871024630a7a598e71eeb4b0dcf03e291623", size = 1162599, upload-time = "2026-08-30T19:52:11.72Z" },
+    { url = "https://files.pythonhosted.org/packages/a1/0d/53d419094e6f8a7e7377c09de15ac23f842ab698ff07241f7b73e19bd559/hypothesis-6.167.1-cp313-cp313-musllinux_1_2_aarch64.whl", hash = "sha256:fb9194f450417cf35f66b6c72737cc6b8f21f20567ba4d39822c64f0d1075784", size = 1292230, upload-time = "2026-08-30T19:52:49.732Z" },
+    { url = "https://files.pythonhosted.org/packages/5c/df/cf4c482323ae4f06b5326b5bdd89cf17d8232fdb3186c9913e0b19a5fa58/hypothesis-6.167.1-cp313-cp313-musllinux_1_2_x86_64.whl", hash = "sha256:c63a0d292a5dde3c0fe999892e76d8375a003ca40c0c00d763f3360f91be5b96", size = 1328899, upload-time = "2026-08-30T19:51:05.366Z" },
+    { url = "https://files.pythonhosted.org/packages/49/05/780c4b0396491d294fda69a541cb1dedb37fb9eb2e3a696e85fe19064c40/hypothesis-6.167.1-cp313-cp313-win_amd64.whl", hash = "sha256:ff07f98a0b230632bb2836b5dad3e94d85c114ae155a316afd251c58760958ae", size = 675927, upload-time = "2026-08-30T19:50:58.472Z" },
+    { url = "https://files.pythonhosted.org/packages/6a/f1/1e602f090dcb7e38655f1f7909482742891332275fc01f241e255cdfa514/hypothesis-6.167.1-cp314-cp314-macosx_10_12_x86_64.whl", hash = "sha256:fcfc2a78fc1025644f889a74684b3201f4652ce8e6694c2a01af0f100d0348cf", size = 787054, upload-time = "2026-08-30T19:50:40.88Z" },
+    { url = "https://files.pythonhosted.org/packages/52/57/cfa930719c7af5a33627abba826a3fa2efa61a5f23e38d4111eace5dfe53/hypothesis-6.167.1-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:769bdd9aa0af08c063327730ab6dc18b7a23837a2912f2aeaab3912f11a7e3ad", size = 778721, upload-time = "2026-08-30T19:50:36.503Z" },
+    { url = "https://files.pythonhosted.org/packages/01/37/4c4bc3d319eac85bfe17515c9786bf49e57181ca8757886110d2cfb13d10/hypothesis-6.167.1-cp314-cp314-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d75d44bdead6679b6ee9a7c90d10207db865ca0c77c5212103b5ff421379f99e", size = 1116972, upload-time = "2026-08-30T19:51:52.472Z" },
+    { url = "https://files.pythonhosted.org/packages/33/45/0d61ceef2739e7b96ea1faa0f3d5aa5917c8156797993bf3acbadfcd7f0a/hypothesis-6.167.1-cp314-cp314-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:742be00d7bb53d10634e6435e5b98f51fcdbe7ed377d473ab7387d9499c87169", size = 1162776, upload-time = "2026-08-30T19:51:09.084Z" },
+    { url = "https://files.pythonhosted.org/packages/4a/cf/756666ce2262e90fd61fec41a95548cceab94b0669381d8f0387cd89af93/hypothesis-6.167.1-cp314-cp314-musllinux_1_2_aarch64.whl", hash = "sha256:b6fcdc8d03b37a902262be13112d113bb4ac87edf3b08afb47f3d1210deb038a", size = 1292748, upload-time = "2026-08-30T19:51:17.22Z" },
+    { url = "https://files.pythonhosted.org/packages/e8/b4/87eb3c695d6c37fb44f4d49f9faa2033af496e24965658942a1706e22620/hypothesis-6.167.1-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:e56e7841514276c308c2bb4d033cf01860d0fc8c76e2b79ce748a9f123eaf83b", size = 1329101, upload-time = "2026-08-30T19:51:36.313Z" },
+    { url = "https://files.pythonhosted.org/packages/a7/3b/87faa4a86533eaaa19037741fb9cdde8647f7ffdf8fd4279828ac9d81f8b/hypothesis-6.167.1-cp314-cp314-pyemscripten_2026_0_wasm32.whl", hash = "sha256:bbf4f0cad201d0b8e821e82ad828b2aec99ce6d9967779eecb2cad4d4a93debd", size = 618079, upload-time = "2026-08-30T19:52:43.226Z" },
+    { url = "https://files.pythonhosted.org/packages/e0/46/96b7ac9605887447d267b4b3a9ecf61c6caaabf39eef667173b0cc9222b3/hypothesis-6.167.1-cp314-cp314-win_amd64.whl", hash = "sha256:3e04f6001299708b6fd4512267b189c0b029ef1e34500deb4e4c9639023598d7", size = 675812, upload-time = "2026-08-30T19:52:52.298Z" },
+    { url = "https://files.pythonhosted.org/packages/6e/f6/0337a91c50ce4be323c3d6aa852fcf08199ffbb1072da09fbe6d602f4dfe/hypothesis-6.167.1-cp314-cp314t-macosx_10_12_x86_64.whl", hash = "sha256:47c99256df28555ecc2aed0e22ca17cd61c63c8c44207a07b4e402cc49661fae", size = 785525, upload-time = "2026-08-30T19:51:03.637Z" },
+    { url = "https://files.pythonhosted.org/packages/5a/08/9bb52de855169d31888c7033ee2f94b94138fde021c1af9dbc7ba5e83cd5/hypothesis-6.167.1-cp314-cp314t-macosx_11_0_arm64.whl", hash = "sha256:5d6614e88fd267bbd870e3ec02f8a5387897d2d573626a02f5ac06d81533afa6", size = 777142, upload-time = "2026-08-30T19:52:18.472Z" },
+    { url = "https://files.pythonhosted.org/packages/85/79/f1a7e088e13a641357abb9b43d75c116c2a0902711b1a25a203864b96c9b/hypothesis-6.167.1-cp314-cp314t-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:27829aa89fe2e47c5c8d13b3ce31e0f53a01a98f76a4f99cfa3369ab35362f33", size = 1115311, upload-time = "2026-08-30T19:52:40.975Z" },
+    { url = "https://files.pythonhosted.org/packages/e0/38/e28b1fc20bd3d67d43cf1aab7a15daa2a24ec01d17a153e82fcf38c882f3/hypothesis-6.167.1-cp314-cp314t-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e936777d92ae27393b4a941839bbb43c1f339b5a0e394c7f3730454cdf091b3a", size = 1161238, upload-time = "2026-08-30T19:52:20.501Z" },
+    { url = "https://files.pythonhosted.org/packages/5a/de/9b4fc7992166299e0fc5c13c8766919ae57d0fb9eed5319b7a3bad4f2f17/hypothesis-6.167.1-cp314-cp314t-musllinux_1_2_aarch64.whl", hash = "sha256:971ce0d8a367a37c4690b83a2e7f6ef832fa3543357eb6da0da27ba078e5088a", size = 1290974, upload-time = "2026-08-30T19:51:15.673Z" },
+    { url = "https://files.pythonhosted.org/packages/cc/79/ca086eea02588212ab796ee4bd7fe6ed514e10d1a99967e478691608e8d9/hypothesis-6.167.1-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:af84ce2416be2a65bc0ea18e64d2dbb9796b7692593b5b2064d60ea1d52ec1e2", size = 1327969, upload-time = "2026-08-30T19:52:35.846Z" },
+    { url = "https://files.pythonhosted.org/packages/26/79/1875380fa30e8411553e76b3e9695aca0845f9f265523f20f879bdea2b82/hypothesis-6.167.1-cp314-cp314t-win_amd64.whl", hash = "sha256:3b596efec5bd714588e3bb269544d993c5258c979f3a26f51fadf62c215d0e68", size = 675735, upload-time = "2026-08-30T19:52:22.821Z" },
+    { url = "https://files.pythonhosted.org/packages/b0/45/59abecd75e52b9dfb5b3eb991276f54954c44917a1c83d148cfb3580bd39/hypothesis-6.167.1-cp315-abi3.abi3t-macosx_10_12_x86_64.whl", hash = "sha256:c25c556d51d55d94988dc0a2c716d471ff19cf9632cd33f5e6db2914d802428a", size = 785097, upload-time = "2026-08-30T19:51:12.528Z" },
+    { url = "https://files.pythonhosted.org/packages/01/7c/e6d978dc9564ba70352da60c00f55f6ad7d66d99ecbf336a228978206cb4/hypothesis-6.167.1-cp315-abi3.abi3t-macosx_11_0_arm64.whl", hash = "sha256:b57e950f9d5c93ca335bc612e8fa8fb49abb187c3fc9d5e7d9966d52eb27d747", size = 776798, upload-time = "2026-08-30T19:50:47.247Z" },
+    { url = "https://files.pythonhosted.org/packages/d1/a4/f8ecedcf96790aab69d750afe3fcbf503229d0bb4e0c32be655385a4fc8c/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:0807ae8d399162827fc1c396ab4c697a41921c2a2baacfada439771e9dc2b867", size = 1115116, upload-time = "2026-08-30T19:52:16.029Z" },
+    { url = "https://files.pythonhosted.org/packages/7e/97/8bca7c262ac4fcb1ee684c04e4ba75f26541d3a30417e3743d19912d257e/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:546fef39c7aadba74bf3e592585694a71340d1775e9b3274bb3f94106dbde4b7", size = 1137812, upload-time = "2026-08-30T19:51:25.507Z" },
+    { url = "https://files.pythonhosted.org/packages/d6/07/9cb4a7fc2aa2b2ad063b446c378f0d7acfa5303e84afd1b1374ba23fd6f3/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:a17a5618b6a5b84f17c8acb3bce37122647cf7a3e48b660a39d68a773bd627dd", size = 1140384, upload-time = "2026-08-30T19:52:05.264Z" },
+    { url = "https://files.pythonhosted.org/packages/1c/18/813bf7efa18f11ce938a0da22a7518a54db46d4a181ebf4cb0a8061c263f/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:0ff5ad833480c1e34ae902cb52fc00802b08ac6c87bd22d7e5f04fd925869608", size = 1160569, upload-time = "2026-08-30T19:51:40.099Z" },
+    { url = "https://files.pythonhosted.org/packages/52/1d/6658d9294ed33bb17da4acf06fe63b010b205ea1861f6921c38060793255/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_31_riscv64.whl", hash = "sha256:630eb37df80b5bc4ec6f391b13caaecc06942ff5da3aadebacf85f53bbc55757", size = 1120605, upload-time = "2026-08-30T19:53:01.848Z" },
+    { url = "https://files.pythonhosted.org/packages/a6/81/a75821c0e879223a2635b8eded84ae874cb6c711b23e9930668008d0b13f/hypothesis-6.167.1-cp315-abi3.abi3t-manylinux_2_5_i686.manylinux1_i686.whl", hash = "sha256:bc4e65f7c43b187f7a40964706b5ded1073e0c1839e9fb5e041d7ed973bb65fe", size = 1149479, upload-time = "2026-08-30T19:51:30.972Z" },
+    { url = "https://files.pythonhosted.org/packages/4e/f2/c8c3faf4ec796d6dbf36b84662806696434aef38616e5a65b46188c04262/hypothesis-6.167.1-cp315-abi3.abi3t-musllinux_1_2_aarch64.whl", hash = "sha256:e849f518cbc4e76ab15f2f1473c60dd3103da8d32399187325ceb84309105976", size = 1290423, upload-time = "2026-08-30T19:50:51.114Z" },
+    { url = "https://files.pythonhosted.org/packages/a4/e3/e0903abe7e8634cedb5414931452e82daacdd3b8b46d6348bbefcaa45f2a/hypothesis-6.167.1-cp315-abi3.abi3t-musllinux_1_2_armv7l.whl", hash = "sha256:5c5a26d4d3dca0c84e01bde41df4cabaa5a373c7393f9eef372d19fe93b07ccd", size = 1415749, upload-time = "2026-08-30T19:51:48.456Z" },
+    { url = "https://files.pythonhosted.org/packages/a6/9b/03f09c1ecac1dfb0f4cd7fcc6dc50d9c6ea8067b295a728e242650bafe32/hypothesis-6.167.1-cp315-abi3.abi3t-musllinux_1_2_ppc64le.whl", hash = "sha256:4819adbc5911648f6bfaeb574f276add184b4b49f54731dbde46fd71256bb157", size = 1272086, upload-time = "2026-08-30T19:52:29.368Z" },
+    { url = "https://files.pythonhosted.org/packages/7c/7c/7a63bfc0bfbaf000f71352c4faac72ff611376330a2ce2e9a1bf4668848a/hypothesis-6.167.1-cp315-abi3.abi3t-musllinux_1_2_riscv64.whl", hash = "sha256:3ad7206de9c398c8da5745b69b5ba2ef45100082eeb174656490bc4f262b112c", size = 1291553, upload-time = "2026-08-30T19:52:07.545Z" },
+    { url = "https://files.pythonhosted.org/packages/e8/5c/8065bdab53bc81743ca68fc76ca53fc7531a5b3f01c0de4ba40467955d6a/hypothesis-6.167.1-cp315-abi3.abi3t-musllinux_1_2_x86_64.whl", hash = "sha256:96d5e8017a9508f06c8a61a6130cb0d0b4810847ed5c76923cb5cfb9952b31af", size = 1327734, upload-time = "2026-08-30T19:51:42.306Z" },
+    { url = "https://files.pythonhosted.org/packages/0d/a0/5c15d480aea3a8e6e5c17c7cb1170171707ac643ffd319473bc194743ad8/hypothesis-6.167.1-cp315-abi3.abi3t-win32.whl", hash = "sha256:a4e4de36a397cba49d949d89cbc26135977c15f9d797caa95317962ceb5b5674", size = 669115, upload-time = "2026-08-30T19:51:14.192Z" },
+    { url = "https://files.pythonhosted.org/packages/1b/36/4cf494bc96384189fedb7d3f272580315f2284a9f8a7f6a59796612eb76d/hypothesis-6.167.1-cp315-abi3.abi3t-win_amd64.whl", hash = "sha256:f6fe9c40ab14def363d9e7ab22863fa31652bd5e08f8495b34ff7bd0062b3f8d", size = 675438, upload-time = "2026-08-30T19:51:06.881Z" },
+    { url = "https://files.pythonhosted.org/packages/b0/7f/db1a37e5f45be32c0e64f9ed1268eba56aeedcb2ef20d195fa60c6610347/hypothesis-6.167.1-cp315-abi3.abi3t-win_arm64.whl", hash = "sha256:627ce3bd166799a6c0ddcf1351049be5b9a772d5bce436216d42b41a935f42c0", size = 673123, upload-time = "2026-08-30T19:52:13.991Z" },
+]
+
 [[package]]
 name = "idna"
 version = "3.18"
@@ -2068,6 +2133,7 @@ dev = [
     { name = "datashader" },
     { name = "doubletdetection" },
     { name = "hdf5plugin" },
+    { name = "hypothesis" },
     { name = "pytest" },
     { name = "pytest-cov" },
     { name = "ruff" },
@@ -2126,6 +2192,7 @@ dev = [
     { name = "datashader", specifier = ">=0.19" },
     { name = "doubletdetection", specifier = ">=4.3" },
     { name = "hdf5plugin", specifier = ">=7.0.0" },
+    { name = "hypothesis", specifier = ">=6.140" },
     { name = "pytest", specifier = ">=9.0" },
     { name = "pytest-cov", specifier = ">=7.0" },
     { name = "ruff", specifier = ">=0.16" },

From 4db09c7b62f2b5f1a62ace8f7d647196bf3a82ab Mon Sep 17 00:00:00 2001
From: Aaron Meyer <git@ameyer.me>
Date: Fri, 4 Sep 2026 10:27:06 -0700
Subject: [PATCH 2/2] Include weights in energy-based component ordering

---
 docs/component_ordering.md              | 24 +++++++++++++-----------
 scrise/factorization.py                 |  8 ++++----
 scrise/tests/test_component_ordering.py | 14 ++++++++------
 scrise/tests/test_golden_regression.py  |  4 ++--
 scrise/tests/test_invariants.py         |  9 +++++++--
 5 files changed, 34 insertions(+), 25 deletions(-)

diff --git a/docs/component_ordering.md b/docs/component_ordering.md
index e6d1bab1..1a441dea 100644
--- a/docs/component_ordering.md
+++ b/docs/component_ordering.md
@@ -39,19 +39,21 @@ footing (i.e. a fixed sign and a fixed ordering rule).
 RISE orders components by their intrinsic energy,
 
 $$
-e_r = \lVert \mathbf{A}[:, r] \rVert \cdot \lVert \mathbf{C}[:, r] \rVert,
+e_r = |w_r| \cdot \lVert \mathbf{A}[:, r] \rVert \cdot \lVert \mathbf{C}[:, r] \rVert,
 $$
 
-the product of the condition-factor and gene-factor column norms for
-component $r$. This quantity is directly determined by the fit: unlike the
-Gini coefficient, it is not distorted by the arbitrary rescaling that can
-occur between $\mathbf{A}$ and $\mathbf{C}$ during optimization (their
-product is fixed by the fit, but how that product is split between the two
-factors is not). Components are sorted from highest to lowest energy, so
-that low-energy components — typically the ones that appear only once the
-rank is increased — land at the high end of the ordering, while
-established, high-energy components stay near the front. This is
-implemented in [`RISE.order_components_by_energy`][RISE.factorization.order_components_by_energy]
+the product of the component weight, condition-factor column norm, and
+gene-factor column norm for component $r$. (Because standard PARAFAC2 factor
+normalization scales the columns of $\mathbf{A}$, $\mathbf{B}$, and
+$\mathbf{C}$ to unit length, this energy corresponds directly to $|w_r|$ while
+remaining invariant to any alternative scale distribution among the factors.)
+This quantity is directly determined by the fit: unlike the Gini coefficient,
+it is not distorted by arbitrary rescaling between factor matrices. Components
+are sorted from highest to lowest energy, so that low-energy components —
+typically the ones that appear only once the rank is increased — land at the
+high end of the ordering, while established, high-energy components stay near
+the front. This is implemented in
+[`RISE.order_components_by_energy`][RISE.factorization.order_components_by_energy]
 and is applied automatically inside [`RISE.pf2`][RISE.factorization.pf2].
 
 ## Sign convention
diff --git a/scrise/factorization.py b/scrise/factorization.py
index cc40f90e..d6983f78 100644
--- a/scrise/factorization.py
+++ b/scrise/factorization.py
@@ -100,9 +100,9 @@ def order_components_by_energy(X: anndata.AnnData) -> anndata.AnnData:
     components of the new fit correspond to the N components of the old fit.
 
     This function instead orders components by their intrinsic energy,
-    ``||A[:, r]|| * ||C[:, r]||`` (the product of the condition-factor and
-    gene-factor column norms), which is directly determined by the fit and
-    is not subject to the arbitrary rescaling that can occur between A and C.
+    ``|weights[r]| * ||A[:, r]|| * ||C[:, r]||`` (the product of component
+    weights, condition-factor, and gene-factor column norms), which is
+    directly determined by the fit and is not subject to arbitrary rescaling.
     Components are ordered from highest to lowest energy, so that low-energy
     components -- which tend to be the ones added when the rank is
     increased -- land at the high end of the ordering. This is consistent
@@ -142,7 +142,7 @@ def order_components_by_energy(X: anndata.AnnData) -> anndata.AnnData:
     A = A * signs
     C = C * signs
 
-    energy = np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+    energy = np.abs(weights) * np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
     order = np.argsort(energy)[::-1]
 
     X.uns["Pf2_A"] = A[:, order]
diff --git a/scrise/tests/test_component_ordering.py b/scrise/tests/test_component_ordering.py
index 04d8524e..e7233418 100644
--- a/scrise/tests/test_component_ordering.py
+++ b/scrise/tests/test_component_ordering.py
@@ -66,15 +66,17 @@ def test_order_components_by_energy_descending():
     scales = np.array([0.01, 10.0, 1.0, 5.0])
     A = A * scales
     expected_energy_order = np.argsort(
-        np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+        np.abs(weights) * np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
     )[::-1]
 
     adata = _mock_adata(A.copy(), B.copy(), C.copy(), weights.copy(), projections)
     ordered = order_components_by_energy(adata)
 
-    new_energy = np.linalg.norm(
-        np.array(ordered.uns["Pf2_A"]), axis=0
-    ) * np.linalg.norm(np.array(ordered.varm["Pf2_C"]), axis=0)
+    new_energy = (
+        np.abs(np.array(ordered.uns["Pf2_weights"]))
+        * np.linalg.norm(np.array(ordered.uns["Pf2_A"]), axis=0)
+        * np.linalg.norm(np.array(ordered.varm["Pf2_C"]), axis=0)
+    )
     assert np.all(np.diff(new_energy) <= 1e-8)
 
     # Check that the columns were permuted as expected (up to sign).
@@ -109,7 +111,7 @@ def test_order_components_by_energy_preserves_reconstruction():
 
     # Column r of after_wp, weighted by A/C for that component, should match
     # some permuted (and consistently signed) column of before_wp.
-    energy = np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+    energy = np.abs(weights) * np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
     order = np.argsort(energy)[::-1]
 
     # B itself is left as the unflipped sign reference, so weighted_projections
@@ -128,7 +130,7 @@ def test_order_components_by_energy_reorders_weights_and_B():
     weights = np.arange(rank, dtype=float)
     projections, _ = np.linalg.qr(rng.normal(size=(n_cells, rank)))
 
-    energy = np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+    energy = np.abs(weights) * np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
     order = np.argsort(energy)[::-1]
 
     adata = _mock_adata(A.copy(), B.copy(), C.copy(), weights.copy(), projections)
diff --git a/scrise/tests/test_golden_regression.py b/scrise/tests/test_golden_regression.py
index 14621de7..431eb3d3 100644
--- a/scrise/tests/test_golden_regression.py
+++ b/scrise/tests/test_golden_regression.py
@@ -57,10 +57,10 @@ def test_pf2_golden_values_on_fixed_synthetic_fixture():
     # Component weights, order, and which gene dominates each component are
     # a much coarser (and thus more robust-to-numerical-noise) fingerprint
     # of the fit than raw factor values.
-    np.testing.assert_allclose(weights, [22.2080, 56.6194, 24.3079], rtol=1e-2)
+    np.testing.assert_allclose(weights, [56.6194, 24.3079, 22.2080], rtol=1e-2)
 
     top_gene_per_component = np.argmax(np.abs(C), axis=0)
-    np.testing.assert_array_equal(top_gene_per_component, [12, 13, 5])
+    np.testing.assert_array_equal(top_gene_per_component, [13, 5, 12])
 
 
 if __name__ == "__main__":
diff --git a/scrise/tests/test_invariants.py b/scrise/tests/test_invariants.py
index cba20083..2b625580 100644
--- a/scrise/tests/test_invariants.py
+++ b/scrise/tests/test_invariants.py
@@ -121,10 +121,15 @@ def test_order_components_by_energy_is_a_permutation(
 
     new_A = np.array(ordered.uns["Pf2_A"])
     new_C = np.array(ordered.varm["Pf2_C"])
-    new_energy = np.linalg.norm(new_A, axis=0) * np.linalg.norm(new_C, axis=0)
+    new_weights = np.array(ordered.uns["Pf2_weights"])
+    new_energy = (
+        np.abs(new_weights)
+        * np.linalg.norm(new_A, axis=0)
+        * np.linalg.norm(new_C, axis=0)
+    )
     assert np.all(np.diff(new_energy) <= 1e-6)
 
-    energy = np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
+    energy = np.abs(weights) * np.linalg.norm(A, axis=0) * np.linalg.norm(C, axis=0)
     order = np.argsort(energy)[::-1]
     np.testing.assert_allclose(
         np.array(ordered.obsm["weighted_projections"]), before_wp[:, order], atol=1e-5