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 diff --git a/CHANGELOG.md b/CHANGELOG.md index a336fd0..e91d7c9 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.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. + +### Removed + + - `sync_parameters` on `DifferentialAbundance.fit` and `da()`; use `param_scheme="pooled"`. + +### Fixed + + - `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 (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 + 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`. + ## [0.8.0] - 2026-07-28 ### Changed — statistics now match the manuscript diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 36dc71a..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 - * - ``sigma``, ``ls``, ``ls_factor``, ``n_landmarks``, ``batch_size``, ``eps``, ``jit_compile``, ``random_state``, ``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 + - 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/anndata/_da_helpers.py b/kompot/anndata/_da_helpers.py index 4c61b6e..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,6 +27,7 @@ def _check_da_overwrites( condition2: str, obsm_key: str, ls_factor: float, + param_scheme: str, ): """Check for existing DA results and handle overwrite logic.""" from .utils import detect_output_field_overwrite @@ -101,9 +104,13 @@ 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. + 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 042c3b0..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,6 +42,7 @@ def _check_overwrites( obsm_key: str, layer: Optional[str], ls_factor: float, + param_scheme: str, ): """Check for existing results and handle overwrite logic. @@ -134,9 +136,15 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", + "param_scheme", ]: curr_val = locals().get(param_name) - prev_val = _pg(prev_params, param_name) + # A run recorded before `param_scheme` existed used the default. + prev_val = _pg( + prev_params, + param_name, + DEFAULT_PARAM_SCHEME if param_name == "param_scheme" else None, + ) if curr_val != prev_val: params_match = False logger.debug( @@ -189,9 +197,15 @@ def _check_overwrites( "obsm_key", "layer", "ls_factor", + "param_scheme", ]: curr_val = locals().get(param_name) - prev_val = _pg(prev_params, param_name) + # A run recorded before `param_scheme` existed used the default. + prev_val = _pg( + prev_params, + param_name, + DEFAULT_PARAM_SCHEME 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..577032d 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, @@ -76,9 +78,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 @@ -90,7 +99,8 @@ def da( ``density_predictor1/2`` and ``variance_predictor1/2`` are used. See :class:`~kompot.ModelSettings`. **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 +117,33 @@ 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 "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 = [ + 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 + + 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 +194,7 @@ def da( condition2=condition2, obsm_key=obsm_key, ls_factor=ls_factor, + param_scheme=param_scheme, ) # ---- 1. Copy if requested ---- @@ -199,6 +237,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 +274,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 +297,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 4315a0d..74eb2c3 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, @@ -97,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 @@ -132,6 +138,12 @@ 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 use_empirical_variance = _gp.use_empirical_variance @@ -305,6 +317,7 @@ def de( obsm_key=obsm_key, layer=layer, ls_factor=ls_factor, + param_scheme=param_scheme, ) # ---- 2. Copy if requested ---- @@ -455,6 +468,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, + param_scheme=param_scheme, landmarks=landmarks, condition1_sample_indices=condition1_sample_indices, condition2_sample_indices=condition2_sample_indices, @@ -635,6 +649,7 @@ def de( sigma=sigma, ls=ls, ls_factor=ls_factor, + param_scheme=param_scheme, n_landmarks=n_landmarks, use_empirical_variance=use_empirical_variance, batch_size=batch_size, diff --git a/kompot/anndata/utils/params.py b/kompot/anndata/utils/params.py index 4a46042..545a434 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"), + "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..ce2b3d9 100644 --- a/kompot/cli/da.py +++ b/kompot/cli/da.py @@ -184,6 +184,22 @@ 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`. +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 +308,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 5e6f80c..5ee951e 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", + "param_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/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/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 80abd92..6697c4a 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 +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 # Memory management: diff --git a/kompot/differential/differential_abundance.py b/kompot/differential/differential_abundance.py index 17ab1ec..acfa3f9 100644 --- a/kompot/differential/differential_abundance.py +++ b/kompot/differential/differential_abundance.py @@ -9,10 +9,88 @@ 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=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"``. ``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)) + # 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, d=None +): + """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, d)[0] + if param_scheme == "condition2": + 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 + # 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, 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). + 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. + 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. + 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 +99,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 +122,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. """ @@ -76,6 +154,9 @@ 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. 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 @@ -157,7 +238,7 @@ 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, + param_scheme: Optional[str] = None, allow_single_condition_variance: bool = False, **density_kwargs, ): @@ -176,6 +257,19 @@ 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. 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. 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. @@ -188,11 +282,58 @@ 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 - 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. + param_scheme : str, optional + Where the density hyperparameters ``d``, ``mu`` and ``ls`` come from. + 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. + * ``"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, ``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 + 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 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. + + ``"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. + 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. @@ -202,6 +343,8 @@ def fit( The fitted instance. """ + 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 @@ -229,11 +372,23 @@ 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 (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." + ) + 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, @@ -246,8 +401,8 @@ 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. + if param_scheme == "pooled": # Combine data from both conditions for parameter estimation X_combined = np.vstack([X_condition1, X_condition2]) logger.info( @@ -278,6 +433,34 @@ 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, + d=density_kwargs.get("d"), + ) + 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 6af2bee..86c0cd7 100644 --- a/kompot/differential/differential_expression.py +++ b/kompot/differential/differential_expression.py @@ -10,10 +10,81 @@ 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") +#: 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: + """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_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 param_scheme == "condition1": + return None, None + 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 + # 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"param_scheme='condition2': shared length scale {shared:.4f}") + return shared, shared + if param_scheme == "separate": + return None, None + if param_scheme == "pooled": + shared = _auto_ls( + np.vstack([np.asarray(X_condition1), np.asarray(X_condition2)]), + ls_factor, + ) + logger.info(f"param_scheme='pooled': shared length scale {shared:.4f}") + return shared, shared + 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) + 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"param_scheme='symmetric': condition length scales {ls1:.4f} / " + f"{ls2:.4f} -> shared {shared:.4f}" + ) + return shared, shared + raise ValueError( + f"Unknown param_scheme {param_scheme!r}. Expected one of {PARAM_SCHEMES}." + ) + class DifferentialExpression: """ @@ -315,6 +386,7 @@ def fit( 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, ): """ @@ -336,7 +408,8 @@ 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 ``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. @@ -352,6 +425,89 @@ def fit( 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 + so that their + fitted surfaces are comparable; the schemes differ in which cells the + 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 + 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, 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:: + + **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 against + ``assert_allclose(rtol=1e-6, atol=1e-8)`` — the tolerance this + 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 + ``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** + ``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 + 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 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. **function_kwargs : dict Additional arguments to pass to the FunctionEstimator. @@ -361,6 +517,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 @@ -394,6 +552,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}") @@ -408,6 +567,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 ( + param_scheme == "condition2" + and not landmarks_were_provided + and self.n_landmarks is not None + and self.n_landmarks > 0 + ): + logger.warning( + "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 " + "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( @@ -420,13 +605,33 @@ def fit( disk_storage_dir=self.disk_storage_dir, ) + # -- Resolve the length scale each condition will be fitted with -- + # 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 `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_scheme_ls( + param_scheme, X_condition1, X_condition2, ls_factor + ) + else: + if param_scheme != "condition1": + logger.info( + f"param_scheme={param_scheme!r} ignored: the length scale is " + "already fixed by an explicit ls." + ) + 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 +642,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 param_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..baf5496 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,8 +60,8 @@ class GPSettings: sigma : float Noise level for the GP. ls : float, optional - Length scale. If *None*, estimated automatically using - ``ls_factor``. + 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. n_landmarks : int, optional @@ -56,6 +78,38 @@ class GPSettings: 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`` + 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")`` 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`. """ sigma: float = 1.0 @@ -68,11 +122,13 @@ 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") validate_positive_float(self.ls, "ls", optional=True) validate_positive_float(self.ls_factor, "ls_factor") + 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_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.""" 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", ) diff --git a/tests/test_da_swap_equivalence.py b/tests/test_da_swap_equivalence.py new file mode 100644 index 0000000..15bd9bc --- /dev/null +++ b/tests/test_da_swap_equivalence.py @@ -0,0 +1,311 @@ +"""Swap equivalence for differential ABUNDANCE at its default ``param_scheme``. + +``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 -- + + density_estimator_condition1 = mellon.DensityEstimator(**estimator_defaults) + density_estimator_condition2 = mellon.DensityEstimator(**estimator_defaults) + +-- 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 ``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 +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 ``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`` +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. +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, "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. + +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 +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, 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, param_scheme=param_scheme, landmarks=landmarks) + + rev = DifferentialAbundance(n_landmarks=n_landmarks, random_state=0) + rev.fit(X2, X1, param_scheme=param_scheme, 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(param_scheme="pooled"), + "d / mu / ls synchronised from the stacked union", + id="pooled", + ), + ], +) +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 ``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. + + 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(): + """Share the landmarks, and the automatic-landmark half is exact again. + + 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 + ``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 + + 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 + ) + + +# -- 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) 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_param_scheme.py b/tests/test_param_scheme.py new file mode 100644 index 0000000..27c8035 --- /dev/null +++ b/tests/test_param_scheme.py @@ -0,0 +1,893 @@ +"""``param_scheme``: where the hyperparameters of a two-condition fit come 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``, 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. +""" + +import contextlib +import logging + +import numpy as np +import pytest + +from kompot.differential.differential_expression import ( + DifferentialExpression, + _auto_ls, + _resolve_scheme_ls, +) +from kompot.settings import PARAM_SCHEMES, 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_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_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 -- + # 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_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] + 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_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_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. + assert fwd == pytest.approx(rev, rel=1e-2) + + +def test_unknown_scheme_is_rejected(): + X1, _, X2, _ = _data() + 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="param_scheme"): + GPSettings(param_scheme="nope") + for scheme in PARAM_SCHEMES: + GPSettings(param_scheme=scheme) + + +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") + + +# -------------------------------------------------------------------------- +# 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, param_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", PARAM_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) + + +# -------------------------------------------------------------------------- +# 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, param_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 `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. + """ + 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) + + +# -------------------------------------------------------------------------- +# "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, 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_scheme_ls("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_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) + + +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, param_scheme="condition1", landmarks=landmarks) + + mir = DifferentialExpression(n_landmarks=n_landmarks, random_state=0) + mir.fit(X2, y2, X1, y1, param_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 `param_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 `param_scheme`, not just the writer +# -------------------------------------------------------------------------- + + +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` + 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 PARAM_SCHEMES: + params = build_params_dict( + {"groupby": "g", "condition1": "A", "condition2": "B"}, + gp=GPSettings(param_scheme=scheme), + ) + assert params["gp"]["param_scheme"] == scheme + assert params_get(params, "param_scheme") == scheme + + +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 + 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 "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"} + + +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 + 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() + + with _captured_kompot_warnings() as messages: + DifferentialExpression(n_landmarks=20, random_state=0).fit( + X1, y1, X2, y2, param_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, param_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, 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 anndata # noqa: E402 + +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 + +# 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) + 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 -------------------------------------------------------------- + + +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]) + + +def test_pooled_and_separate_differ(): + 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_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) + + +_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 --------------------------------------------------------- + + +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 + # 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 + 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)) + assert np.abs(pooled - separate).max() > 1e-2 + + +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) + + +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)) + 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"} + + +# -- 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