Add property-based, contract, regression, and plotting tests

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

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/*"]

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]

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()

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)

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()