feat(de): add ls_scheme to choose how the shared length scale is estimated

CHANGELOG.md

@@ -2,6 +2,40 @@
 
 All notable changes to this project will be documented in this file.
 
+## [Unreleased]
+
+### Added
+
+ - **`GPSettings.ls_scheme`** — chooses which cells the automatic shared length scale is
+   estimated from in differential expression. The default, `"condition1"`, is exactly the
+   existing behaviour and nothing changes unless you ask for something else.
+
+   The two conditions have always been smoothed at a **shared** length scale so that their
+   fitted surfaces are comparable, and that part is right: giving each condition its own
+   (`ls_scheme="separate"`) makes the two surfaces differently smooth, and the mismatch alone
+   manufactures apparent fold changes. What was undocumented is that the shared value is
+   estimated from **condition 1's cells only**, so it is a function of `n_condition1` and
+   `de(adata, condition1=X, condition2=Y)` is *not* equivalent to
+   `de(adata, condition1=Y, condition2=X)`. On an exchangeable null with nothing to find, the
+   two orientations of the same partition differ in false-positive rate by 0.06–0.08 at
+   2:1 and beyond.
+
+   `"symmetric"` shares the size-weighted geometric mean of the two conditions' own estimates
+   — the same estimator, with nearest neighbours looked up *within* each condition — and is
+   invariant under swapping the conditions. `"pooled"` estimates from the union, which is also
+   swap-invariant but yields a systematically smaller value, because the union is denser than
+   either condition and nearest-neighbour distances shrink with cell count alone.
+
+### Changed
+
+ - `ls_scheme` participates in run-parameter matching, so re-running under a different scheme
+   is no longer treated as a matching rerun.
+
+### Documentation
+
+ - `GPSettings.ls` and `DifferentialExpression.fit` now state where the shared length scale
+   comes from and that the default makes the contrast depend on argument order.
+
 ## [0.8.0] - 2026-07-28
 
 ### Changed — statistics now match the manuscript

kompot/anndata/_de_helpers.py

