From 8cd7688080d6e6d78b7b216387e87dfc126e43f5 Mon Sep 17 00:00:00 2001 From: Dominik Date: Fri, 7 Aug 2026 06:56:33 -0700 Subject: [PATCH 01/23] feat(de): add ls_scheme to choose how the shared length scale is estimated The two conditions have always been smoothed at a shared length scale so that their fitted surfaces are comparable. That part is right: giving each condition its own makes the two surfaces differently smooth, and on an exchangeable null the mismatch alone roughly doubles the false-positive rate -- worse than any shared value. What was undocumented is that the shared value is estimated from condition 1's cells only (`ls_for_model2 = self.model1.ls`), so it is a function of n_condition1 and de(adata, condition1=X, condition2=Y) is not equivalent to de(adata, condition1=Y, condition2=X). Measured on a cell-level exchangeable null (2,917 cells, 16,285 genes, truth = 0 differential genes), the two orientations of one partition differ in false-positive rate by 0.059 at 2:1, 0.082 at 3:1 and 0.079 at 6:1, against a floor of ~0.004. `ls_scheme` makes the choice explicit without changing anything by default: "condition1" (default) exactly today's behaviour "symmetric" size-weighted geometric mean of the two conditions' own estimates -- the same estimator, with neighbours looked up within each condition rather than across the union; swap-invariant "pooled" estimated from both conditions' cells together; also swap-invariant, but systematically smaller, since the union is denser than either condition and the estimate tracks nearest-neighbour distance "separate" each condition estimates its own; diagnostics only It is ignored, with a log line, when `ls` or a `cov_func` is supplied. It also joins the run-parameter match list, so re-running under a different scheme is not mistaken for a matching rerun. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 34 ++++ kompot/anndata/_de_helpers.py | 3 + kompot/anndata/differential_expression.py | 6 + .../differential/differential_expression.py | 122 ++++++++++++++- kompot/settings.py | 21 ++- tests/test_ls_scheme.py | 145 ++++++++++++++++++ 6 files changed, 324 insertions(+), 7 deletions(-) create mode 100644 tests/test_ls_scheme.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a336fd0..270c203 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,40 @@ All notable changes to this project will be documented in this file. +## [Unreleased] + +### Added + + - **`GPSettings.ls_scheme`** — chooses which cells the automatic shared length scale is + estimated from in differential expression. The default, `"condition1"`, is exactly the + existing behaviour and nothing changes unless you ask for something else. + + The two conditions have always been smoothed at a **shared** length scale so that their + fitted surfaces are comparable, and that part is right: giving each condition its own + (`ls_scheme="separate"`) makes the two surfaces differently smooth, and the mismatch alone + manufactures apparent fold changes. What was undocumented is that the shared value is + estimated from **condition 1's cells only**, so it is a function of `n_condition1` and + `de(adata, condition1=X, condition2=Y)` is *not* equivalent to + `de(adata, condition1=Y, condition2=X)`. On an exchangeable null with nothing to find, the + two orientations of the same partition differ in false-positive rate by 0.06–0.08 at + 2:1 and beyond. + + `"symmetric"` shares the size-weighted geometric mean of the two conditions' own estimates + — the same estimator, with nearest neighbours looked up *within* each condition — and is + invariant under swapping the conditions. `"pooled"` estimates from the union, which is also + swap-invariant but yields a systematically smaller value, because the union is denser than + either condition and nearest-neighbour distances shrink with cell count alone. + +### Changed + + - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme + is no longer treated as a matching rerun. + +### Documentation + + - `GPSettings.ls` and `DifferentialExpression.fit` now state where the shared length scale + comes from and that the default makes the contrast depend on argument order. + ## [0.8.0] - 2026-07-28 ### Changed — statistics now match the manuscript diff --git a/kompot/anndata/_de_helpers.py b/kompot/anndata/_de_helpers.py index 042c3b0..6548d13 100644 --- a/kompot/anndata/_de_helpers.py +++ b/kompot/anndata/_de_helpers.py @@ -41,6 +41,7 @@ def _check_overwrites( obsm_key: str, layer: Optional[str], ls_factor: float, + ls_scheme: str = "condition1", ): """Check for existing results and handle overwrite logic. @@ -134,6 +135,7 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", + "ls_scheme", ]: curr_val = locals().get(param_name) prev_val = _pg(prev_params, param_name) @@ -189,6 +191,7 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", + "ls_scheme", ]: curr_val = locals().get(param_name) prev_val = _pg(prev_params, param_name) diff --git a/kompot/anndata/differential_expression.py b/kompot/anndata/differential_expression.py index 4315a0d..d0a833b 100644 --- a/kompot/anndata/differential_expression.py +++ b/kompot/anndata/differential_expression.py @@ -132,6 +132,7 @@ def de( sigma = _gp.sigma ls = _gp.ls ls_factor = _gp.ls_factor + ls_scheme = _gp.ls_scheme n_landmarks = _gp.n_landmarks landmarks = _gp.landmarks use_empirical_variance = _gp.use_empirical_variance @@ -305,6 +306,7 @@ def de( obsm_key=obsm_key, layer=layer, ls_factor=ls_factor, + ls_scheme=ls_scheme, ) # ---- 2. Copy if requested ---- @@ -455,6 +457,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, + ls_scheme=ls_scheme, landmarks=landmarks, condition1_sample_indices=condition1_sample_indices, condition2_sample_indices=condition2_sample_indices, @@ -635,6 +638,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, + ls_scheme=ls_scheme, n_landmarks=n_landmarks, use_empirical_variance=use_empirical_variance, batch_size=batch_size, @@ -730,6 +734,7 @@ def compute_differential_expression( sigma: float = 1.0, ls=None, ls_factor: float = 10.0, + ls_scheme: str = "condition1", compute_mahalanobis: bool = True, jit_compile: bool = False, eps: float = 1e-8, @@ -783,6 +788,7 @@ def compute_differential_expression( sigma=sigma, ls=ls, ls_factor=ls_factor, + ls_scheme=ls_scheme, n_landmarks=n_landmarks, landmarks=landmarks, use_empirical_variance=use_empirical_variance, diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 6af2bee..5bfdb9b 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -14,6 +14,61 @@ logger = logging.getLogger("kompot") +#: Recognised values for ``ls_scheme``. See :meth:`DifferentialExpression.fit`. +LS_SCHEMES = ("condition1", "symmetric", "pooled", "separate") + + +def _auto_ls(X: np.ndarray, ls_factor: float) -> float: + """Automatic length scale for one design matrix. + + Identical to what :meth:`ExpressionModel.fit` derives internally when + ``ls`` is None, exposed here so a scheme can combine the two conditions' + values *before* either model is fitted. + """ + from mellon.parameters import compute_ls, compute_nn_distances + + return float(compute_ls(compute_nn_distances(np.asarray(X))) * ls_factor) + + +def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor): + """Return the ``(ls_model1, ls_model2)`` a scheme prescribes. + + ``None`` in a slot means "let that model estimate its own". For + ``"condition1"`` both slots are None and the caller reproduces the + historical behaviour by copying model1's fitted value into model2. + """ + if ls_scheme == "condition1": + return None, None + if ls_scheme == "separate": + return None, None + if ls_scheme == "pooled": + shared = _auto_ls( + np.vstack([np.asarray(X_condition1), np.asarray(X_condition2)]), + ls_factor, + ) + logger.info(f"ls_scheme='pooled': shared length scale {shared:.4f}") + return shared, shared + if ls_scheme == "symmetric": + n1 = np.asarray(X_condition1).shape[0] + n2 = np.asarray(X_condition2).shape[0] + ls1 = _auto_ls(X_condition1, ls_factor) + ls2 = _auto_ls(X_condition2, ls_factor) + # Size-weighted GEOMETRIC mean. kompot's estimator is + # exp(mean(log nn_dist)) * const, so this is exactly that estimator + # applied to the two conditions' WITHIN-condition nearest-neighbour + # distances pooled together -- unlike 'pooled', which looks the + # neighbours up across the union and therefore shrinks purely because + # the union holds more cells. + shared = float(np.exp((n1 * np.log(ls1) + n2 * np.log(ls2)) / (n1 + n2))) + logger.info( + f"ls_scheme='symmetric': condition length scales {ls1:.4f} / " + f"{ls2:.4f} -> shared {shared:.4f}" + ) + return shared, shared + raise ValueError( + f"Unknown ls_scheme {ls_scheme!r}. Expected one of {LS_SCHEMES}." + ) + class DifferentialExpression: """ @@ -310,6 +365,7 @@ def fit( sigma: float = 1.0, ls: Optional[float] = None, ls_factor: float = 10.0, + ls_scheme: str = "condition1", landmarks: Optional[np.ndarray] = None, sample_estimator_ls: Optional[float] = None, condition1_sample_indices: Optional[np.ndarray] = None, @@ -336,10 +392,38 @@ def fit( sigma : float, optional Noise level for function estimator, by default 1.0. ls : float, optional - Length scale for the GP kernel. If None, it will be estimated, by default None. + Length scale for the GP kernel, shared by both conditions. If None it is + estimated from the data according to ``ls_scheme``, by default None. ls_factor : float, optional Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is None. + ls_scheme : str, optional + How the automatic length scale is derived when ``ls`` is None. Both + conditions are normally smoothed at the *same* scale so that their + fitted surfaces are comparable; the schemes differ in which cells the + shared value is estimated from. + + * ``"condition1"`` (default) — estimate from condition 1's cells and + reuse the value for condition 2. **Not symmetric**: the length + scale is a function of ``n_condition1`` alone, so swapping + ``condition1`` and ``condition2`` changes the result even when + nothing else does. + * ``"symmetric"`` — estimate a length scale from each condition + separately and share their size-weighted geometric mean. Invariant + under swapping the two conditions, and it does not inherit the + cell-count artefact of ``"pooled"``. + * ``"pooled"`` — estimate from the two conditions' cells taken + together. Also swap-invariant, but the union is denser than either + condition, so nearest-neighbour distances — and hence the length + scale — shrink purely because there are more cells. + * ``"separate"`` — let each condition estimate its own length scale + and do not share. The two surfaces are then smoothed differently + and the mismatch itself produces apparent fold changes; measured on + an exchangeable null this is markedly worse than any shared value. + Provided for diagnostics, not recommended. + + Ignored when ``ls`` is given explicitly or when ``cov_func`` / + ``cov_func_curry`` is supplied in ``function_kwargs``. landmarks : np.ndarray, optional Pre-computed landmarks to use. If provided, n_landmarks will be ignored. Shape (n_landmarks, n_features). @@ -420,13 +504,35 @@ def fit( disk_storage_dir=self.disk_storage_dir, ) + # -- Resolve the length scale each condition will be fitted with -- + # Only relevant when the caller left `ls` open and did not hand mellon a + # kernel directly; otherwise the explicit value applies to both models + # and the scheme has nothing to decide. + ls_auto = ( + ls is None + and "ls" not in function_kwargs + and "cov_func" not in function_kwargs + and "cov_func_curry" not in function_kwargs + ) + if ls_auto: + ls_model1, ls_model2 = _resolve_ls_scheme( + ls_scheme, X_condition1, X_condition2, ls_factor + ) + else: + if ls_scheme != "condition1": + logger.info( + f"ls_scheme={ls_scheme!r} ignored: the length scale is " + "already determined by an explicit ls or cov_func." + ) + ls_model1 = ls_model2 = ls + if self.model1.predictor is None: logger.info("Fitting expression estimator for condition 1...") self.model1.fit( X_condition1, y_condition1, sigma=sigma, - ls=ls, + ls=ls_model1, ls_factor=ls_factor, landmarks=landmarks, sample_indices=condition1_sample_indices @@ -437,9 +543,15 @@ def fit( **function_kwargs, ) - # Extract ls from model1 for model2 consistency - ls_for_model2 = ls - if ls is None and "ls" not in function_kwargs and self.model1.ls is not None: + # Under the default scheme condition 2 inherits condition 1's fitted + # length scale, so the shared value is a function of n_condition1 alone. + # Every other scheme has already fixed both values above. + ls_for_model2 = ls_model2 + if ( + ls_auto + and ls_scheme == "condition1" + and self.model1.ls is not None + ): ls_for_model2 = self.model1.ls # -- Fit model2 -- diff --git a/kompot/settings.py b/kompot/settings.py index 1f80ca1..3bc062a 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -38,10 +38,19 @@ class GPSettings: sigma : float Noise level for the GP. ls : float, optional - Length scale. If *None*, estimated automatically using - ``ls_factor``. + Length scale, shared by both conditions. If *None*, estimated + automatically using ``ls_factor`` and ``ls_scheme``. ls_factor : float Multiplier applied to the automatically inferred length scale. + ls_scheme : str + Which cells the automatic shared length scale is estimated from. + ``"condition1"`` (default) uses condition 1's cells only, which makes + the result depend on which condition is passed first; ``"symmetric"`` + shares the size-weighted geometric mean of the two per-condition + estimates; ``"pooled"`` estimates from both conditions' cells taken + together; ``"separate"`` gives each condition its own. See + :meth:`kompot.differential.DifferentialExpression.fit`. Ignored when + ``ls`` is given explicitly. n_landmarks : int, optional Number of landmarks for the Nystrom approximation. landmarks : np.ndarray, optional @@ -61,6 +70,7 @@ class GPSettings: sigma: float = 1.0 ls: Optional[float] = None ls_factor: float = 10.0 + ls_scheme: str = "condition1" n_landmarks: Optional[int] = 5000 landmarks: Optional[np.ndarray] = None use_empirical_variance: bool = False @@ -73,6 +83,13 @@ def __post_init__(self): validate_positive_float(self.sigma, "sigma") validate_positive_float(self.ls, "ls", optional=True) validate_positive_float(self.ls_factor, "ls_factor") + from .differential.differential_expression import LS_SCHEMES + + if self.ls_scheme not in LS_SCHEMES: + raise ValueError( + f"'ls_scheme' must be one of {LS_SCHEMES} (got " + f"{self.ls_scheme!r})." + ) validate_positive_int(self.n_landmarks, "n_landmarks", optional=True) validate_bool(self.use_empirical_variance, "use_empirical_variance") validate_positive_int(self.batch_size, "batch_size", optional=True) diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py new file mode 100644 index 0000000..2d327f4 --- /dev/null +++ b/tests/test_ls_scheme.py @@ -0,0 +1,145 @@ +"""Length-scale determination schemes for the two-condition expression fit. + +The default (``ls_scheme="condition1"``) estimates the shared length scale from +condition 1's cells only, which makes the contrast depend on which condition is +passed first. These tests pin the default's behaviour and check that each +alternative does what it says. +""" + +import numpy as np +import pytest + +from kompot.differential.differential_expression import ( + LS_SCHEMES, + DifferentialExpression, + _auto_ls, + _resolve_ls_scheme, +) +from kompot.settings import GPSettings + + +def _data(n1=120, n2=40, n_genes=6, seed=0): + rng = np.random.default_rng(seed) + X = rng.normal(size=(n1 + n2, 3)) + y = rng.normal(size=(n1 + n2, n_genes)) + return X[:n1], y[:n1], X[n1:], y[n1:] + + +# -------------------------------------------------------------------------- +# scheme resolution +# -------------------------------------------------------------------------- + + +def test_condition1_and_separate_defer_to_the_models(): + X1, _, X2, _ = _data() + assert _resolve_ls_scheme("condition1", X1, X2, 10.0) == (None, None) + assert _resolve_ls_scheme("separate", X1, X2, 10.0) == (None, None) + + +def test_pooled_uses_the_union_and_is_smaller_than_either_condition(): + X1, _, X2, _ = _data() + p1, p2 = _resolve_ls_scheme("pooled", X1, X2, 10.0) + assert p1 == p2 + assert p1 == pytest.approx(_auto_ls(np.vstack([X1, X2]), 10.0)) + # The union is denser than either part, so nearest-neighbour distances -- + # and hence the length scale -- shrink. This is the property that makes + # "pooled" a different question from "what scale does this condition + # support", and it is why it is not the default. + assert p1 < _auto_ls(X1, 10.0) + assert p1 < _auto_ls(X2, 10.0) + + +def test_symmetric_is_the_size_weighted_geometric_mean(): + X1, _, X2, _ = _data() + got, got2 = _resolve_ls_scheme("symmetric", X1, X2, 10.0) + assert got == got2 + a, b = _auto_ls(X1, 10.0), _auto_ls(X2, 10.0) + n1, n2 = X1.shape[0], X2.shape[0] + assert got == pytest.approx(np.exp((n1 * np.log(a) + n2 * np.log(b)) / (n1 + n2))) + assert min(a, b) <= got <= max(a, b) + + +def test_symmetric_is_invariant_under_swapping_the_conditions(): + X1, _, X2, _ = _data() + fwd, _ = _resolve_ls_scheme("symmetric", X1, X2, 10.0) + rev, _ = _resolve_ls_scheme("symmetric", X2, X1, 10.0) + assert fwd == pytest.approx(rev, rel=1e-9) + + +def test_pooled_is_invariant_under_swapping_the_conditions(): + X1, _, X2, _ = _data() + fwd, _ = _resolve_ls_scheme("pooled", X1, X2, 10.0) + rev, _ = _resolve_ls_scheme("pooled", X2, X1, 10.0) + # The union is the same point set either way; the only difference is row + # order, which the approximate nearest-neighbour index is mildly + # sensitive to. + assert fwd == pytest.approx(rev, rel=1e-2) + + +def test_unknown_scheme_is_rejected(): + X1, _, X2, _ = _data() + with pytest.raises(ValueError, match="Unknown ls_scheme"): + _resolve_ls_scheme("nope", X1, X2, 10.0) + + +def test_gp_settings_rejects_unknown_scheme(): + with pytest.raises(ValueError, match="ls_scheme"): + GPSettings(ls_scheme="nope") + for scheme in LS_SCHEMES: + GPSettings(ls_scheme=scheme) + + +def test_gp_settings_default_preserves_historical_behaviour(): + assert GPSettings().ls_scheme == "condition1" + + +# -------------------------------------------------------------------------- +# end-to-end through DifferentialExpression.fit +# -------------------------------------------------------------------------- + + +def _fit(scheme, swap=False, **kw): + X1, y1, X2, y2 = _data() + if swap: + X1, y1, X2, y2 = X2, y2, X1, y1 + de = DifferentialExpression(n_landmarks=0) + de.fit(X1, y1, X2, y2, ls_scheme=scheme, **kw) + return float(de.model1.ls), float(de.model2.ls) + + +def test_default_gives_condition2_condition1s_length_scale(): + ls1, ls2 = _fit("condition1") + assert ls1 == pytest.approx(ls2) + X1, _, _, _ = _data() + assert ls1 == pytest.approx(_auto_ls(X1, 10.0), rel=1e-6) + + +def test_default_is_not_invariant_under_swapping_the_conditions(): + """The defect #310 reports: de(A, B) and de(B, A) smooth differently.""" + fwd, _ = _fit("condition1") + rev, _ = _fit("condition1", swap=True) + assert fwd != pytest.approx(rev, rel=1e-3) + + +def test_separate_gives_each_condition_its_own(): + X1, _, X2, _ = _data() + ls1, ls2 = _fit("separate") + assert ls1 == pytest.approx(_auto_ls(X1, 10.0), rel=1e-6) + assert ls2 == pytest.approx(_auto_ls(X2, 10.0), rel=1e-6) + assert ls1 != pytest.approx(ls2, rel=1e-3) + + +@pytest.mark.parametrize("scheme", ["symmetric", "pooled"]) +def test_shared_schemes_are_swap_invariant_end_to_end(scheme): + fwd = _fit(scheme) + rev = _fit(scheme, swap=True) + assert fwd[0] == pytest.approx(fwd[1]) + assert rev[0] == pytest.approx(rev[1]) + assert fwd[0] == pytest.approx(rev[0], rel=1e-2) + + +@pytest.mark.parametrize("scheme", LS_SCHEMES) +def test_explicit_ls_overrides_every_scheme(scheme): + ls1, ls2 = _fit(scheme, ls=0.5) + assert ls1 == pytest.approx(0.5) + assert ls2 == pytest.approx(0.5) From 487a02e984fb374743b96cf59d4cab133e424aac Mon Sep 17 00:00:00 2001 From: Dominik Date: Fri, 7 Aug 2026 07:43:48 -0700 Subject: [PATCH 02/23] fix(de): a custom cov_func must not turn the default into `separate` The `ls_scheme` gate excluded `cov_func` / `cov_func_curry` from the automatic path. Under the default scheme that dropped the historical `ls_for_model2 = self.model1.ls` inheritance, so a caller supplying a custom kernel and never touching `ls_scheme` silently got `separate` behaviour -- the worst-measured scheme in the comparison this branch is based on (tail mass 0.322 vs 0.172 at 2:1 on an exchangeable null). `function_kwargs` is public and the CLI forwards into it, so the path is reachable from the shipped API. The gate is now exactly the historical one -- `ls is None and "ls" not in function_kwargs`. A supplied kernel no longer disables scheme resolution: the resolved value is handed to the custom kernel via `estimator_defaults["ls"]`, which is what the old inheritance did too, so both conditions keep a shared scale. Verified against 19ee1a1 in a clean worktree, default scheme, `ls=None`: no function_kwargs old 92.589588/92.589588 new 92.589588/92.589588 cov_func_curry=Matern52 old 92.589588/92.589588 new 92.589588/92.589588 Three tests added covering the curry path for the default, the shared schemes and `separate`. Caught by an independent skeptic review, which also noted that 17/17 green missed it because no test exercised `function_kwargs`. --- .../differential/differential_expression.py | 24 +++++----- tests/test_ls_scheme.py | 44 +++++++++++++++++++ 2 files changed, 56 insertions(+), 12 deletions(-) diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 5bfdb9b..8624269 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -422,8 +422,10 @@ def fit( an exchangeable null this is markedly worse than any shared value. Provided for diagnostics, not recommended. - Ignored when ``ls`` is given explicitly or when ``cov_func`` / - ``cov_func_curry`` is supplied in ``function_kwargs``. + Ignored when ``ls`` is given explicitly, or when a length scale is + passed through ``function_kwargs``. Supplying a ``cov_func`` or + ``cov_func_curry`` does *not* disable it: the resolved value is + handed to the custom kernel, so both conditions keep a shared scale. landmarks : np.ndarray, optional Pre-computed landmarks to use. If provided, n_landmarks will be ignored. Shape (n_landmarks, n_features). @@ -505,15 +507,13 @@ def fit( ) # -- Resolve the length scale each condition will be fitted with -- - # Only relevant when the caller left `ls` open and did not hand mellon a - # kernel directly; otherwise the explicit value applies to both models - # and the scheme has nothing to decide. - ls_auto = ( - ls is None - and "ls" not in function_kwargs - and "cov_func" not in function_kwargs - and "cov_func_curry" not in function_kwargs - ) + # The gate is deliberately the historical one: `ls` left open and no + # length scale smuggled in through function_kwargs. A supplied + # `cov_func`/`cov_func_curry` must NOT disable it -- the two conditions + # still need a shared scale, and gating on the kernel would drop the + # condition-1 inheritance for callers who never touched `ls_scheme`, + # silently giving them "separate" behaviour. + ls_auto = ls is None and "ls" not in function_kwargs if ls_auto: ls_model1, ls_model2 = _resolve_ls_scheme( ls_scheme, X_condition1, X_condition2, ls_factor @@ -522,7 +522,7 @@ def fit( if ls_scheme != "condition1": logger.info( f"ls_scheme={ls_scheme!r} ignored: the length scale is " - "already determined by an explicit ls or cov_func." + "already fixed by an explicit ls." ) ls_model1 = ls_model2 = ls diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py index 2d327f4..2d3a2ad 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_ls_scheme.py @@ -143,3 +143,47 @@ def test_explicit_ls_overrides_every_scheme(scheme): ls1, ls2 = _fit(scheme, ls=0.5) assert ls1 == pytest.approx(0.5) assert ls2 == pytest.approx(0.5) + + +# -------------------------------------------------------------------------- +# a custom kernel must not silently change the default's sharing behaviour +# -------------------------------------------------------------------------- + + +def _fit_with_curry(scheme): + from mellon.cov import Matern52 + + X1, y1, X2, y2 = _data() + de = DifferentialExpression(n_landmarks=0) + de.fit(X1, y1, X2, y2, ls_scheme=scheme, cov_func_curry=Matern52) + return float(de.model1.ls), float(de.model2.ls) + + +def test_default_still_shares_when_a_cov_func_curry_is_supplied(): + """Regression: `cov_func_curry` must not turn the default into `separate`. + + `function_kwargs` is public and the CLI forwards into it, so a caller who + never touched `ls_scheme` can reach this path. Before the fix, supplying a + curry dropped the condition-1 inheritance and let condition 2 estimate its + own -- i.e. `separate` behaviour, silently, at the default. + """ + ls1, ls2 = _fit_with_curry("condition1") + assert ls1 == pytest.approx(ls2), ( + "condition 2 stopped inheriting condition 1's length scale when a " + "cov_func_curry was supplied" + ) + # and it is condition 1's own estimate that both are using + plain1, plain2 = _fit("condition1") + assert ls1 == pytest.approx(plain1, rel=1e-6) + + +def test_shared_schemes_still_share_when_a_cov_func_curry_is_supplied(): + for scheme in ("symmetric", "pooled"): + ls1, ls2 = _fit_with_curry(scheme) + assert ls1 == pytest.approx(ls2), scheme + assert ls1 == pytest.approx(_fit(scheme)[0], rel=1e-6), scheme + + +def test_separate_still_unshares_when_a_cov_func_curry_is_supplied(): + ls1, ls2 = _fit_with_curry("separate") + assert ls1 != pytest.approx(ls2, rel=1e-3) From 3cfc90417caa204f772ac53fca4ede0c3317e168 Mon Sep 17 00:00:00 2001 From: Dominik Date: Mon, 21 Sep 2026 16:57:38 -0700 Subject: [PATCH 03/23] feat(de): add the "condition2" ls_scheme, the mirror of the default "condition1" estimates the shared length scale from condition 1's cells and gives it to model 2. "condition2" does the reverse. It is not a fifth alternative -- it is as asymmetric as the default, deliberately -- and it earns its place because of an exact equivalence: de(X, Y, ls_scheme="condition1") == de(Y, X, ls_scheme="condition2") as the same computation with the labels exchanged. That makes the default's swap-dependence demonstrable from a single call site, without swapping condition1/condition2 at the call site and re-deriving which sign is which. Measured, not argued (120+40 cells, seed 0, identical landmarks; max_rel = max |u-v| / max(|u|,|v|)): condition1_smoothed <-> condition2_smoothed 0.000e+00 condition1_std <-> condition2_std 0.000e+00 mahalanobis_distances 0.000e+00 neg_log10_ptp 0.000e+00 fold_change vs NEGATED mirror 0.000e+00 mean_log_fold_change vs NEGATED mirror 0.000e+00 fold_change_zscores vs NEGATED mirror 0.000e+00 control: unflipped fold_change max_rel = 2.000e+00 (the signature of u = -v), max|fold_change| = 0.5054, so the arrays are not near zero. The signed quantities flip sign because fold_change is condition2 - condition1. The unsigned ones do not, because the fold change enters the Mahalanobis form quadratically over a cov1 + cov2 denominator that is symmetric in the two models. WHAT DOES NOT MIRROR, and it is not about ls_scheme: fit() derives automatic landmarks from np.vstack([X_condition1, X_condition2]), whose ROW ORDER differs between the two orientations. compute_landmarks is order-sensitive (measured: different sets at random_state=42, and not even run-to-run stable at random_state=None, the DifferentialExpression default). With automatic landmarks that is the ONLY remaining disagreement between the orientations -- max_rel 1.5e-05 on the smoothed surfaces, 1.8e-06 on the Mahalanobis distances -- so the equivalence is stated as "given identical landmarks". A pre-existing, scheme-independent swap asymmetry that no ls_scheme removes. Implementation note. "condition1" resolves lazily -- model1 estimates internally, model2 inherits self.model1.ls -- and "condition2" cannot, because model1 is fitted first. It resolves eagerly through _auto_ls, so the mirror holds only while the eager and lazy estimates agree exactly. They do, to the last bit, and that is now pinned by a test rather than assumed. The cov_func gate that 487a02e fixed is value-blind (ls is None and "ls" not in function_kwargs), so condition2 inherits the fix rather than needing its own -- now a measurement, not a reading. Also fixed: two readers of ls_scheme that could not see it. - params_get could not read it back. build_params_dict stores it nested at params["gp"]["ls_scheme"], _check_overwrites compares it through params_get, and params_get resolves nested keys through _LEGACY_MAP, which did not list it. So prev_val was always None and params_match was False for EVERY rerun, including a byte-identical one. - The DE CLI silently misrouted a flat ls_scheme config key. An unrecognised key is not rejected; it falls through to de(**function_kwargs) and is handed to mellon. gp_keys is now the module-level GP_CONFIG_KEYS, asserted against GPSettings' own fields, so a future field added without a CLI route fails a test instead of going quietly into mellon. tests/test_ls_scheme.py: 20 -> 33 tests. Negative control: planting _auto_ls(X_condition2) -> _auto_ls(X_condition1) flips exactly the 7 mirror tests and correctly leaves the other two condition2 tests green. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 19 ++ docs/source/cli.rst | 2 +- kompot/anndata/utils/params.py | 1 + kompot/cli/de.py | 31 ++- kompot/cli/templates/de_config_template.yaml | 2 + .../differential/differential_expression.py | 22 +- kompot/settings.py | 12 +- tests/test_ls_scheme.py | 248 +++++++++++++++++- 8 files changed, 318 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 270c203..124e7b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,15 @@ All notable changes to this project will be documented in this file. two orientations of the same partition differ in false-positive rate by 0.06–0.08 at 2:1 and beyond. + `"condition2"` is the default's mirror — it estimates from condition 2's cells and hands the + value to condition 1. It is not an alternative to the default but a **diagnostic**: + `de(X, Y, ls_scheme="condition1")` and `de(Y, X, ls_scheme="condition2")` are the same + computation with the labels exchanged (measured bit-identical for the smoothed surfaces, the + posterior standard deviations, the Mahalanobis distances and `neg_log10_ptp`; the fold change + and its z-scores flip sign, as the definition `condition2 - condition1` requires). So the + default's swap-dependence can be exhibited from a single call site, without swapping + `condition1` and `condition2` and re-deriving which sign is which. + `"symmetric"` shares the size-weighted geometric mean of the two conditions' own estimates — the same estimator, with nearest neighbours looked up *within* each condition — and is invariant under swapping the conditions. `"pooled"` estimates from the union, which is also @@ -31,6 +40,16 @@ All notable changes to this project will be documented in this file. - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme is no longer treated as a matching rerun. +### Fixed + + - `ls_scheme` is now readable back out of a stored `run_info` params dict. It is written nested + under `params["gp"]`, but `params_get` resolves nested keys through `_LEGACY_MAP` and the + field was missing from it — so it read back as `None` and *every* rerun, including an + identical one, compared unequal and was reported as a parameter change. + - The DE CLI now routes a flat `ls_scheme` config key into `GPSettings`. An unrecognised key + is not rejected: it falls through to `**function_kwargs` and is handed to mellon, so the + documented knob was unreachable from a config file and the omission was silent. + ### Documentation - `GPSettings.ls` and `DifferentialExpression.fit` now state where the shared length scale diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 36dc71a..975f3b8 100644 --- a/docs/source/cli.rst +++ b/docs/source/cli.rst @@ -443,7 +443,7 @@ The table below shows how config keys correspond to Python Settings: * - Config key - Python equivalent - Description - * - ``sigma``, ``ls``, ``ls_factor``, ``n_landmarks``, ``batch_size``, ``eps``, ``jit_compile``, ``random_state``, ``use_empirical_variance`` + * - ``sigma``, ``ls``, ``ls_factor``, ``ls_scheme`` (DE only), ``n_landmarks``, ``batch_size``, ``eps``, ``jit_compile``, ``random_state``, ``use_empirical_variance`` - ``gp=GPSettings(...)`` - GP model parameters * - ``null_genes``, ``null_seed``, ``fdr_threshold`` diff --git a/kompot/anndata/utils/params.py b/kompot/anndata/utils/params.py index 4a46042..6732f14 100644 --- a/kompot/anndata/utils/params.py +++ b/kompot/anndata/utils/params.py @@ -197,6 +197,7 @@ def params_to_call_args( "sigma": ("gp", "sigma"), "ls": ("gp", "ls"), "ls_factor": ("gp", "ls_factor"), + "ls_scheme": ("gp", "ls_scheme"), "n_landmarks": ("gp", "n_landmarks"), "use_empirical_variance": ("gp", "use_empirical_variance"), "batch_size": ("gp", "batch_size"), diff --git a/kompot/cli/de.py b/kompot/cli/de.py index 5e6f80c..9a865ae 100644 --- a/kompot/cli/de.py +++ b/kompot/cli/de.py @@ -205,6 +205,24 @@ def add_de_parser(subparsers) -> argparse.ArgumentParser: return parser +#: Flat config keys the DE CLI routes into :class:`~kompot.settings.GPSettings`. +#: Module level so a test can assert it against the dataclass's own fields -- +#: a key missing here is not rejected, it falls through to ``**function_kwargs`` +#: and is handed to mellon, so the omission is silent at the config surface. +GP_CONFIG_KEYS = { + "sigma", + "ls", + "ls_factor", + "ls_scheme", + "n_landmarks", + "use_empirical_variance", + "batch_size", + "eps", + "jit_compile", + "random_state", +} + + def run_de(args): """ Run differential expression analysis. @@ -317,18 +335,7 @@ def run_de(args): sample_col = params.pop("sample_col", None) # Build Settings from remaining params - gp_keys = { - "sigma", - "ls", - "ls_factor", - "n_landmarks", - "use_empirical_variance", - "batch_size", - "eps", - "jit_compile", - "random_state", - } - gp_kwargs = {k: params.pop(k) for k in list(params) if k in gp_keys} + gp_kwargs = {k: params.pop(k) for k in list(params) if k in GP_CONFIG_KEYS} gp = GPSettings(**gp_kwargs) if gp_kwargs else None fdr_kwargs = {} diff --git a/kompot/cli/templates/de_config_template.yaml b/kompot/cli/templates/de_config_template.yaml index 80abd92..c2d1f18 100644 --- a/kompot/cli/templates/de_config_template.yaml +++ b/kompot/cli/templates/de_config_template.yaml @@ -43,6 +43,8 @@ jit_compile: false # Use JAX JIT compilation sigma: 1.0 # Noise level for function estimator ls: null # Length scale (null = auto-estimate) ls_factor: 10.0 # Length scale multiplication factor +ls_scheme: condition1 # Which cells the shared length scale is estimated from: + # condition1 | condition2 | symmetric | pooled | separate eps: 1.0e-8 # Numerical stability constant # Memory management: diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 8624269..76b5bd8 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -15,7 +15,7 @@ logger = logging.getLogger("kompot") #: Recognised values for ``ls_scheme``. See :meth:`DifferentialExpression.fit`. -LS_SCHEMES = ("condition1", "symmetric", "pooled", "separate") +LS_SCHEMES = ("condition1", "condition2", "symmetric", "pooled", "separate") def _auto_ls(X: np.ndarray, ls_factor: float) -> float: @@ -39,6 +39,17 @@ def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor): """ if ls_scheme == "condition1": return None, None + if ls_scheme == "condition2": + # The mirror of "condition1": estimate from condition 2's cells and give + # the value to condition 1. Resolved eagerly rather than by inheritance + # because model1 is fitted first; `_auto_ls` reproduces what a model + # derives internally bit for bit, which is what makes + # de(X, Y, "condition1") and de(Y, X, "condition2") the same computation + # with the labels exchanged. Pinned by + # test_auto_ls_is_bit_identical_to_the_models_internal_estimate. + shared = _auto_ls(X_condition2, ls_factor) + logger.info(f"ls_scheme='condition2': shared length scale {shared:.4f}") + return shared, shared if ls_scheme == "separate": return None, None if ls_scheme == "pooled": @@ -408,6 +419,15 @@ def fit( scale is a function of ``n_condition1`` alone, so swapping ``condition1`` and ``condition2`` changes the result even when nothing else does. + * ``"condition2"`` — the mirror of the default: estimate from + condition 2's cells and reuse the value for condition 1. Equally + asymmetric, and deliberately so. Its use is **diagnostic**: + ``de(X, Y, ls_scheme="condition1")`` and + ``de(Y, X, ls_scheme="condition2")`` are the same computation with + the labels exchanged, so running both orientations of a contrast + under the two schemes exposes the default's swap-dependence from a + single call site, without having to swap ``condition1`` and + ``condition2`` at the call site and re-derive which sign is which. * ``"symmetric"`` — estimate a length scale from each condition separately and share their size-weighted geometric mean. Invariant under swapping the two conditions, and it does not inherit the diff --git a/kompot/settings.py b/kompot/settings.py index 3bc062a..a0e2236 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -45,10 +45,14 @@ class GPSettings: ls_scheme : str Which cells the automatic shared length scale is estimated from. ``"condition1"`` (default) uses condition 1's cells only, which makes - the result depend on which condition is passed first; ``"symmetric"`` - shares the size-weighted geometric mean of the two per-condition - estimates; ``"pooled"`` estimates from both conditions' cells taken - together; ``"separate"`` gives each condition its own. See + the result depend on which condition is passed first; ``"condition2"`` + is its mirror and exists to make that dependence visible in a single + run, since ``de(X, Y, ls_scheme="condition1")`` and + ``de(Y, X, ls_scheme="condition2")`` are the same computation with the + labels exchanged; ``"symmetric"`` shares the size-weighted geometric + mean of the two per-condition estimates; ``"pooled"`` estimates from + both conditions' cells taken together; ``"separate"`` gives each + condition its own. See :meth:`kompot.differential.DifferentialExpression.fit`. Ignored when ``ls`` is given explicitly. n_landmarks : int, optional diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py index 2d3a2ad..8aa1359 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_ls_scheme.py @@ -3,7 +3,10 @@ The default (``ls_scheme="condition1"``) estimates the shared length scale from condition 1's cells only, which makes the contrast depend on which condition is passed first. These tests pin the default's behaviour and check that each -alternative does what it says. +alternative does what it says -- including ``"condition2"``, the default's +mirror, for which the statement is an exact equivalence: +``de(X, Y, ls_scheme="condition1")`` and ``de(Y, X, ls_scheme="condition2")`` +are the same computation with the labels exchanged. """ import numpy as np @@ -187,3 +190,246 @@ def test_shared_schemes_still_share_when_a_cov_func_curry_is_supplied(): def test_separate_still_unshares_when_a_cov_func_curry_is_supplied(): ls1, ls2 = _fit_with_curry("separate") assert ls1 != pytest.approx(ls2, rel=1e-3) + + +# -------------------------------------------------------------------------- +# "condition2" -- the mirror of the default +# -------------------------------------------------------------------------- + + +def test_auto_ls_is_bit_identical_to_the_models_internal_estimate(): + """`_auto_ls` must reproduce what `ExpressionModel.fit` derives for itself. + + This is the load-bearing assumption behind the `condition1` / `condition2` + mirror. `condition1` resolves lazily -- model1 estimates internally and + model2 inherits the fitted value -- while `condition2` has to resolve + *eagerly*, because model1 is fitted first and cannot inherit from a model + that does not exist yet. The two schemes are the same computation only for + as long as the eager and lazy estimates agree exactly. Measured equal to + the last bit on numpy 2.4 / jax 0.10 / mellon 1.7.1. + """ + X1, y1, X2, y2 = _data() + de = DifferentialExpression(n_landmarks=0) + de.fit(X1, y1, X2, y2, ls_scheme="condition1") + assert float(de.model1.ls) == _auto_ls(X1, 10.0) + + +def test_condition2_shares_condition2s_length_scale(): + X1, _, X2, _ = _data() + a, b = _resolve_ls_scheme("condition2", X1, X2, 10.0) + assert a == b + assert a == pytest.approx(_auto_ls(X2, 10.0)) + + +def test_condition2_end_to_end_gives_condition1_condition2s_length_scale(): + X1, _, X2, _ = _data() + ls1, ls2 = _fit("condition2") + assert ls1 == pytest.approx(ls2) + assert ls1 == pytest.approx(_auto_ls(X2, 10.0), rel=1e-6) + + +def test_condition2_is_not_invariant_under_swapping_the_conditions(): + """Characterisation: `condition2` is as asymmetric as the default, by design. + + It is not an alternative to the default -- it is the default's mirror, and + its purpose is diagnostic. If this ever starts passing, the scheme has + stopped being a mirror and the equivalence below is no longer meaningful. + """ + fwd, _ = _fit("condition2") + rev, _ = _fit("condition2", swap=True) + assert fwd != pytest.approx(rev, rel=1e-3) + + +def test_condition2_resolves_to_what_condition1_resolves_to_on_the_swapped_pair(): + """The resolver-level statement of the mirror.""" + X1, _, X2, _ = _data() + c2, _ = _resolve_ls_scheme("condition2", X1, X2, 10.0) + # `condition1` defers to the model, so its value is `_auto_ls` of whatever + # is passed first -- which after the swap is X2. + assert c2 == _auto_ls(X2, 10.0) + + +def test_condition2_still_shares_when_a_cov_func_curry_is_supplied(): + """Regression for the `487a02e` defect, re-checked at the new value. + + The gate is `ls is None and "ls" not in function_kwargs` -- blind to the + scheme's value -- so `condition2` inherits the fix rather than needing its + own. This test is what makes that "inherits" a measurement. + """ + ls1, ls2 = _fit_with_curry("condition2") + assert ls1 == pytest.approx(ls2) + assert ls1 == pytest.approx(_fit("condition2")[0], rel=1e-6) + + +# -------------------------------------------------------------------------- +# the equivalence: de(X, Y, "condition1") == de(Y, X, "condition2") +# -------------------------------------------------------------------------- + + +def _mirror_pair(n_landmarks=0, landmarks=None, n_genes=6): + """Fit both orientations and predict both on the same points. + + Returns ``(de_forward, de_mirror, result_forward, result_mirror)``. + Forward is ``de(X, Y, "condition1")``; mirror is ``de(Y, X, "condition2")``. + """ + X1, y1, X2, y2 = _data(n_genes=n_genes) + X_new = np.vstack([X1, X2]) + + fwd = DifferentialExpression(n_landmarks=n_landmarks, random_state=0) + fwd.fit(X1, y1, X2, y2, ls_scheme="condition1", landmarks=landmarks) + + mir = DifferentialExpression(n_landmarks=n_landmarks, random_state=0) + mir.fit(X2, y2, X1, y1, ls_scheme="condition2", landmarks=landmarks) + + return ( + fwd, + mir, + fwd.predict(X_new, compute_mahalanobis=True, progress=False), + mir.predict(X_new, compute_mahalanobis=True, progress=False), + ) + + +# Measured bit-identical on the reference platform (numpy 2.4.4, jax 0.10.0, +# mellon 1.7.1, CPU); the tolerance is here so a different BLAS or an +# accelerator does not turn an equivalence into a flake. +_MIRROR_TOL = dict(rtol=1e-6, atol=1e-8) + + +def test_the_mirror_is_a_relabelling_of_the_same_fit(): + """`de(X, Y, "condition1")` and `de(Y, X, "condition2")` fit the same GPs. + + Both schemes put a single shared length scale on both models; forward takes + it from X (condition 1 there) and the mirror takes it from X as well (X is + condition 2 there). So the same GP is fitted to X and the same GP to Y, + with the model1 / model2 *roles* exchanged. + """ + fwd, mir, rf, rm = _mirror_pair() + + assert float(fwd.model1.ls) == pytest.approx(float(fwd.model2.ls)) + assert float(mir.model1.ls) == pytest.approx(float(mir.model2.ls)) + assert float(fwd.model1.ls) == pytest.approx(float(mir.model1.ls)) + + # the per-condition surfaces, with the roles exchanged + np.testing.assert_allclose( + rf["condition1_smoothed"], rm["condition2_smoothed"], **_MIRROR_TOL + ) + np.testing.assert_allclose( + rf["condition2_smoothed"], rm["condition1_smoothed"], **_MIRROR_TOL + ) + np.testing.assert_allclose( + rf["condition1_std"], rm["condition2_std"], **_MIRROR_TOL + ) + np.testing.assert_allclose( + rf["condition2_std"], rm["condition1_std"], **_MIRROR_TOL + ) + + +def test_the_mirror_flips_the_sign_of_every_signed_quantity(): + """Fold change is `condition2 - condition1`, so exchanging the labels negates it. + + This is the half of the equivalence that is *not* an equality, and stating + it as a sign flip rather than weakening the test to "both ran" is the point. + """ + _, _, rf, rm = _mirror_pair() + + # positive control: the fold change is not ~0, so the flip below is a real + # statement and not two zero arrays agreeing. + assert np.abs(np.asarray(rf["fold_change"])).max() > 1e-2 + assert not np.allclose(rf["fold_change"], rm["fold_change"]) + + for key in ("fold_change", "mean_log_fold_change", "fold_change_zscores"): + np.testing.assert_allclose( + np.asarray(rf[key]), -np.asarray(rm[key]), **_MIRROR_TOL + ) + + +def test_the_mirror_leaves_the_unsigned_statistics_unchanged(): + """Mahalanobis distance and its tail probability are blind to the exchange. + + The fold change enters the Mahalanobis form quadratically and the + denominator is `cov1 + cov2`, which is symmetric in the two models -- so + the sign flip above cancels and the gene-level statistics are identical. + This is what licenses reading a published `condition1` swap-invariance + number as a `condition2` number under label exchange. + """ + _, _, rf, rm = _mirror_pair() + + assert np.asarray(rf["mahalanobis_distances"]).max() > 0 + np.testing.assert_allclose( + rf["mahalanobis_distances"], rm["mahalanobis_distances"], **_MIRROR_TOL + ) + np.testing.assert_allclose( + rf["neg_log10_ptp"], rm["neg_log10_ptp"], **_MIRROR_TOL + ) + + +def test_the_mirror_holds_with_landmarks_when_the_landmarks_are_shared(): + """The equivalence survives the Nystrom path -- given the *same* landmarks. + + `fit()` derives automatic landmarks from `np.vstack([X_condition1, + X_condition2])`, whose ROW ORDER differs between the two orientations, so + the automatic path evaluates the two runs at different points. That is a + pre-existing, scheme-independent asymmetry; no `ls_scheme` removes it. Pass + one landmark array to both runs and the equivalence is exact again. + """ + from mellon.parameters import compute_landmarks + + X1, _, X2, _ = _data() + shared = np.asarray( + compute_landmarks( + np.vstack([X1, X2]), gp_type="fixed", n_landmarks=40, random_state=7 + ) + ) + _, _, rf, rm = _mirror_pair(n_landmarks=40, landmarks=shared) + + np.testing.assert_allclose( + rf["condition1_smoothed"], rm["condition2_smoothed"], **_MIRROR_TOL + ) + np.testing.assert_allclose( + np.asarray(rf["fold_change"]), -np.asarray(rm["fold_change"]), **_MIRROR_TOL + ) + np.testing.assert_allclose( + rf["mahalanobis_distances"], rm["mahalanobis_distances"], **_MIRROR_TOL + ) + + +# -------------------------------------------------------------------------- +# every reader of `ls_scheme`, not just the writer +# -------------------------------------------------------------------------- + + +def test_ls_scheme_is_readable_back_out_of_a_stored_params_dict(): + """`ls_scheme` is in the run-parameter match list; it must also be READABLE. + + `build_params_dict` stores it nested under `params["gp"]`, and + `_check_overwrites` compares it through `params_get`. `params_get` + resolves a nested key via `_LEGACY_MAP`; a field missing from that map + reads back as `None`, so every rerun -- including an identical one -- + compared unequal and was reported as a parameter change. + """ + from kompot.anndata.utils.params import build_params_dict, params_get + + for scheme in LS_SCHEMES: + params = build_params_dict( + {"groupby": "g", "condition1": "A", "condition2": "B"}, + gp=GPSettings(ls_scheme=scheme), + ) + assert params["gp"]["ls_scheme"] == scheme + assert params_get(params, "ls_scheme") == scheme + + +def test_the_de_cli_routes_ls_scheme_into_gp_settings(): + """A flat config key the CLI does not recognise is not rejected. + + It falls through to `**function_kwargs` and is handed to mellon, so an + omission from the routing set is silent at the config surface. `landmarks` + is deliberately absent -- it is an array, not a scalar config value. + """ + import dataclasses + + from kompot.cli.de import GP_CONFIG_KEYS + + assert "ls_scheme" in GP_CONFIG_KEYS + fields = {f.name for f in dataclasses.fields(GPSettings)} + assert GP_CONFIG_KEYS <= fields, GP_CONFIG_KEYS - fields + assert fields - GP_CONFIG_KEYS == {"landmarks"} From 082a702d13e5f07b020e642b8782de4860339924 Mon Sep 17 00:00:00 2001 From: Dominik Date: Mon, 21 Sep 2026 18:42:57 -0700 Subject: [PATCH 04/23] docs(de): the condition1/condition2 equivalence has a precondition the defaults violate Skeptic finding, and it is the right one. The equivalence was asserted UNCONDITIONALLY in settings.py, in the fit() docstring and in the CHANGELOG -- the "given identical landmarks" clause existed only in the PR body, the commit message and a TEST NAME. A boundary that lives in a test name is invisible to everyone the boundary protects. And the configuration that violates it is the LIBRARY DEFAULT. GPSettings(n_landmarks=5000, landmarks=None) computes landmarks from np.vstack([X_condition1, X_condition2]), whose ROW ORDER differs between the two orientations, so the two runs select different landmark sets. Measured at the toy scale with landmarks=None: 6 of 9 mirrored comparisons MISS assert_allclose(rtol=1e-6, atol=1e-8) -- the tolerance this package's own tests use -- and 0 of 9 miss it at n_landmarks=0. Silently. No error, no log line, nothing downstream that signals it: the numbers come back close instead of equal. "IT IS SMALL" WAS ONE CONFIGURATION GENERALISED. Landmark-fraction sweep, 800 cells, RMS-relative gap ||u-v||_2 / ||u||_2: frac n_lm smoothed mahalanobis 0.25 200 1.662e-06 1.624e-06 0.125 100 6.056e-06 4.135e-06 0.05 40 3.587e-05 1.677e-05 0.025 20 7.577e-05 4.136e-05 control, ONE shared landmark array: 0.000e+00 at both ends of the sweep So it grows ~45x (smoothed) and ~25x (Mahalanobis) as the fraction falls 10x, and the control pins the growth on landmark SELECTION rather than on the fraction itself. A default 5,000-landmark run on a large dataset sits at the LOW-fraction end, which is the bad end. An independent measurement in a separate clone agreed in direction with a somewhat larger factor. Three doc sites now carry the clause and say what happens when it does NOT hold, not only when it does: - kompot/settings.py GPSettings.ls_scheme - kompot/differential/differential_expression.py fit(), as a .. warning:: - CHANGELOG.md the [Unreleased] Added entry Plus enforcement, because docs are not a mechanism: fit() logs a warning when ls_scheme="condition2" runs with automatic landmarks. It is scoped to that scheme deliberately -- condition2's ONLY purpose is the cross-orientation comparison, so anyone who reaches that branch is relying on the equivalence. Four tests, and the positive one earned its place immediately: it FAILED on the first run, because kompot/__init__.py configures the logger with propagate=False, so pytest's caplog sees an EMPTY record list. That does not fail loudly -- it makes every "did not warn" assertion pass for a reason unrelated to the code. The three negative tests were vacuous until the positive one flushed it out. caplog is now replaced by a handler attached directly to the kompot logger, with a positive control for the helper itself. Mutation control: forcing the warning to fire unconditionally flips exactly the three negative tests and leaves the two positives green. tests/test_ls_scheme.py: 33 -> 38. No behaviour change beyond the added log line. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 22 +++- .../differential/differential_expression.py | 54 ++++++++++ kompot/settings.py | 6 +- tests/test_ls_scheme.py | 102 ++++++++++++++++++ 4 files changed, 178 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 124e7b5..aee84e9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,11 +23,23 @@ All notable changes to this project will be documented in this file. `"condition2"` is the default's mirror — it estimates from condition 2's cells and hands the value to condition 1. It is not an alternative to the default but a **diagnostic**: `de(X, Y, ls_scheme="condition1")` and `de(Y, X, ls_scheme="condition2")` are the same - computation with the labels exchanged (measured bit-identical for the smoothed surfaces, the - posterior standard deviations, the Mahalanobis distances and `neg_log10_ptp`; the fold change - and its z-scores flip sign, as the definition `condition2 - condition1` requires). So the - default's swap-dependence can be exhibited from a single call site, without swapping - `condition1` and `condition2` and re-deriving which sign is which. + computation with the labels exchanged, so the default's swap-dependence can be exhibited from + a single call site without swapping `condition1` and `condition2` and re-deriving which sign + is which. **Given the same landmarks in both runs** the equivalence is exact — measured + bit-identical for the smoothed surfaces, the posterior standard deviations, the Mahalanobis + distances and `neg_log10_ptp`, with the fold change and its z-scores flipping sign as the + definition `condition2 - condition1` requires. + + **That precondition is not satisfied by the defaults.** Automatic landmarks are computed from + the two conditions' cells stacked in the order given, so the two orientations select + *different* landmark sets; with `n_landmarks=5000, landmarks=None` the equivalence is only + approximate and nothing downstream signals it. Measured at a small scale, 6 of 9 mirrored + comparisons then miss `rtol=1e-6, atol=1e-8`, and the gap **grows as the landmark fraction + falls** — roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as the + fraction goes from 0.25 to 0.025 — so a default run on a large dataset sits at the worse end. + Pass one `landmarks` array to both runs, or use `n_landmarks=0`; both are exact. `fit()` now + logs a warning when `"condition2"` is used with automatic landmarks. The order dependence is + pre-existing and independent of `ls_scheme`: no scheme removes it. `"symmetric"` shares the size-weighted geometric mean of the two conditions' own estimates — the same estimator, with nearest neighbours looked up *within* each condition — and is diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 76b5bd8..87fa435 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -428,6 +428,33 @@ def fit( under the two schemes exposes the default's swap-dependence from a single call site, without having to swap ``condition1`` and ``condition2`` at the call site and re-derive which sign is which. + + .. warning:: + + **That equivalence holds only if both runs use the same + landmarks, and the library defaults do not.** Automatic + landmarks are computed from + ``np.vstack([X_condition1, X_condition2])``, whose row order + differs between the two orientations, so the two runs select + *different* landmark sets and evaluate at different points. + With ``n_landmarks=5000`` and ``landmarks=None`` — the + :class:`~kompot.settings.GPSettings` defaults — the mirror is + therefore only approximate, and nothing downstream signals it: + the numbers are silently close rather than equal. Measured at + a small scale, **6 of 9 mirrored comparisons miss** + ``assert_allclose(rtol=1e-6, atol=1e-8)`` — the tolerance this + package's own tests use — and the gap **grows as the landmark + fraction falls**, so a default 5,000-landmark run on a large + dataset sits at the worse end, not the better one. + + To rely on the equivalence, make the landmarks common to both + runs: pass one ``landmarks`` array to both, or set + ``n_landmarks=0``. Both are exact. ``fit()`` logs a warning + when ``"condition2"`` is used with automatic landmarks. + + The order dependence is **pre-existing and independent of** + ``ls_scheme`` — no scheme removes it, and it affects any + comparison of two orientations, not only this one. * ``"symmetric"`` — estimate a length scale from each condition separately and share their size-weighted geometric mean. Invariant under swapping the two conditions, and it does not inherit the @@ -500,6 +527,7 @@ def fit( ) # Compute shared landmarks (needs both conditions) + landmarks_were_provided = landmarks is not None if self.function_predictor1 is None or self.function_predictor2 is None: if landmarks is not None: logger.info(f"Using provided landmarks with shape {landmarks.shape}") @@ -514,6 +542,32 @@ def fit( ) self.computed_landmarks = landmarks + # `"condition2"` exists ONLY to be compared against `"condition1"` on the + # swapped orientation, and that comparison is exact only when both runs + # evaluate at the same landmarks. Automatic landmarks are derived from + # `np.vstack([X_condition1, X_condition2])`, whose ROW ORDER differs + # between the two orientations, so the two runs pick different sets and + # the equivalence silently degrades to approximate -- no error, no + # downstream signal, just numbers that are close instead of equal. The + # scheme's whole purpose is the comparison, so anyone who reaches here is + # relying on it; say so rather than leaving the boundary in the docs. + if ( + ls_scheme == "condition2" + and not landmarks_were_provided + and self.n_landmarks is not None + and self.n_landmarks > 0 + ): + logger.warning( + "ls_scheme='condition2' with automatic landmarks: the " + "condition1/condition2 equivalence is EXACT only when both " + "orientations use the same landmarks, and automatic landmarks " + "are order-dependent (they are computed from the two " + "conditions' cells stacked in the order given). Pass one " + "`landmarks` array to both runs, or use n_landmarks=0, if you " + "intend to compare the two orientations. The gap grows as the " + "landmark fraction falls." + ) + # -- Fit model1 -- if self.model1 is None: self.model1 = ExpressionModel( diff --git a/kompot/settings.py b/kompot/settings.py index a0e2236..c6f84e0 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -49,7 +49,11 @@ class GPSettings: is its mirror and exists to make that dependence visible in a single run, since ``de(X, Y, ls_scheme="condition1")`` and ``de(Y, X, ls_scheme="condition2")`` are the same computation with the - labels exchanged; ``"symmetric"`` shares the size-weighted geometric + labels exchanged **provided both runs use the same landmarks** — see + the warning under :meth:`kompot.differential.DifferentialExpression.fit`, + because the defaults here (``n_landmarks=5000``, ``landmarks=None``) do + *not* satisfy that and the equivalence then holds only approximately; + ``"symmetric"`` shares the size-weighted geometric mean of the two per-condition estimates; ``"pooled"`` estimates from both conditions' cells taken together; ``"separate"`` gives each condition its own. See diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py index 8aa1359..9cf42bc 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_ls_scheme.py @@ -9,6 +9,9 @@ are the same computation with the labels exchanged. """ +import contextlib +import logging + import numpy as np import pytest @@ -433,3 +436,102 @@ def test_the_de_cli_routes_ls_scheme_into_gp_settings(): fields = {f.name for f in dataclasses.fields(GPSettings)} assert GP_CONFIG_KEYS <= fields, GP_CONFIG_KEYS - fields assert fields - GP_CONFIG_KEYS == {"landmarks"} + + +class _CaptureWarnings(logging.Handler): + """Collect WARNING records straight off the `kompot` logger. + + `caplog` cannot see them: `kompot/__init__.py` configures the logger with + `"propagate": False`, so records never reach the root handler pytest + installs. Using `caplog` here does not fail loudly -- it yields an EMPTY + record list, which makes a "did not warn" assertion pass for a reason that + has nothing to do with the code under test. The three negative tests below + were vacuous until the positive one flushed this out. + """ + + def __init__(self): + super().__init__(level=logging.WARNING) + self.messages = [] + + def emit(self, record): + self.messages.append(record.getMessage()) + + +@contextlib.contextmanager +def _captured_kompot_warnings(): + logger = logging.getLogger("kompot") + handler = _CaptureWarnings() + previous = logger.level + logger.addHandler(handler) + logger.setLevel(min(previous or logging.WARNING, logging.WARNING)) + try: + yield handler.messages + finally: + logger.removeHandler(handler) + logger.setLevel(previous) + + +def test_the_warning_capture_helper_actually_captures(): + """Positive control for the helper itself, so the negatives below mean something.""" + with _captured_kompot_warnings() as messages: + logging.getLogger("kompot").warning("canary: same landmarks") + assert any("same landmarks" in m for m in messages) + + +def test_condition2_warns_when_the_landmarks_are_not_shared(): + """The equivalence's precondition must be audible at runtime, not just in docs. + + `"condition2"` exists only to be compared against `"condition1"` on the + swapped orientation, and that comparison is exact only when both runs + evaluate at the same landmarks. Automatic landmarks are order-dependent, so + the default configuration (`n_landmarks=5000`, `landmarks=None`) degrades the + equivalence to approximate with nothing downstream to signal it -- measured, + 6 of 9 mirrored comparisons then miss the tolerance asserted above. A + boundary that lives only in prose is invisible to the people it protects. + """ + X1, y1, X2, y2 = _data() + + with _captured_kompot_warnings() as messages: + DifferentialExpression(n_landmarks=20, random_state=0).fit( + X1, y1, X2, y2, ls_scheme="condition2" + ) + assert any("same landmarks" in m for m in messages), messages + + +@pytest.mark.parametrize( + "kwargs", + [ + pytest.param(dict(n_landmarks=0), id="no-landmarks"), + pytest.param(dict(n_landmarks=20, shared=True), id="landmarks-shared"), + ], +) +def test_condition2_does_not_warn_when_the_equivalence_is_exact(kwargs): + """The other half: a warning that always fires teaches nothing.""" + from mellon.parameters import compute_landmarks + + X1, y1, X2, y2 = _data() + kwargs = dict(kwargs) + landmarks = None + if kwargs.pop("shared", False): + landmarks = np.asarray( + compute_landmarks( + np.vstack([X1, X2]), gp_type="fixed", n_landmarks=20, random_state=7 + ) + ) + + with _captured_kompot_warnings() as messages: + DifferentialExpression(random_state=0, **kwargs).fit( + X1, y1, X2, y2, ls_scheme="condition2", landmarks=landmarks + ) + assert not [m for m in messages if "same landmarks" in m] + + +def test_the_other_schemes_do_not_warn_about_landmarks(): + """Only `condition2` makes the cross-orientation comparison its purpose.""" + X1, y1, X2, y2 = _data() + for scheme in ("condition1", "symmetric", "pooled", "separate"): + with _captured_kompot_warnings() as messages: + DifferentialExpression(n_landmarks=20, random_state=0).fit( + X1, y1, X2, y2, ls_scheme=scheme + ) + assert not [m for m in messages if "same landmarks" in m], scheme From 55d4ef279ad642a32541f4a81d83da141bb98b6d Mon Sep 17 00:00:00 2001 From: Dominik Date: Mon, 21 Sep 2026 19:07:16 -0700 Subject: [PATCH 05/23] docs(de): name the seed beside the 6-of-9 count, and say the default is worse Skeptic follow-up, and it is right that a bare count invites a reader to reproduce it and conclude one of us is wrong. `6 of 9` is a `random_state=0` measurement. `GPSettings.random_state` defaults to None, and the landmark draw is then not reproducible from one run to the next -- so the count at the documented default is not a fixed number at all. Measured, three consecutive runs at `random_state=None`, same toy data, same `assert_allclose(rtol=1e-6, atol=1e-8)`: run 1 7 of 9 MISS run 2 8 of 9 MISS run 3 7 of 9 MISS The skeptic reported a stable 7 across its three runs; mine were 7, 8, 7. Both observations are of the same thing -- at `random_state=None` the count is a draw, not a constant -- so the disclosure now states the seed beside the count AND that the default is worse and unstable, rather than quoting a second bare number that a reader could also fail to reproduce. `mean_log_fold_change` is the only comparison that never missed in any run. Three sites updated in step: the fit() warning block, the CHANGELOG entry, and the docstring of test_condition2_warns_when_the_landmarks_are_not_shared. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 9 ++++++--- kompot/differential/differential_expression.py | 11 ++++++++--- tests/test_ls_scheme.py | 6 ++++-- 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aee84e9..fbe3759 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,9 +34,12 @@ All notable changes to this project will be documented in this file. the two conditions' cells stacked in the order given, so the two orientations select *different* landmark sets; with `n_landmarks=5000, landmarks=None` the equivalence is only approximate and nothing downstream signals it. Measured at a small scale, 6 of 9 mirrored - comparisons then miss `rtol=1e-6, atol=1e-8`, and the gap **grows as the landmark fraction - falls** — roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as the - fraction goes from 0.25 to 0.025 — so a default run on a large dataset sits at the worse end. + comparisons then miss `rtol=1e-6, atol=1e-8` **at `random_state=0`** — and at the default + `random_state=None` the landmark draw is not reproducible between runs, so the count is worse + and unstable (7, 8, 7 over three consecutive runs). The gap also **grows as the landmark + fraction falls** — roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as + the fraction goes from 0.25 to 0.025 — so a default run on a large dataset sits at the worse + end. Pass one `landmarks` array to both runs, or use `n_landmarks=0`; both are exact. `fit()` now logs a warning when `"condition2"` is used with automatic landmarks. The order dependence is pre-existing and independent of `ls_scheme`: no scheme removes it. diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 87fa435..0ad6abb 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -443,9 +443,14 @@ def fit( the numbers are silently close rather than equal. Measured at a small scale, **6 of 9 mirrored comparisons miss** ``assert_allclose(rtol=1e-6, atol=1e-8)`` — the tolerance this - package's own tests use — and the gap **grows as the landmark - fraction falls**, so a default 5,000-landmark run on a large - dataset sits at the worse end, not the better one. + package's own tests use. That count is taken at + ``random_state=0``; at the :class:`~kompot.settings.GPSettings` + default ``random_state=None`` the landmark draw is not + reproducible from one run to the next, so the count is **worse + and unstable** — 7, 8 and 7 over three consecutive runs. The + gap also **grows as the landmark fraction falls**, so a default + 5,000-landmark run on a large dataset sits at the worse end, + not the better one. To rely on the equivalence, make the landmarks common to both runs: pass one ``landmarks`` array to both, or set diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py index 9cf42bc..02abfa1 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_ls_scheme.py @@ -486,8 +486,10 @@ def test_condition2_warns_when_the_landmarks_are_not_shared(): evaluate at the same landmarks. Automatic landmarks are order-dependent, so the default configuration (`n_landmarks=5000`, `landmarks=None`) degrades the equivalence to approximate with nothing downstream to signal it -- measured, - 6 of 9 mirrored comparisons then miss the tolerance asserted above. A - boundary that lives only in prose is invisible to the people it protects. + 6 of 9 mirrored comparisons then miss the tolerance asserted above at + `random_state=0`, and 7-8 of 9 at the `random_state=None` default, where the + landmark draw is not even reproducible between runs. A boundary that lives + only in prose is invisible to the people it protects. """ X1, y1, X2, y2 = _data() From 2bff8e34cf49e6307cff75b8d6e6dce62a4c398a Mon Sep 17 00:00:00 2001 From: Dominik Date: Mon, 21 Sep 2026 19:20:09 -0700 Subject: [PATCH 06/23] docs(de): replace the miss-count range with a per-comparison localisation My previous commit published "7, 8 and 7 over three consecutive runs" at random_state=None. That is refuted, and the lower bound was wrong in the direction that undermines the disclosure: a 6 occurs at the default too. Twelve independent draws, random_state=None, n_landmarks=40, toy scale, against the suite's own assert_allclose(rtol=1e-6, atol=1e-8): 12/12 condition1_smoothed <-> condition2_smoothed 12/12 condition2_smoothed <-> condition1_smoothed 8/12 condition1_std <-> condition2_std 2/12 condition2_std <-> condition1_std 12/12 fold_change vs NEGATED 1/12 mean_log_fold_change vs NEGATED 12/12 fold_change_zscores vs NEGATED 12/12 mahalanobis_distances 12/12 neg_log10_ptp per-draw totals: 6 x4, 7 x6, 8 x1, 9 x1 (observed range 6-9, not 7-8) random_state=0 reference: total 6, and THE SAME SIX comparisons So a range quoted from three samples was three samples dressed as a range. What is actually stable is WHICH comparisons miss, and that is what is now published: six always-miss, the two posterior standard-deviation comparisons carrying nearly all the variance, and mean_log_fold_change nearly exempt because averaging over cells cancels most of the landmark noise -- nearly, not immune: it missed once in twelve, so "never" would have been the same mistake one level down. Six is therefore a FLOOR, not a typical value, and it is also the seeded figure with the identical comparison set -- so a reader who runs once at the default and sees 6 learns nothing about the seed either way. Stating the localisation rather than a count makes the seed parenthetical nearly unnecessary, and it is reproducible where a range is not. An independent review reported a stable 7 across three draws and 0/10 for mean_log_fold_change; both are consistent with a draw whose modal value is 7 and whose tails are thin, and neither supports a range. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 17 ++++++----- .../differential/differential_expression.py | 30 ++++++++++++------- tests/test_ls_scheme.py | 11 +++---- 3 files changed, 36 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fbe3759..1bda902 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,13 +33,16 @@ All notable changes to this project will be documented in this file. **That precondition is not satisfied by the defaults.** Automatic landmarks are computed from the two conditions' cells stacked in the order given, so the two orientations select *different* landmark sets; with `n_landmarks=5000, landmarks=None` the equivalence is only - approximate and nothing downstream signals it. Measured at a small scale, 6 of 9 mirrored - comparisons then miss `rtol=1e-6, atol=1e-8` **at `random_state=0`** — and at the default - `random_state=None` the landmark draw is not reproducible between runs, so the count is worse - and unstable (7, 8, 7 over three consecutive runs). The gap also **grows as the landmark - fraction falls** — roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as - the fraction goes from 0.25 to 0.025 — so a default run on a large dataset sits at the worse - end. + approximate and nothing downstream signals it. Measured at a small scale against + `rtol=1e-6, atol=1e-8`, **six of the nine mirrored comparisons miss in every draw** — both + smoothed surfaces, `fold_change`, `fold_change_zscores`, `mahalanobis_distances` and + `neg_log10_ptp`. The rest are intermittent: the two posterior standard-deviation comparisons + missed 8/12 and 2/12 draws at the default `random_state=None`, and `mean_log_fold_change` 1/12 + (averaging over cells cancels most of the landmark noise). Totals ran 6 to 9 of 9; 0 of 9 miss + at `n_landmarks=0`. Six is a floor rather than a typical value, and it is also the + `random_state=0` figure. The gap additionally **grows as the landmark fraction falls** — + roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as the fraction goes + from 0.25 to 0.025 — so a default run on a large dataset sits at the worse end. Pass one `landmarks` array to both runs, or use `n_landmarks=0`; both are exact. `fit()` now logs a warning when `"condition2"` is used with automatic landmarks. The order dependence is pre-existing and independent of `ls_scheme`: no scheme removes it. diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 0ad6abb..5e11b4a 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -440,17 +440,27 @@ def fit( With ``n_landmarks=5000`` and ``landmarks=None`` — the :class:`~kompot.settings.GPSettings` defaults — the mirror is therefore only approximate, and nothing downstream signals it: - the numbers are silently close rather than equal. Measured at - a small scale, **6 of 9 mirrored comparisons miss** + the numbers are silently close rather than equal. + + Measured at a small scale against ``assert_allclose(rtol=1e-6, atol=1e-8)`` — the tolerance this - package's own tests use. That count is taken at - ``random_state=0``; at the :class:`~kompot.settings.GPSettings` - default ``random_state=None`` the landmark draw is not - reproducible from one run to the next, so the count is **worse - and unstable** — 7, 8 and 7 over three consecutive runs. The - gap also **grows as the landmark fraction falls**, so a default - 5,000-landmark run on a large dataset sits at the worse end, - not the better one. + package's own tests use — **six of the nine mirrored + comparisons miss in every draw**: both smoothed surfaces, + ``fold_change``, ``fold_change_zscores``, + ``mahalanobis_distances`` and ``neg_log10_ptp`` (12/12 draws + each). The other three are intermittent — the two posterior + standard-deviation comparisons missed 8/12 and 2/12 draws at + the default ``random_state=None``, and ``mean_log_fold_change`` + 1/12, since averaging over cells cancels most of the landmark + noise. Totals ran 6 to 9 of 9. **Six is a floor, not a typical + value**, and it is also what ``random_state=0`` gives, so a + single run at the default can coincide with the seeded number + while nothing about the seed is thereby established. At + ``n_landmarks=0``, 0 of 9 miss. + + The gap also **grows as the landmark fraction falls**, so a + default 5,000-landmark run on a large dataset sits at the worse + end, not the better one. To rely on the equivalence, make the landmarks common to both runs: pass one ``landmarks`` array to both, or set diff --git a/tests/test_ls_scheme.py b/tests/test_ls_scheme.py index 02abfa1..b0e5745 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_ls_scheme.py @@ -485,11 +485,12 @@ def test_condition2_warns_when_the_landmarks_are_not_shared(): swapped orientation, and that comparison is exact only when both runs evaluate at the same landmarks. Automatic landmarks are order-dependent, so the default configuration (`n_landmarks=5000`, `landmarks=None`) degrades the - equivalence to approximate with nothing downstream to signal it -- measured, - 6 of 9 mirrored comparisons then miss the tolerance asserted above at - `random_state=0`, and 7-8 of 9 at the `random_state=None` default, where the - landmark draw is not even reproducible between runs. A boundary that lives - only in prose is invisible to the people it protects. + equivalence to approximate with nothing downstream to signal it. Measured + against the tolerance asserted above, six of the nine mirrored comparisons + miss in EVERY draw and the totals ran 6 to 9 of 9 at the + `random_state=None` default, where the landmark draw is not even + reproducible between runs. A boundary that lives only in prose is + invisible to the people it protects. """ X1, y1, X2, y2 = _data() From c907da9d431b515bdcbb816c8c3d0578d35d704a Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 09:40:14 -0700 Subject: [PATCH 07/23] test(da): differential abundance is swap-equivalent at its defaults The operator asked to make sure differential abundance works equivalently. It already does, structurally, and no ls_scheme analogue is warranted -- establishing that, rather than adding a knob, is the answer. WHY IT HOLDS. DifferentialAbundance.fit builds BOTH density estimators from the same estimator_defaults dict, with no `ls_for_model2 = self.model1.ls` analogue anywhere in the module. There is no condition-1 inheritance to be asymmetric about, so the defect ls_scheme exists to expose in differential EXPRESSION has no counterpart in abundance. Read at 082a702. The argument rests on one line, so it was measured rather than asserted: `**d` unpacks into a fresh dict, but the VALUES are shared references, and an in-place mutation of the landmarks array by estimator 1's fit would make the result luck rather than structure. Measured -- estimator_defaults is byte-identical after CONSTRUCTING estimator 1 and again after FITTING it, landmark object identity stable, and both estimators' landmarks equal the dict's. MEASURED at the DA defaults (n_landmarks=None, sync_parameters=False; 400+400 cells, 2-D, seed 0). The configuration is named because it is the REASON: with no landmarks there is no np.vstack row-order dependence to be sensitive to. log_density_c1(A) vs log_density_c2(B) 0.000e+00 log_density_c2(A) vs log_density_c1(B) 0.000e+00 log_fold_change_uncertainty 0.000e+00 neg_log10_fold_change_ptp 0.000e+00 log_fold_change vs NEGATED 0.000e+00 log_fold_change_zscore vs NEGATED 0.000e+00 direction: A vs SWAPPED(B) 800/800 agree direction: A vs UNSWAPPED(B) 0/800 agree That is one dataset at one configuration. It generalises through the mechanism above, not through the zero. TWO CONFIGURATIONS ARE ONLY APPROXIMATE, and they are findings, not footnotes -- both relatively WORSE than the 1.5e-05 landmark gap measured for DE: sync_parameters=True log_fold_change 1.220e-02, zscore 1.223e-02 n_landmarks=40 (automatic) uncertainty 5.634e-02, zscore 2.858e-02 Same root cause, wider blast radius: np.vstack([X_condition1, X_condition2]) is row-order dependent, and under sync_parameters it feeds compute_d_factal / compute_mu / compute_ls as well as compute_landmarks. So the landmark precondition documented for DE is in scope for DA and has a SECOND instance there. Sharing one landmark array restores exactness; that is tested. THE POSITIVE CONTROL IS WHY THESE TESTS ARE WORTH ANYTHING. My first two direction checks read 160/160 and 600/600 agree and were VACUOUS: every cell was 'neutral' (max|lfc| 0.07 then 0.42 against the package's own default threshold of 1.0), and two all-neutral arrays agree perfectly whether or not the labels are swapped. A single arm reading 600/600 looks like a triumph; only running the UNSWAPPED arm exposed it. The generalisable part is why a ported test could not have caught it: differential expression has no categorical output, so a control built for the DE family structurally under-covers the member being added. A ported shape inherits the green without earning it. The direction test therefore asserts both arms and requires 'up' and 'down' to actually occur; separation was increased until |log_fold_change| clears the DEFAULT threshold, and no threshold was tuned. tests/test_da_swap_equivalence.py: 8 tests, all passing. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 88 +++++------- tests/test_da_swap_equivalence.py | 222 ++++++++++++++++++++++++++++++ 2 files changed, 255 insertions(+), 55 deletions(-) create mode 100644 tests/test_da_swap_equivalence.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 1bda902..0bee9a6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,72 +6,50 @@ All notable changes to this project will be documented in this file. ### Added - - **`GPSettings.ls_scheme`** — chooses which cells the automatic shared length scale is - estimated from in differential expression. The default, `"condition1"`, is exactly the - existing behaviour and nothing changes unless you ask for something else. - - The two conditions have always been smoothed at a **shared** length scale so that their - fitted surfaces are comparable, and that part is right: giving each condition its own - (`ls_scheme="separate"`) makes the two surfaces differently smooth, and the mismatch alone - manufactures apparent fold changes. What was undocumented is that the shared value is - estimated from **condition 1's cells only**, so it is a function of `n_condition1` and - `de(adata, condition1=X, condition2=Y)` is *not* equivalent to - `de(adata, condition1=Y, condition2=X)`. On an exchangeable null with nothing to find, the - two orientations of the same partition differ in false-positive rate by 0.06–0.08 at - 2:1 and beyond. - - `"condition2"` is the default's mirror — it estimates from condition 2's cells and hands the - value to condition 1. It is not an alternative to the default but a **diagnostic**: - `de(X, Y, ls_scheme="condition1")` and `de(Y, X, ls_scheme="condition2")` are the same - computation with the labels exchanged, so the default's swap-dependence can be exhibited from - a single call site without swapping `condition1` and `condition2` and re-deriving which sign - is which. **Given the same landmarks in both runs** the equivalence is exact — measured - bit-identical for the smoothed surfaces, the posterior standard deviations, the Mahalanobis - distances and `neg_log10_ptp`, with the fold change and its z-scores flipping sign as the - definition `condition2 - condition1` requires. - - **That precondition is not satisfied by the defaults.** Automatic landmarks are computed from - the two conditions' cells stacked in the order given, so the two orientations select - *different* landmark sets; with `n_landmarks=5000, landmarks=None` the equivalence is only - approximate and nothing downstream signals it. Measured at a small scale against - `rtol=1e-6, atol=1e-8`, **six of the nine mirrored comparisons miss in every draw** — both - smoothed surfaces, `fold_change`, `fold_change_zscores`, `mahalanobis_distances` and - `neg_log10_ptp`. The rest are intermittent: the two posterior standard-deviation comparisons - missed 8/12 and 2/12 draws at the default `random_state=None`, and `mean_log_fold_change` 1/12 - (averaging over cells cancels most of the landmark noise). Totals ran 6 to 9 of 9; 0 of 9 miss - at `n_landmarks=0`. Six is a floor rather than a typical value, and it is also the - `random_state=0` figure. The gap additionally **grows as the landmark fraction falls** — - roughly 45x on the smoothed surfaces and 25x on the Mahalanobis distances as the fraction goes - from 0.25 to 0.025 — so a default run on a large dataset sits at the worse end. - Pass one `landmarks` array to both runs, or use `n_landmarks=0`; both are exact. `fit()` now - logs a warning when `"condition2"` is used with automatic landmarks. The order dependence is - pre-existing and independent of `ls_scheme`: no scheme removes it. - - `"symmetric"` shares the size-weighted geometric mean of the two conditions' own estimates - — the same estimator, with nearest neighbours looked up *within* each condition — and is - invariant under swapping the conditions. `"pooled"` estimates from the union, which is also - swap-invariant but yields a systematically smaller value, because the union is denser than - either condition and nearest-neighbour distances shrink with cell count alone. + - **`GPSettings.ls_scheme`**: chooses which cells the automatic shared length scale is estimated + from in differential expression — `"condition1"` (the default, exactly the existing behaviour), + `"condition2"`, `"symmetric"`, `"pooled"`, `"separate"`. The shared value has always been + estimated from **condition 1's cells only**, making it a function of `n_condition1`, so + `de(condition1=X, condition2=Y)` is *not* equivalent to `de(condition1=Y, condition2=X)`; on an + exchangeable null the two orientations differ in false-positive rate by 0.06–0.08 beyond 2:1. + `"symmetric"` and `"pooled"` are swap-invariant, `"separate"` is diagnostics-only. See + `DifferentialExpression.fit` for which to pick and why. + - **`ls_scheme="condition2"`**: the default's mirror — estimated from condition 2's cells — for + diagnosis rather than analysis. `de(X, Y, ls_scheme="condition1")` and + `de(Y, X, ls_scheme="condition2")` are the same computation with the labels exchanged, so the + default's swap-dependence is visible from a single call site. The equivalence is exact only + when both runs use the same landmarks; automatic landmarks are order-dependent, and `fit()` + warns when the scheme is used with them. Measurements and the precondition: + `DifferentialExpression.fit`. ### Changed - - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme - is no longer treated as a matching rerun. + - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme is + no longer treated as a matching rerun. ### Fixed - - `ls_scheme` is now readable back out of a stored `run_info` params dict. It is written nested - under `params["gp"]`, but `params_get` resolves nested keys through `_LEGACY_MAP` and the - field was missing from it — so it read back as `None` and *every* rerun, including an - identical one, compared unequal and was reported as a parameter change. - - The DE CLI now routes a flat `ls_scheme` config key into `GPSettings`. An unrecognised key + - **`ls_scheme` is now readable back out of a stored `run_info`.** It is written nested under + `params["gp"]` but was missing from `params_get`'s key map, so it read back as `None` and + *every* rerun — including a byte-identical one — compared unequal and was reported as a + parameter change. + - **The DE CLI now routes a flat `ls_scheme` config key into `GPSettings`.** An unrecognised key is not rejected: it falls through to `**function_kwargs` and is handed to mellon, so the documented knob was unreachable from a config file and the omission was silent. ### Documentation - - `GPSettings.ls` and `DifferentialExpression.fit` now state where the shared length scale - comes from and that the default makes the contrast depend on argument order. + - `GPSettings.ls_scheme` and `DifferentialExpression.fit` now state where the shared length scale + comes from, that the default makes the contrast depend on argument order, and what the + `"condition2"` equivalence requires of landmarks. + - **Differential abundance is swap-equivalent at its defaults**, and now has tests saying so: + `da(X, Y)` and `da(Y, X)` give the same densities with the roles exchanged, negated log fold + changes and z-scores, identical uncertainty, PTP and significant-cell set, and exchanged + `'up'`/`'down'` labels. It needs no `ls_scheme` analogue — abundance builds both density + estimators from one shared configuration and never had the condition-1 inheritance + `ls_scheme` exists to expose. The equivalence is approximate rather than exact under + `sync_parameters=True` or automatic landmarks, both because the stacked union is row-order + dependent; pass one `landmarks` array to both runs to restore it. ## [0.8.0] - 2026-07-28 diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py new file mode 100644 index 0000000..9239c7a --- /dev/null +++ b/tests/test_da_swap_equivalence.py @@ -0,0 +1,222 @@ +"""Swap equivalence for differential ABUNDANCE. + +`ls_scheme` and `tests/test_ls_scheme.py` are about differential EXPRESSION, where +the shared length scale was estimated from condition 1's cells alone and +`"condition2"` exists to expose that. Abundance is a separate code path and +inherits none of it, so the property is established here rather than assumed. + +It holds for a structural reason worth stating: `DifferentialAbundance.fit` builds +BOTH density estimators from the same ``estimator_defaults`` dict -- + + density_estimator_condition1 = mellon.DensityEstimator(**estimator_defaults) + density_estimator_condition2 = mellon.DensityEstimator(**estimator_defaults) + +-- with no ``ls_for_model2 = self.model1.ls`` analogue anywhere. There is no +condition-1 inheritance to be asymmetric about, so DA never had the defect +`ls_scheme` was added for, and no DA counterpart of `"condition2"` is warranted. + +Two configurations DO break the exactness, and both are covered below: automatic +landmarks, and ``sync_parameters=True``. Both trace to the same construction as +the DE case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent +-- but under ``sync_parameters`` it feeds ``d``, ``mu`` and ``ls`` as well as the +landmarks, so DA has that hazard in two places rather than one. +""" + +import numpy as np +import pytest + +from kompot.differential.differential_abundance import DifferentialAbundance + +# Well separated so that |log_fold_change| clears the package's OWN default +# `log_fold_change_threshold = 1.0`. No threshold is tuned anywhere in this file: +# an earlier version of these tests used weaker separation, every cell came back +# 'neutral', and the direction check below passed while being vacuous. +_SEPARATION = 2.5 +_SCALE = 0.45 + + +def _data(n=400, seed=0): + rng = np.random.default_rng(seed) + X1 = rng.normal(loc=[-_SEPARATION, 0.0], scale=_SCALE, size=(n, 2)) + X2 = rng.normal(loc=[+_SEPARATION, 0.0], scale=_SCALE, size=(n, 2)) + return X1, X2 + + +def _swap_pair(n_landmarks=None, sync_parameters=False, landmarks=None): + """Fit both orientations and predict both on the same points.""" + X1, X2 = _data() + X_new = np.vstack([X1, X2]) + + fwd = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) + fwd.fit(X1, X2, sync_parameters=sync_parameters, landmarks=landmarks) + + rev = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) + rev.fit(X2, X1, sync_parameters=sync_parameters, landmarks=landmarks) + + return ( + fwd.predict(X_new, progress=False), + rev.predict(X_new, progress=False), + ) + + +# Measured bit-identical on the reference platform at the DA defaults; the +# tolerance is here so a different BLAS or an accelerator does not turn an +# equivalence into a flake. +_TOL = dict(rtol=1e-6, atol=1e-8) + + +@pytest.fixture(scope="module") +def default_pair(): + return _swap_pair() + + +def test_the_swap_is_a_relabelling_of_the_same_densities(default_pair): + """`da(X, Y)` and `da(Y, X)` fit the same density to each condition. + + Both estimators are constructed from one shared ``estimator_defaults``, and at + the defaults nothing in that dict depends on the ORDER the conditions were + passed in -- so each condition's density follows the condition, not the slot. + """ + fwd, rev = default_pair + np.testing.assert_allclose( + fwd["log_density_condition1"], rev["log_density_condition2"], **_TOL + ) + np.testing.assert_allclose( + fwd["log_density_condition2"], rev["log_density_condition1"], **_TOL + ) + + +def test_the_swap_flips_the_sign_of_every_signed_quantity(default_pair): + """`log_fold_change` is `condition2 - condition1`, so exchanging labels negates it.""" + fwd, rev = default_pair + + # positive control: the fold changes are not ~0, so the flip is a real claim + assert np.abs(np.asarray(fwd["log_fold_change"])).max() > 1.0 + assert not np.allclose(fwd["log_fold_change"], rev["log_fold_change"]) + + for key in ("log_fold_change", "log_fold_change_zscore"): + np.testing.assert_allclose( + np.asarray(fwd[key]), -np.asarray(rev[key]), **_TOL + ) + + +def test_the_swap_leaves_the_unsigned_statistics_unchanged(default_pair): + """Uncertainty is `u1 + u2` and the PTP is `min(logcdf(z), logcdf(-z))`. + + Both are symmetric in the two conditions -- the first because addition is, the + second because it is an even function of the z-score -- so the sign flip above + cancels and these come back identical. + """ + fwd, rev = default_pair + np.testing.assert_allclose( + fwd["log_fold_change_uncertainty"], rev["log_fold_change_uncertainty"], **_TOL + ) + np.testing.assert_allclose( + fwd["neg_log10_fold_change_ptp"], rev["neg_log10_fold_change_ptp"], **_TOL + ) + + +def test_the_swap_exchanges_the_direction_labels(default_pair): + """`log_fold_change_direction` is categorical: 'up' and 'down' must EXCHANGE. + + This is the one abundance output with no differential-expression counterpart, + so porting the DE test shape does not cover it. The control is load-bearing + rather than decorative: with a weaker density difference every cell comes back + 'neutral', and two all-'neutral' arrays agree perfectly *whether or not* the + labels are swapped -- so the swap assertion passes while establishing nothing. + Both halves are asserted here, and the unswapped comparison must DISAGREE. + """ + fwd, rev = default_pair + forward = np.asarray(fwd["log_fold_change_direction"]) + reverse = np.asarray(rev["log_fold_change_direction"]) + + # control 1: both labels actually occur, so the swap has something to exchange + assert (forward == "up").sum() > 0 + assert (forward == "down").sum() > 0 + + swap = {"up": "down", "down": "up", "neutral": "neutral"} + expected = np.array([swap[d] for d in reverse], dtype=object) + assert np.array_equal(forward, expected) + + # control 2: the swap is not a no-op -- without it the two must disagree + assert not np.array_equal(forward, reverse) + + +def test_the_significant_set_is_the_same_cells_either_way(default_pair): + """Significance is `|log_fold_change| > threshold` and a PTP, both sign-blind. + + So the swap may relabel a cell's direction but must never change WHICH cells + are called. + """ + fwd, rev = default_pair + fwd_sig = np.asarray(fwd["log_fold_change_direction"]) != "neutral" + rev_sig = np.asarray(rev["log_fold_change_direction"]) != "neutral" + assert fwd_sig.sum() > 0 + assert np.array_equal(fwd_sig, rev_sig) + + +# -------------------------------------------------------------------------- +# where the exactness STOPS -- reported, not hidden +# -------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "kwargs, why", + [ + pytest.param( + dict(n_landmarks=40), + "automatic landmarks", + id="automatic-landmarks", + ), + pytest.param( + dict(sync_parameters=True), + "d / mu / ls synchronised from the stacked union", + id="sync-parameters", + ), + ], +) +def test_the_swap_is_only_approximate_when_the_union_order_matters(kwargs, why): + """Characterisation: two configurations break the exactness, for one reason. + + ``np.vstack([X_condition1, X_condition2])`` is row-order dependent. It feeds + ``compute_landmarks`` always, and under ``sync_parameters=True`` it also feeds + ``compute_d_factal`` / ``compute_mu`` / ``compute_ls``. Either way the two + orientations are then computed from differently-ordered unions and agreement + degrades -- silently, which is why it is pinned here. + + This is NOT asserting the gap is small. It asserts the direction of the + finding: exact at the defaults, inexact here. + """ + fwd, rev = _swap_pair(**kwargs) + with pytest.raises(AssertionError): + np.testing.assert_allclose( + np.asarray(fwd["log_fold_change"]), + -np.asarray(rev["log_fold_change"]), + **_TOL, + ) + + # and the defaults, by contrast, ARE exact -- the comparison that gives the + # statement above its meaning + d_fwd, d_rev = _swap_pair() + np.testing.assert_allclose( + np.asarray(d_fwd["log_fold_change"]), -np.asarray(d_rev["log_fold_change"]), **_TOL + ) + + +def test_shared_landmarks_restore_the_exact_swap(): + """The remedy is the same as for differential expression: share the landmarks.""" + from mellon.parameters import compute_landmarks + + X1, X2 = _data() + shared = np.asarray( + compute_landmarks( + np.vstack([X1, X2]), gp_type="fixed", n_landmarks=40, random_state=7 + ) + ) + fwd, rev = _swap_pair(n_landmarks=40, landmarks=shared) + np.testing.assert_allclose( + np.asarray(fwd["log_fold_change"]), -np.asarray(rev["log_fold_change"]), **_TOL + ) + np.testing.assert_allclose( + fwd["neg_log10_fold_change_ptp"], rev["neg_log10_fold_change_ptp"], **_TOL + ) From f72573e7bfd8ad4466c46da94735a46426efcbd7 Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 10:26:29 -0700 Subject: [PATCH 08/23] fix(docs): correct two over-reaching claims about differential abundance Both corrections are the same move: a correct local measurement generalised past the condition that made it true. Found by an independent review pass. CORRECTION A -- "at its defaults" named no entry point, and the entry point users configure through does not have those defaults. anndata/differential_abundance.py:104 takes n_landmarks=None only when gp is None, and settings.py:82 defaults GPSettings.n_landmarks to 5000. So da(adata, gp=GPSettings(...)) lands in the automatic-landmark branch the same paragraph calls approximate. GPSettings being the documented way to configure da(), the exact configuration is arguably the one a user is LEAST likely to be in. Said at that strength rather than softened to a caveat. CORRECTION B -- "pass one landmarks array to both runs to restore it" was applied to BOTH approximate configurations. It is measured for the landmark half only, and it is now measured FALSE for the other: sync_parameters=True with 40 shared landmarks reads 4.784e-03 against a 4.720e-03 reference with no landmarks at all. Indistinguishable, and the reference arm is what makes that a number rather than an impression. mu and ls come from the union's nearest-neighbour distances, which landmarks do not enter. THE REMEDY, WITH ITS SIZE BOUNDARY IN THE SAME SENTENCE, because a boundary in a following clause is one that gets quoted away: pass mu and ls explicitly -- and d as well when the combined data exceeds 500 cells, because compute_d_factal subsamples 500 INDICES above that and becomes row-order dependent. Below 500, mu and ls suffice. Passing all three is always correct. measured: d+mu+ls explicit, 400+400 -> 0.000e+00 mu+ls only, 400+400 -> 6.085e-04 (d 1.4972 vs 1.4983) mu+ls only, 200+200 -> 0.000e+00 AND THE BOUNDARY WAS ONE ARM FROM BEING MISSED. My own earlier note that "d is identical across orientations" was measured at 160 combined cells. True there; false above 500. Had the mu+ls arm been run at a single size, "pass mu and ls" would have shipped as the remedy and been wrong for nearly every real dataset. That is a THIRD instance of the same move in one PR, and the test module docstring now records it beside the boundary so the next reader meets the trap with the fact. Sites corrected: CHANGELOG.md (both), tests/test_da_swap_equivalence.py module docstring and test_shared_landmarks_restore_the_exact_swap (scoped to the landmark half). The PR body and a retraction comment follow; the previously posted comment is not silently editable, so its retraction is stated in the new one rather than edited around. Docstrings and CHANGELOG only -- test file verified AST-identical to HEAD with docstrings stripped. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 32 ++++++++++++++++++++++-------- tests/test_da_swap_equivalence.py | 33 ++++++++++++++++++++++++++----- 2 files changed, 52 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bee9a6..9b0ab24 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,14 +42,30 @@ All notable changes to this project will be documented in this file. - `GPSettings.ls_scheme` and `DifferentialExpression.fit` now state where the shared length scale comes from, that the default makes the contrast depend on argument order, and what the `"condition2"` equivalence requires of landmarks. - - **Differential abundance is swap-equivalent at its defaults**, and now has tests saying so: - `da(X, Y)` and `da(Y, X)` give the same densities with the roles exchanged, negated log fold - changes and z-scores, identical uncertainty, PTP and significant-cell set, and exchanged - `'up'`/`'down'` labels. It needs no `ls_scheme` analogue — abundance builds both density - estimators from one shared configuration and never had the condition-1 inheritance - `ls_scheme` exists to expose. The equivalence is approximate rather than exact under - `sync_parameters=True` or automatic landmarks, both because the stacked union is row-order - dependent; pass one `landmarks` array to both runs to restore it. + - **Differential abundance is swap-equivalent when it uses no landmarks** — the + `DifferentialAbundance` class default (`n_landmarks=None`, `sync_parameters=False`) — and now + has tests saying so: `da(X, Y)` and `da(Y, X)` give the same densities with the roles + exchanged, negated log fold changes and z-scores, identical uncertainty, PTP and + significant-cell set, and exchanged `'up'`/`'down'` labels. It needs no `ls_scheme` analogue — + abundance builds both density estimators from one shared configuration and never had the + condition-1 inheritance `ls_scheme` exists to expose. + + Note which entry point that is: `da()` takes `n_landmarks=None` only when `gp is None`, and any + `GPSettings` carries `n_landmarks=5000`, so `da(adata, gp=GPSettings(...))` lands in the + automatic-landmark branch instead. Since `GPSettings` is the documented way to configure + `da()`, the exact configuration is arguably the one you are least likely to be in — and at that + configuration the uncertainty-derived quantities are not merely approximate, see + settylab/kompot#32, which is filed separately and is not caused by this change. + + The equivalence is approximate rather than exact under `sync_parameters=True` or automatic + landmarks — both because the stacked union is row-order dependent, but at different points, so + they need different remedies. For automatic landmarks, pass one `landmarks` array to both runs; + that restores exactness. It does **not** fix `sync_parameters=True`: `mu` and `ls` come from + the union's nearest-neighbour distances, which landmarks do not enter — measured, sharing + landmarks leaves that half at 4.784e-03 against a 4.720e-03 reference with no landmarks at all. + Pass `mu` and `ls` explicitly — **and `d` as well when the combined data exceeds 500 cells**, + because `compute_d_factal` subsamples 500 indices above that and becomes row-order dependent. + Below 500, `mu` and `ls` suffice. **Passing all three is always correct.** ## [0.8.0] - 2026-07-28 diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index 9239c7a..b3df6fe 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -16,10 +16,28 @@ `ls_scheme` was added for, and no DA counterpart of `"condition2"` is warranted. Two configurations DO break the exactness, and both are covered below: automatic -landmarks, and ``sync_parameters=True``. Both trace to the same construction as -the DE case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent --- but under ``sync_parameters`` it feeds ``d``, ``mu`` and ``ls`` as well as the -landmarks, so DA has that hazard in two places rather than one. +landmarks, and ``sync_parameters=True``. Both trace to the same ROOT as the DE +case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent -- but +they act through DIFFERENT proximate mechanisms, and **a remedy applies at the +proximate level, not the root**: ``compute_landmarks`` (k-means init) for the +first, ``compute_nn_distances`` (which sets ``mu`` and ``ls``) for the second. +So sharing a landmark array fixes the first and cannot reach the second -- +landmarks are not an input to ``compute_nn_distances``, and sharing them leaves +the ``sync_parameters`` half at 4.784e-03 against a 4.720e-03 no-landmark +reference. + +The remedy for that half is to pass ``mu`` and ``ls`` explicitly -- **and ``d`` +as well once the combined data exceeds 500 cells**, because +``compute_d_factal`` subsamples 500 INDICES above that (mellon +``parameters.py``) and so becomes row-order dependent itself. Measured: +``mu``+``ls`` alone is exact at 200+200 and NOT exact at 400+400 (6.085e-04, +with ``d`` reading 1.4972 vs 1.4983); all three explicit is exact at 400+400. +Passing all three is always correct. + +Note the trap in that boundary. An earlier version of this file said ``d`` is +identical across orientations. That was measured at 160 combined cells and is +TRUE THERE -- and false above 500, which is most real datasets. A correct local +measurement generalised past the condition that made it true. """ import numpy as np @@ -204,7 +222,12 @@ def test_the_swap_is_only_approximate_when_the_union_order_matters(kwargs, why): def test_shared_landmarks_restore_the_exact_swap(): - """The remedy is the same as for differential expression: share the landmarks.""" + """Share the landmarks, and the automatic-landmark half is exact again. + + Scope: this covers the LANDMARK half only -- ``sync_parameters`` is left at + its default ``False`` here. The same remedy does not reach the + ``sync_parameters`` half, and no test asserts that it does. + """ from mellon.parameters import compute_landmarks X1, X2 = _data() From c962036c93b47d81bba25d457776dd4894b8b328 Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 12:05:33 -0700 Subject: [PATCH 09/23] docs(da): the two swap remedies compose; move the precondition into the docstrings Skeptic round on the corrections themselves (lsc2sk3). Four fixes, one bundle. (1) "PASSING ALL THREE IS ALWAYS CORRECT" WAS FALSE WITH AUTOMATIC LANDMARKS. It fixes the PARAMETER half and leaves the LANDMARK half untouched. Measured at 400+400: sync + AUTO landmarks + d,mu,ls explicit -> uncertainty 8.4e-01 sync + SHARED landmarks + d,mu,ls explicit -> 0 So the two remedies COMPOSE and neither substitutes for the other. Every site now says so: with sync_parameters=True passing all three is always correct, AND if landmarks are automatic pass one landmarks array as well. This is the fourth instance on this PR of one habit -- a correct local result generalised past the condition that made it true -- and the first where the measurement I over-generalised was not even my own. I transcribed a peer's number, which had been taken with landmarks controlled, and stated it unconditionally. Re-reading it would not have helped; the condition is not visible in the number. (2) The CHANGELOG had grown 49 -> 65 lines against a request to keep it concise. The DA precondition and both remedies now live in DifferentialAbundance.fit, under `sync_parameters` and `landmarks`, which said NOTHING about order dependence before -- so this is a move into the place a user meets it, not a deletion. The bullet is 9 lines and points there; [Unreleased] is back to 50. (3)(4) belong to the PR body: the suite section still called the DA whole-suite run outstanding and carried a spliced duplicate sentence, and the 1.220e-02 table was unlabelled and read as the 400+400 run when it is 120+40 in 3-D. Verified: kompot/differential/differential_abundance.py and tests/test_da_swap_equivalence.py are both AST-identical to f72573e with docstrings stripped. tests/test_da_swap_equivalence.py + tests/test_ls_scheme.py 46 passed. Docstrings and CHANGELOG only, so no full-suite re-run. Refs settylab/dotto-nexus#310 --- CHANGELOG.md | 29 +++++------------ kompot/differential/differential_abundance.py | 31 +++++++++++++++++++ tests/test_da_swap_equivalence.py | 14 +++++++-- 3 files changed, 50 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9b0ab24..1efc011 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,28 +44,13 @@ All notable changes to this project will be documented in this file. `"condition2"` equivalence requires of landmarks. - **Differential abundance is swap-equivalent when it uses no landmarks** — the `DifferentialAbundance` class default (`n_landmarks=None`, `sync_parameters=False`) — and now - has tests saying so: `da(X, Y)` and `da(Y, X)` give the same densities with the roles - exchanged, negated log fold changes and z-scores, identical uncertainty, PTP and - significant-cell set, and exchanged `'up'`/`'down'` labels. It needs no `ls_scheme` analogue — - abundance builds both density estimators from one shared configuration and never had the - condition-1 inheritance `ls_scheme` exists to expose. - - Note which entry point that is: `da()` takes `n_landmarks=None` only when `gp is None`, and any - `GPSettings` carries `n_landmarks=5000`, so `da(adata, gp=GPSettings(...))` lands in the - automatic-landmark branch instead. Since `GPSettings` is the documented way to configure - `da()`, the exact configuration is arguably the one you are least likely to be in — and at that - configuration the uncertainty-derived quantities are not merely approximate, see - settylab/kompot#32, which is filed separately and is not caused by this change. - - The equivalence is approximate rather than exact under `sync_parameters=True` or automatic - landmarks — both because the stacked union is row-order dependent, but at different points, so - they need different remedies. For automatic landmarks, pass one `landmarks` array to both runs; - that restores exactness. It does **not** fix `sync_parameters=True`: `mu` and `ls` come from - the union's nearest-neighbour distances, which landmarks do not enter — measured, sharing - landmarks leaves that half at 4.784e-03 against a 4.720e-03 reference with no landmarks at all. - Pass `mu` and `ls` explicitly — **and `d` as well when the combined data exceeds 500 cells**, - because `compute_d_factal` subsamples 500 indices above that and becomes row-order dependent. - Below 500, `mu` and `ls` suffice. **Passing all three is always correct.** + has tests saying so. It needs no `ls_scheme` analogue: abundance shares one configuration + between both density estimators and never had the condition-1 inheritance `ls_scheme` exists + to expose. `sync_parameters=True` and automatic landmarks each break the equivalence, by + different mechanisms needing different remedies that do not substitute for one another — see + `DifferentialAbundance.fit` for both, and note that `da(adata, gp=GPSettings(...))` carries + `n_landmarks=5000` and so is not at the class default. A pre-existing instability at + `n_landmarks >= n_combined` is tracked at settylab/kompot#32. ## [0.8.0] - 2026-07-28 diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 17ab1ec..f575311 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -176,6 +176,15 @@ def fit( landmarks : np.ndarray, optional Pre-computed landmarks to use. If provided, n_landmarks will be ignored. Shape (n_landmarks, n_features). + + Automatic landmarks are computed from ``np.vstack([X_condition1, + X_condition2])``, whose **row order differs between the two + orientations** of a contrast, so ``da(X, Y)`` and ``da(Y, X)`` select + different landmark sets and their results stop being exact mirrors of + one another. Passing one array here to both runs removes that; so + does ``n_landmarks=None``, which builds no landmarks at all. With + ``random_state=None`` the selection is not even reproducible between + two runs of the same orientation. ls_factor : float, optional Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is not explicitly provided in density_kwargs. @@ -193,6 +202,28 @@ def fit( the combined dataset. When True, parameters are computed once from the combined data to ensure models for both conditions use identical parameter values. This is especially important for consistent density estimation across conditions. Default is False. + + **Order dependence.** The combined dataset is + ``np.vstack([X_condition1, X_condition2])``, so the parameters it + derives depend on the order the conditions were passed in, and + ``da(X, Y)`` and ``da(Y, X)`` stop agreeing. ``mu`` and ``ls`` come + from the union's nearest-neighbour distances; ``d`` additionally + subsamples 500 indices once the **combined** cell count is **above** + 500, and so becomes order-dependent there too. + + To make the two orientations agree, pass ``d``, ``mu`` and ``ls`` + explicitly through ``density_kwargs`` rather than relying on + synchronisation. **Two remedies are needed and neither substitutes + for the other**: passing all three fixes the parameter half, and it + leaves an uncertainty discrepancy behind whenever landmarks are still + chosen automatically — so pass one ``landmarks`` array to both runs + as well (or use ``n_landmarks=None``). Measured on 400 + 400 cells: + all three parameters with automatic landmarks still differs between + orientations, while all three with a shared landmark array agrees + exactly. + + Below 500 combined cells ``mu`` and ``ls`` alone suffice for the + parameter half, but passing ``d`` as well is always safe. **density_kwargs : dict Additional arguments to pass to the DensityEstimator. diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index b3df6fe..9844109 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -32,7 +32,14 @@ ``parameters.py``) and so becomes row-order dependent itself. Measured: ``mu``+``ls`` alone is exact at 200+200 and NOT exact at 400+400 (6.085e-04, with ``d`` reading 1.4972 vs 1.4983); all three explicit is exact at 400+400. -Passing all three is always correct. +With ``sync_parameters=True``, passing all three is always correct. + +**THE TWO REMEDIES COMPOSE AND NEITHER SUBSTITUTES FOR THE OTHER.** Passing +``d``, ``mu`` and ``ls`` fixes the PARAMETER half and leaves the LANDMARK half +untouched: measured at 400+400, sync with automatic landmarks and all three +parameters explicit still differs between orientations on the uncertainty +(8.4e-01), while the same run with one shared landmark array agrees exactly. +So if landmarks are automatic, pass one ``landmarks`` array as well. Note the trap in that boundary. An earlier version of this file said ``d`` is identical across orientations. That was measured at 160 combined cells and is @@ -226,7 +233,10 @@ def test_shared_landmarks_restore_the_exact_swap(): Scope: this covers the LANDMARK half only -- ``sync_parameters`` is left at its default ``False`` here. The same remedy does not reach the - ``sync_parameters`` half, and no test asserts that it does. + ``sync_parameters`` half, and no test asserts that it does. Nor does the + parameter remedy reach THIS half: the two compose, and a run with + ``sync_parameters=True`` needs both a shared landmark array and explicit + ``d``/``mu``/``ls`` before the two orientations agree. """ from mellon.parameters import compute_landmarks From 9fc37948e68d520961957cfaa9debb84587ba209 Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 12:58:35 -0700 Subject: [PATCH 10/23] feat: one GPSettings.param_scheme for de() and da(); remove ls_scheme Replaces the unreleased `ls_scheme` and the shipped `sync_parameters` with a single `GPSettings.param_scheme` ("condition1" | "condition2" | "symmetric" | "pooled" | "separate"), read by both entry points. The model decides what it covers: `ls` in differential expression, `d`, `mu` and `ls` in differential abundance. `None` keeps each entry point's default ("condition1" for de, "separate" for da), so no existing call changes its numbers. - `ls_scheme` is removed outright (never released): parameter, field, CLI key, template, docs, CHANGELOG, tests (test_ls_scheme.py -> test_param_scheme.py). - `sync_parameters` is a deprecated alias on DifferentialAbundance.fit and a named, deprecated parameter on da() (it used to arrive only via **density_kwargs): True -> "pooled", False -> "separate"; both at once raises. - da() gains "condition1", "condition2" and "symmetric"; one-sided schemes share exactly the values that condition's estimator derives on its own, and "symmetric" is bit-identical under a swap by construction. - Fix: da() silently ignored GPSettings.ls. It now honours ls and param_scheme, warns on the expression-only sigma/eps/use_empirical_variance, and the DA CLI routes ls and param_scheme into GPSettings. --- CHANGELOG.md | 65 ++- docs/source/cli.rst | 4 +- kompot/anndata/_da_helpers.py | 11 + kompot/anndata/_de_helpers.py | 20 +- kompot/anndata/differential_abundance.py | 70 +++- kompot/anndata/differential_expression.py | 12 +- kompot/anndata/utils/params.py | 2 +- kompot/cli/da.py | 27 +- kompot/cli/de.py | 2 +- kompot/cli/templates/da_config_template.yaml | 4 +- kompot/cli/templates/de_config_template.yaml | 2 +- kompot/differential/differential_abundance.py | 202 ++++++++-- .../differential/differential_expression.py | 78 ++-- kompot/settings.py | 87 ++-- tests/test_da_swap_equivalence.py | 43 +- tests/test_differential.py | 10 +- tests/test_differential_extended.py | 6 +- ...test_ls_scheme.py => test_param_scheme.py} | 376 ++++++++++++++++-- 18 files changed, 792 insertions(+), 229 deletions(-) rename tests/{test_ls_scheme.py => test_param_scheme.py} (57%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1efc011..113dab4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,51 +6,44 @@ All notable changes to this project will be documented in this file. ### Added - - **`GPSettings.ls_scheme`**: chooses which cells the automatic shared length scale is estimated - from in differential expression — `"condition1"` (the default, exactly the existing behaviour), - `"condition2"`, `"symmetric"`, `"pooled"`, `"separate"`. The shared value has always been - estimated from **condition 1's cells only**, making it a function of `n_condition1`, so - `de(condition1=X, condition2=Y)` is *not* equivalent to `de(condition1=Y, condition2=X)`; on an - exchangeable null the two orientations differ in false-positive rate by 0.06–0.08 beyond 2:1. - `"symmetric"` and `"pooled"` are swap-invariant, `"separate"` is diagnostics-only. See - `DifferentialExpression.fit` for which to pick and why. - - **`ls_scheme="condition2"`**: the default's mirror — estimated from condition 2's cells — for - diagnosis rather than analysis. `de(X, Y, ls_scheme="condition1")` and - `de(Y, X, ls_scheme="condition2")` are the same computation with the labels exchanged, so the - default's swap-dependence is visible from a single call site. The equivalence is exact only - when both runs use the same landmarks; automatic landmarks are order-dependent, and `fit()` - warns when the scheme is used with them. Measurements and the precondition: - `DifferentialExpression.fit`. + - **`GPSettings.param_scheme`**: one field, read by both `kompot.de()` and `kompot.da()`, that + says where the hyperparameters estimated from cells come from — `"condition1"`, + `"condition2"`, `"symmetric"`, `"pooled"` or `"separate"`. It covers `ls` in differential + expression and `d`, `mu` and `ls` in differential abundance; an explicit value pins that + parameter. The default, `None`, keeps each entry point's existing behaviour: `"condition1"` for + `de()` and `"separate"` for `da()`. + - **Expression results depend on argument order at the default.** The shared length scale is + estimated from condition 1's cells only, so `de(condition1=X, condition2=Y)` is not equivalent + to `de(condition1=Y, condition2=X)`; on an exchangeable null the two orientations differ in + false-positive rate by 0.06–0.08 beyond 2:1. `"symmetric"` and `"pooled"` are swap-invariant, + `"separate"` is for diagnostics, and `"condition2"` mirrors the default so its order dependence + is visible from one call site. See `DifferentialExpression.fit`. + - **Abundance gains `"condition1"`, `"condition2"` and `"symmetric"`.** `"symmetric"` is + swap-invariant by construction. See `DifferentialAbundance.fit`. ### Changed - - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme is - no longer treated as a matching rerun. + - `param_scheme` participates in run-parameter matching for both `de()` and `da()`. + +### Deprecated + + - **`sync_parameters`** on `DifferentialAbundance.fit` and `da()`: `True` is + `param_scheme="pooled"` and `False` is `"separate"`, with identical results. It previously + reached `da()` only through `**density_kwargs`. ### Fixed - - **`ls_scheme` is now readable back out of a stored `run_info`.** It is written nested under - `params["gp"]` but was missing from `params_get`'s key map, so it read back as `None` and - *every* rerun — including a byte-identical one — compared unequal and was reported as a - parameter change. - - **The DE CLI now routes a flat `ls_scheme` config key into `GPSettings`.** An unrecognised key - is not rejected: it falls through to `**function_kwargs` and is handed to mellon, so the - documented knob was unreachable from a config file and the omission was silent. + - **`da()` now honours `GPSettings.ls`.** It read six `GPSettings` fields and silently ignored the + rest. The expression-only fields `sigma`, `eps` and `use_empirical_variance` now log a warning + when set, and the DA CLI routes `ls` and `param_scheme` config keys into `GPSettings`. ### Documentation - - `GPSettings.ls_scheme` and `DifferentialExpression.fit` now state where the shared length scale - comes from, that the default makes the contrast depend on argument order, and what the - `"condition2"` equivalence requires of landmarks. - - **Differential abundance is swap-equivalent when it uses no landmarks** — the - `DifferentialAbundance` class default (`n_landmarks=None`, `sync_parameters=False`) — and now - has tests saying so. It needs no `ls_scheme` analogue: abundance shares one configuration - between both density estimators and never had the condition-1 inheritance `ls_scheme` exists - to expose. `sync_parameters=True` and automatic landmarks each break the equivalence, by - different mechanisms needing different remedies that do not substitute for one another — see - `DifferentialAbundance.fit` for both, and note that `da(adata, gp=GPSettings(...))` carries - `n_landmarks=5000` and so is not at the class default. A pre-existing instability at - `n_landmarks >= n_combined` is tracked at settylab/kompot#32. + - Abundance is swap-equivalent at `"separate"` without landmarks, and tested. `"pooled"` and + automatic landmarks each break that, and need different remedies that do not substitute for + one another: see `DifferentialAbundance.fit`. `da(adata, gp=GPSettings(...))` carries + `n_landmarks=5000`. A pre-existing instability at `n_landmarks >= n_combined` is tracked at + settylab/kompot#32. ## [0.8.0] - 2026-07-28 diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 975f3b8..1f84983 100644 --- a/docs/source/cli.rst +++ b/docs/source/cli.rst @@ -443,9 +443,9 @@ The table below shows how config keys correspond to Python Settings: * - Config key - Python equivalent - Description - * - ``sigma``, ``ls``, ``ls_factor``, ``ls_scheme`` (DE only), ``n_landmarks``, ``batch_size``, ``eps``, ``jit_compile``, ``random_state``, ``use_empirical_variance`` + * - ``ls``, ``ls_factor``, ``param_scheme``, ``n_landmarks``, ``batch_size``, ``jit_compile``, ``random_state``; DE also ``sigma``, ``eps``, ``use_empirical_variance`` - ``gp=GPSettings(...)`` - - GP model parameters + - GP model parameters. ``param_scheme`` defaults to ``condition1`` in DE and ``separate`` in DA * - ``null_genes``, ``null_seed``, ``fdr_threshold`` - ``fdr=FDRSettings(...)`` - FDR / null distribution (DE) diff --git a/kompot/anndata/_da_helpers.py b/kompot/anndata/_da_helpers.py index 4c61b6e..ba9b263 100644 --- a/kompot/anndata/_da_helpers.py +++ b/kompot/anndata/_da_helpers.py @@ -25,6 +25,7 @@ def _check_da_overwrites( condition2: str, obsm_key: str, ls_factor: float, + param_scheme: str = "separate", ): """Check for existing DA results and handle overwrite logic.""" from .utils import detect_output_field_overwrite @@ -101,9 +102,19 @@ def _check_da_overwrites( "condition2", "obsm_key", "ls_factor", + "param_scheme", ]: curr_val = locals().get(param_name) prev_val = _pg(prev_params, param_name) + if param_name == "param_scheme" and prev_val is None: + # Recorded before `param_scheme` existed: sharing was + # then only reachable as a `sync_parameters` kwarg. + prev_extra = prev_params.get("extra_kwargs") or {} + prev_val = ( + "pooled" + if prev_extra.get("sync_parameters") + else "separate" + ) if curr_val != prev_val: params_match = False logger.debug( diff --git a/kompot/anndata/_de_helpers.py b/kompot/anndata/_de_helpers.py index 6548d13..f17190a 100644 --- a/kompot/anndata/_de_helpers.py +++ b/kompot/anndata/_de_helpers.py @@ -41,7 +41,7 @@ def _check_overwrites( obsm_key: str, layer: Optional[str], ls_factor: float, - ls_scheme: str = "condition1", + param_scheme: str = "condition1", ): """Check for existing results and handle overwrite logic. @@ -135,10 +135,15 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", - "ls_scheme", + "param_scheme", ]: curr_val = locals().get(param_name) - prev_val = _pg(prev_params, param_name) + # A run recorded before `param_scheme` existed used "condition1". + prev_val = _pg( + prev_params, + param_name, + "condition1" if param_name == "param_scheme" else None, + ) if curr_val != prev_val: params_match = False logger.debug( @@ -191,10 +196,15 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", - "ls_scheme", + "param_scheme", ]: curr_val = locals().get(param_name) - prev_val = _pg(prev_params, param_name) + # A run recorded before `param_scheme` existed used "condition1". + prev_val = _pg( + prev_params, + param_name, + "condition1" if param_name == "param_scheme" else None, + ) if curr_val != prev_val: params_match = False diff --git a/kompot/anndata/differential_abundance.py b/kompot/anndata/differential_abundance.py index a2de64a..d8f1872 100644 --- a/kompot/anndata/differential_abundance.py +++ b/kompot/anndata/differential_abundance.py @@ -9,7 +9,9 @@ from typing import Union, Dict, Any from ..differential import DifferentialAbundance +from ..differential.differential_abundance import DEFAULT_PARAM_SCHEME from ..settings import ( + resolve_param_scheme, GPSettings, DAThresholdSettings, StorageSettings, @@ -48,6 +50,7 @@ def da( storage: "StorageSettings | None" = None, output: "OutputSettings | None" = None, model: "ModelSettings | None" = None, + sync_parameters: "bool | None" = None, **density_kwargs, ) -> "Union[Dict[str, np.ndarray], Any]": """Run differential abundance analysis on an AnnData object. @@ -76,9 +79,16 @@ def da( sample_col : str, optional Column with biological-replicate labels. gp : GPSettings, optional - GP model parameters (``ls_factor``, ``n_landmarks``, - ``landmarks``, ``batch_size``, ``jit_compile``, - ``random_state``). + GP model parameters. Used here: ``ls``, ``ls_factor``, + ``param_scheme``, ``n_landmarks``, ``landmarks``, ``batch_size``, + ``jit_compile`` and ``random_state``. ``param_scheme`` decides + whether the two density estimators share ``d``, ``mu`` and ``ls`` + and where the shared values come from; ``None`` means + ``"separate"``. See + :meth:`~kompot.differential.DifferentialAbundance.fit`. + ``sigma``, ``eps`` and ``use_empirical_variance`` are + expression-only and are ignored with a warning when set to + anything but their defaults. threshold : DAThresholdSettings, optional Significance thresholds for abundance changes. storage : StorageSettings, optional @@ -89,8 +99,14 @@ def da( Pre-fitted models or predictors to inject. For DA, only ``density_predictor1/2`` and ``variance_predictor1/2`` are used. See :class:`~kompot.ModelSettings`. + sync_parameters : bool, optional + .. deprecated:: + Use ``gp=GPSettings(param_scheme=...)``. ``True`` is + ``"pooled"`` and ``False`` is ``"separate"``; passing both + raises. **density_kwargs - Forwarded to :class:`~mellon.DensityEstimator`. + Forwarded to :class:`~mellon.DensityEstimator`. ``d``, ``mu`` and + ``ls`` given here pin that parameter under any ``param_scheme``. Returns ------- @@ -107,6 +123,42 @@ def da( batch_size = gp.batch_size if gp is not None else None jit_compile = gp.jit_compile if gp is not None else False random_state = gp.random_state if gp is not None else None + param_scheme = gp.param_scheme if gp is not None else None + + if gp is not None: + _unused = [ + name + for name in ("sigma", "eps", "use_empirical_variance") + if getattr(gp, name) != getattr(GPSettings, name) + ] + if _unused: + logger.warning( + f"GPSettings field(s) {_unused} have no effect in da(): they " + "configure the expression model only." + ) + if gp.ls is not None: + if "ls" in density_kwargs: + raise ValueError( + "`ls` given both as GPSettings.ls and as a density keyword " + "argument; pass it once." + ) + density_kwargs["ls"] = gp.ls + + if sync_parameters is not None: + if param_scheme is not None: + raise ValueError( + "Pass either GPSettings.param_scheme or the deprecated " + "`sync_parameters`, not both." + ) + warnings.warn( + "`sync_parameters` is deprecated; use " + "gp=GPSettings(param_scheme='pooled') for True or " + "param_scheme='separate' for False.", + DeprecationWarning, + stacklevel=2, + ) + param_scheme = "pooled" if sync_parameters else "separate" + param_scheme = resolve_param_scheme(param_scheme, DEFAULT_PARAM_SCHEME) _threshold = threshold if threshold is not None else DAThresholdSettings() log_fold_change_threshold = _threshold.lfc_threshold @@ -157,6 +209,7 @@ def da( condition2=condition2, obsm_key=obsm_key, ls_factor=ls_factor, + param_scheme=param_scheme, ) # ---- 1. Copy if requested ---- @@ -199,6 +252,7 @@ def da( condition1_sample_indices=data["condition1_sample_indices"], condition2_sample_indices=data["condition2_sample_indices"], allow_single_condition_variance=allow_single_condition_variance, + param_scheme=param_scheme, **density_kwargs, ) @@ -235,7 +289,9 @@ def da( sample_col=sample_col, ), gp=GPSettings( + ls=gp.ls if gp is not None else None, ls_factor=ls_factor, + param_scheme=param_scheme, n_landmarks=n_landmarks, batch_size=batch_size, jit_compile=jit_compile, @@ -256,7 +312,11 @@ def da( allow_single_condition_variance=allow_single_condition_variance, progress=progress, ), - extra_kwargs=density_kwargs, + extra_kwargs={ + k: v + for k, v in density_kwargs.items() + if not (k == "ls" and gp is not None and gp.ls is not None) + }, landmarks_provided=landmarks is not None, ) diff --git a/kompot/anndata/differential_expression.py b/kompot/anndata/differential_expression.py index d0a833b..bcd6e88 100644 --- a/kompot/anndata/differential_expression.py +++ b/kompot/anndata/differential_expression.py @@ -9,6 +9,8 @@ from typing import Union, Dict, Any from ..differential import DifferentialExpression +from ..differential.differential_expression import DEFAULT_PARAM_SCHEME +from ..settings import resolve_param_scheme from ..settings import ( GPSettings, FDRSettings, @@ -132,7 +134,7 @@ def de( sigma = _gp.sigma ls = _gp.ls ls_factor = _gp.ls_factor - ls_scheme = _gp.ls_scheme + param_scheme = resolve_param_scheme(_gp.param_scheme, DEFAULT_PARAM_SCHEME) n_landmarks = _gp.n_landmarks landmarks = _gp.landmarks use_empirical_variance = _gp.use_empirical_variance @@ -306,7 +308,7 @@ def de( obsm_key=obsm_key, layer=layer, ls_factor=ls_factor, - ls_scheme=ls_scheme, + param_scheme=param_scheme, ) # ---- 2. Copy if requested ---- @@ -457,7 +459,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, - ls_scheme=ls_scheme, + param_scheme=param_scheme, landmarks=landmarks, condition1_sample_indices=condition1_sample_indices, condition2_sample_indices=condition2_sample_indices, @@ -638,7 +640,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, - ls_scheme=ls_scheme, + param_scheme=param_scheme, n_landmarks=n_landmarks, use_empirical_variance=use_empirical_variance, batch_size=batch_size, @@ -734,7 +736,6 @@ def compute_differential_expression( sigma: float = 1.0, ls=None, ls_factor: float = 10.0, - ls_scheme: str = "condition1", compute_mahalanobis: bool = True, jit_compile: bool = False, eps: float = 1e-8, @@ -788,7 +789,6 @@ def compute_differential_expression( sigma=sigma, ls=ls, ls_factor=ls_factor, - ls_scheme=ls_scheme, n_landmarks=n_landmarks, landmarks=landmarks, use_empirical_variance=use_empirical_variance, diff --git a/kompot/anndata/utils/params.py b/kompot/anndata/utils/params.py index 6732f14..545a434 100644 --- a/kompot/anndata/utils/params.py +++ b/kompot/anndata/utils/params.py @@ -197,7 +197,7 @@ def params_to_call_args( "sigma": ("gp", "sigma"), "ls": ("gp", "ls"), "ls_factor": ("gp", "ls_factor"), - "ls_scheme": ("gp", "ls_scheme"), + "param_scheme": ("gp", "param_scheme"), "n_landmarks": ("gp", "n_landmarks"), "use_empirical_variance": ("gp", "use_empirical_variance"), "batch_size": ("gp", "batch_size"), diff --git a/kompot/cli/da.py b/kompot/cli/da.py index aa77da8..2c82f47 100644 --- a/kompot/cli/da.py +++ b/kompot/cli/da.py @@ -184,6 +184,23 @@ def add_da_parser(subparsers) -> argparse.ArgumentParser: return parser +#: Flat config keys the DA CLI routes into :class:`~kompot.settings.GPSettings`. +#: Every other key -- including ``sigma``, which here means mellon's density +#: noise rather than the expression model's -- is forwarded to ``da()`` and +#: from there to :class:`mellon.DensityEstimator`. The deprecated +#: ``sync_parameters`` key reaches ``da()``'s named, deprecated parameter. +GP_CONFIG_KEYS = { + "ls", + "ls_factor", + "param_scheme", + "n_landmarks", + "landmarks", + "batch_size", + "jit_compile", + "random_state", +} + + def run_da(args): """ Run differential abundance analysis. @@ -292,15 +309,7 @@ def run_da(args): sample_col = params.pop("sample_col", None) # Build Settings from remaining params - gp_keys = { - "n_landmarks", - "landmarks", - "ls_factor", - "batch_size", - "jit_compile", - "random_state", - } - gp_kwargs = {k: params.pop(k) for k in list(params) if k in gp_keys} + gp_kwargs = {k: params.pop(k) for k in list(params) if k in GP_CONFIG_KEYS} gp = GPSettings(**gp_kwargs) if gp_kwargs else None threshold_kwargs = {} diff --git a/kompot/cli/de.py b/kompot/cli/de.py index 9a865ae..5ee951e 100644 --- a/kompot/cli/de.py +++ b/kompot/cli/de.py @@ -213,7 +213,7 @@ def add_de_parser(subparsers) -> argparse.ArgumentParser: "sigma", "ls", "ls_factor", - "ls_scheme", + "param_scheme", "n_landmarks", "use_empirical_variance", "batch_size", diff --git a/kompot/cli/templates/da_config_template.yaml b/kompot/cli/templates/da_config_template.yaml index df1febd..275db47 100644 --- a/kompot/cli/templates/da_config_template.yaml +++ b/kompot/cli/templates/da_config_template.yaml @@ -33,7 +33,10 @@ random_state: null # Random seed for landmark selection jit_compile: false # Use JAX JIT compilation # GP kernel parameters: +ls: null # Length scale (null = auto-estimate) ls_factor: 10.0 # Length scale multiplication factor +param_scheme: null # Where d, mu and ls are estimated from (null = separate): + # separate | pooled | condition1 | condition2 | symmetric # Memory management: store_arrays_on_disk: null # Store large arrays on disk (null = auto) @@ -57,7 +60,6 @@ return_full_results: false # Return full results dict (not used in CLI) # # d_method: "fractal" # Density estimation method # sigma: null # Noise level (null = auto) -# ls: null # Length scale (null = auto) # gp_type: "vonkarman" # GP kernel type # cov_func_curry: null # Custom covariance function # rank: null # Rank for low-rank approximation diff --git a/kompot/cli/templates/de_config_template.yaml b/kompot/cli/templates/de_config_template.yaml index c2d1f18..6697c4a 100644 --- a/kompot/cli/templates/de_config_template.yaml +++ b/kompot/cli/templates/de_config_template.yaml @@ -43,7 +43,7 @@ jit_compile: false # Use JAX JIT compilation sigma: 1.0 # Noise level for function estimator ls: null # Length scale (null = auto-estimate) ls_factor: 10.0 # Length scale multiplication factor -ls_scheme: condition1 # Which cells the shared length scale is estimated from: +param_scheme: null # Where the shared length scale is estimated from (null = condition1): # condition1 | condition2 | symmetric | pooled | separate eps: 1.0e-8 # Numerical stability constant diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index f575311..cf0ad7c 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -3,16 +3,90 @@ import numpy as np from typing import Optional, Dict, Any import logging +import warnings from scipy.stats import norm as normal import mellon from mellon.parameters import compute_landmarks from ..batch_utils import apply_batched +from ..settings import PARAM_SCHEMES, resolve_param_scheme from .sample_variance_estimator import SampleVarianceEstimator logger = logging.getLogger("kompot") +#: The ``param_scheme`` an abundance fit uses when none is given. +DEFAULT_PARAM_SCHEME = "separate" + +#: The density hyperparameters a ``param_scheme`` covers. +_DENSITY_PARAMS = ("d", "mu", "ls") + +# ``compute_d_factal`` averages local dimensionality over at most this many +# query cells (mellon ``parameters.py``); the ``"symmetric"`` scheme weights +# each condition's ``d`` by the number of cells it actually averaged over. +_D_FACTAL_QUERY_CELLS = 500 + + +def _condition_density_params(X, ls_factor, seed): + """``d``, ``mu`` and ``ls`` exactly as a ``mellon.DensityEstimator`` derives them. + + Mirrors ``DensityEstimator._compute_d`` (``d_method="fractal"``), + ``_compute_nn_distances``, ``_compute_mu`` and ``_compute_ls``, so that a + value estimated here *before* either estimator is fitted is the value the + estimator would have reached on its own. Returns the validated + nearest-neighbour distances as well, for ``"symmetric"``. + """ + from mellon.validation import validate_nn_distances + + X = np.asarray(X) + nn = validate_nn_distances(mellon.parameters.compute_nn_distances(X, seed=seed)) + d = mellon.parameters.compute_d_factal(X) + mu = mellon.parameters.compute_mu(nn, d) + ls = mellon.parameters.compute_ls(nn) * ls_factor + return {"d": d, "mu": mu, "ls": ls}, nn + + +def _resolve_scheme_density_params( + param_scheme, X_condition1, X_condition2, ls_factor, seed +): + """Return the shared ``{d, mu, ls}`` a one-sided or symmetric scheme prescribes. + + ``"separate"`` and ``"pooled"`` are handled by the caller; this covers the + three schemes that estimate from each condition on its own. + """ + if param_scheme == "condition1": + return _condition_density_params(X_condition1, ls_factor, seed)[0] + if param_scheme == "condition2": + return _condition_density_params(X_condition2, ls_factor, seed)[0] + if param_scheme == "symmetric": + # Each shared value is the model's own estimator applied to the two + # conditions' WITHIN-condition statistics pooled together -- never to + # the stacked union, whose row order depends on the orientation. Every + # combination below is written so that exchanging the conditions + # exchanges two commutative terms, which makes the result bit-identical + # under a swap by construction, not merely close. + p1, nn1 = _condition_density_params(X_condition1, ls_factor, seed) + p2, nn2 = _condition_density_params(X_condition2, ls_factor, seed) + n1, n2 = np.asarray(X_condition1).shape[0], np.asarray(X_condition2).shape[0] + # ls: geometric mean of nn distances pooled == size-weighted geometric + # mean of the per-condition values (the differential-expression rule). + ls = float( + np.exp((n1 * np.log(p1["ls"]) + n2 * np.log(p2["ls"])) / (n1 + n2)) + ) + # d: mean local dimensionality pooled over the query cells each + # condition's estimate averaged. + q1 = min(n1, _D_FACTAL_QUERY_CELLS) + q2 = min(n2, _D_FACTAL_QUERY_CELLS) + d = float((q1 * p1["d"] + q2 * p2["d"]) / (q1 + q2)) + # mu: the 1st-percentile rule over the pooled within-condition nn + # distances, at the shared d. A quantile sorts, so the order the two + # conditions are concatenated in does not reach the result. + mu = mellon.parameters.compute_mu(np.concatenate([nn1, nn2]), d) + return {"d": d, "mu": mu, "ls": ls} + raise ValueError( + f"Unknown param_scheme {param_scheme!r}. Expected one of {PARAM_SCHEMES}." + ) + class DifferentialAbundance: """ @@ -21,9 +95,9 @@ class DifferentialAbundance: This class analyzes the differences in cell density between two conditions (e.g., control to treatment) using density estimation and fold change analysis. - The analysis can be performed with synchronized parameters between conditions - by setting sync_parameters=True in the fit method, which ensures consistent - density estimation across both conditions. + Whether the two conditions' density estimators share their ``d``, ``mu`` + and ``ls``, and which cells the shared values come from, is set by + ``param_scheme`` in :meth:`fit`. Attributes ---------- @@ -44,8 +118,8 @@ class DifferentialAbundance: Methods ------- - fit(X_condition1, X_condition2, sync_parameters=False, **density_kwargs) - Fit density estimators for both conditions, optionally with synchronized parameters. + fit(X_condition1, X_condition2, param_scheme=None, **density_kwargs) + Fit density estimators for both conditions, optionally sharing parameters. predict(X_new) Predict log density and log fold change for new points. """ @@ -157,8 +231,9 @@ def fit( condition1_sample_indices: Optional[np.ndarray] = None, condition2_sample_indices: Optional[np.ndarray] = None, sample_estimator_ls: Optional[float] = None, - sync_parameters: bool = False, + sync_parameters: Optional[bool] = None, allow_single_condition_variance: bool = False, + param_scheme: Optional[str] = None, **density_kwargs, ): """ @@ -183,8 +258,9 @@ def fit( different landmark sets and their results stop being exact mirrors of one another. Passing one array here to both runs removes that; so does ``n_landmarks=None``, which builds no landmarks at all. With - ``random_state=None`` the selection is not even reproducible between - two runs of the same orientation. + ``random_state=None`` (the default) the selection is not even + reproducible between two runs of the same orientation. This is + independent of ``param_scheme``: see there for the other half. ls_factor : float, optional Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is not explicitly provided in density_kwargs. @@ -198,12 +274,38 @@ def fit( Length scale for the sample-specific variance estimators. If None, will use the same value as ls or it will be estimated, by default None. sync_parameters : bool, optional - Whether to synchronize model parameters (d, mu, ls) between both conditions using - the combined dataset. When True, parameters are computed once from the combined data - to ensure models for both conditions use identical parameter values. This is especially - important for consistent density estimation across conditions. Default is False. - - **Order dependence.** The combined dataset is + .. deprecated:: + Use ``param_scheme``. ``True`` is ``param_scheme="pooled"`` and + ``False`` is ``param_scheme="separate"``; passing both raises. + allow_single_condition_variance : bool, optional + Allow sample-variance estimation when only one condition has + multiple samples, by default False. + param_scheme : str, optional + Where the density hyperparameters ``d``, ``mu`` and ``ls`` come from. + The same field drives :meth:`DifferentialExpression.fit`, where it + covers ``ls`` only. A value passed explicitly through + ``density_kwargs`` (``d=``, ``mu=``, ``ls=``) pins that parameter + and takes precedence. ``None`` (default) means ``"separate"``. + + * ``"separate"`` (default) — each density estimator derives its own + ``d``, ``mu`` and ``ls`` from its own condition; nothing is shared. + * ``"pooled"`` — estimate all three once from both conditions' cells + taken together and share them (formerly ``sync_parameters=True``). + * ``"condition1"`` / ``"condition2"`` — estimate all three from one + condition's cells, exactly as that condition's own estimator + would, and share them. As asymmetric as the names say. + * ``"symmetric"`` — estimate from each condition separately and share + the pooled estimate: ``ls`` is the size-weighted geometric mean of + the two (as in differential expression), ``d`` the mean local + dimensionality over both conditions' query cells, and ``mu`` the + 1st-percentile rule over both conditions' within-condition + nearest-neighbour distances at that ``d``. Bit-identical under + swapping the conditions. + + A shared ``d`` is always the fractal estimate; pass ``d`` to use + another. + + **Order dependence.** ``"pooled"`` works on ``np.vstack([X_condition1, X_condition2])``, so the parameters it derives depend on the order the conditions were passed in, and ``da(X, Y)`` and ``da(Y, X)`` stop agreeing. ``mu`` and ``ls`` come @@ -211,19 +313,23 @@ def fit( subsamples 500 indices once the **combined** cell count is **above** 500, and so becomes order-dependent there too. - To make the two orientations agree, pass ``d``, ``mu`` and ``ls`` - explicitly through ``density_kwargs`` rather than relying on - synchronisation. **Two remedies are needed and neither substitutes - for the other**: passing all three fixes the parameter half, and it - leaves an uncertainty discrepancy behind whenever landmarks are still - chosen automatically — so pass one ``landmarks`` array to both runs - as well (or use ``n_landmarks=None``). Measured on 400 + 400 cells: - all three parameters with automatic landmarks still differs between + To make the two orientations agree under ``"pooled"``, pass ``d``, + ``mu`` and ``ls`` explicitly through ``density_kwargs``. **Two + remedies are needed and neither substitutes for the other**: passing + all three fixes the parameter half, and it leaves an uncertainty + discrepancy behind whenever landmarks are still chosen + automatically — so pass one ``landmarks`` array to both runs as well + (or use ``n_landmarks=None``). Measured on 400 + 400 cells: all + three parameters with automatic landmarks still differs between orientations, while all three with a shared landmark array agrees - exactly. - - Below 500 combined cells ``mu`` and ``ls`` alone suffice for the - parameter half, but passing ``d`` as well is always safe. + exactly. Below 500 combined cells ``mu`` and ``ls`` alone suffice + for the parameter half, but passing ``d`` as well is always safe. + + ``"separate"``, ``"symmetric"`` and the ``"condition1"`` / + ``"condition2"`` mirror never touch the union, so only the landmark + half applies to them: ``da(X, Y, "condition1")`` and + ``da(Y, X, "condition2")`` are the same computation with the labels + exchanged when both runs use the same landmarks. **density_kwargs : dict Additional arguments to pass to the DensityEstimator. @@ -233,6 +339,22 @@ def fit( The fitted instance. """ + if sync_parameters is not None: + if param_scheme is not None: + raise ValueError( + "Pass either `param_scheme` or the deprecated " + "`sync_parameters`, not both." + ) + warnings.warn( + "`sync_parameters` is deprecated; use " + "`param_scheme='pooled'` (for True) or " + "`param_scheme='separate'` (for False).", + DeprecationWarning, + stacklevel=2, + ) + param_scheme = "pooled" if sync_parameters else "separate" + param_scheme = resolve_param_scheme(param_scheme, DEFAULT_PARAM_SCHEME) + # Create or use density predictors if self.density_predictor1 is None or self.density_predictor2 is None: # Configure density estimator defaults @@ -277,8 +399,9 @@ def fit( # Store computed landmarks for future use self.computed_landmarks = computed_landmarks - # Handle synchronization of parameters when sync_parameters is True - if sync_parameters: + # "pooled": estimate d / mu / ls from the stacked union and share. + # This block is the former `sync_parameters=True` path, unchanged. + if param_scheme == "pooled": # Combine data from both conditions for parameter estimation X_combined = np.vstack([X_condition1, X_condition2]) logger.info( @@ -309,6 +432,29 @@ def fit( estimator_defaults["ls"] = ls logger.info(f"Synchronizing parameter ls to {ls:.4f}") + # "condition1" / "condition2" / "symmetric": estimate from each + # condition on its own, before either estimator is fitted, and share. + elif param_scheme != "separate": + unpinned = [p for p in _DENSITY_PARAMS if p not in density_kwargs] + if unpinned: + # the seed a DensityEstimator uses for its own nn search + from mellon.parameters import DEFAULT_RANDOM_SEED + + seed = density_kwargs.get("random_state") + if seed is None: + seed = DEFAULT_RANDOM_SEED + shared = _resolve_scheme_density_params( + param_scheme, X_condition1, X_condition2, ls_factor, seed + ) + for name in unpinned: + estimator_defaults[name] = shared[name] + logger.info( + f"param_scheme={param_scheme!r}: sharing " + + ", ".join( + f"{name}={float(shared[name]):.4f}" for name in unpinned + ) + ) + # Fit density estimators for both conditions logger.info("Fitting density estimator for condition 1...") density_estimator_condition1 = mellon.DensityEstimator(**estimator_defaults) diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index 5e11b4a..b9ef5bd 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -10,12 +10,13 @@ from ..utils import compute_mahalanobis_distances from ..batch_utils import apply_batched, is_jax_memory_error # noqa: F401 +from ..settings import PARAM_SCHEMES, resolve_param_scheme from .expression_model import ExpressionModel logger = logging.getLogger("kompot") -#: Recognised values for ``ls_scheme``. See :meth:`DifferentialExpression.fit`. -LS_SCHEMES = ("condition1", "condition2", "symmetric", "pooled", "separate") +#: The ``param_scheme`` an expression fit uses when none is given. +DEFAULT_PARAM_SCHEME = "condition1" def _auto_ls(X: np.ndarray, ls_factor: float) -> float: @@ -30,16 +31,20 @@ def _auto_ls(X: np.ndarray, ls_factor: float) -> float: return float(compute_ls(compute_nn_distances(np.asarray(X))) * ls_factor) -def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor): - """Return the ``(ls_model1, ls_model2)`` a scheme prescribes. +def _resolve_scheme_ls(param_scheme, X_condition1, X_condition2, ls_factor): + """Return the ``(ls_model1, ls_model2)`` a ``param_scheme`` prescribes. + + ``ls`` is the only hyperparameter a :class:`mellon.FunctionEstimator` + estimates from cells (it has no ``d``, and its ``mu`` is a constant prior + mean), so in differential expression the scheme covers ``ls`` alone. ``None`` in a slot means "let that model estimate its own". For ``"condition1"`` both slots are None and the caller reproduces the historical behaviour by copying model1's fitted value into model2. """ - if ls_scheme == "condition1": + if param_scheme == "condition1": return None, None - if ls_scheme == "condition2": + if param_scheme == "condition2": # The mirror of "condition1": estimate from condition 2's cells and give # the value to condition 1. Resolved eagerly rather than by inheritance # because model1 is fitted first; `_auto_ls` reproduces what a model @@ -48,18 +53,18 @@ def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor): # with the labels exchanged. Pinned by # test_auto_ls_is_bit_identical_to_the_models_internal_estimate. shared = _auto_ls(X_condition2, ls_factor) - logger.info(f"ls_scheme='condition2': shared length scale {shared:.4f}") + logger.info(f"param_scheme='condition2': shared length scale {shared:.4f}") return shared, shared - if ls_scheme == "separate": + if param_scheme == "separate": return None, None - if ls_scheme == "pooled": + if param_scheme == "pooled": shared = _auto_ls( np.vstack([np.asarray(X_condition1), np.asarray(X_condition2)]), ls_factor, ) - logger.info(f"ls_scheme='pooled': shared length scale {shared:.4f}") + logger.info(f"param_scheme='pooled': shared length scale {shared:.4f}") return shared, shared - if ls_scheme == "symmetric": + if param_scheme == "symmetric": n1 = np.asarray(X_condition1).shape[0] n2 = np.asarray(X_condition2).shape[0] ls1 = _auto_ls(X_condition1, ls_factor) @@ -72,12 +77,12 @@ def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor): # the union holds more cells. shared = float(np.exp((n1 * np.log(ls1) + n2 * np.log(ls2)) / (n1 + n2))) logger.info( - f"ls_scheme='symmetric': condition length scales {ls1:.4f} / " + f"param_scheme='symmetric': condition length scales {ls1:.4f} / " f"{ls2:.4f} -> shared {shared:.4f}" ) return shared, shared raise ValueError( - f"Unknown ls_scheme {ls_scheme!r}. Expected one of {LS_SCHEMES}." + f"Unknown param_scheme {param_scheme!r}. Expected one of {PARAM_SCHEMES}." ) @@ -376,7 +381,7 @@ def fit( sigma: float = 1.0, ls: Optional[float] = None, ls_factor: float = 10.0, - ls_scheme: str = "condition1", + param_scheme: Optional[str] = None, landmarks: Optional[np.ndarray] = None, sample_estimator_ls: Optional[float] = None, condition1_sample_indices: Optional[np.ndarray] = None, @@ -404,15 +409,19 @@ def fit( Noise level for function estimator, by default 1.0. ls : float, optional Length scale for the GP kernel, shared by both conditions. If None it is - estimated from the data according to ``ls_scheme``, by default None. + estimated from the data according to ``param_scheme``, by default None. ls_factor : float, optional Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is None. - ls_scheme : str, optional - How the automatic length scale is derived when ``ls`` is None. Both + param_scheme : str, optional + Where the hyperparameters estimated from cells come from. The same + field drives :meth:`DifferentialAbundance.fit`; here the only such + hyperparameter is the length scale, so the scheme decides which cells + the shared ``ls`` is estimated from when ``ls`` is None. Both conditions are normally smoothed at the *same* scale so that their fitted surfaces are comparable; the schemes differ in which cells the - shared value is estimated from. + shared value is estimated from. ``None`` (default) means + ``"condition1"``. * ``"condition1"`` (default) — estimate from condition 1's cells and reuse the value for condition 2. **Not symmetric**: the length @@ -422,12 +431,13 @@ def fit( * ``"condition2"`` — the mirror of the default: estimate from condition 2's cells and reuse the value for condition 1. Equally asymmetric, and deliberately so. Its use is **diagnostic**: - ``de(X, Y, ls_scheme="condition1")`` and - ``de(Y, X, ls_scheme="condition2")`` are the same computation with - the labels exchanged, so running both orientations of a contrast - under the two schemes exposes the default's swap-dependence from a - single call site, without having to swap ``condition1`` and - ``condition2`` at the call site and re-derive which sign is which. + ``de(X, Y, param_scheme="condition1")`` and + ``de(Y, X, param_scheme="condition2")`` are the same computation + with the labels exchanged, so running both orientations of a + contrast under the two schemes exposes the default's + swap-dependence from a single call site, without having to swap + ``condition1`` and ``condition2`` at the call site and re-derive + which sign is which. .. warning:: @@ -468,7 +478,7 @@ def fit( when ``"condition2"`` is used with automatic landmarks. The order dependence is **pre-existing and independent of** - ``ls_scheme`` — no scheme removes it, and it affects any + ``param_scheme`` — no scheme removes it, and it affects any comparison of two orientations, not only this one. * ``"symmetric"`` — estimate a length scale from each condition separately and share their size-weighted geometric mean. Invariant @@ -509,6 +519,8 @@ def fit( The fitted instance. """ + param_scheme = resolve_param_scheme(param_scheme, DEFAULT_PARAM_SCHEME) + # Check if sample indices are provided have_sample_indices = ( condition1_sample_indices is not None @@ -567,13 +579,13 @@ def fit( # scheme's whole purpose is the comparison, so anyone who reaches here is # relying on it; say so rather than leaving the boundary in the docs. if ( - ls_scheme == "condition2" + param_scheme == "condition2" and not landmarks_were_provided and self.n_landmarks is not None and self.n_landmarks > 0 ): logger.warning( - "ls_scheme='condition2' with automatic landmarks: the " + "param_scheme='condition2' with automatic landmarks: the " "condition1/condition2 equivalence is EXACT only when both " "orientations use the same landmarks, and automatic landmarks " "are order-dependent (they are computed from the two " @@ -600,17 +612,17 @@ def fit( # length scale smuggled in through function_kwargs. A supplied # `cov_func`/`cov_func_curry` must NOT disable it -- the two conditions # still need a shared scale, and gating on the kernel would drop the - # condition-1 inheritance for callers who never touched `ls_scheme`, + # condition-1 inheritance for callers who never touched `param_scheme`, # silently giving them "separate" behaviour. ls_auto = ls is None and "ls" not in function_kwargs if ls_auto: - ls_model1, ls_model2 = _resolve_ls_scheme( - ls_scheme, X_condition1, X_condition2, ls_factor + ls_model1, ls_model2 = _resolve_scheme_ls( + param_scheme, X_condition1, X_condition2, ls_factor ) else: - if ls_scheme != "condition1": + if param_scheme != "condition1": logger.info( - f"ls_scheme={ls_scheme!r} ignored: the length scale is " + f"param_scheme={param_scheme!r} ignored: the length scale is " "already fixed by an explicit ls." ) ls_model1 = ls_model2 = ls @@ -638,7 +650,7 @@ def fit( ls_for_model2 = ls_model2 if ( ls_auto - and ls_scheme == "condition1" + and param_scheme == "condition1" and self.model1.ls is not None ): ls_for_model2 = self.model1.ls diff --git a/kompot/settings.py b/kompot/settings.py index c6f84e0..19ff328 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -28,6 +28,28 @@ validate_bool, ) +#: Recognised values for ``GPSettings.param_scheme``. ``None`` is also +#: accepted and means "this entry point's default". +PARAM_SCHEMES = ("condition1", "condition2", "symmetric", "pooled", "separate") + + +def validate_param_scheme(value, optional: bool = False): + """Return ``value`` if it is a recognised ``param_scheme``, else raise.""" + if value is None: + if optional: + return None + raise ValueError(f"'param_scheme' must be one of {PARAM_SCHEMES} (got None).") + if value not in PARAM_SCHEMES: + raise ValueError( + f"'param_scheme' must be one of {PARAM_SCHEMES} (got {value!r})." + ) + return value + + +def resolve_param_scheme(value, default: str) -> str: + """Map ``None`` onto an entry point's default and validate the result.""" + return validate_param_scheme(default if value is None else value) + @dataclass class GPSettings: @@ -38,27 +60,44 @@ class GPSettings: sigma : float Noise level for the GP. ls : float, optional - Length scale, shared by both conditions. If *None*, estimated - automatically using ``ls_factor`` and ``ls_scheme``. + Length scale of the GP kernel, shared by both conditions. If *None* + it is estimated from the cells, from where ``param_scheme`` says. ls_factor : float Multiplier applied to the automatically inferred length scale. - ls_scheme : str - Which cells the automatic shared length scale is estimated from. - ``"condition1"`` (default) uses condition 1's cells only, which makes - the result depend on which condition is passed first; ``"condition2"`` - is its mirror and exists to make that dependence visible in a single - run, since ``de(X, Y, ls_scheme="condition1")`` and - ``de(Y, X, ls_scheme="condition2")`` are the same computation with the - labels exchanged **provided both runs use the same landmarks** — see - the warning under :meth:`kompot.differential.DifferentialExpression.fit`, - because the defaults here (``n_landmarks=5000``, ``landmarks=None``) do - *not* satisfy that and the equivalence then holds only approximately; - ``"symmetric"`` shares the size-weighted geometric - mean of the two per-condition estimates; ``"pooled"`` estimates from - both conditions' cells taken together; ``"separate"`` gives each - condition its own. See - :meth:`kompot.differential.DifferentialExpression.fit`. Ignored when - ``ls`` is given explicitly. + param_scheme : str, optional + Where the hyperparameters a model estimates from cells come from. + Read by both :func:`kompot.de` and :func:`kompot.da`; the model + decides which parameters it covers — ``ls`` in differential + expression, ``d``, ``mu`` and ``ls`` in differential abundance. A + value given explicitly (``ls`` here, or ``d`` / ``mu`` passed to + ``da()``) pins that parameter and takes precedence. + + ``None`` (default) keeps each entry point's own default: + ``"condition1"`` for ``de()`` and ``"separate"`` for ``da()``. + + * ``"condition1"`` — estimate from condition 1's cells and reuse for + condition 2. Makes the result depend on which condition is passed + first. + * ``"condition2"`` — the mirror: estimate from condition 2's cells and + reuse for condition 1. Equally asymmetric, by design; it exists so + that ``de(X, Y, "condition1")`` and ``de(Y, X, "condition2")`` (and + likewise ``da``) are the same computation with the labels exchanged + **provided both runs use the same landmarks**. The defaults here + (``n_landmarks=5000``, ``landmarks=None``) do *not* satisfy that and + the equivalence then holds only approximately — see + :meth:`kompot.differential.DifferentialExpression.fit`. + * ``"symmetric"`` — estimate from each condition separately and share + the size-weighted combination (geometric mean for ``ls``). + Invariant under swapping the conditions. + * ``"pooled"`` — estimate from both conditions' cells taken together. + Swap-invariant up to the row order of the stacked union, and the + union is denser than either condition, so the length scale shrinks + with cell count alone. + * ``"separate"`` — each condition estimates its own; nothing shared. + + Which to pick, and the measurements behind each: + :meth:`kompot.differential.DifferentialExpression.fit` and + :meth:`kompot.differential.DifferentialAbundance.fit`. n_landmarks : int, optional Number of landmarks for the Nystrom approximation. landmarks : np.ndarray, optional @@ -78,7 +117,7 @@ class GPSettings: sigma: float = 1.0 ls: Optional[float] = None ls_factor: float = 10.0 - ls_scheme: str = "condition1" + param_scheme: Optional[str] = None n_landmarks: Optional[int] = 5000 landmarks: Optional[np.ndarray] = None use_empirical_variance: bool = False @@ -91,13 +130,7 @@ def __post_init__(self): validate_positive_float(self.sigma, "sigma") validate_positive_float(self.ls, "ls", optional=True) validate_positive_float(self.ls_factor, "ls_factor") - from .differential.differential_expression import LS_SCHEMES - - if self.ls_scheme not in LS_SCHEMES: - raise ValueError( - f"'ls_scheme' must be one of {LS_SCHEMES} (got " - f"{self.ls_scheme!r})." - ) + validate_param_scheme(self.param_scheme, optional=True) validate_positive_int(self.n_landmarks, "n_landmarks", optional=True) validate_bool(self.use_empirical_variance, "use_empirical_variance") validate_positive_int(self.batch_size, "batch_size", optional=True) diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index 9844109..9e28222 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -1,9 +1,8 @@ -"""Swap equivalence for differential ABUNDANCE. +"""Swap equivalence for differential ABUNDANCE at its default ``param_scheme``. -`ls_scheme` and `tests/test_ls_scheme.py` are about differential EXPRESSION, where -the shared length scale was estimated from condition 1's cells alone and -`"condition2"` exists to expose that. Abundance is a separate code path and -inherits none of it, so the property is established here rather than assumed. +``tests/test_param_scheme.py`` covers what each ``param_scheme`` shares, +including the ``"condition1"`` / ``"condition2"`` mirror. This file establishes +the property for the DEFAULT, ``"separate"``, rather than assuming it. It holds for a structural reason worth stating: `DifferentialAbundance.fit` builds BOTH density estimators from the same ``estimator_defaults`` dict -- @@ -11,19 +10,19 @@ density_estimator_condition1 = mellon.DensityEstimator(**estimator_defaults) density_estimator_condition2 = mellon.DensityEstimator(**estimator_defaults) --- with no ``ls_for_model2 = self.model1.ls`` analogue anywhere. There is no -condition-1 inheritance to be asymmetric about, so DA never had the defect -`ls_scheme` was added for, and no DA counterpart of `"condition2"` is warranted. +-- and under ``"separate"`` that dict carries no ``d``, ``mu`` or ``ls``, so each +estimator derives its own from its own cells. There is no condition-1 +inheritance to be asymmetric about. Two configurations DO break the exactness, and both are covered below: automatic -landmarks, and ``sync_parameters=True``. Both trace to the same ROOT as the DE +landmarks, and ``param_scheme="pooled"`` (formerly ``sync_parameters=True``). Both trace to the same ROOT as the DE case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent -- but they act through DIFFERENT proximate mechanisms, and **a remedy applies at the proximate level, not the root**: ``compute_landmarks`` (k-means init) for the first, ``compute_nn_distances`` (which sets ``mu`` and ``ls``) for the second. So sharing a landmark array fixes the first and cannot reach the second -- landmarks are not an input to ``compute_nn_distances``, and sharing them leaves -the ``sync_parameters`` half at 4.784e-03 against a 4.720e-03 no-landmark +the ``pooled`` half at 4.784e-03 against a 4.720e-03 no-landmark reference. The remedy for that half is to pass ``mu`` and ``ls`` explicitly -- **and ``d`` @@ -32,11 +31,11 @@ ``parameters.py``) and so becomes row-order dependent itself. Measured: ``mu``+``ls`` alone is exact at 200+200 and NOT exact at 400+400 (6.085e-04, with ``d`` reading 1.4972 vs 1.4983); all three explicit is exact at 400+400. -With ``sync_parameters=True``, passing all three is always correct. +Under ``"pooled"``, passing all three is always correct. **THE TWO REMEDIES COMPOSE AND NEITHER SUBSTITUTES FOR THE OTHER.** Passing ``d``, ``mu`` and ``ls`` fixes the PARAMETER half and leaves the LANDMARK half -untouched: measured at 400+400, sync with automatic landmarks and all three +untouched: measured at 400+400, "pooled" with automatic landmarks and all three parameters explicit still differs between orientations on the uncertainty (8.4e-01), while the same run with one shared landmark array agrees exactly. So if landmarks are automatic, pass one ``landmarks`` array as well. @@ -67,16 +66,16 @@ def _data(n=400, seed=0): return X1, X2 -def _swap_pair(n_landmarks=None, sync_parameters=False, landmarks=None): +def _swap_pair(n_landmarks=None, param_scheme=None, landmarks=None): """Fit both orientations and predict both on the same points.""" X1, X2 = _data() X_new = np.vstack([X1, X2]) fwd = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) - fwd.fit(X1, X2, sync_parameters=sync_parameters, landmarks=landmarks) + fwd.fit(X1, X2, param_scheme=param_scheme, landmarks=landmarks) rev = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) - rev.fit(X2, X1, sync_parameters=sync_parameters, landmarks=landmarks) + rev.fit(X2, X1, param_scheme=param_scheme, landmarks=landmarks) return ( fwd.predict(X_new, progress=False), @@ -194,9 +193,9 @@ def test_the_significant_set_is_the_same_cells_either_way(default_pair): id="automatic-landmarks", ), pytest.param( - dict(sync_parameters=True), + dict(param_scheme="pooled"), "d / mu / ls synchronised from the stacked union", - id="sync-parameters", + id="pooled", ), ], ) @@ -204,7 +203,7 @@ def test_the_swap_is_only_approximate_when_the_union_order_matters(kwargs, why): """Characterisation: two configurations break the exactness, for one reason. ``np.vstack([X_condition1, X_condition2])`` is row-order dependent. It feeds - ``compute_landmarks`` always, and under ``sync_parameters=True`` it also feeds + ``compute_landmarks`` always, and under ``param_scheme="pooled"`` it also feeds ``compute_d_factal`` / ``compute_mu`` / ``compute_ls``. Either way the two orientations are then computed from differently-ordered unions and agreement degrades -- silently, which is why it is pinned here. @@ -231,11 +230,11 @@ def test_the_swap_is_only_approximate_when_the_union_order_matters(kwargs, why): def test_shared_landmarks_restore_the_exact_swap(): """Share the landmarks, and the automatic-landmark half is exact again. - Scope: this covers the LANDMARK half only -- ``sync_parameters`` is left at - its default ``False`` here. The same remedy does not reach the - ``sync_parameters`` half, and no test asserts that it does. Nor does the + Scope: this covers the LANDMARK half only -- ``param_scheme`` is left at + its default ``"separate"`` here. The same remedy does not reach the + ``pooled`` half, and no test asserts that it does. Nor does the parameter remedy reach THIS half: the two compose, and a run with - ``sync_parameters=True`` needs both a shared landmark array and explicit + ``param_scheme="pooled"`` needs both a shared landmark array and explicit ``d``/``mu``/``ls`` before the two orientations agree. """ from mellon.parameters import compute_landmarks diff --git a/tests/test_differential.py b/tests/test_differential.py index 945a412..a7a7630 100644 --- a/tests/test_differential.py +++ b/tests/test_differential.py @@ -49,8 +49,8 @@ def test_differential_abundance_fit(): # No longer checking for class attributes - they're not updated in the new version -def test_differential_abundance_sync_parameters(): - """Test synchronizing parameters between conditions with DifferentialAbundance.""" +def test_differential_abundance_pooled_parameters(): + """Test sharing pooled parameters between conditions with DifferentialAbundance.""" # Generate data with different distributions np.random.seed(42) # Set seed for reproducibility X_condition1 = np.random.randn(100, 5) @@ -58,11 +58,11 @@ def test_differential_abundance_sync_parameters(): # First fit without parameter synchronization diff_abundance_nosync = DifferentialAbundance() - diff_abundance_nosync.fit(X_condition1, X_condition2, sync_parameters=False) + diff_abundance_nosync.fit(X_condition1, X_condition2, param_scheme="separate") # Then fit with parameter synchronization diff_abundance_sync = DifferentialAbundance() - diff_abundance_sync.fit(X_condition1, X_condition2, sync_parameters=True) + diff_abundance_sync.fit(X_condition1, X_condition2, param_scheme="pooled") # Both should have valid predictors assert diff_abundance_nosync.density_predictor1 is not None @@ -114,7 +114,7 @@ def test_differential_abundance_sync_parameters(): diff_abundance_custom.fit( X_condition1, X_condition2, - sync_parameters=True, # Should sync but these values take precedence + param_scheme="pooled", # Should share but these values take precedence d=custom_d, mu=custom_mu, ls=custom_ls, diff --git a/tests/test_differential_extended.py b/tests/test_differential_extended.py index a20cc2e..0859910 100644 --- a/tests/test_differential_extended.py +++ b/tests/test_differential_extended.py @@ -120,8 +120,8 @@ def test_differential_abundance_fit_with_landmarks(self): assert hasattr(da, "computed_landmarks") np.testing.assert_array_equal(da.computed_landmarks, landmarks) - def test_differential_abundance_fit_sync_parameters(self): - """Test DifferentialAbundance fit with parameter synchronization.""" + def test_differential_abundance_fit_pooled_parameters(self): + """Test DifferentialAbundance fit with pooled shared parameters.""" try: from kompot.differential import DifferentialAbundance except ImportError as e: @@ -148,7 +148,7 @@ def test_differential_abundance_fit_sync_parameters(self): mock_estimator.predict = mock_predictor mock_mellon.DensityEstimator.return_value = mock_estimator - da.fit(X1, X2, sync_parameters=True, n_landmarks=10) + da.fit(X1, X2, param_scheme="pooled", n_landmarks=10) # Should call parameter computation functions mock_mellon.parameters.compute_d_factal.assert_called_once() diff --git a/tests/test_ls_scheme.py b/tests/test_param_scheme.py similarity index 57% rename from tests/test_ls_scheme.py rename to tests/test_param_scheme.py index b0e5745..51b4354 100644 --- a/tests/test_ls_scheme.py +++ b/tests/test_param_scheme.py @@ -1,11 +1,15 @@ -"""Length-scale determination schemes for the two-condition expression fit. +"""``param_scheme``: where the hyperparameters of a two-condition fit come from. -The default (``ls_scheme="condition1"``) estimates the shared length scale from +One ``GPSettings`` field drives both entry points. In differential EXPRESSION +it covers the length scale; in differential ABUNDANCE it covers ``d``, ``mu`` +and ``ls`` (the second half of this file). + +The expression default (``param_scheme=None`` -> ``"condition1"``) estimates the shared length scale from condition 1's cells only, which makes the contrast depend on which condition is passed first. These tests pin the default's behaviour and check that each alternative does what it says -- including ``"condition2"``, the default's mirror, for which the statement is an exact equivalence: -``de(X, Y, ls_scheme="condition1")`` and ``de(Y, X, ls_scheme="condition2")`` +``de(X, Y, param_scheme="condition1")`` and ``de(Y, X, param_scheme="condition2")`` are the same computation with the labels exchanged. """ @@ -16,12 +20,11 @@ import pytest from kompot.differential.differential_expression import ( - LS_SCHEMES, DifferentialExpression, _auto_ls, - _resolve_ls_scheme, + _resolve_scheme_ls, ) -from kompot.settings import GPSettings +from kompot.settings import PARAM_SCHEMES, GPSettings def _data(n1=120, n2=40, n_genes=6, seed=0): @@ -38,13 +41,13 @@ def _data(n1=120, n2=40, n_genes=6, seed=0): def test_condition1_and_separate_defer_to_the_models(): X1, _, X2, _ = _data() - assert _resolve_ls_scheme("condition1", X1, X2, 10.0) == (None, None) - assert _resolve_ls_scheme("separate", X1, X2, 10.0) == (None, None) + assert _resolve_scheme_ls("condition1", X1, X2, 10.0) == (None, None) + assert _resolve_scheme_ls("separate", X1, X2, 10.0) == (None, None) def test_pooled_uses_the_union_and_is_smaller_than_either_condition(): X1, _, X2, _ = _data() - p1, p2 = _resolve_ls_scheme("pooled", X1, X2, 10.0) + p1, p2 = _resolve_scheme_ls("pooled", X1, X2, 10.0) assert p1 == p2 assert p1 == pytest.approx(_auto_ls(np.vstack([X1, X2]), 10.0)) # The union is denser than either part, so nearest-neighbour distances -- @@ -57,7 +60,7 @@ def test_pooled_uses_the_union_and_is_smaller_than_either_condition(): def test_symmetric_is_the_size_weighted_geometric_mean(): X1, _, X2, _ = _data() - got, got2 = _resolve_ls_scheme("symmetric", X1, X2, 10.0) + got, got2 = _resolve_scheme_ls("symmetric", X1, X2, 10.0) assert got == got2 a, b = _auto_ls(X1, 10.0), _auto_ls(X2, 10.0) n1, n2 = X1.shape[0], X2.shape[0] @@ -67,15 +70,15 @@ def test_symmetric_is_the_size_weighted_geometric_mean(): def test_symmetric_is_invariant_under_swapping_the_conditions(): X1, _, X2, _ = _data() - fwd, _ = _resolve_ls_scheme("symmetric", X1, X2, 10.0) - rev, _ = _resolve_ls_scheme("symmetric", X2, X1, 10.0) + fwd, _ = _resolve_scheme_ls("symmetric", X1, X2, 10.0) + rev, _ = _resolve_scheme_ls("symmetric", X2, X1, 10.0) assert fwd == pytest.approx(rev, rel=1e-9) def test_pooled_is_invariant_under_swapping_the_conditions(): X1, _, X2, _ = _data() - fwd, _ = _resolve_ls_scheme("pooled", X1, X2, 10.0) - rev, _ = _resolve_ls_scheme("pooled", X2, X1, 10.0) + fwd, _ = _resolve_scheme_ls("pooled", X1, X2, 10.0) + rev, _ = _resolve_scheme_ls("pooled", X2, X1, 10.0) # The union is the same point set either way; the only difference is row # order, which the approximate nearest-neighbour index is mildly # sensitive to. @@ -84,19 +87,23 @@ def test_pooled_is_invariant_under_swapping_the_conditions(): def test_unknown_scheme_is_rejected(): X1, _, X2, _ = _data() - with pytest.raises(ValueError, match="Unknown ls_scheme"): - _resolve_ls_scheme("nope", X1, X2, 10.0) + with pytest.raises(ValueError, match="Unknown param_scheme"): + _resolve_scheme_ls("nope", X1, X2, 10.0) + with pytest.raises(ValueError, match="param_scheme"): + DifferentialExpression(n_landmarks=0).fit(*_data(), param_scheme="nope") def test_gp_settings_rejects_unknown_scheme(): - with pytest.raises(ValueError, match="ls_scheme"): - GPSettings(ls_scheme="nope") - for scheme in LS_SCHEMES: - GPSettings(ls_scheme=scheme) + with pytest.raises(ValueError, match="param_scheme"): + GPSettings(param_scheme="nope") + for scheme in PARAM_SCHEMES: + GPSettings(param_scheme=scheme) -def test_gp_settings_default_preserves_historical_behaviour(): - assert GPSettings().ls_scheme == "condition1" +def test_gp_settings_default_is_each_entry_points_own(): + """``None`` defers to the entry point; for expression that is condition1.""" + assert GPSettings().param_scheme is None + assert _fit(None) == _fit("condition1") # -------------------------------------------------------------------------- @@ -109,7 +116,7 @@ def _fit(scheme, swap=False, **kw): if swap: X1, y1, X2, y2 = X2, y2, X1, y1 de = DifferentialExpression(n_landmarks=0) - de.fit(X1, y1, X2, y2, ls_scheme=scheme, **kw) + de.fit(X1, y1, X2, y2, param_scheme=scheme, **kw) return float(de.model1.ls), float(de.model2.ls) @@ -144,7 +151,7 @@ def test_shared_schemes_are_swap_invariant_end_to_end(scheme): assert fwd[0] == pytest.approx(rev[0], rel=1e-2) -@pytest.mark.parametrize("scheme", LS_SCHEMES) +@pytest.mark.parametrize("scheme", PARAM_SCHEMES) def test_explicit_ls_overrides_every_scheme(scheme): ls1, ls2 = _fit(scheme, ls=0.5) assert ls1 == pytest.approx(0.5) @@ -161,7 +168,7 @@ def _fit_with_curry(scheme): X1, y1, X2, y2 = _data() de = DifferentialExpression(n_landmarks=0) - de.fit(X1, y1, X2, y2, ls_scheme=scheme, cov_func_curry=Matern52) + de.fit(X1, y1, X2, y2, param_scheme=scheme, cov_func_curry=Matern52) return float(de.model1.ls), float(de.model2.ls) @@ -169,7 +176,7 @@ def test_default_still_shares_when_a_cov_func_curry_is_supplied(): """Regression: `cov_func_curry` must not turn the default into `separate`. `function_kwargs` is public and the CLI forwards into it, so a caller who - never touched `ls_scheme` can reach this path. Before the fix, supplying a + never touched `param_scheme` can reach this path. Before the fix, supplying a curry dropped the condition-1 inheritance and let condition 2 estimate its own -- i.e. `separate` behaviour, silently, at the default. """ @@ -213,13 +220,13 @@ def test_auto_ls_is_bit_identical_to_the_models_internal_estimate(): """ X1, y1, X2, y2 = _data() de = DifferentialExpression(n_landmarks=0) - de.fit(X1, y1, X2, y2, ls_scheme="condition1") + de.fit(X1, y1, X2, y2, param_scheme="condition1") assert float(de.model1.ls) == _auto_ls(X1, 10.0) def test_condition2_shares_condition2s_length_scale(): X1, _, X2, _ = _data() - a, b = _resolve_ls_scheme("condition2", X1, X2, 10.0) + a, b = _resolve_scheme_ls("condition2", X1, X2, 10.0) assert a == b assert a == pytest.approx(_auto_ls(X2, 10.0)) @@ -246,7 +253,7 @@ def test_condition2_is_not_invariant_under_swapping_the_conditions(): def test_condition2_resolves_to_what_condition1_resolves_to_on_the_swapped_pair(): """The resolver-level statement of the mirror.""" X1, _, X2, _ = _data() - c2, _ = _resolve_ls_scheme("condition2", X1, X2, 10.0) + c2, _ = _resolve_scheme_ls("condition2", X1, X2, 10.0) # `condition1` defers to the model, so its value is `_auto_ls` of whatever # is passed first -- which after the swap is X2. assert c2 == _auto_ls(X2, 10.0) @@ -279,10 +286,10 @@ def _mirror_pair(n_landmarks=0, landmarks=None, n_genes=6): X_new = np.vstack([X1, X2]) fwd = DifferentialExpression(n_landmarks=n_landmarks, random_state=0) - fwd.fit(X1, y1, X2, y2, ls_scheme="condition1", landmarks=landmarks) + fwd.fit(X1, y1, X2, y2, param_scheme="condition1", landmarks=landmarks) mir = DifferentialExpression(n_landmarks=n_landmarks, random_state=0) - mir.fit(X2, y2, X1, y1, ls_scheme="condition2", landmarks=landmarks) + mir.fit(X2, y2, X1, y1, param_scheme="condition2", landmarks=landmarks) return ( fwd, @@ -372,7 +379,7 @@ def test_the_mirror_holds_with_landmarks_when_the_landmarks_are_shared(): `fit()` derives automatic landmarks from `np.vstack([X_condition1, X_condition2])`, whose ROW ORDER differs between the two orientations, so the automatic path evaluates the two runs at different points. That is a - pre-existing, scheme-independent asymmetry; no `ls_scheme` removes it. Pass + pre-existing, scheme-independent asymmetry; no `param_scheme` removes it. Pass one landmark array to both runs and the equivalence is exact again. """ from mellon.parameters import compute_landmarks @@ -397,12 +404,12 @@ def test_the_mirror_holds_with_landmarks_when_the_landmarks_are_shared(): # -------------------------------------------------------------------------- -# every reader of `ls_scheme`, not just the writer +# every reader of `param_scheme`, not just the writer # -------------------------------------------------------------------------- -def test_ls_scheme_is_readable_back_out_of_a_stored_params_dict(): - """`ls_scheme` is in the run-parameter match list; it must also be READABLE. +def test_param_scheme_is_readable_back_out_of_a_stored_params_dict(): + """`param_scheme` is in the run-parameter match list; it must also be READABLE. `build_params_dict` stores it nested under `params["gp"]`, and `_check_overwrites` compares it through `params_get`. `params_get` @@ -412,16 +419,16 @@ def test_ls_scheme_is_readable_back_out_of_a_stored_params_dict(): """ from kompot.anndata.utils.params import build_params_dict, params_get - for scheme in LS_SCHEMES: + for scheme in PARAM_SCHEMES: params = build_params_dict( {"groupby": "g", "condition1": "A", "condition2": "B"}, - gp=GPSettings(ls_scheme=scheme), + gp=GPSettings(param_scheme=scheme), ) - assert params["gp"]["ls_scheme"] == scheme - assert params_get(params, "ls_scheme") == scheme + assert params["gp"]["param_scheme"] == scheme + assert params_get(params, "param_scheme") == scheme -def test_the_de_cli_routes_ls_scheme_into_gp_settings(): +def test_the_de_cli_routes_param_scheme_into_gp_settings(): """A flat config key the CLI does not recognise is not rejected. It falls through to `**function_kwargs` and is handed to mellon, so an @@ -432,7 +439,7 @@ def test_the_de_cli_routes_ls_scheme_into_gp_settings(): from kompot.cli.de import GP_CONFIG_KEYS - assert "ls_scheme" in GP_CONFIG_KEYS + assert "param_scheme" in GP_CONFIG_KEYS fields = {f.name for f in dataclasses.fields(GPSettings)} assert GP_CONFIG_KEYS <= fields, GP_CONFIG_KEYS - fields assert fields - GP_CONFIG_KEYS == {"landmarks"} @@ -496,7 +503,7 @@ def test_condition2_warns_when_the_landmarks_are_not_shared(): with _captured_kompot_warnings() as messages: DifferentialExpression(n_landmarks=20, random_state=0).fit( - X1, y1, X2, y2, ls_scheme="condition2" + X1, y1, X2, y2, param_scheme="condition2" ) assert any("same landmarks" in m for m in messages), messages @@ -524,7 +531,7 @@ def test_condition2_does_not_warn_when_the_equivalence_is_exact(kwargs): with _captured_kompot_warnings() as messages: DifferentialExpression(random_state=0, **kwargs).fit( - X1, y1, X2, y2, ls_scheme="condition2", landmarks=landmarks + X1, y1, X2, y2, param_scheme="condition2", landmarks=landmarks ) assert not [m for m in messages if "same landmarks" in m] @@ -535,6 +542,287 @@ def test_the_other_schemes_do_not_warn_about_landmarks(): for scheme in ("condition1", "symmetric", "pooled", "separate"): with _captured_kompot_warnings() as messages: DifferentialExpression(n_landmarks=20, random_state=0).fit( - X1, y1, X2, y2, ls_scheme=scheme + X1, y1, X2, y2, param_scheme=scheme ) assert not [m for m in messages if "same landmarks" in m], scheme + + +# ========================================================================== +# differential ABUNDANCE: the same field covers d, mu and ls +# ========================================================================== + +import warnings # noqa: E402 + +import anndata # noqa: E402 + +import kompot # noqa: E402 +from kompot.differential.differential_abundance import ( # noqa: E402 + DifferentialAbundance, + _resolve_scheme_density_params, +) +from kompot.settings import OutputSettings # noqa: E402 + +# Separated enough that the fold change is far from zero, so every equality +# below compares arrays that COULD differ. +_DA_SEP, _DA_SCALE = 2.5, 0.45 + + +def _da_data(n1=90, n2=60, seed=0): + rng = np.random.default_rng(seed) + X1 = rng.normal(loc=[-_DA_SEP, 0.0], scale=_DA_SCALE, size=(n1, 2)) + X2 = rng.normal(loc=[+_DA_SEP, 0.0], scale=_DA_SCALE, size=(n2, 2)) + return X1, X2 + + +def _da_fit(X1, X2, **fit_kwargs): + """Fit with no landmarks, so the only orientation effect is the scheme's.""" + model = DifferentialAbundance(n_landmarks=None, random_state=0) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", DeprecationWarning) + model.fit(X1, X2, **fit_kwargs) + return model + + +def _da_predict(model, X_new): + return model.predict(X_new, progress=False) + + +def _shared(model): + p1, p2 = model.density_predictor1, model.density_predictor2 + return ( + (float(p1.d), float(p1.mu), float(p1.cov_func.ls)), + (float(p2.d), float(p2.mu), float(p2.cov_func.ls)), + ) + + +# -- defaults and the deprecated alias -------------------------------------- + + +def test_da_default_is_separate(): + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + a = _da_predict(_da_fit(X1, X2), Xn) + b = _da_predict(_da_fit(X1, X2, param_scheme="separate"), Xn) + for key in ("log_density_condition1", "log_density_condition2", "log_fold_change"): + np.testing.assert_array_equal(a[key], b[key]) + + +@pytest.mark.parametrize("sync, scheme", [(True, "pooled"), (False, "separate")]) +def test_sync_parameters_is_an_exact_alias(sync, scheme): + """The shipped spelling and the new one are the same computation, to the bit.""" + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + model = DifferentialAbundance(n_landmarks=None, random_state=0) + with pytest.warns(DeprecationWarning, match="param_scheme"): + model.fit(X1, X2, sync_parameters=sync) + old = _da_predict(model, Xn) + new = _da_predict(_da_fit(X1, X2, param_scheme=scheme), Xn) + for key in ( + "log_density_condition1", + "log_density_condition2", + "log_fold_change", + "log_fold_change_uncertainty", + ): + np.testing.assert_array_equal(old[key], new[key]) + + +def test_pooled_and_separate_differ(): + """Positive control for the alias test: the two targets are distinguishable.""" + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + a = _da_predict(_da_fit(X1, X2, param_scheme="pooled"), Xn) + b = _da_predict(_da_fit(X1, X2, param_scheme="separate"), Xn) + assert np.abs(a["log_fold_change"] - b["log_fold_change"]).max() > 1e-2 + + +def test_sync_parameters_and_param_scheme_together_raise(): + X1, X2 = _da_data() + with pytest.raises(ValueError, match="not both"): + DifferentialAbundance().fit(X1, X2, sync_parameters=True, param_scheme="pooled") + + +def test_da_rejects_unknown_scheme(): + X1, X2 = _da_data() + with pytest.raises(ValueError, match="param_scheme"): + DifferentialAbundance().fit(X1, X2, param_scheme="nope") + + +# -- what each scheme shares ------------------------------------------------ + + +def test_da_separate_does_not_share(): + X1, X2 = _da_data() + s1, s2 = _shared(_da_fit(X1, X2, param_scheme="separate")) + assert s1[0] != s2[0] and s1[1] != s2[1] and s1[2] != s2[2] + + +@pytest.mark.parametrize("scheme", ["pooled", "condition1", "condition2", "symmetric"]) +def test_da_sharing_schemes_share_d_mu_and_ls(scheme): + X1, X2 = _da_data() + s1, s2 = _shared(_da_fit(X1, X2, param_scheme=scheme)) + assert s1 == s2 + + +@pytest.mark.parametrize("scheme, side", [("condition1", 0), ("condition2", 1)]) +def test_da_one_sided_scheme_shares_that_conditions_own_estimate(scheme, side): + """The shared values are what that condition's estimator reaches on its own. + + Load-bearing for the mirror below: the scheme resolves eagerly, and it is + the same computation as ``"separate"`` for the condition it names only if + the eager estimate equals the estimator's internal one to the last bit. + """ + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + sep = _da_fit(X1, X2, param_scheme="separate") + one = _da_fit(X1, X2, param_scheme=scheme) + assert _shared(one)[side] == _shared(sep)[side] + key = ("log_density_condition1", "log_density_condition2")[side] + np.testing.assert_array_equal(_da_predict(one, Xn)[key], _da_predict(sep, Xn)[key]) + + +@pytest.mark.parametrize("scheme", PARAM_SCHEMES) +def test_da_explicit_values_pin_every_scheme(scheme): + X1, X2 = _da_data() + s1, s2 = _shared(_da_fit(X1, X2, param_scheme=scheme, d=2.0, mu=-5.0, ls=3.0)) + assert s1 == s2 == (2.0, -5.0, 3.0) + + +# -- swap behaviour --------------------------------------------------------- + + +def test_da_mirror_is_a_relabelling_of_the_same_fit(): + """``da(X, Y, "condition1")`` and ``da(Y, X, "condition2")``, labels exchanged.""" + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + fwd = _da_predict(_da_fit(X1, X2, param_scheme="condition1"), Xn) + mir = _da_predict(_da_fit(X2, X1, param_scheme="condition2"), Xn) + assert np.abs(fwd["log_fold_change"]).max() > 1e-2 + np.testing.assert_array_equal(fwd["log_density_condition1"], mir["log_density_condition2"]) + np.testing.assert_array_equal(fwd["log_density_condition2"], mir["log_density_condition1"]) + np.testing.assert_array_equal(fwd["log_fold_change"], -mir["log_fold_change"]) + np.testing.assert_array_equal( + fwd["log_fold_change_uncertainty"], mir["log_fold_change_uncertainty"] + ) + + +def test_da_one_sided_schemes_are_not_swap_invariant(): + """Characterisation: ``condition1`` is asymmetric by design.""" + X1, X2 = _da_data() + fwd = _shared(_da_fit(X1, X2, param_scheme="condition1"))[0] + rev = _shared(_da_fit(X2, X1, param_scheme="condition1"))[0] + assert fwd != rev + + +@pytest.mark.parametrize("n1, n2", [(90, 60), (700, 300)]) +def test_da_symmetric_is_bit_identical_under_a_swap(n1, n2): + """Swap-invariant by construction, including above the 500-cell ``d`` subsample. + + Positive control: ``"pooled"``'s ``d`` on the same data IS order-dependent + above 500 combined cells, so the size is one where a union-based rule fails. + """ + X1, X2 = _da_data(n1, n2) + fwd = _resolve_scheme_density_params("symmetric", X1, X2, 10.0, 42) + rev = _resolve_scheme_density_params("symmetric", X2, X1, 10.0, 42) + assert fwd == rev + one = _resolve_scheme_density_params("condition1", X1, X2, 10.0, 42) + two = _resolve_scheme_density_params("condition2", X1, X2, 10.0, 42) + for name in ("d", "mu", "ls"): + assert fwd[name] != one[name] and fwd[name] != two[name], name + assert min(one[name], two[name]) <= fwd[name] <= max(one[name], two[name]), name + if n1 + n2 > 500: + from mellon.parameters import compute_d_factal + + assert compute_d_factal(np.vstack([X1, X2])) != compute_d_factal( + np.vstack([X2, X1]) + ) + + +def test_da_symmetric_swap_is_exact_end_to_end(): + X1, X2 = _da_data() + Xn = np.vstack([X1, X2]) + fwd = _da_predict(_da_fit(X1, X2, param_scheme="symmetric"), Xn) + rev = _da_predict(_da_fit(X2, X1, param_scheme="symmetric"), Xn) + assert np.abs(fwd["log_fold_change"]).max() > 1e-2 + np.testing.assert_array_equal(fwd["log_fold_change"], -rev["log_fold_change"]) + np.testing.assert_array_equal( + fwd["log_fold_change_uncertainty"], rev["log_fold_change_uncertainty"] + ) + + +# -- kompot.da(): GPSettings is honoured ------------------------------------ + + +def _adata(n1=90, n2=60): + X1, X2 = _da_data(n1, n2) + X = np.vstack([X1, X2]) + ad = anndata.AnnData(X=np.zeros((X.shape[0], 1), dtype=np.float32)) + ad.obs["cond"] = ["A"] * n1 + ["B"] * n2 + ad.obsm["X"] = X + return ad + + +def _da_lfc(gp=None, **kwargs): + result, _ = kompot.da( + _adata(), + "cond", + "A", + "B", + obsm_key="X", + gp=gp, + output=OutputSettings(copy=True, return_full_results=True, progress=False), + **kwargs, + ) + return result["table"]["lfc"].to_numpy() + + +_NO_LM = dict(n_landmarks=None, random_state=0) + + +def test_da_honours_gp_ls(): + """Regression: ``da()`` used to read six GPSettings fields and drop ``ls``. + + The control is the same length scale passed as a density keyword, which has + always reached mellon. If ``gp.ls`` were still ignored, the first result + would equal the default run, not the control. + """ + default = _da_lfc(GPSettings(**_NO_LM)) + control = _da_lfc(GPSettings(**_NO_LM), ls=0.05) + via_gp = _da_lfc(GPSettings(ls=0.05, **_NO_LM)) + assert np.abs(control - default).max() > 1e-1 + np.testing.assert_array_equal(via_gp, control) + + +def test_da_honours_gp_param_scheme(): + separate = _da_lfc(GPSettings(**_NO_LM)) + pooled = _da_lfc(GPSettings(param_scheme="pooled", **_NO_LM)) + with pytest.warns(DeprecationWarning, match="param_scheme"): + legacy = _da_lfc(GPSettings(**_NO_LM), sync_parameters=True) + assert np.abs(pooled - separate).max() > 1e-2 + np.testing.assert_array_equal(pooled, legacy) + + +def test_da_rejects_ls_given_twice_and_scheme_given_twice(): + with pytest.raises(ValueError, match="pass it once"): + _da_lfc(GPSettings(ls=0.05, **_NO_LM), ls=0.05) + with pytest.raises(ValueError, match="not both"): + _da_lfc(GPSettings(param_scheme="pooled", **_NO_LM), sync_parameters=True) + + +def test_da_warns_about_expression_only_gp_fields(): + with _captured_kompot_warnings() as messages: + _da_lfc(GPSettings(sigma=0.5, **_NO_LM)) + assert any("have no effect in da()" in m for m in messages), messages + with _captured_kompot_warnings() as messages: + _da_lfc(GPSettings(**_NO_LM)) + assert not [m for m in messages if "have no effect in da()" in m] + + +def test_the_da_cli_routes_param_scheme_and_ls_into_gp_settings(): + """Every GPSettings field the DA model uses is routed; the rest fall through.""" + import dataclasses + + from kompot.cli.da import GP_CONFIG_KEYS + + fields = {f.name for f in dataclasses.fields(GPSettings)} + assert {"param_scheme", "ls"} <= GP_CONFIG_KEYS <= fields + assert fields - GP_CONFIG_KEYS == {"sigma", "eps", "use_empirical_variance"} From 241b36db9936c5ca177c20f678ab51cc0954fe48 Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 13:06:17 -0700 Subject: [PATCH 11/23] docs(tests): wrap two docstring lines to the 100-column limit --- tests/test_da_swap_equivalence.py | 7 ++++--- tests/test_param_scheme.py | 11 ++++++----- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index 9e28222..4e864a7 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -15,9 +15,10 @@ inheritance to be asymmetric about. Two configurations DO break the exactness, and both are covered below: automatic -landmarks, and ``param_scheme="pooled"`` (formerly ``sync_parameters=True``). Both trace to the same ROOT as the DE -case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent -- but -they act through DIFFERENT proximate mechanisms, and **a remedy applies at the +landmarks, and ``param_scheme="pooled"`` (formerly ``sync_parameters=True``). +Both trace to the same ROOT as the DE case -- +``np.vstack([X_condition1, X_condition2])`` is row-order dependent -- but they +act through DIFFERENT proximate mechanisms, and **a remedy applies at the proximate level, not the root**: ``compute_landmarks`` (k-means init) for the first, ``compute_nn_distances`` (which sets ``mu`` and ``ls``) for the second. So sharing a landmark array fixes the first and cannot reach the second -- diff --git a/tests/test_param_scheme.py b/tests/test_param_scheme.py index 51b4354..549d238 100644 --- a/tests/test_param_scheme.py +++ b/tests/test_param_scheme.py @@ -4,11 +4,12 @@ it covers the length scale; in differential ABUNDANCE it covers ``d``, ``mu`` and ``ls`` (the second half of this file). -The expression default (``param_scheme=None`` -> ``"condition1"``) estimates the shared length scale from -condition 1's cells only, which makes the contrast depend on which condition is -passed first. These tests pin the default's behaviour and check that each -alternative does what it says -- including ``"condition2"``, the default's -mirror, for which the statement is an exact equivalence: +The expression default (``param_scheme=None``, i.e. ``"condition1"``) estimates +the shared length scale from condition 1's cells only, which makes the contrast +depend on which condition is passed first. These tests pin the default's +behaviour and check that each alternative does what it says -- including +``"condition2"``, the default's mirror, for which the statement is an exact +equivalence: ``de(X, Y, param_scheme="condition1")`` and ``de(Y, X, param_scheme="condition2")`` are the same computation with the labels exchanged. """ From 87708797fd9de85f0e35f4584d25d518613868f9 Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 13:20:13 -0700 Subject: [PATCH 12/23] test(da): take symmetric's reference bounds from the estimator recipe, not the condition2 scheme --- tests/test_param_scheme.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/test_param_scheme.py b/tests/test_param_scheme.py index 549d238..5e854bf 100644 --- a/tests/test_param_scheme.py +++ b/tests/test_param_scheme.py @@ -559,6 +559,7 @@ def test_the_other_schemes_do_not_warn_about_landmarks(): import kompot # noqa: E402 from kompot.differential.differential_abundance import ( # noqa: E402 DifferentialAbundance, + _condition_density_params, _resolve_scheme_density_params, ) from kompot.settings import OutputSettings # noqa: E402 @@ -725,8 +726,10 @@ def test_da_symmetric_is_bit_identical_under_a_swap(n1, n2): fwd = _resolve_scheme_density_params("symmetric", X1, X2, 10.0, 42) rev = _resolve_scheme_density_params("symmetric", X2, X1, 10.0, 42) assert fwd == rev - one = _resolve_scheme_density_params("condition1", X1, X2, 10.0, 42) - two = _resolve_scheme_density_params("condition2", X1, X2, 10.0, 42) + # per-condition references straight from the estimator recipe, so this test + # reads nothing of the one-sided schemes it is not about + one = _condition_density_params(X1, 10.0, 42)[0] + two = _condition_density_params(X2, 10.0, 42)[0] for name in ("d", "mu", "ls"): assert fwd[name] != one[name] and fwd[name] != two[name], name assert min(one[name], two[name]) <= fwd[name] <= max(one[name], two[name]), name From e9f798bd552f31651b1691c86502261631dfa65c Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 13:23:30 -0700 Subject: [PATCH 13/23] docs(de): point the de() gp entry at param_scheme and its default --- kompot/anndata/differential_expression.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/kompot/anndata/differential_expression.py b/kompot/anndata/differential_expression.py index bcd6e88..b106319 100644 --- a/kompot/anndata/differential_expression.py +++ b/kompot/anndata/differential_expression.py @@ -99,6 +99,10 @@ def de( Column with biological-replicate labels. gp : GPSettings, optional GP model parameters (sigma, ls, n_landmarks, etc.). + ``param_scheme`` sets which cells the shared length scale is + estimated from; ``None`` means ``"condition1"``, which makes the + result depend on which condition is passed first. See + :meth:`~kompot.differential.DifferentialExpression.fit`. fdr : FDRSettings, optional FDR / null-distribution parameters. filter : FilterSettings, optional From 298e032af4c4e3dd227be83ab85149638e41ea5f Mon Sep 17 00:00:00 2001 From: Dominik Date: Tue, 22 Sep 2026 14:46:18 -0700 Subject: [PATCH 14/23] fix(da): derive mu at a pinned d under every param_scheme; refuse flat param_scheme - With `d` pinned (alone), "condition1", "condition2" and "symmetric" computed the shared `mu` at the fractal `d` rather than the pinned one, so the shared value was not what the estimator derives. A pinned `d` now feeds `mu` (and, for "symmetric", replaces the pooled d). "pooled" was already correct. - Tests pin `d` alone and `ls` alone for each one-sided/symmetric scheme; the d-alone tests fail at e9f798b and pass here. - `param_scheme` passed to de()/da() as a flat keyword now raises a ValueError pointing at gp=GPSettings(param_scheme=...) instead of a TypeError. - CHANGELOG: the argument-order bullet moves from Added to Documentation. --- CHANGELOG.md | 12 ++-- kompot/anndata/differential_abundance.py | 5 ++ kompot/anndata/differential_expression.py | 5 ++ kompot/differential/differential_abundance.py | 34 ++++++---- tests/test_param_scheme.py | 62 +++++++++++++++++++ 5 files changed, 100 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 113dab4..d353a95 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,12 +12,6 @@ All notable changes to this project will be documented in this file. expression and `d`, `mu` and `ls` in differential abundance; an explicit value pins that parameter. The default, `None`, keeps each entry point's existing behaviour: `"condition1"` for `de()` and `"separate"` for `da()`. - - **Expression results depend on argument order at the default.** The shared length scale is - estimated from condition 1's cells only, so `de(condition1=X, condition2=Y)` is not equivalent - to `de(condition1=Y, condition2=X)`; on an exchangeable null the two orientations differ in - false-positive rate by 0.06–0.08 beyond 2:1. `"symmetric"` and `"pooled"` are swap-invariant, - `"separate"` is for diagnostics, and `"condition2"` mirrors the default so its order dependence - is visible from one call site. See `DifferentialExpression.fit`. - **Abundance gains `"condition1"`, `"condition2"` and `"symmetric"`.** `"symmetric"` is swap-invariant by construction. See `DifferentialAbundance.fit`. @@ -39,6 +33,12 @@ All notable changes to this project will be documented in this file. ### Documentation + - **Expression results depend on argument order at the default.** The shared length scale is + estimated from condition 1's cells only, so `de(condition1=X, condition2=Y)` is not equivalent + to `de(condition1=Y, condition2=X)`; on an exchangeable null the two orientations differ in + false-positive rate by 0.06–0.08 beyond 2:1. `"symmetric"` and `"pooled"` are swap-invariant, + `"separate"` is for diagnostics, and `"condition2"` mirrors the default so its order dependence + is visible from one call site. See `DifferentialExpression.fit`. - Abundance is swap-equivalent at `"separate"` without landmarks, and tested. `"pooled"` and automatic landmarks each break that, and need different remedies that do not substitute for one another: see `DifferentialAbundance.fit`. `da(adata, gp=GPSettings(...))` carries diff --git a/kompot/anndata/differential_abundance.py b/kompot/anndata/differential_abundance.py index d8f1872..f747909 100644 --- a/kompot/anndata/differential_abundance.py +++ b/kompot/anndata/differential_abundance.py @@ -124,6 +124,11 @@ def da( jit_compile = gp.jit_compile if gp is not None else False random_state = gp.random_state if gp is not None else None param_scheme = gp.param_scheme if gp is not None else None + if "param_scheme" in density_kwargs: + raise ValueError( + "`param_scheme` is a GPSettings field: pass " + "gp=GPSettings(param_scheme=...), not a keyword argument of da()." + ) if gp is not None: _unused = [ diff --git a/kompot/anndata/differential_expression.py b/kompot/anndata/differential_expression.py index b106319..74eb2c3 100644 --- a/kompot/anndata/differential_expression.py +++ b/kompot/anndata/differential_expression.py @@ -138,6 +138,11 @@ def de( sigma = _gp.sigma ls = _gp.ls ls_factor = _gp.ls_factor + if "param_scheme" in function_kwargs: + raise ValueError( + "`param_scheme` is a GPSettings field: pass " + "gp=GPSettings(param_scheme=...), not a keyword argument of de()." + ) param_scheme = resolve_param_scheme(_gp.param_scheme, DEFAULT_PARAM_SCHEME) n_landmarks = _gp.n_landmarks landmarks = _gp.landmarks diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index cf0ad7c..171dd63 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -27,27 +27,31 @@ _D_FACTAL_QUERY_CELLS = 500 -def _condition_density_params(X, ls_factor, seed): +def _condition_density_params(X, ls_factor, seed, d=None): """``d``, ``mu`` and ``ls`` exactly as a ``mellon.DensityEstimator`` derives them. Mirrors ``DensityEstimator._compute_d`` (``d_method="fractal"``), ``_compute_nn_distances``, ``_compute_mu`` and ``_compute_ls``, so that a value estimated here *before* either estimator is fitted is the value the estimator would have reached on its own. Returns the validated - nearest-neighbour distances as well, for ``"symmetric"``. + nearest-neighbour distances as well, for ``"symmetric"``. ``d``, when + given, is a pinned value and replaces the fractal estimate everywhere. """ from mellon.validation import validate_nn_distances X = np.asarray(X) nn = validate_nn_distances(mellon.parameters.compute_nn_distances(X, seed=seed)) - d = mellon.parameters.compute_d_factal(X) + # A pinned ``d`` is what the estimator will use, so ``mu`` is derived at it + # -- never at the fractal estimate the pin replaces. + if d is None: + d = mellon.parameters.compute_d_factal(X) mu = mellon.parameters.compute_mu(nn, d) ls = mellon.parameters.compute_ls(nn) * ls_factor return {"d": d, "mu": mu, "ls": ls}, nn def _resolve_scheme_density_params( - param_scheme, X_condition1, X_condition2, ls_factor, seed + param_scheme, X_condition1, X_condition2, ls_factor, seed, d=None ): """Return the shared ``{d, mu, ls}`` a one-sided or symmetric scheme prescribes. @@ -55,9 +59,9 @@ def _resolve_scheme_density_params( three schemes that estimate from each condition on its own. """ if param_scheme == "condition1": - return _condition_density_params(X_condition1, ls_factor, seed)[0] + return _condition_density_params(X_condition1, ls_factor, seed, d)[0] if param_scheme == "condition2": - return _condition_density_params(X_condition2, ls_factor, seed)[0] + return _condition_density_params(X_condition2, ls_factor, seed, d)[0] if param_scheme == "symmetric": # Each shared value is the model's own estimator applied to the two # conditions' WITHIN-condition statistics pooled together -- never to @@ -65,8 +69,8 @@ def _resolve_scheme_density_params( # combination below is written so that exchanging the conditions # exchanges two commutative terms, which makes the result bit-identical # under a swap by construction, not merely close. - p1, nn1 = _condition_density_params(X_condition1, ls_factor, seed) - p2, nn2 = _condition_density_params(X_condition2, ls_factor, seed) + p1, nn1 = _condition_density_params(X_condition1, ls_factor, seed, d) + p2, nn2 = _condition_density_params(X_condition2, ls_factor, seed, d) n1, n2 = np.asarray(X_condition1).shape[0], np.asarray(X_condition2).shape[0] # ls: geometric mean of nn distances pooled == size-weighted geometric # mean of the per-condition values (the differential-expression rule). @@ -75,9 +79,10 @@ def _resolve_scheme_density_params( ) # d: mean local dimensionality pooled over the query cells each # condition's estimate averaged. - q1 = min(n1, _D_FACTAL_QUERY_CELLS) - q2 = min(n2, _D_FACTAL_QUERY_CELLS) - d = float((q1 * p1["d"] + q2 * p2["d"]) / (q1 + q2)) + if d is None: + q1 = min(n1, _D_FACTAL_QUERY_CELLS) + q2 = min(n2, _D_FACTAL_QUERY_CELLS) + d = float((q1 * p1["d"] + q2 * p2["d"]) / (q1 + q2)) # mu: the 1st-percentile rule over the pooled within-condition nn # distances, at the shared d. A quantile sorts, so the order the two # conditions are concatenated in does not reach the result. @@ -444,7 +449,12 @@ def fit( if seed is None: seed = DEFAULT_RANDOM_SEED shared = _resolve_scheme_density_params( - param_scheme, X_condition1, X_condition2, ls_factor, seed + param_scheme, + X_condition1, + X_condition2, + ls_factor, + seed, + d=density_kwargs.get("d"), ) for name in unpinned: estimator_defaults[name] = shared[name] diff --git a/tests/test_param_scheme.py b/tests/test_param_scheme.py index 5e854bf..0ea2f9f 100644 --- a/tests/test_param_scheme.py +++ b/tests/test_param_scheme.py @@ -689,6 +689,60 @@ def test_da_explicit_values_pin_every_scheme(scheme): assert s1 == s2 == (2.0, -5.0, 3.0) +_ONE_SIDED = [("condition1", 0), ("condition2", 1)] + + +@pytest.mark.parametrize("scheme, side", _ONE_SIDED) +def test_da_one_sided_scheme_derives_mu_at_a_pinned_d(scheme, side): + """Pinning ``d`` ALONE must reach ``mu``, which the estimator derives from ``d``. + + Reference: the unshared fit with the same pin, i.e. what that condition's + estimator does on its own. Pinning d, mu and ls together (the test above) + cannot see this, because the pinned mu hides the derived one. + """ + X1, X2 = _da_data() + sep = _shared(_da_fit(X1, X2, param_scheme="separate", d=2.0))[side] + got = _shared(_da_fit(X1, X2, param_scheme=scheme, d=2.0)) + assert got[0] == got[1] + assert got[side] == sep + assert got[side][0] == 2.0 + # control: the pin moves mu, so equality above is not the unpinned value + assert got[side][1] != _shared(_da_fit(X1, X2, param_scheme=scheme))[side][1] + + +def test_da_symmetric_derives_mu_at_a_pinned_d(): + from mellon.parameters import compute_mu + + X1, X2 = _da_data() + s1, s2 = _shared(_da_fit(X1, X2, param_scheme="symmetric", d=2.0)) + assert s1 == s2 + nn1 = _condition_density_params(X1, 10.0, 42)[1] + nn2 = _condition_density_params(X2, 10.0, 42)[1] + assert s1[0] == 2.0 + assert s1[1] == float(compute_mu(np.concatenate([nn1, nn2]), 2.0)) + unpinned = _resolve_scheme_density_params("symmetric", X1, X2, 10.0, 42) + assert s1[2] == float(unpinned["ls"]) # ls does not depend on d + assert s1[1] != float(unpinned["mu"]) # control: the pin moves mu + + +@pytest.mark.parametrize("scheme, side", _ONE_SIDED) +def test_da_one_sided_scheme_with_ls_pinned_alone(scheme, side): + """Control for the d-pin tests: ``ls`` does not enter d or mu.""" + X1, X2 = _da_data() + sep = _shared(_da_fit(X1, X2, param_scheme="separate", ls=3.0))[side] + got = _shared(_da_fit(X1, X2, param_scheme=scheme, ls=3.0)) + assert got[0] == got[1] + assert got[side] == sep + assert got[side][2] == 3.0 + + +def test_da_symmetric_with_ls_pinned_alone(): + X1, X2 = _da_data() + s1, s2 = _shared(_da_fit(X1, X2, param_scheme="symmetric", ls=3.0)) + unpinned = _resolve_scheme_density_params("symmetric", X1, X2, 10.0, 42) + assert s1 == s2 == (float(unpinned["d"]), float(unpinned["mu"]), 3.0) + + # -- swap behaviour --------------------------------------------------------- @@ -812,6 +866,14 @@ def test_da_rejects_ls_given_twice_and_scheme_given_twice(): _da_lfc(GPSettings(param_scheme="pooled", **_NO_LM), sync_parameters=True) +def test_param_scheme_as_a_flat_keyword_is_refused_clearly(): + """It is a GPSettings field; as a keyword it used to die in a TypeError.""" + with pytest.raises(ValueError, match=r"gp=GPSettings\(param_scheme"): + _da_lfc(GPSettings(**_NO_LM), param_scheme="pooled") + with pytest.raises(ValueError, match=r"gp=GPSettings\(param_scheme"): + kompot.de(_adata(), "cond", "A", "B", obsm_key="X", param_scheme="pooled") + + def test_da_warns_about_expression_only_gp_fields(): with _captured_kompot_warnings() as messages: _da_lfc(GPSettings(sigma=0.5, **_NO_LM)) From 7e37ac65cf4acd7598b8bc202f4f4bb6a812cf6b Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 10:34:19 -0700 Subject: [PATCH 15/23] refactor(da)!: remove sync_parameters; param_scheme takes its place No alias and no deprecation shim. param_scheme="pooled" is the former sync_parameters=True. param_scheme occupies the removed argument's position in DifferentialAbundance.fit, so allow_single_condition_variance keeps the position it has on main. --- kompot/anndata/differential_abundance.py | 20 --------- kompot/cli/da.py | 3 +- kompot/differential/differential_abundance.py | 42 +++++-------------- tests/test_da_swap_equivalence.py | 2 +- tests/test_param_scheme.py | 41 ++---------------- 5 files changed, 15 insertions(+), 93 deletions(-) diff --git a/kompot/anndata/differential_abundance.py b/kompot/anndata/differential_abundance.py index f747909..577032d 100644 --- a/kompot/anndata/differential_abundance.py +++ b/kompot/anndata/differential_abundance.py @@ -50,7 +50,6 @@ def da( storage: "StorageSettings | None" = None, output: "OutputSettings | None" = None, model: "ModelSettings | None" = None, - sync_parameters: "bool | None" = None, **density_kwargs, ) -> "Union[Dict[str, np.ndarray], Any]": """Run differential abundance analysis on an AnnData object. @@ -99,11 +98,6 @@ def da( Pre-fitted models or predictors to inject. For DA, only ``density_predictor1/2`` and ``variance_predictor1/2`` are used. See :class:`~kompot.ModelSettings`. - sync_parameters : bool, optional - .. deprecated:: - Use ``gp=GPSettings(param_scheme=...)``. ``True`` is - ``"pooled"`` and ``False`` is ``"separate"``; passing both - raises. **density_kwargs Forwarded to :class:`~mellon.DensityEstimator`. ``d``, ``mu`` and ``ls`` given here pin that parameter under any ``param_scheme``. @@ -149,20 +143,6 @@ def da( ) density_kwargs["ls"] = gp.ls - if sync_parameters is not None: - if param_scheme is not None: - raise ValueError( - "Pass either GPSettings.param_scheme or the deprecated " - "`sync_parameters`, not both." - ) - warnings.warn( - "`sync_parameters` is deprecated; use " - "gp=GPSettings(param_scheme='pooled') for True or " - "param_scheme='separate' for False.", - DeprecationWarning, - stacklevel=2, - ) - param_scheme = "pooled" if sync_parameters else "separate" param_scheme = resolve_param_scheme(param_scheme, DEFAULT_PARAM_SCHEME) _threshold = threshold if threshold is not None else DAThresholdSettings() diff --git a/kompot/cli/da.py b/kompot/cli/da.py index 2c82f47..ce2b3d9 100644 --- a/kompot/cli/da.py +++ b/kompot/cli/da.py @@ -187,8 +187,7 @@ def add_da_parser(subparsers) -> argparse.ArgumentParser: #: Flat config keys the DA CLI routes into :class:`~kompot.settings.GPSettings`. #: Every other key -- including ``sigma``, which here means mellon's density #: noise rather than the expression model's -- is forwarded to ``da()`` and -#: from there to :class:`mellon.DensityEstimator`. The deprecated -#: ``sync_parameters`` key reaches ``da()``'s named, deprecated parameter. +#: from there to :class:`mellon.DensityEstimator`. GP_CONFIG_KEYS = { "ls", "ls_factor", diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 171dd63..929a615 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -3,7 +3,6 @@ import numpy as np from typing import Optional, Dict, Any import logging -import warnings from scipy.stats import norm as normal import mellon @@ -236,9 +235,8 @@ def fit( condition1_sample_indices: Optional[np.ndarray] = None, condition2_sample_indices: Optional[np.ndarray] = None, sample_estimator_ls: Optional[float] = None, - sync_parameters: Optional[bool] = None, - allow_single_condition_variance: bool = False, param_scheme: Optional[str] = None, + allow_single_condition_variance: bool = False, **density_kwargs, ): """ @@ -278,31 +276,23 @@ def fit( sample_estimator_ls : float, optional Length scale for the sample-specific variance estimators. If None, will use the same value as ls or it will be estimated, by default None. - sync_parameters : bool, optional - .. deprecated:: - Use ``param_scheme``. ``True`` is ``param_scheme="pooled"`` and - ``False`` is ``param_scheme="separate"``; passing both raises. - allow_single_condition_variance : bool, optional - Allow sample-variance estimation when only one condition has - multiple samples, by default False. param_scheme : str, optional Where the density hyperparameters ``d``, ``mu`` and ``ls`` come from. - The same field drives :meth:`DifferentialExpression.fit`, where it - covers ``ls`` only. A value passed explicitly through - ``density_kwargs`` (``d=``, ``mu=``, ``ls=``) pins that parameter - and takes precedence. ``None`` (default) means ``"separate"``. + A value passed explicitly through ``density_kwargs`` (``d=``, + ``mu=``, ``ls=``) pins that parameter and takes precedence. ``None`` + (default) means ``"separate"``. * ``"separate"`` (default) — each density estimator derives its own ``d``, ``mu`` and ``ls`` from its own condition; nothing is shared. * ``"pooled"`` — estimate all three once from both conditions' cells - taken together and share them (formerly ``sync_parameters=True``). + taken together and share them. * ``"condition1"`` / ``"condition2"`` — estimate all three from one condition's cells, exactly as that condition's own estimator would, and share them. As asymmetric as the names say. * ``"symmetric"`` — estimate from each condition separately and share the pooled estimate: ``ls`` is the size-weighted geometric mean of - the two (as in differential expression), ``d`` the mean local - dimensionality over both conditions' query cells, and ``mu`` the + the two, ``d`` the mean local dimensionality over both + conditions' query cells, and ``mu`` the 1st-percentile rule over both conditions' within-condition nearest-neighbour distances at that ``d``. Bit-identical under swapping the conditions. @@ -335,6 +325,9 @@ def fit( half applies to them: ``da(X, Y, "condition1")`` and ``da(Y, X, "condition2")`` are the same computation with the labels exchanged when both runs use the same landmarks. + allow_single_condition_variance : bool, optional + Allow sample-variance estimation when only one condition has + multiple samples, by default False. **density_kwargs : dict Additional arguments to pass to the DensityEstimator. @@ -344,20 +337,6 @@ def fit( The fitted instance. """ - if sync_parameters is not None: - if param_scheme is not None: - raise ValueError( - "Pass either `param_scheme` or the deprecated " - "`sync_parameters`, not both." - ) - warnings.warn( - "`sync_parameters` is deprecated; use " - "`param_scheme='pooled'` (for True) or " - "`param_scheme='separate'` (for False).", - DeprecationWarning, - stacklevel=2, - ) - param_scheme = "pooled" if sync_parameters else "separate" param_scheme = resolve_param_scheme(param_scheme, DEFAULT_PARAM_SCHEME) # Create or use density predictors @@ -405,7 +384,6 @@ def fit( self.computed_landmarks = computed_landmarks # "pooled": estimate d / mu / ls from the stacked union and share. - # This block is the former `sync_parameters=True` path, unchanged. if param_scheme == "pooled": # Combine data from both conditions for parameter estimation X_combined = np.vstack([X_condition1, X_condition2]) diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index 4e864a7..c73b769 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -15,7 +15,7 @@ inheritance to be asymmetric about. Two configurations DO break the exactness, and both are covered below: automatic -landmarks, and ``param_scheme="pooled"`` (formerly ``sync_parameters=True``). +landmarks, and ``param_scheme="pooled"``. Both trace to the same ROOT as the DE case -- ``np.vstack([X_condition1, X_condition2])`` is row-order dependent -- but they act through DIFFERENT proximate mechanisms, and **a remedy applies at the diff --git a/tests/test_param_scheme.py b/tests/test_param_scheme.py index 0ea2f9f..2c87fd0 100644 --- a/tests/test_param_scheme.py +++ b/tests/test_param_scheme.py @@ -552,8 +552,6 @@ def test_the_other_schemes_do_not_warn_about_landmarks(): # differential ABUNDANCE: the same field covers d, mu and ls # ========================================================================== -import warnings # noqa: E402 - import anndata # noqa: E402 import kompot # noqa: E402 @@ -579,9 +577,7 @@ def _da_data(n1=90, n2=60, seed=0): def _da_fit(X1, X2, **fit_kwargs): """Fit with no landmarks, so the only orientation effect is the scheme's.""" model = DifferentialAbundance(n_landmarks=None, random_state=0) - with warnings.catch_warnings(): - warnings.simplefilter("ignore", DeprecationWarning) - model.fit(X1, X2, **fit_kwargs) + model.fit(X1, X2, **fit_kwargs) return model @@ -597,7 +593,7 @@ def _shared(model): ) -# -- defaults and the deprecated alias -------------------------------------- +# -- defaults -------------------------------------------------------------- def test_da_default_is_separate(): @@ -609,27 +605,7 @@ def test_da_default_is_separate(): np.testing.assert_array_equal(a[key], b[key]) -@pytest.mark.parametrize("sync, scheme", [(True, "pooled"), (False, "separate")]) -def test_sync_parameters_is_an_exact_alias(sync, scheme): - """The shipped spelling and the new one are the same computation, to the bit.""" - X1, X2 = _da_data() - Xn = np.vstack([X1, X2]) - model = DifferentialAbundance(n_landmarks=None, random_state=0) - with pytest.warns(DeprecationWarning, match="param_scheme"): - model.fit(X1, X2, sync_parameters=sync) - old = _da_predict(model, Xn) - new = _da_predict(_da_fit(X1, X2, param_scheme=scheme), Xn) - for key in ( - "log_density_condition1", - "log_density_condition2", - "log_fold_change", - "log_fold_change_uncertainty", - ): - np.testing.assert_array_equal(old[key], new[key]) - - def test_pooled_and_separate_differ(): - """Positive control for the alias test: the two targets are distinguishable.""" X1, X2 = _da_data() Xn = np.vstack([X1, X2]) a = _da_predict(_da_fit(X1, X2, param_scheme="pooled"), Xn) @@ -637,12 +613,6 @@ def test_pooled_and_separate_differ(): assert np.abs(a["log_fold_change"] - b["log_fold_change"]).max() > 1e-2 -def test_sync_parameters_and_param_scheme_together_raise(): - X1, X2 = _da_data() - with pytest.raises(ValueError, match="not both"): - DifferentialAbundance().fit(X1, X2, sync_parameters=True, param_scheme="pooled") - - def test_da_rejects_unknown_scheme(): X1, X2 = _da_data() with pytest.raises(ValueError, match="param_scheme"): @@ -853,17 +823,12 @@ def test_da_honours_gp_ls(): def test_da_honours_gp_param_scheme(): separate = _da_lfc(GPSettings(**_NO_LM)) pooled = _da_lfc(GPSettings(param_scheme="pooled", **_NO_LM)) - with pytest.warns(DeprecationWarning, match="param_scheme"): - legacy = _da_lfc(GPSettings(**_NO_LM), sync_parameters=True) assert np.abs(pooled - separate).max() > 1e-2 - np.testing.assert_array_equal(pooled, legacy) -def test_da_rejects_ls_given_twice_and_scheme_given_twice(): +def test_da_rejects_ls_given_twice(): with pytest.raises(ValueError, match="pass it once"): _da_lfc(GPSettings(ls=0.05, **_NO_LM), ls=0.05) - with pytest.raises(ValueError, match="not both"): - _da_lfc(GPSettings(param_scheme="pooled", **_NO_LM), sync_parameters=True) def test_param_scheme_as_a_flat_keyword_is_refused_clearly(): From 12e1253e5b5bede0dd2b63d3aa8999f68e6c6df6 Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 10:34:19 -0700 Subject: [PATCH 16/23] refactor(anndata): overwrite checks take param_scheme from their caller _check_da_overwrites and _check_overwrites carried their own hard-coded defaults ("separate", "condition1"), a second copy of each entry point's default that could drift from it. param_scheme is now a required argument, and a run recorded before param_scheme existed reads back as DEFAULT_PARAM_SCHEME. --- kompot/anndata/_da_helpers.py | 14 +++++--------- kompot/anndata/_de_helpers.py | 11 ++++++----- tests/test_da_helpers_overwrite.py | 6 ++++++ 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/kompot/anndata/_da_helpers.py b/kompot/anndata/_da_helpers.py index ba9b263..f2f833b 100644 --- a/kompot/anndata/_da_helpers.py +++ b/kompot/anndata/_da_helpers.py @@ -11,6 +11,8 @@ import pandas as pd from typing import Optional +from ..differential.differential_abundance import DEFAULT_PARAM_SCHEME + logger = logging.getLogger("kompot") @@ -25,7 +27,7 @@ def _check_da_overwrites( condition2: str, obsm_key: str, ls_factor: float, - param_scheme: str = "separate", + param_scheme: str, ): """Check for existing DA results and handle overwrite logic.""" from .utils import detect_output_field_overwrite @@ -107,14 +109,8 @@ def _check_da_overwrites( curr_val = locals().get(param_name) prev_val = _pg(prev_params, param_name) if param_name == "param_scheme" and prev_val is None: - # Recorded before `param_scheme` existed: sharing was - # then only reachable as a `sync_parameters` kwarg. - prev_extra = prev_params.get("extra_kwargs") or {} - prev_val = ( - "pooled" - if prev_extra.get("sync_parameters") - else "separate" - ) + # Recorded before `param_scheme` existed. + prev_val = DEFAULT_PARAM_SCHEME if curr_val != prev_val: params_match = False logger.debug( diff --git a/kompot/anndata/_de_helpers.py b/kompot/anndata/_de_helpers.py index f17190a..4b8c539 100644 --- a/kompot/anndata/_de_helpers.py +++ b/kompot/anndata/_de_helpers.py @@ -13,6 +13,7 @@ from typing import Optional, Dict, List, Tuple from scipy import sparse +from ..differential.differential_expression import DEFAULT_PARAM_SCHEME from .utils import check_underrepresentation, apply_cell_filter logger = logging.getLogger("kompot") @@ -41,7 +42,7 @@ def _check_overwrites( obsm_key: str, layer: Optional[str], ls_factor: float, - param_scheme: str = "condition1", + param_scheme: str, ): """Check for existing results and handle overwrite logic. @@ -138,11 +139,11 @@ def _check_overwrites( "param_scheme", ]: curr_val = locals().get(param_name) - # A run recorded before `param_scheme` existed used "condition1". + # A run recorded before `param_scheme` existed used the default. prev_val = _pg( prev_params, param_name, - "condition1" if param_name == "param_scheme" else None, + DEFAULT_PARAM_SCHEME if param_name == "param_scheme" else None, ) if curr_val != prev_val: params_match = False @@ -199,11 +200,11 @@ def _check_overwrites( "param_scheme", ]: curr_val = locals().get(param_name) - # A run recorded before `param_scheme` existed used "condition1". + # A run recorded before `param_scheme` existed used the default. prev_val = _pg( prev_params, param_name, - "condition1" if param_name == "param_scheme" else None, + DEFAULT_PARAM_SCHEME if param_name == "param_scheme" else None, ) if curr_val != prev_val: params_match = False diff --git a/tests/test_da_helpers_overwrite.py b/tests/test_da_helpers_overwrite.py index 2ab0df2..7d81a03 100644 --- a/tests/test_da_helpers_overwrite.py +++ b/tests/test_da_helpers_overwrite.py @@ -109,6 +109,7 @@ def test_overwrite_false_raises(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) def test_overwrite_none_with_sample_var_change_params_match(self): @@ -131,6 +132,7 @@ def test_overwrite_none_with_sample_var_change_params_match(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) def test_overwrite_none_with_sample_var_change_params_mismatch(self): @@ -153,6 +155,7 @@ def test_overwrite_none_with_sample_var_change_params_mismatch(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) def test_overwrite_false_with_sample_var_fields_listed(self): @@ -175,6 +178,7 @@ def test_overwrite_false_with_sample_var_fields_listed(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) def test_overwrite_false_prev_sample_var_true_current_false(self): @@ -195,6 +199,7 @@ def test_overwrite_false_prev_sample_var_true_current_false(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) def test_overwrite_none_no_sample_var_change(self): @@ -214,6 +219,7 @@ def test_overwrite_none_no_sample_var_change(self): condition2="B", obsm_key="X_pca", ls_factor=1.0, + param_scheme="separate", ) From 812c15b9eb3b2fe3ec394e2e1c2e8f8fa1555e51 Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 10:34:19 -0700 Subject: [PATCH 17/23] docs: state only where DA differs; shorten the changelog --- CHANGELOG.md | 42 ++++++------------- docs/source/cli.rst | 4 +- .../differential/differential_expression.py | 8 ++-- kompot/settings.py | 14 +++---- 4 files changed, 23 insertions(+), 45 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d353a95..1126187 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,44 +6,26 @@ All notable changes to this project will be documented in this file. ### Added - - **`GPSettings.param_scheme`**: one field, read by both `kompot.de()` and `kompot.da()`, that - says where the hyperparameters estimated from cells come from — `"condition1"`, - `"condition2"`, `"symmetric"`, `"pooled"` or `"separate"`. It covers `ls` in differential - expression and `d`, `mu` and `ls` in differential abundance; an explicit value pins that - parameter. The default, `None`, keeps each entry point's existing behaviour: `"condition1"` for - `de()` and `"separate"` for `da()`. - - **Abundance gains `"condition1"`, `"condition2"` and `"symmetric"`.** `"symmetric"` is - swap-invariant by construction. See `DifferentialAbundance.fit`. + - **`GPSettings.param_scheme`** chooses which cells the length scale is estimated from: + `"condition1"`, `"condition2"`, `"symmetric"`, `"pooled"` or `"separate"`. In `da()` it covers + `d` and `mu` as well. The default, `None`, keeps current behaviour: `"condition1"` in `de()`, + `"separate"` in `da()`. `"symmetric"` is swap-invariant. Both CLIs accept it as a config key. -### Changed - - - `param_scheme` participates in run-parameter matching for both `de()` and `da()`. - -### Deprecated +### Removed - - **`sync_parameters`** on `DifferentialAbundance.fit` and `da()`: `True` is - `param_scheme="pooled"` and `False` is `"separate"`, with identical results. It previously - reached `da()` only through `**density_kwargs`. + - `sync_parameters` on `DifferentialAbundance.fit` and `da()`; use `param_scheme="pooled"`. ### Fixed - - **`da()` now honours `GPSettings.ls`.** It read six `GPSettings` fields and silently ignored the - rest. The expression-only fields `sigma`, `eps` and `use_empirical_variance` now log a warning - when set, and the DA CLI routes `ls` and `param_scheme` config keys into `GPSettings`. + - `da()` now honours `GPSettings.ls`, and warns when the expression-only fields `sigma`, `eps` or + `use_empirical_variance` are set. ### Documentation - - **Expression results depend on argument order at the default.** The shared length scale is - estimated from condition 1's cells only, so `de(condition1=X, condition2=Y)` is not equivalent - to `de(condition1=Y, condition2=X)`; on an exchangeable null the two orientations differ in - false-positive rate by 0.06–0.08 beyond 2:1. `"symmetric"` and `"pooled"` are swap-invariant, - `"separate"` is for diagnostics, and `"condition2"` mirrors the default so its order dependence - is visible from one call site. See `DifferentialExpression.fit`. - - Abundance is swap-equivalent at `"separate"` without landmarks, and tested. `"pooled"` and - automatic landmarks each break that, and need different remedies that do not substitute for - one another: see `DifferentialAbundance.fit`. `da(adata, gp=GPSettings(...))` carries - `n_landmarks=5000`. A pre-existing instability at `n_landmarks >= n_combined` is tracked at - settylab/kompot#32. + - At the default, `de()` depends on argument order: the shared length scale comes from condition + 1 alone. See `DifferentialExpression.fit`. + - Automatic landmarks and `"pooled"` make `da()` depend on argument order; see + `DifferentialAbundance.fit`. The instability at `n_landmarks >= n_combined` is settylab/kompot#32. ## [0.8.0] - 2026-07-28 diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 1f84983..6929d90 100644 --- a/docs/source/cli.rst +++ b/docs/source/cli.rst @@ -443,9 +443,9 @@ The table below shows how config keys correspond to Python Settings: * - Config key - Python equivalent - Description - * - ``ls``, ``ls_factor``, ``param_scheme``, ``n_landmarks``, ``batch_size``, ``jit_compile``, ``random_state``; DE also ``sigma``, ``eps``, ``use_empirical_variance`` + * - ``sigma``, ``ls``, ``ls_factor``, ``param_scheme``, ``n_landmarks``, ``batch_size``, ``eps``, ``jit_compile``, ``random_state``, ``use_empirical_variance`` - ``gp=GPSettings(...)`` - - GP model parameters. ``param_scheme`` defaults to ``condition1`` in DE and ``separate`` in DA + - GP model parameters. In DA, ``sigma`` is the density estimator's and ``eps`` and ``use_empirical_variance`` do not apply. ``param_scheme`` defaults to ``condition1`` in DE and ``separate`` in DA * - ``null_genes``, ``null_seed``, ``fdr_threshold`` - ``fdr=FDRSettings(...)`` - FDR / null distribution (DE) diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index b9ef5bd..c34d5fa 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -414,11 +414,9 @@ def fit( Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is None. param_scheme : str, optional - Where the hyperparameters estimated from cells come from. The same - field drives :meth:`DifferentialAbundance.fit`; here the only such - hyperparameter is the length scale, so the scheme decides which cells - the shared ``ls`` is estimated from when ``ls`` is None. Both - conditions are normally smoothed at the *same* scale so that their + Which cells the shared length scale is estimated from when ``ls`` + is None. Both conditions are normally smoothed at the *same* scale + so that their fitted surfaces are comparable; the schemes differ in which cells the shared value is estimated from. ``None`` (default) means ``"condition1"``. diff --git a/kompot/settings.py b/kompot/settings.py index 19ff328..b34e073 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -65,12 +65,10 @@ class GPSettings: ls_factor : float Multiplier applied to the automatically inferred length scale. param_scheme : str, optional - Where the hyperparameters a model estimates from cells come from. - Read by both :func:`kompot.de` and :func:`kompot.da`; the model - decides which parameters it covers — ``ls`` in differential - expression, ``d``, ``mu`` and ``ls`` in differential abundance. A - value given explicitly (``ls`` here, or ``d`` / ``mu`` passed to - ``da()``) pins that parameter and takes precedence. + Which cells the length scale is estimated from; in :func:`kompot.da` + it covers ``d`` and ``mu`` as well. A value given explicitly (``ls`` + here, or ``d`` / ``mu`` passed to ``da()``) pins that parameter and + takes precedence. ``None`` (default) keeps each entry point's own default: ``"condition1"`` for ``de()`` and ``"separate"`` for ``da()``. @@ -80,8 +78,8 @@ class GPSettings: first. * ``"condition2"`` — the mirror: estimate from condition 2's cells and reuse for condition 1. Equally asymmetric, by design; it exists so - that ``de(X, Y, "condition1")`` and ``de(Y, X, "condition2")`` (and - likewise ``da``) are the same computation with the labels exchanged + that ``de(X, Y, "condition1")`` and ``de(Y, X, "condition2")`` are + the same computation with the labels exchanged **provided both runs use the same landmarks**. The defaults here (``n_landmarks=5000``, ``landmarks=None``) do *not* satisfy that and the equivalence then holds only approximately — see From 52c84843bdf366451cfe35b9c81750b57dc5d5c2 Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 11:52:05 -0700 Subject: [PATCH 18/23] fix(da): n_landmarks >= n_combined fits the full GP, not the union as landmarks With as many landmarks as cells, DifferentialAbundance used the stacked union [X_condition1; X_condition2] as the landmark array. mellon's diagonal Laplace uncertainty depends on the landmarks' row order, and that order is set by which condition is passed first, so da(X, Y) and da(Y, X) disagreed on the uncertainty (median 97%) while the fold changes agreed to 1e-10. Any GPSettings routes da() here on data under 5000 cells. No approximation was asked for, so each condition now gets its full GP on its own cells, as n_landmarks=None does. Fixes #32 --- kompot/differential/differential_abundance.py | 19 ++++++- tests/test_da_swap_equivalence.py | 56 +++++++++++++++++++ 2 files changed, 72 insertions(+), 3 deletions(-) diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 929a615..5ddcf51 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -154,6 +154,8 @@ def __init__( Threshold for considering a PTP significant, by default 1e-2. n_landmarks : int, optional Number of landmarks to use for approximation. If None, use all points, by default None. + At or above the number of cells in both conditions together, no landmarks + are built and each condition gets its full GP, as with None. use_sample_variance : bool, optional Whether to use sample variance for uncertainty estimation. By default None. - If None (recommended): Automatically determined based on variance_predictor1/2 @@ -366,11 +368,22 @@ def fit( estimator_defaults["gp_type"] = "fixed" # Store provided landmarks for future use self.computed_landmarks = landmarks + elif self.n_landmarks is not None and self.n_landmarks >= ( + len(X_condition1) + len(X_condition2) + ): + # As many landmarks as cells would make the stacked union itself + # the landmark array, ordered by which condition was passed first. + # mellon's Laplace uncertainty depends on the landmarks' row order, + # so da(X, Y) and da(Y, X) then disagree on the uncertainty + # (settylab/kompot#32). No approximation was asked for: fit each + # condition's exact GP, on its own cells, as n_landmarks=None does. + logger.info( + f"n_landmarks={self.n_landmarks:,} >= {len(X_condition1) + len(X_condition2):,} " + "cells: fitting each condition's full GP instead of landmarks." + ) + estimator_defaults.setdefault("n_landmarks", 0) elif self.n_landmarks is not None: # Use mellon's compute_landmarks function to get properly distributed landmarks. - # gp_type="fixed" guarantees a populated landmark array even when - # n_landmarks >= n_combined, so both per-condition density estimators - # share the same inducing points. X_combined = np.vstack([X_condition1, X_condition2]) computed_landmarks = compute_landmarks( X_combined, diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py index c73b769..15bd9bc 100644 --- a/tests/test_da_swap_equivalence.py +++ b/tests/test_da_swap_equivalence.py @@ -253,3 +253,59 @@ def test_shared_landmarks_restore_the_exact_swap(): np.testing.assert_allclose( fwd["neg_log10_fold_change_ptp"], rev["neg_log10_fold_change_ptp"], **_TOL ) + + +# -- settylab/kompot#32: as many landmarks as cells --------------------------- + +_OUTPUTS = ( + "log_density_condition1", + "log_density_condition2", + "log_fold_change", + "log_fold_change_uncertainty", + "log_fold_change_zscore", + "neg_log10_fold_change_ptp", +) + + +@pytest.mark.parametrize("n_landmarks", [5000, 800], ids=["above", "equal"]) +def test_landmarks_covering_every_cell_fit_the_full_gp(n_landmarks): + """``n_landmarks >= n_combined`` is the full GP, so it is swap-symmetric. + + It used to make the stacked union the landmark array. mellon's Laplace + uncertainty depends on the landmarks' ROW ORDER, and that order is set by + which condition is passed first, so the two orientations disagreed on the + uncertainty by a median of 97% while the fold changes agreed to 1e-10. + ``any GPSettings`` routes ``da()`` here on data under 5000 cells. + """ + X1, X2 = _data() + X_new = np.vstack([X1, X2]) + assert n_landmarks >= len(X_new) + + full = DifferentialAbundance(n_landmarks=None, random_state=0) + full.fit(X1, X2) + ref = full.predict(X_new, progress=False) + + model = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) + model.fit(X1, X2) + assert getattr(model, "computed_landmarks", None) is None + got = model.predict(X_new, progress=False) + for key in _OUTPUTS: + np.testing.assert_array_equal(np.asarray(got[key]), np.asarray(ref[key]), err_msg=key) + + fwd, rev = _swap_pair(n_landmarks=n_landmarks) + np.testing.assert_allclose( + fwd["log_fold_change_uncertainty"], rev["log_fold_change_uncertainty"], **_TOL + ) + np.testing.assert_allclose( + np.asarray(fwd["log_fold_change_zscore"]), + -np.asarray(rev["log_fold_change_zscore"]), + **_TOL, + ) + + +def test_fewer_landmarks_than_cells_still_builds_landmarks(): + """Control for the test above: one landmark short of every cell is unchanged.""" + X1, X2 = _data() + model = DifferentialAbundance(n_landmarks=len(X1) + len(X2) - 1, random_state=0) + model.fit(X1, X2) + assert model.computed_landmarks.shape == (len(X1) + len(X2) - 1, 2) From 6606ee2aa0ac12d734f752dfb3631d704b4f0aff Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 12:13:57 -0700 Subject: [PATCH 19/23] docs(da): scope the #32 change in the changelog; landmark order reaches the uncertainty --- CHANGELOG.md | 10 +++++++++- kompot/differential/differential_abundance.py | 7 +++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1126187..139af23 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,13 +19,21 @@ All notable changes to this project will be documented in this file. - `da()` now honours `GPSettings.ls`, and warns when the expression-only fields `sigma`, `eps` or `use_empirical_variance` are set. + - **Abundance uncertainty no longer depends on argument order when `n_landmarks` is at least the + number of cells** (settylab/kompot#32). That is `da(gp=GPSettings(...))` on data under 5000 + cells, and the CLI whenever a GP key is set without `n_landmarks`. Each condition now gets its + full GP, as with `n_landmarks=None`. **This changes results there**: z-scores and PTPs move to the + full-GP values (on 400 + 400 cells, z by a median of 34%; the two orientations now agree exactly, + where they differed by a median of 82%), and fold changes and densities by at most 6e-4. + `da()` without `gp`, `DifferentialAbundance()`, and every run with fewer landmarks than cells are + unchanged. ### Documentation - At the default, `de()` depends on argument order: the shared length scale comes from condition 1 alone. See `DifferentialExpression.fit`. - Automatic landmarks and `"pooled"` make `da()` depend on argument order; see - `DifferentialAbundance.fit`. The instability at `n_landmarks >= n_combined` is settylab/kompot#32. + `DifferentialAbundance.fit`. ## [0.8.0] - 2026-07-28 diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 5ddcf51..041d997 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -261,8 +261,11 @@ def fit( X_condition2])``, whose **row order differs between the two orientations** of a contrast, so ``da(X, Y)`` and ``da(Y, X)`` select different landmark sets and their results stop being exact mirrors of - one another. Passing one array here to both runs removes that; so - does ``n_landmarks=None``, which builds no landmarks at all. With + one another. The uncertainty is hit hardest, because mellon's + Laplace uncertainty depends on the order of the landmark rows, not + only on the set. Passing one array here to both runs removes that; so + does ``n_landmarks=None``, or any ``n_landmarks`` at least the number + of cells, which build no landmarks at all. With ``random_state=None`` (the default) the selection is not even reproducible between two runs of the same orientation. This is independent of ``param_scheme``: see there for the other half. From 794bf7f9b817208a6da668b126522191e7c32b1b Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 12:22:17 -0700 Subject: [PATCH 20/23] ci: bound and diagnose the intermittent 3.11 hang A job timeout of 60 min, fail-fast off, and pytest-timeout at 300 s with the thread method, which dumps every thread's stack on a hang; plus --durations=25. This bounds the hang and captures evidence; it does not cure it. --- .github/workflows/tests.yml | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index e96858f..d8ff5b1 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -9,7 +9,11 @@ on: jobs: test: runs-on: ubuntu-22.04 + # A healthy job takes 20-31 min. Without this, a hung pytest holds the runner + # for GitHub's 360-min default and is reported as "cancelled". + timeout-minutes: 60 strategy: + fail-fast: false matrix: python-version: ['3.10', '3.11', '3.12'] @@ -27,7 +31,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install pytest pytest-cov + pip install pytest pytest-cov pytest-timeout # Install JAX CPU version to avoid GPU dependencies in CI pip install --upgrade "jax[cpu]" # Install package in development mode with all dependencies including dask @@ -35,7 +39,10 @@ jobs: - name: Test with pytest and generate coverage report run: | - pytest --cov=kompot --cov-report=xml --cov-report=term-missing:skip-covered + # The slowest test takes about 70 s. A test that exceeds 300 s is hung, so fail it + # and dump every thread's stack (thread method), so the log shows where it hung. + pytest --timeout=300 --timeout-method=thread --durations=25 \ + --cov=kompot --cov-report=xml --cov-report=term-missing:skip-covered - name: Upload coverage to Codecov uses: codecov/codecov-action@v5 From 31315eab627b8316fe0e3acc8eab8efd1b7aab93 Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 13:41:23 -0700 Subject: [PATCH 21/23] fix: keep main's positional order in DifferentialExpression.fit and GPSettings param_scheme was inserted ahead of landmarks in DifferentialExpression.fit and ahead of n_landmarks in GPSettings, so positional calls that work on main bound the wrong parameter. It is now the last parameter of each. DifferentialAbundance.fit keeps param_scheme in the removed sync_parameters slot, so allow_single_condition_variance keeps its position. --- .../differential/differential_expression.py | 26 +++++++------- kompot/settings.py | 30 ++++++++-------- tests/test_param_scheme.py | 34 +++++++++++++++++++ 3 files changed, 62 insertions(+), 28 deletions(-) diff --git a/kompot/differential/differential_expression.py b/kompot/differential/differential_expression.py index c34d5fa..86c0cd7 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -381,12 +381,12 @@ def fit( sigma: float = 1.0, ls: Optional[float] = None, ls_factor: float = 10.0, - param_scheme: Optional[str] = None, landmarks: Optional[np.ndarray] = None, sample_estimator_ls: Optional[float] = None, condition1_sample_indices: Optional[np.ndarray] = None, condition2_sample_indices: Optional[np.ndarray] = None, allow_single_condition_variance: bool = False, + param_scheme: Optional[str] = None, **function_kwargs, ): """ @@ -413,6 +413,18 @@ def fit( ls_factor : float, optional Multiplication factor to apply to length scale when it's automatically inferred, by default 10.0. Only used when ls is None. + landmarks : np.ndarray, optional + Pre-computed landmarks to use. If provided, n_landmarks will be ignored. + Shape (n_landmarks, n_features). + sample_estimator_ls : float, optional + Length scale for the sample-specific variance estimators. If None, will use + the same value as ls or it will be estimated, by default None. + condition1_sample_indices : np.ndarray, optional + Sample indices for first condition. Used for sample variance estimation. + Unique values in this array define different sample groups. + condition2_sample_indices : np.ndarray, optional + Sample indices for second condition. Used for sample variance estimation. + Unique values in this array define different sample groups. param_scheme : str, optional Which cells the shared length scale is estimated from when ``ls`` is None. Both conditions are normally smoothed at the *same* scale @@ -496,18 +508,6 @@ def fit( passed through ``function_kwargs``. Supplying a ``cov_func`` or ``cov_func_curry`` does *not* disable it: the resolved value is handed to the custom kernel, so both conditions keep a shared scale. - landmarks : np.ndarray, optional - Pre-computed landmarks to use. If provided, n_landmarks will be ignored. - Shape (n_landmarks, n_features). - sample_estimator_ls : float, optional - Length scale for the sample-specific variance estimators. If None, will use - the same value as ls or it will be estimated, by default None. - condition1_sample_indices : np.ndarray, optional - Sample indices for first condition. Used for sample variance estimation. - Unique values in this array define different sample groups. - condition2_sample_indices : np.ndarray, optional - Sample indices for second condition. Used for sample variance estimation. - Unique values in this array define different sample groups. **function_kwargs : dict Additional arguments to pass to the FunctionEstimator. diff --git a/kompot/settings.py b/kompot/settings.py index b34e073..baf5496 100644 --- a/kompot/settings.py +++ b/kompot/settings.py @@ -64,6 +64,20 @@ class GPSettings: it is estimated from the cells, from where ``param_scheme`` says. ls_factor : float Multiplier applied to the automatically inferred length scale. + n_landmarks : int, optional + Number of landmarks for the Nystrom approximation. + landmarks : np.ndarray, optional + Pre-computed landmark coordinates. + use_empirical_variance : bool + Estimate per-gene heteroscedastic noise from GP residuals. + batch_size : int, optional + Number of cells processed at once during prediction. + eps : float + Small constant for numerical stability. + jit_compile : bool + Use JAX JIT compilation. + random_state : int, optional + Random seed for landmark selection. param_scheme : str, optional Which cells the length scale is estimated from; in :func:`kompot.da` it covers ``d`` and ``mu`` as well. A value given explicitly (``ls`` @@ -96,26 +110,11 @@ class GPSettings: Which to pick, and the measurements behind each: :meth:`kompot.differential.DifferentialExpression.fit` and :meth:`kompot.differential.DifferentialAbundance.fit`. - n_landmarks : int, optional - Number of landmarks for the Nystrom approximation. - landmarks : np.ndarray, optional - Pre-computed landmark coordinates. - use_empirical_variance : bool - Estimate per-gene heteroscedastic noise from GP residuals. - batch_size : int, optional - Number of cells processed at once during prediction. - eps : float - Small constant for numerical stability. - jit_compile : bool - Use JAX JIT compilation. - random_state : int, optional - Random seed for landmark selection. """ sigma: float = 1.0 ls: Optional[float] = None ls_factor: float = 10.0 - param_scheme: Optional[str] = None n_landmarks: Optional[int] = 5000 landmarks: Optional[np.ndarray] = None use_empirical_variance: bool = False @@ -123,6 +122,7 @@ class GPSettings: eps: float = 1e-8 jit_compile: bool = False random_state: Optional[int] = None + param_scheme: Optional[str] = None def __post_init__(self): validate_positive_float(self.sigma, "sigma") diff --git a/tests/test_param_scheme.py b/tests/test_param_scheme.py index 2c87fd0..27c8035 100644 --- a/tests/test_param_scheme.py +++ b/tests/test_param_scheme.py @@ -857,3 +857,37 @@ def test_the_da_cli_routes_param_scheme_and_ls_into_gp_settings(): fields = {f.name for f in dataclasses.fields(GPSettings)} assert {"param_scheme", "ls"} <= GP_CONFIG_KEYS <= fields assert fields - GP_CONFIG_KEYS == {"sigma", "eps", "use_empirical_variance"} + + +# -- positional order is main's ---------------------------------------------- + + +def test_de_fit_keeps_the_positional_order_of_main(): + """``param_scheme`` is appended, so ``landmarks`` is still the 8th argument. + + Inserted ahead of ``landmarks`` it received the landmark array, and a call + that works on main died with "truth value of an array is ambiguous". + """ + X1, y1, X2, y2 = _data() + L = np.vstack([X1, X2])[::8] + + def fit(*args, **kwargs): + de = DifferentialExpression() + de.fit(*args, **kwargs) + return de.predict(np.vstack([X1, X2]), progress=False) + + positional = fit(X1, y1, X2, y2, 1.0, None, 10.0, L) + keyword = fit(X1, y1, X2, y2, sigma=1.0, ls=None, ls_factor=10.0, landmarks=L) + for key in ("fold_change", "fold_change_zscores"): + np.testing.assert_array_equal(positional[key], keyword[key]) + + +def test_gpsettings_keeps_the_positional_order_of_main(): + """``param_scheme`` is the LAST field, so main's positional order still binds.""" + g = GPSettings(1.0, None, 10.0, 100) + assert g.n_landmarks == 100 and g.param_scheme is None + import dataclasses + + main = ["sigma", "ls", "ls_factor", "n_landmarks", "landmarks", + "use_empirical_variance", "batch_size", "eps", "jit_compile", "random_state"] + assert [f.name for f in dataclasses.fields(GPSettings)][: len(main)] == main From 343486fd6ec85d1a24b773b550a745880d058535 Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 13:41:24 -0700 Subject: [PATCH 22/23] docs(da): the full-GP fallback equals n_landmarks=None only up to 5000 cells per condition --- CHANGELOG.md | 3 ++- kompot/differential/differential_abundance.py | 6 ++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 139af23..e91d7c9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,7 +22,8 @@ All notable changes to this project will be documented in this file. - **Abundance uncertainty no longer depends on argument order when `n_landmarks` is at least the number of cells** (settylab/kompot#32). That is `da(gp=GPSettings(...))` on data under 5000 cells, and the CLI whenever a GP key is set without `n_landmarks`. Each condition now gets its - full GP, as with `n_landmarks=None`. **This changes results there**: z-scores and PTPs move to the + full GP (the same as `n_landmarks=None` up to 5000 cells per condition, where mellon's own + landmark default starts). **This changes results there**: z-scores and PTPs move to the full-GP values (on 400 + 400 cells, z by a median of 34%; the two orientations now agree exactly, where they differed by a median of 82%), and fold changes and densities by at most 6e-4. `da()` without `gp`, `DifferentialAbundance()`, and every run with fewer landmarks than cells are diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 041d997..acfa3f9 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -155,7 +155,8 @@ def __init__( n_landmarks : int, optional Number of landmarks to use for approximation. If None, use all points, by default None. At or above the number of cells in both conditions together, no landmarks - are built and each condition gets its full GP, as with None. + are built and each condition gets its full GP. That equals None up to 5000 cells + per condition; above that, None leaves mellon at its default of 5000 landmarks. use_sample_variance : bool, optional Whether to use sample variance for uncertainty estimation. By default None. - If None (recommended): Automatically determined based on variance_predictor1/2 @@ -379,7 +380,8 @@ def fit( # mellon's Laplace uncertainty depends on the landmarks' row order, # so da(X, Y) and da(Y, X) then disagree on the uncertainty # (settylab/kompot#32). No approximation was asked for: fit each - # condition's exact GP, on its own cells, as n_landmarks=None does. + # condition's exact GP, on its own cells (n_landmarks=0 to mellon; the + # same as n_landmarks=None only up to mellon's 5000-landmark default). logger.info( f"n_landmarks={self.n_landmarks:,} >= {len(X_condition1) + len(X_condition2):,} " "cells: fitting each condition's full GP instead of landmarks." From da733b46c093c6baff68988dcda399ee7674bbed Mon Sep 17 00:00:00 2001 From: Dominik Date: Wed, 23 Sep 2026 13:56:23 -0700 Subject: [PATCH 23/23] fix(cli): kompot dm validates its input before requiring Palantir run_dm imported Palantir first, so without it installed a missing PCA key was never reported and test_dm_missing_pca failed in any environment without the optional dependency. CI installs .[all], which includes Palantir, so it passed there. The import now follows the input checks; a new test covers the Palantir-missing error with valid input. --- kompot/cli/dm.py | 17 +++++++++-------- tests/test_cli.py | 16 ++++++++++++++++ 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/kompot/cli/dm.py b/kompot/cli/dm.py index b3515b0..6712247 100644 --- a/kompot/cli/dm.py +++ b/kompot/cli/dm.py @@ -131,14 +131,6 @@ def run_dm(args): args Parsed arguments from argparse """ - # Check if palantir is available - try: - import palantir - except ImportError: - logger.error("Palantir is not installed. Install it with: pip install palantir") - logger.error("Or install kompot with: pip install kompot[recommended]") - sys.exit(1) - # Validate input file input_path = validate_anndata_path(args.input) @@ -188,6 +180,15 @@ def run_dm(args): logger.error("Run PCA first using scanpy: sc.pp.pca(adata)") sys.exit(1) + # Palantir is checked only now, so an input problem is reported whether or not + # the optional dependency is installed. + try: + import palantir + except ImportError: + logger.error("Palantir is not installed. Install it with: pip install palantir") + logger.error("Or install kompot with: pip install kompot[recommended]") + sys.exit(1) + logger.info("Computing diffusion maps") logger.info(f" PCA key: {pca_key}") logger.info(f" n_components: {params.get('n_components', 10)}") diff --git a/tests/test_cli.py b/tests/test_cli.py index 694301f..eb3a6fa 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -763,6 +763,22 @@ def test_dm_missing_pca(self, sample_adata, temp_dir): or "PCA coordinates not found" in result.stdout ) + def test_dm_without_palantir_says_so(self, sample_adata, temp_dir): + """With valid input and no Palantir, the error names the missing package.""" + try: + import palantir # noqa: F401 + except ImportError: + pass + else: + pytest.skip("palantir is installed") + input_file = Path(temp_dir) / "input.h5ad" + sample_adata.write_h5ad(input_file) + cmd = [sys.executable, "-m", "kompot.cli", "dm", str(input_file), + "-o", str(Path(temp_dir) / "output.h5ad")] + result = subprocess.run(cmd, capture_output=True, text=True) + assert result.returncode != 0 + assert "Palantir is not installed" in result.stderr + result.stdout + class TestCLIMainEntry: """Test main CLI entry point."""