@@ -41,6 +41,7 @@ def _check_overwrites(
     obsm_key: str,
     layer: Optional[str],
     ls_factor: float,
+    ls_scheme: str = "condition1",
 ):
     """Check for existing results and handle overwrite logic.
 
@@ -134,6 +135,7 @@ def _check_overwrites(
                 "obsm_key",
                 "layer",
                 "ls_factor",
+                "ls_scheme",
             ]:
                 curr_val = locals().get(param_name)
                 prev_val = _pg(prev_params, param_name)
@@ -189,6 +191,7 @@ def _check_overwrites(
                 "obsm_key",
                 "layer",
                 "ls_factor",
+                "ls_scheme",
             ]:
                 curr_val = locals().get(param_name)
                 prev_val = _pg(prev_params, param_name)

kompot/anndata/differential_expression.py

@@ -132,6 +132,7 @@ def de(
     sigma = _gp.sigma
     ls = _gp.ls
     ls_factor = _gp.ls_factor
+    ls_scheme = _gp.ls_scheme
     n_landmarks = _gp.n_landmarks
     landmarks = _gp.landmarks
     use_empirical_variance = _gp.use_empirical_variance
@@ -305,6 +306,7 @@ def de(
         obsm_key=obsm_key,
         layer=layer,
         ls_factor=ls_factor,
+        ls_scheme=ls_scheme,
     )
 
     # ---- 2. Copy if requested ----
@@ -455,6 +457,7 @@ def de(
         sigma=sigma,
         ls=ls,
         ls_factor=ls_factor,
+        ls_scheme=ls_scheme,
         landmarks=landmarks,
         condition1_sample_indices=condition1_sample_indices,
         condition2_sample_indices=condition2_sample_indices,
@@ -635,6 +638,7 @@ def de(
                 sigma=sigma,
                 ls=ls,
                 ls_factor=ls_factor,
+                ls_scheme=ls_scheme,
                 n_landmarks=n_landmarks,
                 use_empirical_variance=use_empirical_variance,
                 batch_size=batch_size,
@@ -730,6 +734,7 @@ def compute_differential_expression(
     sigma: float = 1.0,
     ls=None,
     ls_factor: float = 10.0,
+    ls_scheme: str = "condition1",
     compute_mahalanobis: bool = True,
     jit_compile: bool = False,
     eps: float = 1e-8,
@@ -783,6 +788,7 @@ def compute_differential_expression(
             sigma=sigma,
             ls=ls,
             ls_factor=ls_factor,
+            ls_scheme=ls_scheme,
             n_landmarks=n_landmarks,
             landmarks=landmarks,
             use_empirical_variance=use_empirical_variance,

kompot/differential/differential_expression.py

@@ -14,6 +14,61 @@
 
 logger = logging.getLogger("kompot")
 
+#: Recognised values for ``ls_scheme``.  See :meth:`DifferentialExpression.fit`.
+LS_SCHEMES = ("condition1", "symmetric", "pooled", "separate")
+
+
+def _auto_ls(X: np.ndarray, ls_factor: float) -> float:
+    """Automatic length scale for one design matrix.
+
+    Identical to what :meth:`ExpressionModel.fit` derives internally when
+    ``ls`` is None, exposed here so a scheme can combine the two conditions'
+    values *before* either model is fitted.
+    """
+    from mellon.parameters import compute_ls, compute_nn_distances
+
+    return float(compute_ls(compute_nn_distances(np.asarray(X))) * ls_factor)
+
+
+def _resolve_ls_scheme(ls_scheme, X_condition1, X_condition2, ls_factor):
+    """Return the ``(ls_model1, ls_model2)`` a scheme prescribes.
+
+    ``None`` in a slot means "let that model estimate its own".  For
+    ``"condition1"`` both slots are None and the caller reproduces the
+    historical behaviour by copying model1's fitted value into model2.
+    """
+    if ls_scheme == "condition1":
+        return None, None
+    if ls_scheme == "separate":
+        return None, None
+    if ls_scheme == "pooled":
+        shared = _auto_ls(
+            np.vstack([np.asarray(X_condition1), np.asarray(X_condition2)]),
+            ls_factor,
+        )
+        logger.info(f"ls_scheme='pooled': shared length scale {shared:.4f}")
+        return shared, shared
+    if ls_scheme == "symmetric":
+        n1 = np.asarray(X_condition1).shape[0]
+        n2 = np.asarray(X_condition2).shape[0]
+        ls1 = _auto_ls(X_condition1, ls_factor)
+        ls2 = _auto_ls(X_condition2, ls_factor)
+        # Size-weighted GEOMETRIC mean.  kompot's estimator is
+        # exp(mean(log nn_dist)) * const, so this is exactly that estimator
+        # applied to the two conditions' WITHIN-condition nearest-neighbour
+        # distances pooled together -- unlike 'pooled', which looks the
+        # neighbours up across the union and therefore shrinks purely because
+        # the union holds more cells.
+        shared = float(np.exp((n1 * np.log(ls1) + n2 * np.log(ls2)) / (n1 + n2)))
+        logger.info(
+            f"ls_scheme='symmetric': condition length scales {ls1:.4f} / "
+            f"{ls2:.4f} -> shared {shared:.4f}"
+        )
+        return shared, shared
+    raise ValueError(
+        f"Unknown ls_scheme {ls_scheme!r}. Expected one of {LS_SCHEMES}."
+    )
+
 
 class DifferentialExpression:
     """
@@ -310,6 +365,7 @@ def fit(
         sigma: float = 1.0,
         ls: Optional[float] = None,
         ls_factor: float = 10.0,
+        ls_scheme: str = "condition1",
         landmarks: Optional[np.ndarray] = None,
         sample_estimator_ls: Optional[float] = None,
         condition1_sample_indices: Optional[np.ndarray] = None,
@@ -336,10 +392,40 @@ def fit(
         sigma : float, optional
             Noise level for function estimator, by default 1.0.
         ls : float, optional
-            Length scale for the GP kernel. If None, it will be estimated, by default None.
+            Length scale for the GP kernel, shared by both conditions. If None it is
+            estimated from the data according to ``ls_scheme``, by default None.
         ls_factor : float, optional
             Multiplication factor to apply to length scale when it's automatically inferred,
             by default 10.0. Only used when ls is None.
+        ls_scheme : str, optional
+            How the automatic length scale is derived when ``ls`` is None. Both
+            conditions are normally smoothed at the *same* scale so that their
+            fitted surfaces are comparable; the schemes differ in which cells the
+            shared value is estimated from.
+
+            * ``"condition1"`` (default) — estimate from condition 1's cells and
+              reuse the value for condition 2. **Not symmetric**: the length
+              scale is a function of ``n_condition1`` alone, so swapping
+              ``condition1`` and ``condition2`` changes the result even when
+              nothing else does.
+            * ``"symmetric"`` — estimate a length scale from each condition
+              separately and share their size-weighted geometric mean. Invariant
+              under swapping the two conditions, and it does not inherit the
+              cell-count artefact of ``"pooled"``.
+            * ``"pooled"`` — estimate from the two conditions' cells taken
+              together. Also swap-invariant, but the union is denser than either
+              condition, so nearest-neighbour distances — and hence the length
+              scale — shrink purely because there are more cells.
+            * ``"separate"`` — let each condition estimate its own length scale
+              and do not share. The two surfaces are then smoothed differently
+              and the mismatch itself produces apparent fold changes; measured on
+              an exchangeable null this is markedly worse than any shared value.
+              Provided for diagnostics, not recommended.
+
+            Ignored when ``ls`` is given explicitly, or when a length scale is
+            passed through ``function_kwargs``. Supplying a ``cov_func`` or
+            ``cov_func_curry`` does *not* disable it: the resolved value is
+            handed to the custom kernel, so both conditions keep a shared scale.
         landmarks : np.ndarray, optional
             Pre-computed landmarks to use. If provided, n_landmarks will be ignored.
             Shape (n_landmarks, n_features).
@@ -420,13 +506,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 `ls_scheme`,
+        # silently giving them "separate" behaviour.
+        ls_auto = ls is None and "ls" not in function_kwargs
+        if ls_auto:
+            ls_model1, ls_model2 = _resolve_ls_scheme(
+                ls_scheme, X_condition1, X_condition2, ls_factor
+            )
+        else:
+            if ls_scheme != "condition1":
+                logger.info(
+                    f"ls_scheme={ls_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 +543,15 @@ def fit(
                 **function_kwargs,
             )
 
-        # Extract ls from model1 for model2 consistency
-        ls_for_model2 = ls
-        if ls is None and "ls" not in function_kwargs and self.model1.ls is not None:
+        # Under the default scheme condition 2 inherits condition 1's fitted
+        # length scale, so the shared value is a function of n_condition1 alone.
+        # Every other scheme has already fixed both values above.
+        ls_for_model2 = ls_model2
+        if (
+            ls_auto
+            and ls_scheme == "condition1"
+            and self.model1.ls is not None
+        ):
             ls_for_model2 = self.model1.ls
 
         # -- Fit model2 --

kompot/settings.py

@@ -38,10 +38,19 @@ class GPSettings:
     sigma : float
         Noise level for the GP.
     ls : float, optional
-        Length scale.  If *None*, estimated automatically using
-        ``ls_factor``.
+        Length scale, shared by both conditions.  If *None*, estimated
+        automatically using ``ls_factor`` and ``ls_scheme``.
     ls_factor : float
         Multiplier applied to the automatically inferred length scale.
+    ls_scheme : str
+        Which cells the automatic shared length scale is estimated from.
+        ``"condition1"`` (default) uses condition 1's cells only, which makes
+        the result depend on which condition is passed first; ``"symmetric"``
+        shares the size-weighted geometric mean of the two per-condition
+        estimates; ``"pooled"`` estimates from both conditions' cells taken
+        together; ``"separate"`` gives each condition its own.  See
+        :meth:`kompot.differential.DifferentialExpression.fit`.  Ignored when
+        ``ls`` is given explicitly.
     n_landmarks : int, optional
         Number of landmarks for the Nystrom approximation.
     landmarks : np.ndarray, optional
@@ -61,6 +70,7 @@ class GPSettings:
     sigma: float = 1.0
     ls: Optional[float] = None
     ls_factor: float = 10.0
+    ls_scheme: str = "condition1"
     n_landmarks: Optional[int] = 5000
     landmarks: Optional[np.ndarray] = None
     use_empirical_variance: bool = False
@@ -73,6 +83,13 @@ def __post_init__(self):
         validate_positive_float(self.sigma, "sigma")
         validate_positive_float(self.ls, "ls", optional=True)
         validate_positive_float(self.ls_factor, "ls_factor")
+        from .differential.differential_expression import LS_SCHEMES
+
+        if self.ls_scheme not in LS_SCHEMES:
+            raise ValueError(
+                f"'ls_scheme' must be one of {LS_SCHEMES} (got "
+                f"{self.ls_scheme!r})."
+            )
         validate_positive_int(self.n_landmarks, "n_landmarks", optional=True)
         validate_bool(self.use_empirical_variance, "use_empirical_variance")
         validate_positive_int(self.batch_size, "batch_size", optional=True)